Identifiability of ModelingToolkit models about structuralidentifiability.jl HOT 16 CLOSED

Hello Gleb, thanks for trying to target my problem. Indeed it seems like a complicated problem, would be nice if it would work some day :)

from structuralidentifiability.jl.

Hi @MaAl13 , thank you for the question! Can you share an example of the model you're working with? The code should support the ODE problems defined by ModelingToolkit generally.

from structuralidentifiability.jl.

Hello Ilia, thank you so much for your response!!
Maybe it is trivial for you, basically i thought it would be nice to do sth like below. I want to import the model with the SBMLToolkit and put it in the StructuralIdentifiability/SIAN. But for me it is unclear how to add the observables to the ODE system and if the ModelingToolkit system is okay as an input for assess_identifiability.

using SBMLToolkit
using StructuralIdentifiability
mdl = readSBML("Wodarz2018_1.xml", doc -> begin
            set_level_and_version(3, 1)(doc) # c.f. Format on BioModels homepage, i.e. SBML (L2V1) 
            convert_simplify_math(doc)
        end)

rs = ReactionSystem(mdl)
odesys_SBML = convert(ODESystem, rs)

measured_quantities = string.(collect(ModelingToolkit.get_states(odesys_SBML)))[1:2]

assess_identifiability(odesys_SBML, measured_quantities=measured_quantities)

from structuralidentifiability.jl.

Thank you! One that is expected by assess_identifiability is that measured_quantities is passed as an array of Equation objects. For example, if your measured quantity is x1(t), you would pass [y1~x1].

I could not run the code without the XML file, but try something like this:

using ModelingToolkit

y_functions = @variables t y1(t) y2(t)
y_functions = y_functions[2:end]
measured_quantities = collect(ModelingToolkit.get_states(odesys_SBML))[1:2]

measured_quantities = [y_functions[i] ~ measured_quantities[i] for i in 1:2]
assess_identifiability(odesys_SBML, measured_quantities=measured_quantities)

from structuralidentifiability.jl.

Thank you very much for your input! It definietly makes sene to add the Array of Equations, I modified the code with your suggestions, however there still occurs an error.
The code is now the following

using SBMLToolkit
using StructuralIdentifiability
using ModelingToolkit
mdl = readSBML("Wodarz2018_1.xml", doc -> begin
            set_level_and_version(3, 1)(doc) # c.f. Format on BioModels homepage, i.e. SBML (L2V1) 
            convert_simplify_math(doc)
        end)

rs = ReactionSystem(mdl)
odesys_SBML = convert(ODESystem, rs)

y_functions = @variables t y1(t) y2(t)
y_functions = y_functions[2:end]

measured_quantities = string.(collect(ModelingToolkit.get_states(odesys_SBML)))[1:2]

measured_quantities = [y_functions[i] ~ measured_quantities[i] for i in 1:2]
assess_identifiability(odesys_SBML, measured_quantities=measured_quantities)

This leads to the error
[ Info: Preproccessing ModelingToolkit.ODESystem object
ERROR: MethodError: no method matching getmetadata(::Int64, ::Type{ModelingToolkit.VariableOutput}, ::Nothing)
Closest candidates are:
getmetadata(::SymbolicUtils.Symbolic, ::Any, ::Any) at ~/.julia/packages/SymbolicUtils/qulQp/src/types.jl:28
getmetadata(::SymbolicUtils.Symbolic, ::Any) at ~/.julia/packages/SymbolicUtils/qulQp/src/types.jl:19
getmetadata(::Complex{Num}, ::Any) at ~/.julia/packages/Symbolics/FGTCH/src/Symbolics.jl:162
...
Stacktrace:
[1] isvarkind(m::Type, x::Int64)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/iWH2V/src/variables.jl:23
[2] isoutput(x::Int64)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/iWH2V/src/variables.jl:33
[3] (::StructuralIdentifiability.var"#151#157")(eq::Equation)
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/egNv7/src/ODE.jl:378
[4] filter(f::StructuralIdentifiability.var"#151#157", a::Vector{Equation})
@ Base ./array.jl:2484
[5] PreprocessODE(de::ODESystem, measured_quantities::Vector{Equation})
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/egNv7/src/ODE.jl:378
[6] assess_identifiability(ode::ODESystem; measured_quantities::Vector{Equation}, funcs_to_check::Vector{Any}, p::Float64)
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/egNv7/src/StructuralIdentifiability.jl:166
[7] top-level scope
@ ~/Documents/julia_GPU_solver/SBML/Test_Identifiability.jl:18

By the way, the xml file is downloaded from SBML repository found under the following link, if you want to execute the code :)
https://www.ebi.ac.uk/biomodels/MODEL1908010003#Files

from structuralidentifiability.jl.

Thanks for linking to the model. I checked the equations:

ModelingToolkit.equations(odesys_SBML)

returns

5-element Vector{Equation}:
 Differential(t)(S(t)) ~ D(t)*g(t) + (2.0p(t) - 1.0)*S(t)*r(t)
 Differential(t)(D(t)) ~ 2.0(1.0 - p(t))*S(t)*r(t) - a*D(t) - D(t)*g(t)
 0 ~ p_0 / (1.0 + h2*((D(t) / tme)^k2)) - p(t)
 0 ~ g_0 / (1.0 + h3*((S(t) / tme)^k3)) - g(t)
 0 ~ r_0 / (1.0 + h1*((D(t) / tme)^k1)) - r(t)

The last three are causing the errors issues, since these are not ODEs and assess_identifiability does not know how to deal with them.

I think plugging the last three equations into the first two would work, but this might need to be done manually.

from structuralidentifiability.jl.

Do structural_simplify(odesys_SBML) to eliminate the algebraic equations.

from structuralidentifiability.jl.

@ChrisRackauckas thank you! This works, however, it appears the eval_at_nemo does not detect Float coefficients in this particular ODE, I will look into it.

from structuralidentifiability.jl.

Thank you both! Indeed i also get an error if i use the structural_simplify beforehand. Right now i am checking if the system contains no algebraic equations with the condition that the right hand side of each equation has a differential at the beginning. I use now a pure ODE system but i am still running in an error

using SBMLToolkit
using StructuralIdentifiability
using ModelingToolkit
mdl = readSBML("Varusai2018.xml", doc -> begin
            set_level_and_version(2, 4)(doc) # c.f. Format on BioModels homepage, i.e. SBML (L2V1) 
            convert_simplify_math(doc)
        end)

rs = ReactionSystem(mdl)
odesys_SBML = convert(ODESystem, rs)

#odesys_SBML_new = structural_simplify(odesys_SBML)
eqs = ModelingToolkit.equations(odesys_SBML)
num_real_diff_eqs = 0
for i in 1:length(eqs)
    global num_real_diff_eqs
    if length(split(string(eqs[i]), "~")[1]) > 12
        if "Differential" == split(string(eqs[i]), "~")[1][1:12]
            num_real_diff_eqs+=1
        end
    end
end


if num_real_diff_eqs == length(eqs)
    string_y_func = "y_functions = @variables t"
    for i in 1:num_real_diff_eqs
        global string_y_func
        string_y_func = string_y_func * " y$i(t)"
    end
    eval(Meta.parse(string_y_func))
    #y_functions = @variables t y1(t) y2(t)
    #y_functions = y_functions[2:end]

    measured_quantities = string.(collect(ModelingToolkit.get_states(odesys_SBML)))[1:num_real_diff_eqs]

    measured_quantities = [y_functions[i] ~ measured_quantities[i] for i in 1:num_real_diff_eqs]
    assess_identifiability(odesys_SBML, measured_quantities=measured_quantities)
else
    println("System is not a pure ODE and contains algebraic equations!")
end

[ Info: Preproccessing ModelingToolkit.ODESystem object
ERROR: MethodError: Cannot convert an object of type
Nothing to an object of type
Union{AbstractAlgebra.Generic.Frac{Nemo.fmpq_mpoly}, Nemo.fmpq_mpoly}
Closest candidates are:
convert(::Type{T}, ::T) where T at /Applications/Julia-1.7.3.app/Contents/Resources/julia/share/julia/base/essentials.jl:218
Stacktrace:
[1] setindex!(h::Dict{Nemo.fmpq_mpoly, Union{AbstractAlgebra.Generic.Frac{Nemo.fmpq_mpoly}, Nemo.fmpq_mpoly}}, v0::Nothing, key::Nemo.fmpq_mpoly)
@ Base ./dict.jl:381
[2] PreprocessODE(de::ODESystem, measured_quantities::Vector{Equation})
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/egNv7/src/ODE.jl:402
[3] assess_identifiability(ode::ODESystem; measured_quantities::Vector{Equation}, funcs_to_check::Vector{Any}, p::Float64)
@ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/egNv7/src/StructuralIdentifiability.jl:166
[4] top-level scope
@ ~/Documents/julia_GPU_solver/SBML/Test_Identifiability.jl:36

The file can be downloaded at https://www.ebi.ac.uk/biomodels/BIOMD0000000823#Files

from structuralidentifiability.jl.

@MaAl13 for the code you attached the issue was this line:

measured_quantities = string.(collect(ModelingToolkit.get_states(odesys_SBML)))[1:num_real_diff_eqs]

The string conversion prevented the StructuralIdentifiability code from properly parsing the right-hand side of output equation. removing the string conversion works:

using SBMLToolkit
using StructuralIdentifiability
using ModelingToolkit
mdl = readSBML("./Varusai2018.xml", doc -> begin
    set_level_and_version(2, 4)(doc) # c.f. Format on BioModels homepage, i.e. SBML (L2V1)
    convert_simplify_math(doc)
end)

rs = ReactionSystem(mdl)
odesys_SBML = convert(ODESystem, rs)

#odesys_SBML_new = structural_simplify(odesys_SBML)
eqs = ModelingToolkit.equations(odesys_SBML)
num_real_diff_eqs = 0
for i in 1:length(eqs)
    global num_real_diff_eqs
    if length(split(string(eqs[i]), "~")[1]) > 12
        if "Differential" == split(string(eqs[i]), "~")[1][1:12]
            num_real_diff_eqs += 1
        end
    end
end


if num_real_diff_eqs == length(eqs)
    string_y_func = "y_functions = @variables t"
    for i in 1:num_real_diff_eqs
        global string_y_func
        string_y_func = string_y_func * " y$i(t)"
    end
    eval(Meta.parse(string_y_func))
    #y_functions = @variables t y1(t) y2(t)
    #y_functions = y_functions[2:end]

    measured_quantities = collect(ModelingToolkit.get_states(odesys_SBML))

    measured_quantities = [y_functions[i+1] ~ measured_quantities[i] for i in 1:num_real_diff_eqs]
    assess_identifiability(odesys_SBML; measured_quantities=measured_quantities)
else
    println("System is not a pure ODE and contains algebraic equations!")
end

from structuralidentifiability.jl.

Thank you very much for your input! It definietly makes sene to add the Array of Equations, I modified the code with your suggestions, however there still occurs an error. The code is now the following

using SBMLToolkit
using StructuralIdentifiability
using ModelingToolkit
mdl = readSBML("Wodarz2018_1.xml", doc -> begin
            set_level_and_version(3, 1)(doc) # c.f. Format on BioModels homepage, i.e. SBML (L2V1) 
            convert_simplify_math(doc)
        end)

rs = ReactionSystem(mdl)
odesys_SBML = convert(ODESystem, rs)

y_functions = @variables t y1(t) y2(t)
y_functions = y_functions[2:end]

measured_quantities = string.(collect(ModelingToolkit.get_states(odesys_SBML)))[1:2]

measured_quantities = [y_functions[i] ~ measured_quantities[i] for i in 1:2]
assess_identifiability(odesys_SBML, measured_quantities=measured_quantities)

In this model, after structural_simplify there are functions in the right-hand side that cause issues: p(t), r(t), g(t):

 Differential(t)(S(t)) ~ D(t)*g(t) + (2.0p(t) - 1.0)*S(t)*r(t)
 Differential(t)(D(t)) ~ (2.0 - 2.0p(t))*S(t)*r(t) - a*D(t) - D(t)*g(t)

StructuralIdentifiability does not support variable coefficients in ODE so they would be treated as input variables. Are these the time-varying parameters in the model?

from structuralidentifiability.jl.

Ah, thank you for pointing out my error!

from structuralidentifiability.jl.

@MaAl13 @iliailmer
It seems like the question has been resolved. If not, feel free to reopen.

from structuralidentifiability.jl.

So i have again a problem when trying to run an SBML model. This time the whole program is killed, can you tell me what is causing this error and how i can catch it? I am going through the SBML databse and search for appropriate models for Structural identifiability. Also the warning/errors for float seems to come up, is there a way to make them work aswell?

using SBMLToolkit
using StructuralIdentifiability
using ModelingToolkit
mdl = readSBML("models/BIOMD0000000070_url.xml", doc -> begin
    set_level_and_version(2, 4)(doc) # c.f. Format on BioModels homepage, i.e. SBML (L2V1)
    convert_simplify_math(doc)
end)

rs = ReactionSystem(mdl)
odesys_SBML = convert(ODESystem, rs)

#odesys_SBML_new = structural_simplify(odesys_SBML)
eqs = ModelingToolkit.equations(odesys_SBML)
num_real_diff_eqs = 0
for i in 1:length(eqs)
    global num_real_diff_eqs
    if length(split(string(eqs[i]), "~")[1]) > 12
        if "Differential" == split(string(eqs[i]), "~")[1][1:12]
            num_real_diff_eqs += 1
        end
    end
end


if num_real_diff_eqs == length(eqs)
    string_y_func = "y_functions = @variables t"
    for i in 1:num_real_diff_eqs
        global string_y_func
        string_y_func = string_y_func * " y$i(t)"
    end
    eval(Meta.parse(string_y_func))
    #y_functions = @variables t y1(t) y2(t)
    #y_functions = y_functions[2:end]

    measured_quantities = collect(ModelingToolkit.get_states(odesys_SBML))

    measured_quantities = [y_functions[i+1] ~ measured_quantities[i] for i in 1:num_real_diff_eqs]
    assess_identifiability(odesys_SBML; measured_quantities=measured_quantities)
else
    println("System is not a pure ODE and contains algebraic equations!")
end

[ Info: Preproccessing ModelingToolkit.ODESystem object
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 1.55 will be converted to 31//20.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 2.0 will be converted to 2.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 2.0 will be converted to 2.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 2.0 will be converted to 2.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 2.0 will be converted to 2.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 2.0 will be converted to 2.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 4.0 will be converted to 4.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
┌ Warning: Floating points are not allowed, value 1.55 will be converted to 31//20.
└ @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/m9hcN/src/util.jl:323
Killed

The model causing this error can be found here
https://www.ebi.ac.uk/biomodels/BIOMD0000000070

from structuralidentifiability.jl.

Some preliminary results: it runs out of memory during converting a model from MTK to the internal format. The expressions do not look scary at all, I am looking into it.

from structuralidentifiability.jl.

@MaAl13
Here is what happens with the model. The right-hand sides are expressions involving several levels of nested fractions. StructuralIdentifiability unfolds this to the form of rational function (that is, polynomial divided by a polynomial) and, because of the depth of nestedness, this blows up the memory. I do not have a good solution for now (even if the conversion was finished, I do not think the further computations would go far). It would be good to have the local identifiability part being able to work directly on the expressions using symbolic AD but this is a long story, I will create an issue for this.

from structuralidentifiability.jl.