From 86f47d5eccd2d9dfcfd03a9ec5314bc6e9794d7a Mon Sep 17 00:00:00 2001 From: paulflang Date: Thu, 4 Jul 2024 00:27:57 +0200 Subject: [PATCH 1/6] omit conversion to mass action for stoichiometries > 100 to avoid StackOverflow in simplify_fractions --- README.md | 39 +++++++++++++++++++++++++++------------ src/reactions.jl | 13 +++++++++++++ test/reactions.jl | 1 + 3 files changed, 41 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index 4d1f83a..355b430 100644 --- a/README.md +++ b/README.md @@ -27,28 +27,43 @@ Pkg.add("SBMLToolkit") SBML models can be simulated with the following steps (note that `sol` is in absolute quantities rather than concentrations): ```julia -using SBMLToolkit, ModelingToolkit, OrdinaryDiffEq +using SBMLToolkit, OrdinaryDiffEq -SBMLToolkit.checksupport_file("my_model.xml") -mdl = readSBML("my_model.xml", doc -> begin - set_level_and_version(3, 2)(doc) - convert_promotelocals_expandfuns(doc) -end) - -rs = ReactionSystem(mdl) # If you want to create a reaction system -odesys = convert(ODESystem, rs) # Alternatively: ODESystem(mdl) +odesys = readSBML("my_model.xml", ODESystemImporter()) tspan = (0.0, 1.0) prob = ODEProblem(odesys, [], tspan, []) sol = solve(prob, Tsit5()) ``` -Alternatively, SBMLToolkit also provides more concise methods to import `SBML.Models`, `Catalyst.ReactionSystems`, and `ModelingToolkit.ODESystems`. +While this imports an `ODESystem` directly, you can also import a Catalyst.jl `ReactionSystem`: ```julia -mdl = readSBML("my_model.xml", DefaultImporter()) +using SBMLToolkit + rs = readSBML("my_model.xml", ReactionSystemImporter()) -odesys = readSBML("my_model.xml", ODESystemImporter()) +``` + +One common case where this is useful is if you want to run stochastic instead of ODE simulations. + +In the very unlikely case that you need fine-grained control over the SBML import, you can create an SBML.jl `Model` (we strongly recommend manually running `SBMLToolkit.checksupport_file("my_model.xml")` before) + +```julia +using SBML + +mdl = readSBML("my_model.xml", doc -> begin + set_level_and_version(3, 2)(doc) + convert_promotelocals_expandfuns(doc) +end) +``` + +The conversion to SBML level 3 version 2 is necessary, because older versions are not well supported in SBMLToolkit. `convert_promotelocals_expandfuns` basically flattens the SBML before the import. Once you have obtained the `Model`, you can convert it to a `ReactionSystem` and `ODESystem`. + +```julia +using SBMLToolkit + +rs = ReactionSystem(mdl) +odesys = convert(ODESystem, rs) ``` ## License diff --git a/src/reactions.jl b/src/reactions.jl index 86e3409..4de5bc6 100644 --- a/src/reactions.jl +++ b/src/reactions.jl @@ -19,6 +19,7 @@ function get_reactions(model::SBML.Model) enforce_rate = enforce_rate) else kl = substitute(symbolic_math, subsdict) # Todo: use SUBSDICT + println(reaction) add_reaction!(rxs, kl, reactant_references, product_references, model) end end @@ -61,7 +62,9 @@ function add_reaction!(rxs::Vector{Reaction}, reactants, products, rstoichvals, pstoichvals = get_reagents(reactant_references, product_references, model) isnothing(reactants) && isnothing(products) && return + println("rstoichvals: $rstoichvals") rstoichvals = stoich_convert_to_ints(rstoichvals) + println("rstoichvals: $rstoichvals") pstoichvals = stoich_convert_to_ints(pstoichvals) kl, our = use_rate(kl, reactants, rstoichvals) our = enforce_rate ? true : our @@ -183,7 +186,17 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing}, rate_const = kl elseif isnothing(reactants) | isnothing(stoich) throw(ErrorException("`reactants` and `stoich` are inconsistent: `reactants` are $(reactants) and `stoich` is $(stoich).")) + elseif max(stoich...) > 100 # simplify_fractions might StackOverflow + rate_const = kl else + println("kl") + println(kl) + println("reactants") + println(reactants) + println(stoich) + println(*((.^(reactants, stoich))...)) + println("fraction") + println(kl / *((.^(reactants, stoich))...)) rate_const = SymbolicUtils.simplify_fractions(kl / *((.^(reactants, stoich))...)) end diff --git a/test/reactions.jl b/test/reactions.jl index 0dafb76..1fb47f3 100644 --- a/test/reactions.jl +++ b/test/reactions.jl @@ -135,6 +135,7 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r)) @test isequal(SBMLToolkit.get_massaction(k1 + c1, nothing, nothing), k1 + c1) # Case zero order kineticLaw @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2 / (c1 + s2), [s1], [1])) # Case Michaelis-Menten kinetics @test isnan(SBMLToolkit.get_massaction(k1 * s1 * IV, [s1], [1])) # Case kineticLaw with time +@test isequal(SBMLToolkit.get_massaction(k1 * (s1 + s1), [s1], [101]), k1 * (s1 + s1)) # Case no simplication performed due to large rstoich @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1])) @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1])) From c8fd129da81ccbc88d2c8fe6a42021d8ec5623b9 Mon Sep 17 00:00:00 2001 From: paulflang Date: Thu, 4 Jul 2024 00:35:03 +0200 Subject: [PATCH 2/6] fix typo --- test/reactions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/reactions.jl b/test/reactions.jl index 1fb47f3..8f1ee9a 100644 --- a/test/reactions.jl +++ b/test/reactions.jl @@ -135,7 +135,7 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r)) @test isequal(SBMLToolkit.get_massaction(k1 + c1, nothing, nothing), k1 + c1) # Case zero order kineticLaw @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2 / (c1 + s2), [s1], [1])) # Case Michaelis-Menten kinetics @test isnan(SBMLToolkit.get_massaction(k1 * s1 * IV, [s1], [1])) # Case kineticLaw with time -@test isequal(SBMLToolkit.get_massaction(k1 * (s1 + s1), [s1], [101]), k1 * (s1 + s1)) # Case no simplication performed due to large rstoich +@test isequal(SBMLToolkit.get_massaction(k1 * (s1 + s1), [s1], [101]), k1 * (s1 + s1)) # Case no simplification performed due to large rstoich @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1])) @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1])) From 938a5bd31d9941b635dcd0db6ceeaadbe9798afc Mon Sep 17 00:00:00 2001 From: paulflang Date: Thu, 4 Jul 2024 00:52:26 +0200 Subject: [PATCH 3/6] improve warning when SBML file contains SpeciesReference with undefined stoichiometry --- src/reactions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactions.jl b/src/reactions.jl index 4de5bc6..f38534a 100644 --- a/src/reactions.jl +++ b/src/reactions.jl @@ -92,7 +92,7 @@ function get_reagents(reactant_references::Vector{SBML.SpeciesReference}, sn = rr.species stoich = rr.stoichiometry if isnothing(stoich) - @warn "Stoichiometries of SpeciesReferences are not defined. Setting to 1." maxlog=1 + @warn "SBML SpeciesReferences does not contain stoichiometries. Assuming stoichiometry of 1." maxlog=1 stoich = 1.0 end iszero(stoich) && @error("Stoichiometry of $sn must be non-zero") From 36f3609dd36b62cfcb17e8af251de164e355bb8b Mon Sep 17 00:00:00 2001 From: paulflang Date: Thu, 4 Jul 2024 01:04:03 +0200 Subject: [PATCH 4/6] fix test case for high stoich --- test/reactions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/reactions.jl b/test/reactions.jl index 8f1ee9a..a39147e 100644 --- a/test/reactions.jl +++ b/test/reactions.jl @@ -135,7 +135,7 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r)) @test isequal(SBMLToolkit.get_massaction(k1 + c1, nothing, nothing), k1 + c1) # Case zero order kineticLaw @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2 / (c1 + s2), [s1], [1])) # Case Michaelis-Menten kinetics @test isnan(SBMLToolkit.get_massaction(k1 * s1 * IV, [s1], [1])) # Case kineticLaw with time -@test isequal(SBMLToolkit.get_massaction(k1 * (s1 + s1), [s1], [101]), k1 * (s1 + s1)) # Case no simplification performed due to large rstoich +@test isnan(SBMLToolkit.get_massaction(k1 * s1, [s1], [101])) # Case no simplification performed due to large rstoich @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1])) @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1])) From 6c2a84bf5f6d15f059e7201a187f237036f6b866 Mon Sep 17 00:00:00 2001 From: paulflang Date: Thu, 4 Jul 2024 01:04:55 +0200 Subject: [PATCH 5/6] clean up reactions.jl --- src/reactions.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/reactions.jl b/src/reactions.jl index f38534a..86ff18c 100644 --- a/src/reactions.jl +++ b/src/reactions.jl @@ -19,7 +19,6 @@ function get_reactions(model::SBML.Model) enforce_rate = enforce_rate) else kl = substitute(symbolic_math, subsdict) # Todo: use SUBSDICT - println(reaction) add_reaction!(rxs, kl, reactant_references, product_references, model) end end @@ -62,9 +61,7 @@ function add_reaction!(rxs::Vector{Reaction}, reactants, products, rstoichvals, pstoichvals = get_reagents(reactant_references, product_references, model) isnothing(reactants) && isnothing(products) && return - println("rstoichvals: $rstoichvals") rstoichvals = stoich_convert_to_ints(rstoichvals) - println("rstoichvals: $rstoichvals") pstoichvals = stoich_convert_to_ints(pstoichvals) kl, our = use_rate(kl, reactants, rstoichvals) our = enforce_rate ? true : our @@ -189,14 +186,6 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing}, elseif max(stoich...) > 100 # simplify_fractions might StackOverflow rate_const = kl else - println("kl") - println(kl) - println("reactants") - println(reactants) - println(stoich) - println(*((.^(reactants, stoich))...)) - println("fraction") - println(kl / *((.^(reactants, stoich))...)) rate_const = SymbolicUtils.simplify_fractions(kl / *((.^(reactants, stoich))...)) end From c909d5aba9375f8bf52917d7485a4c1b4332dce5 Mon Sep 17 00:00:00 2001 From: paulflang Date: Thu, 4 Jul 2024 12:51:25 +0200 Subject: [PATCH 6/6] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d298d37..6ba044c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SBMLToolkit" uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e" authors = ["paulflang", "anandijain"] -version = "0.1.25" +version = "0.1.26" [deps] Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"