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

ModelingToolkit Workflow only returns Parameter Identifiability #200

Closed lucakollmer closed 1 year ago

lucakollmer commented 1 year ago

Hello,

I am writing to ask if it is intended functionality for assess_identifiability and assess_local_identifiability to only return the structural identifiability of parameters in the ODEs defined using the ModelingToolkit workflow, as opposed to the @ODEmodel macro workflow, which additionally returns the observability of the states.

Analysis of the compartmental dynamic model given in using_modeling_toolkit.md using both methods yields the same set of identifiable parameters, but the @ODEmodel method also returns which states are observable. Is there any way to get the observability of states out of the ModelingToolkit method? It would be very useful to be able to use the ModelingToolkit method and also analyse state observability. Thank you.

ModelingToolkit method:



@parameters b q N_inv k r alpha g1 g2
@variables t S(t) E(t) A(t) I(t) J(t) C(t) y1(t) y2(t)

D = Differential(t)

eqs = [
    D(S) ~ -b * S * (I + J + q * A) * N_inv,
    D(E) ~ b * S * (I + J + q * A) * N_inv - k * E,
    D(A) ~ k * (1 - r) * E - g1 * A,
    D(I) ~ k * r * E - (alpha + g1) * I,
    D(J) ~ alpha * I - g2 * J,
    D(C) ~ alpha * I,
]

ode = ODESystem(eqs, t, name = :SEIAJRCmodel)

measured_quantities = [y1 ~ C, y2 ~ N_inv]
@time assess_local_identifiability(ode, measured_quantities = measured_quantities)

Dict{SymbolicUtils.BasicSymbolic{Real}, Bool} with 8 entries:
  g1    => 1
  q     => 0
  r     => 0
  k     => 1
  b     => 1
  g2    => 1
  alpha => 1
  N_inv => 1

@ODEmodel method:

ode = @ODEmodel(
    S'(t) = -b * S(t) * (I(t) + J(t) + q * A(t)) * N_inv,
    E'(t) = b * S(t) * (I(t) + J(t) + q * A(t)) * N_inv - k * E(t),
    A'(t) = k * (1 - r) * E(t) - g1 * A(t),
    I'(t) = k * r * E(t) - (alpha + g1) * I(t),
    J'(t) = alpha * I(t) - g2 * J(t),
    C'(t) = alpha * I(t),
    y1(t) = C(t),
    y2(t) = N_inv
)

assess_local_identifiability(ode)

Dict{Nemo.fmpq_mpoly, Bool} with 14 entries:
  k     => 1
  g1    => 1
  g2    => 1
  q     => 0
  alpha => 1
  b     => 1
  N_inv => 1
  E     => 0
  J     => 1
  I     => 1
  r     => 0
  C     => 1
  S     => 0
  A     => 0
pogudingleb commented 1 year ago

@lucakollmer Thanks for bringing this up! This is not intended, and it has been fixed in the latest release 0.4.10. Which version do you use?

lucakollmer commented 1 year ago

@pogudingleb Thank you for your reply. I am on v0.4.9 which I installed a few weeks ago using Pkg.add("StructuralIdentifiability"). I will update using the GitHub repo.