Flux balance analysis

Flux balance analysis (FBA) predicts steady-state fluxes in a metabolic network. For a stoichiometric matrix \(S\) and reaction-flux vector \(v\), FBA imposes

\[ S v = 0, \qquad v_{\mathrm{lower}} \leq v \leq v_{\mathrm{upper}}, \]

and optimizes an objective such as biomass production. These models are commonly used to predict growth, compare media, simulate knockouts, and identify metabolic-engineering strategies.

CORNETO brings FBA into the same optimization framework used for other network-inference problems. It interoperates with COBRApy models and supports standard FBA, while exposing the resulting optimization problem so that constraints, conditions, sparsity, and omics-derived evidence can be combined through a common interface.

COBRApy remains a natural tool for curating metabolic models and running established COBRA workflows. CORNETO is useful when the formulation itself must be extended or several network-inference problems must be expressed together.

From COBRApy to CORNETO

COBRApy concept	CORNETO concept
cobra.Model	metabolic hypergraph
metabolite	vertex
reaction and stoichiometry	directed hyperedge and its coefficients
reaction bounds	edge-flow bounds
objective coefficients	reaction objectives
solution.fluxes	problem.expr.flow.value
optimize a configured model	build and extend an optimization problem

The graph representation preserves reaction identifiers, stoichiometry, default bounds, and gene–protein–reaction annotations imported from the COBRApy model.

Solve the COBRApy model

We start with the familiar E. coli core model distributed with COBRApy.

from contextlib import redirect_stdout
from io import StringIO

import numpy as np
import pandas as pd
from cobra.io import load_model

from corneto.io import cobra_model_to_graph
from corneto.methods import MultiSampleFBA

with redirect_stdout(StringIO()):
    model = load_model("textbook")
biomass_id = "Biomass_Ecoli_core"
cobra_solution = model.optimize()

pd.Series(
    {
        "metabolites": len(model.metabolites),
        "reactions": len(model.reactions),
        "biomass flux": cobra_solution.fluxes[biomass_id],
    },
    name="COBRApy model",
)

metabolites     72.000000
reactions       95.000000
biomass flux     0.873922
Name: COBRApy model, dtype: float64

Convert the model and solve with CORNETO

cobra_model_to_graph converts each reaction into a stoichiometric hyperedge. MultiSampleFBA.build then constructs the mass-balance constraints, applies the imported bounds, and adds the requested reaction objective.

MultiSampleFBA minimizes a weighted objective. Therefore, the coefficient -1 below maximizes biomass. CORNETO’s general optimization API also supports an explicit Direction.MAX when constructing a problem directly, but direction is not an argument of MultiSampleFBA.build.

G = cobra_model_to_graph(model)
reaction_ids = list(G.get_attr_from_edges("id"))
reaction_index = {reaction_id: i for i, reaction_id in enumerate(reaction_ids)}
biomass_idx = reaction_index[biomass_id]

problem = MultiSampleFBA().build(
    G,
    objectives={biomass_id: -1},
)
problem.solve(solver="scipy")

corneto_biomass = float(problem.expr.flow[biomass_idx].value)
comparison = pd.Series(
    {
        "COBRApy": cobra_solution.fluxes[biomass_id],
        "CORNETO": corneto_biomass,
    },
    name="optimal biomass flux",
)

assert np.isclose(comparison["COBRApy"], comparison["CORNETO"], atol=1e-7)
comparison

COBRApy    0.873922
CORNETO    0.873922
Name: optimal biomass flux, dtype: float64

Fluxes follow the graph’s edge order. Mapping that order back to reaction identifiers gives a labeled result analogous to solution.fluxes.

fluxes = pd.Series(
    np.asarray(problem.expr.flow.value),
    index=reaction_ids,
    name="flux",
)
fluxes.loc[[biomass_id, "EX_glc__D_e", "EX_o2_e", "ATPM"]]

Biomass_Ecoli_core     0.873922
EX_glc__D_e          -10.000000
EX_o2_e              -21.799493
ATPM                   8.390000
Name: flux, dtype: float64

Extend the optimization problem

A built CORNETO problem is still editable. For example, the following constraint limits oxygen uptake to 10 mmol gDW\(^{-1}\) h\(^{-1}\). In this model, exchange uptake is a negative flux, so a less-negative lower limit permits less uptake.

oxygen_idx = reaction_index["EX_o2_e"]

oxygen_limited = MultiSampleFBA().build(
    G,
    objectives={biomass_id: -1},
)
oxygen_limited += oxygen_limited.expr.flow[oxygen_idx] >= -10
oxygen_limited.solve(solver="scipy")

oxygen_limited_biomass = float(oxygen_limited.expr.flow[biomass_idx].value)
assert oxygen_limited_biomass < corneto_biomass

pd.Series(
    {
        "default medium": corneto_biomass,
        "oxygen-limited": oxygen_limited_biomass,
    },
    name="optimal biomass flux",
)

default medium    0.873922
oxygen-limited    0.559051
Name: optimal biomass flux, dtype: float64

Solve multiple media conditions

build_many creates one flux vector per named condition. Objectives and reaction bounds use an outer mapping with the same condition names. Conditions are solved in one problem, which also makes it possible for advanced formulations to couple them.

condition_objectives = {
    "glucose_rich": {biomass_id: -1},
    "glucose_limited": {biomass_id: -1},
}
condition_bounds = {
    "glucose_rich": {
        "EX_glc__D_e": (-10.0, 1000.0),
        "EX_o2_e": (-20.0, 1000.0),
    },
    "glucose_limited": {
        "EX_glc__D_e": (-2.0, 1000.0),
        "EX_o2_e": (-20.0, 1000.0),
    },
}

media_problem = MultiSampleFBA().build_many(
    G,
    objectives=condition_objectives,
    reaction_bounds=condition_bounds,
)
media_problem.solve(solver="scipy")

condition_names = list(condition_objectives)
reported_reactions = [biomass_id, "EX_glc__D_e", "EX_o2_e"]
media_fluxes = pd.DataFrame(
    {
        condition: {
            reaction_id: float(media_problem.expr.flow[reaction_index[reaction_id], column].value)
            for reaction_id in reported_reactions
        }
        for column, condition in enumerate(condition_names)
    }
)

assert media_problem.expr.flow.value.shape == (G.num_edges, len(condition_names))
assert media_fluxes.loc[biomass_id, "glucose_limited"] < media_fluxes.loc[biomass_id, "glucose_rich"]
media_fluxes

	glucose_rich	glucose_limited
Biomass_Ecoli_core	0.832614	0.140604
EX_glc__D_e	-10.000000	-2.000000
EX_o2_e	-20.000000	-5.853994

Sparse FBA

CORNETO can associate each reaction with a binary indicator that records whether its flux is nonzero. With lambda_reg > 0, it penalizes the number of active reactions. Because of these indicators, sparse FBA is a mixed-integer linear program (MILP), rather than the linear program used by the standard mathematical FBA formulation.

Here we require at least 90% of the optimal biomass and optimize biomass together with a penalty on the number of active reactions. lambda_reg controls that trade-off; larger values give reaction count more influence. The biomass floor prevents sparsity from being achieved by suppressing growth.

biomass_floor = 0.90 * corneto_biomass
sparse_problem = MultiSampleFBA(lambda_reg=0.1).build(
    G,
    objectives={biomass_id: -1},
    reaction_bounds={biomass_id: (biomass_floor, None)},
)
sparse_problem.solve(solver="highs")

sparse_biomass = float(sparse_problem.expr.flow[biomass_idx].value)
standard_active = int(np.count_nonzero(np.abs(problem.expr.flow.value) > 1e-6))
sparse_active = int(np.count_nonzero(np.abs(sparse_problem.expr.flow.value) > 1e-6))

assert sparse_biomass >= biomass_floor - 1e-7
assert sparse_active <= standard_active

pd.Series(
    {
        "biomass floor": biomass_floor,
        "sparse biomass": sparse_biomass,
        "active reactions in standard solution": standard_active,
        "active reactions in sparse solution": sparse_active,
    }
)

biomass floor                             0.786529
sparse biomass                            0.814298
active reactions in standard solution    48.000000
active reactions in sparse solution      45.000000
dtype: float64

This reaction-count criterion is not the same as parsimonious FBA (pFBA), which conventionally minimizes total flux after fixing the primary objective. Sparse FBA instead penalizes how many reactions are active.

For several conditions, CORNETO can apply structured sparsity to the union of active reactions. A reaction shared by several conditions is then counted once, encouraging compact shared metabolic programs while retaining a separate flux vector for every condition. Continue with Multi-condition FBA for the general formulation and a worked example, followed by gene expression integration for context-specific metabolic inference.

The examples above use the method-specific build and build_many interfaces. The general Data interface remains available for workflows that require custom feature metadata or already represent measurements as CORNETO data objects.