You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
During parameter fitting, parameter values are inferred from data. Parameter identifiability refers to whether inferring parameter values for a given model is mathematically feasible. Ideally, parameter fitting should always be accompanied with an identifiability analysis of the problem.
3
3
4
-
Identifiability can be divided into *structural* and *practical* identifiability[^1]. Structural identifiability considers only the mathematical model, and which parameters are and are not inherently identifiable due to model structure. Practical identifiability also considers the available data, and determines what system quantities can be inferred from it. In the idealized case of an infinite amount of noise-less data, practical identifiability becomes identical to structural identifiability. Generally, structural identifiability is assessed before parameters are fitted, while practical identifiability is assessed afterwards.
4
+
Identifiability can be divided into *structural* and *practical* identifiability[^1]. Structural identifiability considers only the mathematical model, and which parameters are and are not inherently identifiable due to model structure. Practical identifiability also considers the available data, and determines what system quantities can be inferred from it. In the idealised case of an infinite amount of non-noisy data, practical identifiability converges to structural identifiability. Generally, structural identifiability is assessed before parameters are fitted, while practical identifiability is assessed afterwards.
5
5
6
6
Structural identifiability (which is what this tutorial considers) can be illustrated by the following differential equation:
7
7
${dx \over dt} = p1*p2*x(t)$
8
8
where, however much data is collected on $x$, it is impossible to determine the distinct values of $p1$ and $p2$. Hence, these parameters are these are non-identifiable (however, their product, $p1*p2$, *is* identifiable).
9
9
10
-
Catalyst contains a special extension for carrying out structural identifiability analysis using the [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) package. This enables StructuralIdentifiability's `assess_identifiability`, `assess_local_identifiability`, and `find_identifiable_functions` functions to be called directly on Catalyst `ReactionSystem`s. It also implements specialised routines to make these more efficient when applied to reaction network models. How to use these functions is described in the following tutorial, with [StructuralIdentifiability providing a more extensive documentation](https://docs.sciml.ai/StructuralIdentifiability/stable/). If you use this in your research, please [cite the StructuralIdentifiability.jl](@ref structural_identifiability_citation) and [Catalyst.jl](@ref catalyst_citation) publications.
10
+
Catalyst contains a special extension for carrying out structural identifiability analysis using the [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) package. This enables StructuralIdentifiability's `assess_identifiability`, `assess_local_identifiability`, and `find_identifiable_functions` functions to be called directly on Catalyst `ReactionSystem`s. It also implements specialised routines to make these more efficient when applied to reaction network models. How to use these functions is described in the following tutorial, with [StructuralIdentifiability providing a more extensive documentation](https://docs.sciml.ai/StructuralIdentifiability/stable/).
11
11
12
-
Structural identifiability can be divided into *local* and *global* identifiability. If a model quantity is locally identifiable, it means that its true value can be determined down to a finite-number of possible options. This also means that there is some limited region around its true value, where the true value is the only possible value. Globally identifiable quantities' values, on the other hand, can be uniquely determined. Again, while identifiability can be confirmed structurally for a model, it does not necessarily mean that they are practically identifiable for some given data.
12
+
Structural identifiability can be divided into *local* and *global* identifiability. If a model quantity is locally identifiable, it means that its true value can be determined down to a finite-number of possible options. This also means that there is some limited region around the quantity's true value where this true value is the only possible value (and hence, within this region, the quantity is fully identifiable). Globally identifiable quantities' values, on the other hand, can be uniquely determined. Again, while identifiability can be confirmed structurally for a model, it does not necessarily mean that they are practically identifiable for some given data.
13
13
14
14
Generally, there are three types of quantities for which identifiability can be assessed.
15
15
- Parameters (e.g. $p1$ and $p2$).
16
16
- Full variable trajectories (e.g. $x(t)$).
17
17
- Variable initial conditions (e.g. $x(0)$).
18
18
19
-
StructutralIdentifiability currently assesses identifiability for the first two (however, if $x(t)$ is identifiable, $x(0)$ will be as well).
19
+
StructuralIdentifiability currently assesses identifiability for the first two only (however, if $x(t)$ is identifiable, then $x(0)$ will be as well).
20
20
21
21
## Global identifiability analysis
22
22
23
23
### Basic example
24
24
25
25
Local identifiability can be assessed using the `assess_identifiability` function. For each model quantity (parameters and variables), it will assess whether they are:
26
-
-globally identifiable.
27
-
-locally identifiable.
28
-
-unidentifiable.
26
+
-Globally identifiable.
27
+
-Locally identifiable.
28
+
-Unidentifiable.
29
29
30
-
To it, we provide our `ReactionSystem` model and a list of quantities that we are able to measure. Here, we consider a Goodwind oscillator (a simple 3-component model, where the three species $M$, $E$, and $P$ are produced and degraded, which may exhibit oscillations)[^2]. Let us say that we are able to measure the concentration of $M$, we then provide designate this using the `measured_quantities` argument. We can now assess identifiability in the following way:
30
+
To it, we provide our `ReactionSystem` model and a list of quantities that we are able to measure. Here, we consider a Goodwind oscillator (a simple 3-component model, where the three species $M$, $E$, and $P$ are produced and degraded, which may exhibit oscillations)[^2]. Let us say that we are able to measure the concentration of $M$, we then designate this using the `measured_quantities` argument. We can now assess identifiability in the following way:
31
31
```example si1
32
32
using Catalyst, Logging, StructuralIdentifiability
33
33
goodwind_oscillator = @reaction_network begin
@@ -37,15 +37,13 @@ goodwind_oscillator = @reaction_network begin
From the output, we find that `E(t)`, `pₑ`, and `pₚ` (the trajectory of $E$, and the production rates of $E$ and $P$, respectively) are non-identifiable. Next, `dₑ` and `dₚ` (the degradation rates of $E$ and $P$, respectively) are locally identifiable. Finally, `P(t)`, `M(t)`, `pₘ`, and `dₘ` (the trajectories of `P` and `M`, and the production and degradation rate of `M`, respectively) are all globally identifiable. We note that we also imported the Logging package, and provided the `loglevel=Logging.Error` input argument. StructuralIdentifiability functions generally provide a large number of output messages. Hence, we will use this argument (which requires the Logging package) throughout this tutorial to decrease the amount of printed text.
40
+
From the output, we find that `E(t)`, `pₑ`, and `pₚ` (the trajectory of $E$, and the production rates of $E$ and $P$, respectively) are non-identifiable. Next, `dₑ` and `dₚ` (the degradation rates of $E$ and $P$, respectively) are locally identifiable. Finally, `P(t)`, `M(t)`, `pₘ`, and `dₘ` (the trajectories of `P` and `M`, and the production and degradation rate of `M`, respectively) are all globally identifiable. We note that we also imported the Logging.jl package, and provided the `loglevel=Logging.Error` input argument. StructuralIdentifiability functions generally provide a large number of output messages. Hence, we will use this argument (which requires the Logging package) throughout this tutorial to decrease the amount of printed text.
41
41
42
-
43
-
Next, we can assess identifiability in the case where we can measure all three species concentrations:
42
+
Next, we also assess identifiability in the case where we can measure all three species concentrations:
in which case all initial conditions and parameters become identifiable.
48
-
46
+
in which case all species trajectories and parameters become identifiable.
49
47
50
48
### Indicating known parameters
51
49
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ₑ$):
@@ -57,7 +55,7 @@ Not only does this turn the previously non-identifiable `pₑ` (globally) identi
57
55
To, in a similar manner, indicate that certain initial conditions are known is a work in progress. Hopefully the feature should be an available in the near future.
58
56
59
57
### Providing non-trivial measured quantities
60
-
Sometimes, we are not actually measuring species species, but rather some combinations of species (or possibly parameters). Here, any algebraic expression can be used in `measured_quantities`. Species and parameters have to first be `@unpack`'ed from the model to be used in such expressions. 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$:
58
+
Sometimes, we are not actually measuring species species, but rather some combinations of species (or possibly parameters). To account for this, `measured_quantities` accepts any algebraic expressions (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$:
61
59
```example si2
62
60
using Catalyst, StructuralIdentifiability # hide
63
61
enzyme_activation = @reaction_network begin
@@ -81,12 +79,12 @@ nothing # hide
81
79
```
82
80
83
81
### Probability of correctness
84
-
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 `p` (by default set to `0.99`, that is, at least a $99%$ chance of correctness). We can e.g. increase the bound through:
82
+
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 `p` (by default set to `0.99`, that is, at least a $99\%$ chance of correctness). We can e.g. increase the bound through:
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 `p` increases the certainty of correctness, it will also increase the time required to assess identifiability.
87
+
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 `p` increases the certainty of correctness, it will also increase the time required to assess identifiability.
90
88
91
89
## Local identifiability analysis
92
90
Local identifiability can be assessed through the `assess_local_identifiability` function. While this is already determined by `assess_identifiability`, assessing local identifiability only have the advantage that it is easier to compute. Hence, there might be models where global identifiability analysis fails (or takes prohibitively long time), where instead `assess_local_identifiability` can be used. This functions 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:
@@ -100,12 +98,12 @@ Finally, StructuralIdentifiability provides the `find_identifiable_functions` fu
Again, these results are consistent with those of produced by `assess_identifiability`. There, `pₑ` and `pₚ` where found to be globally identifiable. Here, they correspond directly to identifiable expressions. The remaining four parameters (`pₘ`, `dₘ`, `dₑ`, and `dₚ`) occurs as part of more complicated composite expressions.
101
+
Again, these results are consistent with those produced by `assess_identifiability`. There, `pₑ` and `pₚ` where found to be globally identifiable. Here, they correspond directly to identifiable expressions. The remaining four parameters (`pₘ`, `dₘ`, `dₑ`, and `dₚ`) occurs as part of more complicated composite expressions.
104
102
105
-
The `find_identifiable_functions` functions tries to simplify its output functions to create nice expression. 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 runtimes).
103
+
The `find_identifiable_functions` 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 runtimes).
106
104
107
105
## Creating StructuralIdentifiability compatible ODE models from Catalyst `ReactionSystem`s
108
-
While the functionality described above covers the vast majority of analysis as 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. Similarly to the previous functions, it takes a `ReactionSystem` and 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:
106
+
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. Similarly to the previous functions, it takes a `ReactionSystem` and 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:
we see that `Γ[1]` (`= X1 + X2`) is detected as an identifiable expression. If we want to disable this feature for any function, we can use the `remove_conserved = false` option:
133
+
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:
0 commit comments