4. Debugging

It is often observed in genome scale models that adding thermodynamic constraints make the model infeasible. This might be due to variety of reasons, like ill-formatted covariance between groups/components, lumping of reactions etc. It is always not easy to find the exact cause of the problem. But, if you are using Gurobi/Cplex, you can find the constraints that render the model infeasible.

Please note that these methods give you random set of constraints that represent infeasibility. Also, a model can contain multiple infeasible sets. While modifying/removing these constraints, it is advisable to look for biological meaning.

Let us demonstrate the usage by breaking the model.

import sys
sys.path.insert(0, "../../examples")
from paper_data_example import build_core_model

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

tfa_model = build_core_model()

Using license file /home/vishnu/gurobi.lic
Academic license - for non-commercial use only
Read LP format model from file /tmp/tmpdzvydb7y.lp
Reading time = 0.00 seconds
: 72 rows, 190 columns, 720 nonzeros

Now lets change the concentration atp and find the objective value. Then, we set the lower bound of objective variable higher than the solution to break the model.

atp = tfa_model.metabolites.get_by_id('atp_c')
atp.concentration_min = 2e-3
atp.concentration_max = 5e-2
solution = tfa_model.optimize()
print(solution)

<Solution 0.374 at 0x7f3c6ad4dda0>

Now lets adjust the biomass lower bound to 0.5 and find out the IIS set of reactions using Gurobi

biomass = tfa_model.reactions.get_by_id('BIOMASS_Ecoli_core_w_GAM')
biomass.lower_bound = 0.5
solution = tfa_model.optimize()
tfa_model.solver.problem.computeIIS() # model.solver.problem represents the solver interface model
for cons in tfa_model.solver.problem.getConstrs():
    if cons.IISConstr:
        print(cons)

<gurobi.Constr glc__D_e>
<gurobi.Constr gln__L_c>
<gurobi.Constr gln__L_e>
<gurobi.Constr glu__L_c>
<gurobi.Constr glu__L_e>
<gurobi.Constr icit_c>
<gurobi.Constr lac__D_c>
<gurobi.Constr lac__D_e>
<gurobi.Constr mal__L_c>
<gurobi.Constr mal__L_e>
<gurobi.Constr 13dpg_c>
<gurobi.Constr 2pg_c>
<gurobi.Constr 3pg_c>
<gurobi.Constr oaa_c>
<gurobi.Constr pep_c>
<gurobi.Constr 6pgc_c>
<gurobi.Constr 6pgl_c>
<gurobi.Constr pyr_c>
<gurobi.Constr pyr_e>
<gurobi.Constr r5p_c>
<gurobi.Constr ru5p__D_c>
<gurobi.Constr acald_c>
<gurobi.Constr s7p_c>
<gurobi.Constr acald_e>
<gurobi.Constr accoa_c>
<gurobi.Constr succ_c>
<gurobi.Constr succ_e>
<gurobi.Constr succoa_c>
<gurobi.Constr acon_C_c>
<gurobi.Constr xu5p__D_c>
<gurobi.Constr actp_c>
<gurobi.Constr adp_c>
<gurobi.Constr akg_c>
<gurobi.Constr akg_e>
<gurobi.Constr amp_c>
<gurobi.Constr cit_c>
<gurobi.Constr dhap_c>
<gurobi.Constr e4p_c>
<gurobi.Constr etoh_c>
<gurobi.Constr etoh_e>
<gurobi.Constr f6p_c>
<gurobi.Constr fdp_c>
<gurobi.Constr fru_e>
<gurobi.Constr fum_c>
<gurobi.Constr fum_e>
<gurobi.Constr g3p_c>
<gurobi.Constr g6p_c>
<gurobi.Constr directionality_ACONTb>
<gurobi.Constr ind_ACONTb>
<gurobi.Constr delG_ACONTb>
<gurobi.Constr directionality_ATPS4r>
<gurobi.Constr ind_ATPS4r>
<gurobi.Constr delG_ATPS4r>

This returned us a list of constraints that can be modified to make the model status optimal. One has to choose and modify the list of constraints to manipulate based on biological meaning.