Skip to content

Commit 392b8aa

Browse files
iliailmeriliailmer
committed
add identifiability tutorial
1 parent 98ebb0c commit 392b8aa

File tree

2 files changed

+138
-6
lines changed

2 files changed

+138
-6
lines changed

‎docs/make.jl

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -5,8 +5,8 @@ makedocs(
55
authors="Chris Rackauckas",
66
modules=[ModelingToolkit],
77
clean=true,doctest=false,
8-
format = Documenter.HTML(#analytics = "UA-90474609-3",
9-
assets = ["assets/favicon.ico"],
8+
format=Documenter.HTML(# analytics = "UA-90474609-3",
9+
assets=["assets/favicon.ico"],
1010
canonical="https://mtk.sciml.ai/stable/"),
1111
pages=[
1212
"Home" => "index.md",
@@ -19,12 +19,13 @@ makedocs(
1919
"tutorials/nonlinear.md",
2020
"tutorials/optimization.md",
2121
"tutorials/stochastic_diffeq.md",
22-
"tutorials/nonlinear_optimal_control.md"
22+
"tutorials/nonlinear_optimal_control.md",
23+
"tutorials/parameter_identifiability.md"
2324
],
2425
"ModelingToolkitize Tutorials" => Any[
2526
"mtkitize_tutorials/modelingtoolkitize.md",
2627
"mtkitize_tutorials/modelingtoolkitize_index_reduction.md",
27-
#"mtkitize_tutorials/sparse_jacobians",
28+
# "mtkitize_tutorials/sparse_jacobians",
2829
],
2930
"Basics" => Any[
3031
"basics/AbstractSystem.md",
@@ -49,6 +50,6 @@ makedocs(
4950
)
5051

5152
deploydocs(
52-
repo = "github.com/SciML/ModelingToolkit.jl.git";
53-
push_preview = true
53+
repo="github.com/SciML/ModelingToolkit.jl.git";
54+
push_preview=true
5455
)

‎docs/src/tutorials/parameter_identifiability.md

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
# Parameter Identifiability in ODE Models
2+
3+
Using ordinary differential equations in modeling processes is commonplace and the challenge of parameter identifiability is one of the key design challenges. In this tutorial, we will show how to use `StructuralIdentifiability.jl` with `ModelingToolkit.jl` to assess parameter identifiability.
4+
5+
We will start with determining local identifiability, where a parameter is known up to finitely many values, and then proceed to determining global identifiability properties, that is, which parameters can be identified uniquely.
6+
7+
## Local Identifiability
8+
### Input System
9+
10+
We will consider a simple two-species competition model
11+
12+
$$\begin{cases}
13+
\frac{d\,x_4}{d\,t} = - \frac{k_5 x_4}{k_6 + x_4},\\
14+
\frac{d\,x_5}{d\,t} = \frac{k_5 x_4}{k_6 + x_4} - \frac{k_7 x_5}{(k_8 + x_5 + x_6)},\\
15+
\frac{d\,x_6}{d\,t} = \frac{k_7 x_5}{(k_8 + x_5 + x_6)} - \frac{k_9 x_6 (k_{10} - x_6) }{k_{10}},\\
16+
\frac{d\,x_7}{d\,t} = \frac{k_9 x_6 (k_{10} - x_6)}{ k_{10}},\\
17+
y_1 = x_4,\\
18+
y_2 = x_5\end{cases}$$
19+
20+
This model describes the biohydrogenation[^1] process[^2] with unknown initial conditions.
21+
22+
### Using the `ODESystem` object
23+
To define the system in Julia, we use `ModelingToolkit.jl`.
24+
25+
We first define the parameters, variables, differential equations and the output equations. Notice that the system does not have any input functions, so inputs will be an empty array.
26+
27+
```@example
28+
using StructuralIdentifiability, ModelingToolkit
29+
30+
# define parameters and variables
31+
@variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t)
32+
@parameters k5 k6 k7 k8 k9 k10
33+
D = Differential(t)
34+
35+
# define equations
36+
eqs = [
37+
D(x4) ~ - k5 * x4 / (k6 + x4),
38+
D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5/(k8 + x5 + x6),
39+
D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10,
40+
D(x7) ~ k9 * x6 * (k10 - x6) / k10
41+
]
42+
43+
# define observed functions
44+
observed = [
45+
y1 ~ x4,
46+
y2 ~ x5
47+
]
48+
49+
# define the system
50+
de = ODESystem(eqs, t, [x4, x5, x6, x7], [k5, k6, k7, k8, k9, k10], observed=observed, name=:Biohydrogenation)
51+
52+
# no input functions:
53+
inputs = []
54+
55+
# we want to check everything
56+
to_check = []
57+
58+
# query local identifiability
59+
# we pass the ode-system
60+
local_id_all = assess_local_identifiability(de, inputs, to_check, 0.99)
61+
62+
# let's try to check specific parameters and their combinations
63+
to_check = [k5, k7, k10/k9, k5+k6]
64+
local_id_some = assess_local_identifiability(de, inputs, to_check, 0.99)
65+
66+
```
67+
68+
Notice that in this case, everything (except the state variable $x_7$) is locally identifiable, including combinations such as $k_{10}/k_9, k_5+k_6$
69+
70+
## Global Identifiability
71+
72+
In this tutorial, let us cover an example problem of querying the ODE for globally identifiable parameters.
73+
74+
### Input System
75+
76+
Let us consider the following four-dimensional model with two outputs:
77+
78+
$\begin{cases}x'(t) = lm - d \, x(t) - \beta \, x(t) \, v(t),\\
79+
y'(t) = \beta \, x(t) \, v(t) - a \, y(t),\\
80+
v'(t) = k \, y(t) - u \, v(t),\\
81+
w'(t) = c \, x(t) \, y(t) \, w(t) - c \, q \, y(t) \, w(t) - b \, w(t),\\
82+
z'(t) = c \, q \, y(t) \, w(t) - h \, z(t),\\
83+
y_1(t) = w(t),\\
84+
y_2(t) = z(t)\end{cases}$
85+
86+
This model describes HIV dynamics[^1]. Let us run a global identifiability check on this model to get the result with probability of correctness being `p=0.999`. To do this, we will use `assess_identifiability(ode, p)` function.
87+
88+
Global identifiability needs information about local identifiability first, hence the function we chose here will take care of that extra step for us.
89+
90+
```@repl
91+
using StructuralIdentifiability
92+
93+
ode = @ODEmodel(
94+
x'(t) = lm - d * x(t) - beta * x(t) * v(t),
95+
y'(t) = beta * x(t) * v(t) - a * y(t),
96+
v'(t) = k * y(t) - u * v(t),
97+
w'(t) = c * x(t) * y(t) * w(t) - c * q * y(t) * w(t) - b * w(t),
98+
z'(t) = c * q * y(t) * w(t) - h * z(t),
99+
y1(t) = w(t),
100+
y2(t) = z(t)
101+
)
102+
@time global_id = assess_identifiability(ode, 0.999)
103+
```
104+
105+
Now let us compare the same system but with probability being `p=0.99`. We will see a reduction in runtime:
106+
107+
```@repl
108+
using StructuralIdentifiability
109+
110+
ode = @ODEmodel(
111+
x'(t) = lm - d * x(t) - beta * x(t) * v(t),
112+
y'(t) = beta * x(t) * v(t) - a * y(t),
113+
v'(t) = k * y(t) - u * v(t),
114+
w'(t) = c * x(t) * y(t) * w(t) - c * q * y(t) * w(t) - b * w(t),
115+
z'(t) = c * q * y(t) * w(t) - h * z(t),
116+
y1(t) = w(t),
117+
y2(t) = z(t)
118+
)
119+
@time global_id = assess_identifiability(ode, 0.99)
120+
```
121+
122+
Indeed, notice how much quicker we obtained the result with 99% correctness guarantee! This illustrates the fact that you may sometimes sacrifice probability slightly to get results much faster.
123+
124+
[^1]:
125+
> R. Munoz-Tamayo, L. Puillet, J.B. Daniel, D. Sauvant, O. Martin, M. Taghipoor, P. Blavy [*Review: To be or not to be an identifiable model. Is this a relevant question in animal science modelling?*](https://doi.org/10.1017/S1751731117002774), Animal, Vol 12 (4), 701-712, 2018. The model is the ODE system (3) in Supplementary Material 2, initial conditions are assumed to be unknown.
126+
127+
[^2]:
128+
> Moate P.J., Boston R.C., Jenkins T.C. and Lean I.J., [*Kinetics of Ruminal Lipolysis of Triacylglycerol and Biohydrogenationof Long-Chain Fatty Acids: New Insights from Old Data*](doi:10.3168/jds.2007-0398), Journal of Dairy Science 91, 731–742, 2008
129+
130+
[^3]:
131+
> D. Wodarz, M. Nowak, [*Specific therapy regimes could lead to long-term immunological control of HIV*](https://doi.org/10.1073/pnas.96.25.14464), PNAS December 7, 1999 96 (25) 14464-14469;

0 commit comments

Comments
 (0)