Efficient Simulation of Dynamic Flux Balance Models

Dynamic flux balance modeling is a powerful tool for the analysis and optimization of biochemical reactors. Dynamic flux balance analysis (DFBA) assumes that species within a cell rapidly reach equilibrium with the extracellular environment. Intracellular metabolism is then coupled with the dynamic mass balance equations for the extracellular environment via expressions for the rates of substrate uptake and product excretion. The model for the intracellular metabolism is a stoichiometric model that relates species and fluxes via a reaction network. However, this linear system is underdetermined since there are more fluxes than species. This is overcome by assuming that the cell will maximize some objective within the imposed flux constraints. If this objective, usually cell growth, can be modeled as a linear function, the result is a linear program (LP). The overall dynamic flux balance model is an ordinary differential equation (ODE) with a LP "embedded''.

The general form is

x'(t) = f(t,x(t),u(t,x(t))),  x(t₀) = x₀

with {u(t,z)} = arg min_v c^Tv                                                                             (1)

                        s.t. Av = b(t,z),

                                    v ≥ 0.

These models have proven useful and accurate, and one could use these models to perform dynamic optimization of bioreactor systems. However, solving (1) by embedding a LP solver in a code for numerically integrating ODEs leads to quite an inefficient and inaccurate implementation, especially if the dynamics are stiff, the LP is large and/or several LPs modeling several micro-organisms are embedded. The issue at hand is that the vector field, which depends on the optimal solution to the LP, is nonsmooth. In addition, solving a LP is normally a more computationally complex task than a simple function evaluation and this induces numerical "noise''. On the other hand, the algorithm described in this work overcomes both of these issues, and the result is an efficient solution procedure for problems of the form (1).

A code implementing the proposed algorithm has been written, called DSL48LPR. In the most general sense, DSL48LPR solves the LP once using the simplex algorithm and records the basis set B, an index set which partitions the primal variables v. The solution set (while this basis is optimal) is obtained very simply. The basic variables v_B of the solution set are given by the solution to A_Bv_B = b(t,z) (where A_B is a submatrix constructed from the columns of A corresponding to B) and the nonbasic variables are simply zero. For the LP in (1) a basis set is optimal until it becomes infeasible;  since the LP is in standard from, this is indicated by a variable crossing zero. Using the event detection capabilities of DAEPACK component DSL48E, integration stops when a variable crosses zero. The basis set is updated by re-solving the LP before continuing integration. Thus, for most integration steps, a system of linear equations is solved, rather than a linear program. Overall, (1) is reformulated as differential-algebraic equations (DAEs) with discrete modes. This framework provides an efficient and accurate method for solving a system of ODEs with embedded LPs.