Advice on getting global identifiabiltiy to work for a relatively small (?) system

[TorkelE]

I have a system for which I want to assess identifiability. While local identifiability works, global does not:

using StructuralIdentifiability, ModelingToolkit

@parameters k1 k2 k3
@variables t LPd(t) ArBr(t) LPdArBr(t) RB(t) LPdArR(t) BBr(t) ArR(t)
D = Differential(t)

eqs = [
    D(LPd) ~ k3*LPdArR - k1*ArBr*LPd,
    D(ArBr) ~ -k1*ArBr*LPd,
    D(LPdArBr) ~ k1*ArBr*LPd - k2*LPdArBr*RB,
    D(RB) ~ -k2*LPdArBr*RB,
    D(LPdArR) ~ k2*LPdArBr*RB - k3*LPdArR,
    D(BBr) ~ k2*LPdArBr*RB,
    D(ArR) ~ k3*LPdArR
]

@named ode = ODESystem(eqs, t)

measured_quantities = [ArR]
@time assess_local_identifiability(ode, measured_quantities=measured_quantities)
@time global_id = assess_identifiability(ode, measured_quantities=measured_quantities)

it gets stuck on:

[ Info: Preproccessing `ModelingToolkit.AbstractTimeDependentSystem` object
[ Info: Assessing local identifiability
[ Info: Local identifiability assessed in 0.117122279 seconds
[ Info: Assessing global identifiability
[ Info: Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`
[ Info: Computing IO-equations

I understand that global identifiability might not be feasible for larger system, but I figured this one might be small enough? I tried a simplified version, and that was really fast. I have left it for like 2 hours, with no result. If I put it to an HPC for like 3 days, might that help, or at this point, are you probably screwed?

[pogudingleb]

@TorkelE Thanks for the nice example! The model is indeed not so complicated but it turn out to have some tricky structural property which complicates the computation. I would suggest to use a reduction with first integrals as described in #63 (we have a code for automatising this, hopefully will be finalised and merged over the summer). More precisely, one can observe that LPd - ArBr - ArR = const and LPdArR - RB - ArR = const. We denote these constants by C and C2 and use them to eliminate LPd and LPdArR from the model, so we get an equivalent model

using StructuralIdentifiability, ModelingToolkit

@parameters k1 k2 k3 C C2
@variables t LPd(t) ArBr(t) LPdArBr(t) RB(t) LPdArR(t) BBr(t) ArR(t)
D = Differential(t)

eqs = [
    D(ArBr) ~ -k1*ArBr*(C + ArBr + ArR),
    D(LPdArBr) ~ k1*ArBr*(C + ArBr + ArR) - k2*LPdArBr*RB,
    D(RB) ~ -k2*(C2 + RB + ArR)*RB,
    D(ArR) ~ k3*(C2 + RB + ArR)
]

@named ode = ODESystem(eqs, t)
measured_quantities = [ArR]
@time assess_local_identifiability(ode, measured_quantities=measured_quantities)
@time global_id = assess_identifiability(ode, measured_quantities=measured_quantities)

This takes 30 sec on my laptop. Let me know if you have any questions/comments.

[TorkelE]

Thanks :)