Recall that we consider ODE models in the state-space form
\[\begin{cases}
+\mathbf{x}'(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}, \mathbf{u}(t)),\\
+\mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t), \mathbf{p}, \mathbf{u(t)}),
+\end{cases}\]
where $\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively, and $\mathbf{p}$ are scalar parameters. We will call that a parameter or a states (or a function of them) is identifiable if its value can be recovered from time series for inputs and outputs. Typically, two types of identifiability are distinguished
Note that in the case of states, term observability it typically used. We will use identifiability for both states and parameters for brevity and uniformity. While the notion of global identifiability is much more precise, assessing local identifiability is typically much faster, and can be performed for the models whose global identifiability analysis is out of reach.
We consider the Lotka-Volterra model:
\[\begin{cases}
+x_1'(t) = a x_1(t) - b x_1(t) x_2(t) + u(t),\\
+x_2'(t) = -c x_2(t) + d x_1(t) x_2(t),\\
+y(t) = x_1(t)
+\end{cases}\]
The local identifiability of all parameters and states in this model can be assessed as follows
using StructuralIdentifiability
+
+ode = @ODEmodel(
+ x1'(t) = a * x1(t) - b * x1(t) * x2(t) + u(t),
+ x2'(t) = -c * x2(t) + d * x1(t) * x2(t),
+ y(t) = x1(t)
+)
+
+assess_local_identifiability(ode)
OrderedCollections.OrderedDict{Any, Bool} with 6 entries:
+ x1(t) => 1
+ x2(t) => 0
+ a => 1
+ b => 0
+ c => 1
+ d => 1
We see that $x_1(t)$ is locally identifiable (no surprises, this is an output), $a, c,$ and $d$ are identifiable as well. On the other hand, $x_2(t)$ and $b$ are nonidentifiable. This can be explained by the following scaling symmetry
\[x_2(t) \to \lambda x_2(t), \quad b \to \frac{b}{\lambda}\]
which preserves input and output for every nonzero $\lambda$. The algorithm behind this call is the one due to Sedoglavic.
Function assess_local_identifiability has several optional parameters
funcs_to_check a list of specific functions of parameters and states to check identifiability for (see an example below). If not provided, the identifiability is assessed for all parameters and states.
prob_threshold (default $0.99$, i.e. 99%) is the probability of correctness. The algorithm can, in theory, produce wrong result, but the probability that it is correct is guaranteed to be at least prob_threshold. However, the probability bounds we use are quite conservative, so the actual probability of correctness is likely to be much higher.
type (default :SE). By default, the algorithm checks the standard single-experiment identifiability. If one sets type = :ME, then the algorithm checks multi-experiment identifiability, that is, identifiability from several experiments with independent initial conditions (the algorithm from is used).
loglevel (default Logging.Info). The minimal level of logging messages to be displayed. Available options: Logging.Debug, Logging.Info, Logging.Warn, and Logging.Error.
Note that the scaling symmetry given above suggests that $b x_2(t)$ may in fact be identifiable. This can be checked using funcs_to_check parameter:
assess_local_identifiability(ode, funcs_to_check = [b * x2])
OrderedCollections.OrderedDict{Any, Bool} with 1 entry:
+ x2(t)*b => 1
Indeed!
One can obtain more refined information about a model using assess_identifiability function. We will showcase it using the Goodwin oscillator model .
using StructuralIdentifiability
+
+ode = @ODEmodel(
+ x1'(t) = -b * x1(t) + 1 / (c + x4(t)),
+ x2'(t) = alpha * x1(t) - beta * x2(t),
+ x3'(t) = gama * x2(t) - delta * x3(t),
+ x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t),
+ y(t) = x1(t)
+)
+
+assess_identifiability(ode)
OrderedCollections.OrderedDict{Any, Symbol} with 11 entries:
+ x1(t) => :globally
+ x2(t) => :nonidentifiable
+ x3(t) => :nonidentifiable
+ x4(t) => :globally
+ alpha => :nonidentifiable
+ b => :globally
+ beta => :locally
+ c => :globally
+ delta => :locally
+ gama => :nonidentifiable
+ sigma => :globally
As a result, each parameter/state is assigned one of the labels :globally (globally identifiable), :locally (locally but not globally identifiable), or :nonidentifiable (not identifiable, even locally). The algorithm behind this computation follows .
Similarly to assess_local_identifiability, this function has optional parameters:
funcs_to_check a list of specific functions of parameters and states to check identifiability for (see an example below). If not provided, the identifiability is assessed for all parameters and states. Note that the computations for states may be more involved than for the parameters, so one may want to call the function with funcs_to_check = ode.parameters if the call assess_identifiability(ode) takes too long.
prob_threshold (default $0.99$, i.e. 99%) is the probability of correctness. Same story as above: the probability estimates are very conservative, so the actual error probability is much lower than 1%. Also, currently, the probability of correctness does not include the probability of correctness of the modular reconstruction for Groebner bases. This probability is ensured by an additional check modulo a large prime, and can be neglected for practical purposes.
loglevel (default Logging.Info). The minimal level of logging messages to be displayed. Available options: Logging.Debug, Logging.Info, Logging.Warn, and Logging.Error.
Using funcs_to_check parameter, one can further inverstigate the nature of the lack of identifiability in the model at hand. For example, for the Goodwin oscillator, we can check if beta + delta and beta * delta are identifiable:
assess_identifiability(ode, funcs_to_check = [beta + delta, beta * delta])
OrderedCollections.OrderedDict{Any, Symbol} with 2 entries:
+ beta + delta => :globally
+ beta*delta => :globally
And we see that they indeed are. This means, in particular, that the reason why beta and delta are not identifiable is because their values can be exchanged. One may wonder how could we guess these functions beta + delta, beta * delta. In fact, they can be just computed using find_identifiable_functions function as we will explain in the next tutorial. Stay tuned!
An experimental feature allows to provide an additional keyword argument known_ic to inidcate functions of states and parameters for which the initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be generic). In this case, the identifiability will be assessed for parameters and all the initial conditions or for the initial conditions of funcs_to_check. Let us add an assumption that the initial conditions x2(0) and x3(0) are known:
assess_identifiability(ode, known_ic = [x2, x3])
OrderedCollections.OrderedDict{Nemo.QQMPolyRingElem, Symbol} with 11 entries:
+ x1(0) => :globally
+ x2(0) => :globally
+ x3(0) => :globally
+ x4(0) => :globally
+ alpha => :locally
+ b => :globally
+ beta => :locally
+ c => :globally
+ delta => :locally
+ gama => :locally
+ sigma => :globally
And we see that now alpha and gama become locally identifiable.