GitHub

Making adjustments to the constraint system

In the previous example about model adjustments, we noted that some constraint systems may be too complex to be changed within the limits of the usual FBC model view, and we may require a sharper tool to do the changes we need. This example shows how to do that by modifying the constraint systems that are generated within COBREXA to represent the metabolic model contents.

using COBREXA

download_model(
    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
    "e_coli_core.json",
    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)

import JSONFBCModels
import HiGHS

model = load_model("e_coli_core.json") # flux balance type model
JSONFBCModels.JSONFBCModel(#= 95 reactions, 72 metabolites =#)

Background: Constraint trees

COBREXA uses ConstraintTrees to represent model structures internally. This framework provides a powerful unified interface over all constraints and variables in the model, making its manipulation much more convenient.

import ConstraintTrees as C

In general, constraint-based models use fluxes as variables, and all the constraints are in terms of them (or derived quantities). We can get a constraint tree for the usual flux-balance-style models quite easily:

ct = flux_balance_constraints(model)
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 4 elements:
  :coupling           => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 0 …
  :flux_stoichiometry => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 72…
  :fluxes             => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 95…
  :objective          => ConstraintTrees.Constraint(ConstraintTrees.LinearValue…

The fluxes are represented by constraints for individual variables:

ct.fluxes
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 95 elements:
  :ACALD                    => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ACALDt                   => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ACKr                     => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ACONTa                   => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ACONTb                   => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ACt2r                    => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ADK1                     => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :AKGDH                    => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :AKGt2r                   => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ALCD2x                   => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ATPM                     => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ATPS4r                   => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :BIOMASS_Ecoli_core_w_GAM => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :CO2t                     => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :CS                       => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :CYTBD                    => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :D_LACt2                  => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ENO                      => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  :ETOHt2r                  => ConstraintTrees.Constraint(ConstraintTrees.Linea…
  ⋮                         => ⋮

The "mass balance" is represented as equality constraints:

ct.flux_stoichiometry
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 72 elements:
  Symbol("13dpg_c") => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  Symbol("2pg_c")   => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  Symbol("3pg_c")   => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  Symbol("6pgc_c")  => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  Symbol("6pgl_c")  => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :ac_c             => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :ac_e             => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :acald_c          => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :acald_e          => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :accoa_c          => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :acon_C_c         => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :actp_c           => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :adp_c            => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :akg_c            => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :akg_e            => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :amp_c            => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :atp_c            => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :cit_c            => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  :co2_c            => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#…
  ⋮                 => ⋮

The objective is represented as a "transparent reference" to the variables that specify the biomass. Notice that it has no bound (thus it's technically not a constraint, just a "label" for something that has a sensible semantic and can be constrained or optimized).

ct.objective
ConstraintTrees.Constraint(ConstraintTrees.LinearValue([25], [1.0]), nothing)

Customizing the model

New values and constraints can be easily created from the existing ones. For example, this is a total flux through exchanges of typical fermentation products:

total_fermentation = ct.fluxes.EX_ac_e.value + ct.fluxes.EX_etoh_e.value
ConstraintTrees.LinearValue([44, 48], [1.0, 1.0])

With the value in hand, we can constraint it (enforcing that the model outputs at least some fermentation products):

fermentation_constraint = C.Constraint(total_fermentation, (10.0, 1000.0))
ConstraintTrees.Constraint(ConstraintTrees.LinearValue([44, 48], [1.0, 1.0]), ConstraintTrees.Between(10.0, 1000.0))

We can assign a name to the constraint, creating a small (singleton) constraint tree:

:fermentation^fermentation_constraint
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 1 element:
  :fermentation => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ..…

Named constraints can be freely combined, and we combine our new constraint with the whole original constraint tree, getting a differently constrained system:

fermenting_ct = ct * :fermentation^fermentation_constraint
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 5 elements:
  :coupling           => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 0 …
  :fermentation       => ConstraintTrees.Constraint(ConstraintTrees.LinearValue…
  :flux_stoichiometry => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 72…
  :fluxes             => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 95…
  :objective          => ConstraintTrees.Constraint(ConstraintTrees.LinearValue…
What if I need more complicated changes?

Almost all analysis functions have an associated constraint-building function that internally builds the constraint system, for example the parsimonious_flux_balance_analysis is implemented with parsimonious_flux_balance_constraints, which can be used just as flux_balance_constraints here. Additionally, to reach various custom goals, it is recommended to re-use and modify the source of the functions – use the macro @edit, such as @edit parsimonious_flux_balance_constraints, to get a working source code that serves well as a base for implementing new constraint systems.

Constraint trees can be "solved", simply by choosing the objective and sending them to the appropriate function. Here, optimized_values rewrites the constraints into a JuMP model, which is subsequently solved and the solved variables are transformed back into semantically labeled values, in the same structure as the original constraint tree.

solution = optimized_values(
    fermenting_ct,
    objective = fermenting_ct.objective.value,
    optimizer = HiGHS.Optimizer,
)
ConstraintTrees.Tree{Float64} with 5 elements:
  :coupling           => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
  :fermentation       => 10.0
  :flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
  :fluxes             => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :objective          => 0.633738

Models that can not be solved (for any reason) would instead return nothing. We demonstrate that by breaking the bounds of the original constraint trees to an unsolvable state:

ct.fluxes.ATPM.bound = C.Between(1000.0, 10000.0)

solution = optimized_values(ct, objective = ct.objective.value, optimizer = HiGHS.Optimizer)

print(solution)
nothing

Several functions exist to simplify the construction of more complicated constraints. See the reference documentation for generic constraint builders for details.

This page was generated using Literate.jl.