Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PR to sync with latest breaking SBML change to Reaction #26

Closed
wants to merge 4 commits into from
Closed

PR to sync with latest breaking SBML change to Reaction #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e664b13
PR to sync with latest breaking SBML change to Reaction
Aug 6, 2021
441f3a1
fix in reversible
Aug 8, 2021
045a42b
fix test stoich
Aug 8, 2021
0eccdde
one test left to fix
Aug 8, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 42 additions & 22 deletions src/reactionsystem.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(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))
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)
Expand Down Expand Up @@ -134,36 +137,53 @@ 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)

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)
Comment on lines -154 to -156
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@paulflang okay yeah I think i see. what is the reason we don't push to reactants if bc is true here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i guess because that would double things up right

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if bc is true we need to make sure that the stoichiometry of the species in question is the same on the right and left side of the reaction equation. however, we never want to modify the left hand side. this is because the left hand side can be used to determine if a reaction is mass action kinetics (i.e. A + B -> C gives mass action kinetic kAB, A -> B gives mass action kinetic k*A, etc.). This leaves us with the option to modify the product side. But this is already taken care of in line 149. In line 154 we only need to take care of the default case of bc = false, which is that reaction products are actually produced.

end
else
@error("Stoichiometry of $k must be non-zero")
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)
#
# 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`."
Expand Down
17 changes: 8 additions & 9 deletions test/reactionsystem.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
anandijain marked this conversation as resolved.
Show resolved Hide resolved
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())
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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("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))

# Test getunidirectionalcomponents
km = SBML.MathApply("-", SBML.Math[KINETICMATH1, SBML.MathIdent("c1")])
Expand Down