add support for coupled-reaction bounds

.github/workflows/ci.yml

@@ -37,6 +37,7 @@ jobs:
       - uses: julia-actions/julia-runtest@latest
         continue-on-error: ${{ matrix.version == 'nightly' }}
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v4
         with:
-          file: lcov.info
+          files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}

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

@@ -64,7 +64,7 @@ A.run_fbcmodel_type_tests(Model);
 # For testing the values, you need to provide an existing file that contains
 # the model. Let's create some contents first:
 
-import AbstractFBCModels.CanonicalModel: Reaction, Metabolite, Gene
+import AbstractFBCModels.CanonicalModel: Reaction, Metabolite, Gene, Coupling
 
 m = Model()
 m.metabolites["m1"] = Metabolite(compartment = "inside")
@@ -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

docs/src/utilities.jl

@@ -12,7 +12,7 @@ import AbstractFBCModels as A
 
 A.accessors()
 
-@test length(A.accessors()) == 27 #src
+@test length(A.accessors()) == 35 #src
 
 # You do not need to overload all of them (e.g., if you model does not have any
 # genes you can completely omit all gene-related functions). The main required

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_bounds(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,Coupling} = 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(
+            c => 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