|
9 | 9 | [](http://pkg.julialang.org/?pkg=DiffEqBiological) |
10 | 10 | [](http://pkg.julialang.org/?pkg=DiffEqBiological) |
11 | 11 |
|
12 | | -This repository holds functionality for easily building systems biological models |
13 | | -for the DifferentialEquations ecosystem. This has tools for developing chemical |
14 | | -reaction networks and simulating them with Gillespie, ODE, SDE, DDE, and DAE |
15 | | -algorithms. |
16 | | - |
17 | 12 | Full documentation is in the |
18 | 13 | [DifferentialEquations.jl models documentation](http://docs.juliadiffeq.org/latest/models/biological.html) |
| 14 | + |
| 15 | +##The Reaction DSL |
| 16 | + |
| 17 | +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. |
| 18 | + |
| 19 | +The basic syntax is |
| 20 | +```julia |
| 21 | +rn = @reaction_network rType begin |
| 22 | + 2.0, X + Y --> XY |
| 23 | + 1.0, XY --> Z |
| 24 | +end |
| 25 | +``` |
| 26 | +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`). |
| 27 | + |
| 28 | +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 |
| 29 | +```julia |
| 30 | +rn1 = @reaction_network rType begin |
| 31 | + 2.0, X + Y → XY |
| 32 | + 1.0, XY > Z |
| 33 | + 1.0, X + Y ← XY |
| 34 | + 0.5, XY < Z |
| 35 | +end |
| 36 | +rn1 = @reaction_network rType begin |
| 37 | + (2.0,1.0), X + Y ↔ XY |
| 38 | + (1.0, 0.5), XY ⟷ Z |
| 39 | +end |
| 40 | +``` |
| 41 | +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. |
| 42 | +```julia |
| 43 | +rn1 = @reaction_network rType begin |
| 44 | + 2.0, 2X --> 0 |
| 45 | + 2.0, ∅ --> X |
| 46 | +end |
| 47 | +``` |
| 48 | +Multiple reactions can be declared in a single line |
| 49 | +```julia |
| 50 | +rn = @reaction_network rType begin |
| 51 | + 2.0, (X,Y) --> 0 #Identical to reactions [2.0, X --> 0] and [2.0, Y --> 0] |
| 52 | + (2.0, 1.0), (X,Y) --> 0 #Identical to reactions [2.0, X --> 0] and [1.0, X --> 0] |
| 53 | + 2.0, (X1,Y1) --> (X2,Y2) #Identical to reactions [2.0, X1 --> X2] and [2.0, Y1 --> Y2] |
| 54 | + (2.0,1.0), X + Y ↔ XY #Identical to reactions [2.0, X + Y --> XY] and [1.0, XY --> X + Y]. |
| 55 | + ((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]. |
| 56 | +end |
| 57 | + ``` |
| 58 | +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. |
| 59 | +```julia |
| 60 | +rn = @reaction_network rType begin |
| 61 | + (kB, kD), X + Y ↔ XY |
| 62 | +end kB, kD |
| 63 | +p = [2.0, 1.0] |
| 64 | +``` |
| 65 | +The parameter set `p` must be passed to the problem constructor. The parameter values can be changed after the reaction network is defined. |
| 66 | + |
| 67 | +The reaction rate do not need to be constant, but maybe depend on the concentration of the reactants. |
| 68 | +```julia |
| 69 | +rn = @reaction_network rType begin |
| 70 | + (1.0,2XY), X + Y ↔ XY |
| 71 | +end |
| 72 | +``` |
| 73 | +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`). |
| 74 | +```julia |
| 75 | +rn = @reaction_network rType begin |
| 76 | + (1.0,hill(XY,1.5,2.0,2)), X + Y ↔ XY |
| 77 | +end |
| 78 | +``` |
| 79 | +By using the `@reaction_func` macro it is possible to define your own functions, which may then be used when creating new reaction networks. |
| 80 | +```julia |
| 81 | +@reaction_func hill2(x, v, k) = v*x^2/(k^2+x^2) |
| 82 | +@reaction_network macro can see. |
| 83 | +rn = @reaction_network rType begin |
| 84 | + (1.0,hill2(XY,1.5,2.0)), X + Y ↔ XY |
| 85 | +end |
| 86 | +``` |
| 87 | + |
| 88 | +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 |
| 89 | +rn = @reaction_network rType begin |
| 90 | + 2.0, X + Y --> XY |
| 91 | + 2.0*X*Y X + Y ⟾ XY |
| 92 | +end |
| 93 | +will both have reaction rate equal to 2[X][Y]. |
| 94 | + |
| 95 | +Once a reaction network has been created it can be passed as input to either one of the `ODEProblem`, `SDEProblem` or `JumpProblem` constructors. |
| 96 | +```julia |
| 97 | + probODE = ODEProblem(rn, args...; kwargs...) |
| 98 | + probSDE = SDEProblem(rn, args...; kwargs...) |
| 99 | + probJump = JumpProblem(prob,aggregator::Direct,rn) |
| 100 | +``` |
| 101 | +the output problems may then be used as normal input to the solvers of the `DifferentialEquations` package. |
| 102 | + |
| 103 | +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. |
| 104 | +```julia |
| 105 | +rn = @reaction_network \eta begin |
| 106 | + 2.0, X + Y ↔ XY |
| 107 | +end |
| 108 | +p = [0.5] |
| 109 | +``` |
| 110 | +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). |
| 111 | + |
| 112 | +It is possible to access expressions corresponding to the functions determining the deterministic and stochastic development of the network using. |
| 113 | +```julia |
| 114 | + f_expr = rn.f_func |
| 115 | + g_expr = rn.g_func |
| 116 | +``` |
| 117 | +This can e.g. be used to generate LaTeX code corresponding to the system. |
0 commit comments