Tutorial

This overarching tutorial of SBMLImporter will cover how to import an SBML model into a JumpProblem for Gillespie simulations, a SDEProblem for chemical Langevin simulations, or an ODEProblem for deterministic simulations.

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.

load_SBML returns two outputs: a ParsedReactionSystem (prn) and a CallbackSet (cb). The ParsedReactionSystem includes a Catalyst ReactionSystem (prn.rn), a map for the initial condition values of each species (prn.u0), and a map for setting the model parameter values (prn.p). The CallbackSet holds any potential SBML events, along with SBML piecewise functions that have been parsed into events.

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; lw = 3, xlabel = "Time [s]", ylabel = "Protein amount")

For more information on jump simulations, see JumpProcesses.jl's documentation.

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)

WARNING: could not import OrdinaryDiffEqCore.has_lazy_interpolation into OrdinaryDiffEqLowOrderRK

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; lw = 3, xlabel = "Time [s]", ylabel = "Protein amount")

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

ODE simulations

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

using ModelingToolkit, OrdinaryDiffEq
tspan = (0.0, 10.0)
sys = structural_simplify(convert(ODESystem, prn.rn))
oprob = ODEProblem(sys, prn.u0, tspan, prn.p, jac=true)

Here we call structural_simplify to inline any potential assignment rules into the model equations. Additionally, 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; lw = 3, xlabel = "Time [s]", ylabel = "Protein amount")

For more information on ODE simulations, see OrdinaryDiffEq.jl's documentation.

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.