Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

Identifiability of ModelingToolkit models #128

Closed
MaAl13 opened this issue Nov 8, 2022 · 16 comments
Closed

Identifiability of ModelingToolkit models #128

MaAl13 opened this issue Nov 8, 2022 · 16 comments

Comments

@MaAl13
Copy link

MaAl13 commented Nov 8, 2022

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....

@iliailmer
Copy link
Collaborator

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.

@MaAl13
Copy link
Author

MaAl13 commented Nov 8, 2022

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)

@iliailmer
Copy link
Collaborator

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)

@MaAl13
Copy link
Author

MaAl13 commented Nov 9, 2022

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

@iliailmer
Copy link
Collaborator

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.

@ChrisRackauckas
Copy link
Member

Do structural_simplify(odesys_SBML) to eliminate the algebraic equations.

@iliailmer
Copy link
Collaborator

iliailmer commented Nov 10, 2022

@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.

@MaAl13
Copy link
Author

MaAl13 commented Nov 14, 2022

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

@iliailmer
Copy link
Collaborator

@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

@iliailmer
Copy link
Collaborator

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?

@MaAl13
Copy link
Author

MaAl13 commented Nov 28, 2022

Ah, thank you for pointing out my error!

@pogudingleb
Copy link
Collaborator

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

@MaAl13
Copy link
Author

MaAl13 commented Feb 7, 2023

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

@pogudingleb pogudingleb reopened this Feb 17, 2023
@pogudingleb
Copy link
Collaborator

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.

@pogudingleb
Copy link
Collaborator

@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.

@MaAl13
Copy link
Author

MaAl13 commented Apr 17, 2023

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 :)

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

No branches or pull requests

4 participants