DFBAlab: a fast and reliable MATLAB code for dynamic flux balance analysis

BMC Bioinformatics DFBAlab: a fast and reliable MATLAB code for dynamic flux balance analysis Jose A Gomez 0 Kai Hffner 0 Paul I Barton 0 0 Process Systems Engineering Laboratory, Massachusetts Institute of Technology , Cambridge, MA 02139 , USA Background: Dynamic Flux Balance Analysis (DFBA) is a dynamic simulation framework for biochemical processes. DFBA can be performed using different approaches such as static optimization (SOA), dynamic optimization (DOA), and direct approaches (DA). Few existing simulators address the theoretical and practical challenges of nonunique exchange fluxes or infeasible linear programs (LPs). Both are common sources of failure and inefficiencies for these simulators. Results: DFBAlab, a MATLAB-based simulator that uses the LP feasibility problem to obtain an extended system and lexicographic optimization to yield unique exchange fluxes, is presented. DFBAlab is able to simulate complex dynamic cultures with multiple species rapidly and reliably, including differential-algebraic equation (DAE) systems. In addition, DFBAlab?s running time scales linearly with the number of species models. Three examples are presented where the performance of COBRA, DyMMM and DFBAlab are compared. Conclusions: Lexicographic optimization is used to determine unique exchange fluxes which are necessary for a well-defined dynamic system. DFBAlab does not fail during numerical integration due to infeasible LPs. The extended system obtained through the LP feasibility problem in DFBAlab provides a penalty function that can be used in optimization algorithms. Dynamic flux balance analysis; Nonsmooth dynamic systems; Linear programming; Lexicographic optimization - Background The acceleration in the process of genome sequencing in recent years has increased the availability of genomescale metabolic network reconstructions for a variety of species. These genome-based networks can be used within the framework of flux balance analysis (FBA) to predict steady-state growth and uptake rates accurately [1]. Dynamic flux balance analysis (DFBA) enables the simulation of dynamic biological systems by assuming organisms reach steady state rapidly in response to changes in the extracellular environment. Then, the rates predicted by FBA are used to update the extracellular environment. There exist three approaches to simulate DFBA models: the static optimization approach (SOA) [2], the dynamic optimization approach [2] (DOA), and the direct approach (DA). The static optimization approach uses the Euler forward method, solving the embedded LPs at each time step. Since most DFBA models are stiff, small time steps are required for stability, making this approach computationally expensive. Meanwhile, the DOA approach discretizes the time horizon and optimizes simultaneously over the entire time period of interest by solving a nonlinear programming problem (NLP). The dimension of this NLP increases with time discretization, therefore it is limited to small-scale metabolic models [3]. Finally, a DA has been proposed recently by including the LP solver in the right-hand side evaluator for the ordinary differential equations (ODEs) and taking advantage of reliable implicit ODE integrators with adaptive step size for error control. At present, the DOA is rarely used due to the intractability of the resulting NLP. DFBA can be easily performed on MATLAB using the constraint-based reconstruction and analysis (COBRA) toolbox [1,4], which implements the SOA. Recently, the DA has been implemented by Hanly and Henson [5], Mao and Verwoerd in the ORCA toolbox [6], Zhuang et al. in the dynamic multispecies metabolic modeling (DyMMM) framework [7,8], and others. A comprehensive list of DFBA implementations can be found in Table I of [3]. COBRA, DyMMM and ORCA codes are available on the Web. Of these, only DyMMM allows community simulations. Since ORCA and DyMMM are extremely similar, only COBRA and DyMMM were implemented in the case studies presented. These implementations present several shortcomings. The COBRA Toolbox uses a fixed time step and does not take advantage of the high quality built-in integrators provided by MATLAB. Simulation stability and accuracy are closely linked to a uniformly small step size which can greatly increase simulation time. It can fail if the extracellular conditions are close to the FBA model becoming infeasible. In addition, it uses a simple exchange flux bounding scheme that does not allow the implementation of Michaelis-Menten kinetics or other more complex dynamic behaviors such as day/night shifts for photosynthetic organisms or system feed and discharge rates. It does not allow community simulations. The ORCA toolbox and the DyMMM framework use the MATLAB built-in integrators. ORCA simulates monocultures only, whereas DyMMM can simulate cocultures. The ORCA toolbox allows the implementation of Michaelis-Menten and Hill kinetics only, whereas DyMMM provides the flexibility to implement more complex dynamics such as day/night shifts for photosynthetic organisms or system feed and discharge rates. Both attempt to carry on with simulations when the FBA model is infeasible by setting the growth rate and exchange fluxes equal to zero and displaying a death phase message. This message may be displayed even though the system is not infeasible. None of these implementations (COBRA, ORCA, and DyMMM) account for the solution of a linear program (LP) being a nonsingleton set. Therefore, exchange fluxes are not necessarily unique and the dynamic system is not well-defined. Nonunique optimal fluxes have been discussed elsewhere in [9] and [5]. If no effort is made to obtain unique fluxes, different integrators could yield different results. Hffner et al. have designed a fast and reliable community simulator that has the flexibility of implementing complex dynamics, does not fail, identifies precisely when a system becomes infeasible, and performs lexicographic optimization to render unique exchange fluxes [3]. In particular, it avoids numerical failure by reformulating the LP as an algebraic system and integrating an index-1 differential-algebraic equation (DAE) system. Despite these advantages, this simulator has not been widely used due to being coded in FORTRAN. In this paper, we implement the LP feasibility problem combined with lexicographic optimization in our Dynamic vkUB(x(t)) v vkLB(x(t)), where ck Rnrk is the cost vector that maximizes growth fluxes, and vkLB, vkUB are lower and upper bounds as functions of the extracellular concentrations. The vector hk then takes the solution of this linear program to find the values of the exchange fluxes (e.g. biomass production rate, O2 consumption rate, ethanol production rate, etc.). This definition of DFBA has a serious problem: the solution set of the LP (2) can be nonunique (e.g. different flux distributions v can attain the maximum growth rate) and Flux Balance Analysis laboratory (DFBAlab), a MATLAB c (...truncated)