7. Tutorial: Create iML1515 TFA model¶
The thermodynamics (TD) model imposes additional constraints on the metabolic model, with flux directions being coupled to the transformed Gibbs energy of reaction. These transformations are with respect to compartmental pH and ionic strength and aligned to the calculations performed in the pyTFA package (Salvy et al., 2019). The standard Gibbs energies of formation, which are prerequisites for the calculations, are obtained from the thermodynamics database (Alberty, 2003; Jankowski et al., 2008) utilized by the pyTFA package (https://github.com/EPFL-LCSB/pytfa/blob/master/data/thermo_data.thermodb). This TD database contains data records for more than 18,000 molecules, which are referenced by SEED identifiers (Seaver et al., 2021). The records include molecule names, chemical formula, electrical charge, Gibbs energy of formation estimated using the group contribution method (Mavrovouniotis, 1990), the composition of cues (groups) making up the molecule, and acid dissociation constants. The TD database also includes the list of cues used in estimating the Gibbs energy of formation, together with their estimated Gibbs energy contribution and the estimation error.
The implementation of TD constraints is contingent upon the availability of valid TD data records for all participating metabolites. Consequently, the modeler must ensure that the majority of the metabolites in the model have a corresponding TD data record in the TD database. The mapping of model metabolites to TD database entries can be achieved through various methods. One approach involves manual mapping via the table modify_td_sids in the TFA configuration file. This approach
necessitates the definition of a regular expression pattern (mid_regex_pattern of the general table) to extract metabolite identifiers from model species identifiers. Alternatively, the name of the model species, in lower case, can be searched in the TD database. In the event of a failure to locate a match, a sorted list of SEED compound identifiers, extracted from the model species MIRIAM annotation data, is examined against entries in the TD database. In the event that a matching TD
data record is identified, a final verification is conducted to ensure that the TD data record contains valid data and that the chemical formula and charge are consistent with the data configured for the model species. A discrepancy in protonation state with adjusted charge is permissible.
The workflow for generating TD extended models and optimizing these models aligns with the workflows for enzyme constraint and resource balanced constraint models. TD extended models will be exported in SBML-compliant format, thereby facilitating model sharing and downstream processing by SBML-compliant tools. The optimization of the models will be done using cobrapy and alternatively using the gurobipy interface. When employing cobrapy, it is imperative to utilize a solver that supports mixed integer linear programming (MILP), and it is recommended to use a commercial solver such as IBM CPLEX or Gurobi.
The extension of a model through the incorporation of thermodynamic constraints necessitates the establishment of a TFA configuration file. The process commences with the creation of an XbaModel, derived from a genome-scale metabolic model and configured without additional data. The XbaModel is then provided as a reference to a TfaModel, which, based on a preliminary TFA configuration file, adds thermodynamic constraints to yield a preliminary TD constraint model. It is important to note that the optimization of this preliminary TD constraint model may result in the identification of infeasible solutions or a substantial reduction in solution spaces. Constraints must be relaxed in such cases, as error estimates may be too narrow, estimated energies may be incorrectly determined, and thermodynamics calculations may be inexact. Additionally, issues inherent to the metabolic model may arise if multistep enzymatic reactions with energy conservation inside the enzyme are split into individual, independent reactions. Employing a model relaxation step that utilizes a temporary TD constraint “slack” model can ensure that predicted growth rates are in proximity to the predicted growth rates for the model devoid of TD constraints. In this workflow, we also demonstrate how to constrain selected metabolite concentrations to measured values (Buckstein Michael et al., 2008).
Upon completion of this tutorial, the reader will find a description of the variables and constraints that have been incorporated into the genome-scale metabolic model, as well as the technical details about the calculations performed.
Peter Schubert, Heinrich-Heine University Duesseldorf, Institute for Computational Cell Biology (Prof. Dr. M. Lercher), March, 2025
Step 1: Initial setup¶
In order to facilitate thermodynamics modeling, it is necessary to import the TfaModel. Furthermore, in order to optimize the TFA model, we use FbaOptimization. To analyze the results, we use FbaResults. ] Furthermore, the dictionary metabolite_molar_concs is created, with minimum and maximum values based on average nucleotide concentrations determined in E. coli (Buckstein Michael et al., 2008). This data will be used to constrain selected metabolite concentrations in the TFA
model. As the focus remains on the level of FBA modeling, proteomics data is not imported.
[1]:
# Required imports and model names
import os
import re
import time
from collections import defaultdict
import numpy as np
import pandas as pd
import cobra
from f2xba import XbaModel, TfaModel
from f2xba import FbaOptimization, FbaResults
from f2xba.utils.mapping_utils import load_parameter_file, write_parameter_file
fba_model = 'iML1515'
target_model = 'iML1515_TFA'
reference_cond = 'Glucose'
biomass_rid = 'BIOMASS_Ec_iML1515_core_75p37M'
# Create media conditions
media_grs = {'Acetate': ['ac', 0.29], 'Glycerol': ['glyc', 0.47], 'Fructose': ['fru', 0.54],
'L-Malate': ['mal__L', 0.55], 'Glucose': ['glc__D', 0.66], 'Glucose 6-Phosphate': ['g6p', 0.78]}
base_medium = ['ca2', 'cbl1', 'cl', 'co2', 'cobalt2', 'cu2', 'fe2', 'fe3', 'h2o', 'h', 'k', 'mg2',
'mn2', 'mobd', 'na1', 'nh4', 'ni2', 'o2', 'pi', 'sel', 'slnt', 'so4', 'tungs', 'zn2']
conditions = {}
exp_grs = {}
for cond, (carbon_sid, exp_gr )in media_grs.items():
conditions[cond] = {f'EX_{sidx}_e': 1000.0 for sidx in base_medium}
conditions[cond][f'EX_{carbon_sid}_e'] = 10.0 # << carbon source uptake in FBA limited at 10.0 mmol/gDWh
exp_grs[cond] = exp_gr
print(f'{len(conditions)} minimal media conditions created for simulation')
# Metabolite concentrations to be fixed with a margin (Buckstein, 2008)
nt_concs_µmol_per_l = {
'atp_c': 3560, 'adp_c': 116, 'datp_c': 181, 'ctp_c': 325, 'dctp_c': 184,
'gtp_c': 1660, 'gdp_c': 203, 'dgtp_c': 92, 'utp_c': 667, 'udp_c': 54,
'dttp_c': 256, 'ppgpp_c': 113, 'accoa_c': 1390}
factor = 2.0
metabolite_molar_concs = {sid: (µmol_per_l * 1e-6/factor, µmol_per_l*1e-6*factor)
for sid, µmol_per_l in nt_concs_µmol_per_l.items()}
6 minimal media conditions created for simulation
download thermodynamics database¶
In order to incorporate thermodynamic constraints into the model, it is necessary to obtain information on Gibbs energies of formation for metabolites. This information is extracted from a thermodynamics database obtained from the pyTFA (Salvy et al., 2019) repository (https://github.com/EPFL-LCSB/pytfa/raw/refs/heads/master/data/thermo_data.thermodb). This database contains Gibbs energies of formation for metabolites, which have been calculated using the group contribution method (Alberty, 2003; Jankowski et al., 2008).
[2]:
# Download thermodynamics database
import wget
base_url = 'https://github.com/EPFL-LCSB/pytfa/raw/refs/heads/master/data/'
fname = 'thermo_data.thermodb'
full_fname = os.path.join('data', fname)
if not os.path.exists(full_fname):
url = base_url + fname
print(f'\ndownloading {fname} from {url}')
wget.download(url, out=full_fname)
print(f'thermodynamics database (pyTFA): {full_fname}')
thermodynamics database (pyTFA): data/thermo_data.thermodb
Step 2: Create XBA and preliminary TFA configuration files¶
It is not necessary to apply any changes to the genome-scale metabolic model; therefore, a XBA configuration file will not be provided.
The configuration data pertaining to thermodynamics constraints is collated in multiple tables and exported to a TFA configuration file. The tables designated as general and td_compartments are obligatory, while the remaining tables are optional. Within the general table, the file name of the TD database must be specified in the thermo_data_fname parameter. The corresponding file must have the same structure as the file ‘thermo_data.thermodb’ used in the pyTFA package. As we
intend to hard link several metabolite identifiers to TD data records using the table ‘modify_td_sids’, we have to assign a regular expression pattern to the parameter mid_regex_pattern. This pattern will be used to derive metabolite identifiers from model species identifiers. The table td_compartments contains compartment-specific configuration parameters, including the compartment pH, the ionic strength (ionic_strength_M) in mol/L, the minimum and maximum metabolite concentrations
(c_min_M, c_max_M) in mol/L, and the membrane potentials (<other_cid>_mV) in mV. The latter are calculated by subtracting the local from the <other_cid> compartment electrical potential. The table modify_thermo_data is optional; its utilization is for modifying data in the TD database, which was found to be inconsistent. A specific TD data record is identified by its identifier (td_sid), and the value of a certain attribute can be changed. The notes field is an
optional field that is used to annotate modifications. The optional table modify_td_sids facilitates the hard linking of selected metabolites to specific TD data records, thereby overruling the automated matching procedure. Another important configuration table, designated modify_drg0_bounds, will be added after model relaxation.
[3]:
# Data in thermo-db, which requires update
mod_thermo_data_cols = ['td_sid', 'attribute', 'value', 'notes']
mod_thermo_data = [
['cpd00637', 'charge_std', 1, 'met__D'],
['cpd15461', 'charge_std', 2, 'fe3hox_un'],
['cpd19833', 'charge_std', 0, '3amac'],
['cpd19835', 'charge_std', 1, 'athtp'],
['cpd19836', 'charge_std', 0, 'cdg'],
['cpd19837', 'charge_std', 1, 'chtbs6p'],
['cpd19838', 'charge_std', 0, 'cph4'],
['cpd19839', 'charge_std', 1, 'ddcap'],
['cpd19844', 'charge_std', 1, 'hdcap'],
['cpd19845', 'charge_std', 1, 'hdceap'],
['cpd19854', 'charge_std', 1, 'ocdcap'],
['cpd19855', 'charge_std', 1, 'ocdceap'],
['cpd19861', 'charge_std', 1, 'rephaccoa'],
['cpd19863', 'charge_std', 1, 'ttdcap'],
['cpd19864', 'charge_std', 1, 'ttdceap'],
['cpd19848', 'charge_std', 1, 'malcoame']
]
# Model annotations can contain several seed ids, here we hard link selected metabolites to a specific seed id
modify_td_sids_data = [['h2o', 'cpd00001'], ['lcts', 'cpd00208'], ['feoxam', 'cpd04761'],
['ddcap', 'cpd19839'], ['pheme', 'cpd00028'], ['malcoame', 'cpd19848']]
modify_td_sids_cols = ['mid', 'td_sid']
[4]:
# preliminary TFA configuration file
tfa_params = {}
data = [['thermo_data_fname', os.path.join('data', 'thermo_data.thermodb')],
['mid_regex_pattern', r'^M_(\w+)_\w+$']]
tfa_params['general'] = pd.DataFrame(data, columns=['parameter', 'value']).set_index('parameter')
cols = ['cid', 'ph', 'ionic_strength_M', 'c_min_M', 'c_max_M', 'c_mV', 'p_mV', 'e_mV']
data = [['c', 7.5, 0.25, 1.0e-8, 0.05, 0.0, 150.0, 150.0],
['p', 7.0, 0.0, 1.0e-8, 0.05, -150.0, 0.0, 0.0],
['e', 7.0, 0.0, 1.0e-8, 0.10, -150.0, 0.0, 0.0]]
tfa_params['td_compartments'] = pd.DataFrame(data, columns=cols).set_index('cid')
# try models without 'modify_td_sids' and/or without 'modify_thermo_data'
tfa_params['modify_td_sids'] = pd.DataFrame(modify_td_sids_data, columns=modify_td_sids_cols).set_index('mid')
tfa_params['modify_thermo_data'] = pd.DataFrame(mod_thermo_data, columns=mod_thermo_data_cols).set_index('td_sid')
write_parameter_file(os.path.join('data', 'preliminary_tfa_parameters.xlsx'), tfa_params)
4 table(s) with parameters written to data/preliminary_tfa_parameters.xlsx
Step 3: Create preliminary TFA model¶
The preliminary TFA model will be utilized for constraints relaxation. The creation of the TFA model is analogous to the creation of other extended models, except that we do not provide additional configuration data to the Xba model.
[5]:
# Create preliminary TFA model
xba_model = XbaModel(os.path.join('SBML_models', f'{fba_model}.xml'))
xba_model.configure()
tfa_model = TfaModel(xba_model)
tfa_model.configure(os.path.join('data', 'preliminary_tfa_parameters.xlsx'))
loading: SBML_models/iML1515.xml (last modified: Thu Dec 5 10:03:46 2024)
1877 constraints (+0); 2712 variables (+0); 1516 genes (+0); 5 parameters (+0)
>>> BASELINE XBA model configured!
4 table(s) with parameters loaded from data/preliminary_tfa_parameters.xlsx (Fri Nov 21 11:15:26 2025)
16 thermo data attributes updated
938 metabolites with TD data covering 1546 model species; (331 not supported)
1899 reactions supported by TD data; (476 not supported)
7221 constraints to add
7221 constraint ids added to the model (9098 total constraints)
3668 ∆rG'/∆rG'˚ variables to add
3668 variable ids added to the model (6380 total variables)
2407 forward/reverse use variables to add
2407 variable ids added to the model (8787 total variables)
1525 log concentration variables to add
1525 variable ids added to the model (10312 total variables)
1834 TD reactions split in forward/reverse, 0 opened reverse direction
2407 fwd/rev reactions to couple with flux direction
2407 attributes on reaction instances updated
2 RHS variables to add
2 variable ids added to the model (10887 total variables)
3682 parameters
1525 species (929 metabolites) with TD data from total 1877 (1169)
1834 metabolic/transporter reactions with TD data from total 2375
9098 constraints (+7221); 10887 variables (+8175); 1516 genes (+0); 3682 parameters (+3677)
[5]:
True
Step 4: Model relaxation (cobrapy)¶
The optimization of our preliminary TFA model may result in infeasible solutions or solutions with significantly reduced solution space. By relaxing constraints, namely the bounds of selected standard-transformed Gibbs energy of reaction variables, we can recover feasible optimization solutions for our selected media conditions. As a reference, we use the predicted growth rates for the original FBA model. Though not required, we perform the relaxation considering stiffer bounds on selected metabolites.
The model relaxation workflow initiates with the construction of a slack model, which is subsequently loaded with cobrapy. The FBA growth rates are then determined for utilization as a reference. Consequently, the slack model undergoes optimization across the designated conditions, facilitating the identification of variables necessitating relaxation.
4.1 create slack model¶
The utilization of the function tfa_model.export_slack_model() facilitates the generation of a slack model, incorporating positive and negative slack variables within the constraints that determine the transformed Gibbs free energies or reactions. The optimization objective of the slack model is modified to minimize the sum of slack.
[6]:
# Create slack model
slack_fname = os.path.join('SBML_models', f'{target_model}_slack.xml')
tfa_model.export_slack_model(slack_fname)
3668 slack variables for ∆rG'˙ to add
3668 variable ids added to the model (14555 total variables)
model exported to SBML format: SBML_models/iML1515_TFA_slack.xml
4.2 loading of slack model¶
The slack is loaded with cobrapy, and the solver parameter ‘presolve’ is activated for model relaxation. Selected optimization variables are reconfigured from ‘continuous’ to ‘binary’ by FbaOptimization, and the sign of selected constraints is changed from ‘equality’ to ‘inequality’. The value ranges of selected nucleotide concentrations in the model are constrained by fo.set_tfa_metab_concentrations(), which is optional.
[7]:
# Load slack model (cobrapy)
tfa_slack_model = cobra.io.read_sbml_model(slack_fname)
print(f'TD Slack model loaded from {slack_fname}, last modified: {time.ctime(os.path.getmtime(slack_fname))}')
tfa_slack_model.solver.configuration.presolve = True
# Reconfigure variables and constraints
fo = FbaOptimization(slack_fname, tfa_slack_model)
# Set nucleotide concentrations (optional)
for sid, (lb, ub) in metabolite_molar_concs.items():
tfa_slack_model.reactions.get_by_id('V_LC_' + sid).bounds = (np.log(lb), np.log(ub))
print(f'{len(metabolite_molar_concs)} metabolite concentrations constrained')
Set parameter Username
Set parameter LicenseID to value 2731209
Academic license - for non-commercial use only - expires 2026-10-31
TD Slack model loaded from SBML_models/iML1515_TFA_slack.xml, last modified: Fri Nov 21 11:15:58 2025
SBML model loaded by sbmlxdf: SBML_models/iML1515_TFA_slack.xml (Fri Nov 21 11:15:58 2025)
Thermodynamic use variables (V_FU_xxx and V_RU_xxx) as binary
Thermodynamic constraints (C_F[FR]C_xxx, C_G[FR]C_xxx, C_SU_xxx) ≤ 0
1834 TD reaction constraints
13 metabolite concentrations constrained
4.3 FBA growth rates¶
We employ cobrapy to predict the growth rates of the FBA model under our specific media conditions. These predictions subsequently serve as targets for the slack model.
[8]:
# Determine FBA growth rates (using cobrapy)
fba_fname = os.path.join('SBML_models', f'{fba_model}.xml')
gem_model = cobra.io.read_sbml_model(fba_fname)
gem_model.objective = biomass_rid
fba_grs = {}
for condition, medium in conditions.items():
gem_model.medium = medium
fba_grs[condition] = gem_model.slim_optimize()
print(f'{len(fba_grs)} growth conditions with FBA growth rates in range '
f'[{min(fba_grs.values()):.2f} - {max(fba_grs.values()):.2f}] h-1.')
6 growth conditions with FBA growth rates in range [0.21 - 0.91] h-1.
4.4 relax variable bounds¶
In this process, instances of standard-transformed Gibbs energy of reaction variables where bounds require modification are identified. The slack model is optimized under specific media conditions, while the lower bound of the biomass reactions is configured at a level of 95% of the respective predicted growth rate of the FBA model. The objective function of the slack model is to minimize the total slack. Any slack encountered is used to modify the bound of the respective standard-transformed
Gibbs energy of reaction variable using fo.update_relaxation(). The modifications are collected in drg0_relaxations and subsequently used to modify the TFA configuration file.
[9]:
# relax drg0 bounds across several conditions
start = time.time()
pct_fba_gr = 0.95
drg0_relaxations = {}
for condition, medium in conditions.items():
tfa_slack_model.medium = medium
tfa_slack_model.reactions.get_by_id(biomass_rid).lower_bound = pct_fba_gr * fba_grs[condition]
if np.isfinite(tfa_slack_model.slim_optimize()) is False:
print(f'{condition} no solution during slack minimization')
else:
solution = tfa_slack_model.optimize()
gr = solution.fluxes[biomass_rid]
req_relaxations = fo.update_relaxation(solution.fluxes)
drg0_relaxations.update(req_relaxations)
print(f'{condition:20s} ({gr:5.3f} h-1): slack of {solution.objective_value:7.2f}'
f"({solution.status}) -> {len(req_relaxations):2d} ∆rG'˚ updates")
tfa_slack_model.reactions.get_by_id(biomass_rid).lower_bound = 0.0
print(f'{len(drg0_relaxations)} dgr0 relaxations required - update ')
print(f'Duration: {time.time()-start:.1f} s')
print(f'{len(drg0_relaxations)} variables that need relaxation:')
print(drg0_relaxations)
Acetate (0.200 h-1): slack of 167.18(optimal) -> 3 ∆rG'˚ updates
Glycerol (0.470 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
Fructose (0.833 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
L-Malate (0.391 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
Glucose (0.833 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
Glucose 6-Phosphate (0.860 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
3 dgr0 relaxations required - update
Duration: 6.5 s
3 variables that need relaxation:
{'V_DRG0_MECDPS': {'fbc_lower_bound': 83.39185422440065}, 'V_DRG0_ATPS4rpp': {'fbc_lower_bound': -29.10627304266519}, 'V_DRG0_MALCOAMT': {'fbc_lower_bound': 75.96556319818774}}
Step 5: Update TFA parameters and create final TFA model¶
The TFA configuration file is extended with the table designated as modify_drg0_bounds, which encompasses the parameters determined during the process of model relaxation. For the purpose of model inspection, the option of exporting the model to spreadsheet format is available.
[10]:
# update TFA configuration file
tfa_params = load_parameter_file(os.path.join('data', 'preliminary_tfa_parameters.xlsx'))
data = []
for var_id, bounds in drg0_relaxations.items():
for bound_id, val in bounds.items():
data.append([var_id, 'reaction', bound_id, val])
cols = ['id', 'component', 'attribute', 'value']
tfa_params['modify_drg0_bounds'] = pd.DataFrame(data, columns=cols).set_index('id')
write_parameter_file(os.path.join('data', f'{target_model}_tfa_parameters.xlsx'), tfa_params)
# create TFA model
xba_model = XbaModel(os.path.join('SBML_models', f'{fba_model}.xml'))
xba_model.configure()
tfa_model = TfaModel(xba_model)
tfa_model.configure(os.path.join('data', f'{target_model}_tfa_parameters.xlsx'))
tfa_model.export(os.path.join('SBML_models', f'{target_model}.xml'))
tfa_model.export(os.path.join('xlsx_models', f'{target_model}.xlsx')) # optional export to spreadsheet format
4 table(s) with parameters loaded from data/preliminary_tfa_parameters.xlsx (Fri Nov 21 11:15:26 2025)
5 table(s) with parameters written to data/iML1515_TFA_tfa_parameters.xlsx
loading: SBML_models/iML1515.xml (last modified: Thu Dec 5 10:03:46 2024)
1877 constraints (+0); 2712 variables (+0); 1516 genes (+0); 5 parameters (+0)
>>> BASELINE XBA model configured!
5 table(s) with parameters loaded from data/iML1515_TFA_tfa_parameters.xlsx (Fri Nov 21 11:16:34 2025)
16 thermo data attributes updated
938 metabolites with TD data covering 1546 model species; (331 not supported)
1899 reactions supported by TD data; (476 not supported)
7221 constraints to add
7221 constraint ids added to the model (9098 total constraints)
3668 ∆rG'/∆rG'˚ variables to add
3668 variable ids added to the model (6380 total variables)
2407 forward/reverse use variables to add
2407 variable ids added to the model (8787 total variables)
1525 log concentration variables to add
1525 variable ids added to the model (10312 total variables)
1834 TD reactions split in forward/reverse, 0 opened reverse direction
2407 fwd/rev reactions to couple with flux direction
2407 attributes on reaction instances updated
2 RHS variables to add
2 variable ids added to the model (10887 total variables)
3682 parameters
3 ∆Gr'˚ variables need relaxation.
3 attributes on parameter instances updated
1525 species (929 metabolites) with TD data from total 1877 (1169)
1834 metabolic/transporter reactions with TD data from total 2375
9098 constraints (+7221); 10887 variables (+8175); 1516 genes (+0); 3682 parameters (+3677)
model exported to SBML format: SBML_models/iML1515_TFA.xml
model exported to Excel Spreadsheet: xlsx_models/iML1515_TFA.xlsx
[10]:
True
Step 6. Load and optimize TFA model (cobrapy)¶
In this section, we will demonstrate the utilization of cobrapy for the purpose of TFA model loading and optimization. In the following section, we will show the optimization using the gurobipy interface. The optimization problem is a mixed-integer linear optimization (MILP) problem. While MILP optimization is supported by cobrapy, some variables and constraints need to be modified. This is done automatically by FbaOptimization. As an alternative option, in a non-Python environment, these
variables and constraints could be executed separately after model loading. The following Python example code is provided:
tfam = cobra.io.read_sbml_model(os.path.join('SBML_models', f'{target_model}.xml'))
for var in tfam.variables:
for vprefix in ['V_FU_', 'V_RU_']:
if re.match(vprefix, var.name) and 'reverse' not in var.name:
var.type = 'binary'
for constr in tfam.constraints:
for cprefix in ['C_FFC_', 'C_FRC_', 'C_GFC_', 'C_GRC_', 'C_SU_']:
if re.match(cprefix, constr.name):
constr.lb = None
Constraining the concentrations of nucleotides is an optional process that can be carried out by utilizing the function fo.set_tfa_metab_concentrations(). In instances where the ‘FbaOptimization’ class is unavailable, the concentrations can be set using the Python code provided below.
import numpy as np
for sid, (lb, ub) in metabolite_molar_concs.items():
tfam.reactions.get_by_id('V_LC_' + sid).bounds = (np.log(lb), np.log(ub))
[11]:
# Load model using cobrapy
fname = os.path.join('SBML_models', f'{target_model}.xml')
tfam = cobra.io.read_sbml_model(fname)
for sid, (lb, ub) in metabolite_molar_concs.items():
tfam.reactions.get_by_id('V_LC_' + sid).bounds = (np.log(lb), np.log(ub))
# Reconfigure variables and constraints
fo = FbaOptimization(fname, tfam)
# Set nucleotide concentrations (optional)
for sid, (lb, ub) in metabolite_molar_concs.items():
tfam.reactions.get_by_id('V_LC_' + sid).bounds = (np.log(lb), np.log(ub))
print(f'{len(metabolite_molar_concs)} metabolite concentrations constrained')
SBML model loaded by sbmlxdf: SBML_models/iML1515_TFA.xml (Fri Nov 21 11:16:59 2025)
Thermodynamic use variables (V_FU_xxx and V_RU_xxx) as binary
Thermodynamic constraints (C_F[FR]C_xxx, C_G[FR]C_xxx, C_SU_xxx) ≤ 0
1834 TD reaction constraints
13 metabolite concentrations constrained
[12]:
# Optimize model using cobrapy and analyze results
pred_results = {}
for cond, medium in conditions.items():
with tfam as model:
model.medium = medium
solution = model.optimize()
if solution.status == 'optimal':
gr = solution.objective_value
pred_results[cond] = solution
print(f'{cond:25s}: TFA gr: {gr:.3f} h-1 vs. FBA {fba_grs[cond]:.3f}, '
f'diff: {gr - fba_grs[cond]:6.3f}')
else:
print(f'{cond} ended with status {solution.status}')
fr = FbaResults(fo, pred_results)
df_fluxes = fr.collect_fluxes()
df_net_fluxes = fr.collect_fluxes(net=True)
df_species_conc = fr.collect_species_conc()
fr.save_to_escher(df_net_fluxes[reference_cond], os.path.join('escher', target_model))
fr.save_to_escher(df_species_conc[reference_cond], os.path.join('escher', target_model))
Acetate : TFA gr: 0.210 h-1 vs. FBA 0.210, diff: -0.000
Glycerol : TFA gr: 0.495 h-1 vs. FBA 0.495, diff: -0.000
Fructose : TFA gr: 0.877 h-1 vs. FBA 0.877, diff: -0.000
L-Malate : TFA gr: 0.411 h-1 vs. FBA 0.411, diff: -0.000
Glucose : TFA gr: 0.877 h-1 vs. FBA 0.877, diff: -0.000
Glucose 6-Phosphate : TFA gr: 0.905 h-1 vs. FBA 0.905, diff: -0.000
1 file(s) exported for "Load reaction data" into Escher maps
1 file(s) exported for "Load metabolite data" into Escher maps
species concentrations¶
The function fr.collect_species_conc() enables the retrieval of a table that contains metabolite concentrations.
[13]:
print(f'{len(df_species_conc)} species concentrations')
df_species_conc.loc[list(metabolite_molar_concs)].head()
1525 species concentrations
[13]:
| name | rank | mean mmol_per_l | stdev | Acetate | Glycerol | Fructose | L-Malate | Glucose | Glucose 6-Phosphate | |
|---|---|---|---|---|---|---|---|---|---|---|
| sid | ||||||||||
| atp_c | ATP C10H12N5O13P3 | 1046 | 1.780000 | 0.000000e+00 | 1.7800 | 1.7800 | 1.7800 | 1.7800 | 1.780000 | 1.780000 |
| adp_c | ADP C10H12N5O10P2 | 1130 | 0.203000 | 7.103520e-02 | 0.2320 | 0.2320 | 0.0580 | 0.2320 | 0.232000 | 0.232000 |
| datp_c | DATP C10H12N5O12P3 | 1154 | 0.090500 | 1.520235e-17 | 0.0905 | 0.0905 | 0.0905 | 0.0905 | 0.090500 | 0.090500 |
| ctp_c | CTP C9H12N3O14P3 | 1069 | 0.645782 | 6.534249e-03 | 0.6500 | 0.6500 | 0.6500 | 0.6500 | 0.637346 | 0.637346 |
| dctp_c | DCTP C9H12N3O13P3 | 1153 | 0.092000 | 1.520235e-17 | 0.0920 | 0.0920 | 0.0920 | 0.0920 | 0.092000 | 0.092000 |
Closing remarks¶
The generation of a thermodynamics constraint model of iML1515 is a straightforward process requiring minimal effort in configuration. A preliminary TFA model was used to relax some of the thermodynamics constraints, and the final TFA model predicts similar growth rates as the FBA model. It also determines a consistent set of metabolite concentrations that ensure that flux directions align with negative values of transformed Gibbs energies of reaction.
In the subsequent tutorials, we will integrate the TFA-specific configuration with ECM, and the RBA-specific configuration data, to generate thermodynamics constraint TGECKO and TRBA models of iML1515.
(Alternative) gurobipy – TFA model creation¶
The workflow of generating a TFA model using the gurobipy interface involves similar steps as presented above; however, the workflow deviates with the loading of the slack model using FbaOptimization. In this section, we will demonstrate the use of gurobipy to predict the growth rates of the FBA model.During parameter relaxation of the slack model, the minimal expected growth rates are set using fo.set_variable_bounds(). The update for the TFA configuration file and the creation of the
final TFA model remain unchanged.
[14]:
# Load slack model (gurobipy)
fo = FbaOptimization(slack_fname)
# Set nucleotide concentrations (optional)
orig_concs = fo.set_tfa_metab_concentrations(metabolite_molar_concs)
print(f'{len(metabolite_molar_concs)} metabolite concentrations constrained')
SBML model loaded by sbmlxdf: SBML_models/iML1515_TFA_slack.xml (Fri Nov 21 11:15:58 2025)
MILP Model of iML1515_TFA
14555 variables, 9098 constraints, 39194 non-zero matrix coefficients
1834 TD reaction constraints
13 metabolite concentrations constrained
[15]:
# Determine FBA growth rates (gurobipy)
fba_fo = FbaOptimization(os.path.join('SBML_models', f'{fba_model}.xml'))
fba_fo.set_objective({biomass_rid: 1.0}, 'max')
fba_grs = {}
for condition, medium in conditions.items():
fba_fo.medium = medium
fba_grs[condition] = fba_fo.optimize().objective_value
print(f'{len(fba_grs)} growth conditions with FBA growth rates in range '
f'[{min(fba_grs.values()):.2f} - {max(fba_grs.values()):.2f}] h-1.')
SBML model loaded by sbmlxdf: SBML_models/iML1515.xml (Thu Dec 5 10:03:46 2024)
LP Model of iML1515
2712 variables, 1877 constraints, 10575 non-zero matrix coefficients
6 growth conditions with FBA growth rates in range [0.21 - 0.91] h-1.
[16]:
# relax drg0 bounds across several conditions (gurobipy)
start = time.time()
pct_fba_gr = 0.95
drg0_relaxations = {}
for condition, medium in conditions.items():
fo.medium = medium
orig_bounds = fo.set_variable_bounds({biomass_rid:(pct_fba_gr * fba_grs[condition], None)})
solution = fo.optimize()
if solution.status == 'optimal':
gr = solution.fluxes[biomass_rid]
req_relaxations = fo.update_relaxation(solution.fluxes)
drg0_relaxations.update(req_relaxations)
print(f'{condition:20s} ({gr:5.3f} h-1): slack of {solution.objective_value:7.2f}'
f"({solution.status}) -> {len(req_relaxations):2d} ∆rG'˚ updates")
else:
print(f'{condition} no solution during slack minimization')
fo.set_variable_bounds({biomass_rid: (0.0, 1000.0)})
print(f'{len(drg0_relaxations)} dgr0 relaxations required - update ')
print(f'Duration: {time.time()-start:.1f} s')
print(f'{len(drg0_relaxations)} variables that need relaxation:')
print(drg0_relaxations)
Acetate (0.200 h-1): slack of 167.18(optimal) -> 4 ∆rG'˚ updates
Glycerol (0.470 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
Fructose (0.833 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
L-Malate (0.391 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
Glucose (0.833 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
Glucose 6-Phosphate (0.860 h-1): slack of 0.00(optimal) -> 0 ∆rG'˚ updates
4 dgr0 relaxations required - update
Duration: 4.7 s
4 variables that need relaxation:
{'V_DRG0_MECDPS': {'fbc_lower_bound': 83.39185422440065}, 'V_DRG0_ATPS4rpp': {'fbc_lower_bound': -25.456120671846797}, 'V_DRG0_PRAIS': {'fbc_lower_bound': 62.67390227094072}, 'V_DRG0_MALCOAMT': {'fbc_lower_bound': 75.96556319818777}}
[17]:
# update TFA configuration file
tfa_params = load_parameter_file(os.path.join('data', 'preliminary_tfa_parameters.xlsx'))
data = []
for var_id, bounds in drg0_relaxations.items():
for bound_id, val in bounds.items():
data.append([var_id, 'reaction', bound_id, val])
cols = ['id', 'component', 'attribute', 'value']
tfa_params['modify_drg0_bounds'] = pd.DataFrame(data, columns=cols).set_index('id')
write_parameter_file(os.path.join('data', f'{target_model}_tfa_parameters.xlsx'), tfa_params)
# create TFA model
xba_model = XbaModel(os.path.join('SBML_models', f'{fba_model}.xml'))
xba_model.configure()
tfa_model = TfaModel(xba_model)
tfa_model.configure(os.path.join('data', f'{target_model}_tfa_parameters.xlsx'))
tfa_model.export(os.path.join('SBML_models', f'{target_model}.xml'))
4 table(s) with parameters loaded from data/preliminary_tfa_parameters.xlsx (Fri Apr 11 15:06:56 2025)
5 table(s) with parameters written to data/iML1515_TFA_tfa_parameters.xlsx
loading: SBML_models/iML1515.xml (last modified: Thu Dec 5 10:03:46 2024)
1877 constraints (+0); 2712 variables (+0); 1516 genes (+0); 5 parameters (+0)
>>> BASELINE XBA model configured!
5 table(s) with parameters loaded from data/iML1515_TFA_tfa_parameters.xlsx (Fri Apr 11 15:09:45 2025)
16 thermo data attributes updated
938 metabolites with TD data covering 1546 model species; (331 not supported)
1899 reactions supported by TD data; (476 not supported)
7221 constraints to add
7221 constraint ids added to the model (9098 total constraints)
3668 ∆rG'/∆rG'˚ variables to add
3668 variable ids added to the model (6380 total variables)
2407 forward/reverse use variables to add
2407 variable ids added to the model (8787 total variables)
1525 log concentration variables to add
1525 variable ids added to the model (10312 total variables)
1834 TD reactions split in forward/reverse, 0 opened reverse direction
2407 fwd/rev reactions to couple with flux direction
2407 attributes on reaction instances updated
2 RHS variables to add
2 variable ids added to the model (10887 total variables)
3682 parameters
3 ∆Gr'˚ variables need relaxation.
3 attributes on parameter instances updated
1525 species (929 metabolites) with TD data from total 1877 (1169)
1834 metabolic/transporter reactions with TD data from total 2375
9098 constraints (+7221); 10887 variables (+8175); 1516 genes (+0); 3682 parameters (+3677)
model exported to SBML: SBML_models/iML1515_TFA.xml
[17]:
True
(Alternative) gurobipy – TFA model optimization¶
[17]:
# Load TFA model using gurobipy
fname = os.path.join('SBML_models', f'{target_model}.xml')
fo = FbaOptimization(fname)
# Set nucleotide concentrations (optional)
orig_concs = fo.set_tfa_metab_concentrations(metabolite_molar_concs)
print(f'{len(orig_concs)} metabolite concentrations constrained')
SBML model loaded by sbmlxdf: SBML_models/iML1515_TFA.xml (Fri Nov 21 11:16:59 2025)
MILP Model of iML1515_TFA
10887 variables, 9098 constraints, 35526 non-zero matrix coefficients
1834 TD reaction constraints
13 metabolite concentrations constrained
[18]:
# Optimize model using gurobipy and analyze results
pred_results = {}
for cond, medium in conditions.items():
fo.medium = medium
solution = fo.optimize()
if solution.status == 'optimal':
gr = solution.objective_value
pred_results[cond] = solution
print(f'{cond:25s}: TFA gr: {gr:.3f} h-1 vs. FBA {fba_grs[cond]:.3f}, '
f'diff: {gr - fba_grs[cond]:6.3f}')
else:
print(f'{cond} ended with status {solution.status}')
fr = FbaResults(fo, pred_results)
df_fluxes = fr.collect_fluxes()
df_net_fluxes = fr.collect_fluxes(net=True)
df_species_conc = fr.collect_species_conc()
fr.save_to_escher(df_net_fluxes[reference_cond], os.path.join('escher', target_model))
fr.save_to_escher(df_species_conc[reference_cond], os.path.join('escher', target_model))
Acetate : TFA gr: 0.210 h-1 vs. FBA 0.210, diff: -0.000
Glycerol : TFA gr: 0.495 h-1 vs. FBA 0.495, diff: -0.000
Fructose : TFA gr: 0.877 h-1 vs. FBA 0.877, diff: -0.000
L-Malate : TFA gr: 0.411 h-1 vs. FBA 0.411, diff: -0.000
Glucose : TFA gr: 0.877 h-1 vs. FBA 0.877, diff: -0.000
Glucose 6-Phosphate : TFA gr: 0.905 h-1 vs. FBA 0.905, diff: -0.000
1 file(s) exported for "Load reaction data" into Escher maps
1 file(s) exported for "Load metabolite data" into Escher maps
Implementation details¶
The ensuing sections delineate the variables and constraints that have been integrated into the genome-scale metabolic model during its conversion into a thermodynamics constraint model. These sections also elucidate the formulas employed to derive the numerical values.
Variables and Constraints¶
The following paragraphs detail the variables and constraints that have been incorporated into the genome-scale metabolic model. To illustrate this, the reaction of fructose-bisphosphate aldolase is examined, which is implemented as a reversible reaction in the iML1515 model with the identifier R_FBA: ‘fdp_c -> dhap_c + g3p_c’. During the extension of the model, the reaction is rendered irreversible, designated as R_FBA: ‘fdp_c => dhap_c + g3p_c’, and a new reaction catalyzing the reverse
direction is incorporated. This reverse reaction is designated as R_FBA_REV: ‘dhap_c + g3p_c => fdp_c’. It is noteworthy that the reverse reaction is not included for reactions that have been configured as irreversible in the original model.
Two additional variables are incorporated to represent the transformed Gibbs energy of reaction, designated as V_DRG_FBA and V_DRG0_FBA, respectively. The former is assigned unlimited bounds, while the latter is constrained to the calculated standard transformed Gibbs energy of reaction plus or minus the estimation error.
The log concentration variables, designated as V_LC_fdp_c, V_LC_dhap_c, and V_LC_g3p_c, are incorporated with default bounds as per compartmental configuration, and the concentrations are expressed in units of mol/L. The calculation formula for the transformed Gibbs energy of reaction is implemented by the constraint C_DRG_FBA: ‘V_DRG_FBA = V_DRG0_FBA - 2.48 V_LC_fdp_c + 2.48 V_LC_dhap_c + 2.48 V_LC_g3p_c’ (RT = 2.48 kJ/mol).
Two binary variables (values 0 or 1) designated “V_FU_FBA” (“forward use”) and “V_RU_FBA” (“reverse use”) are introduced to couple the transformed Gibbs energy of reaction to the flux direction. The implementation of a “simultaneous use” constraint, denoted as C_SU_FBA, ensures that only one of the use variables can take the value “1.” This is expressed as “V_FU_FBA + V_RU_FBA ≤ 1.” Two Gibbs energy coupling constraints, designated as C_GFC_FBA and C_GRC_FBA, couple the forward and
reverse use variables to the transformed Gibbs energy of reaction with inequalities ‘V_DRG_FBA ≤ 999.99 - 1000 V_FU_FBA’, thereby forcing V_DRG_FBA ≤ 0.01 kJ/mol when V_FU_FBA is active, and ‘V_DRG_FBA ≥ 1000 V_RU_FBA - 999.99’, forcing V_DRG_FBA ≥ 0.01 kJ/mol when V_RU_FBA is active. In a similar fashion, reactions fluxes in forward and reverse direction are coupled to the forward and reverse use variables via the flux coupling constraints C_FFC_FBA: R_FBA ≤ 1000 V_FU_FBA and C_FRC_FBA:
R_FBA_REV ≤ 1000 V_FB_FBA. This configuration can be readily verified by exporting the TFA model to the ‘.xlsx’ format.
Calculation Details¶
Thermodynamic constraints couple reaction flux directionality with Gibbs energy of reaction. The transformed Gibbs energy values employed in this context are derived through transformations with respect to compartmental pH and ionic strength at the default temperature of 298.15 K (25˚C). Negative values of the transformed Gibbs energy of reaction will drive the reaction in the forward direction, while positive values will drive it in the reverse direction. The incorporation of thermodynamic constraints into a model is a straightforward process that necessitates minimal configuration input. However, it should be noted that the underlying calculations are of a highly complex nature. The f2xba model is aligned with the formulation implemented in the pyTFA package (Salvy et al., 2019), with the formulas being based on the book by Alberty (Alberty, 2003).
The natural log of metabolite concentrations is a variable in the optimization problem, with lower and upper bounds defined in the TFA configuration file and potentially further limited prior to optimization. The factor of gas constant times temperature, ‘RT’, used in subsequent equations, is evaluated to 2.48 kJ/mol at T = 298.15 K.
The transformed Gibbs energy of reaction, denoted by \(\Delta_r \text G^{'}\), is calculated from the standard transformed Gibbs energy of reaction, denoted by \(\Delta_r \text G^{'0}\), the metabolite concentrations \(c_i\) [mol/L], and the stoichiometric coefficients \(\nu_i\) of reactants (negative) and products (positive) using the following equation 4.5-10 (Alberty, 2003):
\(\Delta_r \text G^{'} = \Delta_r \text G^{'0} + \text{RT} \sum_i {\nu_i \ln c_i}\)
The group contribution method involves the decomposition of a molecule into cues (groups) for which the Gibbs free energy of formation has been estimated with high confidence. During a reaction, only some of the cues of reactants and products, the net cues, undergo modification, while the majority of the other cues remain unmodified. The TD database contains estimated errors for each of the cues. The estimated errors of the net cues, denoted by \(\text{cue_est_error}_j\), are then utilized to ascertain the estimation error for the Gibbs energy of reaction. This estimation error is subsequently employed to establish the bounds of the variable \(\Delta_r \text G^{'0}\). The calculation is performed as per pyTFA:
\(\text{estimation_error} = \sqrt {\sum_j (\nu_j \, \text{cue_est_error}_j)^2}\)
The standard-transformed Gibbs energy of reaction, denoted by \(\Delta_r \text G^{'0}\), is derived from the standard-transformed Gibbs energies of formation, denoted by \(\Delta_f \text G_i^{'0}\), of the reactants and products. This derivation is expressed through the following equation 4.4-2 (Alberty, 2003):
\(\Delta_r \text G^{'0} = \sum_i \nu_i \Delta_f \text G_i^{'0}\)
In the context of a transport process, it is imperative to incorporate electrical work terms, which are calculated from the membrane potentials, denoted by the symbol \(\Delta \varphi_{sd}\) (destination minus source potential), and the transported charges, denoted by the symbol \(z_{sd}\) (source to destination compartment), which can assume positive or negative values. \(F\) is the Faraday constant (96.485 kJ mol-1 V-1). To derive the equation, we have consulted the work of Jol (Jol et al., 2010).
\(\Delta_r \text G^{'0} = \sum_i \nu_i \Delta_f \text G_i^{'0} + F \Delta \sum_{sd} \varphi_{ds} z_{sd}\)
The reaction network implemented by the genome-scale metabolic model consists of biochemical reactions and biochemical reactants (metabolites). At the molecular level, a biochemical reactant, such as ATP, can be regarded as a group (pseudoisomer group) of related chemical species in different protonation states and different complexations with metal ions, such as \(\text{ATP}^{4-}\), \(\text{HATP}^{3-}\), \(\text{H}_2\text{ATP}^{2-}\), \(\text{MgATP}^{2-}\) or \(\text{MgHATP}^-\). The f2xba modeling framework considers different protonation states, but not different complexations with metal ions.
The Gibbs energy of formation, denoted as \(\Delta_f \text G^{'0}(I)\), of a biochemical reactant can be determined from the Gibbs energy of formation of the least protonated chemical species, denoted as \(\Delta_f \text {G1}^{'0}(I)\), in the pseudoisomeric group and the contribution of the other chemical species in the pseudoisomeric group, which are considered by the binding polynomial, denoted as \(\text P(I)\). It is imperative to note that these values are contingent on the isomeric strength \(I\) [mol/L] of the compartment, as elucidated in equations 4.5-6 (Alberty, 2003).
\(\Delta_f \text G^{'0}(I) = \Delta_f \text{G1}^{'0}(I) - \text{RT} \ln(\text P(I))\)
The determination of the least protonated chemical species is contingent upon the compartmental pH, as defined in the TFA configuration file, in conjunction with the \(\text{pKa}_j\) values and the electrical charge extracted from the pertinent TD data record. The standard-transformed Gibbs energy of formation for the least protonated species, denoted as \(\Delta_f \text{G1}^{'0}\), is derived from the standard Gibbs energy of formation, denoted as \(\Delta_f \text G^{0}\), extracted from the corresponding TD database record. The latter is transformed to the compartmental \(pH\) using equation 4.10-12 (Alberty, 2003):
\(\Delta_f \text {G1}^{'0} = \Delta_f \text G^0 + RT \ln{(10)} \sum_j {\text{pKa}_j}\)
The standard transformed Gibbs energy of formation for this least protonated state, denoted as \(\Delta_f \text{G1}^{'0}(I)\) at a given ionic strength, is calculated from its value at zero ionic strength, denoted as \(\Delta_f \text{G1}^{'0}\), and adjustments with respect to the energy contribution by the number of protons, denoted as \(nH\), in the structure and electrical charge, denoted as \(z\). Ionic strength \(I\), defined as the measure of ion concentration, exerts a significant influence on the activity coefficients employed in equilibrium equations by means of shielding charges. This adjustment is achieved through the utilization of equation 4.4-10 (Alberty, 2003), which is founded on the extended Debye-Hueckel equation with constants \(A = 0.51065 \sqrt\frac{l}{mol}\) and \(B = 1.6\sqrt\frac{l}{mol}\).
\(\Delta_f \text{G1}^{'0}(I) = \Delta_f \text{G1}^{'0} + \text{RT} \cdot \ln(10) \cdot \text{nH} \cdot \text{pH} - \text{RT} \cdot \ln(10) \cdot (z^2 - \text{nH}) \cdot \frac{A \sqrt I}{1 + B \sqrt I}\)
Ionic strength-adjusted acid dissociation constants, denoted by the symbol \(\text{pKa}^{'}(I)\), are essential for determining the binding polynomial. These constants can be calculated from the corresponding \(\text{pKa}\) values using the following equation 4.10-11 (Alberty, 2003):
\(\text{pKa}^{'}(I) = \text{pKa} - \sum_j \nu_j z_j^2 \cdot \frac{A \sqrt I}{1 + B \sqrt I}\)
The binding polynomial, denoted by \(P(I)\), which accounts for the energy contribution of the other chemical species in the pseudoisomer group, is calculated from the equilibrium constants \(K_x\), which relate to the ionic strength-adjusted acid dissociation constants (\(K_x = 10^{-\text{pKa}_x^{'}}\)), with \(K_1\) having the smallest value (highest pKa) and the compartmental proton concentration \([H^+] = 10^{-pH}\), using equation 4.5-7 (Alberty, 2003).
\(P(I) = 1 + \frac {[H^+]} {K_1} + \frac {[H^+]^2} {K_1 \cdot K_2} + \dotsb\)
The mean number of bound hydrogens, denoted by \(\text{avg_h_binding}\), in the pseudoisomeric group can be calculated from the binding polynomial, as outlined in equation 1.3-9 (Alberty, 2003).
\(\text{avg_h_binding} = \frac {[H+]} {\text P} \frac {d \text P} {d[H+]} = \frac{[H+]/K_1 + 2[H+]^2/K_1K_2 + \dotsb } {\text P(I)}\)
References¶
Alberty, R. A. (2003). Thermodynamics of Biochemical Reactions. Massachusetts Institute of Technology Press, Cambridge, MA.
Buckstein Michael, H., He, J., & Rubin, H. (2008). Characterization of Nucleotide Pools as a Function of Physiological State in Escherichia coli. Journal of Bacteriology, 190(2), 718-726. https://doi.org/10.1128/jb.01020-07
Jankowski, M. D., Henry, C. S., Broadbelt, L. J., & Hatzimanikatis, V. (2008). Group contribution method for thermodynamic analysis of complex metabolic networks. Biophys J, 95(3), 1487-1499. https://doi.org/10.1529/biophysj.107.124784
Jol, S. J., Kümmel, A., Hatzimanikatis, V., Beard, D. A., & Heinemann, M. (2010). Thermodynamic calculations for biochemical transport and reaction processes in metabolic networks. Biophys J, 99(10), 3139-3144. https://doi.org/10.1016/j.bpj.2010.09.043
Mavrovouniotis, M. L. (1990). Group contributions for estimating standard gibbs energies of formation of biochemical compounds in aqueous solution. Biotechnol Bioeng, 36(10), 1070-1082. https://doi.org/10.1002/bit.260361013
Salvy, P., Fengos, G., Ataman, M., Pathier, T., Soh, K. C., & Hatzimanikatis, V. (2019). pyTFA and matTFA: a Python package and a Matlab toolbox for Thermodynamics-based Flux Analysis. Bioinformatics, 35(1), 167-169. https://doi.org/10.1093/bioinformatics/bty499
Seaver, S. M. D., Liu, F., Zhang, Q., Jeffryes, J., Faria, J. P., Edirisinghe, J. N., Mundy, M., Chia, N., Noor, E., Beber, M. E., Best, A. A., DeJongh, M., Kimbrel, J. A., D’Haeseleer, P., McCorkle, S. R., Bolton, J. R., Pearson, E., Canon, S., Wood-Charlson, E. M.,…Henry, C. S. (2021). The ModelSEED Biochemistry Database for the integration of metabolic annotations and the reconstruction, comparison and analysis of metabolic models for plants, fungi and microbes. Nucleic Acids Res, 49(D1), D575-d588. https://doi.org/10.1093/nar/gkaa746