4. Variability Analysis

4.1. Flux Variability Analysis

Flux variability analysis (FVA) can be used to assess the possible range of metabolic fluxes within a metabolic network while still achieving a certain level of optimality in the objective function. It essentially determines the flexibility of the metabolic network conditioned by the network’s objective. The etfba package allows for FVA along with additional constraints on enzyme protein allocation and thermodynamics, implemented separately or together. The complete form (ETFVA) can be expressed as:

Here, \(\gamma\) in the range [0, 1] determines the fraction of the objective’s optimum that should be achieved. Try a smaller \(\gamma\) if the solver cannot find a feasible solution.

Let’s begin with the basic flux variability analysis, focusing solely on the mass balance constraints of metabolites. You can perform FVA using the following code:

from etfba import Model

model_file = '../../models/e_coli/etfba_iML1515.bin'
model = Model.load(model_file)
print(model)

model iML1515 with 2712 reactions and 1877 metabolites

objective = {'BIOMASS_Ec_iML1515_core_75p37M': 1}
flux_bound = (0, 1000)
spec_flux_bound = {'ATPM': (6.86, 1000)}
preset_flux = {'EX_glc__D_e_b': 10, 'FHL': 0}

res = model.evaluate_variability(
    'fva',
    objective=objective,
    obj_value=0.877,
    gamma=0.99,
    flux_bound=flux_bound,
    spec_flux_bound=spec_flux_bound,
    preset_flux=preset_flux
).solve(solver='gurobi', n_jobs=36)

Note: It is highly recommended to run variability analysis in parallel jobs by specifying the “n_jobs” argument, especially when dealing with large-scale models..

The estimated feasible range of fluxes can be accessed with the flux_ranges attribute, which can be further saved directly using its save method.

for rxnid, flux_range in list(res.flux_ranges.items())[:10]:
    print(f'{rxnid:>10} [{flux_range[0]:>6.3f} {flux_range[1]:>6.3f}]')

    CYTDK2 [ 0.000  2.331]
      XPPT [ 0.000  1.165]
     HXPRT [ 0.000  1.165]
     NDPK5 [-1.576  1.186]
    SHK3Dr [ 0.331  0.388]
     NDPK6 [ 0.000  0.643]
     NDPK8 [-1.576  1.185]
    DHORTS [-0.445 -0.287]
     OMPDC [ 0.287  0.445]
    PYNP2r [-2.331  1.188]

Next, we can extend the analysis to include constraints on enzyme protein allocation and thermodynamics, namely EFVA and TFVA, respectively. Here’s how you can perform EFVA:

import pandas as pd

kcat_file = '../../models/e_coli/kcats.xlsx'
mw_file = '../../models/e_coli/mws.xlsx'
dgpm_file = '../../models/e_coli/dgpms.xlsx'
kcats = pd.read_excel(kcat_file, header=None, index_col=0).squeeze()
mws = pd.read_excel(mw_file, header=None, index_col=0).squeeze()
dgpms = pd.read_excel(dgpm_file, header=None, index_col=0).squeeze()

# set reactions with available kcats and MWs
eff_rxns = []
for rxnid, rxn in model.reactions.items():
    if rxn.rev:
        if rxnid+'_f' in kcats.index and rxnid+'_b' in kcats.index and rxnid in mws.index:
            rxn.forward_kcat = kcats[rxnid+'_f']
            rxn.backward_kcat = kcats[rxnid+'_b']
            rxn.molecular_weight = mws[rxnid]
            eff_rxns.append(rxnid)
    else:
        if rxnid in kcats.index and rxnid in mws.index:
            rxn.forward_kcat = kcats[rxnid]
            eff_rxns.append(rxnid)
enz_ub = 0.15

# set reactions with available ΔG'm
ex_rxns = []
for rxnid, rxn in model.reactions.items():
    if rxnid in dgpms.index:
        rxn.standard_gibbs_energy = dgpms[rxnid]
    else:
        ex_rxns.append(rxnid)

# perform EFVA
preset_flux = {'FHL': 0}
res = model.evaluate_variability(
    'efva',
    objective=objective,
    obj_value=0.866,
    gamma=0.99,
    flux_bound=flux_bound,
    spec_flux_bound=spec_flux_bound,
    preset_flux=preset_flux,
    inc_enz_cons=eff_rxns,
    enz_prot_lb=enz_ub
).solve(solver='gurobi', n_jobs=36)

And TFVA:

# perform TFVA
preset_flux = {'EX_glc__D_e_b': 10, 'FHL': 0}
conc_bound = (0.0001, 100)
spec_conc_bound = {
    'o2_c': (0.0001, 0.0082),
    'co2_c': (0.1, 100)
}
preset_conc = {
    'glc__D_p': 20,
    'pi_p': 56,
    'so4_p': 3,
    'nh4_p': 19,
    'na1_p': 160,
    'k_p': 22,
    'fe2_p': 62
}
res = model.evaluate_variability(
    'tfva',
    objective=objective,
    obj_value=0.877,
    gamma=0.99,
    flux_bound=flux_bound,
    conc_bound=conc_bound,
    spec_flux_bound=spec_flux_bound,
    spec_conc_bound=spec_conc_bound,
    preset_flux=preset_flux,
    preset_conc=preset_conc,
    ex_thermo_cons=ex_rxns
).solve(solver='gurobi', n_jobs=36)

Moreover, we can conduct flux variability analysis with comprehensive constraints on enzyme protein allocation and thermodynamics, i.e., ETFVA.

# perform ETFVA
preset_flux = {'FHL': 0}
res = model.evaluate_variability(
    'etfva',
    objective=objective,
    obj_value=0.864,
    gamma=0.99,
    flux_bound=flux_bound,
    conc_bound=conc_bound,
    spec_flux_bound=spec_flux_bound,
    spec_conc_bound=spec_conc_bound,
    preset_flux=preset_flux,
    preset_conc=preset_conc,
    ex_thermo_cons=ex_rxns,
    inc_enz_cons=eff_rxns,
    enz_prot_lb=enz_ub
).solve(solver='gurobi', n_jobs=36)

4.2. Enzyme Protein Cost Variability Analysis

Similar to flux variability analysis, the variability of enzyme protein costs in enzymatic reactions can also be evaluated. Below is the complete form of enzyme protein variability analysis (TEVA), where the constraints are identical to those in ETFVA.

If only the constraint of enzyme protein is considered, basic EVA can be performed using the following code:

preset_flux = {'FHL': 0}
res = model.evaluate_variability(
    'eva',
    objective=objective,
    obj_value=0.866,
    gamma=0.99,
    flux_bound=flux_bound,
    spec_flux_bound=spec_flux_bound,
    preset_flux=preset_flux,
    inc_enz_cons=eff_rxns,
    enz_prot_lb=enz_ub
).solve(solver='gurobi', n_jobs=36)

The estimated feasible ranges of enzyme protein costs can be found in the protein_cost_ranges attribute. You can use the save method of this attribute to save the results.

pro_cost_sorted = sorted(
    res.protein_cost_ranges.items(),
    key=lambda item: item[1][0],
    reverse=True
)
for rxnid, pro_cost_range in pro_cost_sorted[:10]:
    print(f'{rxnid:>10} [{pro_cost_range[0]:>6.5f} {pro_cost_range[1]:>6.5f}]')

       ENO [0.00875 0.01220]
      GAPD [0.00612 0.00885]
     KARA2 [0.00499 0.00641]
     KARA1 [0.00490 0.00632]
       PGK [0.00489 0.00736]
      METS [0.00474 0.00548]
       FBA [0.00357 0.01068]
       PGM [0.00289 0.00498]
    GHMT2r [0.00249 0.00504]
  ATPS4rpp [0.00196 0.00346]

Below shows the code of enzyme protein variability analysis under the complete constraints of enzyme protein allocation and thermodynamics, i.e., TEVA.

res = model.evaluate_variability(
    'teva',
    objective=objective,
    obj_value=0.864,
    gamma=0.98,
    flux_bound=flux_bound,
    conc_bound=conc_bound,
    spec_flux_bound=spec_flux_bound,
    spec_conc_bound=spec_conc_bound,
    preset_flux=preset_flux,
    preset_conc=preset_conc,
    ex_thermo_cons=ex_rxns,
    inc_enz_cons=eff_rxns,
    enz_prot_lb=enz_ub
).solve(solver='gurobi', n_jobs=36)

4.3. Thermodynamic Variability Analysis

Thermodynamic variability analysis can be used to evaluate the feasible ranges of reaction Gibbs energy change. Its complete form (ETVA) can be represented as below, with constraints identical to those in ETFVA.

The basic version of TVA without enzyme protein allocation constraint can be conducted with the following code:

preset_flux = {'EX_glc__D_e_b': 10, 'FHL': 0}
res = model.evaluate_variability(
    'tva',
    objective=objective,
    obj_value=0.877,
    gamma=0.99,
    flux_bound=flux_bound,
    conc_bound=conc_bound,
    spec_flux_bound=spec_flux_bound,
    spec_conc_bound=spec_conc_bound,
    preset_flux=preset_flux,
    preset_conc=preset_conc,
    ex_thermo_cons=ex_rxns
).solve(solver='gurobi', n_jobs=36)

To inspect the bounds of feasible \(\Delta G'\), one can use the gibbs_energy_ranges attribute, which also has a save method for direct saving.

for rxnid, dgp_range in list(res.gibbs_energy_ranges.items())[:10]:
    print(f'{rxnid:>10} [{dgp_range[0]:>6.1f} {dgp_range[1]:>6.1f}]')

    CYTDK2 [-102.0   44.3]
      XPPT [ -93.9   34.7]
     HXPRT [ -97.7   31.0]
     NDPK5 [ -60.0   41.8]
    SHK3Dr [ -56.7   -0.0]
     NDPK6 [ -63.4   38.4]
     NDPK8 [ -59.6   42.3]
    DHORTS [   0.0   17.1]
     OMPDC [ -69.4   -0.0]
    PYNP2r [ -68.5   67.5]

Variability of reactions Gibbs energy under complete constraints of enzyme protein allocation and thermodynamics, i.e., ETVA, can be conducted as follows:

preset_flux = {'FHL': 0}
res = model.evaluate_variability(
    'etva',
    objective=objective,
    obj_value=0.864,
    gamma=0.98,
    flux_bound=flux_bound,
    conc_bound=conc_bound,
    spec_flux_bound=spec_flux_bound,
    spec_conc_bound=spec_conc_bound,
    preset_flux=preset_flux,
    preset_conc=preset_conc,
    ex_thermo_cons=ex_rxns,
    inc_enz_cons=eff_rxns,
    enz_prot_lb=enz_ub
).solve(solver='gurobi', n_jobs=36)