sciml / structuralidentifiability.jl Goto Github PK

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

Hi.

When I use StructuralIdentifiability to define an ODE-model, and then check its identifiability, these big "Info:" blocks are printed after either action. I'm not that interested in them, and they take up a lot of space. Is there a way to suppress or hide these/to have them not appear?
Any help would appreciated.

Here is an example of what I'm talking about:
Input:

ode_fhn = @ODEmodel(
    v'(t) = v - v*v*v - w + I_ext,
    w'(t) = (v - a - b*w)/tau,
    y(t) = v(t)
)

Output:

┌ Info: Summary of the model:
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/ODE.jl:303
┌ Info: State variables: w, v
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/ODE.jl:304
┌ Info: Parameters: I_ext, a, b, tau
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/ODE.jl:305
┌ Info: Inputs: 
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/ODE.jl:306
┌ Info: Outputs: y
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/ODE.jl:307

v'(t) = I_ext - w(t) - v(t)^3 + v(t)
w'(t) = (-w(t)*b - a + v(t))//tau
y(t) = v(t)

Input:

print(assess_identifiability(ode_fhn))

Output (only the first line here is relevant for me):

Dict{Any, Symbol}(I_ext => :nonidentifiable, a => :nonidentifiable, b => :globally, tau => :globally)

┌ Info: Assessing local identifiability
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/StructuralIdentifiability.jl:99
┌ Info: Local identifiability assessed in 0.0034953 seconds
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/StructuralIdentifiability.jl:101
┌ Info: Assessing global identifiability
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/StructuralIdentifiability.jl:117
┌ Info: Computing IO-equations
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:200
┌ Info: Computed in 0.0015555 seconds
│   :ioeq_time = ioeq_time
│   ioeq_time = 0.0015555
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:203
┌ Info: Computing Wronskians
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:206
┌ Info: Computed in 0.0016289 seconds
│   :wrnsk_time = wrnsk_time
│   wrnsk_time = 0.0016289
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:208
┌ Info: Dimensions of the wronskians [4]
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:212
┌ Info: Ranks of the Wronskians computed in 1.38e-5 seconds
│   :rank_time = rank_time
│   rank_times = 1.38e-5
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:217
┌ Info: Assessing global identifiability using the coefficients of the io-equations
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:224
┌ Info: Computed in 0.0017633 seconds
│   :check_time = check_time
│   check_time = 0.0017633
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/global_identifiability.jl:226
┌ Info: Global identifiability assessed in 0.0053466 seconds
└ @ StructuralIdentifiability /home/jovyan/.julia/packages/StructuralIdentifiability/bRtkP/src/StructuralIdentifiability.jl:119

Hello people,

The package is throwing errors when I use it in Julia version 1.7.3. But surprisingly, in Julia version 1.6.3, the same piece of code runs smoothly without any errors. Even the example given in the Git page's readme is throwing an error!

Output from v 1.6.3:

[ Info: Global identifiability assessed in 15.339882808 seconds
Dict{Any, Symbol} with 4 entries:
  a12 => :locally
  a21 => :globally
  a01 => :locally
  b   => :nonidentifiable

Output from v 1.7.3:

[ Info: Computing Wronskians
ERROR: MethodError: Cannot `convert` an object of type
  Nemo.gfp_abs_series to an object of type
  AbstractAlgebra.Generic.AbsSeries
Closest candidates are:
  convert(::Type{T1}, ::CxxWrap.CxxWrapCore.SmartPointer{DerivedT}) where {BaseT, T1<:BaseT, DerivedT<:BaseT} at ~/.julia/packages/CxxWrap/OcN1Z/src/CxxWrap.jl:276
  convert(::Type{T}, ::T) where T at /Apps/conda/share/julia/base/essentials.jl:218
Stacktrace:
 [1] setindex!(h::Dict{AbstractAlgebra.MPolyElem, AbstractAlgebra.Generic.AbsSeries}, v0::Nemo.gfp_abs_series, key::Nemo.gfp_mpoly)
   @ Base ./dict.jl:381
 [2] wronskian(io_equations::Dict{Nemo.fmpq_mpoly, Nemo.fmpq_mpoly}, ode::ODE{Nemo.fmpq_mpoly})
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/wronskian.jl:223
 [3] macro expansion
   @ ./timing.jl:299 [inlined]
 [4] assess_global_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Any}, p::Float64; var_change::Symbol, gb_method::Symbol)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:207
 [5] assess_global_identifiability
   @ ~/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:200 [inlined]
 [6] macro expansion
   @ ./timing.jl:299 [inlined]
 [7] assess_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Nemo.fmpq_mpoly}, p::Float64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:110
 [8] assess_identifiability (repeats 2 times)
   @ ~/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:68 [inlined]
 [9] top-level scope
   @ REPL[3]:1

Also surprisingly, the version of the package is v0.3.9 in Julia 1.6.3 while it is v0.1.0 in Julia 1.7.3, and I am unable to update to higher version in the 1.7 Julia.

Please help me out!

The following code

using StructuralIdentifiability, OrdinaryDiffEq, ModelingToolkit
function f(du,u,p,t)
  b,c,a,beta,g,delta,sigma = p
  du[1] = -b * u[1] + 1/(c + u[4])
  du[2] = a * u[1] - beta * u[2]
  du[3] = g * u[2] - delta * u[3]
  du[4] = sigma * u[4] * (g * u[2] - delta * u[3])/u[3]
end
prob = ODEProblem(f,ones(4),(0.0,1.0),ones(7))
sys = modelingtoolkitize(prob)
st = states(sys)
params = parameters(sys)
@variables t y(t) x[1:length(st)](t) 
@parameters mu[1:length(params)]

states_map = Dict(s=>x[i] for (i, s) ∈ enumerate(st))
param_map = Dict(p=>mu[i] for (i, p) ∈ enumerate(params))

eqs = substitute(substitute(equations(sys), states_map), param_map)

sys = ODESystem(eqs, t, name=:Name)
measured_quantities = [y~states_map[us[1]]]
@time global_id = assess_identifiability(sys, measured_quantities=measured_quantities)

Produces this error at the very end

ERROR: ArgumentError: The polynomial contains a variable mu[2] not present in the new ring 9733449230697594235070667*mu[2]^2*mu[1]^5*mu[4]*mu[7]^2 + 9733449230697594235070667*mu[2]^2*mu[1]^5*mu[6]*mu[7]^2 - 48667246153487971175353335*mu[2]^2*mu[1]^4*mu[4]*mu[6]*mu[7]^3 - 12977932307596792313427556*mu[2]^2*mu[1]^4*mu[4]*mu[6]*mu[7]^2 + 520606133531898376008304868289618073659180126368695604752803844382800307368205982265758604268378386364287597134108550309252380302258675*mu[7]^2

This is due to a Singular issue, as in #43. I could not find a workaround so far. If it's possible to automate this part

@variables t y(t) x[1:length(st)](t) 
@parameters mu[1:length(params)]

avoiding indexed notation, i.e. x[1].

In addition to knowing the ODE system and time series data for some outputs, one frequently also knows the initial conditions for some of the states. It would be great to use this information to clarify the identifiability analysis.
For an algorithm based on input-output equations, an interesting discussion of this problem can be found in this paper. To the best of my understanding, the paper does not give a complete algorithm but some ideas can be of interest.

Some interesting specific things we could do:

• fixed initial conditions can be used when verifying the rank of the Wronskian. If it is still of full rank, then identifiability with these fixed initial conditions is at least as good as without them (that is, we do not loose identifiability because of hitting a singular point).
• As suggested in the paper mentioned above, if the system is controllable from the fixed initial conditions, then again the result is at least as good.
• What would be interesting to see is if we could get a different result of elimination by making some initial conditions fixed in the GenericPointGenerator used in the elimination process...

Nemo has now multivariate factorization, so there is no need to do conversion to Singular for this purpose anymore in fast_factor.

Sometimes (see #128) the equations have the nested structure which, after unfolding into a rational function, becomes prohibitively large. For these cases it would be cool to have the local identifiability test working directly on expressions via symbols automatic differentiation, I believe that the original Sedoglavic algorithm we use can be extended to do this.

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?

It isn't an issue, but more of a functionality. Is there a way to impose constraint like, the bounds for parameters, or the bounds for the state-variables, in the identifiability analysis?

I wish to create a StructuralIdentifaibiltiy ODE from a ModelingToolkit model. However, it fails if the measured quantities are a combination of known parameters and measured quantities:

# Fetch packages.
using ModelingToolkit, StructuralIdentifiability

# Create model
@parameters p1 p2
@variables t x1(t) x2(t) y(t)
D = Differential(t)
eqs = [D(x1) ~ p1* x2 - x1,
       D(x2) ~ p2 - p1* x2]
@named osys = ODESystem(eqs)

# Create StructuralIdentifiability ODE.
known_parameter = [p1]
measured_quantities = [y ~ x1]
known_stuff = [Num.(measured_quantities);]
si_ode = preprocess_ode(osys, known_stuff)

# This, providing either only known parameters or equations work:
si_ode = preprocess_ode(osys, known_parameter)
si_ode = preprocess_ode(osys, measured_quantities)

Is there a way to create a system with this kind of information?

At the moment all the algorithms in the package require input models to have rational right-hand side (that is, the rhs must be a quotient of two polynomials). One way to lift this limitation is to transform the system into a rational or even polynomial one. It is know that this can be done by introducing new variables (detailed algorithm with implementation in Python).

However, there is an important subtlety. Consider an example:

x'(t) = a * e^x(t)
y(t) = x(t)

The standard idea would be to introduce a new variable z(t) := e^x(t), this will recast the system into

x'(t) = a * z(t)
z'(t) = a * z(t)^2
y(t) = x(t)

In the latter model, parameter a will be non identifiable. However, if we add a different new variable w'(t) := a * e^x(t), the system will become

x'(t) = w(t)
w'(t) = w(t)^2
y(t) = x(t)

Although a has disappeared from the system, knowing x(0) and w(0) is sufficient to find it. Since these initial conditions are identifiable (and they are), a is in fact identifiable as well.

The difference between the two transformations above is that the former increases the number of degrees of freedom in the system and the latter does not. Because of this, the former yield an incorrect result. Therefore, the overall goal would be to have a procedure for transforming the model into a rational one while keeping the number of degrees of freedom the same. For practical examples of "careful" transformation, see sections A.2 and A.3 in this paper.

One plausible approach is:

1. Implement any reasonable procedure for transforming a system into a rational one (e.g. adding all non rational terms as new variables until we are done). This may increase the number of degrees of freedom.
2. Reduce back the number of degrees of freedom by applying an algorithm for finding scaling transformation from this paper.

This approach will work nicely for the example above. Implementing any of the two steps above separately would be a nice and welcome contribution on its own.

UPD: A simpler procedure can be established for parameters only #190

I am trying to write an extension for Catalyst/StructuralIdentifiability that would allow user to input by Catalyst defined ReactionSystems to assess_identifiability and assess_local_identifiability.

Basically, ReactionSystems can be converted to ODESystems, so I would simply use this conversion. However, generally, Catalyst does not have observables directly defined, so I wanted the user to input these using measured_quantities. Typically, they would just be species of the system, so the user would do something like

rs = @reaction_network begin
    (p,d), 0 <--> X
    (k1,k2), X <--> Y
end
assess_identifiability(rs; measured_quantities=[:Y])

and the appropriate conversion would be done under the hood. Also, we could do some automatic removal of conserved quantities to e.g. deal with cases like: #189

However, I was trying to prototype things following something like: #133 but fails. Here is my minimal example:

# Fectch packages.
using Catalyst, StructuralIdentifiability

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

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

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

which gives

[ Info: Preproccessing `ModelingToolkit.AbstractTimeDependentSystem` object
[ Info: Assessing local identifiability
[ Info: Local identifiability assessed in 0.002107838 seconds
[ Info: Assessing global identifiability
[ Info: Functions to check involve states
[ Info: Computing IO-equations
┌ Info: Computed in 0.001850315 seconds
│   :ioeq_time = :ioeq_time
└   ioeq_time = 0.001850315
[ Info: Computing Wronskians
┌ Info: Computed in 0.001188107 seconds
│   :wrnsk_time = :wrnsk_time
└   wrnsk_time = 0.001188107
[ Info: Dimensions of the Wronskians [2]
┌ Info: Ranks of the Wronskians computed in 1.259e-5 seconds
│   :rank_time = :rank_time
└   rank_times = 1.259e-5
[ Info: Simplifying identifiable functions
ERROR: type VSCodeLogger has no field min_level
Stacktrace:
  [1] getproperty
    @ Base ./Base.jl:37 [inlined]
  [2] groebner_loglevel()
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/utils.jl:28
  [3] discover_shape!(state::ParamPunPam.GroebnerState{…}, modular::ParamPunPam.ModularTracker{…}; η::Int64)
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:181
  [4] discover_shape!
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:163 [inlined]
  [5] _paramgb(blackbox::StructuralIdentifiability.IdealMQS{…}, ordering::Groebner.InputOrdering, up_to_degree::Tuple{…}, estimate_degrees::Bool, assess_correctness::Bool, rational_interpolator::Symbol, polynomial_interpolator::Symbol)
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:123
  [6] paramgb(blackbox::StructuralIdentifiability.IdealMQS{…}; kwargs::@Kwargs{…})
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:97
  [7] paramgb
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:58 [inlined]
  [8] macro expansion
    @ StructuralIdentifiability ./timing.jl:395 [inlined]
  [9] groebner_basis_coeffs(rff::StructuralIdentifiability.RationalFunctionField{…}; seed::Int64, ordering::Groebner.InputOrdering, up_to_degree::Tuple{…}, rational_interpolator::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/RationalFunctionFields/RationalFunctionField.jl:328
 [10] simplified_generating_set(rff::StructuralIdentifiability.RationalFunctionField{…}; p::Float64, seed::Int64, simplify::Symbol, check_variables::Bool, rational_interpolator::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/RationalFunctionFields/RationalFunctionField.jl:529
 [11] macro expansion
    @ StructuralIdentifiability ./timing.jl:395 [inlined]
 [12] initial_identifiable_functions(ode::ODE{…}; p::Float64, known::Vector{…}, with_states::Bool, var_change_policy::Symbol, rational_interpolator::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:139
 [13] initial_identifiable_functions
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:83 [inlined]
 [14] check_identifiability(ode::ODE{…}, funcs_to_check::Vector{…}; known::Vector{…}, p::Float64, var_change_policy::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:181
 [15] check_identifiability
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:163 [inlined]
 [16] assess_global_identifiability(ode::ODE{…}, funcs_to_check::Vector{…}, known::Vector{…}, p::Float64; var_change::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:280
 [17] assess_global_identifiability
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:268 [inlined]
 [18] macro expansion
    @ StructuralIdentifiability ./timing.jl:395 [inlined]
 [19] assess_identifiability(ode::ODE{Nemo.QQMPolyRingElem}; funcs_to_check::Vector{Nemo.QQMPolyRingElem}, p::Float64)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/StructuralIdentifiability.jl:132
 [20] assess_identifiability(ode::ODESystem; measured_quantities::Vector{Equation}, funcs_to_check::Vector{Any}, p::Float64)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/StructuralIdentifiability.jl:184
 [21] top-level scope
    @ ~/Desktop/Julia Playground/Catalyst playground/playground.jl:115
Some type information was truncated. Use `show(err)` to see complete types.

I am using StructuralIdentifiability.jl package to check the Identifiability of ODE by fixing some model parameters. The model includes seven variable states and 22 parameters, of which 12 are fixed. (I defined fixed parameters by param_values as Dict.)
I have googled it and found nothing more except using StructuralIdentifiability.set_parameter_values(ode, param_values) function, but when I use it and define param_values by param_values = Dict(ηψᵥ=>19.03,γc=>1.31,ψₛ=>410e-3,ξ=>510e-9,
μₛ=>0.14,μᵢ=>0.15,ψᵣ=>0.02,αᵣ=>0.001,μᵣ=>0.02,
κ=>0.05,μₐ=>0.02,νₐ=>0.07
)
, face the following error: Any suggestion or other solutions??
MethodError: no method matching set_parameter_values(::ODE{Nemo.fmpq_mpoly}, ::Dict{Nemo.fmpq_mpoly, Float64})
Closest candidates are:
set_parameter_values(::ODE{P}, ::Dict{P, T}) where {T<:AbstractAlgebra.FieldElem, P<:AbstractAlgebra.MPolyElem{T}} at ~/.julia/packages/StructuralIdentifiability/ooOlO/src/ODE.jl:62

Main idea

Many models of interest (especially, in life sciences) have linear first integrals (for example, conservation of mass, etc.). If there is such a first integral a_1 * x_1 + ... + a_n * x_n, where x_1, ..., x_n are the states and a_1, ...., a_n are constants (maybe involving parameters), one can do the following:

1. introduce a new parameter C = a_1 * x_1 + ... + a_n * x_n
2. use x_1 = 1 / a_1 * (C - a_2 * x_2 - ... - a_n * x_n) to eliminate x_1 from the system

As a result, we reduce the number of states by one at the expense of having one more parameter. Since the typical bottleneck of the algorithm is differential elimination, this will very likely speed up the computation.

Motivating example

As a motivating/illustrating example, we consider this CRN model but with only one output:

ode = @ODEmodel(
    x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t),
    x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t),
    x3'(t) = k3 * x4(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
    x4'(t) = k1 * x1(t) * x2(t) - k2 * x4(t) - k3 * x4(t),
    x5'(t) = k4 * x6(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
    x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * x3(t) * x5(t),
    y1(t) = x3(t)
)

This model is quite hard, it does not finish in an hour (and I remember trying it earlier, I think it just took all my memory in a couple of hours).

One can observe that x5(t) + x6(t) is constant. Thus, we can introduce C1 = x5(t) + x6(t) and eliminate x5(t):

ode = @ODEmodel(
    x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t),
    x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t),
    x3'(t) = k3 * x4(t) + k5 * x6(t) - k6 * x3(t) * (C1 - x6(t)),
    x4'(t) = k1 * x1(t) * x2(t) - k2 * x4(t) - k3 * x4(t),
    x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * x3(t) * (C1 - x6(t)),
    y1(t) = x3(t)
)

Doing the same with first integrals C2 = x2(t) + x4(t) and C3 = x1(t) - x2(t) + x3(t) + x6(t), we arrive at the following model with only three states:

ode = @ODEmodel(
    x1'(t) = -k1 * x1(t) * x2(t) + k2 * (C2 - x2(t)) + k4 * x6(t),
    x2'(t) = -k1 * x1(t) * x2(t) + k2 * (C2 - x2(t)) + k3 * (C2 - x2(t)),
    x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * (C3 - x1(t) + x2(t) - x6(t)) * (C1 - x6(t)),
    y1(t) = (C3 - x1(t) + x2(t) - x6(t))
)

Now global identifiability is successfully analyzed in 60 seconds (and computing the input-output equation takes only 15 seconds)

What to do

The substantial speedup above can be achieved automatically if we can

1. find all linear first integrals
2. perform a reduction using this integrals
3. add this as a preprocessing in assess_global_identifiability

If a have just a expression x(t), which is not an ODE, should I caluculate the derivative of it in order to be able to test structural identifiability or os there a way to implemet just a funtion?

Hi,

I am just curious and testing your library.

Now I have one question:
Which type of system does your library support? Only polynomial ODEs?
I am asking because I am getting the following error, when parsing a MWE model:
ode = @ODEmodel( x1'(t) = a*x1(t)+sin(u(t)), y(t) = x1(t) )
This leads to the following error:
[ Info: Summary of the model: [ Info: State variables: x1 [ Info: Parameters: a, sin [ Info: Inputs: u [ Info: Outputs: y ERROR: Number of variables does not match number o f values
where sin is recognized as parameter. Am I doing something wrong or is this a limitation of your toolbox?
Thx.

Hi,

I am new to both structural identifiability analysis and Julia. So I hope the answer to my question is not completely obvious. I am able to obtain the local SI results but when commuting the global ones, the output in the terminal always stalls at "Computing IO-equations". I have left it running on an HPC over night.

Here's the code:

#import Pkg; 
#using Pkg
#Pkg.add("StructuralIdentifiability")
#Pkg.add("OrdinaryDiffEq")

using StructuralIdentifiability, OrdinaryDiffEq

odeMyModel = @ODEmodel(
	Complex'(t) = 1/C*(((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) - (ke_Complex*Complex(t)*C) - ((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C)),
	Complex_2'(t) = 1/C*(((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C) - (ke_Complex_2*Complex_2(t)*C)),
        Drug'(t) = 1/C*(-(ke_Drug*Drug(t)*C) - ((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) + (45/100*ka*Drug_SC(t))),
        free_receptor'(t) = 1/C*(-((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) + (66/2500*C) - (kdeg_free_receptor*free_receptor(t)*C) - ((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C)),
	Drug_SC'(t) = -(45/100*ka*Drug_SC(t)) - ((1-45/100)*ka*Drug_SC(t)) + u_SC(t),
        y1(t) = Drug(t),                                   # output free Drug
        y2(t) = free_receptor(t) + Complex(t) + 2*Complex_2(t)      # output total receptor
)

local_id_all = assess_local_identifiability(odeMyModel, 0.99)

assess_identifiability(odeMyModel,0.99) 

The parameters, state variables, inputs and outputs seem to be correctly recognized:

[ Info: Summary of the model:
[ Info: State variables: Drug, free_receptor, Complex, Complex_2, Drug_SC
[ Info: Parameters: ke_Complex_2, kdeg_free_receptor, ka, kon, ke_Drug, ke_Complex, kon_2, koff, C
[ Info: Inputs: u_SC
[ Info: Outputs: y1, y2
[ Info: Assessing local identifiability
[ Info: Local identifiability assessed in 0.091845932 seconds
[ Info: Assessing global identifiability
[ Info: Computing IO-equations

But it doesn't move past there.

Also, I was looking to include initial conditions but have failed in the attempt in multiple ways. I think it might not be possible yet. If that is the case, can you recommend an alternative for SI - maybe even additionally allowing bounds for parameters?

Thanks a lot in advance for your help!

Some biological models are represented via difference equations of the form: $x(t + 1) = f\big( x(t), \theta \big)$. According to Chris Rackauckas on Twitter, an extension of the identifiability functions could be made for those systems using the DiscreteSystem from ModelingToolkit.

Some references on the subject are:

Nõmm, Sven, and C. H. Moog. "Identifiability of discrete-time nonlinear systems." IFAC Proceedings Volumes 37.13 (2004): 333-338.

Nõmm, Sven, and Claude H. Moog. "Further results on identifiability of discrete-time nonlinear systems." Automatica 68 (2016): 69-74.

Hello,
When I run this system of ODEs on StructuralIdentifiability vs SIAN, I get different results,

solver = :SIAN

if solver == :StructuralIdentifiability
    using StructuralIdentifiability
elseif solver == :SIAN
    using SIAN
end

ode = @ODEmodel(
    n'(t) = (co(t)/(co(t)+Ko)) * (cg(t)/(cg(t)+Kg)) * (Kl/(cl(t)+Kl)) * beta * n(t) * (1-n(t)/nmax) - delta * n(t),
    co'(t) = 0,
    cg'(t) = -Vg * n(t) * (cg(t) / (cg(t) + cbarg)),
    cl'(t) = - 2 * (-Vo * n(t) * (co(t) / (co(t) + cbaro))) + (-Vg * n(t) * (cg(t) / (cg(t) + cbaro)))/3,
    y1(t) = n(t),
    y2(t) = cg(t),
    y3(t) = cl(t),
    y4(t) = co(t)
)

if solver == :StructuralIdentifiability
    output = assess_identifiability(ode, 0.99)
elseif solver == :SIAN
    output_single = identifiability_ode(ode, get_parameters(ode; initial_conditions=false); p=0.99)
end

I assumed that the reason might be because StructuralIdentifiability assumes the initial conditions to be known while SIAN assumes them to be unknown. Is there a way for me to set the initial conditions to know/unknown on either package?

I just installed & tried this out.

julia> using StructuralIdentifiability
julia> ode = @ODEmodel(
           x1'(t) = -(a01 + a21) * x1(t) + a12 * x2(t) + u(t),
           x2'(t) = a21 * x1(t) - a12 * x2(t) - x3(t) / b,
           x3'(t) = x3(t),
           y(t) = x2(t)
       )
julia> assess_identifiability(ode)
[ Info: Assessing local identifiability
[ Info: Local identifiability assessed in 3.320626708 seconds
[ Info: Assessing global identifiability
[ Info: Computing IO-equations
┌ Info: Computed in 4.313144799 seconds
│   :ioeq_time = :ioeq_time
└   ioeq_time = 4.313144799
[ Info: Computing Wronskians
ERROR: MethodError: Cannot `convert` an object of type 
  Nemo.gfp_abs_series to an object of type 
  AbstractAlgebra.Generic.AbsSeries
Closest candidates are:
  convert(::Type{T1}, ::CxxWrap.CxxWrapCore.SmartPointer{DerivedT}) where {BaseT, T1<:BaseT, DerivedT<:BaseT} at /Users/lmorton/.julia/packages/CxxWrap/OcN1Z/src/CxxWrap.jl:276
  convert(::Type{T}, ::T) where T at essentials.jl:205
Stacktrace:
 [1] setindex!(h::Dict{AbstractAlgebra.MPolyElem, AbstractAlgebra.Generic.AbsSeries}, v0::Nemo.gfp_abs_series, key::Nemo.gfp_mpoly)
   @ Base ./dict.jl:382
 [2] wronskian(io_equations::Dict{Nemo.fmpq_mpoly, Nemo.fmpq_mpoly}, ode::ODE{Nemo.fmpq_mpoly})
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/wronskian.jl:223
 [3] macro expansion
   @ ./timing.jl:287 [inlined]
 [4] assess_global_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Any}, p::Float64; var_change::Symbol, gb_method::Symbol)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:207
 [5] assess_global_identifiability
   @ ~/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:200 [inlined]
 [6] macro expansion
   @ ./timing.jl:287 [inlined]
 [7] assess_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Nemo.fmpq_mpoly}, p::Float64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:110
 [8] assess_identifiability (repeats 2 times)
   @ ~/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:68 [inlined]
 [9] top-level scope
   @ REPL[7]:1

My environment:

julia> versioninfo()
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i9-9980HK CPU @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 
  JULIA_PKG_DEVDIR = /Users/lmorton/Code/Jdev/

Packages:

(Julia) pkg> status
      Status `~/Code/Julia/Project.toml`
  [39de3d68] AxisArrays v0.4.4
  [a134a8b2] BlackBoxOptim v0.6.0
  [82cc6244] DataInterpolations v3.6.0
  [aae7a2af] DiffEqFlux v1.43.0
  [0c46a032] DifferentialEquations v6.18.0
  [587475ba] Flux v0.12.6
  [f6369f11] ForwardDiff v0.10.19
  [a75be94c] GalacticOptim v2.0.3
  [bb4c363b] GridInterpolations v1.1.2
  [f67ccb44] HDF5 v0.15.6
  [615f187c] IfElse v0.1.0
  [b964fa9f] LaTeXStrings v1.2.1
  [2ee39098] LabelledArrays v1.6.4
  [23fbe1c1] Latexify v0.15.6
  [961ee093] ModelingToolkit v6.3.0 `~/Code/Jdev/ModelingToolkit`
  [8913a72c] NonlinearSolve v0.3.9
  [429524aa] Optim v1.4.1
  [1dea7af3] OrdinaryDiffEq v5.61.1
  [a03496cd] PlotlyBase v0.8.7
  [91a5bcdd] Plots v1.20.1
  [438e738f] PyCall v1.92.3
  [d330b81b] PyPlot v2.9.0
  [0bca4576] SciMLBase v1.18.4
  [220ca800] StructuralIdentifiability v0.1.0
  [0c5d862f] Symbolics v3.0.0
  [1986cc42] Unitful v1.9.0
  [e88e6eb3] Zygote v0.6.17
  [700de1a5] ZygoteRules v0.2.1

Unable to load StructuralIdentifiability.jl after installation without errors. Tried multiple clean installations, but the error persists. Following is the StackTrace error and version numbers:

[220ca800] StructuralIdentifiability v0.3.6
[961ee093] ModelingToolkit v8.4.1

[ Info: Precompiling StructuralIdentifiability [220ca800-aa68-49bb-acd8-6037fa93a544]
ERROR: LoadError: UndefVarError: libsingular not defined
Stacktrace:
[1] getproperty(x::Module, f::Symbol)
@ Base .\Base.jl:26
[2] top-level scope
@ C:\Users\Bharadwaj.julia\packages\Singular\nGiup\src\Singular.jl:66
[3] include
@ .\Base.jl:386 [inlined]
[4] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base .\loading.jl:1235
[5] top-level scope
@ none:1
[6] eval
@ .\boot.jl:360 [inlined]
[7] eval(x::Expr)
@ Base.MainInclude .\client.jl:446
[8] top-level scope
@ none:1
in expression starting at C:\Users\Bharadwaj.julia\packages\Singular\nGiup\src\Singular.jl:1
ERROR: LoadError: Failed to precompile Singular [bcd08a7b-43d2-5ff7-b6d4-c458787f915c] to C:\Users\Bharadwaj.julia\compiled\v1.6\Singular\jl_9131.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::Base.TTY, internal_stdout::Base.TTY, ignore_loaded_modules::Bool)
@ Base .\loading.jl:1385
[3] compilecache(pkg::Base.PkgId, path::String)
@ Base .\loading.jl:1329
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1043
[5] require(uuidkey::Base.PkgId)
@ Base .\loading.jl:936
[6] require(into::Module, mod::Symbol)
@ Base .\loading.jl:923
[7] include
@ .\Base.jl:386 [inlined]
[8] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base .\loading.jl:1235
[9] top-level scope
@ none:1
[10] eval
@ .\boot.jl:360 [inlined]
[11] eval(x::Expr)
@ Base.MainInclude .\client.jl:446
[12] top-level scope
@ none:1
in expression starting at C:\Users\Bharadwaj.julia\packages\GroebnerBasis\89ud7\src\GroebnerBasis.jl:1
ERROR: LoadError: Failed to precompile GroebnerBasis [e39c9192-ea4d-5e15-9584-a488c6d614bd] to C:\Users\Bharadwaj.julia\compiled\v1.6\GroebnerBasis\jl_8BC2.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::Base.TTY, internal_stdout::Base.TTY, ignore_loaded_modules::Bool)
@ Base .\loading.jl:1385
[3] compilecache(pkg::Base.PkgId, path::String)
@ Base .\loading.jl:1329
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1043
[5] require(uuidkey::Base.PkgId)
@ Base .\loading.jl:936
[6] require(into::Module, mod::Symbol)
@ Base .\loading.jl:923
[7] include
@ .\Base.jl:386 [inlined]
[8] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::Nothing)
@ Base .\loading.jl:1235
[9] top-level scope
@ none:1
[10] eval
@ .\boot.jl:360 [inlined]
[11] eval(x::Expr)
@ Base.MainInclude .\client.jl:446
[12] top-level scope
@ none:1
in expression starting at C:\Users\Bharadwaj.julia\packages\StructuralIdentifiability\0qdZi\src\StructuralIdentifiability.jl:1
ERROR: Failed to precompile StructuralIdentifiability [220ca800-aa68-49bb-acd8-6037fa93a544] to C:\Users\Bharadwaj.julia\compiled\v1.6\StructuralIdentifiability\jl_84FE.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::Base.TTY, internal_stdout::Base.TTY, ignore_loaded_modules::Bool)
@ Base .\loading.jl:1385
[3] compilecache(pkg::Base.PkgId, path::String)
@ Base .\loading.jl:1329
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1043
[5] require(uuidkey::Base.PkgId)
@ Base .\loading.jl:936
[6] require(into::Module, mod::Symbol)
@ Base .\loading.jl:923

I was just trying to run the example code copied directly from here: https://github.com/SciML/StructuralIdentifiability.jl/blob/master/examples/CRN.jl

and am getting the output pasted below. I also tried running the library on an ODE model of my own and other examples, always getting the same output. Is this a known bug? I tried running `Pkg.update("StructuralIndentifiability") to see if that helped, but no luck.

`Singular.jl, based on
SINGULAR /
A Computer Algebra System for Polynomial Computations / Singular.jl: 0.5.7
0< Singular : 4.2.0p1
by: W. Decker, G.-M. Greuel, G. Pfister, H. Schoenemann
FB Mathematik der Universitaet, D-67653 Kaiserslautern \

┌ Info: Summary of the model:
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/ODE.jl:299
┌ Info: State variables: x2, x5, x6, x3, x4, x1
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/ODE.jl:300
┌ Info: Parameters: k5, k3, k4, k2, k6, k1
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/ODE.jl:301
┌ Info: Inputs:
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/ODE.jl:302
┌ Info: Outputs: y1, y2
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/ODE.jl:303
┌ Info: Assessing local identifiability
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:92
┌ Info: Local identifiability assessed in 4.943092853 seconds
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:94
┌ Info: Assessing global identifiability
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:109
┌ Info: Computing IO-equations
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:200
┌ Info: Computed in 8.872011279 seconds
│ :ioeq_time = ioeq_time
│ ioeq_time = 8.872011279
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:203
┌ Info: Computing Wronskians
└ @ StructuralIdentifiability /home/nao5/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:206
MethodError: Cannot convert an object of type
Nemo.gfp_abs_series to an object of type
AbstractAlgebra.Generic.AbsSeries
Closest candidates are:
convert(::Type{T1}, ::CxxWrap.CxxWrapCore.SmartPointer{DerivedT}) where {BaseT, T1<:BaseT, DerivedT<:BaseT} at /home/nao5/.julia/packages/CxxWrap/OcN1Z/src/CxxWrap.jl:276
convert(::Type{T}, ::T) where T at essentials.jl:205

Stacktrace:
[1] setindex!(h::Dict{AbstractAlgebra.MPolyElem, AbstractAlgebra.Generic.AbsSeries}, v0::Nemo.gfp_abs_series, key::Nemo.gfp_mpoly)
@ Base ./dict.jl:382
[2] wronskian(io_equations::Dict{Nemo.fmpq_mpoly, Nemo.fmpq_mpoly}, ode::ODE{Nemo.fmpq_mpoly})
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/wronskian.jl:223
[3] macro expansion
@ ./timing.jl:287 [inlined]
[4] assess_global_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Any}, p::Float64; var_change::Symbol, gb_method::Symbol)
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:207
[5] assess_global_identifiability
@ ~/.julia/packages/StructuralIdentifiability/wOUIO/src/global_identifiability.jl:200 [inlined]
[6] macro expansion
@ ./timing.jl:287 [inlined]
[7] assess_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Nemo.fmpq_mpoly}, p::Float64)
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:110
[8] assess_identifiability (repeats 2 times)
@ ~/.julia/packages/StructuralIdentifiability/wOUIO/src/StructuralIdentifiability.jl:68 [inlined]
[9] top-level scope
@ ./timing.jl:210 [inlined]
[10] top-level scope
@ ./In[1]:0
[11] eval
@ ./boot.jl:360 [inlined]
[12] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1116`

Hello,

I have implemented the minimal glucose model (Bergman et al 1979) with Catalyst.jl, from which I get an ODESystem object. As described in my paper of reference, the small model includes interpolated insulin data (to avoid regulation modeling):

using Catalyst, DifferentialEquations, StructuralIdentifiability, Interpolations

tticks = [0.0, 15.0, 30.0, 60.0, 90.0, 120.0, 180.0, 240.0, 300.0]
datamean = [6.12,18.10, 32.57, 25.77, 27.94, 33.47, 35.18, 39.88, 46.20] 

# minimal model
Iintpol(t) = linear_interpolation(tticks,datamean)(t)
@register_symbolic Iintpol(t)
minimalModel_rn = @reaction_network minimalmodel begin
    B_0,     ∅   ⇒   G_p              # hepatic glucose production
    X,       G_p --> ∅                   # glucose insulin-dependent regulation
    k_greg,  G_p --> 2*G_p      # glucose "self-regulation"
    k_ipan*Iintpol(t), ∅ ⇒ X    # blood insulin --> tissue insulin            
    k_ireg, X --> 2*X                  # insulin "self-regulation"
end B_0 k_ipan k_ireg k_greg

minimalModel_ode = convert(ODESytem, minimalModel_rn)

When I try to perform a full structural analysis, I obtain an issue because of Iintpol(t):

@variables t y(t)
assess_identifiability(minimalModel_ode,measured_quantities = [y ~ states(minimalModel_ode)[1]]) # G_p is the observable

output:

[ Info: Preproccessing `ModelingToolkit.ODESystem` object
ArgumentError: Function Iintpol is not supported

Is there any workaround for this case?

In addition to knowing the ODE system and time series data for some outputs, one frequently also knows the initial conditions for some of the states. It would be great to use this information to clarify the identifiability analysis.
In the case of local identifiability, the algorithm computes the observability matrix at a random point using the algorithm by Sedoglavic.

Instead of using a completely random point one could use the user-specified initial conditions for some of the states. The issue with such an algorithm is that it may produce wrong result as shown here. I am not aware of any algorithm resolving the issue with reasonable complexity. What we could still do is:

• One can always get a correct result from an observability matrix if it will be computed up to high enough order. The problem is that we do not know how far to go but we could use some heuristics, say, double the order until the rank stabilises.
• In such an algorithm the only way things can go wrong is when a locally identifiable parameter is classified as non identifiable. Therefore, we can show a warning but only in the case there are non identifiable parameters in the result.

With Catalyst (and other packages) we generate models ranging from small to large that we can encode as ModeligToolkit ODESystems (or other stuff). I want to assess the identifiability of these systems. However, I am unsure how to apply certain functions in the system to these (which seems primarily designed for StructuralIdentifiabiltiy ODEs).

Say that I have an identifiability problem, but know the value of certain initial conditions and parameters, how do I give StructuralIdentifiabiltiy this information? My only lead, the set_parameter_values function, seems to only work on ODEs defined by StructutralIdentifiability.

Hi,
First, thank for this cool package.

I am trying to compose model based on several submodels written in MTK. To do I have input variable in each submodel which allows me to define input and output from other states variables in the composed model.

However, this seems to interfere with the identifiability analysis. I have written a MWE with a Lotka-Volterra model exemplifying the issue. The non-composed models works fine but the composed version fails (see error below). Is this expected? There seem to similar issues posted here.

using DifferentialEquations
using CairoMakie
using ModelingToolkit
using StructuralIdentifiability


function lotka_volterra!(du, u, p, t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α * x - β * x * y ## rabbits
    du[2] = dy = -δ * y + γ * x * y ## wolves
end

# Initial condition
u0 = [1.0, 2.0]

# Simulation interval
tspan = (0.0, 10.0)

# LV equation parameter. p = [α, β, δ, γ]
p = [2.5, 3.0, 7.0, 2.0]

# Setup the ODE problem, then solve
prob = ODEProblem(lotka_volterra!, u0, tspan, p)
sol = solve(prob)
datasol = solve(prob, saveat = 1)
data = Array(datasol)

f = Figure()
ax = Axis(f[1,1])
tplot = collect(tspan[1]:0.01:tspan[2])
lines!(ax, tplot, [sol(j)[1] for j in tplot])
lines!(ax, tplot, [sol(j)[2] for j in tplot])
scatter!(ax, collect(0:1:10), data[1,:])
scatter!(ax, collect(0:1:10), data[2,:])
f

@variables t
function rabbits_creator(;name)
    ps = @parameters α=1.5
    vars = @variables x(t)=1.0 z(t)=0.0 [input = true]
    D = Differential(t)
    equs = [
        D(x) ~ α *x + z
    ]

    ODESystem(equs, t, vars, ps; name = name)
end

function wolves_creator(;name)
    ps = @parameters δ = 3.0
    vars = @variables y(t)=1.0 q(t)=0.0 [input = true]
    D = Differential(t)
    equs = [
        D(y) ~ -δ * y + q
    ]

    ODESystem(equs, t, vars, ps; name = name)
end

function lotka_volterra_creator(;name)
    @named wolves = wolves_creator()
    @named rabbits = rabbits_creator()

    ps = @parameters β=1.0 γ=1.0
    D = Differential(t)
    
    eqs = [ 
        rabbits.z ~  - β *  wolves.y * rabbits.x, 
        wolves.q ~ γ *  wolves.y *rabbits.x]
           
    compose(ODESystem(eqs, t, [],ps; name= name), wolves, rabbits)
end

@named ltk_mtk = lotka_volterra_creator()
simp_ltk_mtk = structural_simplify(ltk_mtk)
equations(simp_ltk_mtk)
parameters(simp_ltk_mtk)
@unpack wolves₊δ, rabbits₊α , β, γ,  wolves₊y, rabbits₊x = simp_ltk_mtk

u0_sim = [
    rabbits₊x => u0[1],
    wolves₊y => u0[2]
]

pars_sim = [
    rabbits₊α => p[1],
    β =>  p[2],
    wolves₊δ => p[3],
    γ => p[4]
]

sol_mtk = solve(ODEProblem(simp_ltk_mtk, u0_sim, tspan, pars_sim))

prob = ODEProblem(lotka_volterra!, u0, tspan, p)
sol = solve(prob)

f = Figure()
ax = Axis(f[1,1])
tplot = collect(tspan[1]:0.01:tspan[2])
lines!(ax, tplot, [sol_mtk(j)[1] for j in tplot])
lines!(ax, tplot, [sol_mtk(j)[2] for j in tplot])
lines!(ax, tplot, [sol(j)[1] for j in tplot], linestyle=:dash)
lines!(ax, tplot, [sol(j)[2] for j in tplot], linestyle=:dash)
f

ModelingToolkit.observed(simp_ltk_mtk)
original_sys = modelingtoolkitize(prob)

@unpack wolves₊y, rabbits₊x = simp_ltk_mtk
@unpack x₁, x₂ = original_sys
@variables y1(t) y2(t)

This works

measured_quantities = [y1 ~ x₁, y2 ~ x₂]
assess_identifiability(original_sys, measured_quantities = measured_quantities)

but this fails with ERROR: MethodError: no method matching +(::Nemo.fmpq_mpoly, ::Term{Real, Base.ImmutableDict{DataType, Any}})

measured_quantities = [y1 ~ rabbits₊x, y2 ~ wolves₊y]
assess_identifiability(simp_ltk_mtk, measured_quantities = measured_quantities)

Hello,
I have a model for which I want to perform identifiability analysis. Unfortunately, running it simply crashes Julia, without any errors or other useful information to help figure out what is going on. I am not really sure how to deal with this.

Here's the code:

using Catalyst, DifferentialEquations, DelimitedFiles, StructuralIdentifiability

jakobs_network = @reaction_network begin
    k0, Rf3PPd + C6H5Br -->  Rf3PPdC6H5Br
    k1, Rf3PPdC6H5Br + Re3SiONa --> NaBr + Rf3PPdC6H5Re3SiO
    k2, Rf3PPdC6H5Re3SiO + Rg3SiH --> Rf3PPd + C6H6 + Re3Si2ORg3
    k3, Rf3PPdC6H5Br + RaP2CuC2H4Rh -->  Rf3PPdC8H9Rh + RaP2CuBr
    k4, RaP2CuBr + Re3SiONa --> NaBr + RaP2CuRe3SiO
    k5, RaP2CuRe3SiO + Rg3SiH --> RaP2CuH + Re3Si2ORg3
    k6, RaP2CuH + C2H3Rh --> RaP2CuC2H4Rh
    k7, Rf3PPdC8H9Rh --> Rf3PPd + C8H9Rh
    k8, Rf3PPd + Rg3SiH --> PdH + Rg3SiRf3P
    k9, PdH + PdH --> Pd2 + H2
    k10, RaP2CuC2H4Rh + H2O --> C2H5Rh + RaP2CuOH
    k11, RaP2CuOH + Rg3SiH --> RaP2CuRg3SiO + H2
    k12, RaP2CuRg3SiO + Rg3SiH --> RaP2CuH + Rg6Si2O
    k13, C2H3Rh --> C2H3 + Rh
end k0 k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13;

osys = convert(ODESystem, jakobs_network)
meassured_species_idxs = [8,15,16,22]
meassured_species = map(s -> states(osys)[s], meassured_species_idxs)

@variables t y1(t) y2(t) y3(t) y4(t)
eq1 = y1 ~ meassured_species[1]
eq2 = y2 ~ meassured_species[2]
eq3 = y3 ~ meassured_species[3]
eq4 = y4 ~ meassured_species[4]

assess_identifiability(osys; measured_quantities=[eq1, eq2, eq3, eq4])

and here's the environment:

(Scripts) pkg> status
Status `~/Desktop/Research Projects/Jakobs-Problem/Scripts/Project.toml`
  [479239e8] Catalyst v12.3.1
  [1130ab10] DiffEqParamEstim v1.28.0
  [0c46a032] DifferentialEquations v7.6.0
  [86b6b26d] Evolutionary v0.11.1
  [7073ff75] IJulia v1.23.3
  [429524aa] Optim v1.7.4
⌃ [1dea7af3] OrdinaryDiffEq v6.35.0
⌃ [91a5bcdd] Plots v1.37.0
  [220ca800] StructuralIdentifiability v0.4.2
  [8bb1440f] DelimitedFiles
Info Packages marked with ⌃ have new versions available and may be upgradable.

I had a go at defining a system with state h, inputs u and outputs y,

using StructuralIdentifiability

ode = @ODEmodel(
    h1'(t) = -a/A * sqrt(18*h1(t)) + a/A*sqrt(18*h3(t)) +     γ*k/A * u1(t),
    h2'(t) = -a/A * sqrt(18*h2(t)) + a/A*sqrt(18*h4(t)) +     γ*k/A * u2(t),
    h3'(t) = -a/A*sqrt(18*h3(t))                          + (1-γ)*k/A * u2(t),
    h4'(t) = -a/A*sqrt(18*h4(t))                          + (1-γ)*k/A * u1(t),
	y1(t) = h1(t),
    y2(t) = h2(t),
)

local_id = assess_local_identifiability(ode, 0.99)

it looks like parameters, states, inputs and outputs were correctly identified, except it thinks that sqrt is a parameter, and as a consequence it fails:

[ Info: Summary of the model:
[ Info: State variables: h1, h2, h3, h4
[ Info: Parameters: A, k, a, γ, sqrt
[ Info: Inputs: u1, u2
[ Info: Outputs: y1, y2
ERROR: Number of variables does not match number of values

If I remove sqrt, it accepts the definition. I couldn't figure out from the docs whether or not this is a known limitation or a bug?

Running the MTK tutorial and testing this package fail because of this libsingular not defined error. It has been there since the beginning. Would be good if this issue could finally be resolved so that one can give a real test of this excellent package.

           _

_ _ ()_ | Documentation: https://docs.julialang.org
() | () () |
_ _ | | __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ ` | |
| | || | | | (| | | Version 1.7.1 (2021-12-22)
/ |_'|||_'_| | Official https://julialang.org/ release
|__/ |

(@v1.7) pkg> test StructuralIdentifiability
Testing StructuralIdentifiability
Status C:\Users\Denis\AppData\Local\Temp\jl_IbOUUA\Project.toml
[c3fe647b] AbstractAlgebra v0.20.1
[e39c9192] GroebnerBasis v0.3.3
[c8e1da08] IterTools v1.4.0
[1914dd2f] MacroTools v0.5.9
[2edaba10] Nemo v0.25.1
[27ebfcd6] Primes v0.5.1
[bcd08a7b] Singular v0.5.8
[220ca800] StructuralIdentifiability v0.1.0
[98d24dd4] TestSetExtensions v2.0.0
[ade2ca70] Dates @stdlib/Dates
[37e2e46d] LinearAlgebra @stdlib/LinearAlgebra
[56ddb016] Logging @stdlib/Logging
[8dfed614] Test @stdlib/Test
Status C:\Users\Denis\AppData\Local\Temp\jl_IbOUUA\Manifest.toml
[c3fe647b] AbstractAlgebra v0.20.1
[b99e7846] BinaryProvider v0.5.10
[1f15a43c] CxxWrap v0.11.2
[ab62b9b5] DeepDiffs v1.2.0
[e39c9192] GroebnerBasis v0.3.3
[3e1990a7] Hecke v0.10.15
[c8e1da08] IterTools v1.4.0
[692b3bcd] JLLWrappers v1.3.0
[472f376f] LoadFlint v0.5.1
[1914dd2f] MacroTools v0.5.9
[2edaba10] Nemo v0.25.1
[21216c6a] Preferences v1.2.3
[27ebfcd6] Primes v0.5.1
[fb686558] RandomExtensions v0.4.3
[ae029012] Requires v1.2.0
[bcd08a7b] Singular v0.5.8
[90137ffa] StaticArrays v1.3.1
[220ca800] StructuralIdentifiability v0.1.0
[98d24dd4] TestSetExtensions v2.0.0
[e21ec000] Antic_jll v0.200.400+0
[d9960996] Arb_jll v200.1900.1+0
[fcfa6d1b] Calcium_jll v0.400.0+0
[e134572f] FLINT_jll v200.700.100+0
[e8aa6df9] GLPK_jll v5.0.1+0
[43d676ae] Singular_jll v402.0.104+0
[f07e07eb] cddlib_jll v0.94.13+0
[006bdf2e] gb_jll v0.17.0+0
[1493ae25] lib4ti2_jll v1.6.10+0
[3eaa8342] libcxxwrap_julia_jll v0.8.8+1
[ae4fbd8f] libsingular_julia_jll v0.15.7+0
[0dad84c5] ArgTools @stdlib/ArgTools
[56f22d72] Artifacts @stdlib/Artifacts
[2a0f44e3] Base64 @stdlib/Base64
[ade2ca70] Dates @stdlib/Dates
[8ba89e20] Distributed @stdlib/Distributed
[f43a241f] Downloads @stdlib/Downloads
[b77e0a4c] InteractiveUtils @stdlib/InteractiveUtils
[b27032c2] LibCURL @stdlib/LibCURL
[76f85450] LibGit2 @stdlib/LibGit2
[8f399da3] Libdl @stdlib/Libdl
[37e2e46d] LinearAlgebra @stdlib/LinearAlgebra
[56ddb016] Logging @stdlib/Logging
[d6f4376e] Markdown @stdlib/Markdown
[ca575930] NetworkOptions @stdlib/NetworkOptions
[44cfe95a] Pkg @stdlib/Pkg
[de0858da] Printf @stdlib/Printf
[9abbd945] Profile @stdlib/Profile
[3fa0cd96] REPL @stdlib/REPL
[9a3f8284] Random @stdlib/Random
[ea8e919c] SHA @stdlib/SHA
[9e88b42a] Serialization @stdlib/Serialization
[6462fe0b] Sockets @stdlib/Sockets
[2f01184e] SparseArrays @stdlib/SparseArrays
[10745b16] Statistics @stdlib/Statistics
[fa267f1f] TOML @stdlib/TOML
[a4e569a6] Tar @stdlib/Tar
[8dfed614] Test @stdlib/Test
[cf7118a7] UUIDs @stdlib/UUIDs
[4ec0a83e] Unicode @stdlib/Unicode
[e66e0078] CompilerSupportLibraries_jll @stdlib/CompilerSupportLibraries_jll
[781609d7] GMP_jll @stdlib/GMP_jll
[deac9b47] LibCURL_jll @stdlib/LibCURL_jll
[29816b5a] LibSSH2_jll @stdlib/LibSSH2_jll
[3a97d323] MPFR_jll @stdlib/MPFR_jll
[c8ffd9c3] MbedTLS_jll @stdlib/MbedTLS_jll
[14a3606d] MozillaCACerts_jll @stdlib/MozillaCACerts_jll
[4536629a] OpenBLAS_jll @stdlib/OpenBLAS_jll
[83775a58] Zlib_jll @stdlib/Zlib_jll
[8e850b90] libblastrampoline_jll @stdlib/libblastrampoline_jll
[8e850ede] nghttp2_jll @stdlib/nghttp2_jll
[3f19e933] p7zip_jll @stdlib/p7zip_jll
Precompiling project...
✗ Singular
✗ GroebnerBasis
✗ StructuralIdentifiability
0 dependencies successfully precompiled in 15 seconds (35 already precompiled)
3 dependencies errored. To see a full report either run import Pkg; Pkg.precompile() or load the packages
Testing Running tests...
ERROR: LoadError: UndefVarError: libsingular not defined
Stacktrace:
[1] getproperty(x::Module, f::Symbol)
@ Base .\Base.jl:35
[2] top-level scope
@ C:\Users\Denis.julia\packages\Singular\nGiup\src\Singular.jl:66
[3] include
@ .\Base.jl:418 [inlined]
[4] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base .\loading.jl:1318
[5] top-level scope
@ none:1
[6] eval
@ .\boot.jl:373 [inlined]
[7] eval(x::Expr)
@ Base.MainInclude .\client.jl:453
[8] top-level scope
@ none:1
in expression starting at C:\Users\Denis.julia\packages\Singular\nGiup\src\Singular.jl:1
ERROR: LoadError: Failed to precompile Singular [bcd08a7b-43d2-5ff7-b6d4-c458787f915c] to C:\Users\Denis.julia\compiled\v1.7\Singular\jl_B9D.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, ignore_loaded_modules::Bool)
@ Base .\loading.jl:1466
[3] compilecache(pkg::Base.PkgId, path::String)
@ Base .\loading.jl:1410
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1120
[5] require(uuidkey::Base.PkgId)
@ Base .\loading.jl:1013
[6] require(into::Module, mod::Symbol)
@ Base .\loading.jl:997
[7] include
@ .\Base.jl:418 [inlined]
[8] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base .\loading.jl:1318
[9] top-level scope
@ none:1
[10] eval
@ .\boot.jl:373 [inlined]
[11] eval(x::Expr)
@ Base.MainInclude .\client.jl:453
[12] top-level scope
@ none:1
in expression starting at C:\Users\Denis.julia\packages\GroebnerBasis\89ud7\src\GroebnerBasis.jl:1
ERROR: LoadError: Failed to precompile GroebnerBasis [e39c9192-ea4d-5e15-9584-a488c6d614bd] to C:\Users\Denis.julia\compiled\v1.7\GroebnerBasis\jl_6AB.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, ignore_loaded_modules::Bool)
@ Base .\loading.jl:1466
[3] compilecache(pkg::Base.PkgId, path::String)
@ Base .\loading.jl:1410
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1120
[5] require(uuidkey::Base.PkgId)
@ Base .\loading.jl:1013
[6] require(into::Module, mod::Symbol)
@ Base .\loading.jl:997
[7] include
@ .\Base.jl:418 [inlined]
[8] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base .\loading.jl:1318
[9] top-level scope
@ none:1
[10] eval
@ .\boot.jl:373 [inlined]
[11] eval(x::Expr)
@ Base.MainInclude .\client.jl:453
[12] top-level scope
@ none:1
in expression starting at C:\Users\Denis.julia\packages\StructuralIdentifiability\wOUIO\src\StructuralIdentifiability.jl:1
ERROR: LoadError: Failed to precompile StructuralIdentifiability [220ca800-aa68-49bb-acd8-6037fa93a544] to C:\Users\Denis.julia\compiled\v1.7\StructuralIdentifiability\jl_81.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, ignore_loaded_modules::Bool)
@ Base .\loading.jl:1466
[3] compilecache(pkg::Base.PkgId, path::String)
@ Base .\loading.jl:1410
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1120
[5] require(uuidkey::Base.PkgId)
@ Base .\loading.jl:1013
[6] require(into::Module, mod::Symbol)
@ Base .\loading.jl:997
[7] include(fname::String)
@ Base.MainInclude .\client.jl:451
[8] top-level scope
@ none:6
in expression starting at C:\Users\Denis.julia\packages\StructuralIdentifiability\wOUIO\test\runtests.jl:1
ERROR: Package StructuralIdentifiability errored during testing

(@v1.7) pkg>

Hi, thank you for the great package! I was trying to do the structural global identifiability on a simple SIR model with a nonlinear term in the transmission rate, it works perfectly

sir_ode = @ODEmodel(
    S'(t) = -τ*S(t)*I(t)/(1+ξ*I(t)),
    I'(t) = τ*S(t)*I(t)/(1+ξ*I(t)) - γi*I(t),
    R'(t) = γi*I(t),
    y1(t) = I(t),
)

However, when I tried to extend it to an SEIR model it doesn't work, here is a MWE

seir_ode = @ODEmodel(
    S'(t) = -τ*S(t)*I(t)/(1+ξ*I(t)),
    E'(t) = τ*S(t)*I(t)/(1+ξ*I(t)) - σ*E(t),
    I'(t) = σ*E(t) - γi*I(t),
    R'(t) = γi*I(t),
    y1(t) = σ*E(t),
)
res_seir = assess_local_identifiability(seir_ode, seir_ode.parameters) # works
res_seir = assess_identifiability(seir_ode) # error

ERROR: MethodError: no method matching ODE{Nemo.fmpq_mpoly}(::Dict{Nemo.fmpq_mpoly, RingElem}, ::Dict{Nemo.fmpq_mpoly, Nemo.fmpq_mpoly}, ::Vector{Nemo.fmpq_mpoly})
Closest candidates are:
  ODE{P}(::Dict{P, <:Union{AbstractAlgebra.Generic.Frac{P}, P}}, ::Dict{P, <:Union{AbstractAlgebra.Generic.Frac{P}, P}}, ::Vector{P}) where P<:(MPolyElem{<:FieldElem}) at ~/.julia/packages/StructuralIdentifiability/gEpsR/src/ODE.jl:17
Stacktrace:
  [1] ode_aux(ode::ODE{Nemo.fmpq_mpoly}, submodel::Set{Nemo.fmpq_mpoly})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/gEpsR/src/submodels.jl:120
  [2] #431
    @ ./array.jl:0 [inlined]
  [3] iterate
    @ ./generator.jl:47 [inlined]
  [4] collect(itr::Base.Generator{Vector{Set{Nemo.fmpq_mpoly}}, StructuralIdentifiability.var"#431#432"{ODE{Nemo.fmpq_mpoly}}})
    @ Base ./array.jl:787
  [5] submodel2ode
    @ ~/.julia/packages/StructuralIdentifiability/gEpsR/src/submodels.jl:140 [inlined]
  [6] find_submodels(ode::ODE{Nemo.fmpq_mpoly})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/gEpsR/src/submodels.jl:179
  [7] assess_global_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Any}, known::Vector{Nemo.fmpq_mpoly}, p::Float64; var_change::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/gEpsR/src/global_identifiability.jl:215
  [8] assess_global_identifiability
    @ ~/.julia/packages/StructuralIdentifiability/gEpsR/src/global_identifiability.jl:207 [inlined]
  [9] macro expansion
    @ ./timing.jl:382 [inlined]
 [10] assess_identifiability(ode::ODE{Nemo.fmpq_mpoly}, funcs_to_check::Vector{Nemo.fmpq_mpoly}, p::Float64)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/gEpsR/src/StructuralIdentifiability.jl:123
 [11] assess_identifiability (repeats 2 times)
    @ ~/.julia/packages/StructuralIdentifiability/gEpsR/src/StructuralIdentifiability.jl:78 [inlined]
 [12] top-level scope
    @ ~/Documents/GitHub/Post-doc_Germany/Project_EpidInference/EpidInference/scripts/Examples_identifiability/identifiability_full.jl:56

If one removes the nonlinear term, then again it works. I tried to looked into AbstractAlgebra.jl a little bit, but couldn't find a workaround.

@JuliaRegistrator register()

Greek letters yield an error, MWE:

using StructuralIdentifiability

ode = @ODEmodel(
    x'(t) = α * x(t),
    y(t) = x(t)
)
assess_identifiability(ode)

This is because they are not supported by Singular. Conversion to Singular is currently used in two places

• factorization, to be eliminated by this issue
• GB computation at the end

How are initial conditions handled?
I can imagine some chemical system where all the initial concentrations are known,
but then only some concentrations are measured during the reaction.
I think identifiability of the reaction kinetics might depend on whether those initial conditions are known or not?
I'm not sure if it is possible to pass that information.

Hi

I want to check StructuralIdentifiability some differential equation which include "exp" but it didnt work with message "Function exp is not supported"

Is there any solutiong for that? or internally not supported?

Thank you
Park

Did the test twice, and got very different results. Rare events?

Here is the model:

using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
S'(t) = -beta_SA * S(t)(A(t) + P(t) ) - beta_SI S(t)In(t),
E'(t) = beta_SA S(t) * (A(t) + P(t) ) + beta_SI * S(t)In(t) - sigma_1E(t),
P'(t) = (1-phi)sigma_1E(t) - sigma_2P(t),
A'(t) = phisigma_1E(t) - M_ar A(t),
In'(t) = sigma_2P(t) - MIn(t),
H'(t) = M*In(t) - kappa_s *H(t),
D'(t) = H(t),
y(t) = D(t)
)

@time println(assess_identifiability(ode))

1st run: all are globally identifiable

Dict{Any, Symbol}(sigma_2 => :globally, M_ar => :globally, kappa_s => :globally, sigma_1 => :globally,
beta_SI => :globally, phi => :globally, M => :globally, beta_SA => :globally)
26381.402992 seconds (214.93 M allocations: 13.735 GiB, 0.04% gc time, 0.31% compilation time)

2nd run: took much longer to finish, and only 2 are globally identifiable

Dict{Any, Symbol}(sigma_2 => :locally, M_ar => :globally, kappa_s => :globally, sigma_1 => :locally,
beta_SI => :locally, phi => :locally, M => :locally, beta_SA => :locally)
67551.049371 seconds (214.83 M allocations: 13.644 GiB, 0.02% gc time, 0.12% compilation time)

I want to convert a ModelingToolkit ODE to the same ODE form as used in StructuralIdentifiabiltiy (and use it for various input). The function is documented, but does not seem to exist?

using StructuralIdentifiability
PreprocessODE
StructuralIdentifiability.PreprocessODE

Both the two bottom lines gives UnderVarError.

GroebnerBasis.jl is now depreciated, so it would be great to move the code to msolve. At the moment, F4 works only over finite fields, so one should either

• refine the sampling bound from the paper (Theorem 3.3 in https://arxiv.org/abs/2111.00991)
• reduce to the zero-dimensional solving

Hello,

i have the following problem. My goal is to get the identifiability of parameters in SBML models. These are imported via SBMLToolkit into ModelingToolkit format. Is there an easy way to run your package on this ModelingToolkit ODE problem? The observables that i have are currently in a Vector of Strings. Currently i am building your ODESystem macro the hard way with the eval() function, based upon a string. But this is quite nasty....

Hi,

I have an ode for the input in my model. See below as an example. Is there any way to define u(t) as an input in such a case? Thanks.

ode = @ODEmodel(
x1'(t) = -(a01 + a21) * x1(t) + a12 * x2(t) + u(t),
x2'(t) = a21 * x1(t) - a12 * x2(t) - x3(t) / b,
x3'(t) = x3(t),
u'(t) = -kd*u(t),
y(t) = x2(t)
)

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

Hello,
I encountered an error while playing around with ModelingToolkit.jl and StructuralIdentifiability.jl. I think StructuralIdentifiability.jl has trouble digesting @variables x(t)[1:2], the supported way of defining tensors in Symbolics.jl.

Here is a minimal working example with StructuralIdentifiability.jl:

ode_1 = @ODEmodel(
    x[1]'(t) = -k1 * x[2],
    x[2]'(t) = -k2 * x[1],
    y[1](t) = x[1](t),
    y[2](t) = x[2](t)
)

This snippet results in ERROR: UndefVarError: x not defined.

I encountered this while using StructuralIdentifiability.jl with ModelingToolkit.jl:

@variables t, x(t)[1:2], y(t)[1:2]
@parameters k1, k2
D = Differential(t)

eqs_1 = [
    D(x[1]) ~ -k1 * x[2],
    D(x[2]) ~ -k2 * x[1]
]

sys_1 = ODESystem(eqs_1, t, name=:example_1)
assess_identifiability(sys_1, measured_quantities=[y[1]~x[1], y[2]~x[2]])

The last line throws

[ Info: Preproccessing `ModelingToolkit.ODESystem` object
ERROR: MethodError: no method matching occursin(::Vector{SymbolicUtils.BasicSymbolic{Real}}, ::Int64)

Closest candidates are:
  occursin(::SymbolicUtils.Symbolic, ::Any)
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/H684H/src/substitute.jl:52
  occursin(::Any)
   @ Base strings/search.jl:713
  occursin(::Any, ::SymbolicUtils.Symbolic)
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/H684H/src/substitute.jl:51

Stacktrace:
 [1] _occursin(needle::Vector{SymbolicUtils.BasicSymbolic{Real}}, haystack::SymbolicUtils.BasicSymbolic{Real})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/H684H/src/substitute.jl:59
 [2] occursin
   @ ~/.julia/packages/SymbolicUtils/H684H/src/substitute.jl:51 [inlined]
 [3] check_variables(dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, iv::Vector{SymbolicUtils.BasicSymbolic{Real}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/FbXPg/src/utils.jl:156
 [4] ODESystem(tag::UInt64, deqs::Vector{Equation}, iv::Vector{SymbolicUtils.BasicSymbolic{Real}}, dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, ps::Vector{Any}, tspan::Nothing, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{ODESystem}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, preface::Nothing, cevents::Vector{ModelingToolkit.SymbolicContinuousCallback}, devents::Vector{ModelingToolkit.SymbolicDiscreteCallback}, metadata::Nothing, gui_metadata::Nothing, tearing_state::Nothing, substitutions::Nothing, complete::Bool, discrete_subsystems::Nothing, unknown_states::Nothing; checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/FbXPg/src/systems/diffeqs/odesystem.jl:152
 [5] ODESystem(deqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Vector{Real}}, dvs::Vector{Any}, ps::OrderedCollections.OrderedSet{Any}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, tspan::Nothing, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Any, Any}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, discrete_events::Nothing, checks::Bool, metadata::Nothing, gui_metadata::Nothing)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/FbXPg/src/systems/diffeqs/odesystem.jl:217
 [6] ODESystem(eqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Vector{Real}}; kwargs::Base.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:name,), Tuple{Symbol}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/FbXPg/src/systems/diffeqs/odesystem.jl:263
 [7] preprocess_ode(de::ODESystem, measured_quantities::Vector{Equation})
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/mGoiQ/src/ODE.jl:458
 [8] assess_identifiability(ode::ODESystem; measured_quantities::Vector{Equation}, funcs_to_check::Vector{Any}, p::Float64)
   @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/mGoiQ/src/StructuralIdentifiability.jl:187
 [9] top-level scope

⌃ [0c5d862f] Symbolics v5.3.1
⌃ [961ee093] ModelingToolkit v8.55.1
  [220ca800] StructuralIdentifiability v0.4.9 

Julia Version 1.9.1
Commit 147bdf428cd (2023-06-07 08:27 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, apple-m1)
  Threads: 4 on 4 virtual cores

SymbolicUtils v1 is a breaking release. Please follow https://github.com/JuliaSymbolics/SymbolicUtils.jl/blob/master/page/upgrade.md to upgrade.

Can we have a new release of StructuralIdentifiability?

Groebner compat is bumped to 0.4; but that is unreleased¹. This results in following incompatibility while building MTK docs:

ERROR: Unsatisfiable requirements detected for package MultivariatePolynomials [102ac46a]:
 MultivariatePolynomials [102ac46a] log:
 ├─possible versions are: 0.2.0-0.5.2 or uninstalled
:
:
 └─restricted by compatibility requirements with Groebner [0b43b601] to versions: 0.4.3-0.4.7 — no versions left
   └─Groebner [0b43b601] log:
     ├─possible versions are: 0.1.0-0.4.4 or uninstalled
     ├─restricted by compatibility requirements with StructuralIdentifiability [220ca800] to versions: 0.2.0-0.3.6
     │ └─StructuralIdentifiability [220ca800] log:
     │   ├─possible versions are: 0.1.0-0.4.10 or uninstalled
     │   └─restricted to versions 0.4 by an explicit requirement, leaving only versions: 0.4.0-0.4.10
     └─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 0.3.0-0.4.4, leaving only versions: 0.3.0-0.3.6
       └─Symbolics [0c5d862f] log: see above

See full log here²

Temporary fix: just restarting the test

also reported here: ederc/GroebnerBasis.jl#42

Current version of tutorials describe only one way of defining a model, using the @ODEmodel macro. A good text about using MTK for this (and maybe even in conjunction with Catalyst) would be great to have.

UPD: Thanks to @saumil-sh - there is such a text in the MTK tutorials but it would be great to have this in the StructuralIdentifiability.jl docs as well

I have a case where I, instead of writing p1 writes 10^p1. The expressions are equivalent, but the rescaling makes sense in some cases (e.g. in the second case the value cannot be negative). However, the second case causes and error, see the following example:

using StructuralIdentifiability, ModelingToolkit

@parameters k1
@variables t LPd(t) ArBr(t) ArR(t)
D = Differential(t)
measured_quantities = [ArR]

@named ode1 = ODESystem([
    D(LPd) ~ -k1*ArBr*LPd,
    D(ArBr) ~ -k1*ArBr*LPd,
    D(ArR) ~ k1*ArBr*LPd
], t)
global_id = assess_identifiability(ode1, measured_quantities=measured_quantities)   # Works fine.
assess_local_identifiability(ode1, measured_quantities=measured_quantities)   # Works fine.

@named ode2 = ODESystem([
    D(LPd) ~ -10^k1*ArBr*LPd,
    D(ArBr) ~ -10^k1*ArBr*LPd,
    D(ArR) ~ 10^k1*ArBr*LPd
], t)
global_id = assess_identifiability(ode2, measured_quantities=measured_quantities)   # Throws an error.
assess_local_identifiability(ode2, measured_quantities=measured_quantities)   # Throws an error.

the error message is:

ERROR: MethodError: no method matching isless(::Int64, ::Nemo.QQMPolyRingElem)

Closest candidates are:
  isless(::Number, ::Unitful.AbstractQuantity)
   @ Unitful ~/.julia/packages/Unitful/orvol/src/quantities.jl:268
  isless(::Integer, ::Nemo.QQFieldElem)
   @ Nemo ~/.julia/packages/Nemo/EuCgH/src/flint/fmpq.jl:367
  isless(::Integer, ::Nemo.ZZRingElem)
   @ Nemo ~/.julia/packages/Nemo/EuCgH/src/flint/fmpz.jl:753
  ...

Stacktrace:
  [1] <(x::Int64, y::Nemo.QQMPolyRingElem)
    @ Base ./operators.jl:343
  [2] <=(x::Int64, y::Nemo.QQMPolyRingElem)
    @ Base ./operators.jl:392
  [3] >=(x::Nemo.QQMPolyRingElem, y::Int64)
    @ Base ./operators.jl:416
  [4] eval_at_nemo(e::SymbolicUtils.BasicSymbolic{Real}, vals::Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/util.jl:348
  [5] (::StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}})(a::SymbolicUtils.BasicSymbolic{Real})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/util.jl:342
  [6] iterate
    @ ./generator.jl:47 [inlined]
  [7] collect_to!(dest::Vector{Int64}, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}}, offs::Int64, st::Int64)
    @ Base ./array.jl:840
  [8] collect_to_with_first!(dest::Vector{Int64}, v1::Int64, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}}, st::Int64)
    @ Base ./array.jl:818
  [9] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:812
 [10] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}})
    @ Base ./array.jl:711
 [11] map(f::Function, A::Vector{Any})
    @ Base ./abstractarray.jl:3261
 [12] eval_at_nemo(e::SymbolicUtils.BasicSymbolic{Real}, vals::Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/util.jl:342
 [13] __preprocess_ode(de::ODESystem, measured_quantities::Vector{Tuple{String, SymbolicUtils.BasicSymbolic{Real}}})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/ODE.jl:549
 [14] preprocess_ode(de::ODESystem, measured_quantities::Vector{Num})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/ODE.jl:465
 [15] assess_local_identifiability(ode::ODESystem; measured_quantities::Vector{Num}, funcs_to_check::Vector{Array}, p::Float64, type::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/local_identifiability.jl:194
 [16] top-level scope
    @ ~/Desktop/Research Projects/New-Active-CRN-Learning/playground.jl:135

julia> ode = @ODEmodel(
                  S'(t) = -β*S(t)*I(t),
                             I'(t) = β*S(t)*I(t) - γ*I(t),
                  R'(t) = γ*I(t),
                  H(t) = 0.05*I(t)
                         )
[ Info: Summary of the model:
[ Info: State variables: I, R, S
[ Info: Parameters: γ, β
[ Info: Inputs:
[ Info: Outputs: H
ERROR: MethodError: no method matching (::Nemo.FmpqMPolyRing)(::Float64)
Closest candidates are:
  (::Nemo.FmpqMPolyRing)() at C:\Users\arno\.julia\packages\Nemo\ITL8k\src\flint\fmpq_mpoly.jl:1170
  (::Nemo.FmpqMPolyRing)(::Vector{Any}, ::Array{Vector{T}, 1}) where T at C:\Users\arno\.julia\packages\Nemo\ITL8k\src\flint\fmpq_mpoly.jl:1233
  (::Nemo.FmpqMPolyRing)(::Rational) at C:\Users\arno\.julia\packages\Nemo\ITL8k\src\flint\fmpq_mpoly.jl:1199
  ...
Stacktrace:
 [1] *(x::Float64, y::Nemo.fmpq_mpoly)
   @ AbstractAlgebra C:\Users\arno\.julia\packages\AbstractAlgebra\Oe2Uj\src\NCRings.jl:82
 [2] top-level scope
   @ REPL[32]:5

5/100 works luckily.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

https://docs.sciml.ai/StructuralIdentifiability/stable/#StructuralIdentifiability.jl on left side navigation bar it says 'chekc' rather than 'check'

Hello, i have the problem that i have a interpolation function on the righthand-side of one of my ODE terms. In order to have an analytical expression for it i used Lagrange polynomials. However, it seem that the StructuralIdentifiability doesn't really work in this case. Here is some example code:

using DifferentialEquations
using StructuralIdentifiability
using ModelingToolkit

days = [0, 2, 4, 6, 8, 10, 12, 14]
surface = [0, 2, 4, 6, 8, 10, 12, 14] 
diam = [0, 2, 4, 6, 8, 10, 12, 14].*2  

@variables t IB(t)  # Define the variables
@parameters k gamma thres pow  # Define the parameters
D = Differential(t)

ode = @ODEmodel(
    IB'(t) = k * sum(
        (
            begin
                (
                    prod((t .- days[1:end .!= i])) / 
                    prod((days[i] .- days[1:end .!= i]))
                ) * diam[i]
            end
        ) ^ pow / (
            begin
                (
                    prod((t .- days[1:end .!= i])) / 
                    prod((days[i] .- days[1:end .!= i]))
                ) * diam[i]
            end
        ) ^ pow + thres ^ pow
        * (
            begin
                (
                    prod((t .- days[1:end .!= i])) / 
                    prod((days[i] .- days[1:end .!= i]))
                ) * surface[i]
            end
        )
        for i in 1:length(days)
    ) - gamma * IB,
    y(t) = IB(t)
)

measured = [IB]

# Assess structural identifiability
identifiability_result = assess_identifiability(ode, measured_quantities=measured)