You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: README.md
+94-82Lines changed: 94 additions & 82 deletions
Original file line number
Diff line number
Diff line change
@@ -3,15 +3,16 @@
3
3
[](https://gitter.im/JuliaDiffEq/Lobby?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)
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.
15
+
The `@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 `ODEProblem`, `SteadyStateProblem`, `SDEProblem` and `JumpProblem` constructors.
15
16
16
17
The basic syntax is
17
18
```julia
@@ -20,74 +21,70 @@ rn = @reaction_network rType begin
20
21
1.0, XY --> Z
21
22
end
22
23
```
23
-
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`).
24
+
where each line corresponds to a chemical reaction. The (optional) input `rType` designates the type of this instance (all instances will inherit from the abstract type `AbstractReactionNetwork`).
24
25
25
-
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
26
-
```julia
27
-
rn1 =@reaction_network rType begin
28
-
2.0, X + Y → XY
29
-
1.0, XY > Z
30
-
1.0, X + Y ← XY
31
-
0.5, XY < Z
32
-
end
33
-
rn1 =@reaction_network rType begin
34
-
(2.0,1.0), X + Y ↔ XY
35
-
(1.0, 0.5), XY ⟷ Z
36
-
end
37
-
```
38
-
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.
39
-
```julia
40
-
rn1 =@reaction_network rType begin
41
-
2.0, 2X -->0
42
-
2.0, ∅ --> X
43
-
end
44
-
```
45
-
Multiple reactions can be declared in a single line
46
-
```julia
47
-
rn =@reaction_network rType begin
48
-
2.0, (X,Y) -->0#Identical to reactions [2.0, X --> 0] and [2.0, Y --> 0]
49
-
(2.0, 1.0), (X,Y) -->0#Identical to reactions [2.0, X --> 0] and [1.0, X --> 0]
50
-
2.0, (X1,Y1) --> (X2,Y2) #Identical to reactions [2.0, X1 --> X2] and [2.0, Y1 --> Y2]
51
-
(2.0,1.0), X + Y ↔ XY #Identical to reactions [2.0, X + Y --> XY] and [1.0, XY --> X + Y].
52
-
((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].
53
-
end
26
+
The DSL has many features:
27
+
* It supports many different arrow types, corresponding to different directions of reactions and different rate laws:
28
+
```julia
29
+
rn =@reaction_networkbegin
30
+
1.0, X + Y --> XY
31
+
1.0, X + Y → XY
32
+
1.0, XY ← X + Y
33
+
2.0, X + Y ↔ XY
34
+
end
54
35
```
55
-
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.
56
-
```julia
57
-
rn =@reaction_network rType begin
58
-
(kB, kD), X + Y ↔ XY
59
-
end kB, kD
60
-
p = [2.0, 1.0]
61
-
```
62
-
The parameter set `p` must be passed to the problem constructor. The parameter values can be changed after the reaction network is defined.
36
+
* It allows multiple reactions to be defined simultaneously on one line. The following two networks
37
+
are equivalent:
38
+
```julia
39
+
rn1 =@reaction_networkbegin
40
+
(1.0,2.0), (S1,S2) → P
41
+
end
42
+
rn2 =@reaction_networkbegin
43
+
1.0, S1 → P
44
+
2.0, S2 → P
45
+
end
46
+
```
47
+
* It allows the use of symbols to represent reaction rate parameters, with their numeric values specified during problem construction. i.e., the previous example could be given by
48
+
```julia
49
+
rn2 =@reaction_networkbegin
50
+
k1, S1 → P
51
+
k2, S2 → P
52
+
end k1 k2
53
+
```
54
+
with `k1` and `k2` corresponding to the reaction rates.
55
+
* Rate law functions can be pre-defined and used within the DSL:
56
+
```julia
57
+
@reaction_funcmyHill(x) =2.0*x^3/(x^3+1.5^3)
58
+
rn =@reaction_networkbegin
59
+
myHill(X), ∅ → X
60
+
end
61
+
```
62
+
* Pre-made rate laws for Hill and Michaelis-Menten reactions are provided:
63
+
```julia
64
+
rn1 =@reaction_networkbegin
65
+
hill(X,v,K,n), ∅ → X
66
+
mm(X,v,K), ∅ → X
67
+
end v K n
68
+
```
69
+
* Simple rate law functions of the species populations can be used within the DSL:
70
+
```julia
71
+
rn =@reaction_networkbegin
72
+
2.0*X^2, 0→ X + Y
73
+
gamma(Y)/5, X → ∅
74
+
pi*X/Y, Y → ∅
75
+
end
76
+
```
77
+
* It is possible to access expressions corresponding to the functions determining the deterministic and stochastic terms within resulting ODE, SDE or jump models using
78
+
```julia
79
+
f_expr = rn.f_func
80
+
g_expr = rn.g_func
81
+
affects = rn.jump_affect_expr
82
+
rates = rn.jump_rate_expr
83
+
```
84
+
These can be used to generate LaTeX code corresponding to the system using packages such as [`Latexify`](https://github.com/korsbo/Latexify.jl).
63
85
64
-
The reaction rate do not need to be constant, but maybe depend on the concentration of the reactants.
65
-
```julia
66
-
rn =@reaction_network rType begin
67
-
(1.0,2XY), X + Y ↔ XY
68
-
end
69
-
```
70
-
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`).
71
-
```julia
72
-
rn =@reaction_network rType begin
73
-
(1.0,hill(XY,1.5,2.0,2)), X + Y ↔ XY
74
-
end
75
-
```
76
-
By using the `@reaction_func` macro it is possible to define your own functions, which may then be used when creating new reaction networks.
77
-
```julia
78
-
@reaction_funchill2(x, v, k) = v*x^2/(k^2+x^2)
79
-
@reaction_networkmacro can see.
80
-
rn =@reaction_network rType begin
81
-
(1.0,hill2(XY,1.5,2.0)), X + Y ↔ XY
82
-
end
83
-
```
84
86
85
-
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
86
-
rn = @reaction_network rType begin
87
-
2.0, X + Y --> XY
88
-
2.0*X*Y X + Y ⟾ XY
89
-
end
90
-
will both have reaction rate equal to 2[X][Y].
87
+
## Simulating ODE, Steady-State, SDE and Jump Problems
91
88
92
89
Once a reaction network has been created it can be passed as input to either one of the `ODEProblem`, `SteadyStateProblem`, `SDEProblem` or `JumpProblem` constructors.
93
90
```julia
@@ -96,20 +93,35 @@ Once a reaction network has been created it can be passed as input to either one
96
93
probSDE =SDEProblem(rn, args...; kwargs...)
97
94
probJump =JumpProblem(prob, Direct(), rn)
98
95
```
99
-
the output problems may then be used as normal input to the solvers of the `DifferentialEquations` package.
96
+
The output problems may then be used as input to the solvers of [DifferentialEquations.jl](http://juliadiffeq.org/). *Note*, the noise used by the `SDEProblem` will correspond to the Chemical Langevin Equations.
100
97
101
-
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.
98
+
As an example, consider models for a simple birth-death process:
102
99
```julia
103
-
rn =@reaction_network\eta begin
104
-
2.0, X + Y ↔ XY
105
-
end
106
-
p = [0.5]
107
-
```
108
-
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).
100
+
rs =@reaction_networkbegin
101
+
c1, X -->2X
102
+
c2, X -->0
103
+
c3, 0--> X
104
+
end c1 c2 c3
105
+
params = (1.0,2.0,50.)
106
+
tspan = (0.,4.)
107
+
u0 = [5.]
109
108
110
-
It is possible to access expressions corresponding to the functions determining the deterministic and stochastic development of the network using.
111
-
```julia
112
-
f_expr = rn.f_func
113
-
g_expr = rn.g_func
109
+
# solve ODEs
110
+
oprob =ODEProblem(rs, u0, tspan, params)
111
+
osol =solve(oprob, Tsit5())
112
+
113
+
# solve for Steady-States
114
+
ssprob =SteadyStateProblem(rs, u0, params)
115
+
sssol =solve(ssprob, SSRootfind())
116
+
117
+
# solve Chemical Langevin SDEs
118
+
sprob =SDEProblem(rs, u0, tspan, params)
119
+
ssol =solve(sprob, EM(), dt=.01)
120
+
121
+
# solve JumpProblem
122
+
u0 = [5]
123
+
dprob =DiscreteProblem(u0, tspan, params)
124
+
jprob =JumpProblem(dprob, Direct(), rs)
125
+
jsol =solve(jprob, SSAStepper())
114
126
```
115
-
This can e.g. be used to generate LaTeX code corresponding to the system.
0 commit comments