Skip to content

Commit 7969d00

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #1243 from iliailmer/id-tutorial
Identifiability tutorial
2 parents ee767dd + d607cc4 commit 7969d00

File tree

2 files changed

+186
-6
lines changed

2 files changed

+186
-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: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,179 @@
1+
# Parameter Identifiability in ODE Models
2+
3+
Ordinary differential equations are commonly used for modeling real-world processes. The problem of parameter identifiability is one of the key design challenges for mathematical models. A parameter is said to be _identifiable_ if one can recover its value from experimental data. _Structural_ identifiabiliy is a theoretical property of a model that answers this question. In this tutorial, we will show how to use `StructuralIdentifiability.jl` with `ModelingToolkit.jl` to assess identifiability of parameters in ODE models. The theory behind `StructuralIdentifiability.jl` is presented in paper [^4].
4+
5+
We will start by illutrating **local identifiability** in which a parameter is known up to _finitely many values_, and then proceed to determining **global identifiability**, that is, which parameters can be identified _uniquely_.
6+
7+
To install `StructuralIdentifiability.jl`, simply run
8+
```julia
9+
using Pkg
10+
Pkg.add("StructuralIdentifiability")
11+
```
12+
13+
The package has a standalone data structure for ordinary differential equations but is also compatible with `ODESystem` type from `ModelingToolkit.jl`.
14+
15+
## Local Identifiability
16+
### Input System
17+
18+
We will consider the following model:
19+
20+
$$\begin{cases}
21+
\frac{d\,x_4}{d\,t} = - \frac{k_5 x_4}{k_6 + x_4},\\
22+
\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)},\\
23+
\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}},\\
24+
\frac{d\,x_7}{d\,t} = \frac{k_9 x_6 (k_{10} - x_6)}{ k_{10}},\\
25+
y_1 = x_4,\\
26+
y_2 = x_5\end{cases}$$
27+
28+
This model describes the biohydrogenation[^1] process[^2] with unknown initial conditions.
29+
30+
### Using the `ODESystem` object
31+
To define the ode system in Julia, we use `ModelingToolkit.jl`.
32+
33+
We first define the parameters, variables, differential equations and the output equations.
34+
```julia
35+
using StructuralIdentifiability, ModelingToolkit
36+
37+
# define parameters and variables
38+
@variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t)
39+
@parameters k5 k6 k7 k8 k9 k10
40+
D = Differential(t)
41+
42+
# define equations
43+
eqs = [
44+
D(x4) ~ - k5 * x4 / (k6 + x4),
45+
D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5/(k8 + x5 + x6),
46+
D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10,
47+
D(x7) ~ k9 * x6 * (k10 - x6) / k10
48+
]
49+
50+
# define the output functions (quantities that can be measured)
51+
measured_quantities = [y1 ~ x4, y2 ~ x5]
52+
53+
# define the system
54+
de = ODESystem(eqs, t, name=:Biohydrogenation)
55+
56+
```
57+
58+
After that we are ready to check the system for local identifiability:
59+
```julia
60+
# query local identifiability
61+
# we pass the ode-system
62+
local_id_all = assess_local_identifiability(de, measured_quantities=measured_quantities, p=0.99)
63+
# [ Info: Preproccessing `ModelingToolkit.ODESystem` object
64+
# 6-element Vector{Bool}:
65+
# 1
66+
# 1
67+
# 1
68+
# 1
69+
# 1
70+
# 1
71+
```
72+
We can see that all states (except $x_7$) and all parameters are locally identifiable with probability 0.99.
73+
74+
Let's try to check specific parameters and their combinations
75+
```julia
76+
to_check = [k5, k7, k10/k9, k5+k6]
77+
local_id_some = assess_local_identifiability(de, funcs_to_check=to_check, p=0.99)
78+
# 4-element Vector{Bool}:
79+
# 1
80+
# 1
81+
# 1
82+
# 1
83+
```
84+
85+
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$
86+
87+
## Global Identifiability
88+
89+
In this part tutorial, let us cover an example problem of querying the ODE for globally identifiable parameters.
90+
91+
### Input System
92+
93+
Let us consider the following four-dimensional model with two outputs:
94+
95+
$$\begin{cases}
96+
x_1'(t) = -b x_1(t) + \frac{1 }{ c + x_4(t)},\\
97+
x_2'(t) = \alpha x_1(t) - \beta x_2(t),\\
98+
x_3'(t) = \gamma x_2(t) - \delta x_3(t),\\
99+
x_4'(t) = \sigma x_4(t) \frac{(\gamma x_2(t) - \delta x_3(t))}{ x_3(t)},\\
100+
y(t) = x_1(t)
101+
\end{cases}$$
102+
103+
We will run a global identifiability check on this enzyme dynamics[^3] model. We will use the default settings: the probability of correctness will be `p=0.99` and we are interested in identifiability of all possible parameters
104+
105+
Global identifiability needs information about local identifiability first, but the function we chose here will take care of that extra step for us.
106+
107+
__Note__: as of writing this tutorial, UTF-symbols such as Greek characters are not supported by one of the project's dependencies, see (this issue)[https://github.com/SciML/StructuralIdentifiability.jl/issues/43].
108+
109+
```julia
110+
using StructuralIdentifiability, ModelingToolkit
111+
@parameters b c a beta g delta sigma
112+
@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t)
113+
D = Differential(t)
114+
115+
eqs = [
116+
D(x1) ~ -b * x1 + 1/(c + x4),
117+
D(x2) ~ a * x1 - beta * x2,
118+
D(x3) ~ g * x2 - delta * x3,
119+
D(x4) ~ sigma * x4 * (g * x2 - delta * x3)/x3
120+
]
121+
122+
measured_quantities = [y~x1+x2, y2~x2]
123+
124+
125+
ode = ODESystem(eqs, t, name=:GoodwinOsc)
126+
127+
@time global_id = assess_identifiability(ode, measured_quantities=measured_quantities)
128+
# 30.672594 seconds (100.97 M allocations: 6.219 GiB, 3.15% gc time, 0.01% compilation time)
129+
# Dict{Num, Symbol} with 7 entries:
130+
# a => :globally
131+
# b => :globally
132+
# beta => :globally
133+
# c => :globally
134+
# sigma => :globally
135+
# g => :nonidentifiable
136+
# delta => :globally
137+
```
138+
We can see that only parameters `a, g` are unidentifiable and everything else can be uniquely recovered.
139+
140+
Let us consider the same system but with two inputs and we will try to find out identifiability with probability `0.9` for parameters `c` and `b`:
141+
142+
```julia
143+
using StructuralIdentifiability, ModelingToolkit
144+
@parameters b c a beta g delta sigma
145+
@variables t x1(t) x2(t) x3(t) x4(t) y(t) u1(t) [input=true] u2(t) [input=true]
146+
D = Differential(t)
147+
148+
eqs = [
149+
D(x1) ~ -b * x1 + 1/(c + x4),
150+
D(x2) ~ a * x1 - beta * x2 - u1,
151+
D(x3) ~ g * x2 - delta * x3 + u2,
152+
D(x4) ~ sigma * x4 * (g * x2 - delta * x3)/x3
153+
]
154+
measured_quantities = [y~x1+x2, y2~x2]
155+
156+
# check only 2 parameters
157+
to_check = [b, c]
158+
159+
ode = ODESystem(eqs, t, name=:GoodwinOsc)
160+
161+
global_id = assess_identifiability(ode, measured_quantities=measured_quantities, funcs_to_check=to_check, p=0.9)
162+
# Dict{Num, Symbol} with 2 entries:
163+
# b => :globally
164+
# c => :globally
165+
```
166+
167+
Both parameters `b, c` are globally identifiable with probability `0.9` in this case.
168+
169+
[^1]:
170+
> 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.
171+
172+
[^2]:
173+
> 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
174+
175+
[^3]:
176+
> Goodwin, B.C. [*Oscillatory behavior in enzymatic control processes*](https://doi.org/10.1016/0065-2571(65)90067-1), Advances in Enzyme Regulation, Vol 3 (C), 425-437, 1965
177+
178+
[^4]:
179+
> Dong, R., Goodbrake, C., Harrington, H. A., & Pogudin, G. [*Computing input-output projections of dynamical models with applications to structural identifiability*](https://arxiv.org/pdf/2111.00991). arXiv preprint arXiv:2111.00991.

0 commit comments

Comments
 (0)