Creating ODE System

Most of the algorithms in StructuralIdentifiability.jl take as input system of ordinary differential equations (ODEs) in the state space form, that is:

\[\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}\]

which involves

• a vector $\mathbf{x}(t)$ of the state variables of the system,

• a vector $\mathbf{u}(t)$ of external inputs,

• a vector $\mathbf{p}$ of scalar parameters,

• a vector $\mathbf{y}(t)$ of outputs (i.e., observations),

• and vectors of rational functions $\mathbf{f}$ and $\mathbf{g}$ (for discussion of the non-rational case, see this issue).

In the standard setup, inputs and outputs are assumed to be known, and the goal is to assess identifiability of parameters and/or states from the input-output data. In the case of states, this property is also called observability.

There are two ways to define such a system to be processed using StructuralIdentifiability.jl. We will demonstrate them using the following example system (Wright's population model of two mutualist species with control^[1]):

\[\begin{cases} x_1'(t) = r_1 x_1(t)(1 - c_1 x_1(t)) + \frac{\beta_1 x_1(t)x_2(t)}{\chi_1 + x_2(t)} + u(t),\\ x_2'(t) = r_2 x_2(t)(1 - c_2 x_2(t)) + \frac{\beta_2 x_1(t)x_2(t)}{\chi_2 + x_1(t)},\\ y(t) = x_1(t). \end{cases}\]

Defining the model using @ODEmodel macro

One way to define the model is to use the @ODEmodel macro provided by the StructuralIdentifiability.jl package.

using StructuralIdentifiability

ode = @ODEmodel(
    x1'(t) =
        r1 * x1(t) * (1 - c1 * x1(t)) + beta1 * x1(t) * x2(t) / (chi1 + x2(t)) + u(t),
    x2'(t) = r2 * x2(t) * (1 - c2 * x2(t)) + beta2 * x1(t) * x2(t) / (chi2 + x1(t)),
    y(t) = x1(t)
)

x1'(t) = (-x1(t)^2*x2(t)*c1*r1 - x1(t)^2*c1*chi1*r1 + x1(t)*x2(t)*beta1 + x1(t)*x2(t)*r1 + x1(t)*chi1*r1 + x2(t)*u(t) + u(t)*chi1)//(x2(t) + chi1)
x2'(t) = (-x1(t)*x2(t)^2*c2*r2 + x1(t)*x2(t)*beta2 + x1(t)*x2(t)*r2 - x2(t)^2*c2*chi2*r2 + x2(t)*chi2*r2)//(x1(t) + chi2)
y(t) = x1(t)

Then one can, for example, assess identifiability of the parameters and states by

assess_identifiability(ode)

OrderedCollections.OrderedDict{Any, Symbol} with 10 entries:
  x1(t) => :globally
  x2(t) => :nonidentifiable
  beta1 => :globally
  beta2 => :globally
  c1    => :globally
  c2    => :nonidentifiable
  chi1  => :nonidentifiable
  chi2  => :globally
  r1    => :globally
  r2    => :globally

Defining using ModelingToolkit

StructuralIdentifiability has an extension ModelingToolkitSIExt which allows to use System from ModelingToolkit to describe a model. The extension is loaded automatically once ModelingToolkit is loaded via using ModelingToolkit. In this case, one should encode the equations for the states as System and specify the outputs separately. In order to do this, we first introduce all functions and scalars:

using StructuralIdentifiability, ModelingToolkit

@independent_variables t
@parameters r1, r2, c1, c2, beta1, beta2, chi1, chi2
@variables x1(t), x2(t), y(t), u(t)

D = Differential(t)

And then defined the system and separately the outputs:

eqs = [
    D(x1) ~ r1 * x1 * (1 - c1 * x1) + beta1 * x1 * x2 / (chi1 + x2) + u,
    D(x2) ~ r2 * x2 * (1 - c2 * x2) + beta2 * x1 * x2 / (chi2 + x1),
]

measured_quantities = [y ~ x1]

ode_mtk = System(eqs, t, name = :mutualist)

\[ \begin{align} \frac{\mathrm{d} \mathtt{x1}\left( t \right)}{\mathrm{d}t} &= u\left( t \right) + \frac{\mathtt{beta1} \mathtt{x1}\left( t \right) \mathtt{x2}\left( t \right)}{\mathtt{chi1} + \mathtt{x2}\left( t \right)} + \left( 1 - \mathtt{c1} \mathtt{x1}\left( t \right) \right) \mathtt{r1} \mathtt{x1}\left( t \right) \\ \frac{\mathrm{d} \mathtt{x2}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{beta2} \mathtt{x1}\left( t \right) \mathtt{x2}\left( t \right)}{\mathtt{chi2} + \mathtt{x1}\left( t \right)} + \left( 1 - \mathtt{c2} \mathtt{x2}\left( t \right) \right) \mathtt{r2} \mathtt{x2}\left( t \right) \end{align} \]

Then, for example, the identifiability of parameters and states can be assessed as follows:

assess_identifiability(ode_mtk, measured_quantities = measured_quantities)

OrderedCollections.OrderedDict{SymbolicUtils.BasicSymbolic{Real}, Symbol} with 10 entries:
  x1(t) => :globally
  x2(t) => :nonidentifiable
  c2    => :nonidentifiable
  chi1  => :nonidentifiable
  r1    => :globally
  chi2  => :globally
  beta1 => :globally
  c1    => :globally
  r2    => :globally
  beta2 => :globally