Merge branch 'main' into compathelper/new_version/2021-06-18-00-24-45…

…-340-744550298

Project.toml

@@ -7,9 +7,11 @@ version = "0.1.0"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
 Catalyst = "6"
+SBML = "0.5.3"
 julia = "1.6"
 
 [extras]

src/SBMLToolkit.jl

@@ -2,6 +2,7 @@ module SBMLToolkit
 
 using Catalyst, ModelingToolkit
 using SBML
+using SymbolicUtils
 
 include("reactionsystem.jl")
 

src/reactionsystem.jl

@@ -1,7 +1,6 @@
 """ ReactionSystem constructor """
 function ModelingToolkit.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires unique parameters (i.e. SBML must have been imported with localParameter promotion in libSBML)
     checksupport(model)
-    model = make_extensive(model)
     rxs = mtk_reactions(model)
     species = []
     for k in keys(model.species)
@@ -14,45 +13,20 @@ end
 """ ODESystem constructor """
 function ModelingToolkit.ODESystem(model::SBML.Model; kwargs...)
     rs = ReactionSystem(model; kwargs...)
-    model = make_extensive(model)  # PL: consider making `make_extensive!` to avoid duplicate calling in ReactionSystem and here
-    u0map = get_u0(model)
+    u0map = [create_var(k,Catalyst.DEFAULT_IV) => v for (k,v) in SBML.initial_amounts(model, convert_concentrations = true)]
     parammap = get_paramap(model)
     defaults = Dict(vcat(u0map, parammap))
     convert(ODESystem, rs, defaults=defaults)
 end
 
 """ Check if conversion to ReactionSystem is possible """
 function checksupport(model::SBML.Model)
-    for s in values(model.species)
-        if s.boundary_condition
-            @warn "Species $(s.name) has `boundaryCondition` or is `constant`. This will lead to wrong results when simulating the `ReactionSystem`."
-        end
-    end
-end
-
-""" Convert intensive to extensive expressions """
-function make_extensive(model)
-    model = to_initial_amounts(model)
-    model = to_extensive_math!(model)
-    model
-end
-
-""" Convert initial_concentration to initial_amount """
-function to_initial_amounts(model::SBML.Model)
-    model = deepcopy(model)
-    for specie in values(model.species)
-        if isnothing(specie.initial_amount)
-            compartment = model.compartments[specie.compartment]
-            if !isnothing(compartment.size)
-                specie.initial_amount = (specie.initial_concentration[1] * compartment.size, "")
-            else
-                @warn "Compartment $(compartment.name) has no `size`. Cannot calculate `initial_amount` of Species $(specie.name). Setting `initial_amount` to `initial_concentration`."
-                specie.initial_amount = (specie.initial_concentration[1], "")
-            end
-            specie.initial_concentration = nothing
-        end
-    end
-    model
+    # for s in values(model.species)
+    #     if s.boundary_condition
+    #         @warn "Species $(s.name) has `boundaryCondition` or is `constant`. This will lead to wrong results when simulating the `ReactionSystem`."
+    #     end
+    # end
+    return nothing
 end
 
 """ Convert intensive to extensive mathematical expression """
@@ -106,47 +80,64 @@ function mtk_reactions(model::SBML.Model)
     subsdict = _get_substitutions(model)
     rxs = []
     if length(model.reactions) == 0
-        throw(ErrorException("Model contains no reactions."))
+        throw(ErrorException("SBML.Model contains no reactions."))
     end
     for reaction in values(model.reactions)
-        reactants = Num[]
-        rstoich = Num[]
-        products = Num[]
-        pstoich = Num[]
-        for (k,v) in reaction.stoichiometry
-            if v < 0
-                push!(reactants, create_var(k,Catalyst.DEFAULT_IV))
-                push!(rstoich, -v)
-            elseif v > 0
-                push!(products, create_var(k,Catalyst.DEFAULT_IV))
-                push!(pstoich, v)
-            else
-                @error("Stoichiometry of $k must be non-zero")
-            end
-        end
-        if (length(reactants)==0) reactants = nothing; rstoich = nothing end
-        if (length(products)==0) products = nothing; pstoich = nothing end
-        symbolic_math = convert(Num, reaction.kinetic_math)
-
+        extensive_math = SBML.extensive_kinetic_math(model, reaction.kinetic_math)
+        symbolic_math = convert(Num, extensive_math)
         if reaction.reversible
             symbolic_math = getunidirectionalcomponents(symbolic_math)
-            kl = [substitute(x, subsdict) for x in symbolic_math]
-            push!(rxs, ModelingToolkit.Reaction(kl[1],reactants,products,rstoich,pstoich;only_use_rate=true))
-            push!(rxs, ModelingToolkit.Reaction(kl[2],products,reactants,pstoich,rstoich;only_use_rate=true))
+            kl_fw, kl_rv = [substitute(x, subsdict) for x in symbolic_math]
+            reagents = getreagents(reaction.stoichiometry, model)
+            push!(rxs, ModelingToolkit.Reaction(kl_fw, reagents...; only_use_rate=true))
+            reagents = getreagents(reaction.stoichiometry, model; rev=true)
+            push!(rxs, ModelingToolkit.Reaction(kl_rv, reagents...; only_use_rate=true))
         else
             kl = substitute(symbolic_math, subsdict)
-            push!(rxs, ModelingToolkit.Reaction(kl,reactants,products,rstoich,pstoich;only_use_rate=true))
+            reagents = getreagents(reaction.stoichiometry, model)
+            push!(rxs, ModelingToolkit.Reaction(kl, reagents...; only_use_rate=true))
         end
     end
     rxs
 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
+    if (length(reactants)==0) reactants = nothing; rstoich = nothing end
+    if (length(products)==0) products = nothing; pstoich = nothing end
+    (reactants, products, rstoich, pstoich)
+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`."
     bidirectional_math = Symbolics.tosymbol(bidirectional_math)
     bidirectional_math = simplify(bidirectional_math; expand=true)
-    if SymbolicUtils.operation(bidirectional_math) != +
+    if (bidirectional_math isa Real) || (SymbolicUtils.operation(bidirectional_math) != +)
         throw(ErrorException(err))
     end
     terms = SymbolicUtils.arguments(bidirectional_math)
@@ -165,15 +156,6 @@ function getunidirectionalcomponents(bidirectional_math)
     return (fw_terms[1], rv_terms[1])
 end
 
-""" Extract u0map from Model """
-function get_u0(model)
-    u0map = []
-    for (k,v) in model.species
-        push!(u0map,Pair(create_var(k,Catalyst.DEFAULT_IV), v.initial_amount[1]))
-    end
-    u0map
-end
-
 """ Extract paramap from Model """
 function get_paramap(model)
     paramap = Pair{Num, Float64}[]
@@ -188,6 +170,30 @@ function get_paramap(model)
     paramap
 end
 
+""" Get rate constant of mass action kineticLaws """
+function getmassaction(kl::Num, reactants::Union{Vector{Num},Nothing}, stoich::Union{Vector{<:Real},Nothing})
+    function check_args(x::SymbolicUtils.Symbolic{Real})
+        for arg in SymbolicUtils.arguments(x)
+            if isnan(check_args(arg))
+                return NaN
+            end
+        end
+        return 0
+    end
+    check_args(x::Term{Real, Nothing}) = NaN  # Variable leaf node
+    check_args(x::Sym{Real, Base.ImmutableDict{DataType, Any}}) = 0  # Parameter leaf node
+    check_args(x::Real) = 0  # Real leaf node
+    check_args(x) = throw(ErrorException("Cannot handle $(typeof(x)) types."))  # Unknow leaf node
+    if isnothing(reactants) && isnothing(stoich)
+        rate_const = kl
+    elseif isnothing(reactants) | isnothing(stoich)
+        throw(ErrorException("`reactants` and `stoich` are incosistent: `reactants` are $(reactants) and `stoich` is $(stoich)."))
+    else
+        rate_const = kl / *((.^(reactants, stoich))...)
+    end
+    isnan(check_args(rate_const.val)) ? NaN : rate_const
+end
+
 create_var(x, iv) = Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(Symbol(x)))(iv).val
 create_var(x) = Num(Variable(Symbol(x))).val
 function create_param(x)

test/reactionsystem.jl

@@ -10,10 +10,16 @@
     KINETICMATH1 = SBML.MathIdent("k1")
     KINETICMATH2 = SBML.MathApply("*", SBML.Math[
         SBML.MathIdent("k1"), SBML.MathIdent("s2")])
-    REACTION1 = SBML.Reaction(Dict("s1" => 1), nothing, nothing, nothing, nothing, KINETICMATH1, false)
-    REACTION2 = SBML.Reaction(Dict("s2" => -1), nothing, nothing, nothing, nothing, KINETICMATH2, false)
+    KINETICMATH3 = SBML.MathApply("-", SBML.Math[KINETICMATH2, KINETICMATH1])
+    REACTION1 = SBML.Reaction(Dict("s1" => 1), (NaN, ""), (NaN, ""), NaN,
+                              nothing, KINETICMATH1, false)
+    REACTION2 = SBML.Reaction(Dict("s2" => -1), (NaN, ""), (NaN, ""), NaN,
+                              nothing, KINETICMATH2, false)
+    REACTION3 = SBML.Reaction(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())
+    MODEL3 = SBML.Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("s2" => SPECIES2), Dict("r3" => REACTION3), Dict(), Dict())
 
     # Test ReactionSystem constructor
     rs = ReactionSystem(MODEL1)
@@ -34,16 +40,16 @@
     @named rs = ReactionSystem(MODEL1)
     isequal(nameof(rs), :rs)
 
-    rs = ReactionSystem(readSBML(joinpath("data","reactionsystem_05.xml")))  # Contains reversible reaction
+    rs = ReactionSystem(MODEL3)  # Contains reversible reaction
     @test isequal(Catalyst.get_eqs(rs),
                   ModelingToolkit.Reaction[
-                      ModelingToolkit.Reaction(c1*k1*s1*s2, [s1,s2], [s1s2],
-                          [1.,1.], [1.]; use_only_rate=true),
-                      ModelingToolkit.Reaction(c1*k1*s1s2, [s1s2], [s1,s2],
-                          [1.], [1.,1.]; use_only_rate=true)])
+                      ModelingToolkit.Reaction(0.5k1*s2, [s2], nothing,
+                          [1.], nothing; use_only_rate=true),
+                      ModelingToolkit.Reaction(k1, nothing, [s2],
+                          nothing, [1.]; use_only_rate=true)])
     @test isequal(Catalyst.get_iv(rs), t)
-    @test isequal(Catalyst.get_states(rs), [s1, s1s2, s2])
-    @test isequal(Catalyst.get_ps(rs), [k1,c1])
+    @test isequal(Catalyst.get_states(rs), [s2])
+    @test isequal(Catalyst.get_ps(rs), [k1, c1])
 
     @test_nowarn convert(ModelingToolkit.ODESystem, rs)
     @test_nowarn structural_simplify(convert(ModelingToolkit.ODESystem, rs))
@@ -63,6 +69,7 @@
     @test_nowarn structural_simplify(odesys)
 
     odesys = ODESystem(readSBML(sbmlfile))
+    m = readSBML(sbmlfile)
     trueeqs = Equation[Differential(t)(s1) ~ -0.25c1 * k1 * s1 * s2,
                        Differential(t)(s1s2) ~ 0.25c1 * k1 * s1 * s2,
                        Differential(t)(s2) ~ -0.25c1 * k1 * s1 * s2]
@@ -79,86 +86,76 @@
     @test_nowarn ODEProblem(odesys, [], [0., 1.], [])
 
     # Test ODEProblem
-    oprob = ODEProblem(ODESystem(MODEL1), [0., 1.])
+    oprob = ODEProblem(ODESystem(MODEL1), [], [0., 1.], [])
     sol = solve(oprob, Tsit5())
     @test isapprox(sol.u, [[1.], [2.]])
 
-    @test_nowarn ODEProblem(ODESystem(readSBML(sbmlfile)), [0., 1.])
+    @test_nowarn ODEProblem(ODESystem(readSBML(sbmlfile)), [], [0., 1.], [])
 
     # Test checksupport
-    @test_nowarn SBML.checksupport(MODEL1)
-    sbc = SBML.Species("sbc", "c1", true, nothing, nothing, (1., "substance"), nothing, true)   
-    mod = SBML.Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("sbc" => sbc), Dict("r1" => REACTION1), Dict(), Dict())
-    @test_logs (:warn, "Species sbc has `boundaryCondition` or is `constant`. This will lead to wrong results when simulating the `ReactionSystem`.") SBML.checksupport(mod)
-
-    # Test make_extensive
-    model = SBML.make_extensive(MODEL2)
-    @test isequal(model.species["s2"].initial_amount, (2., ""))
-    @test isequal(model.species["s2"].initial_concentration, nothing)
-
-    kineticmath2_true = SBML.MathApply("*", SBML.Math[
-        SBML.MathIdent("k1"),
-        SBML.MathApply("/", SBML.Math[
-            SBML.MathIdent("s2"),
-            SBML.MathVal(2.0)])
-        ])
-    @test isequal(repr(model.reactions["r2"].kinetic_math), repr(kineticmath2_true))
-    @test model.species["s2"].only_substance_units
-
-    # Test to_initial_amounts
-    model = SBML.to_initial_amounts(MODEL1)
-    @test isequal(model.species["s1"].initial_amount, (1., "substance"))
-    @test isequal(model.species["s1"].initial_concentration, nothing)
-    model = SBML.to_initial_amounts(MODEL2)
-    @test isequal(model.species["s2"].initial_amount, (2., ""))
-    @test isequal(model.species["s2"].initial_concentration, nothing)
-
-    # Test to_extensive_math!
-    model = SBML.to_extensive_math!(MODEL1)
-    @test isequal(model.reactions["r1"].kinetic_math, KINETICMATH1)
-    model = SBML.to_extensive_math!(MODEL2)
-    @test isequal(repr(model.reactions["r2"].kinetic_math), repr(kineticmath2_true))
-    #PL: Todo @Mirek/Anand: Why does this not work without the `repr()`?
-    # Do we need to write a isequal(x::SBML.Math, y::SBML.Math) method to get this work?
-    @test model.species["s2"].only_substance_units
+    # @test_nowarn SBMLToolkit.checksupport(MODEL1)
+    # sbc = SBML.Species("sbc", "c1", true, nothing, nothing, (1., "substance"), nothing, true)   
+    # mod = SBML.Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("sbc" => sbc), Dict("r1" => REACTION1), Dict(), Dict())
+    # @test_logs (:warn, "Species sbc has `boundaryCondition` or is `constant`. This will lead to wrong results when simulating the `ReactionSystem`.") SBMLToolkit.checksupport(mod)
 
     # Test _get_substitutions
     truesubs = Dict(Num(Variable(:c1)) => c1,
                 Num(Variable(:k1)) => k1,
                 Num(Variable(:s1)) => s1)
-    subs = SBML._get_substitutions(MODEL1)
+    subs = SBMLToolkit._get_substitutions(MODEL1)
     @test isequal(subs, truesubs)
 
     # Test mtk_reactions
-    reaction = SBML.mtk_reactions(MODEL1)[1]
+    reaction = SBMLToolkit.mtk_reactions(MODEL1)[1]
     truereaction = ModelingToolkit.Reaction(k1, nothing, [s1], nothing, [1]; only_use_rate=true)  # Todo: implement Sam's suggestion on mass action kinetics
     @test isequal(reaction, truereaction)
 
+    # Test getreagents
+    @test isequal((nothing, [s1], nothing, [1.]), SBMLToolkit.getreagents(REACTION1.stoichiometry, MODEL1))
+
+    constspec = SBML.Species("constspec", "c1", true, nothing, nothing, (1., "substance"), nothing, true)
+    ncs = SBMLToolkit.create_var("constspec",Catalyst.DEFAULT_IV)
+    kineticmath = SBML.MathApply("-", SBML.Math[
+        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)
+    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 getunidirectionalcomponents
     km = SBML.MathApply("-", SBML.Math[KINETICMATH1, SBML.MathIdent("c1")])
     sm = convert(Num, km)
-    kl = SBML.getunidirectionalcomponents(sm)
+    kl = SBMLToolkit.getunidirectionalcomponents(sm)
     @test isequal(kl, (k1, c1))
 
     km = SBML.MathApply("-", SBML.Math[KINETICMATH1, KINETICMATH2])
     sm = convert(Num, km)
-    fw, rv = SBML.getunidirectionalcomponents(sm)
-    rv = substitute(rv, Dict(SBML.create_var("s2") => SBML.create_var("s2",Catalyst.DEFAULT_IV)))
+    fw, rv = SBMLToolkit.getunidirectionalcomponents(sm)
+    rv = substitute(rv, Dict(SBMLToolkit.create_var("s2") => SBMLToolkit.create_var("s2",Catalyst.DEFAULT_IV)))
     @test isequal((fw, rv), (k1, k1*s2))
 
     km = SBML.MathIdent("s1s2")
     sm1 = convert(Num, km)
     sm2 = sm - sm1
-    @test_throws ErrorException SBML.getunidirectionalcomponents(sm2)
-
-    # Test get_u0
-    true_u0map = [s1 => 1.]
-    u0map = SBML.get_u0(MODEL1)
-    @test isequal(true_u0map, u0map)
+    @test_throws ErrorException SBMLToolkit.getunidirectionalcomponents(sm2)
 
     # Test get_paramap
     trueparamap = [k1 => 1., c1 => 2.]
-    paramap = SBML.get_paramap(MODEL1)
+    paramap = SBMLToolkit.get_paramap(MODEL1)
     @test isequal(paramap, trueparamap)
 
+    # Test getmassaction
+    @test isequal(SBMLToolkit.getmassaction(k1*s1, [s1], [1]), k1)  # Case hOSU=true
+    @test isequal(SBMLToolkit.getmassaction(k1*s1/c1, [s1], [1]), k1/c1)  # Case hOSU=false
+    @test isequal(SBMLToolkit.getmassaction(k1+c1, nothing, nothing), k1+c1)  # Case zero order kineticLaw
+    @test isnan(SBMLToolkit.getmassaction(k1*s1*s2/(c1+s2), [s1], [1]))  # Case Michaelis-Menten kinetics
+
+    @test isnan(SBMLToolkit.getmassaction(k1*s1*s2, [s1], [1]))  
+    @test isnan(SBMLToolkit.getmassaction(k1+c1, [s1], [1]))
+    @test_throws ErrorException SBMLToolkit.getmassaction(k1, nothing, [1])
+
 end