Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
TorkelE opened this issue Dec 13, 2022 · 2 comments

Comments

@TorkelE
Copy link
Member

TorkelE commented Dec 13, 2022

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
Copy link
Collaborator

Hi @TorkelE , if you use

@variable t y(t)

eq = y~states(osys)[1]

it will work.

@TorkelE
Copy link
Member Author

TorkelE commented Dec 14, 2022

Thanks a lot, yes it works :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants