Newer
Older
<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Create a model · MarkovProcesses.jl</title><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.045/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.4/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.4/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.4/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.24/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><a class="docs-logo" href="index.html"><img src="assets/logo.png" alt="MarkovProcesses.jl logo"/></a><div class="docs-package-name"><span class="docs-autofit"><a href="index.html">MarkovProcesses.jl</a></span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="index.html">Home</a></li><li><a class="tocitem" href="starting.html">Getting Started</a></li><li class="is-active"><a class="tocitem" href="create_model.html">Create a model</a><ul class="internal"><li><a class="tocitem" href="#Load-a-pre-written-model"><span>Load a pre-written model</span></a></li><li><a class="tocitem" href="#Define-a-Chemical-Reaction-Network"><span>Define a Chemical Reaction Network</span></a></li><li><a class="tocitem" href="#Manually-(advanced)"><span>Manually (advanced)</span></a></li><li><a class="tocitem" href="#List-of-pre-written-models"><span>List of pre-written models</span></a></li></ul></li><li><span class="tocitem">API</span><ul><li><a class="tocitem" href="api/model.html">Model</a></li><li><a class="tocitem" href="api/trajectory.html">Trajectory</a></li><li><a class="tocitem" href="api/abc.html">Approximate Bayesian Computation</a></li><li><a class="tocitem" href="api/plots.html">Plots</a></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="create_model.html">Create a model</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="create_model.html">Create a model</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://github.com//blob/master/docs/src/create_model.md" title="Edit on GitHub"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitHub</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Create-a-model"><a class="docs-heading-anchor" href="#Create-a-model">Create a model</a><a id="Create-a-model-1"></a><a class="docs-heading-anchor-permalink" href="#Create-a-model" title="Permalink"></a></h1><p>The package offers different ways to create models based on CRNs.</p><h2 id="Load-a-pre-written-model"><a class="docs-heading-anchor" href="#Load-a-pre-written-model">Load a pre-written model</a><a id="Load-a-pre-written-model-1"></a><a class="docs-heading-anchor-permalink" href="#Load-a-pre-written-model" title="Permalink"></a></h2><p>A bunch of models are already writtne within the package. If <code>str_model::String</code> is the name of an implemented model, then <code>load_model(str_model)</code> creates a variable with name <code>str_model</code>.</p><pre><code class="language-julia hljs">load_model("poisson")</code></pre><p>Available models are listed below.</p><h2 id="Define-a-Chemical-Reaction-Network"><a class="docs-heading-anchor" href="#Define-a-Chemical-Reaction-Network">Define a Chemical Reaction Network</a><a id="Define-a-Chemical-Reaction-Network-1"></a><a class="docs-heading-anchor-permalink" href="#Define-a-Chemical-Reaction-Network" title="Permalink"></a></h2><p>Let's consider the Chemical Reaction Network of the SIR model:</p><p>`` Infection: S + I \xrightarrow{k<em>i} 2I \
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
Recovery: I \xrightarrow{k</em>r} R ``</p><p>The macro <code>@network_model</code> creates easily a CTMC stored in a <code>ContinuousTimeModel</code> variable based on this formalism.</p><pre><code class="language-julia hljs">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</code></pre><p>In the first reaction, <code>ki*I*S</code> is the reaction rate of the reaction <code>Infection</code>. This model is almost ready to use, we have to set the initial state and the parameters.</p><pre><code class="language-julia hljs">julia> set_param!(easy_sir, [0.0012, 0.05])
set_x0!(easy_sir, [95, 5, 0])
σ = simulate(easy_sir)
load_plots()
plot(σ)</code></pre><p><img src="assets/sir_trajectory.png" alt="Plot of a simulated SIR trajectory"/></p><h2 id="Manually-(advanced)"><a class="docs-heading-anchor" href="#Manually-(advanced)">Manually (advanced)</a><a id="Manually-(advanced)-1"></a><a class="docs-heading-anchor-permalink" href="#Manually-(advanced)" title="Permalink"></a></h2><p>This page is intented to advanced uses of the package, in order to use</p><h3 id="Based-on-an-existing-model"><a class="docs-heading-anchor" href="#Based-on-an-existing-model">Based on an existing model</a><a id="Based-on-an-existing-model-1"></a><a class="docs-heading-anchor-permalink" href="#Based-on-an-existing-model" title="Permalink"></a></h3><h3 id="From-scratch"><a class="docs-heading-anchor" href="#From-scratch">From scratch</a><a id="From-scratch-1"></a><a class="docs-heading-anchor-permalink" href="#From-scratch" title="Permalink"></a></h3><p>When the above cases don't fit your application one can create manually a <code>ContinuousTimeModel</code>. Let's take a look about the signature of the constructor method:</p><pre><code class="language-julia hljs">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)</code></pre><p>Let's construct an SIR model manually. First, one has to specify the dimensions of the state space and the parameter space.</p><pre><code class="language-julia hljs">dim_state_sir, dim_params_sir = 3, 2</code></pre><p><code>map_var_idx</code> is a dictionary that maps each model variable (represented by a <code>Symbol</code>) to an index in the state space.</p><pre><code class="language-julia hljs">map_var_idx_sir = Dict(:S => 1, :I => 2, :R => 3)</code></pre><p><code>map_var_params</code> is the equivalent of <code>map_var_idx</code> for parameters.</p><pre><code class="language-julia hljs">map_params_idx_sir = Dict(:ki => 1, :kr => 2)</code></pre><p><code>transitions</code> are the transitions/reactions of the model (vector of <code>Symbol</code>), <code>p</code>, <code>x0</code> and <code>t0</code> are respectively the parameters, the initial state and initial time of the model.</p><pre><code class="language-julia hljs">transitions_sir = [:Infection, :Recovery]
p_sir = [0.0012, 0.05]
x0_sir = [95, 5, 0]
t0_sir = 0.0</code></pre><p>The two last arguments are functions, the first one, called <code>f!</code> must have the signature:</p><pre><code class="language-julia hljs">function f!(xnplus1::Vector{Int}, ptr_t::Vector{Float64}, ptr_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})</code></pre><p>It should return nothing. <code>xnplus1</code>, <code>ptr_t</code> and <code>ptr_tr</code> are vectors where the next values are stored. <code>ptr_t</code> is of length 1 and stores the next time value (<code>ptr_t[1] = tn + delta_t</code>) whereas <code>ptr_tr</code> stores the name of the next transition/reaction (<code>ptr_tr[1] = :Infection</code> for example). This function is implemented in the package as:</p><pre><code class="language-julia hljs">@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</code></pre><p>i The second function called <code>isaborbing</code> must have the signature:</p><pre><code class="language-julia hljs">isabsorbing(p::Vector{Float64}, xn::Vector{Int})</code></pre><p>This function checks if the state <code>xn</code> is an absorbing state according to the model parametrised by <code>p</code>. It has to return true or false.</p><p>For a CTMC, a state is an absorbing state if the total exit rate is zero. In the case of the SIR model;</p><pre><code class="language-julia hljs">@everywhere sir_isabsorbing(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0</code></pre><p>Finally one sets the observed variables and the model can be created. The following lines creates a new type <code>TryhardSIRModel <: ContinuousTimeModel</code>, and the core of simulation.</p><pre><code class="language-julia hljs">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)</code></pre><h2 id="List-of-pre-written-models"><a class="docs-heading-anchor" href="#List-of-pre-written-models">List of pre-written models</a><a id="List-of-pre-written-models-1"></a><a class="docs-heading-anchor-permalink" href="#List-of-pre-written-models" title="Permalink"></a></h2><ul><li><code>load_model("poisson")</code>: Poisson process</li><li><code>load_model("ER")</code>: Michaelis-Menten kinetics (Enzymatic Reactions)</li><li><code>load_model("SIR")</code>: Susceptible-Infected-Removed</li><li><code>load_model("doping_3way_oscillator")</code>: Three-way oscillator with doping reactions</li><li><code>load_model("repressilator")</code>: A repressilator model</li></ul></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="starting.html">« Getting Started</a><a class="docs-footer-nextpage" href="api/model.html">Model »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.24 on <span class="colophon-date" title="Monday 22 May 2023 14:03">Monday 22 May 2023</span>. Using Julia version 1.7.2.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>