Symbolic-Numeric Analysis of Parameter Identifiability and Model Stability

The mixture of symbolic computing with numeric computing, which we call symbolic-numeric programming, is one of the central features of the SciML ecosystem. With core aspects like the Symbolics.jl Computer Algebra System and its integration via ModelingToolkit.jl, the SciML ecosystem gracefully mixes analytical symbolic computations with the numerical solver processes to accelerate solvers, give additional information (sparsity, identifiability), automatically fix numerical stability issues, and more.

In this showcase, we will highlight two aspects of symbolic-numeric programming.

1. Automated index reduction of DAEs. While arbitrary differential-algebraic equation systems can be written in DifferentialEquations.jl, not all mathematical formulations of a system are equivalent. Some are numerically difficult to solve, or even require special solvers. Some are easy. Can we recognize which formulations are hard and automatically transform them into the easy ones? Yes.
2. Structural parameter identifiability. When fitting parameters to data, there's always assumptions about whether there is a unique parameter set that achieves such a data fit. But is this actually the case? The structural identifiability tooling allows one to analytically determine whether, in the limit of infinite data on a subset of observables, one could in theory uniquely identify the parameters (global identifiability), identify the parameters up to a discrete set (local identifiability), or whether there's an infinite manifold of solutions to the inverse problem (nonidentifiable).

Let's dig into these two cases!

Automated Index Reduction of DAEs

In many cases one may accidentally write down a DAE that is not easily solvable by numerical methods. In this tutorial, we will walk through an example of a pendulum which accidentally generates an index-3 DAE, and show how to use the modelingtoolkitize to correct the model definition before solving.

Copy-Pastable Example

using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using Plots

function pendulum!(du, u, p, t)
    x, dx, y, dy, T = u
    g, L = p
    du[1] = dx
    du[2] = T * x
    du[3] = dy
    du[4] = T * y - g
    du[5] = x^2 + y^2 - L^2
    return nothing
end
pendulum_fun! = ODEFunction(pendulum!, mass_matrix = Diagonal([1, 1, 1, 1, 0]))
u0 = [1.0, 0, 0, 0, 0]
p = [9.8, 1]
tspan = (0, 10.0)
pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
prob = ODEProblem(pendulum_sys, [], tspan)
sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8)
plot(sol, vars = unknowns(traced_sys))

Explanation

Attempting to Solve the Equation

In this tutorial, we will look at the pendulum system:

\[\begin{aligned} x^\prime &= v_x\\ v_x^\prime &= Tx\\ y^\prime &= v_y\\ v_y^\prime &= Ty - g\\ 0 &= x^2 + y^2 - L^2 \end{aligned}\]

As a good DifferentialEquations.jl user, one would follow the mass matrix DAE tutorial to arrive at code for simulating the model:

using OrdinaryDiffEq, LinearAlgebra
function pendulum!(du, u, p, t)
    x, dx, y, dy, T = u
    g, L = p
    du[1] = dx
    du[2] = T * x
    du[3] = dy
    du[4] = T * y - g
    du[5] = x^2 + y^2 - L^2
end
pendulum_fun! = ODEFunction(pendulum!, mass_matrix = Diagonal([1, 1, 1, 1, 0]))
u0 = [1.0, 0, 0, 0, 0];
p = [9.8, 1];
tspan = (0, 10.0);
pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
solve(pendulum_prob, Rodas5P())

retcode: Unstable
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 5-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 2.9844380280222975e-5
 3.916953199555279e-5
u: 5-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0, 0.0, 0.0]
 [1.0, 0.0, -4.900000000000009e-12, -9.800000000000072e-6, 0.0]
 [1.0, -5.282200000000052e-16, -5.929000000000031e-10, -0.00010780000000000034, 1.2105282010215346e-11]
 [1.0, -1.7340838406781296e-13, -4.364366468121788e-9, -0.000292474926746187, 3.6305781792415833e-8]
 [1.0, -1.6753118817314227e-12, -7.517835960078196e-9, -0.0003838614135564203, 1.4383072941635003e-6]

However, one will quickly be greeted with the unfortunate message:

┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq C:\Users\accou\.julia\packages\OrdinaryDiffEq\yCczp\src\initdt.jl:76
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\accou\.julia\packages\OrdinaryDiffEq\yCczp\src\solve.jl:485
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\DrPil\src\integrator_interface.jl:325

Did you implement the DAE incorrectly? No. Is the solver broken? No.

Understanding DAE Index

It turns out that this is a property of the DAE that we are attempting to solve. This kind of DAE is known as an index-3 DAE. For a complete discussion of DAE index, see this article. Essentially, the issue here is that we have 4 differential variables ($x$, $v_x$, $y$, $v_y$) and one algebraic variable $T$ (which we can know because there is no D(T) term in the equations). An index-1 DAE always satisfies that the Jacobian of the algebraic equations is non-singular. Here, the first 4 equations are differential equations, with the last term the algebraic relationship. However, the partial derivative of x^2 + y^2 - L^2 w.r.t. T is zero, and thus the Jacobian of the algebraic equations is the zero matrix, and thus it's singular. This is a quick way to see whether the DAE is index 1!

The problem with higher order DAEs is that the matrices used in Newton solves are singular or close to singular when applied to such problems. Because of this fact, the nonlinear solvers (or Rosenbrock methods) break down, making them difficult to solve. The classic paper DAEs are not ODEs goes into detail on this and shows that many methods are no longer convergent when index is higher than one. So, it's not necessarily the fault of the solver or the implementation: this is known.

But that's not a satisfying answer, so what do you do about it?

Transforming Higher Order DAEs to Index-1 DAEs

It turns out that higher order DAEs can be transformed into lower order DAEs. If you differentiate the last equation two times and perform a substitution, you can arrive at the following set of equations:

\[\begin{aligned} x^\prime =& v_x \\ v_x^\prime =& x T \\ y^\prime =& v_y \\ v_y^\prime =& y T - g \\ 0 =& 2 \left(v_x^{2} + v_y^{2} + y ( y T - g ) + T x^2 \right) \end{aligned}\]

Note that this is mathematically equivalent to the equation that we had before, but the Jacobian w.r.t. T of the algebraic equation is no longer zero because of the substitution. This means that if you wrote down this version of the model, it will be index-1 and solve correctly! In fact, this is how DAE index is commonly defined: the number of differentiations it takes to transform the DAE into an ODE, where an ODE is an index-0 DAE by substituting out all of the algebraic relationships.

Automating the Index Reduction

However, requiring the user to sit there and work through this process on potentially millions of equations is an unfathomable mental overhead. But, we can avoid this by using methods like the Pantelides algorithm for automatically performing this reduction to index 1. While this requires the ModelingToolkit symbolic form, we use modelingtoolkitize to transform the numerical code into symbolic code, run structural_simplify to simplify the system and lower the index, then transform back to numerical code with ODEProblem, and solve with a numerical solver. Let's try that out:

traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(traced_sys)
prob = ODEProblem(pendulum_sys, Pair[], tspan)
sol = solve(prob, Rodas5P())

using Plots
plot(sol, vars = unknowns(traced_sys))

Note that plotting using unknowns(traced_sys) is done so that any variables which are symbolically eliminated, or any variable reordering done for enhanced parallelism/performance, still show up in the resulting plot and the plot is shown in the same order as the original numerical code.

Note that we can even go a bit further. If we use the ODEProblem constructor, we represent the mass matrix DAE of the index-reduced system, which can be solved via:

traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
prob = ODEProblem(pendulum_sys, Pair[], tspan)
sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8)
plot(sol, vars = unknowns(traced_sys))

And there you go: this has transformed the model from being too hard to solve with implicit DAE solvers, to something that is easily solved.

Parameter Identifiability in ODE Models

Ordinary differential equations are commonly used for modeling real-world processes. The problem of parameter identifiability is one of the key design challenges for mathematical models. A parameter is said to be identifiable if one can recover its value from experimental data. Structural identifiability is a theoretical property of a model that answers this question. In this tutorial, we will show how to use StructuralIdentifiability.jl with ModelingToolkit.jl to assess identifiability of parameters in ODE models. The theory behind StructuralIdentifiability.jl is presented in paper ^[4].

We will start by illustrating local identifiability in which a parameter is known up to finitely many values, and then proceed to determining global identifiability, that is, which parameters can be identified uniquely.

To install StructuralIdentifiability.jl, simply run

using Pkg
Pkg.add("StructuralIdentifiability")

The package has a standalone data structure for ordinary differential equations, but is also compatible with ODESystem type from ModelingToolkit.jl.

Local Identifiability

Input System

We will consider the following model:

\[\begin{cases} \frac{d\,x_4}{d\,t} = - \frac{k_5 x_4}{k_6 + x_4},\\ \frac{d\,x_5}{d\,t} = \frac{k_5 x_4}{k_6 + x_4} - \frac{k_7 x_5}{(k_8 + x_5 + x_6)},\\ \frac{d\,x_6}{d\,t} = \frac{k_7 x_5}{(k_8 + x_5 + x_6)} - \frac{k_9 x_6 (k_{10} - x_6) }{k_{10}},\\ \frac{d\,x_7}{d\,t} = \frac{k_9 x_6 (k_{10} - x_6)}{ k_{10}},\\ y_1 = x_4,\\ y_2 = x_5\end{cases}\]

This model describes the biohydrogenation^[1] process^[2] with unknown initial conditions.

Using the ODESystem object

To define the ode system in Julia, we use ModelingToolkit.jl.

We first define the parameters, variables, differential equations and the output equations.

using StructuralIdentifiability, ModelingToolkit

# define parameters and variables
@variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t)
@parameters k5 k6 k7 k8 k9 k10
D = Differential(t)

# define equations
eqs = [
    D(x4) ~ -k5 * x4 / (k6 + x4),
    D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6),
    D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10,
    D(x7) ~ k9 * x6 * (k10 - x6) / k10
]

# define the output functions (quantities that can be measured)
measured_quantities = [y1 ~ x4, y2 ~ x5]

# define the system
de = ODESystem(eqs, t, name = :Biohydrogenation)

After that, we are ready to check the system for local identifiability:

# query local identifiability
# we pass the ode-system
local_id_all = assess_local_identifiability(de, measured_quantities = measured_quantities,
    p = 0.99)
# [ Info: Preproccessing `ModelingToolkit.ODESystem` object
# 6-element Vector{Bool}:
#  1
#  1
#  1
#  1
#  1
#  1

We can see that all unknowns (except $x_7$) and all parameters are locally identifiable with probability 0.99.

Let's try to check specific parameters and their combinations

to_check = [k5, k7, k10 / k9, k5 + k6]
local_id_some = assess_local_identifiability(de, measured_quantities = measured_quantities,
    funcs_to_check = to_check, p = 0.99)
# 4-element Vector{Bool}:
#  1
#  1
#  1
#  1

Notice 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$

Global Identifiability

In this part tutorial, let us cover an example problem of querying the ODE for globally identifiable parameters.

Input System

Let us consider the following four-dimensional model with two outputs:

\[\begin{cases} x_1'(t) = -b x_1(t) + \frac{1 }{ c + x_4(t)},\\ x_2'(t) = \alpha x_1(t) - \beta x_2(t),\\ x_3'(t) = \gamma x_2(t) - \delta x_3(t),\\ x_4'(t) = \sigma x_4(t) \frac{(\gamma x_2(t) - \delta x_3(t))}{ x_3(t)},\\ y_1(t) = x_1(t) + x_2(t) y_2(t) = x_2(t) \end{cases}\]

We will run a global identifiability check on this enzyme dynamics^[3] model. We will use the default settings: the probability of correctness will be p=0.99 and we are interested in identifiability of all possible parameters.

Global identifiability needs information about local identifiability first, but the function we chose here will take care of that extra step for us.

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.

using StructuralIdentifiability, ModelingToolkit
@parameters b c a beta g delta sigma
@variables t x1(t) x2(t) x3(t) x4(t) y1(t) y2(t)
D = Differential(t)

eqs = [
    D(x1) ~ -b * x1 + 1 / (c + x4),
    D(x2) ~ a * x1 - beta * x2,
    D(x3) ~ g * x2 - delta * x3,
    D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3
]

measured_quantities = [y1 ~ x1 + x2, y2 ~ x2]

ode = ODESystem(eqs, t, name = :GoodwinOsc)

@time global_id = assess_identifiability(ode, measured_quantities = measured_quantities)
# 30.672594 seconds (100.97 M allocations: 6.219 GiB, 3.15% gc time, 0.01% compilation time)
# Dict{Num, Symbol} with 7 entries:
#   a     => :globally
#   b     => :globally
#   beta  => :globally
#   c     => :globally
#   sigma => :globally
#   g     => :nonidentifiable
#   delta => :globally

We can see that only parameters a, g are unidentifiable, and everything else can be uniquely recovered.

Let us consider the same system but with two inputs, and we will find out identifiability with probability 0.9 for parameters c and b:

using StructuralIdentifiability, ModelingToolkit
@parameters b c a beta g delta sigma
@variables t x1(t) x2(t) x3(t) x4(t) y(t) u1(t) [input = true] u2(t) [input = true]
D = Differential(t)

eqs = [
    D(x1) ~ -b * x1 + 1 / (c + x4),
    D(x2) ~ a * x1 - beta * x2 - u1,
    D(x3) ~ g * x2 - delta * x3 + u2,
    D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3
]
measured_quantities = [y1 ~ x1 + x2, y2 ~ x2]

# check only 2 parameters
to_check = [b, c]

ode = ODESystem(eqs, t, name = :GoodwinOsc)

global_id = assess_identifiability(ode, measured_quantities = measured_quantities,
    funcs_to_check = to_check, p = 0.9)
# Dict{Num, Symbol} with 2 entries:
#   b => :globally
#   c => :globally

Both parameters b, c are globally identifiable with probability 0.9 in this case.