SciML / SBMLToolkit.jl

SBML differential equation and chemical reaction model (Gillespie simulations) for Julia's SciML ModelingToolkit
https://docs.sciml.ai/SBMLToolkit/stable/
MIT License
39 stars 9 forks source link

Flag to drop units #166

Open oxinabox opened 1 month ago

oxinabox commented 1 month ago

Is your feature request related to a problem? Please describe. I am trying to import a SMBL file. In particular this one

When I do so i get an error:

julia> using OrdinaryDiffEq
       using SBMLToolkit

julia> 
       odesys = readSBML(model_file(), ODESystemImporter())
┌ Error: SBML reported error: The units of the 'math' formula in a <kineticLaw> definition are expected to be the equivalent of _substance per time_.
│  Expected units are mole (exponent = 1, multiplier = 1, scale = -3), second (exponent = -1, multiplier = 1, scale = 0) but the units returned by the <math> expression in the <kineticLaw> (from the <reaction> with id 'vENO') are second (exponent = -1, multiplier = 1, scale = 0), litre (exponent = 1, multiplier = 1, scale = 0).
└ @ SBML ~/.julia/packages/SBML/9eaHi/src/utils.jl:390
┌ Error: SBML reported error: The units of the 'math' formula in a <kineticLaw> definition are expected to be the equivalent of _substance per time_.
│  Expected units are mole (exponent = 1, multiplier = 1, scale = -3), second (exponent = -1, multiplier = 1, scale = 0) but the units returned by the <math> expression in the <kineticLaw> (from the <reaction> with id 'vDAHPS') are mole (exponent = -4.29497e+09, multiplier = 0, scale = 0), second (exponent = -1, multiplier = 1, scale = 0), litre (exponent = -4.29497e+09, multiplier = 1, scale = 0).
└ @ SBML ~/.julia/packages/SBML/9eaHi/src/utils.jl:390
┌ Error: SBML reported error: The units of the 'math' formula in a <kineticLaw> definition are expected to be the equivalent of _substance per time_.
│  Expected units are mole (exponent = 1, multiplier = 1, scale = -3), second (exponent = -1, multiplier = 1, scale = 0) but the units returned by the <math> expression in the <kineticLaw> (from the <reaction> with id 'vPDH') are mole (exponent = -2.14748e+09, multiplier = 0, scale = 0), second (exponent = -1, multiplier = 1, scale = 0), litre (exponent = nan, multiplier = 1, scale = 0).
└ @ SBML ~/.julia/packages/SBML/9eaHi/src/utils.jl:390
ERROR: Setting of level and version did not succeed
Stacktrace:
  [1] get_error_messages(doc::Ptr{…}, error::ErrorException, report_severities::Vector{…}, throw_severities::Vector{…})
    @ SBML ~/.julia/packages/SBML/9eaHi/src/utils.jl:399
  [2] check_errors
    @ ~/.julia/packages/SBML/9eaHi/src/utils.jl:412 [inlined]
  [3] (::SBML.var"#16#17"{Int64, Int64, Vector{String}, Vector{String}})(doc::Ptr{Nothing})
    @ SBML ~/.julia/packages/SBML/9eaHi/src/converters.jl:15
  [4] #50
    @ ~/.julia/packages/SBMLToolkit/LdVcf/src/utils.jl:101 [inlined]
  [5] _readSBML(symbol::Symbol, fn::String, sbml_conversion::SBMLToolkit.var"#50#51", report_severities::Vector{…}, throw_severities::Vector{…})
    @ SBML ~/.julia/packages/SBML/9eaHi/src/readsbml.jl:183
  [6] readSBML(fn::String, sbml_conversion::Function; report_severities::Vector{String}, throw_severities::Vector{String})
    @ SBML ~/.julia/packages/SBML/9eaHi/src/readsbml.jl:229
  [7] readSBML(fn::String, sbml_conversion::Function)
    @ SBML ~/.julia/packages/SBML/9eaHi/src/readsbml.jl:222
  [8] readSBML
    @ ~/.julia/packages/SBMLToolkit/LdVcf/src/systems.jl:19 [inlined]
  [9] #readSBML#1
    @ ~/.julia/packages/SBMLToolkit/LdVcf/src/systems.jl:30 [inlined]
 [10] readSBML
    @ ~/.julia/packages/SBMLToolkit/LdVcf/src/systems.jl:29 [inlined]
 [11] readSBML(sbmlfile::String, ::ODESystemImporter; include_zero_odes::Bool, kwargs::@Kwargs{})
    @ SBMLToolkit ~/.julia/packages/SBMLToolkit/LdVcf/src/systems.jl:42
 [12] readSBML(sbmlfile::String, ::ODESystemImporter)
    @ SBMLToolkit ~/.julia/packages/SBMLToolkit/LdVcf/src/systems.jl:40
 [13] top-level scope
    @ ~/repos/POC_kinetic/test/script.jl:5
Some type information was truncated. Use `show(err)` to see complete types.

Now I know that that file is internally consistent enough to be workable. because it works with MEWpy's load_ODEModel. Which I suspect just discards unit information.

The file quiet likely is a bit messed up. Most bio data is a mess after all

Describe the solution you’d like

I would like a flag to discard units. readSBML(file, ODESystemImporter; discard_units=true)

Describe alternatives We could instead have flag that controls if bad units gives an error, or warns and then falls back to no units, or silently falls back to no units.

paulflang commented 1 month ago

That error looks like it is coming from SBML_jll via SBML.jl

Most bio data is a mess after all

😭 I'm just worried that supporting incorrect models (e.g. by having a discard_units flag) accelerates their proliferation. Imo, the model should be fixed. I cannot access the link tho. Can you send it over again?

oxinabox commented 1 month ago

This one: https://raw.githubusercontent.com/BioSystemsUM/MEWpy/6a5031d6b8c654e15b7b53778d7a8186ad735b23/examples/models/kinetic/chassagnole2002.xml

This file does load up fine with SBMLImporter.jl, so maybe it is something SBMLImporter supports but SBMLToolkit doesn't?

paulflang commented 3 weeks ago

Hi @oxinabox sorry for the delay, swamped with work. I think you need finer handling than the ODESystemImporter does. Can you try if the following produces the right simulation results?

using SBML, SBMLToolkit, OrdinaryDiffEq

fn = "chassagnole2002.xml"
mdl = readSBML(fn, doc -> begin
    set_level_and_version(3, 2, [], [])(doc)
    convert_promotelocals_expandfuns(doc)
end)
rs = ReactionSystem(mdl)
odesys = convert(ODESystem, rs)
odesys = structural_simplify(odesys)

tspan = (0.0, 100.0)
prob = ODEProblem(odesys, [], tspan, [])
sol = solve(prob, Tsit5())

What this does is described here. It's kind of like the flag that you suggested, but a bit more general. If this solves your problem, I will update the README accordingly.

oxinabox commented 3 weeks ago

I had to add using ModelingToolkit but incidence matrix looks (from memory) same as what i saw from SBMLImporter.jl

Even once i change tspan to (0.800) Plot looks pretty different to what MEWpy has (cell 29) in https://github.com/BioSystemsUM/MEWpy/blob/master/examples/03-kinetic.ipynb Where as plot from SBMLImporter.jl looked ideantical Biggest difference is there is something that starts at over 120 and then drops expodentially

Plot of that sol

image

Plot from MEWpy (SBMLImporter near identical but not at fixed timesteps) image

paulflang commented 3 weeks ago

Could you try plotting without cglcex and compare the other species? Reason is that until compartments (and their volumes) are implemented in Catalyst, SBMLToolkit uses and reports only absolute amounts (not concentrations). cglcex is the only species in a compartment with size not equal to 1. This may distort the plot.