SciML / StructuralIdentifiability.jl

Fast and automatic structural identifiability software for ODE systems
https://docs.sciml.ai/StructuralIdentifiability/stable/
MIT License
109 stars 17 forks source link

Cannot (?) use `known_ic` in combination with MTK variables #324

Closed TorkelE closed 2 months ago

TorkelE commented 3 months ago

To designate that some initial conditions are known. E.g. if I do

using StructuralIdentifiability, ModelingToolkit

@parameters r1, r2, c1, c2, beta1, beta2, chi1, chi2
@variables t, x1(t), x2(t), y(t), u(t)

D = Differential(t)
eqs = [
    D(x1) ~ r1 * x1 * (1 - c1 * x1) + beta1 * x1 * x2 / (chi1 + x2) + u,
    D(x2) ~ r2 * x2 * (1 - c2 * x2) + beta2 * x1 * x2 / (chi2 + x1),
]

measured_quantities = [y ~ x1]

ode_mtk = ODESystem(eqs, t, name = :mutualist)

assess_identifiability(ode_mtk, measured_quantities = measured_quantities) # works

assess_identifiability(ode_mtk; measured_quantities = measured_quantities, known_ic = [x1])

I get a

ERROR: MethodError: no method matching assess_identifiability(::ODESystem; measured_quantities::Vector{Equation}, known_ic::Vector{Num})
pogudingleb commented 3 months ago

@TorkelE Thanks for pointing this out, indeed, I think this functionality was not added to the MTK interface. Will do this.

TorkelE commented 3 months ago

Thanks! When you have added it, I will add this to the Catalyst/SI extension and docs as well. Thanks a lot for adding this feature!

pogudingleb commented 2 months ago

@TorkelE I have a PR addressing this issue, here is the current interface. For a model defined as

using StructuralIdentifiability
using ModelingToolkit
using Logging

@parameters r1, r2, c1, c2, beta1, beta2, chi1, chi2
@variables t, x1(t), x2(t), y(t), u(t)

D = Differential(t)

eqs = [
    D(x1) ~ r1 * x1 * (1 - c1 * x1) + beta1 * x1 * x2 / (chi1 + x2) + u,
    D(x2) ~ r2 * x2 * (1 - c2 * x2) + beta2 * x1 * x2 / (chi2 + x1),
]

measured_quantities = [y ~ x1]

ode_mtk = ODESystem(eqs, t, name = :mutualist)

One can provide known_ic argument to assess_identifiability and find_identifiable_functions.

assess_identifiability(ode_mtk, measured_quantities = measured_quantities, known_ic = [x2])

will produce

OrderedCollections.OrderedDict{SymbolicUtils.BasicSymbolic{Real}, Symbol}(
    x1(0) => :globally, 
    x2(0) => :globally, 
    u(0) => :globally, 
    chi1 => :globally, 
    r1 => :globally, 
    beta1 => :globally, 
    c1 => :globally, 
    c2 => :globally, 
    chi2 => :globally, 
    r2 => :globally, 
    beta2 => :globally)

and

find_identifiable_functions(ode_mtk, measured_quantities = measured_quantities, known_ic = [x2])

will produce

Num[beta2, r2, chi2, c2, c1, beta1, r1, chi1, u(0), x2(0), x1(0)]

Note that, if known initial conditions are provided, then the returned identifiability results are valid at t = 0, this is why t = 0 is substituted in the results above.

Let me know if this interface is convenient.

TorkelE commented 2 months ago

Looks good, however, a follow-up question: Does this mean that, if you provide any initial conditions, you are not able to infer identifiability for the full trajectories (i.e. just the other initial conditions)? I think the ideal case would be that you provide the known initial condition and measured quantities, and then receive the identifiability for parameters, variables, and variable initial conditions. If that is not possible/significant amount of work to implement, this is cool as well though (but I think providing all the known information, and receiving all the identifiabilites, would be the ultimate goal).

pogudingleb commented 2 months ago

Does this mean that, if you provide any initial conditions, you are not able to infer identifiability for the full trajectories (i.e. just the other initial conditions)?

Yes, this is what happens now.

I think the ideal case would be that you provide the known initial condition and measured quantities, and then receive the identifiability for parameters, variables, and variable initial conditions.

I completely agree that this would be nice. For the moment, however, I do not have a complete algorithm for this. I will summarize my thoughts and create a separate issue for this improvement.

TorkelE commented 2 months ago

Sounds great, again, thanks a lot for supporting this :)