DiffEqBiological.jl
DiffEqBiological.jl provides a domain specific language (DSL) for defining chemical reaction networks in Julia. It interfaces with the broader DifferentialEquations.jl infrastructure to enable the easy generation and solution of corresponding mass action ODE models, Chemical Langevin SDE models, and stochastic chemical kinetics jump models. These generated models can also be used in higher level DifferentialEquations.jl packages (e.g. for sensitivity analysis, parameter estimation, etc).
Here we give a brief introduction to using the DiffEqBiological package, with a focus on how to define reaction networks, and a minimal example showing how to create and solve ODE, SDE and jump models.
More detailed documentation is available from:
- Several DiffEqBiological tutorials are available as part of the
DiffEqTutorials Modeling
Examples. Both html and
interactive IJulia notebook versions are provided there. These include
- An introductory tutorial showing how to specify and solve both ODE and stochastic versions of the repressilator.
- A tutorial exploring the DiffEqBiological API for querying network properties, which also illustrates how to programmatically and incrementally construct and solve a network model using the API.
- A tutorial showing how to use the wrapped HomotopyContinuation.jl functionality to find steady-states and make bifurcation plots.
- Full documentation of the DSL syntax, with information on the generated rate functions and models is available in the DifferentialEquations.jl Chemical Reaction Models documentation.
- API documentation showing how to retrieve network information from a
generated
reaction_network
is available here.
The Reaction DSL
The @reaction_network
DSL allows for the definition of reaction networks using
a simple format. Its input is a set of chemical reactions, from which it
generates a reaction network object which can be used as input to ODEProblem
,
SteadyStateProblem
, SDEProblem
and JumpProblem
constructors.
The basic syntax is
rn = @reaction_network rType begin
2.0, X + Y --> XY
1.0, XY --> Z
end
where each line corresponds to a chemical reaction. The (optional) input rType
designates the type of this instance (all instances will inherit from the
abstract type AbstractReactionNetwork
).
The DSL has many features:
- It supports many different arrow types, corresponding to different directions
of reactions and different rate laws:
rn = @reaction_network begin 1.0, X + Y --> XY 1.0, X + Y → XY 1.0, XY ← X + Y 2.0, X + Y ↔ XY end
- It allows multiple reactions to be defined simultaneously on one line. The
following two networks are equivalent:
rn1 = @reaction_network begin (1.0,2.0), (S1,S2) → P end rn2 = @reaction_network begin 1.0, S1 → P 2.0, S2 → P end
- It allows the use of symbols to represent reaction rate parameters, with their
numeric values specified during problem construction. i.e., the previous
example could be given by
withrn2 = @reaction_network begin k1, S1 → P k2, S2 → P end k1 k2
k1
andk2
corresponding to the reaction rates. - Rate law functions can be pre-defined and used within the DSL:
@reaction_func myHill(x) = 2.0*x^3/(x^3+1.5^3) rn = @reaction_network begin myHill(X), ∅ → X end
- Pre-made rate laws for Hill and Michaelis-Menten reactions are provided:
rn1 = @reaction_network begin hill(X,v,K,n), ∅ → X mm(X,v,K), ∅ → X end v K n
- Simple rate law functions of the species populations can be used within the DSL:
rn = @reaction_network begin 2.0*X^2, 0 → X + Y gamma(Y)/5, X → ∅ pi*X/Y, Y → ∅ end
For sufficiently large and structured network models it can often be easier to
specify some reactions through a programmatic API. For this reason the
@min_reaction_network
and @empty_reaction_network
macros, along with the
corresponding addspecies!
, addparam!
and addreaction!
modifier functions,
are provided in the
API.
DiffEqBiological API for Querying Network Information
A variety of network information is calculated by the reaction_network
macro,
and can then be retrieved using the DiffEqBiological
API. This includes
- Orderings of species and reactions
speciesmap(rn) paramsmap(rn)
- Reaction stoichiometries
substratestoich(rn, rxidx) productstoich(rn, rxidx) netstoich(rn, rxidx)
- Expressions corresponding to the functions determining the deterministic and
stochastic terms within resulting ODE, SDE or jump models
These can be used to generate LaTeX expressions corresponding to the system using packages such asode_exprs = odeexprs(rn) jacobian_exprs = jacobianexprs(rn) noise_expr = noiseexprs(rn) rate_exprs,affect_exprs = jumpexprs(rn)
Latexify
. - Dependency graphs
rxtospecies_depgraph(rn) speciestorx_depgraph(rn) rxtorx_depgraph(rn)
and more.
Simulating ODE, Steady-State, SDE and Jump Problems
Once a reaction network has been created it can be passed as input to either one
of the ODEProblem
, SteadyStateProblem
, SDEProblem
or JumpProblem
constructors.
probODE = ODEProblem(rn, args...; kwargs...)
probSS = SteadyStateProblem(rn, args...; kwargs...)
probSDE = SDEProblem(rn, args...; kwargs...)
probJump = JumpProblem(prob, Direct(), rn)
The output problems may then be used as input to the solvers of
DifferentialEquations.jl. Note, the noise used by
the SDEProblem
will correspond to the Chemical Langevin Equations.
As an example, consider models for a simple birth-death process:
rs = @reaction_network begin
c1, X --> 2X
c2, X --> 0
c3, 0 --> X
end c1 c2 c3
params = (1.0,2.0,50.)
tspan = (0.,4.)
u0 = [5.]
# solve ODEs
oprob = ODEProblem(rs, u0, tspan, params)
osol = solve(oprob, Tsit5())
# solve for Steady-States
ssprob = SteadyStateProblem(rs, u0, params)
sssol = solve(ssprob, SSRootfind())
# solve Chemical Langevin SDEs
sprob = SDEProblem(rs, u0, tspan, params)
ssol = solve(sprob, EM(), dt=.01)
# solve JumpProblem using Gillespie's Direct Method
u0 = [5]
dprob = DiscreteProblem(rs, u0, tspan, params)
jprob = JumpProblem(dprob, Direct(), rs)
jsol = solve(jprob, SSAStepper())
Importing Predefined Networks
ReactionNetworkImporters.jl
can load several different types of predefined networks into DiffEqBiological
reaction_network
s. These include
- A subset of BioNetGen .net files that can be generated from a BioNetGen language file (.bngl). (.net files can be generated using the
generate_network
command within BioNetGen.) - Reaction networks specified by dense or sparse matrices encoding the stoichiometry of substrates and products within each reaction.
- Networks defined by the basic file format used by the RSSA group at COSBI in their model collection.
Finding steady states
The steady states of a reaction network can be found using homotopy continuation (as implemented by HomotopyContinuation.jl). This method is limited to polynomial systems, which includes reaction network not containing non-polynomial rates in the reaction rates (such as logarithms and non integer exponents). Note, both the steady-state and the bifurcation diagram functionality only fully support Julia 1.1 and greater.
The basic syntax is
rn = @reaction_network begin
(1.,2.), 0 ↔ X
end
ss = steady_states(rn,params)
and with parameters
rn = @reaction_network begin
(p,d), 0 ↔ X
end p d
params = [1., 2.]
ss = steady_states(rn,params)
the stability of a steady state (or a vector of several) can be determined by the stability
function:
stability(ss,rn,params)
Here the @reaction_network
creates a multivariate polynomial and stores in the equilibrate_content
field in the reaction network structure. The steady_state
method the inserts the corresponding parameter values and solves the polynomial system. The exception is if there exists a parameter as an exponent (typically n
in a hill function). In this case the steady state polynomial can first be created in the steady_state
method. If one plans to solve a polynomial a large number of times with the same value of n
, then one can get a speed-up by first fixing that value using the fix_parameters
function:
rn = @reaction_network begin
(hill(X,v,K,n),d), 0 ↔ X
end v K n d
fix_parameters(rn,n=4)
for i = 1:10000
params = [i, 2.5, 4, 0.1] #The value of 'n' here doesn't really matter, however, the field must exist.
ss = steady_states(rn,params)
Some networks may have an infinite set of steady states, and which one is interested in depends on the initial conditions. For these networks some additional information is required (typically some concentrations which sums to a fixed value). This information can be added through the @add_constraint
macro:
rn = @reaction_network begin
(k1,k2), X ↔ Y
end k1 k2
params = [2.,1.]
@add_constraint rn X+Y=2.
steady_states(rn,params)
The @add_constraint
macro may contain parameters, as long as these are declared in the network.
rn = @reaction_network begin
(k1,k2), X ↔ Y
end k1 k2 C_tot
params = [2.,1.,2.]
@add_constraint rn X+Y=C_tot
steady_states(rn,params)
The @add_constraints
macro can be used to add several constraints at the same time.
rn = @reaction_network begin
(k1,k2), X ↔ Y
(k3,k4), V ↔ W
end k1 k2 k3 k4
params = [2.,1.,1.,2.]
@add_constraints rn begin
X + Y = 2.
V + W = 4.
end
steady_states(rn,params)
Making bifurcation diagram
For any system for which we can find steady states, we can also make bifurcation diagrams.
rn = @reaction_network begin
(p,d), 0 ↔ X
end p d
params = [1.,2.] #The value of 'p' here doesn't really matter, however, the field must exist.
bif = bifurcations(rn, params, :p, (0.1,5.))
These can then be plotted.
plot(bif)
In the plot blue values correspond to stable steady states, red to unstable. Also, cyan correspond to stable steady states with imaginary eigen values and orange to unstable steady states with imaginary eigen values.
In addition to the normal bifurcation diagram (varying a single parameter over a continuous range) there are three more types available.
A bifurcation grid varies a single parameter over a set of discrete values
bif_grid = bifurcation_grid(rn, params, :p, 1.:5.)
A two dimensional bifurcation grid varies two different parameters over a grid of discrete values.
bif_grid_2d = bifurcation_grid_2d(rn, params, :p, 1.:5., :d, 2.:10.)
A bifurcation diagram grid first varies a single variable over a discrete grid of values. Then, for each such value, in varies a second variable over a continuous interval to create a bifurcation grid.
bif_grid_dia = bifurcation_grid_diagram(rn, params, :p, 1.:5., :d, (2.,10.))
All of these can be plotted.
plot(bif_grid)
plot(bif_grid_2d)
plot(bif_grid_dia)