MFA data as constraints for a COBRA model¶
Here we jump into the (for now) final step in our metabolic flux analysis (MFA) pipeline. After constructing an INCA script, running it in MATLAB and reimporting the data we’re now here. This example notebook will guide you through different ways to integrate your MFA results into COBRA models and how to make them more reliable.
[1]:
import pickle
import pandas as pd
import escher
import cobra
from cobra import Reaction
from tabulate import tabulate
from gurobipy import Model as GRBModel
from BFAIR.mfa.INCA import INCA_reimport
from BFAIR.mfa.sampling import (
add_constraints,
add_feasible_constraints,
find_biomass_reaction,
get_min_solution_val,
replace_biomass_rxn_name,
bound_relaxation,
)
from BFAIR.mfa.visualization import (
reshape_fluxes_escher,
)
Determination of memory status is not supported on this
platform, measuring for memoryleaks will never fail
Preparation¶
This step is detailed in the MFA_compatibility notebook. Please consult it for more detail, we will go through the steps without any explanation here
[2]:
fittedFluxes = pd.read_pickle("data/MFA_sampling/preprocessed_fittedFluxes.obj")
[3]:
model = cobra.io.load_json_model("data/MFA_sampling/preprocessed_model.json")
Academic license - for non-commercial use only - expires 2021-07-30
Using license file /Users/matmat/gurobi.lic
Let’s copy the model so we won’t have to go through the pre-processing steps again
[4]:
model_original = model.copy()
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpj84wal8o.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpz4icoou0.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Re-integration¶
Now let’s have a look at our model
[5]:
original_solution = model.optimize()
original_solution
[5]:
| fluxes | reduced_costs | |
|---|---|---|
| EX_cm_e | 0.000000 | 0.000000e+00 |
| EX_cmp_e | 0.000000 | -2.965572e-01 |
| EX_co2_e | 19.675223 | 0.000000e+00 |
| EX_cobalt2_e | -0.000025 | -0.000000e+00 |
| DM_4crsol_c | 0.000219 | 0.000000e+00 |
| ... | ... | ... |
| PTAr_reverse | 0.000000 | 0.000000e+00 |
| ACONTb_reverse | 0.000000 | 0.000000e+00 |
| PGK_reverse | 0.000000 | -2.775558e-17 |
| ACKr_reverse | 0.000000 | 6.938894e-17 |
| ACS_reverse | 0.000000 | 5.551115e-17 |
2593 rows × 2 columns
Let’s see what happens now when we add our new bound constraints
[6]:
model = add_constraints(model, fittedFluxes)
model.optimize()
--- start ---
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpg8vvow6j.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpu2uzuwjf.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Did not work for EX_co2_e_unlabeled
Did not work for EX_glc_e
Did not work for FADR
Did not work for GluSYN
Did not work for HisSYN
Did not work for IleSYN
Did not work for LeuSYN
Did not work for MetSYN
Did not work for MlthfSYN
Did not work for MlthfSYN_reverse
Did not work for NADH
Did not work for PheSYN
Did not work for ProSYN
Did not work for SerSYN
Did not work for SUCCOAS
Did not work for ThrSYN
Did not work for TKT1a
Did not work for TKT1b
Did not work for TKT2a
Did not work for TKT2b
Did not work for TrpSYN
Did not work for TyrSYN
Did not work for ValSYN
Did not work for CYTBD
Did not work for HYD
Did not work for ATPS4r
add_constraints takes 0h: 0min: 10sec to run
--- end ---
/Users/matmat/opt/anaconda3/envs/bfair/lib/python3.8/site-packages/cobra/util/solver.py:430: UserWarning: solver status is 'infeasible'
warn("solver status is '{}'".format(status), UserWarning)
[6]:
Oh… That’s not good… Well that sucks, seems like we have to deal with an infeasible solution. There are two straight forward ways: exclusion and relaxation
Dealing with infeasible solutions - exclusion¶
The easier way to deal with this issue is to simply exclude the constraints that render a model infeasible. We can do that by adding the calculated bounds one by one. If we come across a reaction whose bounds cause trouble, we restart the process and skip this one. This might have to be done a few times to exclude all troublemakers. the add_feasible_constraints() functions takes care of that for us. Let’s reset the model first.
[7]:
model = model_original.copy()
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp4a1vc36x.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp21n47dhv.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
[8]:
min_val = get_min_solution_val(fittedFluxes, biomass_string='BIOMASS')
[9]:
model, problems = add_feasible_constraints(model, fittedFluxes, min_val=min_val)
--- start ---
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp9p9o82ui.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpw2qb3oh4.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Solution infeasible if adding ASPTA
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp99_809un.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpaxl5yla7.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Solution infeasible if adding DAPDC
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp0p1ss7qp.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpemru1hi4.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Solution infeasible if adding BIOMASS_Ec_iJO1366_core_53p95M
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpq6c6eaas.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpxo7ffcyi.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Did not work for EX_co2_e_unlabeled
Did not work for EX_glc_e
Solution infeasible if adding EX_nh4_e
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp48jx61rv.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpg8545zd2.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Did not work for EX_co2_e_unlabeled
Did not work for EX_glc_e
Solution infeasible if adding EX_o2_e
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmppsuzp_y1.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp9au2uh20.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Did not work for EX_co2_e_unlabeled
Did not work for EX_glc_e
Solution infeasible if adding EX_so4_e
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp7m3xb_ha.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpq_l26let.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Did not work for EX_co2_e_unlabeled
Did not work for EX_glc_e
Did not work for FADR
Solution infeasible if adding GLNS
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpx8osj2na.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmp1tz0jylq.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Did not work for EX_co2_e_unlabeled
Did not work for EX_glc_e
Did not work for FADR
Did not work for GluSYN
Did not work for HisSYN
Did not work for IleSYN
Did not work for LeuSYN
Did not work for MetSYN
Did not work for MlthfSYN
Did not work for MlthfSYN_reverse
Did not work for NADH
Did not work for PheSYN
Did not work for ProSYN
Did not work for SerSYN
Did not work for SUCCOAS
Did not work for ThrSYN
Did not work for TKT1a
Did not work for TKT1b
Did not work for TKT2a
Did not work for TKT2b
Did not work for TrpSYN
Did not work for TyrSYN
Did not work for ValSYN
Did not work for CYTBD
Did not work for HYD
Did not work for ATPS4r
---------------------------
Total number of restarts: 7
add_feasible_constraints takes 0h: 1min: 39sec to run
--- end ---
The model is our newly constrained model and the problematic reactions can be listed in problems.
[10]:
problems
[10]:
['ASPTA',
'DAPDC',
'BIOMASS_Ec_iJO1366_core_53p95M',
'EX_nh4_e',
'EX_o2_e',
'EX_so4_e',
'GLNS']
Now let’s see what an effect these new bounds had on the predicted growth rate (the objective value) of our model
[11]:
new_bounds_solution = model.optimize()
new_bounds_solution
[11]:
| fluxes | reduced_costs | |
|---|---|---|
| EX_cm_e | 0.000000 | 0.0 |
| EX_cmp_e | 0.000000 | -0.0 |
| EX_co2_e | 15.000000 | -0.0 |
| EX_cobalt2_e | -0.000020 | -0.0 |
| DM_4crsol_c | 0.000180 | 0.0 |
| ... | ... | ... |
| PTAr_reverse | 0.000000 | 0.0 |
| ACONTb_reverse | 0.000000 | 0.0 |
| PGK_reverse | 23.706810 | 0.0 |
| ACKr_reverse | 4.783906 | 0.0 |
| ACS_reverse | 0.000000 | 0.0 |
2593 rows × 2 columns
And here’s the star of the show, our sampling method. We trust our models because… we have to! And because smart people that knew what they were doing set them up. So in order to gain more confidence in our MFA data, we sample the model after adding the calculated bound for some of the reactions and re-calculate the fluxes a number of time. Then, we take the mean and take that as the most trustworthy calculated flux. These fluxes can be visualized, for example in tools like Escher
[12]:
sampled_fluxes = cobra.sampling.sample(model, n=100, processes=2)
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpxzgt2q8z.lp
Reading time = 0.03 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpt_ozl2vr.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
[13]:
sampled_fluxes
[13]:
| EX_cm_e | EX_cmp_e | EX_co2_e | EX_cobalt2_e | DM_4crsol_c | DM_5drib_c | DM_aacald_c | DM_amob_c | DM_mththf_c | EX_colipa_e | ... | ACONTa_reverse | FUM_reverse | GAPD_reverse | ICDHyr_reverse | PGM_reverse | PTAr_reverse | ACONTb_reverse | PGK_reverse | ACKr_reverse | ACS_reverse | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.0 | 0.0 | 14.651387 | -2.169991e-08 | 1.943322e-07 | 0.011272 | 0.0 | 1.742395e-09 | 0.000171 | 0.000890 | ... | 2.670170 | 1.050728 | 0.016231 | 1.627122 | 16.664920 | 26.785111 | 2.670025 | 19.447786 | 1.010911 | 27.001166 |
| 1 | 0.0 | 0.0 | 14.464957 | -1.890907e-08 | 1.687986e-07 | 0.010888 | 0.0 | 1.513428e-09 | 0.001666 | 0.000545 | ... | 2.713958 | 1.062594 | 0.016960 | 1.520385 | 16.611964 | 26.683295 | 2.713823 | 19.619448 | 1.110196 | 27.091069 |
| 2 | 0.0 | 0.0 | 14.478298 | -1.882824e-08 | 1.680482e-07 | 0.010887 | 0.0 | 1.506698e-09 | 0.001312 | 0.000545 | ... | 2.668646 | 1.003182 | 0.017329 | 1.453261 | 16.614507 | 26.648677 | 2.668510 | 19.624369 | 1.136609 | 27.066028 |
| 3 | 0.0 | 0.0 | 14.458260 | -1.812741e-08 | 1.621550e-07 | 0.010887 | 0.0 | 1.453843e-09 | 0.011918 | 0.000534 | ... | 2.511191 | 0.799086 | 0.039586 | 1.260704 | 16.625832 | 26.738789 | 2.511056 | 19.518556 | 1.222329 | 27.185154 |
| 4 | 0.0 | 0.0 | 14.460898 | -1.867411e-08 | 1.672003e-07 | 0.010934 | 0.0 | 1.491510e-09 | 0.012093 | 0.000532 | ... | 2.523171 | 0.807694 | 0.046294 | 1.269027 | 16.625135 | 26.776248 | 2.521616 | 19.518551 | 1.198614 | 27.181841 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 95 | 0.0 | 0.0 | 14.445930 | -1.365700e-06 | 1.253752e-05 | 0.004238 | 0.0 | 1.124099e-07 | 0.003438 | 0.000481 | ... | 14.699352 | 12.161994 | 1.340203 | 8.881692 | 19.998642 | 6.149870 | 14.686487 | 21.464382 | 9.080217 | 8.795291 |
| 96 | 0.0 | 0.0 | 14.518657 | -1.339059e-06 | 1.223887e-05 | 0.004255 | 0.0 | 1.097314e-07 | 0.003518 | 0.000247 | ... | 14.879652 | 12.508729 | 1.532087 | 9.058146 | 20.068367 | 5.958601 | 14.866551 | 21.493044 | 8.840546 | 8.423002 |
| 97 | 0.0 | 0.0 | 14.503939 | -1.333918e-06 | 1.219152e-05 | 0.004234 | 0.0 | 1.093068e-07 | 0.006585 | 0.000204 | ... | 14.946130 | 12.246561 | 1.537072 | 8.865320 | 20.141118 | 5.562198 | 14.933305 | 21.456669 | 8.898356 | 8.222981 |
| 98 | 0.0 | 0.0 | 14.507076 | -1.377036e-06 | 1.259269e-05 | 0.002890 | 0.0 | 1.129048e-07 | 0.000632 | 0.000175 | ... | 15.312657 | 12.335102 | 1.584222 | 8.660440 | 19.940211 | 5.690202 | 15.339280 | 21.301991 | 9.250233 | 8.631202 |
| 99 | 0.0 | 0.0 | 14.530041 | -1.425529e-06 | 1.304998e-05 | 0.002979 | 0.0 | 1.170060e-07 | 0.000446 | 0.000186 | ... | 15.015905 | 12.054240 | 1.472965 | 8.312579 | 19.900007 | 5.858699 | 15.042438 | 21.290107 | 9.207467 | 8.811348 |
100 rows × 2593 columns
[14]:
fluxes_sampling = reshape_fluxes_escher(sampled_fluxes)
[15]:
sampled_flux_map = escher.Builder('e_coli_core.Core metabolism',
reaction_data = fluxes_sampling).display_in_notebook()
sampled_flux_map
[15]:
There are also other ways to cosolidate the MFA calculations and the constraint based flux predictions. One of these is lMOMA (linear Minimization Of Metabolic Adjustment). MOMA assumes that the fluxes before and after adding the new constraints should be similar, so it aims to predicting an optimum for the newly constrined model while keeping the differences to the original model small. We suggest using pFBA (parsimonious Flux Balance Analysis) instead of regular FBA for this step as pFBA aims to keep the overal fluxes low.
Dealing with infeasible solutions - relaxation¶
Another way of dealing with infeasible models is to relax the added constraints to the point that it works again. You will need to have the Gurobi solver installed for this. The same principle is used in the BFAIR thermo tools. For that we add our constraints to a model that will now be infeasible. Have I meantioned that this is much more elegant, better and that you should do that? It is. The other method is fine but you exclude reactions and, in general, it is always better to use as
much as possible of the information that is available to you.
[16]:
model = model_original.copy()
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmprxly7pup.lp
Reading time = 0.03 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpm40nroxp.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
[17]:
model = add_constraints(model, fittedFluxes)
--- start ---
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpc0e4oi9r.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpdcw8mfmh.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Did not work for 26dap_DASH_MSYN
Did not work for ArgSYN
Did not work for EX_co2_e_unlabeled
Did not work for EX_glc_e
Did not work for FADR
Did not work for GluSYN
Did not work for HisSYN
Did not work for IleSYN
Did not work for LeuSYN
Did not work for MetSYN
Did not work for MlthfSYN
Did not work for MlthfSYN_reverse
Did not work for NADH
Did not work for PheSYN
Did not work for ProSYN
Did not work for SerSYN
Did not work for SUCCOAS
Did not work for ThrSYN
Did not work for TKT1a
Did not work for TKT1b
Did not work for TKT2a
Did not work for TKT2b
Did not work for TrpSYN
Did not work for TyrSYN
Did not work for ValSYN
Did not work for CYTBD
Did not work for HYD
Did not work for ATPS4r
add_constraints takes 0h: 0min: 11sec to run
--- end ---
[18]:
model.optimize()
/Users/matmat/opt/anaconda3/envs/bfair/lib/python3.8/site-packages/cobra/util/solver.py:430: UserWarning: solver status is 'infeasible'
warn("solver status is '{}'".format(status), UserWarning)
[18]:
Then we make use of the handy bound_relaxation() function that will test our model to figure out which of these added bounds need to be adjusted and return a DataFrame that describes the affected functions and the gravity of the suggested changes. If we allow this function to be desctructive it will adjust the input model right away.
[19]:
cons_table = bound_relaxation(model, fittedFluxes, destructive=True, fluxes_to_ignore=['BIOMASS_Ec_iJO1366_core_53p95M'])
cons_table
--- start ---
bound_relaxation takes 0h: 0min: 0sec to run
--- end ---
[19]:
| lb_change | ub_change | subsystem | |
|---|---|---|---|
| reaction | |||
| EX_o2_e | -16.610505 | 0.000000 | Extracellular exchange |
| ASPTA | -3.865062 | 0.000000 | Alanine and Aspartate Metabolism |
| GLNS | 0.000000 | 0.767431 | Glutamate Metabolism |
| EX_nh4_e_reverse_f9cc6 | 0.000000 | 12.360577 | Extracellular exchange |
| EX_so4_e_reverse_5c8ed | 0.000000 | 0.376547 | Extracellular exchange |
| DAPDC_reverse_d3ab8 | -0.040213 | 0.000000 | Threonine and Lysine Metabolism |
Let’s see if our model is feasible now.
[20]:
relaxed_solution = model.optimize()
relaxed_solution
[20]:
| fluxes | reduced_costs | |
|---|---|---|
| EX_cm_e | 0.000000 | 0.0 |
| EX_cmp_e | 0.000000 | -0.0 |
| EX_co2_e | 14.698232 | 0.0 |
| EX_cobalt2_e | -0.000017 | -0.0 |
| DM_4crsol_c | 0.000156 | 0.0 |
| ... | ... | ... |
| PTAr_reverse | 0.000000 | 0.0 |
| ACONTb_reverse | 0.000000 | 0.0 |
| PGK_reverse | 23.033030 | 0.0 |
| ACKr_reverse | 4.865707 | 0.0 |
| ACS_reverse | 0.000000 | 0.0 |
2593 rows × 2 columns
[21]:
relaxed_solution.fluxes["BIOMASS_Ec_iJO1366_core_53p95M"]
[21]:
0.7
And here’s the star of the show, our sampling method. We trust our models because… we have to! And because smart people that knew what they were doing set them up. So in order to gain more confidence in our MFA data, we sample the model after adding the calculated bound for some of the reactions and re-calculate the fluxes a number of time. Then, we take the mean and take that as the most trustworthy calculated flux. These fluxes can be visualized, for example in tools like Escher. So here
is the point where we can let the power of constraint based models work for us and make use of the tools mentioned above in order to adjust the fluxes calculated using MFA so that they nicely fit into our model.
[22]:
relaxed_sampled_fluxes = cobra.sampling.sample(model, n=100, processes=2)
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmpfjca0h7c.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
Read LP format model from file /var/folders/mb/7cs2dcbn369_w97fqkfx47brjcxl49/T/tmprv60cd5x.lp
Reading time = 0.02 seconds
: 1805 rows, 5186 columns, 20446 nonzeros
[23]:
relaxed_fluxes_sampling = reshape_fluxes_escher(relaxed_sampled_fluxes)
[24]:
relaxed_flux_map = escher.Builder('e_coli_core.Core metabolism',
reaction_data = relaxed_fluxes_sampling).display_in_notebook()
relaxed_flux_map
[24]:
And let’s have a quick look at the model itself
[25]:
model.summary()
[25]:
Objective
1.0 BIOMASS_Ec_iJO1366_core_53p95M = 0.7
Uptake
| Metabolite | Reaction | Flux | C-Number | C-Flux |
|---|---|---|---|---|
| ca2_e | EX_ca2_e | 0.003644 | 0 | 0.00% |
| cl_e | EX_cl_e | 0.003644 | 0 | 0.00% |
| cobalt2_e | EX_cobalt2_e | 1.75E-05 | 0 | 0.00% |
| cu2_e | EX_cu2_e | 0.0004963 | 0 | 0.00% |
| fe2_e | EX_fe2_e | 0.005777 | 0 | 0.00% |
| fe3_e | EX_fe3_e | 0.005466 | 0 | 0.00% |
| glc__D_e | EX_glc__D_e | 10 | 6 | 100.00% |
| k_e | EX_k_e | 0.1366 | 0 | 0.00% |
| mg2_e | EX_mg2_e | 0.006072 | 0 | 0.00% |
| mn2_e | EX_mn2_e | 0.0004837 | 0 | 0.00% |
| mobd_e | EX_mobd_e | 9.03E-05 | 0 | 0.00% |
| nh4_e | EX_nh4_e | 7.561 | 0 | 0.00% |
| ni2_e | EX_ni2_e | 0.0002261 | 0 | 0.00% |
| o2_e | EX_o2_e | 8.811 | 0 | 0.00% |
| pi_e | EX_pi_e | 0.6752 | 0 | 0.00% |
| so4_e | EX_so4_e | 0.1765 | 0 | 0.00% |
| zn2_e | EX_zn2_e | 0.0002387 | 0 | 0.00% |
Secretion
| Metabolite | Reaction | Flux | C-Number | C-Flux |
|---|---|---|---|---|
| 4crsol_c | DM_4crsol_c | -0.0001561 | 7 | 0.00% |
| 5drib_c | DM_5drib_c | -0.0001575 | 5 | 0.00% |
| amob_c | DM_amob_c | -1.4E-06 | 15 | 0.00% |
| mththf_c | DM_mththf_c | -0.0003136 | 5 | 0.01% |
| ac_e | EX_ac_e | -2.3 | 2 | 14.71% |
| co2_e | EX_co2_e | -12.7 | 1 | 40.61% |
| etoh_e | EX_etoh_e | -2.393 | 2 | 15.30% |
| h2o_e | EX_h2o_e | -28.79 | 0 | 0.00% |
| h_e | EX_h_e | -11.8 | 0 | 0.00% |
| lac__D_e | EX_lac__D_e | -3.061 | 3 | 29.37% |
| meoh_e | EX_meoh_e | -1.4E-06 | 1 | 0.00% |
Alternatives¶
There are also other ways to cosolidate the MFA calculations and the constraint based flux predictions.
MOMA¶
One of these is lMOMA (linear Minimization Of Metabolic Adjustment). MOMA assumes that the fluxes before and after adding the new constraints should be similar, so it aims to predicting an optimum for the newly constrined model while keeping the differences to the original model small. We suggest using pFBA (parsimonious Flux Balance Analysis) instead of regular FBA for this step as pFBA aims to keep the overal fluxes low.
From Volkova, 2020:
“While FBA with the objective of maximizing growth results in reasonable solutions for wild-type cells, it does not so for (unevolved) gene knockout mutants. While in principle wild-type cells optimize their growth, unevolved cells with knockouts do not. Since homeostasis governs metabolic reprogramming, we cannot assume that the cell will follow a common objective, such as maximizing its growth. To acknowledge the requirement for metabolic homeostasis, the minimization of metabolic adjustment (MOMA) approach was proposed. The main idea behind this method is that, to maintain homeostasis, the difference in fluxes before and after the perturbation should be minimal. MOMA predicts the fluxes of a knockout strain by assuming that the cell will have a minimal redistribution of fluxes compared to its ancestor.” We’re using pFBA istead of FBA because, on top of optimizing the growth rate, it also minimizes the total sum of fluxes.
[26]:
pfba_relaxed_solution = cobra.flux_analysis.pfba(model)
[27]:
moma_after_relaxed_MFA = cobra.flux_analysis.moma(
model=model, solution=pfba_relaxed_solution, linear=True)
[28]:
moma_after_relaxed_MFA.fluxes["BIOMASS_Ec_iJO1366_core_53p95M"]
[28]:
0.7
ROOM¶
From Volkova, 2020:
Another similar approach is the regulatory on/off minimization (ROOM), which minimizes the number of fluxes that are significantly different from the wild type. Wild-type fluxes, determined by FBA or other methods, need to be known in order to use these approaches.
[29]:
room_after_relaxed_MFA = cobra.flux_analysis.room(
model=model, solution=pfba_relaxed_solution, linear=True, delta=0.03, epsilon=0.001)
[30]:
room_after_relaxed_MFA.fluxes["BIOMASS_Ec_iJO1366_core_53p95M"]
[30]:
0.7
[31]:
from tabulate import tabulate
results = [['Before adding new constraints:', round(original_solution.objective_value, 2)],
['New constraints relaxed FBA:', round(relaxed_solution.objective_value, 2)],
['New constraints relaxed pFBA:', round(pfba_relaxed_solution.fluxes["BIOMASS_Ec_iJO1366_core_53p95M"], 2)],
['New constraints relaxed MOMA:', round(moma_after_relaxed_MFA.fluxes["BIOMASS_Ec_iJO1366_core_53p95M"], 2)],
['New constraints relaxed ROOM:', round(room_after_relaxed_MFA.fluxes["BIOMASS_Ec_iJO1366_core_53p95M"], 2)]]
print(tabulate(results, headers=["Method", "Biomass function"]))
Method Biomass function
------------------------------ ------------------
Before adding new constraints: 0.98
New constraints relaxed FBA: 0.7
New constraints relaxed pFBA: 0.7
New constraints relaxed MOMA: 0.7
New constraints relaxed ROOM: 0.7
Of course also the MOMA results can be visualized with Escher. The reshape_fluxes_escher() can take both pandas DataFrames or cobra solutions as an input.
[32]:
fluxes_relaxed_moma = reshape_fluxes_escher(moma_after_relaxed_MFA)
moma_relaxed_flux_map = escher.Builder('e_coli_core.Core metabolism',
reaction_data = fluxes_relaxed_moma).display_in_notebook()
moma_relaxed_flux_map
[32]:
[33]:
fluxes_relaxed_room = reshape_fluxes_escher(room_after_relaxed_MFA)
room_relaxed_flux_map = escher.Builder('e_coli_core.Core metabolism',
reaction_data = fluxes_relaxed_room).display_in_notebook()
room_relaxed_flux_map
[33]:
[ ]: