Chemical Reaction Models

The biological models functionality is provided by DiffEqBiological.jl and helps the user 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 modeling package.

The Reaction DSL

The @reaction_network DSL allows you to define reaction networks in a more scientific 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, 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 input rType designates the type of this instance (all instances will inherit from the abstract type AbstractReactionNetwork).

The DSL can handle several types of arrows, in both backwards and forward direction. If a bi-directional arrow is used two reaction rates must be designated. These two reaction networks are identical:

rn1 = @reaction_network rType begin
  2.0, X + Y → XY               
  1.0, XY > Z       
  1.0, X + Y ← XY               
  0.5, XY < Z           
end
rn1 = @reaction_network rType begin
  (2.0,1.0), X + Y ↔ XY               
  (1.0, 0.5), XY ⟷ Z       
end

The empty set can be used for production or degradation and is declared using either 0 or ∅. Integers denote the number of each reactant partaking in the reaction.

rn1 = @reaction_network rType begin
  2.0, 2X --> 0        
  2.0, ∅ --> X  
end

Multiple reactions can be declared in a single line:

rn = @reaction_network rType begin
  2.0, (X,Y) --> 0                   #Identical to reactions [2.0, X --> 0] and [2.0, Y --> 0]
  (2.0, 1.0), (X,Y) --> 0            #Identical to reactions [2.0, X --> 0] and [1.0, X --> 0]
  2.0, (X1,Y1) --> (X2,Y2)           #Identical to reactions [2.0, X1 --> X2] and [2.0, Y1 --> Y2]
  (2.0,1.0), X + Y ↔ XY              #Identical to reactions [2.0, X + Y --> XY] and [1.0, XY --> X + Y].
  ((2.0,1.0),(1.0,2.0)), (X,Y) ↔ 0   #Identical to reactions [(2.0,1.0), X ↔ 0] and [(1.0,2.0), Y ↔ 0].
end

Parameters can be added to the network by declaring them after the reaction network. Parameters can only exist in the reaction rate and not as a part of the reaction.

rn = @reaction_network rType begin
    (kB, kD), X + Y ↔ XY
end kB kD
p = [2.0, 1.0]

The parameter set p must be passed to the problem constructor. The parameter values can be changed after the reaction network is defined.

The reaction rate do not need to be constant, but maybe depend on the concentration of the reactants.

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

The hill function hill(x,v,K,n) = v*(x^n)/(x^n + K^n) can be used, as well as the Michaelis Menten function (the hill function with n = 1).

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

By using the @reaction_func macro it is possible to define your own functions, which may then be used when creating new reaction networks.

@reaction_func hill2(x, v, k) = v*x^2/(k^2+x^2)    
rn = @reaction_network rType begin
  (1.0,hill2(XY,1.5,2.0)), X + Y ↔ XY
end

Reaction rates are automatically adjusted according mass kinetics, including taking special account of higher order terms like 2X -->. This can be disabled using any non-filled arrow (⇐, ⟽, ⇒, ⟾, ⇔, ⟺), in which case the reaction rate will be exactly as input. E.g the two reactions in

rn = @reaction_network rType begin
    2.0, X + Y --> XY
    2.0*X*Y, X + Y ⟾ XY
end

will both have reaction rate equal to 2[X][Y].

Once a reaction network has been created it can be passed as input to either one of the ODEProblem, SDEProblem or JumpProblem constructors.

  probODE = ODEProblem(rn, args...; kwargs...)      
  probSDE = SDEProblem(rn, args...; kwargs...)
  probJump = JumpProblem(prob,aggregator::Direct,rn)

the output problems may then be used as normal input to the solvers of the DifferentialEquations package.

The noise used by the SDEProblem will correspond to the Chemical Langevin Equations. However it is possible to scale the amount of noise be declaring a noise parameter. This will be done after declaring the type but before the network.

rn = @reaction_network BiDir begin
    2.0, X + Y ↔ XY
end
p = [0.5]

The noise term is then added as an additional parameter to the network (by default the last parameter in the parameter array, unless also declared after the reaction network among the other parameters). By reducing (or increasing) the noise term the amount stochastic fluctuations in the system can be reduced (or increased).

It is possible to access expressions corresponding to the functions determining the deterministic and stochastic development of the network using.

  f_expr = rn.f_func
  g_expr = rn.g_func

This can e.g. be used to generate LaTeX code corresponding to the system.

Example: Birth-Death Process

rs = @reaction_network begin
  c1, X --> 2X
  c2, X --> 0
  c3, 0 --> X
end c1 c2 c3
p = (2.0,1.0,0.5)
prob = DiscreteProblem([5], (0.0, 4.0), p)
jump_prob = JumpProblem(prob, Direct(), rs)
sol = solve(jump_prob, Discrete())

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)
# S = 301, E = 100, SE = 0, P = 0
prob = DiscreteProblem([301, 100, 0, 0], (0.0, 100.0), p)
jump_prob = JumpProblem(prob, Direct(), rs)
sol = solve(jump_prob, Discrete())