(313e) 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.