Parsimonious growth medium

Parsimonious flux balance analysis may be easily combined with the unidirectional reaction splitting tools to create an analysis that minimizes the total nutrient uptake flux. This is beneficial as a complementary, less "refined" but much more computationally efficient alternative to the growth medium optimization (which requires a demanding mixed-integer programming).

For the demonstration, we use the toy E. coli model:

using COBREXA

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

import JSONFBCModels, HiGHS
model = load_model("e_coli_core.json")

JSONFBCModels.JSONFBCModel(#= 95 reactions, 72 metabolites =#)

First, we create a model with extra direction-splitted variables and constraints for the negative exchange reactions:

c = flux_balance_constraints(model, interface = :identifier_prefixes)
c += :in_exchanges^unsigned_negative_contribution_variables(c.interface.exchanges)
c *=
    :out_exchanges^unsigned_positive_contribution_constraints(
        c.interface.exchanges,
        c.in_exchanges,
    )

ConstraintTrees.ConstraintTree with 7 elements:
  :coupling           => ConstraintTrees.ConstraintTree(#= 0 elements =#)
  :flux_stoichiometry => ConstraintTrees.ConstraintTree(#= 72 elements =#)
  :fluxes             => ConstraintTrees.ConstraintTree(#= 95 elements =#)
  :in_exchanges       => ConstraintTrees.ConstraintTree(#= 20 elements =#)
  :interface          => ConstraintTrees.ConstraintTree(#= 3 elements =#)
  :objective          => ConstraintTrees.Constraint(ConstraintTrees.LinearValue…
  :out_exchanges      => ConstraintTrees.ConstraintTree(#= 20 elements =#)

For the parsimonious optimization, we need a good guess on the objective value that we will require; which we obtain from an optimal solution of the original model:

objective_value = flux_balance_analysis(model, optimizer = HiGHS.Optimizer).objective

0.8739215069684305

Parsimonious analyses need a clear specification of a solution tolerance; we use a simple absolute difference bound:

tolerances = [absolute_tolerance_bound(0.001)]

1-element Vector{COBREXA.var"#474#475"{Float64}}:
 #474 (generic function with 1 method)

We can now run a linear "parsimonious intake" analysis as follows:

result = parsimonious_optimized_values(
    c;
    objective = c.objective.value,
    parsimonious_objective = sum_value(c.in_exchanges),
    objective_value,
    tolerances,
    optimizer = HiGHS.Optimizer,
    output = c.in_exchanges,
)

ConstraintTrees.Tree{Float64} with 20 elements:
  :EX_ac_e     => 0.0
  :EX_acald_e  => 0.0
  :EX_akg_e    => 0.0
  :EX_co2_e    => 0.0
  :EX_etoh_e   => 0.0
  :EX_for_e    => 0.0
  :EX_fru_e    => 0.0
  :EX_fum_e    => 0.0
  :EX_glc__D_e => 10.0
  :EX_gln__L_e => 0.0
  :EX_glu__L_e => 0.0
  :EX_h2o_e    => 0.0
  :EX_h_e      => 0.0
  :EX_lac__D_e => 0.0
  :EX_mal__L_e => 0.0
  :EX_nh4_e    => 4.75987
  :EX_o2_e     => 21.7559
  :EX_pi_e     => 3.21122
  :EX_pyr_e    => 0.0
  :EX_succ_e   => 0.0

If suitable, one may also choose to run a quadratic (L2) parsimonious analysis, simply by swapping the objective:

l2_result = parsimonious_optimized_values(
    c;
    objective = c.objective.value,
    parsimonious_objective = sum_value(c.in_exchanges),
    objective_value,
    tolerances,
    optimizer = HiGHS.Optimizer,
    output = c.in_exchanges,
)

ConstraintTrees.Tree{Float64} with 20 elements:
  :EX_ac_e     => 0.0
  :EX_acald_e  => 0.0
  :EX_akg_e    => 0.0
  :EX_co2_e    => 0.0
  :EX_etoh_e   => 0.0
  :EX_for_e    => 0.0
  :EX_fru_e    => 0.0
  :EX_fum_e    => 0.0
  :EX_glc__D_e => 10.0
  :EX_gln__L_e => 0.0
  :EX_glu__L_e => 0.0
  :EX_h2o_e    => 0.0
  :EX_h_e      => 0.0
  :EX_lac__D_e => 0.0
  :EX_mal__L_e => 0.0
  :EX_nh4_e    => 4.75987
  :EX_o2_e     => 21.7559
  :EX_pi_e     => 3.21122
  :EX_pyr_e    => 0.0
  :EX_succ_e   => 0.0

Because of the model simplicity, the solution coincides with the linear one.

To better demonstrate the difference between L1 and L2 solutions that more complex models would exhibit, we arbitrarily choose to relax the tolerance (indirectly allowing non-unique solution):

tolerances = [relative_tolerance_bound(0.5)]

1-element Vector{COBREXA.var"#476#477"{Float64}}:
 #476 (generic function with 1 method)

The linear solution is appropriately reduced:

parsimonious_optimized_values(
    c;
    objective = c.objective.value,
    parsimonious_objective = sum_value(c.in_exchanges),
    objective_value,
    tolerances,
    optimizer = HiGHS.Optimizer,
    output = c.in_exchanges,
)

ConstraintTrees.Tree{Float64} with 20 elements:
  :EX_ac_e     => 0.0
  :EX_acald_e  => 0.0
  :EX_akg_e    => 0.0
  :EX_co2_e    => 0.0
  :EX_etoh_e   => 0.0
  :EX_for_e    => 0.0
  :EX_fru_e    => 0.0
  :EX_fum_e    => 0.0
  :EX_glc__D_e => 6.52119
  :EX_gln__L_e => 0.0
  :EX_glu__L_e => 0.0
  :EX_h2o_e    => 0.0
  :EX_h_e      => 0.0
  :EX_lac__D_e => 0.0
  :EX_mal__L_e => 0.0
  :EX_nh4_e    => 2.38266
  :EX_o2_e     => 9.72181
  :EX_pi_e     => 1.60745
  :EX_pyr_e    => 0.0
  :EX_succ_e   => 0.0

The quadratic solution additionally reduces the "extremes" in the flux profile, since the quadratic cost makes them unfavorable:

parsimonious_optimized_values(
    c;
    objective = c.objective.value,
    parsimonious_objective = squared_sum_value(c.in_exchanges),
    objective_value,
    tolerances,
    optimizer = HiGHS.Optimizer,
    output = c.in_exchanges,
)

ConstraintTrees.Tree{Float64} with 20 elements:
  :EX_ac_e     => 0.0
  :EX_acald_e  => 0.0
  :EX_akg_e    => 0.0
  :EX_co2_e    => 0.0
  :EX_etoh_e   => 0.0
  :EX_for_e    => 0.0
  :EX_fru_e    => 0.0
  :EX_fum_e    => 0.0
  :EX_glc__D_e => 8.12149
  :EX_gln__L_e => 0.0
  :EX_glu__L_e => 0.0
  :EX_h2o_e    => 0.0
  :EX_h_e      => 0.0
  :EX_lac__D_e => 0.0
  :EX_mal__L_e => 0.0
  :EX_nh4_e    => 2.38266
  :EX_o2_e     => 8.12151
  :EX_pi_e     => 1.60745
  :EX_pyr_e    => 0.0
  :EX_succ_e   => 0.0

This page was generated using Literate.jl.