From cdee4d8be79891b3f732cbdcfc2cab3fbd51246f Mon Sep 17 00:00:00 2001 From: = <=> Date: Mon, 9 Aug 2021 10:39:07 -0400 Subject: [PATCH] separate out reactants and products for new SBML version --- Project.toml | 6 ++--- src/reactionsystem.jl | 52 ++++++++++++++++++++++++------------------ test/reactionsystem.jl | 18 +++++++-------- 3 files changed, 41 insertions(+), 35 deletions(-) diff --git a/Project.toml b/Project.toml index 56c7fb5..764e929 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.5" +version = "0.1.6" [deps] Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" @@ -11,10 +11,10 @@ SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" [compat] -Catalyst = "6" +Catalyst = "6, 7" MathML = "0.1.7" ModelingToolkit = "5" -SBML = "0.6" +SBML = "0.7" SymbolicUtils = "0.1 - 0.13" julia = "1.6" diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 21c8f3e..69d2f94 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -84,22 +84,25 @@ function mtk_reactions(model::SBML.Model) handle_empty_compartment_size = _ -> 1.0) symbolic_math = Num(convert(Num, extensive_math, convert_time = (x::SBML.MathTime) -> Catalyst.DEFAULT_IV)) + + rstoich = reaction.reactants + pstoich = reaction.products if reaction.reversible symbolic_math = getunidirectionalcomponents(symbolic_math) kl_fw, kl_rv = [substitute(x, subsdict) for x in symbolic_math] - reagents = getreagents(reaction.stoichiometry, model) + reagents = getreagents(rstoich, pstoich, model) isnothing(reagents[1]) && isnothing(reagents[2]) && continue kl_fw, our = use_rate(kl_fw, reagents[1], reagents[3]) kl_rv = from_noncombinatoric(kl_rv, reagents[3], our) push!(rxs, ModelingToolkit.Reaction(kl_fw, reagents...; only_use_rate=our)) - reagents = getreagents(reaction.stoichiometry, model; rev=true) + reagents = getreagents(rstoich, pstoich, model; rev=true) kl_rv, our = use_rate(kl_rv, reagents[1], reagents[3]) kl_rv = from_noncombinatoric(kl_rv, reagents[3], our) push!(rxs, ModelingToolkit.Reaction(kl_rv, reagents...; only_use_rate=our)) else kl = substitute(symbolic_math, subsdict) - reagents = getreagents(reaction.stoichiometry, model) + reagents = getreagents(rstoich, pstoich, model) isnothing(reagents[1]) && isnothing(reagents[2]) && continue kl, our = use_rate(kl, reagents[1], reagents[3]) kl = from_noncombinatoric(kl, reagents[3], our) @@ -134,31 +137,36 @@ function use_rate(kl::Num, react::Union{Vector{Num},Nothing}, stoich::Union{Vect end """ Get reagents """ -function getreagents(stoich::Dict{String,<:Real}, model::SBML.Model; rev=false) - if rev - stoich = Dict(k => -v for (k,v) in stoich) - end +function getreagents(rstoichdict::Dict{String,<:Real}, pstoichdict::Dict{String,<:Real}, model::SBML.Model; rev=false) reactants = Num[] products = Num[] rstoich = Float64[] pstoich = Float64[] - for (k,v) in stoich - if v < 0 - push!(reactants, create_var(k,Catalyst.DEFAULT_IV)) - push!(rstoich, -v) - if model.species[k].boundary_condition == true - push!(products, create_var(k,Catalyst.DEFAULT_IV)) - push!(pstoich, -v) - end - elseif v > 0 - if model.species[k].boundary_condition != true - push!(products, create_var(k,Catalyst.DEFAULT_IV)) - push!(pstoich, v) - end - else - @error("Stoichiometry of $k must be non-zero") + + if rev + tmp = rstoichdict + rstoichdict = pstoichdict + pstoichdict = tmp + end + + for (k,v) in rstoichdict + iszero(v) && @error("Stoichiometry of $k must be non-zero") + push!(reactants, create_var(k,Catalyst.DEFAULT_IV)) + push!(rstoich, v) + if model.species[k].boundary_condition == true + push!(products, create_var(k,Catalyst.DEFAULT_IV)) + push!(pstoich, v) + end + end + + for (k,v) in pstoichdict + iszero(v) && @error("Stoichiometry of $k must be non-zero") + if model.species[k].boundary_condition != true + push!(products, create_var(k,Catalyst.DEFAULT_IV)) + push!(pstoich, v) end end + if (length(reactants)==0) reactants = nothing; rstoich = nothing end if (length(products)==0) products = nothing; pstoich = nothing end (reactants, products, rstoich, pstoich) diff --git a/test/reactionsystem.jl b/test/reactionsystem.jl index ab0f47b..5a972c6 100644 --- a/test/reactionsystem.jl +++ b/test/reactionsystem.jl @@ -9,11 +9,11 @@ KINETICMATH1 = SBML.MathIdent("k1") KINETICMATH2 = SBML.MathApply("*", SBML.Math[ SBML.MathIdent("k1"), SBML.MathIdent("s2")]) KINETICMATH3 = SBML.MathApply("-", SBML.Math[KINETICMATH2, KINETICMATH1]) -REACTION1 = SBML.Reaction(Dict("s1" => 1), (NaN, ""), (NaN, ""), NaN, +REACTION1 = SBML.Reaction(Dict(), Dict("s1" => 1), (NaN, ""), (NaN, ""), NaN, nothing, KINETICMATH1, false) -REACTION2 = SBML.Reaction(Dict("s2" => -1), (NaN, ""), (NaN, ""), NaN, +REACTION2 = SBML.Reaction(Dict("s2" => 1), Dict(), (NaN, ""), (NaN, ""), NaN, nothing, KINETICMATH2, false) -REACTION3 = SBML.Reaction(Dict("s2" => -1), (NaN, ""), (NaN, ""), NaN, +REACTION3 = SBML.Reaction(Dict("s2" => 1), Dict(), (NaN, ""), (NaN, ""), NaN, nothing, KINETICMATH3, true) MODEL1 = SBML.Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("s1" => SPECIES1), Dict("r1" => REACTION1), Dict(), Dict()) # PL: For instance in the compartments dict, we may want to enforce that key and compartment.name are identical MODEL2 = SBML.Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("s2" => SPECIES2), Dict("r2" => REACTION2), Dict(), Dict()) @@ -111,8 +111,7 @@ truereaction = ModelingToolkit.Reaction(k1, nothing, [s1], nothing, [1]; only_us @test isequal(reaction, truereaction) km = SBML.MathTime("x") -reac = SBML.Reaction(Dict("s1" => 1), (NaN, ""), (NaN, ""), NaN, -nothing, km, false) +reac = SBML.Reaction(Dict(), Dict("s1" => 1), (NaN, ""), (NaN, ""), NaN, nothing, km, false) mod = SBML.Model(Dict(), Dict(), Dict("c1" => COMP1), Dict("s1" => SPECIES1), Dict("r1" => reac), Dict(), Dict()) @test isequal(Catalyst.DEFAULT_IV, SBMLToolkit.mtk_reactions(mod)[1].rate) @@ -126,7 +125,7 @@ mod = SBML.Model(Dict(), Dict(), Dict("c1" => COMP1), Dict("s1" => SPECIES1), Di @test isequal(SBMLToolkit.use_rate(k1*s1*s2/(c1+s2), [s1], [1]), (k1*s1*s2/(c1+s2), true)) # Case Michaelis-Menten kinetics # Test getreagents -@test isequal((nothing, [s1], nothing, [1.]), SBMLToolkit.getreagents(REACTION1.stoichiometry, MODEL1)) +@test isequal((nothing, [s1], nothing, [1.]), SBMLToolkit.getreagents(REACTION1.reactants, REACTION1.products, MODEL1)) constspec = SBML.Species("constspec", "c1", true, nothing, nothing, (1., "substance"), nothing, true) ncs = SBMLToolkit.create_var("constspec",Catalyst.DEFAULT_IV) @@ -135,11 +134,10 @@ SBML.MathApply("*", SBML.Math[ SBML.MathIdent("k1"), SBML.MathIdent("constspec")]), SBML.MathIdent("k1")]) -constreac = SBML.Reaction(Dict("constspec" => -1), (NaN, ""), (NaN, ""), NaN, -nothing, kineticmath, false) +constreac = SBML.Reaction(Dict("constspec" => 1), Dict(), (NaN, ""), (NaN, ""), NaN, nothing, kineticmath, false) constmod = SBML.Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("constspec" => constspec), Dict("r1" => constreac), Dict(), Dict()) -@test isequal(([ncs], [ncs], [1.], [1.]), SBMLToolkit.getreagents(constreac.stoichiometry, constmod)) -@test isequal((nothing, nothing, nothing, nothing), SBMLToolkit.getreagents(constreac.stoichiometry, constmod; rev=true)) +@test isequal(([ncs], [ncs], [1.], [1.]), SBMLToolkit.getreagents(constreac.reactants, constreac.products, constmod)) +@test isequal((nothing, nothing, nothing, nothing), SBMLToolkit.getreagents(constreac.reactants, constreac.products, constmod; rev=true)) # Test getunidirectionalcomponents km = SBML.MathApply("-", SBML.Math[KINETICMATH1, SBML.MathIdent("c1")])