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

separate out reactants and products for new SBML version #28

Merged
merged 1 commit into from
Aug 9, 2021
Merged

separate out reactants and products for new SBML version #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
6 changes: 3 additions & 3 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"
Expand All @@ -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"

Expand Down
52 changes: 30 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(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)
Expand Down Expand Up @@ -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)
Expand Down
18 changes: 8 additions & 10 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,
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,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)

Expand All @@ -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)
Expand All @@ -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")])
Expand Down