Tutorial

This overarching tutorial of SBMLImporter covers how to import an SBML model into a JumpProblem for Gillespie simulations, a SDEProblem for chemical Langevin simulations, or an ODEProblem for deterministic simulations. It also covers how to change parameter and/or initial values for an imported model.

Input - a valid SBML model file

SBMLImporter only requires one input: a valid SBML file. SBML files can be created from various sources, such as the graphical interface of COPASI, by converting rule-based BioNetGen models, or by using any other SBML exporting tools. Additionally, a large collection of published SBML models are hosted on BioModels.

In this tutorial, we will use the Brusselator model, whose SBML file can be downloaded from here.

Importing a model

To import an SBML model use load_SBML:

using SBMLImporter
prn, cb = load_SBML(path; massaction=true)

Here, massaction=true informs the importer that the model follows chemical mass action kinetics. This is required for efficient jump (Gillespie type) simulations. Additional keyword arguments can be found in the API documentation.

load_SBML returns two outputs: a ParsedReactionSystem (prn) and a CallbackSet (cb). The CallbackSet holds any SBML events, as well as any piecewise functions parsed into events. The ParsedReactionSystem includes a Catalyst ReactionSystem (prn.rn), a map for the initial values of each species (prn.u0), and a map for setting the model parameter values (prn.p). Many modeling tasks can be performed with a Catalyst ReactionSystem, for example, the system reactions can be inspected by:

using Catalyst
reactions(prn.rn)

\[ \begin{align*} \varnothing &\xrightleftharpoons[C]{A C} \mathrm{X} \\ \mathrm{Y} + 2 \mathrm{X} &\xrightarrow{2 C} 3 \mathrm{X} \\ \mathrm{X} &\xrightarrow{B C} \mathrm{Y} \end{align*} \]

Jump simulations

To perform jump simulations (e.g. using Gillespie's algorithm), convert the reaction system (prn.rn) into a JumpProblem (which requires first creating an intermediary DiscreteProblem).

using JumpProcesses
tspan = (0.0, 10.0)
dprob = DiscreteProblem(prn.rn, prn.u0, tspan, prn.p)
jprob = JumpProblem(prn.rn, dprob, Direct())

The JumpProblem can be simulated with any solver from the JumpProcesses.jl package, such as the SSAStepper:

using Plots
sol = solve(jprob, SSAStepper(), callback=cb)
plot(sol; xlabel = "Time [s]", ylabel = "Protein amount")

For more information on jump simulations, see JumpProcesses.jl's documentation. A guide to choosing simulation algorithm can be found here.

SDE simulations

To perform chemical Langevin SDE simulations, convert the reaction-system prn.rn into a SDEProblem.

using StochasticDiffEq
tspan = (0.0, 10.0)
sprob = SDEProblem(prn.rn, prn.u0, tspan, prn.p)

The SDEProblem can be simulated with any solver from the StochasticDiffEq.jl package, such as the LambaEM solver:

sol = solve(sprob, LambaEM(), callback=cb)
plot(sol; xlabel = "Time [s]", ylabel = "Protein amount")

For more information on SDE simulations, see StochasticDiffEq.jl's documentation.

ODE simulations

To perform ODE simulations, first convert the reaction system (prn.rn) into an ODESystem.

using ModelingToolkit
sys = structural_simplify(convert(ODESystem, prn.rn))

Here we call structural_simplify to inline any potential assignment rules into the model equations. Given an ODESystem, a good sanity check is to inspect the generated ODEs:

equations(sys)

\[ \begin{align} \frac{\mathrm{d} Y\left( t \right)}{\mathrm{d}t} &= B C X\left( t \right) - \left( X\left( t \right) \right)^{2} C Y\left( t \right) \\ \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} &= A C - C X\left( t \right) - B C X\left( t \right) + \left( X\left( t \right) \right)^{2} C Y\left( t \right) \end{align} \]

Once we have an ODESystem, we can build an ODEProblem:

using OrdinaryDiffEq
tspan = (0.0, 10.0)
oprob = ODEProblem(sys, prn.u0, tspan, prn.p, jac=true)

Setting jac=true means that the ODE Jacobian is computed symbolically, which often improves simulation performance. The ODEProblem can be simulated with any solver from the OrdinaryDiffEq.jl package, such as the Rodas5P solver:

sol = solve(oprob, Rodas5P(), callback=cb)
plot(sol; xlabel = "Time [s]", ylabel = "Protein amount")

For more information on ODE simulations, see OrdinaryDiffEq.jl's documentation. A guide to choosing ODE solver can be found here.

Changing Parameter Values and/or Initial Values

One of the most common modeling operations is to modify model parameter values and/or initial values. Because SBMLImporter imports models as a Catalyst ReactionSystem that supports the SymbolicIndexingInterface, these modifications are straightforward to perform at the ODEProblem, SDEProblem, or JumpProblem level.

Parameter values can be modified using prob.ps[:param_id]. For example, for the Brusselator imported as an ODEProblem above, the value of parameter A can be changed with:

oprob.ps[:A] = 2.0

Initial values can be modified with prob[:specie_id]. For example, to change the initial value of the species X to 2.0, do:

oprob[:X] = 2.0

Next Steps

For additional modeling tasks that can be carried out with a Catalyst ReactionSystem (e.g., bifurcation analysis, parameter estimation, etc.), see the Catalyst documentation. If you encounter any problems with importing an SBML model, first consult the FAQ. Additional details on importer options can be found in the API.