Case study iMM904 curation

Using the iMM904 model as a case study, we have used Thermo-Flux and did the few manual curation steps to show how models can be converted into a combined thermodynamic-stoichiometric model. Here, we present the steps that required manual attention when following the protocol described in the publication. “

[ ]:

import cobra
from cobra import Model, Reaction, Metabolite
from thermo_flux.core.metabolite import ThermoMetabolite
import thermo_flux
from thermo_flux.io import load_excel as ex
from thermo_flux.core.model import ThermoModel, ThermoReaction
from equilibrator_api import  Q_
import pandas as pd
from thermo_flux.io import load_gams as gs
from thermo_flux.io import helper_load as hl
from thermo_flux.utils.vis import compare_met
import numpy as np
import gurobipy as  gp
from gurobipy import GRB


import seaborn as sns
import matplotlib.pyplot as plt

Import sbml model from BiGG

[2]:

from cobra.io import load_json_model
model=load_json_model('iMM904.json')

#set objective to biomass
model.objective = model.reactions.BIOMASS_SC5_notrace
model.reactions.EX_co2_e.id='co2_EX'
model.reactions.EX_o2_e.id='o2_EX'
model.reactions.EX_etoh_e.id='etoh_EX'
model.repair()

for rxn in model.reactions:
    if 'EX_' not in rxn.id:
        rxn.lower_bound=-500 if rxn.lower_bound==-999999.0 else rxn.lower_bound
        rxn.upper_bound=500 if rxn.upper_bound==999999.0 else rxn.upper_bound
    else:
        rxn.lower_bound=-500 if rxn.lower_bound==-999999.0 else rxn.lower_bound
        rxn.upper_bound=500 if rxn.upper_bound==999999.0 else rxn.upper_bound
       # print(rxn.id,rxn.lower_bound,rxn.upper_bound)

#reset O2 exchange bound
model.optimize().to_frame().loc['BIOMASS_SC5_notrace']['fluxes']
model.reactions.o2_EX.bounds=(-500,0)

print(model.optimize().to_frame().loc['BIOMASS_SC5_notrace']['fluxes'])
model.optimize().to_frame().loc[[rxn.id for rxn in model.boundary]]

Set parameter Username
Academic license - for non-commercial use only - expires 2025-02-26
0.9748814206733495

[2]:

	fluxes	reduced_costs
EX_epistest_SC_e	0.0	-96.275549
EX_epist_e	0.0	-1.676719
EX_ergst_e	0.0	-1.722387
EX_ergstest_SC_e	0.0	-96.321217
EX_13BDglcn_e	0.0	-0.196130
...	...	...
EX_xyl__D_e	0.0	-0.163442
EX_xylt_e	0.0	-0.176421
EX_rib__D_e	0.0	-0.163442
EX_zymst_e	0.0	-1.620476
EX_zymstest_SC_e	0.0	-96.219306

164 rows × 2 columns

[3]:

#plot biomass yield for stoichiometric model
#GUR_range=[-0.25, -0.57, -0.6 , -1.1 , -1.67, -2.14, -2.22, -2.83, -3.24, -3.7, -4.67, -5.44, -7.62, -8.09, -8.33, -10.23, -13.18, -15.7, -22.42]
GUR_range = np.linspace(-0, -23, 50)

GUR_plot_range=[-x for x in GUR_range]
biomass_flux=[]
etoh_flux=[]
o2_flux=[]
co2_flux=[]
for GUR in GUR_range:
    model.reactions.EX_glc__D_e.upper_bound=0
    model.reactions.EX_glc__D_e.lower_bound=GUR
    sol = model.optimize()
    #add values to list
    biomass_flux.append(sol.fluxes.loc['BIOMASS_SC5_notrace'])
    etoh_flux.append(sol.fluxes.loc['etoh_EX'])
    o2_flux.append(sol.fluxes.loc['o2_EX'])
    co2_flux.append(sol.fluxes.loc['co2_EX'])
#plot biomass vs GUR in a scatter plot
sns.scatterplot(x=GUR_plot_range,y=biomass_flux)
plt.xlabel('Glucose uptake rate (mmol/(gDW.h))')
plt.ylabel('Biomass production (mmol/(gDW.h))')
fba_flux=pd.DataFrame({'biomass_EX':biomass_flux,'etoh_EX':etoh_flux,'o2_EX':o2_flux,'co2_EX':co2_flux},index=GUR_range)
# plt.savefig('Figures/growth_stoichio_only.png',dpi=300)
# model.optimize().to_frame().loc['BIOMASS_SC5_notrace']/model.optimize().to_frame().loc['EX_glc__D_e']

c:\Users\tedns\anaconda3\envs\thermoflux\lib\site-packages\cobra\util\solver.py:554: UserWarning: Solver status is 'infeasible'.
  warn(f"Solver status is '{status}'.", UserWarning)

Step 1 - Physical and Biochemical parameters

[ ]:

tmodel=ThermoModel(model,split_biomass=True, add_charge_exchange= True)

Initializing component contribution object...
No valid license for cxcalc installed, operating in read-only mode. A local cache may be loaded, but no compounds can be created. Please obtain a ChemAxon license to enable compound creation.
Loading compounds from iMM904_compound.sqlite
added reaction:  biomass_ce: biomass_c <=> biomass_e
added reaction:  biomass_EX: biomass_e <=>
added reaction:  charge_ce: charge_e <=> charge_c
added reaction:  EX_charge: charge_e <=>
added reaction:  charge_cm: charge_m <=> charge_c
added reaction:  charge_cx: charge_x <=> charge_c
added reaction:  charge_cr: charge_r <=> charge_c
added reaction:  charge_cv: charge_v <=> charge_c
added reaction:  charge_cg: charge_g <=> charge_c
added reaction:  charge_cn: charge_n <=> charge_c

[ ]:

#define the biochemical parameters of the model

tmodel.phi={'ec': Q_(-0.06,'volt'), 'cm': Q_(-0.16,'volt'),'cx':Q_(0,'volt'),'cr':Q_(0,'volt'),'cv':Q_(0,'volt'),'cg':Q_(0,'volt'),'cn':Q_(0,'volt')}
tmodel.pH={'c':Q_(7),'e':Q_(5),'m':Q_(7.4),'r':Q_(7),'v':Q_(5.5),'g':Q_(7),'x':Q_(7),'n':Q_(7)} #On the intracellular pH of baker’s yeast(r), Martinez-Munoz et al. 2008(v), Juan Llopis et al 1998 (g), 7 is default pH

#change max drGcalc_drGt
tmodel._max_drG=Q_(1e6,'kJ/mol') #max is already 1e8 shouldnt need to change

Step 2 - Definition of metabolites

[ ]:

tmodel.get_compounds(search = True)

[██████████████████......................] 568/1236 gsn_m

C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: dolichol_c has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'C': 15, 'H': 28, 'O': 1}, eQuilibrator:{'C': 25, 'O': 1}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: dolmanp_r has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'C': 21, 'H': 38, 'O': 9, 'P': 1}, eQuilibrator:{'C': 31, 'O': 9, 'P': 1}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: glycogen_c has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'C': 6, 'H': 10, 'O': 5}, eQuilibrator:{'O': 21, 'C': 24}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: glycogen_v has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'C': 6, 'H': 10, 'O': 5}, eQuilibrator:{'O': 21, 'C': 24}. Defining as unknown compound
  tmodel.get_compounds(search = True)

[███████████████████████████.............] 843/1236 pmtcoa_c

C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: mhpglu_c has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'C': 25, 'H': 36, 'N': 8, 'O': 12}, eQuilibrator:{'C': 30, 'N': 9, 'O': 12}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: macchitppdol_g has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'C': 37, 'H': 64, 'N': 2, 'O': 22, 'P': 2}, eQuilibrator:{'C': 47, 'O': 22, 'P': 2, 'N': 2}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: hpglu_c has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'C': 24, 'H': 34, 'N': 8, 'O': 12}, eQuilibrator:{'N': 9, 'C': 29, 'O': 12}. Defining as unknown compound
  tmodel.get_compounds(search = True)

[███████████████████████████████████████.] 1221/1236 zym_int1_c

C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdox_c has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdox_m has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdox_n has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdox_x has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdrd_c has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1, 'H': 2}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdrd_m has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1, 'H': 2}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdrd_n has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1, 'H': 2}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)
C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\2562127799.py:1: UserWarning: trdrd_x has inconsistent formula with eQuilibrator identified metabolite. COBRA:{'X': 1, 'H': 2}, eQuilibrator:{'N': 4, 'C': 10, 'S': 2, 'O': 4}. Defining as unknown compound
  tmodel.get_compounds(search = True)

[████████████████████████████████████████] 1236/1236 charge_n

[]

[ ]:

#define biomass parameters
#if the biomass formula is not automatically calcualted it's important to specifiy the hydrogen atoms in biomass to ensure correct proton balancing
tmodel.metabolites.biomass_c.formula = 'H67'
tmodel.metabolites.biomass_e.formula = 'H67'

#manually define charge of cytochrome C so electron movement is correctly calcualted
tmodel.metabolites.focytc_m.charge = 2
tmodel.metabolites.focytc_m.redox= True #define as redox metabolite
tmodel.metabolites.ficytc_m.charge = 3
tmodel.metabolites.ficytc_m.redox = True

Step 3 - metabolites formation energies

When specifying the value of dfGprime for the biomass metabolite, we should plug in the values in units of J/gDW. This differs from all the other reactions in the model, for which values of kJ/mol are used.

The reason for this discrepancy has to do with the units of the fluxes.

For a normal reaction: $g_{diss, rxn} = \Delta_rG \times v = (kJ/mol) \times (mmol/gDW/h) = (10^3 J/mol) \times (10^{-3} mol/gDW/h) = J/gDW/h$

For the biomass reaction: $g_{diss, bio} = \DeltafG{bio} \times `:nbsphinx-math:mu = (J/gDW) :nbsphinx-math:times `(gDW/gDW/h) = J/gDW/h $

Since thermoflux handles all reactions the same when determining the gibbs energy dissipation rate (multiplying $v$ by $\Delta_rG$), we have to pretend that the biomass formation energy has units of kJ/mol. E.g: tmodel.metabolites.biomass_c.dfGprime = Q_(701.767, "kJ/mol") (where 701.767 is the formation energy of E.coli biomass at cytosolic pH, in units of J/gDW)

[8]:

#define biomass properties
tmodel.reactions.BIOMASS_SC5_notrace.id='biomass_c' #update biomass reaction name for convenience

tmodel.metabolites.biomass_c.dfG0 = Q_(-3.04,'kJ/mol')*1000 #From Battley 1991, *1000 for mmol flux units (actual units here are J/gDW)
tmodel.metabolites.biomass_e.dfG0 = Q_(-3.04,'kJ/mol')*1000

Step 4 - Transporters

[9]:

#idenfity transporters where the transported metabolite is empty
#these could either be redox reactions or proton pumps just moving electrons (charge) and protons
#or they could be transporters with chemcial transformation during the transport process
for reaction in tmodel.reactions:
    if len(reaction.compartments) > 1:
        if thermo_flux.tools.drg_tools.calc_transported_mets(reaction) == {}:
            print(reaction.id, reaction, thermo_flux.tools.drg_tools.calc_transported_mets(reaction))

ASPOcm ASPOcm: asp__L_c + fad_m --> fadh2_m + h_c + iasp_c {}
ATPS ATPS: atp_c + h2o_c --> adp_c + h_e + pi_c {}
DHORD4i DHORD4i: dhor__S_c + q6_m --> orot_c + q6h2_m {}
DOLPMTcer DOLPMTcer: dolp_c + gdpmann_c --> dolmanp_r + gdp_c {}
DXHPScm DXHPScm: h2o_c + q6_m + spmd_c --> 13dampp_c + 4abutn_c + q6h2_m {}
D_LACDcm D_LACDcm: 2.0 ficytc_m + lac__D_c --> 2.0 focytc_m + pyr_c {}
FDNG FDNG: for_c + h_c + q6_m --> co2_c + q6h2_m {}
FRDcm FRDcm: fadh2_m + fum_c --> fad_m + succ_c {}
L_LACD2cm L_LACD2cm: 2.0 ficytc_m + lac__L_c --> 2.0 focytc_m + pyr_c {}
NADH2_u6cm NADH2_u6cm: h_c + nadh_c + q6_m --> nad_c + q6h2_m {}

[10]:

#DOLPMTcer involved chemcial transformation of dolp. We assume transport of dolmanp for thermodynamic calcualtion
tmodel.reactions.DOLPMTcer.transported_mets = {tmodel.metabolites.dolmanp_r: -1}

Step 5 - charge and proton balancing

[11]:

for rxn in tmodel.reactions:
    thermo_flux.tools.drg_tools.reaction_balance(rxn, balance_charge=True, balance_mg=False)

C:\Users\tedns\AppData\Local\Temp\ipykernel_1896\1541778503.py:2: UserWarning: ATPS3v is not balanced and could not be automatically balanced, please check reaction stoichiometry or define reaction.transported_h
  thermo_flux.tools.drg_tools.reaction_balance(rxn, balance_charge=True, balance_mg=False)

[12]:

#ATPS3v is vacuolar ATP synthase transporting three protons per ATP
tmodel.reactions.ATPS3v

[12]:

Reaction identifier	ATPS3v
Name	ATP synthase vacuole
Memory address	0x16d19f17f40
Stoichiometry	adp_v + 3.0 h_c + pi_v --> atp_v + h2o_v + 2.0 h_v ADP C10H12N5O10P2 + 3.0 H+ + Phosphate --> ATP C10H12N5O13P3 + H2O H2O + 2.0 H+
GPR	YBR127C and YDL185W and YEL027W and YEL051W and YGR020C and YHR026W and YHR039C_A and YKL080W and...
Lower bound	0.0
Upper bound	500

[13]:

tmodel.reactions.ATPS3v.transported_h = {'v': 3.0, 'c': -3.0}
thermo_flux.tools.drg_tools.reaction_balance(tmodel.reactions.ATPS3v, balance_charge=True, balance_mg=False)
tmodel.reactions.ATPS3v

[13]:

Reaction identifier	ATPS3v
Name	ATP synthase vacuole
Memory address	0x16d19f17f40
Stoichiometry	adp_v + 3.0 charge_c + 3.0 h_c + pi_v --> atp_v + 3.0 charge_v + h2o_v + 3.021815 h_v ADP C10H12N5O10P2 + 3.0 charge + 3.0 H+ + Phosphate --> ATP C10H12N5O13P3 + 3.0 charge + H2O H2O + 3.021815 H+
GPR	YBR127C and YDL185W and YEL027W and YEL051W and YGR020C and YHR026W and YHR039C_A and YKL080W and...
Lower bound	0.0
Upper bound	500

Step 6 - Calculation of Gibbs energy of reactions

Add thermodynamic constraints and parameters

[14]:

#change max flux to 500 for all reactions

tmodel.reactions.get_by_id("biomass_EX").bounds = (0,500)
tmodel.reactions.get_by_id("biomass_c").bounds = (0,500)
tmodel.reactions.get_by_id("biomass_ce").bounds = (0,500)


for rxn in tmodel.reactions:
    rxn.lower_bound = -500 if rxn.lower_bound == -1000 else rxn.lower_bound
    rxn.upper_bound = 500 if rxn.upper_bound == 1000 else rxn.upper_bound
tmodel.objective = tmodel.reactions.biomass_EX

[15]:

tmodel.update_thermo_info(fit_unknown_dfG0=False,search=True)

Identifying compounds...
[████████████████████████████████████████] 1244/1244 Mg_n

Estimating dfG0'...
ficytc_m 0 kilojoule / mole
focytc_m 0 kilojoule / mole
biomass_c 0 kilojoule / mole
biomass_e 0 kilojoule / mole
[████████████████████████████████████████] 1244/1244 Mg_n

Estimating drG0'...
[████████████████████████████████████████] 1587/1587 charge_cn

[16]:

for met in tmodel.metabolites:
    if met.ignore_conc:
        print(met.id,met.ignore_conc)

h2o_c True
h2o_e True
h2o_g True
h2o_m True
h2o_n True
h2o_r True
h2o_v True
h2o_x True
h_c True
h_e True
h_g True
h_m True
h_n True
h_r True
h_v True
h_x True
oh1_c True
oh1_m True
biomass_c True
biomass_e True
charge_c True
charge_e True
charge_m True
charge_x True
charge_r True
charge_v True
charge_g True
charge_n True

[17]:

for rxn in tmodel.reactions:
    if rxn.ignore_snd:
        print(rxn.id,rxn.ignore_snd)

EX_epistest_SC_e True
EX_epist_e True
EX_ergst_e True
EX_ergstest_SC_e True
EX_13BDglcn_e True
EX_etha_e True
EX_2hb_e True
etoh_EX True
EX_fe2_e True
EX_fecost_e True
EX_2mbac_e True
EX_2mbald_e True
EX_2mbtoh_e True
EX_2mppal_e True
EX_2phetoh_e True
EX_fecostest_SC_e True
EX_fmn_e True
EX_3c3hmp_e True
EX_3mbald_e True
EX_for_e True
EX_fru_e True
EX_3mop_e True
EX_4abut_e True
EX_fum_e True
EX_g3pc_e True
EX_4abz_e True
EX_5aop_e True
EX_g3pi_e True
EX_gal_e True
EX_8aonn_e True
EX_Nbfortyr_e True
EX_abt_e True
EX_ac_e True
EX_acald_e True
EX_aces_e True
EX_galur_e True
EX_gam6p_e True
EX_gcald_e True
EX_glc__D_e True
EX_gln__L_e True
EX_glu__L_e True
EX_glx_e True
EX_ade_e True
EX_adn_e True
EX_akg_e True
EX_ala__L_e True
EX_gly_e True
EX_glyc_e True
EX_gsn_e True
EX_gthox_e True
EX_gthrd_e True
EX_alltn_e True
EX_alltt_e True
EX_amet_e True
EX_arab__D_e True
EX_arab__L_e True
EX_gua_e True
EX_arg__L_e True
EX_h2o_e True
EX_h_e True
EX_hdca_e True
EX_hdcea_e True
EX_asn__L_e True
EX_asp__L_e True
EX_btd_RR_e True
EX_btn_e True
EX_hexc_e True
EX_his__L_e True
EX_hxan_e True
EX_iamac_e True
EX_iamoh_e True
EX_ibutac_e True
EX_chol_e True
EX_cit_e True
co2_EX True
EX_crn_e True
EX_csn_e True
EX_cys__L_e True
EX_cytd_e True
EX_dad_2_e True
EX_dann_e True
EX_ibutoh_e True
EX_id3acald_e True
EX_ile__L_e True
EX_ind3eth_e True
EX_inost_e True
EX_dca_e True
EX_dcyt_e True
EX_ddca_e True
EX_ins_e True
EX_k_e True
EX_dgsn_e True
EX_din_e True
EX_lac__D_e True
EX_lac__L_e True
EX_lanost_e True
EX_lanostest_SC_e True
EX_dttp_e True
EX_duri_e True
EX_leu__L_e True
EX_lys__L_e True
EX_mal__L_e True
EX_malt_e True
EX_ribflv_e True
EX_sbt__D_e True
EX_man_e True
EX_melib_e True
EX_met__L_e True
EX_mmet_e True
EX_na1_e True
EX_nac_e True
EX_sbt__L_e True
EX_ser__L_e True
EX_so3_e True
EX_so4_e True
EX_spmd_e True
EX_sprm_e True
EX_srb__L_e True
EX_succ_e True
EX_nadp_e True
EX_nh4_e True
EX_nmn_e True
EX_sucr_e True
o2_EX True
EX_taur_e True
EX_thm_e True
EX_oaa_e True
EX_ocdca_e True
EX_ocdcea_e True
EX_ocdcya_e True
EX_orn_e True
EX_thmmp_e True
EX_thmpp_e True
EX_thr__L_e True
EX_thym_e True
EX_pacald_e True
EX_pap_e True
EX_thymd_e True
EX_tre_e True
EX_trp__L_e True
EX_pc_SC_e True
EX_pectin_e True
EX_pepd_e True
EX_phe__L_e True
EX_pheac_e True
EX_pi_e True
EX_pnto__R_e True
EX_ttdca_e True
EX_tyr__L_e True
EX_ura_e True
EX_urea_e True
EX_uri_e True
EX_val__L_e True
EX_xan_e True
EX_xtsn_e True
EX_pro__L_e True
EX_ptd1ino_SC_e True
EX_ptrc_e True
EX_pyr_e True
EX_xyl__D_e True
EX_xylt_e True
EX_rib__D_e True
EX_zymst_e True
EX_zymstest_SC_e True
H2Ot True
H2Oter True
H2Otm True
H2Otn True
H2Otp True
H2Otv True
biomass_c True
biomass_ce True
biomass_EX True
EX_charge True

[18]:

#ignore snd for charge ?
for rxn in tmodel.reactions:
    if 'charge' in rxn.id and 'EX' not in rxn.id:#charge reactions are constrained by 2nd law
        print(rxn.id,rxn.ignore_snd)
        rxn.ignore_snd=True
        print(rxn.id,rxn.ignore_snd)


#tmodel.reactions.EX_charge.ignore_snd=True

charge_ce False
charge_ce True
charge_cm False
charge_cm True
charge_cx False
charge_cx True
charge_cr False
charge_cr True
charge_cv False
charge_cv True
charge_cg False
charge_cg True
charge_cn False
charge_cn True

Step 7 - Establishment of the thermodynamic-stoichiometric solution space

Note : In our curation of the iMM904 model, due to the limited information available for the membrane potential of organelles, we chose to allow the free leakage of ions (charges) across the membrane. Allowing free movement of charges between compartments prevents over-constraining the model to potentially incorrect information. For more detailed analysis of charge transport or models where more information is available on membrane physiology then a constraint on the second law could be added for these charge transport reactions.

[25]:

# add the TFBA variables to the model
tmodel.m = None #reset the gurobi model object in case you're re-running this cell
tmodel.add_TFBA_variables(gdiss_constraint = False, sigmac_limit = (5000/tmodel.T.m), error_type = 'covariance',qnorm=1,alpha=0.95)
tmodel.m.update()

Set parameter NonConvex to value 2
Set parameter TimeLimit to value 10

[26]:

thermo_flux.solver.gurobi.model_start(tmodel,'fixed_dir_firstTFBA_ENS_8.sol', ignore_vars=['all'],fix_vars=['qm','fluxes'],fix='start')

[27]:

tmodel.m.params.NonConvex=2
tmodel.m.params.IntFeasTol=1e-6
tmodel.m.params.OptimalityTol=1e-6
tmodel.m.params.TimeLimit=60*1
tmodel.m.optimize()

Set parameter IntFeasTol to value 1e-06
Set parameter TimeLimit to value 60
Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (win64 - Windows 11+.0 (22631.2))

CPU model: Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 15685 rows, 12850 columns and 350518 nonzeros
Model fingerprint: 0x1b1ad44f
Model has 1 general constraint
Variable types: 11263 continuous, 1587 integer (1587 binary)
Coefficient statistics:
  Matrix range     [6e-06, 1e+06]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+06]
  RHS range        [1e-05, 1e+06]

Warning: Completing partial solution with 1587 unfixed non-continuous variables out of 1587
User MIP start produced solution with objective 2.14592 (0.30s)
Loaded user MIP start with objective 2.14592

Presolve removed 8574 rows and 6443 columns
Presolve time: 1.61s
Presolved: 7111 rows, 6407 columns, 277793 nonzeros
Variable types: 5169 continuous, 1238 integer (1238 binary)

Root simplex log...

Iteration    Objective       Primal Inf.    Dual Inf.      Time
    6208    1.2993797e+00   3.446101e+05   0.000000e+00      5s
   12785    2.1805119e+00   0.000000e+00   0.000000e+00      8s

Root relaxation: objective 2.180512e+00, 12785 iterations, 5.59 seconds (7.85 work units)
Total elapsed time = 31.75s (DegenMoves)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    2.18051    0   43    2.14592    2.18051  1.61%     -   32s
     0     0    2.18051    0   47    2.14592    2.18051  1.61%     -   32s
     0     0    2.18051    0   46    2.14592    2.18051  1.61%     -   33s
     0     0    2.18051    0   55    2.14592    2.18051  1.61%     -   36s
     0     0    2.18051    0   49    2.14592    2.18051  1.61%     -   39s
     0     0    2.18051    0   19    2.14592    2.18051  1.61%     -   44s
     0     0    2.18051    0   18    2.14592    2.18051  1.61%     -   45s
     0     0    2.18051    0   14    2.14592    2.18051  1.61%     -   48s
     0     0    2.18051    0   14    2.14592    2.18051  1.61%     -   49s
     0     0    2.18051    0   19    2.14592    2.18051  1.61%     -   52s
     0     0    2.18043    0   20    2.14592    2.18043  1.61%     -   53s
     0     0    2.18043    0   19    2.14592    2.18043  1.61%     -   53s
     0     0    2.18043    0   14    2.14592    2.18043  1.61%     -   55s
     0     0    2.18042    0   10    2.14592    2.18042  1.61%     -   56s
     0     0    2.18039    0   11    2.14592    2.18039  1.61%     -   57s
     0     0    2.18039    0    9    2.14592    2.18039  1.61%     -   58s
     0     0    2.18039    0    9    2.14592    2.18039  1.61%     -   58s
     0     0    2.18038    0    9    2.14592    2.18038  1.61%     -   58s

Cutting planes:
  Learned: 60
  Cover: 3
  Implied bound: 10
  MIR: 13
  Flow cover: 33
  Relax-and-lift: 4

Explored 1 nodes (49500 simplex iterations) in 60.25 seconds (53.86 work units)
Thread count was 8 (of 8 available processors)

Solution count 1: 2.14592

Time limit reached
Best objective 2.145918254336e+00, best bound 2.180382715743e+00, gap 1.6060%

[28]:

tmodel.m.write('fixed_dir_firstTFBA_ENS_9.mps')
tmodel.m.write('fixed_dir_firstTFBA_ENS_9.sol')

[45]:

thermo_flux.solver.gurobi.gdiss_model(tmodel)

drG error term
Ex: 143.54648665524567
Int: -143.54620102128138
Diff: 0.00028563396429603927


RTlnC concentration term
Ex: -257.1979974465446
Int: 257.19799744655256
Diff: 7.958078640513122e-12


drG0' term
Ex: 28160.88767254069
Int: -28160.88767254047
Diff: 2.219167072325945e-10


drG0 term non-transformed
Ex: 27952.820538559
Int: -27952.820538558746
Diff: 2.546585164964199e-10


drG0' transform term
Ex: 0.0
Int: 0.0
Diff: 0.0


drG0' total transport term
Ex: 0.0
Int: -0.0007508681337640155
Diff: -0.0007508681337640155


drG0' charge transport term
Ex: 0.0
Int: 7.466136594302952e-05
Diff: 7.466136594302952e-05


drG0'proton transport term
Ex: 0.0
Int: -0.0008255294998775753
Diff: -0.0008255294998775753

[29]:

# add the TFBA variables to the model
tmodel.m = None #reset the gurobi model object in case you're re-running this cell
tmodel.add_TFBA_variables(gdiss_constraint = True, sigmac_limit = (3700/tmodel.T.m), error_type = 'covariance',qnorm=1,alpha=0.95)
tmodel.m.update()
tmodel.m.params.timelimit=60*5

GUR_range = np.linspace(-0, -23, 50)

thermo_flux.solver.gurobi.variable_scan(tmodel, GUR_range, var = tmodel.mvars['v'][0][tmodel.reactions.index(tmodel.reactions.get_by_id('EX_glc__D_e'))]) # here we get the index for the glucose uptake reaction and use this to index the v variable in the gurobi model
thermo_flux.solver.gurobi.model_start(tmodel,'fixed_dir_firstTFBA_ENS_9.sol',ignore_vars=['all'],fix_vars=['qm'],fix='bound') #fix the drG error term to that from the initial TFBA to speed up optimisations

tmodel.m.optimize()

fluxes_gdiss = thermo_flux.solver.gurobi.multi_scenario_sol(tmodel,var = 'v')
fluxes_gdiss= pd.DataFrame(fluxes_gdiss[0].T, index = [rxn.id for rxn in tmodel.reactions])

Set parameter NonConvex to value 2
Set parameter TimeLimit to value 10
Set parameter TimeLimit to value 300
Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (win64 - Windows 11+.0 (22631.2))

CPU model: Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 15686 rows, 13017 columns and 350685 nonzeros
Model fingerprint: 0x659902df
Model has 166 quadratic constraints
Model has 1 general constraint
Variable types: 11430 continuous, 1587 integer (1587 binary)
Coefficient statistics:
  Matrix range     [6e-06, 1e+06]
  QMatrix range    [1e+00, 1e+00]
  QLMatrix range   [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [5e-05, 1e+06]
  RHS range        [1e-05, 1e+06]

Solving a multi-scenario model with 50 scenarios...

Presolve removed 11125 rows and 8819 columns
Presolve time: 0.28s
Presolved: 4888 rows, 4248 columns, 15024 nonzeros
Presolved model has 81 bilinear constraint(s)
Presolved model has 50 scenario(s)

Solving non-convex MIQCP

Variable types: 3286 continuous, 962 integer (962 binary)

Root relaxation: objective 1.202063e+00, 2935 iterations, 0.20 seconds (0.10 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    1.20206    0   50          -    1.20206      -     -    0s
     0     0    1.20206    0   54          -    1.20206      -     -    1s
     0     0    1.20206    0   54          -    1.20206      -     -    1s
     0     0    1.19820    0   30          -    1.19820      -     -    1s
     0     0    1.19820    0   30          -    1.19820      -     -    1s
     0     0    1.19820    0   31          -    1.19820      -     -    1s
     0     0    1.19820    0   32          -    1.19820      -     -    1s
     0     0    1.19820    0   30          -    1.19820      -     -    1s
     0     0    1.19820    0   30          -    1.19820      -     -    1s
     0     0    1.19820    0   27          -    1.19820      -     -    2s
     0     0    1.19820    0   28          -    1.19820      -     -    2s
     0     0    1.19820    0   27          -    1.19820      -     -    2s
     0     0    1.02710    0   27          -    1.02710      -     -    3s
H    0     0                      -0.0000000    1.02710      -     -    3s
     0     2    1.02710    0   27   -0.00000    1.02710      -     -    4s
    38    45   -0.00000   11    8   -0.00000    1.02700      -  23.6    5s
H  839   783                       0.0000000    1.02570      -  12.8   14s
H  840   783                       0.0000000    1.02570      -  12.8   14s
H  849   811                       0.1766679    1.02570   481%  12.9   14s
H  867   811                       0.3850917    1.02570   166%  12.9   14s
   871   853    0.51994   72   11    0.38509    1.02570   166%  13.0   16s
H 1131  1054                       0.3914907    1.02570   162%  13.1   18s
  1135   421    0.39149   31   35    0.39149    1.02570   162%  13.1   20s
  1145   172    0.40662  149   42    0.39149    0.40886  4.44%  17.3   26s
  1152    72    0.39149   66   40    0.39149    0.40886  4.44%  22.1   30s
Scenario 3 has been solved (gap 0.0020%). 48/50 scenarios left.
  1154    49   -0.00000   38   25    0.39149    0.40886  4.44%  27.2   38s
  1322   169   -0.00000   62   19    0.39149    0.40885  4.43%  25.9   40s
  2838  1705    0.25606  150   32    0.39149    0.40690  3.94%  19.7   45s
H 3486  2294                       0.4003544    0.40690  1.64%  18.6   47s
H 3731  2466                       0.4018819    0.40689  1.25%  18.6   48s
H 3835  2560                       0.4018852    0.40687  1.24%  18.6   48s
H 3954  2693                       0.4018853    0.40687  1.24%  18.4   49s
  4185  3072    0.39547  100   40    0.40189    0.40687  1.24%  17.9   50s
H 5709  4341                       0.4019121    0.40682  1.22%  16.5   54s
  5795  4614    0.38552  108   54    0.40191    0.40682  1.22%  16.3   55s
  7998  6722    0.03982  195   20    0.40191    0.40678  1.21%  15.4   60s
  9842  5587    0.40261  122   12    0.40191    0.40675  1.20%  16.7   65s
*10580  5888             290       0.4019171    0.40675  1.20%  16.4   66s
 12029  7125    0.27138  116   38    0.40192    0.40666  1.18%  16.3   70s
H12935  7608                       0.4022758    0.40657  1.07%  16.8   72s
 14248  8911 infeasible  200         0.40228    0.40652  1.06%  16.8   76s
H15435  9540                       0.4022806    0.40617  0.97%  17.1   80s
H15452  9540                       0.4024703    0.40617  0.92%  17.1   80s
 15506  9577    0.39308  118   29    0.40247    0.40617  0.92%  17.1   85s
Scenario 1 has been solved. 47/50 scenarios left.
Scenario 2 has been solved. 46/50 scenarios left.
Scenario 5 has been solved. 45/50 scenarios left.
Scenario 6 has been solved. 44/50 scenarios left.
 15548  9579    0.39308  119   29    0.40247    0.40617  0.92%  17.1  120s
Scenario 4 has been solved. 43/50 scenarios left.
Scenario 7 has been solved (gap 0.0018%). 42/50 scenarios left.
Scenario 8 has been solved (gap 0.0069%). 41/50 scenarios left.
Scenario 9 has been solved. 40/50 scenarios left.
 16277 10992    0.38948   98   35    0.40247    0.40607  0.89%  17.0  127s
H17122 11379                       0.4024815    0.40567  0.79%  16.6  129s
 17580 11996    0.37760  141   24    0.40248    0.40557  0.77%  16.7  132s
 18324 12778    0.38603  109   43    0.40248    0.40442  0.48%  16.7  135s
 20342 14646    0.35600  143   41    0.40248    0.40355  0.27%  16.3  141s
H21366 14798                       0.4026757    0.40350  0.21%  16.3  143s
H21371 14798                       0.4026888    0.40350  0.20%  16.4  143s
 21587 15652    0.40350   80   36    0.40269    0.40350  0.20%  16.4  146s
 23746 17593    0.35200  129   34    0.40269    0.40350  0.20%  16.3  152s
*24067 17593             171       0.4027681    0.40350  0.18%  16.2  152s
*24759 17593             165       0.4028130    0.40350  0.17%  16.1  152s
 24904 18017    0.36228  164   43    0.40281    0.40350  0.17%  16.1  155s
 26619 19080    0.39076  120   49    0.40281    0.40350  0.17%  16.2  160s
H26693 19080                       0.4029112    0.40350  0.14%  16.2  160s
 27890 20560     cutoff  105         0.40291    0.40345  0.13%  16.3  166s
 28939 21573    0.39425  150   30    0.40291    0.40345  0.13%  16.0  170s
 31013 23127    0.39472  195    4    0.40291    0.40345  0.13%  15.7  175s
 32841 24688    0.40336   86   35    0.40291    0.40341  0.12%  15.6  182s
 33978 25996    0.39189  185   23    0.40291    0.40341  0.12%  15.6  186s
 35466 26987 infeasible  105         0.40291    0.40340  0.12%  15.3  190s
 37764 28344    0.38396  168   25    0.40291    0.40340  0.12%  15.3  195s
 39656 29782 infeasible   95         0.40291    0.40338  0.12%  15.5  202s
 40905 30530    0.40155   87   52    0.40291    0.40338  0.12%  15.5  207s
 41838 30924    0.38574  151   28    0.40291    0.40338  0.12%  15.4  210s
 43650 32430    0.39840  118   40    0.40291    0.40337  0.11%  15.7  216s
 44490 32417    0.40259  133   34    0.40291    0.40337  0.11%  15.7  224s
Scenario 12 has been solved (gap 0.0079%). 39/50 scenarios left.
 44498 33271    0.40259  134   33    0.40291    0.40336  0.11%  15.7  227s
 45548 33460 infeasible  160         0.40291    0.40335  0.11%  15.7  230s
 47272 33928     cutoff  113         0.40291    0.40333  0.10%  16.2  236s
 48281 33922    0.40323  146   39    0.40291    0.40333  0.10%  16.5  241s
 49414 34994    0.40329  139   32    0.40291    0.40332  0.10%  16.6  247s
H49845 34994                       0.4029574    0.40332  0.09%  16.6  247s
 50165 35613    0.39175  104   40    0.40296    0.40332  0.09%  16.6  250s
 52053 36182    0.40326  105   36    0.40296    0.40328  0.08%  16.8  268s
Scenario 11 has been solved (gap 0.0056%). 38/50 scenarios left.
 52092 36781    0.40326  106   37    0.40296    0.40328  0.08%  16.8  272s
 52989 37366     cutoff  113         0.40296    0.40327  0.08%  16.8  276s
 53864 38129    0.39532  203   22    0.40296    0.40324  0.07%  16.8  284s
 55063 38537    0.40320  150   41    0.40296    0.40323  0.07%  16.8  291s
 56072 39425 infeasible  110         0.40296    0.40321  0.06%  16.9  298s

Cutting planes:
  Learned: 18
  Flow cover: 1
  RLT: 1

Explored 57432 nodes (975891 simplex iterations) in 300.49 seconds (56.02 work units)
Thread count was 8 (of 8 available processors)

Solution count 10: 0.402957 0.402957 0.402957 ... 0.402911

Time limit reached
Best objective 4.029573999817e-01, best bound 4.032086474970e-01, gap 0.0624%

[25]:

#plot biomass yield for stoichiometric model with balanced reactions and pH specified
df_data2 = pd.read_csv('plotting_Yeast_ALL_experimental_data.csv',index_col=0)
df_data2['regression'].fillna('Not fitted',inplace=True)
df_data2['regression'].replace('yes', 'Fitted', inplace=True)
df_data2.rename(columns={'regression':'Exp. data'}, inplace=True)
#rename index in Condition
df_data2.index.rename('Condition', inplace=True)

GUR_plot_range=[-x for x in GUR_range]
#plot biomass vs GUR in a scatter plot

sns.scatterplot(x=GUR_plot_range,y=fluxes_gdiss.loc['biomass_EX'],color='b',label='TFBA')
# sns.scatterplot(x=GUR_plot_range,y=tbiomass_flux,color='b',label='FBA and balanced')
# sns.scatterplot(x=GUR_plot_range,y=biomass_flux,color='g',label='FBA')
sns.scatterplot(x=df_data2['glc-D_EX'], y=df_data2['biomass_EX'], color=['red'],label='Experimental data')
plt.xlabel('Glucose uptake rate (mmol/(gDW.h))',fontsize=14)
plt.ylabel('Biomass production rate (mmol/(gDW.h))',fontsize=14)

plt.savefig('Figures/growth_all_TFBA_gdissconstraint.png',dpi=300)