add support for reaction coupling

Project.toml

@@ -1,7 +1,7 @@
 name = "AbstractFBCModels"
 uuid = "5a4f3dfa-1789-40f8-8221-69268c29937c"
 authors = ["Authors of AbstractFBCModels.jl"]
-version = "0.2.2"
+version = "0.3.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

docs/src/canonical.jl

@@ -89,6 +89,11 @@ m.reactions["exchange2"] = Reaction(
     stoichiometry = Dict("m2" => -1.0),
     gene_association_dnf = [], # DNF encoding of a reaction that never has gene products available
 )
+m.couplings["total_exchange_limit"] = Coupling(
+    lower_bound = 0,
+    upper_bound = 10,
+    reaction_weights = Dict("exchange$i" => 1.0 for i = 1:2),
+)
 nothing #hide
 
 show_contains(x, y) = contains(sprint(show, MIME"text/plain"(), x), y) #src
@@ -99,7 +104,9 @@ show_contains(x, y) = contains(sprint(show, MIME"text/plain"(), x), y) #src
 @test show_contains(m.reactions["forward"], "\"g1\"") #src
 @test show_contains(m.metabolites["m1"], "\"inside\"") #src
 @test show_contains(m.genes["g1"], "name = nothing") #src
-@test show_contains(m.genes["g1"], "name = nothing") #src
+@test show_contains(m.genes["g2"], "name = nothing") #src
+@test show_contains(m.couplings["total_exchange_limit"], "name = nothing") #src
+@test show_contains(m.couplings["total_exchange_limit"], "upper_bound = 10") #src
 
 # We should immediately find the basic accessors working:
 A.stoichiometry(m)
@@ -108,6 +115,10 @@ A.stoichiometry(m)
 
 A.objective(m)
 
+#
+
+A.coupling(m)
+
 # We can check various side things, such as which reactions would and would not work given all gene products disappear:
 products_available = [
     A.reaction_gene_products_available(m, rid, _ -> false) for

src/accessors.jl

@@ -65,6 +65,30 @@ n_genes(a::AbstractFBCModel)::Int = unimplemented(typeof(a), :n_genes)
 """
 $(TYPEDSIGNATURES)
 
+Return identifiers of all coupling bounds contained in the model. Empty if
+none.
+
+Coupling bounds are typically not named in models, but should be.
+
+COMPATIBILITY NOTE: Couplings currently default to empty to prevent breakage.
+This behavior will change with next major version.
+"""
+couplings(a::AbstractFBCModel)::Vector{String} = String[]
+
+"""
+$(TYPEDSIGNATURES)
+
+The number of coupling bounds in the model (must be equal to the length of
+vector given by [`couplings`](@ref)).
+
+This may be more efficient than calling [`couplings`](@ref) and measuring the
+array.
+"""
+n_couplings(a::AbstractFBCModel)::Int = length(couplings(a))
+
+"""
+$(TYPEDSIGNATURES)
+
 The sparse stoichiometric matrix of a given model.
 
 This usually corresponds to all the equality constraints in the model. The
@@ -75,23 +99,44 @@ stoichiometry(a::AbstractFBCModel)::SparseMat = unimplemented(typeof(a), :stoich
 """
 $(TYPEDSIGNATURES)
 
-Get the lower and upper bounds of all reactions in a model.
+Sparse matrix that describes the coupling of a given model.
+
+This usually corresponds to all additional constraints in the model, such as
+the ones used for split-direction reactions and community modeling. The matrix
+must be of size [`n_couplings`](@ref) by [`n_reactions`](@ref).
+"""
+coupling(a::AbstractFBCModel)::SparseMat = spzeros(n_couplings(a), n_reactions(a))
+
+"""
+$(TYPEDSIGNATURES)
+
+Lower and upper bounds of all reactions in the model.
 """
 bounds(a::AbstractFBCModel)::Tuple{Vector{Float64},Vector{Float64}} =
     unimplemented(typeof(a), :bounds)
 
 """
 $(TYPEDSIGNATURES)
 
+Lower and upper bounds of all couplings in the model.
+"""
+coupling(a::AbstractFBCModel)::Tuple{Vector{Float64},Vector{Float64}} =
+    (fill(-Inf, n_couplings(a)), fill(Inf, n_couplings(a)))
+
+"""
+$(TYPEDSIGNATURES)
+
 Get the sparse balance vector of a model, which usually corresponds to the
 accumulation term associated with stoichiometric matrix.
+
+By default, the balance is assumed to be exactly zero.
 """
 balance(a::AbstractFBCModel)::SparseVec = spzeros(n_metabolites(a))
 
 """
 $(TYPEDSIGNATURES)
 
-Get the objective vector of the model.
+The objective vector of the model.
 """
 objective(a::AbstractFBCModel)::SparseVec = unimplemented(typeof(a), :objective)
 
@@ -245,3 +290,37 @@ $(TYPEDSIGNATURES)
 The name of the given gene in the model, if recorded.
 """
 gene_name(::AbstractFBCModel, gene_id::String)::Maybe{String} = nothing
+
+"""
+$(TYPEDSIGNATURES)
+
+The weights of reactions in the given coupling bound. Returns a dictionary that
+maps the reaction IDs to their weights.
+
+Using this function may be more efficient in cases than loading the whole
+[`coupling`](@ref).
+"""
+coupling_weights(a::AbstractFBCModel, coupling_id::String)::Dict{String,Float64} = Dict()
+
+"""
+$(TYPEDSIGNATURES)
+
+A dictionary of standardized names that may help to identify the corresponding
+coupling.
+"""
+coupling_annotations(::AbstractFBCModel, coupling_id::String)::Annotations = Dict()
+
+"""
+$(TYPEDSIGNATURES)
+
+Free-text notes organized in a dictionary by topics about the given coupling in
+the model.
+"""
+coupling_notes(::AbstractFBCModel, coupling_id::String)::Notes = Dict()
+
+"""
+$(TYPEDSIGNATURES)
+
+The name of the given coupling in the model, if recorded.
+"""
+coupling_name(::AbstractFBCModel, coupling_id::String)::Maybe{String} = nothing

src/canonical.jl

@@ -67,6 +67,26 @@ Base.show(io::Base.IO, ::MIME"text/plain", x::Gene) = A.pretty_print_kwdef(io, x
 """
 $(TYPEDEF)
 
+A canonical Julia representation of a row in a coupling matrix of the
+`AbstractFBCModels` interface.
+
+# Fields
+$(TYPEDFIELDS)
+"""
+Base.@kwdef mutable struct Coupling
+    name::A.Maybe{String} = nothing
+    reaction_weights::Dict{String,Float64} = Dict()
+    lower_bound::Float64 = -Inf
+    upper_bound::Float64 = Inf
+    annotations::A.Annotations = A.Annotations()
+    notes::A.Notes = A.Notes()
+end
+
+Base.show(io::Base.IO, ::MIME"text/plain", x::Coupling) = A.pretty_print_kwdef(io, x)
+
+"""
+$(TYPEDEF)
+
 A canonical Julia representation of a metabolic model that sotres exactly the
 data represented by `AbstractFBCModels` accessors.
 
@@ -85,25 +105,31 @@ Base.@kwdef struct Model <: A.AbstractFBCModel
     reactions::Dict{String,Reaction} = Dict()
     metabolites::Dict{String,Metabolite} = Dict()
     genes::Dict{String,Gene} = Dict()
+    couplings::Dict{String,CouplingBound} = Dict()
 end
 
 Base.show(io::Base.IO, ::MIME"text/plain", x::Model) = A.pretty_print_kwdef(io, x)
 
 A.reactions(m::Model) = sort(collect(keys(m.reactions)))
 A.metabolites(m::Model) = sort(collect(keys(m.metabolites)))
 A.genes(m::Model) = sort(collect(keys(m.genes)))
+A.couplings(m::Model) = sort(collect(keys(m.couplings)))
 A.n_reactions(m::Model) = length(m.reactions)
 A.n_metabolites(m::Model) = length(m.metabolites)
 A.n_genes(m::Model) = length(m.genes)
+A.n_couplings(m::Model) = length(m.couplings)
 A.reaction_name(m::Model, id::String) = m.reactions[id].name
 A.metabolite_name(m::Model, id::String) = m.metabolites[id].name
 A.gene_name(m::Model, id::String) = m.genes[id].name
+A.coupling_name(m::Model, id::String) = m.couplings[id].name
 A.reaction_annotations(m::Model, id::String) = m.reactions[id].annotations
 A.metabolite_annotations(m::Model, id::String) = m.metabolites[id].annotations
 A.gene_annotations(m::Model, id::String) = m.genes[id].annotations
+A.coupling_annotations(m::Model, id::String) = m.couplings[id].annotations
 A.reaction_notes(m::Model, id::String) = m.reactions[id].notes
 A.metabolite_notes(m::Model, id::String) = m.metabolites[id].notes
 A.gene_notes(m::Model, id::String) = m.genes[id].notes
+A.coupling_notes(m::Model, id::String) = m.couplings[id].notes
 
 function A.stoichiometry(m::Model)
     midxs = Dict(mid => idx for (idx, mid) in enumerate(A.metabolites(m)))
@@ -120,11 +146,31 @@ function A.stoichiometry(m::Model)
     sparse(I, J, V, A.n_metabolites(m), A.n_reactions(m))
 end
 
+function A.coupling(m::Model)
+    ridxs = Dict(rid => idx for (idx, rid) in enumerate(A.reactions(m)))
+    I = Int[]
+    J = Int[]
+    V = Float64[]
+    for (cidx, cid) in enumerate(A.couplings(m))
+        for (rid, v) in m.couplings[cid].reaction_weights
+            push!(I, cidx)
+            push!(J, ridxs[rid])
+            push!(V, v)
+        end
+    end
+    sparse(I, J, V, A.n_couplings(m), A.n_reactions(m))
+end
+
 A.bounds(m::Model) = (
     [m.reactions[rid].lower_bound for rid in A.reactions(m)],
     [m.reactions[rid].upper_bound for rid in A.reactions(m)],
 )
 
+A.coupling_bounds(m::Model) = (
+    [m.couplings[cid].lower_bound for cid in A.couplings(m)],
+    [m.couplings[cid].upper_bound for cid in A.couplings(m)],
+)
+
 A.balance(m::Model) =
     sparse(Float64[m.metabolites[mid].balance for mid in A.metabolites(m)])
 A.objective(m::Model) =
@@ -139,15 +185,15 @@ A.metabolite_formula(m::Model, id::String) = m.metabolites[id].formula
 A.metabolite_charge(m::Model, id::String) = m.metabolites[id].charge
 A.metabolite_compartment(m::Model, id::String) = m.metabolites[id].compartment
 
+A.coupling_weights(m::Model, id::String) = m.couplings[id].reaction_weights
+
 A.load(::Type{Model}, path::String) = S.deserialize(path)
 A.save(m::Model, path::String) = S.serialize(path, m)
 A.filename_extensions(::Type{Model}) = ["canonical-serialized-fbc"]
 
 function Base.convert(::Type{Model}, x::A.AbstractFBCModel)
     (lbs, ubs) = A.bounds(x)
-    os = A.objective(x)
-    bs = A.balance(x)
-    mets = A.metabolites(x)
+    (clbs, cubs) = A.coupling_bounds(x)
     Model(
         reactions = Dict(
             r => Reaction(
@@ -159,7 +205,7 @@ function Base.convert(::Type{Model}, x::A.AbstractFBCModel)
                 gene_association_dnf = A.reaction_gene_association_dnf(x, r),
                 annotations = A.reaction_annotations(x, r),
                 notes = A.reaction_notes(x, r),
-            ) for (r, o, lb, ub) in zip(A.reactions(x), os, lbs, ubs)
+            ) for (r, o, lb, ub) in zip(A.reactions(x), A.objective(x), lbs, ubs)
         ),
         metabolites = Dict(
             m => Metabolite(
@@ -170,7 +216,7 @@ function Base.convert(::Type{Model}, x::A.AbstractFBCModel)
                 compartment = A.metabolite_compartment(x, m),
                 annotations = A.metabolite_annotations(x, m),
                 notes = A.metabolite_notes(x, m),
-            ) for (m, b) in zip(mets, bs)
+            ) for (m, b) in zip(A.metabolites(x), A.balance(x))
         ),
         genes = Dict(
             g => Gene(
@@ -179,6 +225,16 @@ function Base.convert(::Type{Model}, x::A.AbstractFBCModel)
                 notes = A.gene_notes(x, g),
             ) for g in A.genes(x)
         ),
+        couplings = Dict(
+            cid => Coupling(
+                name = A.coupling_name(x, c),
+                lower_bound = lb,
+                upper_bound = ub,
+                reaction_weights = A.coupling_weights(x, c),
+                annotations = A.coupling_annotations(x, c),
+                notes = A.coupling_notes(x, c),
+            ) for (c, lb, ub) in zip(A.couplings(x), clbs, cubs)
+        ),
     )
 end
 

src/testing.jl

@@ -44,22 +44,34 @@ function run_fbcmodel_file_tests(
         @test m2 isa X
 
         S = stoichiometry(model)
+        C = coupling(model)
         rxns = reactions(model)
         mets = metabolites(model)
         gens = genes(model)
+        cpls = couplings(model)
 
         @test n_reactions(model) == length(rxns)
         @test n_metabolites(model) == length(mets)
         @test n_genes(model) == length(gens)
+        @test n_couplings(model) == length(cpls)
 
         # test sizing
         @test size(S) == (length(mets), length(rxns))
+        @test size(C) == (length(cpls), length(rxns))
         @test length(balance(model)) == size(S, 1)
+
         bs = bounds(model)
         @test bs isa Tuple{Vector{Float64},Vector{Float64}}
         lbs, ubs = bs
         @test length(lbs) == size(S, 2)
         @test length(ubs) == size(S, 2)
+
+        cbs = coupling_bounds(model)
+        @test cbs isa Tuple{Vector{Float64},Vector{Float64}}
+        clbs, cubs = cbs
+        @test length(clbs) == size(C, 1)
+        @test length(cubs) == size(C, 1)
+
         obj = objective(model)
         @test length(obj) == size(S, 2)
 
@@ -70,6 +82,18 @@ function run_fbcmodel_file_tests(
                     @atest met in ms "metabolite `$met' in reaction_stoichiometry() of `$rid' is in metabolites()"
                     @atest S[mi[met], ridx] == stoi "reaction_stoichiometry() of reaction `$rid' matches the column in stoichiometry() matrix"
                 end
+                # TODO also test the other direction
+            end
+        end
+
+        let ss = Set(rxns), ri = Dict(rxns .=> 1:length(rxns))
+            for (cidx, cid) in enumerate(cpls)
+                for (rxn, w) in coupling_weights(model, cid)
+                    # test if coupling weights are the same as with the matrix
+                    @atest rxn in rs "reaction `$rxn' in coupling_weights() of `$cid' is in reactions()"
+                    @atest C[cidx, ri[rxn]] == w "coupling_weights() of coupling `$cid' matches the row in coupling() matrix"
+                end
+                # TODO also test the other direction
             end
         end
 
@@ -115,6 +139,7 @@ function run_fbcmodel_file_tests(
             @test issetequal(rxns, reactions(m))
             @test issetequal(mets, metabolites(m))
             @test issetequal(gens, genes(m))
+            @test issetequal(cpls, couplings(m))
 
             @test Dict(rxns .=> collect(obj)) ==
                   Dict(reactions(m) .=> collect(objective(m)))
@@ -148,8 +173,12 @@ function run_fbcmodel_type_tests(::Type{X}) where {X<:AbstractFBCModel}
         rt(n_metabolites, Int, X)
         rt(genes, Vector{String}, X)
         rt(n_genes, Int, X)
+        rt(couplings, Vector{String}, X)
+        rt(n_couplings, Int, X)
         rt(stoichiometry, SparseMat, X)
+        rt(coupling, SparseMat, X)
         rt(bounds, Tuple{Vector{Float64},Vector{Float64}}, X)
+        rt(coupling_bounds, Tuple{Vector{Float64},Vector{Float64}}, X)
         rt(objective, SparseVec, X)
         rt(balance, SparseVec, X)
 
@@ -170,5 +199,10 @@ function run_fbcmodel_type_tests(::Type{X}) where {X<:AbstractFBCModel}
         rt(gene_name, Maybe{String}, X, String)
         rt(gene_annotations, Annotations, X, String)
         rt(gene_notes, Notes, X, String)
+
+        rt(coupling_weights, Dict{String,Float64}, X, String)
+        rt(coupling_name, Maybe{String}, X, String)
+        rt(coupling_annotations, Annotations, X, String)
+        rt(coupling_notes, Notes, X, String)
     end
 end