Structural identifiability analysis fails for composed MTK model

[jonasmac16]

Hi, First, thank for this cool package.

I am trying to compose model based on several submodels written in MTK. To do I have input variable in each submodel which allows me to define input and output from other states variables in the composed model.

However, this seems to interfere with the identifiability analysis. I have written a MWE with a Lotka-Volterra model exemplifying the issue. The non-composed models works fine but the composed version fails (see error below). Is this expected? There seem to similar issues posted here.

using DifferentialEquations
using CairoMakie
using ModelingToolkit
using StructuralIdentifiability

function lotka_volterra!(du, u, p, t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α * x - β * x * y ## rabbits
    du[2] = dy = -δ * y + γ * x * y ## wolves
end

# Initial condition
u0 = [1.0, 2.0]

# Simulation interval
tspan = (0.0, 10.0)

# LV equation parameter. p = [α, β, δ, γ]
p = [2.5, 3.0, 7.0, 2.0]

# Setup the ODE problem, then solve
prob = ODEProblem(lotka_volterra!, u0, tspan, p)
sol = solve(prob)
datasol = solve(prob, saveat = 1)
data = Array(datasol)

f = Figure()
ax = Axis(f[1,1])
tplot = collect(tspan[1]:0.01:tspan[2])
lines!(ax, tplot, [sol(j)[1] for j in tplot])
lines!(ax, tplot, [sol(j)[2] for j in tplot])
scatter!(ax, collect(0:1:10), data[1,:])
scatter!(ax, collect(0:1:10), data[2,:])
f

@variables t
function rabbits_creator(;name)
    ps = @parameters α=1.5
    vars = @variables x(t)=1.0 z(t)=0.0 [input = true]
    D = Differential(t)
    equs = [
        D(x) ~ α *x + z
    ]

    ODESystem(equs, t, vars, ps; name = name)
end

function wolves_creator(;name)
    ps = @parameters δ = 3.0
    vars = @variables y(t)=1.0 q(t)=0.0 [input = true]
    D = Differential(t)
    equs = [
        D(y) ~ -δ * y + q
    ]

    ODESystem(equs, t, vars, ps; name = name)
end

function lotka_volterra_creator(;name)
    @named wolves = wolves_creator()
    @named rabbits = rabbits_creator()

    ps = @parameters β=1.0 γ=1.0
    D = Differential(t)

    eqs = [ 
        rabbits.z ~  - β *  wolves.y * rabbits.x, 
        wolves.q ~ γ *  wolves.y *rabbits.x]

    compose(ODESystem(eqs, t, [],ps; name= name), wolves, rabbits)
end

@named ltk_mtk = lotka_volterra_creator()
simp_ltk_mtk = structural_simplify(ltk_mtk)
equations(simp_ltk_mtk)
parameters(simp_ltk_mtk)
@unpack wolves₊δ, rabbits₊α , β, γ,  wolves₊y, rabbits₊x = simp_ltk_mtk

u0_sim = [
    rabbits₊x => u0[1],
    wolves₊y => u0[2]
]

pars_sim = [
    rabbits₊α => p[1],
    β =>  p[2],
    wolves₊δ => p[3],
    γ => p[4]
]

sol_mtk = solve(ODEProblem(simp_ltk_mtk, u0_sim, tspan, pars_sim))

prob = ODEProblem(lotka_volterra!, u0, tspan, p)
sol = solve(prob)

f = Figure()
ax = Axis(f[1,1])
tplot = collect(tspan[1]:0.01:tspan[2])
lines!(ax, tplot, [sol_mtk(j)[1] for j in tplot])
lines!(ax, tplot, [sol_mtk(j)[2] for j in tplot])
lines!(ax, tplot, [sol(j)[1] for j in tplot], linestyle=:dash)
lines!(ax, tplot, [sol(j)[2] for j in tplot], linestyle=:dash)
f

ModelingToolkit.observed(simp_ltk_mtk)
original_sys = modelingtoolkitize(prob)

@unpack wolves₊y, rabbits₊x = simp_ltk_mtk
@unpack x₁, x₂ = original_sys
@variables y1(t) y2(t)

This works

measured_quantities = [y1 ~ x₁, y2 ~ x₂]
assess_identifiability(original_sys, measured_quantities = measured_quantities)

but this fails with ERROR: MethodError: no method matching +(::Nemo.fmpq_mpoly, ::Term{Real, Base.ImmutableDict{DataType, Any}})

measured_quantities = [y1 ~ rabbits₊x, y2 ~ wolves₊y]
assess_identifiability(simp_ltk_mtk, measured_quantities = measured_quantities)

[pogudingleb]

@jonasmac16 I first tried with version 0.4.5 of StructuralIdentifiability and got the same error, this error comes from a recent update of SymbolicUtils. The new version, 0.4.7, does not have this error but has a different one, looking into this. Which version are you using?

[jonasmac16]

Thanks for looking into this. I am using StructuralIdentifiability v0.4.5.

[pogudingleb]

@jonasmac16 I have studied the error I am getting with the more recent version (StructuralIdentifiability v0.4.7 or higher). It seems to me that for some reason structural_simplify does not work properly: for equations(ltd_mtk), I get

4-element Vector{Equation}:
 rabbits₊z(t) ~ -β*rabbits₊x(t)*wolves₊y(t)
 wolves₊q(t) ~ γ*rabbits₊x(t)*wolves₊y(t)
 Differential(t)(wolves₊y(t)) ~ wolves₊q(t) - wolves₊δ*wolves₊y(t)
 Differential(t)(rabbits₊x(t)) ~ rabbits₊α*rabbits₊x(t) + rabbits₊z(t)

but then equations(simp_ltk_mtk) gives

2-element Vector{Equation}:
 Differential(t)(wolves₊y(t)) ~ wolves₊q(t) - wolves₊δ*wolves₊y(t)
 Differential(t)(rabbits₊x(t)) ~ rabbits₊α*rabbits₊x(t) + rabbits₊z(t)

so it looks like the first two equations in the former model are just deleted and not used for any substitution. This is not what you expected, right?

[jonasmac16]

@pogudingleb I think since z(t) and q(t) are inputs they are moved to observed. So if you do observed(simp_ltk_mtk) ypu get

2-element Vector{Equation}:
 wolves₊q(t) ~ γ*rabbits₊x(t)*wolves₊y(t)
 rabbits₊z(t) ~ -β*rabbits₊x(t)*wolves₊y(t)

and running assess_identifiability on the non-simplified model also throws an error:

assess_identifiability(ltk_mtk, measured_quantities = measured_quantities)
[ Info: Preproccessing `ModelingToolkit.ODESystem` object
ERROR: BoundsError: attempt to access 2-element Vector{Any} at index [3]

[pogudingleb]

@jonasmac16 Thank you for explaining, it seems that I had a too simplistic understanding of the internal structure of MTK models. Will work further on this.

[pogudingleb]

@jonasmac16 Thanks again for bringing up the issue and for explanation. The current github version should work for your example (works for me at least). Could you check?

[jonasmac16]

Just updated to the github version and it works. Thank you so much for looking at this so quickly and solving it. I have more complicated models which have a similar logic to them and I will see if I encounter more edge cases and report back.

[pogudingleb]

Then I close the issue for the time being. Feel free to reopen if you encounter something.

[wsphillips]

Hi, I had the same issue with observed variables throwing errors and the latest tagged version of StructuralIdentifiability.jl does not appear to have fixed it.

If you run the following example from the root directory of Conductor.jl you will hit the error:

using Conductor, StructuralIdentifiability, ModelingToolkit
import Unitful: mV

# Run from root of Conductor.jl
include(joinpath("demo", "prinz_kinetics.jl"));

@variables y(t) [unit=mV]

channels = [NaV(100mS / cm^2),
    CaT(0mS / cm^2),
    CaS(4mS / cm^2),
    KA(0mS / cm^2),
    KCa(15mS / cm^2),
    Kdr(50mS / cm^2),
    H(0.02mS / cm^2),
    leak(0.03mS / cm^2)];

dynamics = HodgkinHuxley(channels, gradients);

@named neuron = CompartmentSystem(Vₘ, dynamics;
                                  geometry = geo,
                                  extensions = [calcium_conversion]);

odesys = ODESystem(neuron) # converts to ODESystem and runs structural_simplify
measured_quantities = [y ~ Vₘ]

# Errors because of a variable moved to observed on latest v0.4.9
global_id = assess_identifiability(odesys, measured_quantities = measured_quantities)

[pogudingleb]

@wsphillips Thanks for the example! There is indeed a bug causing this behaviour, I will fix it today/tomorrow. However, I see that your model involves exponential functions, and these are unfortunately not handled right now out of the box (but you can still make a variable transformation, see #144, let me know if you would like my help on this).

[wsphillips]

Thank you for taking another look at this issue.

Regarding exponential functions, they are common in the specification of these models, but could certainly be substituted (and often are for performance reasons). I would indeed need some help to understand what specifically needs to happen. Could you ping me on JuliaLang Slack or email wsphil@mit.edu ?

[pogudingleb]

I have fixed the issue in the current GitHub version, it will be a part of 0.4.10 release. @wsphillips I will follow up by email on exponential functions (there are also max functions somewhere apparently)