Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
# Create a model
The package offers different ways to create models based on CRNs.
## Load a pre-written model
A bunch of models are already writtne within the package. If `str_model::String` is the name of an implemented model, then `load_model(str_model)` creates a variable with name `str_model`.
```julia
load_model("poisson")
```
Available models are listed below.
## Define a Chemical Reaction Network
Let's consider the Chemical Reaction Network of the SIR model:
``
Infection: S + I \xrightarrow{k_i} 2I \\
Recovery: I \xrightarrow{k_r} R
``
The macro `@network_model` creates easily a CTMC stored in a `ContinuousTimeModel` variable based on this formalism.
```julia
julia> easy_sir = @network_model begin
Infection: (S + I => 2I, ki*I*S)
Recovery: (I => R, kr*I)
end "My awesome SIR"
My_awesome_SIRModel <: ContinuousTimeModel model
- variables :
* I (index = 2 in state space)
* R (index = 3 in state space)
* S (index = 1 in state space)
- parameters :
* ki (index = 1 in parameter space)
* kr (index = 2 in parameter space)
- transitions : Infection,Recovery
- observed variables :
* S (index = 1 in observed state space, index = 1 in state space)
* I (index = 2 in observed state space, index = 2 in state space)
* R (index = 3 in observed state space, index = 3 in state space)
p = [0.0, 0.0]
x0 = [0, 0, 0]
t0 = 0.0
time bound = Inf
```
In the first reaction, `ki*I*S` is the reaction rate of the reaction `Infection`. This model is almost ready to use, we have to set the initial state and the parameters.
```julia
julia> set_param!(easy_sir, [0.0012, 0.05])
set_x0!(easy_sir, [95, 5, 0])
σ = simulate(easy_sir)
load_plots()
plot(σ)
```
## Manually (advanced)
This page is intented to advanced uses of the package, in order to use
### Based on an existing model
### From scratch
When the above cases don't fit your application one can create manually a `ContinuousTimeModel`. Let's take a look about the signature of the constructor method:
```julia
function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict{VariableModel,Int},
map_param_idx::Dict{ParameterModel,Int}, transitions::Vector{<:Transition},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
f!::Function, isabsorbing::Function; kwargs)
```
Let's construct an SIR model manually. First, one has to specify the dimensions of the state space and the parameter space.
```julia
dim_state_sir, dim_params_sir = 3, 2
```
`map_var_idx` is a dictionary that maps each model variable (represented by a `Symbol`) to an index in the state space.
```julia
map_var_idx_sir = Dict(:S => 1, :I => 2, :R => 3)
```
`map_var_params` is the equivalent of `map_var_idx` for parameters.
```julia
map_params_idx_sir = Dict(:ki => 1, :kr => 2)
```
`transitions` are the transitions/reactions of the model (vector of `Symbol`), `p`, `x0` and `t0` are respectively the parameters, the initial state and initial time of the model.
```julia
transitions_sir = [:Infection, :Recovery]
p_sir = [0.0012, 0.05]
x0_sir = [95, 5, 0]
t0_sir = 0.0
```
The two last arguments are functions, the first one, called `f!` must have the signature:
```julia
function f!(xnplus1::Vector{Int}, ptr_t::Vector{Float64}, ptr_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
```
It should return nothing. `xnplus1`, `ptr_t` and `ptr_tr` are vectors where the next values are stored. `ptr_t` is of length 1 and stores the next time value (`ptr_t[1] = tn + delta_t`) whereas `ptr_tr` stores the name of the next transition/reaction (`ptr_tr[1] = :Infection` for example). This function is implemented in the package as:
```julia
@everywhere function sir_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
@inbounds a1 = p[1] * xn[1] * xn[2]
@inbounds a2 = p[2] * xn[2]
l_a = (a1, a2)
asum = sum(l_a)
if asum == 0.0
copyto!(xnplus1, xn)
return nothing
end
nu_1 = (-1, 1, 0)
nu_2 = (0, -1, 1)
l_nu = (nu_1, nu_2)
l_str_R = (:Infection, :Recovery)
u1 = rand()
u2 = rand()
tau = - log(u1) / asum
b_inf = 0.0
b_sup = a1
reaction = 0
for i = 1:2
if b_inf < asum*u2 < b_sup
reaction = i
break
end
@inbounds b_inf += l_a[i]
@inbounds b_sup += l_a[i+1]
end
nu = l_nu[reaction]
for i = 1:3
@inbounds xnplus1[i] = xn[i]+nu[i]
end
@inbounds l_t[1] = tn + tau
@inbounds l_tr[1] = l_str_R[reaction]
end
```
i
The second function called `isaborbing` must have the signature:
```julia
isabsorbing(p::Vector{Float64}, xn::Vector{Int})
```
This function checks if the state `xn` is an absorbing state according to the model parametrised by `p`. It has to return true or false.
For a CTMC, a state is an absorbing state if the total exit rate is zero. In the case of the SIR model;
```julia
@everywhere sir_isabsorbing(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
```
Finally one sets the observed variables and the model can be created. The following lines creates a new type `TryhardSIRModel <: ContinuousTimeModel`, and the core of simulation.
```julia
g_sir = [:I]
# Generates simulate method for the new model
@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:TryhardSIRModel))
@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:TryhardSIRModel))
@everywhere @eval $(MarkovProcesses.generate_code_simulation(:TryhardSIRModel, :sir_f!, :sir_isabsorbing))
tryhard_sir = TryhardSIRModel(dim_state_sir, dim_params_sir,
map_var_idx_sir, map_params_idx_sir,
transitions_sir, p_sir, x0_sir, t0_sir,
:sir_f!, :sir_isabsorbing; g = g_sir)
σ = simulate(tryhard_sir)
```
## List of pre-written models
- `load_model("poisson")`: Poisson process
- `load_model("ER")`: Michaelis-Menten kinetics (Enzymatic Reactions)
- `load_model("SIR")`: Susceptible-Infected-Removed
- `load_model("doping_3way_oscillator")`: Three-way oscillator with doping reactions
- `load_model("repressilator")`: A repressilator model