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

Do structural identifiability for ModelingToolkit ODESystem (where observables are not specified in system) #133

Closed TorkelE closed 1 year ago

TorkelE commented 1 year ago

I am trying to do identifiability analysis for a chemical reaction network model. I specify the model via Catalyst, and then convert to ODESystem. But that means that I cannot declare the y as a part of the model, and I am a bit unsure how to add them afterwards.

This was my attempt (toy example, to model is simplified and not that interesting, but gives the same error):

# Fectch packages.
using Catalyst, DifferentialEquations, StructuralIdentifiability

# Create model and convert to ODESystem.
rn = @reaction_network begin
    (k1,k2), X <--> Y
end k1 k2
osys = convert(ODESystem,rn)

# Create equation for observables (X only).
@variables y
eq = y ~ states(osys)[1]

# Do identifiability analysis.
assess_identifiability(osys; measured_quantities=[eq])

However, this gives a:

ERROR: MethodError: no method matching arguments(::Sym{Real, Base.ImmutableDict{DataType, Any}})
Closest candidates are:
  arguments(::Symbolics.CallWithMetadata) at ~/.julia/packages/Symbolics/FGTCH/src/variable.jl:241
  arguments(::SymbolicUtils.Mul) at ~/.julia/packages/SymbolicUtils/qulQp/src/types.jl:808
  arguments(::Term) at ~/.julia/packages/SymbolicUtils/qulQp/src/types.jl:317
  ...
Stacktrace:
 [1] PreprocessODE(de::ODESystem, measured_quantities::Vector{Equation})
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/egNv7/src/ODE.jl:383
 [2] assess_identifiability(ode::ODESystem; measured_quantities::Vector{Equation}, funcs_to_check::Vector{Any}, p::Float64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/egNv7/src/StructuralIdentifiability.jl:166
 [3] top-level scope
   @ ~/Desktop/Research Projects/Jakobs-Problem/Scripts/test_identifiability.jl:45
iliailmer commented 1 year ago

Hi @TorkelE , if you use

@variable t y(t)

eq = y~states(osys)[1]

it will work.

TorkelE commented 1 year ago

Thanks a lot, yes it works :)