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
@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
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
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
∅. 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
Once a reaction network has been created it can be passed as input to either one of the
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
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(, (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())