Structural identifiability analysis fails for composed MTK model about structuralidentifiability.jl HOT 12 CLOSED

@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?

from structuralidentifiability.jl.

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)

from structuralidentifiability.jl.

@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?

from structuralidentifiability.jl.

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

from structuralidentifiability.jl.

@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?

from structuralidentifiability.jl.

@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]

from structuralidentifiability.jl.

@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.

from structuralidentifiability.jl.

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.

from structuralidentifiability.jl.

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

from structuralidentifiability.jl.

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)

from structuralidentifiability.jl.

@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).

from structuralidentifiability.jl.

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 [email protected] ?

from structuralidentifiability.jl.