add interface functionality directly to the model builder

src/COBREXA.jl

@@ -44,6 +44,7 @@ import LinearAlgebra
 import SparseArrays
 
 include("types.jl")
+include("config.jl")
 
 # core functionality
 include("io.jl")
@@ -57,6 +58,7 @@ include("builders/knockouts.jl")
 include("builders/loopless.jl")
 include("builders/mmdf.jl")
 include("builders/objectives.jl")
+include("builders/scale.jl")
 include("builders/unsigned.jl")
 
 # simplified front-ends for the above

src/builders/communities.jl

@@ -17,59 +17,20 @@
 """
 $(TYPEDSIGNATURES)
 
-Create an "interface" block for constraints created by
-[`flux_balance_constraints`](@ref) (or any compatible constraint tree).
-"""
-fbc_boundary_constraints(
-    constraints::C.ConstraintTree;
-    ignore = _ -> false,
-    bound = _ -> nothing,
-) = network_boundary_constraints(c.fluxes.c.flux_stoichiometry; ignore, bound)
-
-"""
-$(TYPEDSIGNATURES)
-
-Create an "interface" block for a constrained reaction network described by
-variables in `fluxes` and flux balance in `balances`. Boundary reactions are
-assumed to be the ones that either "only create" or "only remove" the balanced
-components (metabolites), i.e, each dot product of the reaction description
-with all balance components is either strictly non-negative or non-positive.
-"""
-function network_boundary_constraints(
-    reactions,
-    balances;
-    ignore = _ -> false,
-    bound = _ -> nothing,
-)
-    #TODO
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Linearly scale all bounds in a constraint tree by the `factor`.
-"""
-function scale_bounds(tree::C.ConstraintTree, factor)
-    C.map(tree) do c
-        isnothing(c.bound) ? c : C.Constraint(value = c.value, bound = factor * c.bound)
-    end
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Join multiple modules into a bigger module.
+Join multiple constraint tree modules with interfaces into a bigger module with
+an interface.
 
 Modules are like usual constraint trees, but contain an explicitly declared
-`interface` part, marked properly a tuple (the parameters should form a
-dictionary constructor that would generally look such as `:module_name =>
-(module, module.interface)`; the second tuple member may also be specified just
-by name as e.g. `:interface`, or omitted while relying on `default_interface`).
+`interface` part, marked properly in arguments using e.g. a tuple (the
+parameters should form a dictionary constructor that would generally look such
+as `:module_name => (module, module.interface)`; the second tuple member may
+also be specified just by name as e.g. `:interface`, or omitted while relying
+on `default_interface`).
 
-Interface parts get merged and constrained to create a new interface; networks
-are copied intact.
+Compatible modules with interfaces may be created e.g. by
+[`flux_balance_constraints`](@ref).
 """
-function join_modules(
+function join_module_constraints(
     ps::Pair...;
     default_in_interface = :interface,
     out_interface = :interface,
@@ -85,23 +46,27 @@ function join_modules(
     prep(id::Symbol, (mod, interface)::Tuple{C.ConstraintTree,C.ConstraintTree}) =
         (id^(:network^mod * :interface^interface))
 
-    # collect everything into one huge network (while also renumbering the interfaces)
+    # first, collect everything into one huge network
+    # (while also renumbering the interfaces)
     modules = sum(prep.(ps))
 
-    # extract a union of all non-ignored interface keys
-    interface_keys =
-        filter(!ignore, collect(union((keys(m.interface) for (_, m) in modules)...)))
+    # fold a union of all non-ignored interface keys
+    interface_sum = foldl(modules, init = C.ConstraintTree()) do accs, (_, ms)
+        C.imerge(accs, ms) do path, acc, m
+            ignore(path) ? missing :
+            ismissing(acc) ? C.Constraint(value = m.value) :
+            C.Constraint(value = acc.value + m.value)
+        end
+    end
 
-    # remove the interface placeholders and add variables for the new interface
+    # extract the plain networks and add variables for the new interfaces
     constraints =
         ConstraintTree(id => (m.network) for (id, m) in modules) +
-        out_interface^C.variables(keys = interface_keys, bounds = bound.(interface_keys))
+        out_interface^C.variables_ifor(bound, interface_sum)
 
     # join everything with the interrace balance and return
-    constraints * C.map(constraints.out_interface) do ic
-        C.Constraint(
-            value = sum(c.value for mod in modules for (k, c) in mod.interface) - ic.value,
-        )
+    constraints * C.zip(interface_sum, constraints.out_interface) do sum, out
+        C.Constraint(value = sum.value - out.value)
     end
 end
 

src/builders/fbc.jl

@@ -24,16 +24,36 @@ The constructed tree contains subtrees `fluxes` (with the reaction-defining
 "variables") and `flux_stoichiometry` (with the metabolite-balance-defining
 constraints), and a single constraint `objective` thad describes the objective
 function of the model.
+
+Optionally if `interface` is specified, an "interface block" will be created
+within the constraint tree for later use as a "module" in creating bigger
+models (such as communities) using [`join_module_constraints`](@ref). The
+possible parameter values include:
+- `nothing` -- default, no interface is created
+- `:sbo` -- the interface gets created from model's SBO annotations)
+- `:identifier_prefixes` -- the interface is guesstimated from commonly
+  occurring adhoc reaction ID prefixes used in contemporary models
+- `:boundary` -- the interface is created from all reactions that either only
+  consume or only produce metabolites
+
+Output interface name can be set via `interface_name`.
+
+See [`Configuration`](@ref) for fine-tuning the default interface creation.
 """
-function flux_balance_constraints(model::A.AbstractFBCModel)
-    rxns = Symbol.(A.reactions(model))
+function flux_balance_constraints(
+    model::A.AbstractFBCModel;
+    interface::Maybe{Symbol} = nothing,
+    interface_name = :interface,
+)
+    rxn_strings = A.reactions(model)
+    rxns = Symbol.(rxn_strings)
     mets = Symbol.(A.metabolites(model))
     lbs, ubs = A.bounds(model)
     stoi = A.stoichiometry(model)
     bal = A.balance(model)
     obj = A.objective(model)
 
-    return C.ConstraintTree(
+    constraints = C.ConstraintTree(
         :fluxes^C.variables(keys = reactions, bounds = zip(lbs, ubs)) *
         :flux_stoichiometry^C.ConstraintTree(
             met => C.Constraint(
@@ -43,6 +63,33 @@ function flux_balance_constraints(model::A.AbstractFBCModel)
         ) *
         :objective^C.Constraint(C.LinearValue(SparseArrays.sparse(obj))),
     )
+
+    add_interface(sym, flt) =
+        any(flt) && (
+            constraints *=
+                interface_name^sym^C.ConstraintTree(
+                    r => constraints.fluxes[r] for r in rxns[flt]
+                )
+        )
+    if interface == :sbo
+        sbod(sbos, rid) = any(in(sbos), get(A.reaction_annotation(model, rid), "sbo", []))
+        add_interface(:exchanges, sbod.(Ref(configuration.exchange_sbos), rxn_strings))
+        add_interface(:biomass, sbod.(Ref(configuration.biomass_sbos), rxn_strings))
+        add_interface(
+            :atp_maintenance,
+            sbod.(Ref(configuration.atp_maintenance_sbos), rxn_strings),
+        )
+        add_interface(:demand, sbod.(Ref(configuration.demand_sbos), rxn_strings))
+    elseif interface == :identifier_prefixes
+        prefixed(ps, s) = any(p -> hasprefix(p, s), ps)
+        add_interface(:exchanges, prefixed.(Ref(configuration.exchange_id_prefixes), s))
+        add_interface(:biomass, prefixed.(Ref(configuration.biomass_id_prefixes), s))
+        add_interface(:atp_maintenance, in.(s, Ref(configuration.atp_maintenance_ids)))
+    elseif interface == :boundary
+        add_interface(:boundary, ((all(s .<= 0) || all(s .>= 0)) for s in eachcol(stoi)))
+    end
+
+    return constraints
 end
 
 export flux_balance_constraints

src/builders/scale.jl

@@ -0,0 +1,41 @@
+
+# Copyright (c) 2021-2024, University of Luxembourg
+# Copyright (c) 2021-2024, Heinrich-Heine University Duesseldorf
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""
+$(TYPEDSIGNATURES)
+
+Linearly scale all bounds in a constraint tree by the `factor`. This actually
+changes the model (and may not work in surprising/improper ways with some
+constraint systems, esp.  the MILP and QP ones).
+
+See also [`scale_constraints`](@ref).
+"""
+scale_bounds(tree::C.ConstraintTree, factor) =
+    C.map(tree) do c
+        isnothing(c.bound) ? c : C.Constraint(value = c.value, bound = factor * c.bound)
+    end
+
+"""
+$(TYPEDSIGNATURES)
+
+Linearly scale all constraints in a constraint tree by the `factor`.
+
+See also [`scale_bounds`](@ref).
+"""
+scale_constraints(tree::C.ConstraintTree, factor) =
+    C.map(tree) do c
+        c * factor
+    end

src/config.jl

@@ -0,0 +1,80 @@
+
+# Copyright (c) 2021-2024, University of Luxembourg
+# Copyright (c) 2021-2024, Heinrich-Heine University Duesseldorf
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""
+$(TYPEDEF)
+
+Global configuration options for various COBREXA functions, mainly for various
+non-interesting function parameters that are too inconvenient to be passed
+around manually.
+
+Changing the configuration values at runtime is possible via the global
+[`configuration`](@ref) variable.
+"""
+Base.@kwdef mutable struct Configuration
+
+    """
+    Prefixes that [`flux_balance_constraints`](@ref) uses for guessing
+    which reactions are exchanges.
+    """
+    exchange_id_prefixes::Vector{String} = ["EX_", "R_EX_"]
+
+    """
+    Prefixes that [`flux_balance_constraints`](@ref) uses for guessing
+    which reactions are biomass reactions.
+    """
+    biomass_id_prefixes::Vector{String} = ["BIOMASS_", "R_BIOMASS_"]
+
+    """
+    Reaction identifiers that [`flux_balance_constraints`](@ref) considers to
+    be ATP maintenance reactions.
+    """
+    atp_maintenance_ids::Vector{String} = ["ATPM", "R_ATPM"]
+
+    """
+    SBO numbers that label exchange reactions for
+    [`flux_balance_constraints`](@ref).
+    """
+    exchange_sbos::Vector{Int} = ["SBO:0000627"]
+
+    """
+    SBO numbers that label biomass production reactions for
+    [`flux_balance_constraints`](@ref).
+    """
+    biomass_sbos::Vector{Int} = ["SBO:0000629"]
+
+    """
+    SBO numbers that label ATP maintenance reactions for
+    [`flux_balance_constraints`](@ref).
+    """
+    atp_maintenance_sbos::Vector{Int} = ["SBO:0000630"]
+
+    """
+    SBO numbers that label metabolite demand reactions for
+    [`flux_balance_constraints`](@ref).
+    """
+    demand_sbos::Vector{Int} = ["SBO:0000628"]
+end
+
+"""
+    const configuration
+
+The configuration object. You can change the contents of configuration to
+override the default behavior of some of the functions.
+
+The available options are described by struct [`Configuration`](@ref).
+"""
+const configuration = Configuration()