How to perform relaxedFBA with cobrapy?

[cdiener]

This is not implemented in cobrapy. But it is pretty straightforward. Here is an example on how to do it only for the bounds:

from cobra.io import load_model
from numpy import Inf

mod = load_model("textbook")
print(mod.summary())

mod.reactions.Biomass_Ecoli_core.lower_bound = 1
print(mod.optimize())  # infeasible

def relaxed_bounds_fba(model, exclude=None):
    """Relax the bounds to achieve feasibility.

    Arguments
    ---------
    model : cobra.Model
        A cobra model with an infeasible constraint.
    exclude : list of str
        Reaction IDs which should not be relaxed.

    Returns
    -------
    A relaxed model solution.
    """
    relaxations = []
    constraints = []

    with model:
        for rxn in mod.reactions:
            if rxn.id in exclude:
                continue
            p = mod.problem.Variable(f"p_{rxn.id}", lb=0)  # lower bounds relaxation
            q = mod.problem.Variable(f"q_{rxn.id}", lb=0)  # upper bound relaxation
            const = mod.problem.Constraint(                # relaxed flux
                rxn.flux_expression + p - q,
                lb=rxn.lower_bound,
                ub=rxn.upper_bound,
                name=f"relaxed_{rxn.id}"
            )
            relaxations.extend([p, q])
            constraints.append(const)
            rxn.bounds = -Inf, Inf                         # deactivate the "old" bounds

        mod.add_cons_vars(relaxations + constraints)
        mod.objective = sum(relaxations)                   # minimize relaxation
        mod.objective_direction = "min"
        relaxed_sol = mod.optimize()
    return relaxed_sol

rs = relaxed_bounds_fba(mod, exclude=["Biomass_Ecoli_core"])
print(mod.summary(solution=rs))

Which would give you the following output:

Objective
=========
1.0 Biomass_Ecoli_core = 0.8739215069684302

Uptake
------
Metabolite    Reaction  Flux  C-Number  C-Flux
  glc__D_e EX_glc__D_e    10         6 100.00%
     nh4_e    EX_nh4_e 4.765         0   0.00%
      o2_e     EX_o2_e  21.8         0   0.00%
      pi_e     EX_pi_e 3.215         0   0.00%

Secretion
---------
Metabolite Reaction   Flux  C-Number  C-Flux
     co2_e EX_co2_e -22.81         1 100.00%
     h2o_e EX_h2o_e -29.18         0   0.00%
       h_e   EX_h_e -17.53         0   0.00%

/home/cdiener/code/cobrapy/src/cobra/util/solver.py:554: UserWarning: Solver status is 'infeasible'.
  warn(f"Solver status is '{status}'.", UserWarning)
<Solution infeasible at 0x7ff84cd846d0>
Objective
=========
1.0 Biomass_Ecoli_core = 1.0

Uptake
------
Metabolite    Reaction  Flux  C-Number  C-Flux
  glc__D_e EX_glc__D_e 11.38         6 100.00%
     nh4_e    EX_nh4_e 5.453         0   0.00%
      o2_e     EX_o2_e 24.54         0   0.00%
      pi_e     EX_pi_e 3.679         0   0.00%

Secretion
---------
Metabolite Reaction   Flux  C-Number  C-Flux
     co2_e EX_co2_e  -25.7         1 100.00%
     h2o_e EX_h2o_e -32.98         0   0.00%
       h_e   EX_h_e -20.06         0   0.00%

So it relaxes the imports (there may be several equally good relaxations so you might not see the same result).

For the full relaxedFBA you would have to add relaxations to the equality constraints as well and allow specifying the weights, but it would be the same strategy.