Structural Analysis for algebraic difference equations (DiscreteSystem)

[RodrigoZepeda]

Some biological models are represented via difference equations of the form: $x(t + 1) = f\big( x(t), \theta \big)$. According to Chris Rackauckas on Twitter, an extension of the identifiability functions could be made for those systems using the DiscreteSystem from ModelingToolkit.

Some references on the subject are:

Nõmm, Sven, and C. H. Moog. "Identifiability of discrete-time nonlinear systems." IFAC Proceedings Volumes 37.13 (2004): 333-338.

Nõmm, Sven, and Claude H. Moog. "Further results on identifiability of discrete-time nonlinear systems." Automatica 68 (2016): 69-74.

[pogudingleb]

@RodrigoZepeda
Thanks for bringing this up! Indeed, at least for local identifiability this does not seem to be complicated.
If you have some models you would like to check, let me know - I can try to prioritise this.

[RodrigoZepeda]

@pogudingleb
I'm working with State Space Models and variations upon the SIR model from epidemiology. I cannot copy the exact example. However, here are some similar toy models:

using ModelingToolkit, DifferentialEquations

#SIR Model (should not be identifiable)
@parameters α β
@variables t S(t) I(t) R(t)
D = Difference(t; dt=1.0)

eqs = [D(S) ~ S - β*S*I,
       D(I) ~ I + β*S*I - α*I,
       D(R) ~ R + α*I]       
@named sir = DiscreteSystem(eqs)

#Linear model (A can be a matrix)
#Identifiability depends upon the matrix
@parameters A , b
@variables t x(t)
D = Difference(t; dt=1.0)

eqs = [D(x) ~ A*x + b]       
@named linear = DiscreteSystem(eqs)

#An example that can be identified whenever x(0) != 0
@parameters θ
@variables t x(t)
D = Difference(t; dt=1.0)

eqs = [D(x) ~ θ*x^3]

@named eqs = DiscreteSystem(eqs)

#An example for local identifiability 
#(different values of θ for R^+ and R^-)
@parameters θ
@variables t x(t)
D = Difference(t; dt=1.0)

eqs = [D(x) ~ θ^2 + x]

@named eqs = DiscreteSystem(eqs)

#An example that cannot be identified 
@parameters θ β
@variables t x(t) y(t)
D = Difference(t; dt=1.0)

eqs = [D(x) ~ x + y,
       D(y) ~ θ + β]

@named eqs = DiscreteSystem(eqs)

I am happy to generate more toy models to test the algorithms and help any way I can within my limited expertise on the subject.

[pogudingleb]

@RodrigoZepeda
Thanks for the examples !
I have implemented local identifiability check for discrete time systems, it is available in the current GitHub version (and can be installed by Pkg.add("https://github.com/SciML/StructuralIdentifiability.jl")).
Here is an example of computation:

using ModelingToolkit
using StructuralIdentifiability

@parameters α β
@variables t S(t) I(t) R(t) y(t)
D = Difference(t; dt = 1.0)

eqs = [D(S) ~ S - β * S * I, D(I) ~ I + β * S * I - α * I, D(R) ~ R + α * I]
@named sir = DiscreteSystem(eqs)

println(assess_local_identifiability(sir; measured_quantities = [y ~ I]))

which will print a dictionary from parameters and states to their identifiability (under assumption that the I compartment is observed).
Any feedback will be appreciated !

[RodrigoZepeda]

Woow. I'll take a look during the week for feedback and trying additional systems. Thanks!

[pogudingleb]

Quick clarification: the algorithm check identifiability for generic initial conditions and exciting enough inputs, if the algorithm reports identifiability there may be still some singular solutions from which the parameters/states of interest cannot be identified.

[RodrigoZepeda]

There seems to be an issue when the parameter is a denominator as the following example shows:

using ModelingToolkit, StructuralIdentifiability

@parameters α β
@variables t x(t)
D = Difference(t; dt=1.0)

eqs = [D(x) ~ α + (1/β)*x] #Doesn't work

@named eqs = DiscreteSystem(eqs)

assess_local_identifiability(eqs; measured_quantities = [x ~ x, β ~ β])

which throws an error:

ERROR: Element not invertible
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] inv(a::Nemo.QQFieldElem)
   @ Nemo ~/.julia/packages/Nemo/xM4S1/src/flint/fmpq.jl:487
 [3] eval_at_dict(rational::AbstractAlgebra.Generic.Frac{Nemo.QQMPolyRingElem}, d::Dict{Nemo.QQMPolyRingElem, Nemo.QQFieldElem})
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ExdSQ/src/util.jl:29
 [4] sequence_solution(dds::ODE{Nemo.QQMPolyRingElem}, param_values::Dict{Nemo.QQMPolyRingElem, Nemo.QQFieldElem}, initial_conditions::Dict{Nemo.QQMPolyRingElem, Nemo.QQFieldElem}, input_values::Dict{Nemo.QQMPolyRingElem, Vector{Nemo.QQFieldElem}}, num_terms::Int64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ExdSQ/src/discrete.jl:30
 [5] differentiate_sequence_solution(dds::ODE{Nemo.QQMPolyRingElem}, params::Dict{Nemo.QQMPolyRingElem, Nemo.QQFieldElem}, ic::Dict{Nemo.QQMPolyRingElem, Nemo.QQFieldElem}, inputs::Dict{Nemo.QQMPolyRingElem, Vector{Nemo.QQFieldElem}}, num_terms::Int64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ExdSQ/src/discrete.jl:76
 [6] differentiate_sequence_output(dds::ODE{Nemo.QQMPolyRingElem}, params::Dict{Nemo.QQMPolyRingElem, Nemo.QQFieldElem}, ic::Dict{Nemo.QQMPolyRingElem, Nemo.QQFieldElem}, inputs::Dict{Nemo.QQMPolyRingElem, Vector{Nemo.QQFieldElem}}, num_terms::Int64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ExdSQ/src/discrete.jl:130
 [7] _assess_local_identifiability_discrete(dds::ODE{Nemo.QQMPolyRingElem}, funcs_to_check::Vector{Nemo.QQMPolyRingElem}, known_ic::Vector{Union{}}, p::Float64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ExdSQ/src/discrete.jl:235
 [8] assess_local_identifiability(dds::DiscreteSystem; measured_quantities::Vector{Equation}, funcs_to_check::Vector{Array}, known_ic::Vector{Array}, p::Float64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ExdSQ/src/discrete.jl:341
 [9] top-level scope
   @ REPL[7]:1

Of course, as a temporal solution I can always change variables and do $\gamma = 1/\beta$ and that works.

[pogudingleb]

@RodrigoZepeda
Thanks for testing!
The thing is that the outputs are expected to have new names, so the following will work fine

using ModelingToolkit, StructuralIdentifiability

@parameters α β
@variables t x(t) y1(t) y2(t)
D = Difference(t; dt=1.0)

eqs = [D(x) ~ α + (1/β)*x] #Doesn't work

@named eqs = DiscreteSystem(eqs)

assess_local_identifiability(eqs; measured_quantities = [y1 ~ x, y2 ~ β])

Let me know if this makes sense.

[RodrigoZepeda]

Ohhhhh ok this solves many of the weird isues I was having

[pogudingleb]

@RodrigoZepeda
I thought about this issue and this is indeed not completely intuitive that one has to introduce extra names for the measured quantities. So I have added (to the current GitHub version) an option not to name these variables. The above example can be now also run as (the change in the last line)

using ModelingToolkit
using StructuralIdentifiability

@parameters α β
@variables t S(t) I(t) R(t) y(t)
D = Difference(t; dt = 1.0)

eqs = [D(S) ~ S - β * S * I, D(I) ~ I + β * S * I - α * I, D(R) ~ R + α * I]
@named sir = DiscreteSystem(eqs)

println(assess_local_identifiability(sir; measured_quantities = [I]))

[RodrigoZepeda]

Awesome! I'll try and add a tutorial on this to the documentation during the following week (June 20-something)

[pogudingleb]

@RodrigoZepeda
I made a simple tutorial for this functionality, let me know what you think. If you have any interesting example(s) to add, it will be greatly appreciated

[RodrigoZepeda]

It looks great! Thank you. I'm out of the office for a while but will definitley come back here in ~ a month. Last time I completely forgot about my comment but feel free to ping me if I don't come back.

[pogudingleb]

It looks great! Thank you. I'm out of the office for a while but will definitley come back here in ~ a month. Last time I completely forgot about my comment but feel free to ping me if I don't come back.

Thanks! Note that I have misinterpreted Difference from MTK in my first implementation, so make sure use use version 0.4.16 or older.