Skip to content

Commit 6ac545d

Browse files
committed
docs(refactor): update structural identifiability tutorial with @mtkmodel
1 parent c5bd821 commit 6ac545d

File tree

1 file changed

+102
-49
lines changed

1 file changed

+102
-49
lines changed

docs/src/tutorials/parameter_identifiability.md

Lines changed: 102 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -31,43 +31,57 @@ We first define the parameters, variables, differential equations and the output
3131
```@example SI
3232
using StructuralIdentifiability, ModelingToolkit
3333
34-
# define parameters and variables
35-
@variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t)
36-
@parameters k5 k6 k7 k8 k9 k10
34+
@variables t
3735
D = Differential(t)
3836
39-
# define equations
40-
eqs = [
41-
D(x4) ~ -k5 * x4 / (k6 + x4),
42-
D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6),
43-
D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10,
44-
D(x7) ~ k9 * x6 * (k10 - x6) / k10,
45-
]
46-
47-
# define the output functions (quantities that can be measured)
48-
measured_quantities = [y1 ~ x4, y2 ~ x5]
37+
@mtkmodel Biohydrogenation begin
38+
@variables begin
39+
x4(t)
40+
x5(t)
41+
x6(t)
42+
x7(t)
43+
y1(t), [output = true]
44+
y2(t), [output = true]
45+
end
46+
@parameters begin
47+
k5
48+
k6
49+
k7
50+
k8
51+
k9
52+
k10
53+
end
54+
# define equations
55+
@equations begin
56+
D(x4) ~ -k5 * x4 / (k6 + x4)
57+
D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6)
58+
D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10
59+
D(x7) ~ k9 * x6 * (k10 - x6) / k10
60+
y1 ~ x4
61+
y2 ~ x5
62+
end
63+
end
4964
5065
# define the system
51-
de = ODESystem(eqs, t, name = :Biohydrogenation)
66+
@named de = Biohydrogenation()
67+
de = complete(de)
5268
```
5369

5470
After that, we are ready to check the system for local identifiability:
5571

5672
```@example SI
5773
# query local identifiability
5874
# we pass the ode-system
59-
local_id_all = assess_local_identifiability(de, measured_quantities = measured_quantities,
60-
p = 0.99)
75+
local_id_all = assess_local_identifiability(de, p = 0.99)
6176
```
6277

6378
We can see that all states (except $x_7$) and all parameters are locally identifiable with probability 0.99.
6479

6580
Let's try to check specific parameters and their combinations
6681

6782
```@example SI
68-
to_check = [k5, k7, k10 / k9, k5 + k6]
69-
local_id_some = assess_local_identifiability(de, measured_quantities = measured_quantities,
70-
funcs_to_check = to_check, p = 0.99)
83+
to_check = [de.k5, de.k7, de.k10 / de.k9, de.k5 + de.k6]
84+
local_id_some = assess_local_identifiability(de, funcs_to_check = to_check, p = 0.99)
7185
```
7286

7387
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$
@@ -92,26 +106,43 @@ We will run a global identifiability check on this enzyme dynamics[^3] model. We
92106

93107
Global identifiability needs information about local identifiability first, but the function we chose here will take care of that extra step for us.
94108

95-
__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).
96-
97109
```@example SI2
98110
using StructuralIdentifiability, ModelingToolkit
99-
@parameters b c a beta g delta sigma
100-
@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t)
101-
D = Differential(t)
102-
103-
eqs = [
104-
D(x1) ~ -b * x1 + 1 / (c + x4),
105-
D(x2) ~ a * x1 - beta * x2,
106-
D(x3) ~ g * x2 - delta * x3,
107-
D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3,
108-
]
109111
110-
measured_quantities = [y ~ x1 + x2, y2 ~ x2]
111-
112-
ode = ODESystem(eqs, t, name = :GoodwinOsc)
112+
@variables t
113+
D = Differential(t)
113114
114-
global_id = assess_identifiability(ode, measured_quantities = measured_quantities)
115+
@mtkmodel GoodwinOsc begin
116+
@parameters begin
117+
b
118+
c
119+
α
120+
β
121+
γ
122+
δ
123+
σ
124+
end
125+
@variables begin
126+
x1(t)
127+
x2(t)
128+
x3(t)
129+
x4(t)
130+
y(t), [output = true]
131+
y2(t), [output = true]
132+
end
133+
@equations begin
134+
D(x1) ~ -b * x1 + 1 / (c + x4)
135+
D(x2) ~ α * x1 - β * x2
136+
D(x3) ~ γ * x2 - δ * x3
137+
D(x4) ~ σ * x4 * (γ * x2 - δ * x3) / x3
138+
y ~ x1 + x2
139+
y2 ~ x2
140+
end
141+
end
142+
143+
@named ode = GoodwinOsc()
144+
145+
global_id = assess_identifiability(ode)
115146
```
116147

117148
We can see that only parameters `a, g` are unidentifiable, and everything else can be uniquely recovered.
@@ -120,25 +151,47 @@ Let us consider the same system but with two inputs, and we will find out identi
120151

121152
```@example SI3
122153
using StructuralIdentifiability, ModelingToolkit
123-
@parameters b c a beta g delta sigma
124-
@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t) u1(t) [input = true] u2(t) [input = true]
154+
155+
@variables t
125156
D = Differential(t)
126157
127-
eqs = [
128-
D(x1) ~ -b * x1 + 1 / (c + x4),
129-
D(x2) ~ a * x1 - beta * x2 - u1,
130-
D(x3) ~ g * x2 - delta * x3 + u2,
131-
D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3,
132-
]
133-
measured_quantities = [y ~ x1 + x2, y2 ~ x2]
158+
@mtkmodel GoodwinOscillator begin
159+
@parameters begin
160+
b
161+
c
162+
α
163+
β
164+
γ
165+
δ
166+
σ
167+
end
168+
@variables begin
169+
x1(t)
170+
x2(t)
171+
x3(t)
172+
x4(t)
173+
y(t), [output = true]
174+
y2(t), [output = true]
175+
u1(t), [input = true]
176+
u2(t), [input = true]
177+
end
178+
@equations begin
179+
D(x1) ~ -b * x1 + 1 / (c + x4)
180+
D(x2) ~ α * x1 - β * x2 - u1
181+
D(x3) ~ γ * x2 - δ * x3 + u2
182+
D(x4) ~ σ * x4 * (γ * x2 - δ * x3) / x3
183+
y ~ x1 + x2
184+
y2 ~ x2
185+
end
186+
end
187+
188+
@named ode = GoodwinOscillator()
189+
ode = complete(ode)
134190
135191
# check only 2 parameters
136-
to_check = [b, c]
137-
138-
ode = ODESystem(eqs, t, name = :GoodwinOsc)
192+
to_check = [ode.b, ode.c]
139193
140-
global_id = assess_identifiability(ode, measured_quantities = measured_quantities,
141-
funcs_to_check = to_check, p = 0.9)
194+
global_id = assess_identifiability(ode, funcs_to_check = to_check, p = 0.9)
142195
```
143196

144197
Both parameters `b, c` are globally identifiable with probability `0.9` in this case.

0 commit comments

Comments
 (0)