Chemical Reaction Models

The biological models functionality is provided by DiffEqBiological.jl and helps the user to build discrete stochastic and differential equation based systems biological models. These tools allow one to define the models at a high level by specifying reactions and rate constants, and the creation of the actual problems is then handled by the modelling package.

Installation

This functionality does not come standard with DifferentialEquations.jl. To use this functionality, you must install DiffEqBiological.jl:

]add DiffEqBiological
using DiffEqBiological

The Reaction DSL - Basic

This section covers some of the basic syntax for building chemical reaction network models. Examples showing how to both construct and solve network models are provided in Chemical Reaction Network Examples.

Basic syntax

The @reaction_network macro allows the (symbolic) specification of reaction networks with a simple format. Its input is a set of chemical reactions, and from them it generates a reaction network object which can be used as input to ODEProblem, SteadyStateProblem, SDEProblem and JumpProblem constructors. The @min_reaction_network macro constructs a more simplified reaction network, deferring construction of all the various functions needed for each of these problem types. It can then be incrementally filled in for specific problem types as needed, which reduces network construction time for very large networks (see The Min Reaction Network Object for a detailed description).

The basic syntax is:

rn = @reaction_network begin
  2.0, X + Y --> XY               
  1.0, XY --> Z1 + Z2            
end

where each line corresponds to a chemical reaction. Each reaction consists of a reaction rate (the expression on the left hand side of ,), a set of substrates (the expression in-between , and -->), and a set of products (the expression on the right hand side of -->). The substrates and the products may contain one or more reactants, separated by +. The naming convention for these are the same as for normal variables in Julia.

The chemical reaction model is generated by the @reaction_network macro and stored in the rn variable (a normal variable, do not need to be called rn). The macro generates a differential equation model according to the law of mass action, in the above example the ODEs become:

\[\frac{d[X]}{dt} = -2 [X] [Y]\\ \frac{d[Y]}{dt} = -2 [X] [Y]\\ \frac{d[XY]}{dt} = 2 [X] [Y] - [XY]\\ \frac{d[Z1]}{dt}= [XY]\\ \frac{d[Z2]}{dt} = [XY]\]

Arrow variants

Several types of arrows are accepted by the DSL and works instead of -->. All of these works: >, → ↣, ↦, ⇾, ⟶, ⟼, ⥟, ⥟, ⇀, ⇁. Backwards arrows can also be used to write the reaction in the opposite direction. Hence all of these three reactions are equivalent:

rn = @reaction_network begin
  1.0, X + Y --> XY               
  1.0, X + Y → XY      
  1.0, XY ← X + Y      
end

(note that due to technical reasons <-- cannot be used)

Using bi-directional arrows

Bi-directional arrows can be used to designate a reaction that goes two ways. These two models are equivalent:

rn = @reaction_network begin
  2.0, X + Y → XY             
  2.0, X + Y ← XY          
end
rn = @reaction_network begin
  2.0, X + Y ↔ XY               
end

If the reaction rate in the backwards and forwards directions are different they can be designated in the following way:

rn = @reaction_network begin
  (2.0,1.0) X + Y ↔ XY               
end

which is identical to

rn = @reaction_network begin
  2.0, X + Y → XY             
  1.0, X + Y ← XY          
end

Combining several reactions in one line

Several similar reactions can be combined in one line by providing a tuple of reaction rates and/or substrates and/or products. If several tuples are provided they much all be of identical length. These pairs of reaction networks are all identical:

rn1 = @reaction_network begin
  1.0, S → (P1,P2)               
end
rn2 = @reaction_network begin
  1.0, S → P1     
  1.0, S → P2
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

rn1 = @reaction_network begin
  (1.0,2.0,3.0), (S1,S2,S3) → (P1,P2,P3)        
end
rn2 = @reaction_network begin
  1.0, S1 → P1
  2.0, S2 → P2   
  3.0, S3 → P3  
end

This can also be combined with bi-directional arrows in which case separate tuples can be provided for the backward and forward reaction rates separately. These reaction networks are identical

rn1 = @reaction_network begin
 (1.0,(1.0,2.0)), S ↔ (P1,P2)  
end
rn2 = @reaction_network begin
  1.0, S → P1
  1.0, S → P2
  1.0, P1 → S   
  2.0, P2 → S
end

Production and Destruction and Stoichiometry

Sometimes reactants are produced/destroyed from/to nothing. This can be designated using either 0 or ∅:

rn = @reaction_network begin
  2.0, 0 → X
  1.0, X → ∅
end

Sometimes several molecules of the same reactant is involved in a reaction, the stoichiometry of a reactant in a reaction can be set using a number. Here two species of X forms the dimer X2:

rn = @reaction_network begin
  1.0, 2X → X2
end

this corresponds to the differential equation:

\[\frac{d[X]}{dt} = -[X]^2\\ \frac{d[X2]}{dt} = \frac{1}{2!} [X]^2\]

Other numbers than 2 can be used and parenthesises can be used to use the same stoichiometry for several reactants:

rn = @reaction_network begin
  1.0, X + 2(Y + Z) → XY2Z2
end

Variable reaction rates

Reaction rates do not need to be constant, but can also depend on the current concentration of the various reactants (when e.g. one reactant activate the production of another one). E.g. this is a valid notation:

rn = @reaction_network begin
  X, Y → ∅
end

and will have Y degraded at rate

\[\frac{d[Y]}{dt} = -[X][Y]\]

Note that this is actually equivalent to the reaction

rn = @reaction_network begin
  1.0, X + Y → X
end

Most expressions and functions are valid reaction rates, e.g:

rn = @reaction_network begin
  2.0*X^2, 0 → X + Y
  gamma(Y)/5, X → ∅
  pi*X/Y, Y → ∅
end

please note that user defined functions cannot be used directly (see later section "User defined functions in reaction rates").

Defining parameters

Just as when defining normal differential equations using DifferentialEquations parameter values does not need to be set when the model is created. Components can be designated as parameters by declaring them at the end:

rn = @reaction_network begin
  p, ∅ → X
  d, X → ∅
end p d

Parameters can only exist in the reaction rates (where they can be mixed with reactants). All variables not declared at the end will be considered a reactant.

Pre-defined functions

Hill functions and a Michaelis-Menten function are pre-defined and can be used as rate laws. Below, the pair of reactions within rn1 are equivalent, as are the pair of reactions within rn2:

rn1 = @reaction_network begin
  hill(X,v,K,n), ∅ → X
  v*X^n/(X^n+K^n), ∅ → X
end v K n
rn2 = @reaction_network begin
  mm(X,v,K), ∅ → X
  v*X/(X+K), ∅ → X
end v K

Repressor Hill (hillr) and Michaelis-Menten (mmr) functions are also provided:

rn1 = @reaction_network begin
  hillr(X,v,K,n), ∅ → X
  v*K^n/(X^n+K^n), ∅ → X
end v K n
rn2 = @reaction_network begin
  mmr(X,v,K), ∅ → X
  v*K/(X+K), ∅ → X
end v K

Model Simulation

Once created, a reaction network can be used as input to various problem types which can be solved by DifferentialEquations.jl.

Deterministic simulations using ODEs

A reaction network can be used as input to an ODEProblem instead of a function, using probODE = ODEProblem(rn, args...; kwargs...) E.g. a model can be created and simulated using:

rn = @reaction_network begin
  p, ∅ → X
  d, X → ∅
end p d
p = [1.0,2.0]
u0 = [0.1]
tspan = (0.,1.)
prob = ODEProblem(rn,u0,tspan,p)
sol = solve(prob)

(if no parameters are given p does not need to be provided)

To solve for a steady-state starting from the guess u0, one can use

prob = SteadyStateProblem(rn,u0,p)
sol = solve(prob, SSRootfind())

prob = SteadyStateProblem(rn,u0,p)
sol = solve(prob, DynamicSS(Tsit5()))

Stochastic simulations using SDEs

In a similar way a SDE can be created using probSDE = SDEProblem(rn, args...; kwargs...). In this case the chemical Langevin equations (as derived in Gillespie 2000) will be used to generate stochastic differential equations.

Stochastic simulations using discrete stochastic simulation algorithms

Instead of solving SDEs one can create a stochastic jump process model using integer copy numbers and a discrete stochastic simulation algorithm. This can be done using:

rn = @reaction_network begin
  p, ∅ → X
  d, X → ∅
end p d
p = [1.0,2.0]
u0 = [10]
tspan = (0.,1.)
discrete_prob = DiscreteProblem(rn, u0,tspan,p)
jump_prob = JumpProblem(discrete_prob,Direct(),rn)
sol = solve(jump_prob,SSAStepper())

Here we used Gillespie's Direct method as the underlying stochastic simulation algorithm.

Reaction rate laws used in simulations

In generating mathematical models from a reaction_network, reaction rates are treated as microscopic rates. That is, for a general mass action reaction of the form $n_1 S_1 + n_2 S_2 + \dots n_M S_M \to \dots$ with stoichiometric substrate coefficients $\{n_i\}_{i=1}^M$ and rate constant $k$, the corresponding ODE rate law is taken to be

\[k \prod_{i=1}^M \frac{(S_i)^{n_i}}{n_i!},\]

while the jump process transition rate (i.e. propensity function) is

\[k \prod_{i=1}^M \frac{S_i (S_i-1) \dots (S_i-n_i+1)}{n_i!}.\]

For example, the ODE model of the reaction $2X + 3Y \to Z$ with rate constant $k$ would be

\[\frac{dX}{dt} = -2 k \frac{X^2}{2!} \frac{Y^3}{3!} = -k \frac{X^2 Y^3}{3!} \\ \frac{dY}{dt} = -3 k \frac{X^2}{2!} \frac{Y^3}{3!} = -k \frac{X^2 Y^3}{4} \\ \frac{dZ}{dt} = k \frac{X^2}{2!} \frac{Y^3}{3!}.\]

The Reaction DSL - Advanced

This section covers some of the more advanced syntax for building chemical reaction network models (still not very complicated!).

User defined functions in reaction rates

The reaction network DSL cannot "see" user defined functions. E.g. this is not correct syntax:

myHill(x) = 2.0*x^3/(x^3+1.5^3)
rn = @reaction_network begin
  myHill(X), ∅ → X
end

However, it is possible to define functions in such a way that the DSL can see them using the @reaction_func macro:

@reaction_func myHill(x) = 2.0*x^3/(x^3+1.5^3)
rn = @reaction_network begin
  myHill(X), ∅ → X
end

Defining a custom reaction network type

While the default type of a reaction network is reaction_network (which inherits from AbstractReactionNetwork) it is possible to define a custom type (which also will inherit from AbstractReactionNetwork) by adding the type name as a first argument to the @reaction_network macro:

rn = @reaction_network my_custom_type begin
  1.0, ∅ → X
end

Scaling noise in the chemical Langevin equations

When making stochastic simulations using SDEs it is possible to scale the amount of noise in the simulations by declaring a noise scaling parameter. This parameter is declared as a second argument to the @reaction_network macro (when scaling the noise one have to declare a custom type).

rn = @reaction_network my_custom_type ns begin
  1.0, ∅ → X
end

The noise scaling parameter is automatically added as a last argument to the parameter array (even if not declared at the end). E.g. this is correct syntax:

rn = @reaction_network my_custom_type ns begin
  1.0, ∅ → X
end
p = [0.1,]
u0 = [0.1]
tspan = (0.,1.)
prob = SDEProblem(rn,u0,tspan,p)
sol = solve(prob)

Here the amount of noise in the stochastic simulation will be reduced by a factor 10.

Ignoring mass kinetics

While one in almost all cases want the reaction rate to take the law of mass action into account, so the reaction

rn = @reaction_network my_custom_type ns begin
  k, X → ∅
end k

occur at the rate $d[X]/dt = -k[X]$, it is possible to ignore this by using any of the following non-filled arrows when declaring the reaction: ⇐, ⟽, ⇒, ⟾, ⇔, ⟺. This means that the reaction

rn = @reaction_network my_custom_type ns begin
  k, X ⇒ ∅
end k

will occur at rate $d[X]/dt = -k$ (which might become a problem since $[X]$ will be degraded at a constant rate even when very small or equal to 0.

The Reaction Network Object

The @reaction_network macro generates a reaction_network object, which has a number of fields which can be accessed.

rn.f is a function encoding the right hand side of the ODEs (i.e. the time derivatives of the chemical species).
rn.f_func is a vector of expressions corresponding to the time derivatives of the chemical species.
rn.f_symfuncs is a vector of SymEngine expressions corresponding to the time derivatives of the chemical species.
rn.g is a function encoding the noise terms for the SDEs (see rn.g_func for details).
rn.g_func is a vector containing expressions corresponding to the noise terms used when creating the SDEs (n*m elements when there are n reactants and m reactions. The first m elements correspond to the noise terms for the first reactant and each reaction, the next m elements for the second reactant and all reactions, and so on).
rn.jac is a function that evaluates the Jacobian of rn.f in place. i.e. has the form rn.jac(dJ,u,p,t), for pre-allocated Jacobian matrix dJ.
rn.jump_affect_expr is a vector of expressions for how each reaction causes the species populations to change.
rn.jump_rate_expr is a vector of expressions for how the transition rate (i.e. propensity) of each reaction is calculated from the species populations.
rn.jumps is a vector storing a jump corresponding to each reaction (i.e. ConstantRateJump, VariableRateJump, etc...)
rn.odefun stores an ODEFunction that can be used to create an ODEProblem corresponding to the reaction network.
rn.p_matrix is a prototype matrix with the same size as the noise term.
rn.paramjac is a function that evaluates the Jacobian of rn.f with respect to the parameters, p, in-place. It has the form rn.paramjac(dpJ,u,p,t) for pre-allocated parameter Jacobian matrix dpJ.
rn.params is a vector containing symbols corresponding to all the parameters of the network.
rn.params_to_ints provides a mapping from parameter symbol to the integer id of the parameter (i.e. where it is stored in the parameter vector passed to ODEProblem, SDEProblem, etc...)
rn.reactions stores a vector of DiffEqBiological.ReactionStructs, which collect info for their corresponding reaction (such as stoichiometric coefficients).
rn.regular_jumps stores a RegularJump representation of the network, for use in $\tau$-leaping methods.
rn.scale_noise is the noise scaling parameter symbol (if provided).
rn.sdefun is a SDEFunction that can be used to create an SDEProblem corresponding to the reaction network.
rn.symjac is the symbolically calculated Jacobian of the ODEs corresponding to the model.
rn.syms is a vector containing symbols for all species of the network.
rn.syms_to_ints is a map from the symbol of a species to its integer index within the solution vector.

The Min Reaction Network Object

The @min_reaction_network macro works similarly to the @reaction_network macro, but initially only fills in fields corresponding to basic reaction network properties (i.e. rn.params, rn.params_to_ints, rn.scale_noise, rn.reactions, rn.syms, and rn.syms_to_ints). To fill in the remaining fields call (in the following [val] denotes the default value of a keyword argument):

addodes!(rn) to complete ODE-related fields, optional keyword arguments include:
- build_jac=[true], is true if rn.jac and rn.symjac should be constructed. (Currently these build a dense Jacobian, so should be set to false for sufficiently large systems.)
- build_symfuncs=[true], is true if symbolic functions should be constructed for each ODE rhs. It is recommended to disable this for larger systems to reduce memory usage and speedup network construction.
addsdes!(rn) to complete SDE-related fields.
addjumps!(rn) to complete jump-related fields. addjumps! accepts several keyword arguments to control which jumps get created:
- build_jumps=[true] is true if rn.jumps should be constructed. This can be set to false for regular jump problems, where only rn.regular_jumps is needed.
- build_regular_jumps=[true] is true if rn.regular_jumps should be constructed. This can be set to false for Gillespie-type jump problems, where regular_jumps are not used.
- minimal_jumps=[false] is false if rn.jumps should contain a jump for each possible reaction. If set to true jumps are only added to rn.jumps for non-mass action jumps. (Note, mass action jumps are still resolved within any jump simulation. This option simply speeds up the construction of the jump problem since entries in rn.jumps that correspond to mass action jumps are never directly called within jump simulations.)

For example, to simulate a jump process (i.e. Gillespie) simulation without constructing any RegularJumps, and only constructing a minimal set of jumps:

rs = @min_reaction_network begin
  c1, X --> 2X
  c2, X --> 0
  c3, 0 --> X
end c1 c2 c3
p = (2.0,1.0,0.5)
addjumps!(rs; build_regular_jumps=false, minimal_jumps=true)
prob = DiscreteProblem(rs, [5], (0.0, 4.0), p)
jump_prob = JumpProblem(prob, Direct(), rs)
sol = solve(jump_prob, SSAStepper())

Chemical Reaction Network Examples

Example: Birth-Death Process

rs = @reaction_network begin
  c1, X --> 2X
  c2, X --> 0
  c3, 0 --> X
end c1 c2 c3
p = (1.0,2.0,50.)
tspan = (0.,4.)
u0 = [5.]

# solve ODEs
oprob = ODEProblem(rs, u0, tspan, p)
osol  = solve(oprob, Tsit5())

# solve for Steady-States
ssprob = SteadyStateProblem(rs, u0, p)
sssol  = solve(ssprob, SSRootfind())

# solve SDEs
sprob = SDEProblem(rs, u0, tspan, p)
ssol  = solve(sprob, EM(), dt=.01)

# solve JumpProblem
u0 = [5]
dprob = DiscreteProblem(rs, u0, tspan, p)
jprob = JumpProblem(dprob, Direct(), rs)
jsol = solve(jprob, SSAStepper())

Example: Michaelis-Menten Enzyme Kinetics

rs = @reaction_network begin
  c1, S + E --> SE
  c2, SE --> S + E
  c3, SE --> P + E
end c1 c2 c3
p = (0.00166,0.0001,0.1)
tspan = (0., 100.)
u0 = [301., 100., 0., 0.]  # S = 301, E = 100, SE = 0, P = 0

# solve ODEs
oprob = ODEProblem(rs, u0, tspan, p)
osol  = solve(oprob, Tsit5())

# solve JumpProblem
u0 = [301, 100, 0, 0] 
dprob = DiscreteProblem(rs, u0, tspan, p)
jprob = JumpProblem(dprob, Direct(), rs)
jsol = solve(jprob, SSAStepper())