Solver options

This tutorial demonstrates how the user can modify the choice of algorithm, solver, and display options to be used for simulating a DFBA model. The example DFBA model is identical to that found in the first tutorial but the direct method is employed for simulation.

from os.path import dirname, join, pardir

from cobra.io import read_sbml_model

from dfba import DfbaModel, ExchangeFlux, KineticVariable

1. Load model and add variables

The DFBA model from the first example is generated.

# DfbaModel instance initialized with cobra model
path_to_model = join(pardir, "sbml-models", "iJR904.xml.gz")
fba_model = read_sbml_model(path_to_model)
fba_model.solver = "glpk"
dfba_model = DfbaModel(fba_model)

# instances of KineticVariable (default initial conditions are 0.0, but can be
# set here if wanted e.g. Oxygen)
X = KineticVariable("Biomass")
Gluc = KineticVariable("Glucose")
Xyl = KineticVariable("Xylose")
Oxy = KineticVariable("Oxygen", initial_condition=0.24)
Eth = KineticVariable("Ethanol")

# add kinetic variables to dfba_model
dfba_model.add_kinetic_variables([X, Gluc, Xyl, Oxy, Eth])

# instances of ExchangeFlux
mu = ExchangeFlux("BiomassEcoli")
v_G = ExchangeFlux("EX_glc(e)")
v_Z = ExchangeFlux("EX_xyl_D(e)")
v_O = ExchangeFlux("EX_o2(e)")
v_E = ExchangeFlux("EX_etoh(e)")

# add exchange fluxes to dfba_model
dfba_model.add_exchange_fluxes([mu, v_G, v_Z, v_O, v_E])

# add rhs expressions for kinetic variables in dfba_model
dfba_model.add_rhs_expression("Biomass", mu * X)
dfba_model.add_rhs_expression("Glucose", v_G * 180.1559 * X / 1000.0)
dfba_model.add_rhs_expression("Xylose", v_Z * 150.13 * X / 1000.0)
dfba_model.add_rhs_expression("Oxygen", v_O * 16.0 * X / 1000.0)
dfba_model.add_rhs_expression("Ethanol", v_E * 46.06844 * X / 1000.0)

# add lower/upper bound expressions for exchange fluxes in dfba_model together
# with expression that must be non-negative for correct evaluation of bounds
dfba_model.add_exchange_flux_lb(
    "EX_glc(e)", 10.5 * (Gluc / (0.0027 + Gluc)) * (1 / (1 + Eth / 20.0)), Gluc
)
dfba_model.add_exchange_flux_lb("EX_o2(e)", 15.0 * (Oxy / (0.024 + Oxy)), Oxy)
dfba_model.add_exchange_flux_lb(
    "EX_xyl_D(e)",
    6.0
    * (Xyl / (0.0165 + Xyl))
    * (1 / (1 + Eth / 20.0))
    * (1 / (1 + Gluc / 0.005)),
    Xyl,
)

# add initial conditions for kinetic variables in dfba_model biomass (gDW/L),
# metabolites (g/L)
dfba_model.add_initial_conditions(
    {
        "Biomass": 0.03,
        "Glucose": 15.5,
        "Xylose": 8.0,
        "Oxygen": 0.0,
        "Ethanol": 0.0,
    }
)

2. Choose solver options

By default, the simulate method implements Algorithm 1 as described in Harwood et al., 2016. Here, this default option is overridden using the set_algorithm method of the solver_data attribute to implement the direct method.

The default display is full, meaning that both results from SUNDIALS and GLPK iterations will be printed to screen. Here, this default option is overridden using the set_display method of the solver_data attribute to set display to none. Alternatively, glpk_only or sundials_only could be selected.

The last command uses the set_ode_method method of the solver_data attribute to select Adams’ method for solving the ODEs. By default, simulation of a DFBA model via the direct method will use BDF.

# set solver options for direct method
dfba_model.solver_data.set_algorithm("direct")
dfba_model.solver_data.set_display("none")
dfba_model.solver_data.set_ode_method("ADAMS")

The model is simulated using the simulate method. This simulation covers the interval \([0.0,25.0]\) hours and returns trajectories of kinetic variables and fluxes as separate pandas.DataFrames. Unlike in the first tutorial however, where the third argument (here \(0.001\) hours) simply specified the frequency of storing results, when simulating using the direct method this argument is also used to specify the frequency of calls the the GLPK linear solver. The value of this argument is model dependent (too small means many calls to the GLKP solver that can slow the simulation, while too large means ODE integration will likely fail) and sometimes the direct method can produce spurious numerical results. Indeed, this was the motivation for developing Algorithm 1 in Harwood et al., 2016 and so users are encouraged to employ the default algorithm whenever possible.

concentrations, trajectories = dfba_model.simulate(
    0.0, 25.0, 0.001, ["EX_glc(e)", "EX_xyl_D(e)", "EX_etoh(e)"]
)

3. Plot the results

from dfba.plot.plotly import *

import plotly.io as pio

# in plotly version 4, default version is different than in 3
pio.templates.default = "plotly_white"

fig = plot_concentrations(concentrations)
fig.show()

fig = plot_trajectories(trajectories)
fig.show()