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.
The massaction
keyword argument is only relevant for jump simulations. If you import the model as a SDEProblem
or an ODEProblem
, you can (and should) ignore this keyword argument.
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.
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, 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; 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.