@@ -31,43 +31,57 @@ We first define the parameters, variables, differential equations and the output
3131``` @example SI
3232using 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
3735D = 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
5470After 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
6378We can see that all states (except $x_7$) and all parameters are locally identifiable with probability 0.99.
6479
6580Let'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
7387Notice 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
93107Global 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
98110using 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
117148We 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
122153using 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
125156D = 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
144197Both parameters ` b, c ` are globally identifiable with probability ` 0.9 ` in this case.
0 commit comments