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:
reaction_networkis available here.
@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
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
designates the type of this instance (all instances will inherit from the
The DSL has many features:
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
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
rn2 = @reaction_network begin k1, S1 → P k2, S2 → P end k1 k2with
k2corresponding to the reaction rates.
@reaction_func myHill(x) = 2.0*x^3/(x^3+1.5^3) rn = @reaction_network begin myHill(X), ∅ → X end
rn1 = @reaction_network begin hill(X,v,K,n), ∅ → X mm(X,v,K), ∅ → X end v K n
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
@empty_reaction_network macros, along with the
addreaction! modifier functions,
are provided in the
A variety of network information is calculated by the
and can then be retrieved using the DiffEqBiological
API. This includes
substratestoich(rn, rxidx) productstoich(rn, rxidx) netstoich(rn, rxidx)
ode_exprs = odeexprs(rn) jacobian_exprs = jacobianexprs(rn) noise_expr = noiseexprs(rn) rate_exprs,affect_exprs = jumpexprs(rn)These can be used to generate LaTeX expressions corresponding to the system using packages such as
rxtospecies_depgraph(rn) speciestorx_depgraph(rn) rxtorx_depgraph(rn)and more.
Once a reaction network has been created it can be passed as input to either one
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
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 =  dprob = DiscreteProblem(rs, u0, tspan, params) jprob = JumpProblem(dprob, Direct(), rs) jsol = solve(jprob, SSAStepper())
can load several different types of predefined networks into DiffEqBiological
reaction_networks. These include
generate_networkcommand within BioNetGen.)
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
@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
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
rn = @reaction_network begin (k1,k2), X ↔ Y end k1 k2 params = [2.,1.] @add_constraint rn X+Y=2. steady_states(rn,params)
@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)
@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)
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.
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)
8 days ago