draft Model to MTK.___System conversion.

Project.toml

@@ -4,11 +4,14 @@ authors = ["Mirek Kratochvil <[email protected]>", "LCSB R3 team <lcsb-
 version = "0.4.0"
 
 [deps]
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SBML_jll = "bb12108a-f4ef-5f88-8ef3-0b33ff7017f1"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]

src/SBML.jl

@@ -8,10 +8,11 @@ include("version.jl")
 include("readsbml.jl")
 include("math.jl")
 include("utils.jl")
+include("reactionsystem.jl")
 
 sbml = (sym::Symbol) -> dlsym(SBML_jll.libsbml_handle, sym)
 
 export SBMLVersion,
-    readSBML, Model, UnitPart, Species, Reaction, getS, getLBs, getUBs, getOCs
+    readSBML, Model, UnitPart, Compartment, Species, Reaction, getS, getLBs, getUBs, getOCs
 
 end # module

src/reactionsystem.jl

@@ -0,0 +1,152 @@
+using Catalyst, ModelingToolkit, Symbolics
+
+# module SBML
+# struct Model end
+# struct Reaction end
+# struct Compartment end
+# struct Species end
+# end
+
+""" ReactionSystem constructor """
+function ModelingToolkit.ReactionSystem(model::Model; kwargs...)  # Todo: requires unique parameters (i.e. SBML must have been imported with localParameter promotion in libSBML)
+    model = make_extensive(model)
+    model = expand_reversible(model)
+    rxs = mtk_reaction.(model.reactions)
+    t = DEFAULT_IV
+    species = [create_var(k) for k in keys(model.species)]
+    params = vcat([create_var(k) for k in keys(model.parameters)], [create_par(k) for k in keys(model.compartments)])
+    ReactionSystem(rxs,t,species,params; kwargs...)
+end
+
+""" ReactionSystem constructor """
+function ModelingToolkit.ReactionSystem(sbmlfile::String; kwargs...)
+    model = readSBML(sbmlfile)
+    ReactionSystem(model; kwargs...)
+end
+
+""" ODESystem constructor """
+function ModelingToolkit.ODESystem(model::Model; kwargs...)
+    rs = ReactionSystem(model, kwargs...)
+    u0map = get_u0(model)
+    parammap = get_paramap(model)
+    defaults = vcat(u0map, parammap)
+    convert(ODESystem, rs, defaults=defaults)
+end
+
+""" ODESystem constructor """
+function ModelingToolkit.ODESystem(sbmlfile::String; kwargs...)
+    model = readSBML(sbmlfile)
+    ODESystem(model; kwargs...)
+end
+
+""" ODEProblem constructor """
+function ModelingToolkit.ODEProblem(model::Model,tspan;kwargs...)  # PL: Todo: add u0 and parameters argument
+    odesys = ODESystem(model,kwargs...)
+    ODEProblem(odesys, [], tspan)
+end
+
+""" ODEProblem constructor """
+function ModelingToolkit.ODEProblem(sbmlfile::String,tspan;kwargs...)  # PL: Todo: add u0 and parameters argument
+    odesys = ODESystem(sbmlfile;kwargs...)
+    ODEProblem(odesys, [], tspan)
+end
+
+""" Convert intensive to extensive expressions """
+function make_extensive(model)
+    model = to_initial_amounts(model)
+    # model = to_extensive_math(model)
+    model  # Todo: For spevies with `hOSU=false` multiply all occurences in mathematical expressions by compartment size.
+           # Also convert species initialConcentrations to initialAmounts
+end
+
+""" Convert initial_concentration to initial_amount """
+function to_initial_amounts(model::Model)
+    model = deepcopy(model)
+    for specie in values(model.species)
+        if isequal(specie.initial_amount, nothing)
+            compartment = model.compartments[specie.compartment]
+            specie.initial_amount = (specie.initial_concentration[1] * compartment.size, "")
+            specie.initial_concentration = nothing
+        end
+    end
+    model
+end
+
+""" Convert intensive to extensive mathematical expression """
+function to_extensive_math(model::Model)
+    model = deepcopy(model)
+    for reaction in model.reactions
+        km = reaction.kinetic_math
+        reaction.km = 1.  # PL: Todo: @Anand can you multiply species with `hOSU=true` with their compartment volume?
+    end
+    reaction
+end
+
+""" Expand reversible reactions to two reactions """
+function expand_reversible(model)
+    model  # Todo: convert all Reactions that are `reversible=true` to a forward and reverse reaction with `reversible=false`.
+end
+
+""" Convert SBML.Reaction to MTK.Reaction """
+function mtk_reaction(reaction::SBML.Reaction)
+    reactants = []
+    rstoich = []
+    products = []
+    pstoich = []
+    for (k,v) in reaction.stoichiometry
+        if v < 0
+            push!(reactants, create_var(k))
+            push!(rstoich, -v)
+        elseif v > 0
+            push!(products, create_var(k))
+            push!(pstoich, -v)
+        else
+            @error("Stoichiometry of $k must be non-zero")
+        end
+    end
+    subsdict = _get_substitutions(model)
+    # PL: Todo: @Anand: can you convert kinetic_math to Symbolic expression. Perhaps it would actually better if kinetic Math would be a Symbolics.jl expression rather than of type `Math`? But Mirek wants `Math`, I think.
+    kl = substitute(reaction.kinetic_math, subsdict)  # PL: Todo: might need conversion of kinetic_math to Symbolic MTK expression
+    ModelingToolkit.Reaction(reaction.kinetic_math,reactants,prodcts,rstoich,pstoich;only_use_rate=true)
+end
+
+""" Get dictonary to change types in kineticLaw """
+function _get_substitution(model)
+    subsdict = Dict()
+    for k in keys(model.species)
+        push!(substict, Pair(Num(Variable(Symbol(k))),create_var(k)))
+    end
+    for k in keys(model.parameters)
+        push!(subsdict, Pair(Num(Variable(Symbol(k))),create_par(k)))
+    end
+    for k in keys(model.compartments)
+        push!(subsdict, Pair(Num(Variable(Symbol(k))),create_par(k)))
+    end
+    subsdict
+end
+
+""" Extract u0map from Model """
+function get_u0(model)
+    u0map = []
+    for (k,v) in model.species
+        push!(Pair(create_var(k),vinitial_amount))
+    end
+    u0map
+end
+
+""" Extract paramap from Model """
+function get_paramap(model)
+    paramap = []
+    for (k,v) in model.parameters
+        push!(Pair(create_par(k),v))
+    end
+    for (k,v) in model.compartments
+        push!(Pair(create_par(k),v.size))
+    end
+    paramap
+end
+
+create_var(x) = Num(Variable(Symbol(x)))
+# # create_var(x, iv) = Num(Sym{FnType{Tuple{Real}}}(Symbol(x))(Variable(Symbol(iv)))).val
+# # create_var(x, iv) = Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(Symbol(x)))(Variable(Symbol(iv)))
+create_param(x) = Num(Sym{ModelingToolkit.Parameter{Real}}(Symbol(x)))

src/readsbml.jl

@@ -228,6 +228,14 @@ function extractModel(mdl::VPtr)::Model
             )
         end
 
+        ic = nothing
+        if ccall(sbml(:Species_isSetInitialConcentration), Cint, (VPtr,), sp) != 0
+            ic = (
+                ccall(sbml(:Species_getInitialConcentration), Cdouble, (VPtr,), sp),
+                get_string(sp, :Species_getSubstanceUnits),
+            )
+        end
+
         species[get_string(sp, :Species_getId)] = Species(
             get_optional_string(sp, :Species_getName),
             get_string(sp, :Species_getCompartment),
@@ -239,6 +247,7 @@ function extractModel(mdl::VPtr)::Model
             formula,
             charge,
             ia,
+            ic,
             get_optional_bool(
                 sp,
                 :Species_isSetHasOnlySubstanceUnits,
@@ -273,6 +282,9 @@ function extractModel(mdl::VPtr)::Model
     reactions = Dict{String,Reaction}()
     for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
         re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
+        if ccall(sbml(:Reaction_getReversible), Cint, (VPtr,), re)
+            throw(AssertionError("Reaction $(get_string(re, :Reaction_getId)) is reversible, but currently only irreversible reactions are supported."))
+        end     
         lb = (-Inf, "") # (bound value, unit id)
         ub = (Inf, "")
         oc = 0.0

src/structs.jl

@@ -110,11 +110,11 @@ consumption/production rates (accessible in field `stoichiometry`), lower/upper
 bounds (in tuples `lb` and `ub`, with unit names), and objective coefficient
 (`oc`). Also may contains `notes` and `annotation`.
 """
-struct Reaction
-    stoichiometry::Dict{String,Float64}
-    lb::Tuple{Float64,String}
-    ub::Tuple{Float64,String}
-    oc::Float64
+mutable struct Reaction
+    stoichiometry::Dict{String,Float64}  # PL: Why do we allow/require non-integer stoichiometries?
+    lb::Maybe{Tuple{Float64,String}}
+    ub::Maybe{Tuple{Float64,String}}
+    oc::Maybe{Float64}
     gene_product_association::Maybe{GeneProductAssociation}
     kinetic_math::Maybe{Math}
     notes::Maybe{String}
@@ -126,13 +126,14 @@ end
 Species metadata -- contains a human-readable `name`, a `compartment`
 identifier, `formula`, `charge`, and additional `notes` and `annotation`.
 """
-struct Species
+mutable struct Species
     name::Maybe{String}
     compartment::String
     boundary_condition::Maybe{Bool}
     formula::Maybe{String}
     charge::Maybe{Int}
     initial_amount::Maybe{Tuple{Float64,String}}
+    initial_concentration::Maybe{Tuple{Float64,String}}
     only_substance_units::Maybe{Bool}
     notes::Maybe{String}
     annotation::Maybe{String}
@@ -175,8 +176,8 @@ struct Model
     compartments::Dict{String,Compartment}
     species::Dict{String,Species}
     reactions::Dict{String,Reaction}
-    gene_products::Dict{String,GeneProduct}
-    function_definitions::Dict{String,FunctionDefinition}
+    gene_products::Maybe{Dict{String,GeneProduct}}
+    function_definitions::Maybe{Dict{String,FunctionDefinition}}
     notes::Maybe{String}
     annotation::Maybe{String}
     Model(p, u, c, s, r, g, f, n = nothing, a = nothing) = new(p, u, c, s, r, g, f, n, a)

test/reactionsystem.jl

@@ -0,0 +1,67 @@
+
+# sbmlfile = "reactionsystem_01.xml"
+
+COMP1 = Compartment("c1", true, 3, 2., "nl") 
+SPECIES1 = Species("s1", "c1", false, nothing, nothing, (1., "substance"), nothing, true)  # Todo: Maybe not support units in initial_concentration?
+SPECIES2 = Species("s2", "c1", false, nothing, nothing, nothing, (1., "substance/nl"), false)
+KINETICMATH1 = SBML.MathVal(1.)  # PL: @anand or @mirek help needed. Can someone create `k1 * SPECIES1` as a `Math` type here?
+REACTION1 = Reaction(Dict("r1" => 1), nothing, nothing, nothing, nothing, KINETICMATH1)
+REACTION2 = Reaction(Dict("r2" => 1), nothing, nothing, nothing, nothing, KINETICMATH1)
+MODEL1 = Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("s1" => SPECIES1), Dict("r1" => REACTION1), nothing, nothing)  # PL: For instance in the compartments dict, we may want to enforce that key and compartment.name are identical
+MODEL2 = Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("s2" => SPECIES2), Dict("r2" => REACTION2), nothing, nothing)
+
+@testset "Model to MTK conversions" begin
+
+    # Test to_initial_amounts
+    model = SBML.to_initial_amounts(MODEL2)
+    @test isequal(model.species["s2"].initial_amount, (2., ""))
+    @test isequal(model.species["s2"].initial_concentration, nothing)
+    model = SBML.to_initial_amounts(MODEL1)
+    @test isequal(model.species["s1"].initial_amount, (1., "substance"))
+    @test isequal(model.species["s1"].initial_concentration, nothing)
+
+
+
+    # @test length(mdl.compartments) == 2
+
+    # mets, rxns, S = getS(mdl)
+
+    # @test typeof(S) <: AbstractMatrix{Float64}
+    # @test typeof(getS(mdl; zeros = spzeros)[3]) <: SparseMatrixCSC{Float64}
+    # @test typeof(getS(mdl; zeros = zeros)[3]) == Matrix{Float64}
+
+    # @test length(mets) == 77
+    # @test length(rxns) == 77
+    # @test size(S) == (length(mets), length(rxns))
+
+    # # totally arbitrary value tests
+    # @test isapprox(sum(S), 42.1479)
+    # @test mets[10:12] == ["M_akg_e", "M_fum_c", "M_pyr_c"]
+    # @test rxns[10:12] == ["R_H2Ot", "R_PGL", "R_EX_glc_e_"]
+
+    # lbs = getLBs(mdl)
+    # ubs = getUBs(mdl)
+    # ocs = getOCs(mdl)
+
+    # @test length(ocs) == length(mets)
+    # @test ocs[40] == 1.0
+    # deleteat!(ocs, 40)
+    # @test all(ocs .== 0.0)
+
+    # @test length(getLBs(mdl)) == length(rxns)
+    # @test length(getUBs(mdl)) == length(rxns)
+
+    # getunit = (val, unit)::Tuple -> unit
+    # @test all([broadcast(getunit, lbs) broadcast(getunit, ubs)] .== "mmol_per_gDW_per_hr")
+
+    # getval = (val, unit)::Tuple -> val
+    # lvals = broadcast(getval, lbs)
+    # uvals = broadcast(getval, ubs)
+    # @test isapprox(lvals[27], uvals[27])
+    # @test isapprox(lvals[27], 7.6)
+    # @test isapprox(lvals[12], -10)
+
+    # @test count(isapprox.(lvals, -999999)) == 40
+    # @test count(isapprox.(lvals, 0)) == 35
+    # @test count(isapprox.(uvals, 999999)) == 76
+end

test/runtests.jl

@@ -4,7 +4,8 @@ using SBML
 import Pkg
 
 @testset "SBML test suite" begin
-    include("version.jl")
-    include("ecoli_flux.jl")
-    include("loadmodels.jl")
+    # include("version.jl")
+    # include("ecoli_flux.jl")
+    # include("loadmodels.jl")
+    include("reactionsystem.jl")
 end