Skip to content

Commit cae7118

Browse files
authored
Merge pull request #984 from SciML/improve_si_find_ident_funcs_display
Improve conservation law handling in `find_identifiable_functions`
2 parents 031dc1b + fa1d8c0 commit cae7118

File tree

2 files changed

+31
-40
lines changed

2 files changed

+31
-40
lines changed

docs/src/inverse_problems/structural_identifiability.md

Lines changed: 12 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -21,9 +21,9 @@ StructuralIdentifiability currently assesses identifiability for the first two o
2121
!!! note
2222
Currently, the StructuralIdentifiability.jl extension only considers structural identifiability for the ODE generated by the reaction rate equation. It is possible that for the SDE model (generated by the chemical Langevin equation) and the jump model (generated by stochastic chemical kinetics) the identifiability of model quantities is different.
2323

24-
## Global identifiability analysis
24+
## [Global identifiability analysis](@id structural_identifiability_gi)
2525

26-
### Basic example
26+
### [Basic example](@id structural_identifiability_gi_example)
2727

2828
Global identifiability can be assessed using the `assess_identifiability` function. For each model quantity (parameters and variables), it will assess whether they are:
2929
- Globally identifiable.
@@ -48,7 +48,7 @@ assess_identifiability(gwo; measured_quantities = [:M, :P, :E], loglevel = Loggi
4848
```
4949
in which case all species trajectories and parameters become identifiable.
5050

51-
### Indicating known parameters
51+
### [Indicating known parameters](@id structural_identifiability_gi_known_ps)
5252
In the previous case we assumed that all parameters are unknown, however, this is not necessarily true. If there are parameters with known values, we can supply these using the `known_p` argument. Providing this additional information might also make other, previously unidentifiable, parameters identifiable. Let us consider the previous example, where we measure the concentration of $M$ only, but now assume we also know the production rate of $E$ ($pₑ$):
5353
```@example si1
5454
assess_identifiability(gwo; measured_quantities = [:M], known_p = [:pₑ], loglevel = Logging.Error)
@@ -57,7 +57,7 @@ Not only does this turn the previously non-identifiable `pₑ` (globally) identi
5757

5858
To, in a similar manner, indicate that certain initial conditions are known is a work in progress. Hopefully this feature should be an available in the near future.
5959

60-
### Providing non-trivial measured quantities
60+
### [Providing non-trivial measured quantities](@id structural_identifiability_gi_nontrivial_mq)
6161
Sometimes, ones may not have measurements of species, but rather some combinations of species (or possibly parameters). To account for this, `measured_quantities` accepts any algebraic expression (and not just single species). To form such expressions, species and parameters have to first be `@unpack`'ed from the model. Say that we have a model where an enzyme ($E$) is converted between an active and inactive form, which in turns activates the production of a product, $P$:
6262
```@example si1
6363
rs = @reaction_network begin
@@ -72,30 +72,30 @@ assess_identifiability(rs; measured_quantities = [Eᵢ + Eₐ, :P], loglevel = L
7272
nothing # hide
7373
```
7474

75-
### Assessing identifiability for specified quantities only
75+
### [Assessing identifiability for specified quantities only](@id structural_identifiability_gi_ftc)
7676
By default, StructuralIdentifiability assesses identifiability for all parameters and variables. It is, however, possible to designate precisely which quantities you want to check using the `funcs_to_check` option. This both includes selecting a smaller subset of parameters and variables to check, or defining customised expressions. Let us consider the Goodwind from previously, and say that we would like to check whether the production parameters ($pₘ$, $pₑ$, and $pₚ$) and the total amount of the three species ($P + M + E$) are identifiable quantities. Here, we would first unpack these (allowing us to form algebraic expressions) and then use the following code:
7777
```@example si1
7878
@unpack pₘ, pₑ, pₚ, M, E, P = gwo
7979
assess_identifiability(gwo; measured_quantities = [:M], funcs_to_check = [pₘ, pₑ, pₚ, M + E + P], loglevel = Logging.Error)
8080
nothing # hide
8181
```
8282

83-
### Probability of correctness
83+
### [Probability of correctness](@id structural_identifiability_gi_probs)
8484
The identifiability methods used can, in theory, produce erroneous results. However, it is possible to adjust the lower bound for the probability of correctness using the argument `prob_threshold` (by default set to `0.99`, that is, at least a $99\%$ chance of correctness). We can e.g. increase the bound through:
8585
```@example si1
8686
assess_identifiability(gwo; measured_quantities=[:M], prob_threshold = 0.999, loglevel = Logging.Error)
8787
nothing # hide
8888
```
8989
giving a minimum bound of $99.9\%$ chance of correctness. In practise, the bounds used by StructuralIdentifiability are very conservative, which means that while the minimum guaranteed probability of correctness in the default case is $99\%$, in practise it is much higher. While increasing the value of `prob_threshold` increases the certainty of correctness, it will also increase the time required to assess identifiability.
9090

91-
## Local identifiability analysis
91+
## [Local identifiability analysis](@id structural_identifiability_lit)
9292
Local identifiability can be assessed through the `assess_local_identifiability` function. While this is already determined by `assess_identifiability`, assessing local identifiability only has the advantage that it is easier to compute. Hence, there might be models where global identifiability analysis fails (or takes a prohibitively long time), where instead `assess_local_identifiability` can be used. This function takes the same inputs as `assess_identifiability` and returns, for each quantity, `true` if it is locally identifiable (or `false` if it is not). Here, for the Goodwind oscillator, we assesses it for local identifiability only:
9393
```@example si1
9494
assess_local_identifiability(gwo; measured_quantities = [:M], loglevel = Logging.Error)
9595
```
9696
We note that the results are consistent with those produced by `assess_identifiability` (with globally or locally identifiable quantities here all being assessed as at least locally identifiable).
9797

98-
## Finding identifiable functions
98+
## [Finding identifiable functions](@id structural_identifiability_identifiabile_funcs)
9999
Finally, StructuralIdentifiability provides the `find_identifiable_functions` function. Rather than determining the identifiability of each parameter and unknown of the model, it finds a set of identifiable functions, such as any other identifiable expression of the model can be generated by these. Let us again consider the Goodwind oscillator, using the `find_identifiable_functions` function we find that identifiability can be reduced to five globally identifiable expressions:
100100
```@example si1
101101
find_identifiable_functions(gwo; measured_quantities = [:M], loglevel = Logging.Error)
@@ -104,7 +104,7 @@ Again, these results are consistent with those produced by `assess_identifiabili
104104

105105
`find_identifiable_functions` tries to simplify its output functions to create nice expressions. The degree to which it does this can be adjusted using the `simplify` keywords. Using the `:weak`, `:standard` (default), and `:strong` arguments, increased simplification can be forced (at the expense of longer runtime).
106106

107-
## Creating StructuralIdentifiability compatible ODE models from Catalyst `ReactionSystem`s
107+
## [Creating StructuralIdentifiability compatible ODE models from Catalyst `ReactionSystem`s](@id structural_identifiability_si_odes)
108108
While the functionality described above covers the vast majority of analysis that user might want to perform, the StructuralIdentifiability package supports several additional features. While these does not have inherent Catalyst support, we do provide the `make_si_ode` function to simplify their use. Similar to the previous functions, it takes a `ReactionSystem`, lists of measured quantities, and known parameter values. The output is a [ODE of the standard form supported by StructuralIdentifiability](https://docs.sciml.ai/StructuralIdentifiability/stable/tutorials/creating_ode/#Defining-the-model-using-@ODEmodel-macro). It can be created using the following syntax:
109109
```@example si1
110110
si_ode = make_si_ode(gwo; measured_quantities = [:M])
@@ -116,28 +116,17 @@ print_for_DAISY(si_ode)
116116
nothing # hide
117117
```
118118

119-
## Notes on systems with conservation laws
119+
## [Notes on systems with conservation laws](@id structural_identifiability_conslaws)
120120
Several reaction network models, such as
121121
```@example si3
122122
using Catalyst, Logging, StructuralIdentifiability # hide
123123
rs = @reaction_network begin
124124
(k1,k2), X1 <--> X2
125125
end
126126
```
127-
contain conservation laws (in this case $Γ = X1 + X2$, where $Γ = X1(0) + X2(0)$ is a constant). Because the presence of such conservation laws makes structural identifiability analysis prohibitively computationally expensive (for all but the simplest of cases), these are automatically eliminated by Catalyst (removing one ODE from the resulting ODE system for each conservation law). For the `assess_identifiability` and `assess_local_identifiability` functions, this will be unnoticed by the user. However, for the `find_identifiable_functions` and `make_si_ode` functions, this may result in one, or several, parameters of the form `Γ[i]` (where `i` is an integer) appearing in the produced expressions. These correspond to the conservation law constants and can be found through
128-
```@example si3
129-
conservedequations(rs)
130-
```
131-
E.g. if you run:
132-
```@example si3
133-
find_identifiable_functions(rs; measured_quantities = [:X1, :X2])
134-
```
135-
we see that `Γ[1]` (`= X1(0) + X2(0)`) is detected as an identifiable expression. If we want to disable this feature for any function, we can use the `remove_conserved = false` option:
136-
```@example si3
137-
find_identifiable_functions(rs; measured_quantities = [:X1, :X2], remove_conserved = false)
138-
```
127+
contain [conservation laws](@ref network_analysis_deficiency) (in this case $Γ = X1 + X2$, where $Γ = X1(0) + X2(0)$ is a constant). Because the presence of such conservation laws makes structural identifiability analysis prohibitively computationally expensive (for all but the simplest of cases), these are automatically eliminated by Catalyst. This is handled internally, and should not be noticeable to the user. The exception is the `make_si_ode` function. For each conservation law, its output will have one ODE removed, and instead have a conservation parameter (of the form `Γ[i]`) added to its equations. This feature can be disabled through the `remove_conserved = false` option.
139128

140-
## Systems with exponent parameters
129+
## [Systems with exponent parameters](@id structural_identifiability_exp_params)
141130
Structural identifiability cannot currently be applied to systems with parameters (or species) in exponents. E.g. this
142131
```julia
143132
rn = @reaction_network begin

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

Lines changed: 19 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@ function Catalyst.make_si_ode(rs::ReactionSystem; measured_quantities = [], know
3030
ignore_no_measured_warn = false, remove_conserved = true)
3131
# Creates a MTK ODESystem, and a list of measured quantities (there are equations).
3232
# Gives these to SI to create an SI ode model of its preferred form.
33-
osys, conseqs, _ = make_osys(rs; remove_conserved)
33+
osys, conseqs, _, _ = make_osys(rs; remove_conserved)
3434
measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
3535
conseqs; ignore_no_measured_warn)
3636
return SI.mtk_to_si(osys, measured_quantities)[1]
@@ -67,15 +67,15 @@ function SI.assess_local_identifiability(rs::ReactionSystem, args...;
6767
measured_quantities = [], known_p = [], funcs_to_check = Vector(),
6868
remove_conserved = true, ignore_no_measured_warn = false, kwargs...)
6969
# Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form.
70-
osys, conseqs, vars = make_osys(rs; remove_conserved)
70+
osys, conseqs, consconsts, vars = make_osys(rs; remove_conserved)
7171
measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
7272
conseqs; ignore_no_measured_warn)
7373
funcs_to_check = make_ftc(funcs_to_check, conseqs, vars)
7474

7575
# Computes identifiability and converts it to a easy to read form.
7676
out = SI.assess_local_identifiability(osys, args...; measured_quantities,
7777
funcs_to_check, kwargs...)
78-
return make_output(out, funcs_to_check, reverse.(conseqs))
78+
return make_output(out, funcs_to_check, consconsts)
7979
end
8080

8181
"""
@@ -107,15 +107,15 @@ function SI.assess_identifiability(rs::ReactionSystem, args...;
107107
measured_quantities = [], known_p = [], funcs_to_check = Vector(),
108108
remove_conserved = true, ignore_no_measured_warn = false, kwargs...)
109109
# Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form.
110-
osys, conseqs, vars = make_osys(rs; remove_conserved)
110+
osys, conseqs, consconsts, vars = make_osys(rs; remove_conserved)
111111
measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
112112
conseqs; ignore_no_measured_warn)
113113
funcs_to_check = make_ftc(funcs_to_check, conseqs, vars)
114114

115115
# Computes identifiability and converts it to a easy to read form.
116116
out = SI.assess_identifiability(osys, args...; measured_quantities,
117117
funcs_to_check, kwargs...)
118-
return make_output(out, funcs_to_check, reverse.(conseqs))
118+
return make_output(out, funcs_to_check, consconsts)
119119
end
120120

121121
"""
@@ -147,13 +147,13 @@ function SI.find_identifiable_functions(rs::ReactionSystem, args...;
147147
measured_quantities = [], known_p = [], remove_conserved = true,
148148
ignore_no_measured_warn = false, kwargs...)
149149
# Creates a ODESystem, and list of measured quantities, of SI's preferred form.
150-
osys, conseqs = make_osys(rs; remove_conserved)
150+
osys, conseqs, consconsts, _ = make_osys(rs; remove_conserved)
151151
measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
152152
conseqs; ignore_no_measured_warn)
153153

154154
# Computes identifiable functions and converts it to a easy to read form.
155155
out = SI.find_identifiable_functions(osys, args...; measured_quantities, kwargs...)
156-
return vector_subs(out, reverse.(conseqs))
156+
return vector_subs(out, consconsts)
157157
end
158158

159159
### Helper Functions ###
@@ -172,22 +172,24 @@ function make_osys(rs::ReactionSystem; remove_conserved = true)
172172

173173
# Computes equations for system conservation laws.
174174
# If there are no conserved equations, the `conseqs` variable must still have the `Vector{Pair{Any, Any}}` type.
175-
if !remove_conserved
176-
conseqs = Vector{Pair{Any, Any}}[]
177-
else
175+
if remove_conserved
178176
conseqs = [ceq.lhs => ceq.rhs for ceq in conservedequations(rs)]
177+
consconsts = [cconst.lhs => cconst.rhs for cconst in conservationlaw_constants(rs)]
179178
isempty(conseqs) && (conseqs = Vector{Pair{Any, Any}}[])
179+
isempty(consconsts) && (consconsts = Vector{Pair{Any, Any}}[])
180+
else
181+
conseqs = Vector{Pair{Any, Any}}[]
182+
consconsts = Vector{Pair{Any, Any}}[]
180183
end
181184

182-
return osys, conseqs, vars
185+
return osys, conseqs, consconsts, vars
183186
end
184187

185188
# Creates a list of measured quantities of a form that SI can read.
186189
# Each measured quantity must have a form like:
187190
# `obs_var ~ X` # (Here, `obs_var` is a variable, and X is whatever we can measure).
188-
function make_measured_quantities(
189-
rs::ReactionSystem, measured_quantities::Vector{T}, known_p::Vector{S},
190-
conseqs; ignore_no_measured_warn = false) where {T, S}
191+
function make_measured_quantities(rs::ReactionSystem, measured_quantities::Vector{T},
192+
known_p::Vector{S}, conseqs; ignore_no_measured_warn = false) where {T, S}
191193
# Warning if the user didn't give any measured quantities.
192194
if !ignore_no_measured_warn && isempty(measured_quantities)
193195
@warn "No measured quantity provided to the `measured_quantities` argument, any further identifiability analysis will likely fail. You can disable this warning by setting `ignore_no_measured_warn = true`."
@@ -218,9 +220,9 @@ end
218220
# Processes the outputs to a better form.
219221
# Replaces conservation law equations back in the output (so that e.g. Γ are not displayed).
220222
# Sorts the output according to their input order (defaults to the `[unknowns; parameters]` order).
221-
function make_output(out, funcs_to_check, conseqs)
222-
funcs_to_check = vector_subs(funcs_to_check, conseqs)
223-
out = OrderedDict(zip(vector_subs(keys(out), conseqs), values(out)))
223+
function make_output(out, funcs_to_check, consconsts)
224+
funcs_to_check = vector_subs(funcs_to_check, consconsts)
225+
out = OrderedDict(zip(vector_subs(keys(out), consconsts), values(out)))
224226
sortdict = Dict(ftc => i for (i, ftc) in enumerate(funcs_to_check))
225227
return sort!(out; by = x -> sortdict[x])
226228
end

0 commit comments

Comments
 (0)