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.
In [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.
In [2]:
import cobra
from cobra.test import create_test_model
model = create_test_model('textbook')
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.
In [7]:
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:
add_dynamic_bounds(model, y)
cobra.util.add_lp_feasibility(model)
feasibility = cobra.util.fix_objective_as_constraint(model)
lex_constraints = cobra.util.add_lexicographic_constraints(
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:
add_dynamic_bounds(model, y)
cobra.util.add_lp_feasibility(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
In [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'
)
Because the culture runs out of glucose, the simulation terminates early. The exact time of this 'cell death' is recorded in sol.t_events.
In [5]:
sol
Out[5]:
In [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')
Out[6]: