1. Getting Started

multiTFA is a python package for accurate thermodynamic flux analysis. multiTFA is built on COBRApy to take advantage of COBRA model structure. Installation instructions can be found in README.rst.

1.1. Building a multiTFA model

To create a multitfa model, we require

1. A COBRA model

2. Metabolite database identifiers, such as BIGG, KEGG, SEED etc. - This is required to match and retrieve thermodyanamic information against our database. Metabolite identifiers can be added to the multiTFA model later. For simplicity, we use the same nomneclature as equilibrator-api.

3. Compartment specific pH and ionic strength information - To calculate Gibbs energy related to transmembrane transport, information about membrane potential is also required in the form of pandas Dataframes.

Optionally, intracellular metabolite concentrations can be provided to calculate the transformed Gibbs energy of reactions. If metabolite concentrations are not available, a broad metabolite concentration range is assumed.

NOTE: All the input information should be in S.I units, e.g: ionic strength, metabolite concentration in Moles.

To demonstrate the usage, we use E. Coli core model from BIGG database.

First step in creating a multitfa model is to get a COBRA model. For this example, we load E. coli core model. This is available in Examples directory.

from cobra.io import load_model

# Load the Cobra model
cobra_model = load_model('e_coli_core')

Let us define the compartment properties such as pH, ionic strength (M) etc. For transport calculations, we will provide membrane potential in mV. For convenience, we use pandas DataFrame to define these properties.

# compartment specific pH and Ionic strength information. Ionic strength units should be in M.
pH_I_dict = {
    "pH": {"c": 7.5, "e": 7, "p": 7},
    "I": {"c": 0.25, "e": 0, "p": 0},
}

# Membrane potentials between different compartments. Units for membrane potential is mV.
del_psi_dict = {
"c": {"c":0,"e": 0, "p": 150},
"e": {"c": 0,"e":0, "p": 0},
"p": {"c": -150, "e": 0, "p":0},
}

import pandas as pd
# Define membrane potential and compartment information in Dataframes
del_psi = pd.DataFrame.from_dict(data=del_psi_dict)
comp_info = pd.DataFrame.from_dict(data=pH_I_dict)

Now, it is time to list any reactions that needs to be excluded from thermodynamic analysis. For example, Biomass, demand reactions etc.

# List of reaction ids that you would like to exclude from thermodynamic analysis, here we exclude demand reactions
Excl = [
    rxn.id
    for rxn in cobra_model.reactions
    if rxn.id.startswith("EX_") or rxn.id.startswith("DM_")
] + ["BIOMASS_Ecoli_core_w_GAM", "O2t", "H2Ot"]

Now, we have all the information we need. We can now create an instance of multitfa model.

Please note, when you use multitfa for the very first time, it downloads datafiles of about 3 GB. This is a one time thing only.

# Now build the multiTFA model from COBRA model
from multitfa.core import tmodel
tfa_model = tmodel(
    cobra_model, Exclude_list=Excl, compartment_info=comp_info, membrane_potential=del_psi
)

Downloading package metadata...
Fragments already downloaded
Downloading package metadata...
Fragments already downloaded
Downloading package metadata...
Fragments already downloaded

In order to access the thermodynamic information from equilibrator api, we need to add a external database identifier against all the metabolites in tfa model. This can be done by populating the ‘Kegg_id’ property of the metabolites. If you do not know the identifier for some metabolites, it is okay, you can simply add NA against it. For simplicity, we use the notation similar to equilibrator api. For the example below in the E. coli core model, we use bigg identifiers.

You can follow equilibrator api notation for other database identifiers, for example, if you are using KEGG identifier, simply use kegg:kegg identifier.

for met in tfa_model.metabolites:
    kegg_id = 'bigg.metabolite:'+met.id[:-2]
    met.Kegg_id = kegg_id

tfa_model.update()

Now that we have setup the model and populated the thermodynamic constraints, its better to check if constraints are added properly.

print("Number of constraints before update {}\n".format(len(cobra_model.constraints)))
print("Number of reactions in the model {}\n".format(len(tfa_model.reactions)))
print("Number of reactions excluded from thermodynamic analysis {}\n".format(len(tfa_model.Exclude_reactions)))
print("Number of mass balance constraints {}\n".format(len(tfa_model.metabolites)))
num_cons = 2*3*(len(tfa_model.reactions) - len(tfa_model.Exclude_reactions)) + len(tfa_model.metabolites)
print("Number of constraints to be present after update {}\n".format(num_cons))
print("Number of constraints present after update {}\n".format(len(tfa_model.constraints)))

Number of constraints before update 72

Number of reactions in the model 95

Number of reactions excluded from thermodynamic analysis 23

Number of mass balance constraints 72

Number of constraints to be present after update 504

Number of constraints present after update 504

1.2 Metabolites

Reactions and metabolites in multiTFA model can be accessed similar to a COBRA model. For demonstration, we will consider cytosolic ATP as our metabolite with atp_c identifier in t_model. We will access ATP and print it’s standard Gibbs energy of formation and Also, change its intracellular concentration lower and upper bounds to ‘2e-3’ and ‘5e-2’ M respectively from it’s default bounds. These changes will reflect in solver variables and constraints.

atp = tfa_model.metabolites.get_by_id('atp_c')
print("Standard Gibbs energy of formation of ATP is {} KJ/mol".format(atp.delG_f))
print("initial concentrations")
print('{}\t{}\t{}'.format(atp.id, atp.concentration_min, atp.concentration_max))
atp.concentration_min = 2e-3
atp.concentration_max = 5e-2
print("Concentrations after modifying")
print('{}\t{}\t{}'.format(atp.id, atp.concentration_min, atp.concentration_max))

Standard Gibbs energy of formation of ATP is -2259.1882733696866 KJ/mol
initial concentrations
atp_c   1e-05   0.02

Selected objective sense:  MINIMIZE
Selected objective  name:  f3d52b45-5dba-11eb-bac2-c1c11d7a5b17
Selected RHS        name:  rhs
Selected bound      name:  bnd
Concentrations after modifying
atp_c   0.002   0.05

You can always change the Gibbs energy of formation value of a metabolite by resetting the delG_f property of the metabolite. Also, you can access additional information like, major pseudoisomer at given compartment conditions (major_ms), equilibrator accession of the metabolite.

1.3 Reactions

Reactions in multiTFA model are of type cobra.DictList. Similar to COBRA model, reactions can be accessed/retrieved by its index or by it’s name and properties/attributes can be accessed. Let us consider ATP synthase raction in our model by ‘ATPS4rpp’ and access it’s attributes.

atps = tfa_model.reactions.get_by_id('ATPS4r')
print(atps.reaction)

adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c

Access and change bounds of atps

print("lower bound is {} and upper bound is {}\n".format(atps.lower_bound, atps.upper_bound))
atps.lower_bound = 0
print("updated lower and upper bounds are {} \t {}\n".format(atps.lower_bound, atps.upper_bound))

lower bound is -1000.0 and upper bound is 1000.0

updated lower and upper bounds are 0     1000.0

We can access the Gibbs energy of transport and transformed Gibbs energy of the reaction

print("Gibbs energy of transport is {} KJ/mol".format(atps.delG_transport))
print("Transform Gibbs energy of is {} KJ/mol".format(atps.delG_prime))

Gibbs energy of transport is -4.605170185988092 KJ/mol
Transform Gibbs energy of is 32.39600409402766 KJ/mol

You can always reset the transport and transformed Gibbs energy of reactions. You can also access the stoichiometric matrix of the reaction, all the thermodynamic constraints associated with the reaction. Please refer to api for all the attributes of the model.