# 15. Dynamic Flux Balance Analysis (dFBA) in COBRApy¶

The following notebook shows a simple, but slow example of implementing dFBA using COBRApy and scipy.integrate.solve_ivp. This notebook shows a static optimization approach (SOA) implementation and should not be considered production ready.

The model considers only basic Michaelis-Menten limited growth on glucose.

[1]:

import numpy as np
from tqdm import tqdm

from scipy.integrate import solve_ivp

import matplotlib.pyplot as plt
%matplotlib inline


Create or load a cobrapy model. Here, we use the ‘textbook’ e-coli core model.

[2]:

import cobra

Scaling...
A: min|aij| =  1.000e+00  max|aij| =  1.000e+00  ratio =  1.000e+00
Problem data seem to be well scaled


## 15.1. Set up the dynamic system¶

Dynamic flux balance analysis couples a dynamic system in external cellular concentrations to a pseudo-steady state metabolic model.

In this notebook, we define the function add_dynamic_bounds(model, y) to convert the external metabolite concentrations into bounds on the boundary fluxes in the metabolic model.

[3]:

def add_dynamic_bounds(model, y):
"""Use external concentrations to bound the uptake flux of glucose."""
biomass, glucose = y  # expand the boundary species
glucose_max_import = -10 * glucose / (5 + glucose)
model.reactions.EX_glc__D_e.lower_bound = glucose_max_import

def dynamic_system(t, y):
"""Calculate the time derivative of external species."""

biomass, glucose = y  # expand the boundary species

# Calculate the specific exchanges fluxes at the given external concentrations.
with model:

feasibility = cobra.util.fix_objective_as_constraint(model)
model, ['Biomass_Ecoli_core', 'EX_glc__D_e'], ['max', 'max'])

# Since the calculated fluxes are specific rates, we multiply them by the
# biomass concentration to get the bulk exchange rates.
fluxes = lex_constraints.values
fluxes *= biomass

# This implementation is **not** efficient, so I display the current
# simulation time using a progress bar.
if dynamic_system.pbar is not None:
dynamic_system.pbar.update(1)
dynamic_system.pbar.set_description('t = {:.3f}'.format(t))

return fluxes

dynamic_system.pbar = None

def infeasible_event(t, y):
"""
Determine solution feasibility.

Avoiding infeasible solutions is handled by solve_ivp's built-in event detection.
This function re-solves the LP to determine whether or not the solution is feasible
(and if not, how far it is from feasibility). When the sign of this function changes
from -epsilon to positive, we know the solution is no longer feasible.

"""

with model:

feasibility = cobra.util.fix_objective_as_constraint(model)

return feasibility - infeasible_event.epsilon

infeasible_event.epsilon = 1E-6
infeasible_event.direction = 1
infeasible_event.terminal = True


## 15.2. Run the dynamic FBA simulation¶

[4]:

ts = np.linspace(0, 15, 100)  # Desired integration resolution and interval
y0 = [0.1, 10]

with tqdm() as pbar:
dynamic_system.pbar = pbar

sol = solve_ivp(
fun=dynamic_system,
events=[infeasible_event],
t_span=(ts.min(), ts.max()),
y0=y0,
t_eval=ts,
rtol=1e-6,
atol=1e-8,
method='BDF'
)

t = 5.804: : 185it [00:05, 34.89it/s]


Because the culture runs out of glucose, the simulation terminates early. The exact time of this ‘cell death’ is recorded in sol.t_events.

[5]:

sol

[5]:

  message: 'A termination event occurred.'
nfev: 179
njev: 2
nlu: 14
sol: None
status: 1
success: True
t: array([0.        , 0.15151515, 0.3030303 , 0.45454545, 0.60606061,
0.75757576, 0.90909091, 1.06060606, 1.21212121, 1.36363636,
1.51515152, 1.66666667, 1.81818182, 1.96969697, 2.12121212,
2.27272727, 2.42424242, 2.57575758, 2.72727273, 2.87878788,
3.03030303, 3.18181818, 3.33333333, 3.48484848, 3.63636364,
3.78787879, 3.93939394, 4.09090909, 4.24242424, 4.39393939,
4.54545455, 4.6969697 , 4.84848485, 5.        , 5.15151515,
5.3030303 , 5.45454545, 5.60606061, 5.75757576])
t_events: [array([5.80191035])]
y: array([[ 0.1       ,  0.10897602,  0.11871674,  0.12927916,  0.14072254,
0.15310825,  0.16649936,  0.18095988,  0.19655403,  0.21334507,
0.23139394,  0.25075753,  0.27148649,  0.29362257,  0.31719545,
0.34221886,  0.36868605,  0.3965646 ,  0.42579062,  0.4562623 ,
0.48783322,  0.52030582,  0.55342574,  0.58687742,  0.62028461,
0.65321433,  0.685188  ,  0.71570065,  0.74425054,  0.77037369,
0.79368263,  0.81390289,  0.83089676,  0.84467165,  0.85535715,
0.8631722 ,  0.86843813,  0.8715096 ,  0.8727423 ],
[10.        ,  9.8947027 ,  9.78040248,  9.65642157,  9.52205334,
9.37656372,  9.21919615,  9.04917892,  8.86573366,  8.6680879 ,
8.45549026,  8.22722915,  7.98265735,  7.72122137,  7.442497  ,
7.14623236,  6.83239879,  6.50124888,  6.15338213,  5.78981735,
5.41206877,  5.02222068,  4.62299297,  4.21779303,  3.81071525,
3.40650104,  3.01042208,  2.6280723 ,  2.26504645,  1.92656158,
1.61703023,  1.33965598,  1.09616507,  0.88670502,  0.70995892,
0.56344028,  0.44387781,  0.34762375,  0.27100065]])
y_events: [array([[0.87280538, 0.25178492]])]


### 15.2.1. Plot timelines of biomass and glucose¶

[6]:

ax = plt.subplot(111)
ax.plot(sol.t, sol.y.T[:, 0])
ax2 = plt.twinx(ax)
ax2.plot(sol.t, sol.y.T[:, 1], color='r')

ax.set_ylabel('Biomass', color='b')
ax2.set_ylabel('Glucose', color='r')

[6]:

Text(0, 0.5, 'Glucose')

[ ]: