From e664b13710c12e95eaf923992bdfceefcb62056d Mon Sep 17 00:00:00 2001 From: = <=> Date: Fri, 6 Aug 2021 13:00:30 -0400 Subject: [PATCH 1/4] PR to sync with latest breaking SBML change to Reaction --- src/reactionsystem.jl | 55 +++++++++++++++++++++--------------------- test/reactionsystem.jl | 17 ++++++------- 2 files changed, 35 insertions(+), 37 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 21c8f3e..8cf6610 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) 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,36 +137,32 @@ 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 - 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") - end - end +function getreagents(rstoichdict::Dict{String,<:Real}, pstoichdict::Dict{String,<:Real}, model::SBML.Model) + + reactants, rstoich = stoichmap_to_varmap(rstoichdict) + products, pstoich = stoichmap_to_varmap(pstoichdict) + if (length(reactants)==0) reactants = nothing; rstoich = nothing end if (length(products)==0) products = nothing; pstoich = nothing end (reactants, products, rstoich, pstoich) end +function stoichmap_to_varmap(stoich::Dict{String,<:Real}) + species = Num[] + stoichs = Float64[] + for (k,v) in stoich + if iszero(v) + @error("Stoichiometry of $k must be non-zero") + else + # if model.species[k].boundary_condition == true + push!(species, create_var(k,Catalyst.DEFAULT_IV)) + push!(stoichs, v) + # end + end + end + (species, stoichs) +end + """ Infer forward and reverse components of bidirectional kineticLaw """ function getunidirectionalcomponents(bidirectional_math) err = "Cannot separate bidirectional kineticLaw `$bidirectional_math` to forward and reverse part. Please make reaction irreversible or rearrange kineticLaw to the form `term1 - term2`." diff --git a/test/reactionsystem.jl b/test/reactionsystem.jl index ab0f47b..9bb325e 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("s1" => 1), Dict(), (NaN, ""), (NaN, ""), NaN, nothing, KINETICMATH1, false) -REACTION2 = SBML.Reaction(Dict("s2" => -1), (NaN, ""), (NaN, ""), NaN, +REACTION2 = SBML.Reaction(Dict(), Dict("s2" => 1), (NaN, ""), (NaN, ""), NaN, nothing, KINETICMATH2, false) -REACTION3 = SBML.Reaction(Dict("s2" => -1), (NaN, ""), (NaN, ""), NaN, +REACTION3 = SBML.Reaction(Dict(), Dict("s2" => 1), (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,7 +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, +reac = SBML.Reaction(Dict("s1" => 1), Dict(), (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 +126,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 +135,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(), Dict("constspec" => 1), (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)) # Test getunidirectionalcomponents km = SBML.MathApply("-", SBML.Math[KINETICMATH1, SBML.MathIdent("c1")]) From 441f3a1e03fa762968aff887b322e9797286fbed Mon Sep 17 00:00:00 2001 From: = <=> Date: Sun, 8 Aug 2021 12:49:25 -0400 Subject: [PATCH 2/4] fix in reversible --- src/reactionsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 8cf6610..a77a22e 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -96,7 +96,7 @@ function mtk_reactions(model::SBML.Model) kl_rv = from_noncombinatoric(kl_rv, reagents[3], our) push!(rxs, ModelingToolkit.Reaction(kl_fw, reagents...; only_use_rate=our)) - reagents = getreagents(rstoich, pstoich, model) + reagents = getreagents(pstoich, rstoich, model) 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)) From 045a42b88d3a64aefcea0fb32d738eb3ef7acf89 Mon Sep 17 00:00:00 2001 From: = <=> Date: Sun, 8 Aug 2021 13:10:49 -0400 Subject: [PATCH 3/4] fix test stoich --- test/reactionsystem.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/reactionsystem.jl b/test/reactionsystem.jl index 9bb325e..409e97d 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), Dict(), (NaN, ""), (NaN, ""), NaN, +REACTION1 = SBML.Reaction(Dict(), Dict("s1" => 1), (NaN, ""), (NaN, ""), NaN, nothing, KINETICMATH1, false) -REACTION2 = SBML.Reaction(Dict(), Dict("s2" => 1), (NaN, ""), (NaN, ""), NaN, +REACTION2 = SBML.Reaction(Dict("s2" => 1), Dict(), (NaN, ""), (NaN, ""), NaN, nothing, KINETICMATH2, false) -REACTION3 = SBML.Reaction(Dict(), 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()) From 0eccdde23e61f32f0bd1ef820a49b28fb0be2c5a Mon Sep 17 00:00:00 2001 From: = <=> Date: Sun, 8 Aug 2021 18:59:20 -0400 Subject: [PATCH 4/4] one test left to fix --- src/reactionsystem.jl | 57 +++++++++++++++++++++++++++++------------- test/reactionsystem.jl | 2 +- 2 files changed, 40 insertions(+), 19 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index a77a22e..73582af 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -139,29 +139,50 @@ end """ Get reagents """ function getreagents(rstoichdict::Dict{String,<:Real}, pstoichdict::Dict{String,<:Real}, model::SBML.Model) - reactants, rstoich = stoichmap_to_varmap(rstoichdict) - products, pstoich = stoichmap_to_varmap(pstoichdict) - + reactants = Num[] + products = Num[] + rstoich = Float64[] + pstoich = Float64[] + 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 + # push!(reactants, create_var(k,Catalyst.DEFAULT_IV)) + # push!(rstoich, v) + end + if (length(reactants)==0) reactants = nothing; rstoich = nothing end if (length(products)==0) products = nothing; pstoich = nothing end (reactants, products, rstoich, pstoich) end -function stoichmap_to_varmap(stoich::Dict{String,<:Real}) - species = Num[] - stoichs = Float64[] - for (k,v) in stoich - if iszero(v) - @error("Stoichiometry of $k must be non-zero") - else - # if model.species[k].boundary_condition == true - push!(species, create_var(k,Catalyst.DEFAULT_IV)) - push!(stoichs, v) - # end - end - end - (species, stoichs) -end +# function stoichmap_to_varmap(stoich::Dict{String,<:Real}) +# species = Num[] +# stoichs = Float64[] +# for (k,v) in stoich +# if iszero(v) +# +# else +# # if model.species[k].boundary_condition == true +# push!(species, create_var(k,Catalyst.DEFAULT_IV)) +# push!(stoichs, v) +# # end +# end +# end +# (species, stoichs) +# end """ Infer forward and reverse components of bidirectional kineticLaw """ function getunidirectionalcomponents(bidirectional_math) diff --git a/test/reactionsystem.jl b/test/reactionsystem.jl index 409e97d..ee39aec 100644 --- a/test/reactionsystem.jl +++ b/test/reactionsystem.jl @@ -135,7 +135,7 @@ SBML.MathApply("*", SBML.Math[ SBML.MathIdent("k1"), SBML.MathIdent("constspec")]), SBML.MathIdent("k1")]) -constreac = SBML.Reaction(Dict(), 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.reactants, constreac.products, constmod)) @test isequal((nothing, nothing, nothing, nothing), SBMLToolkit.getreagents(constreac.reactants, constreac.products, constmod))