Bio Systems Biology Flux Balance Analysis

import cobra

model = cobra.io.load_model('textbook')   # E. coli core (e_coli_core), 95 reactions, ships with cobra
model = cobra.io.load_model('iJO1366')    # genome-scale E. coli, 2583 reactions

model = cobra.io.read_sbml_model('model.xml')   # SBML (the standard exchange format)
model = cobra.io.load_json_model('model.json')  # COBRA JSON

# Curated genome-scale models: http://bigg.ucsd.edu/models (King 2016). Record the model
# version; predictions are only comparable within the same model release.

import cobra

model = cobra.io.load_model('textbook')

solution = model.optimize()
print(f'Objective (growth): {solution.objective_value:.4f} /h  status: {solution.status}')
print('Objective reaction:', str(model.objective.expression).split('*')[1].split()[0])

# A growth rate is interpretable ONLY against a stated medium and biomass reaction.
# Compare to a measured doubling time (mu = ln2 / t_double) rather than to a fixed
# "fast/slow" scale; absolute values are organism- and biomass-definition-specific.

def set_minimal_medium(model, carbon_source='EX_glc__D_e', carbon_uptake=10):
    '''Close all uptake, then open a defined minimal medium.

    carbon_uptake in mmol/gDW/h; 10 is the standard E. coli aerobic glucose rate (iJO1366).
    '''
    for rxn in model.exchanges:
        rxn.lower_bound = 0

    minimal = {'EX_o2_e': 1000, 'EX_h2o_e': 1000, 'EX_h_e': 1000, 'EX_nh4_e': 1000,
               'EX_pi_e': 1000, 'EX_so4_e': 1000, 'EX_k_e': 1000, 'EX_mg2_e': 1000}
    for ex_id, uptake in minimal.items():
        if ex_id in model.reactions:
            model.reactions.get_by_id(ex_id).lower_bound = -uptake
    if carbon_source in model.reactions:
        model.reactions.get_by_id(carbon_source).lower_bound = -carbon_uptake
    return model

# The with-block reverts all bound changes on exit, so comparisons never leak state.
for cs in ['EX_glc__D_e', 'EX_ac_e', 'EX_succ_e']:
    with model:
        set_minimal_medium(model, carbon_source=cs)
        print(f'{cs}: growth = {model.slim_optimize():.4f}')   # slim_optimize returns the objective float only (fast)

from cobra.flux_analysis import flux_variability_analysis

# Range each reaction can carry while holding the objective at (fraction_of_optimum) of the max.
fva = flux_variability_analysis(model, fraction_of_optimum=1.0)

# fraction_of_optimum < 1 relaxes the objective and reveals the alternative-optima span.
fva90 = flux_variability_analysis(model, fraction_of_optimum=0.9)

# loopless=True removes thermodynamically infeasible internal loops from the ranges (slower).
fva_ll = flux_variability_analysis(model, loopless=True)

# A reaction with maximum == minimum is fully determined; a wide range means the single
# FBA value for it was arbitrary. Blocked reactions have min == max == 0.
fva['range'] = fva['maximum'] - fva['minimum']
fva['blocked'] = fva['range'].abs() < 1e-9

from cobra.flux_analysis import pfba

# Among all optima, pFBA returns the flux vector minimizing the sum of absolute fluxes,
# an Occam's-razor proxy for minimal total enzyme cost (Lewis 2010). It is a principled
# single representative of the optimal face; it is NOT more "true" than the face itself,
# and it hugs the polytope boundary (a min-flux vertex), so report the FVA range alongside it
# when the internal flux values matter.
pfba_solution = pfba(model)
print(f'FBA total flux : {model.optimize().fluxes.abs().sum():.1f}')
print(f'pFBA total flux: {pfba_solution.fluxes.abs().sum():.1f}')

from cobra.flux_analysis import loopless_solution

# Projects an FBA solution onto a loopless one: no net flux around a closed cycle that
# lacks a thermodynamic driving force (Schellenberger 2011). Internal cycles are a common
# artifact of reversible reactions and gap-filling and inflate apparent flux magnitudes.
loopless = loopless_solution(model)

from cobra.sampling import sample

# n samples; method 'optgp' (parallel) or 'achr'. Thinning reduces autocorrelation.
samples = sample(model, n=1000, method='optgp', thinning=100, seed=1)
print(samples['PFK'].describe())   # per-reaction distribution, e.g. median and IQR

# To sample only high-growth states, constrain the objective first (e.g. biomass >= 0.9*max)
# inside a `with model:` block, then sample. Check mixing before trusting the distribution
# (multiple chains/seeds should agree); short chains give correlated, misleading samples.