Problem with simulating SBML Model

max-de-rooij · February 29, 2024, 10:55am

I am trying to solve an ODE problem from a downloaded SBML file, but I seem to get an error. I used to be able to run this file (although I do not remember the specific versions of each package I used before).

My code is the following:

using DifferentialEquations
using ModelingToolkit
using Plots
using Downloads
using SBMLToolkit


URL = "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000627.3?filename=BIOMD0000000627_url.xml"
output = Downloads.download(URL)
mdl = readSBML(output, doc -> begin
    set_level_and_version(3, 2)(doc)
    convert_promotelocals_expandfuns(doc)
  end)
sys = ODESystem(mdl)
sys = structural_simplify(sys)
prob = ODEProblem(sys, [], (0, 400), [])

sol = solve(prob, alg_hints=[:stiff])

And the error I get is

ERROR: DomainError with -11.600402288322522:
log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
DomainError detected in the user `f` function. This occurs when the domain of a function is violated.
For example, `log(-1.0)` is undefined because `log` of a real number is defined to only output real
numbers, but `log` of a negative number is complex valued and therefore Julia throws a DomainError
by default. Cases to be aware of include:

* `log(x)`, `sqrt(x)`, `cbrt(x)`, etc. where `x<0`
* `x^y` for `x<0` floating point `y` (example: `(-1.0)^(1/2) == im`)

Within the context of SciML, this error can occur within the solver process even if the domain constraint
would not be violated in the solution due to adaptivity. For example, an ODE solver or optimization
routine may check a step at `new_u` which violates the domain constraint, and if violated reject the
step and use a smaller `dt`. However, the throwing of this error will have halted the solving process.

Thus the recommended fix is to replace this function with the equivalent ones from NaNMath.jl
(https://github.com/JuliaMath/NaNMath.jl) which returns a NaN instead of an error. The solver will then
effectively use the NaN within the error control routines to reject the out of bounds step. Additionally,
one could perform a domain transformation on the variables so that such an issue does not occur in the
definition of `f`.

For more information, check out the following FAQ page:
https://docs.sciml.ai/Optimization/stable/API/FAQ/#The-Solver-Seems-to-Violate-Constraints-During-the-Optimization,-Causing-DomainErrors,-What-Can-I-Do-About-That?

Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:301
  [3] log
    @ ./special/log.jl:267 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/SymbolicUtils/r1pzW/src/code.jl:418 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/Symbolics/GQNQy/src/build_function.jl:537 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/SymbolicUtils/r1pzW/src/code.jl:375 [inlined]
  [7] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [8] macro expansion
    @ ./none:0 [inlined]
  [9] generated_callfunc
    @ ./none:0 [inlined]
 [10] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:150
 [11] k
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Gpzyo/src/systems/diffeqs/abstractodesystem.jl:371 [inlined]
 [12] Void
    @ SciMLBase ~/.julia/packages/SciMLBase/Dwomw/src/utils.jl:482 [inlined]
 [13] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::Vector{…}, arg4::Float64)
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
 [14] macro expansion
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
 [15] do_ccall
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
 [16] FunctionWrapper
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
 [17] _call
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined]
 [18] FunctionWrappersWrapper
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined]
 [19] ODEFunction
    @ ~/.julia/packages/SciMLBase/Dwomw/src/scimlfunctions.jl:2168 [inlined]
 [20] commonfun(t::Float64, y::Ptr{Float64}, yp::Ptr{Float64}, comfun::LSODA.CommonFunction{ODEFunction{…}, Vector{…}})
    @ LSODA ~/.julia/packages/LSODA/4dVJc/src/common.jl:16
 [21] lsoda(ctx::lsoda_context_t, y::Vector{Float64}, t::Vector{Float64}, tout::Float64)
    @ LSODA ~/.julia/packages/LSODA/4dVJc/src/types_and_consts.jl:94
 [22] __solve(prob::ODEProblem{…}, alg::lsoda, timeseries::Vector{…}, ts::Vector{…}, ks::Vector{…}; verbose::Bool, abstol::Float64, reltol::Float64, tstops::Vector{…}, saveat::Vector{…}, maxiter::Int64, callback::Nothing, timeseries_errors::Bool, save_everystep::Bool, save_start::Bool, userdata::Nothing, alias_u0::Bool, kwargs::@Kwargs{})
    @ LSODA ~/.julia/packages/LSODA/4dVJc/src/common.jl:181
 [23] __solve(prob::ODEProblem{…}, alg::lsoda, timeseries::Vector{…}, ts::Vector{…}, ks::Vector{…})
    @ LSODA ~/.julia/packages/LSODA/4dVJc/src/common.jl:20
 [24] __solve
    @ LSODA ~/.julia/packages/LSODA/4dVJc/src/common.jl:20 [inlined]
 [25] #solve_call#34
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:612 [inlined]
 [26] solve_call(_prob::ODEProblem{…}, args::lsoda)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:569
 [27] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::Vector{…}, args::lsoda; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1064
 [28] solve_up
    @ ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1050 [inlined]
 [29] solve(prob::ODEProblem{…}, args::lsoda; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:987
 [30] solve(prob::ODEProblem{…}, args::lsoda)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:977
 [31] top-level scope
    @ ~/TUe/8DBM050-development/Examples/Winter_2017.jl:25
Some type information was truncated. Use `show(err)` to see complete types.

ERROR: DomainError with -0.6708653124046805:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
  [1] throw_exp_domainerror(x::Float64)
    @ Base.Math ./math.jl:41
  [2] ^
    @ Base.Math ./math.jl:1206 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/r1pzW/src/code.jl:375 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [5] macro expansion
    @ ./none:0 [inlined]
  [6] generated_callfunc
    @ ./none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:150
  [8] k
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Gpzyo/src/systems/diffeqs/abstractodesystem.jl:371 [inlined]
  [9] Void
    @ SciMLBase ~/.julia/packages/SciMLBase/Dwomw/src/utils.jl:482 [inlined]
 [10] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::Vector{…}, arg4::Float64)
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
 [11] macro expansion
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
 [12] do_ccall
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
 [13] FunctionWrapper
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
 [14] _call
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined]
 [15] FunctionWrappersWrapper
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined]
 [16] Void
    @ SciMLBase ~/.julia/packages/SciMLBase/Dwomw/src/utils.jl:482 [inlined]
 [17] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::Vector{…}, arg4::Float64)
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
 [18] macro expansion
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
 [19] do_ccall
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
 [20] FunctionWrapper
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
 [21] _call
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined]
 [22] FunctionWrappersWrapper
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined]
 [23] ODEFunction
    @ ~/.julia/packages/SciMLBase/Dwomw/src/scimlfunctions.jl:2168 [inlined]
 [24] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Rosenbrock5Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uQcyc/src/perform_step/rosenbrock_perform_step.jl:2560
 [25] perform_step!
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uQcyc/src/perform_step/rosenbrock_perform_step.jl:2427 [inlined]
 [26] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uQcyc/src/solve.jl:543
 [27] __solve(::ODEProblem{…}, ::TRBDF2{…}; kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uQcyc/src/solve.jl:7 [inlined]
 [28] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/uQcyc/src/solve.jl:1 [inlined]
 [29] solve_call(_prob::ODEProblem{…}, args::Rodas5P{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:612
 [30] solve_call
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:569 [inlined]
 [31] #solve_up#42
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1064 [inlined]
 [32] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1050 [inlined]
 [33] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:987 [inlined]
 [34] __solve(::ODEProblem{…}, ::Nothing; default_set::Bool, kwargs::@Kwargs{…})
    @ DifferentialEquations ~/.julia/packages/DifferentialEquations/v7yof/src/default_solve.jl:14
 [35] __solve
    @ DifferentialEquations ~/.julia/packages/DifferentialEquations/v7yof/src/default_solve.jl:1 [inlined]
 [36] #__solve#61
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1378 [inlined]
 [37] __solve
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1370 [inlined]
 [38] #solve_call#34
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:612 [inlined]
 [39] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::Vector{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1056
 [40] solve_up
    @ ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:1050 [inlined]
 [41] solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/solve.jl:987
 [42] top-level scope
    @ ~/TUe/8DBM050-development/Examples/Winter_2017.jl:18
Some type information was truncated. Use `show(err)` to see complete types.

I am able to solve the system using the TRBDF2() solver, but then my results don’t look at all like the ones I get when using COPASI. Was there an update in the past year that could have caused this issue? And would anyone have any pointers to help?

ChrisRackauckas · February 29, 2024, 12:51pm

What about with reduced tolerances?

isaacsas · February 29, 2024, 1:41pm

You can also try GitHub - sebapersson/SBMLImporter.jl: Import dynamic models in the SBML format into a ReactionSystem for Gillespe, SDE and ODE simulations, I think it is more actively developed at this point than SBMLToolkit.

max-de-rooij · March 7, 2024, 10:56am

Seems that the initial conditions are loaded incorrectly when using SBMLToolkit, but unfortunately SBMLImporter gives me an error. I use the following code to load the model:

URL = "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000627.3?filename=BIOMD0000000627_url.xml"
output = Downloads.download(URL)
mdl, cb = load_SBML(output)

But this gives me the error:

ERROR: UndefVarError: `F_in` not defined
Stacktrace:
 [1] top-level scope
   @ none:1

Where F_in is one of the variables in the model.

isaacsas · March 7, 2024, 12:10pm

I’d suggest opening issues on the SBMLToolkit and/or SBMLImporter Github sites. The developers for both packages aren’t that active on Discourse, but will likely see issues and be able to help on Github.