flux balance analysis and optimal growth without N input

Michel_Lavoie

unread,

Sep 11, 2018, 11:40:16 AM9/11/18

to COBRA Toolbox

Hi, I am running a flux balance analysis using a metabolic network of a phytoplanktonic algae. The objective function is simply a biomass function seeking to be maximized. I found that the network can grow without nitrogen inputs and I am not able to find why. In some way nitrogen is created in the network since when I run a FBA with a constraint on N (NH4) uptake rate (lb = -1000, ub = 0), I can find an optimal growth rate, but with a modeled N uptake rate equals to 0. There are no other N exchange reactions or N input in the model. However, when I try different (and negative) bounds on N exchange (lb = -1000 and ub = - 0.001), I cannot find an optimal growth rate and the solution is ‘suboptimal’ indicating that the model always find an optimal solution without taking up N from the medium. How can I resolve this and how can I find an optimal growth rate with a realistic significant N uptake rate (i.e., a negative bound for N exchange)?

Thanks for your help,

You can see below the 'SystemConfigReport'.

>> generateSystemConfigReport

> ---------------------------------- SYSTEM CONFIGURATION REPORT ----------------------------------

_____ _____ _____ _____ _____ |

/ ___| / _ \ | _ \ | _ \ / ___ \ | COnstraint-Based Reconstruction and Analysis

| | | | | | | |_| | | |_| | | |___| | | The COBRA Toolbox - 2018

| | | | | | | _ { | _ / | ___ | |

| |___ | |_| | | |_| | | | \ \ | | | | | Documentation:

\_____| \_____/ |_____/ |_| \_\ |_| |_| | http://opencobra.github.io/cobratoolbox

|

> Checking if git is installed ... Done.

> Checking if the repository is tracked using git ... Done.

> Checking if curl is installed ... Done.

> Checking if remote can be reached ... Done.

> Initializing and updating submodules (this may take a while)... Done.

> Adding all the files of The COBRA Toolbox ... Done.

> Define CB map output... set to svg.

> TranslateSBML is installed and working properly.

> Configuring solver environment variables ...

- [----] ILOG_CPLEX_PATH: --> set this path manually after installing the solver ( see instructions )

- [*---] GUROBI_PATH: /Library/gurobi752/mac64/matlab

- [----] TOMLAB_PATH: --> set this path manually after installing the solver ( see instructions )

- [----] MOSEK_PATH: --> set this path manually after installing the solver ( see instructions )

Done.

> Checking available solvers and solver interfaces ... Done.

> Setting default solvers ... Done.

> Saving the MATLAB path ... Done.

- The MATLAB path was saved in the default location.

> Summary of available solvers and solver interfaces

Support LP MILP QP MIQP NLP

----------------------------------------------------------------------

gurobi active 1 1 1 1 -

ibm_cplex active 0 0 0 - -

tomlab_cplex active 0 0 0 0 -

glpk active 1 1 - - -

mosek active 0 - 0 - -

matlab active 1 - - - 1

cplex_direct active 0 0 0 0 -

dqqMinos active 1 - - - -

pdco active 1 - 1 - -

quadMinos active 1 - - - -

qpng passive - - 1 - -

tomlab_snopt passive - - - - 0

gurobi_mex legacy 0 0 0 0 -

lindo_old legacy 0 - - - -

lindo_legacy legacy 0 - - - -

lp_solve legacy 1 - - - -

opti legacy 0 0 0 0 0

----------------------------------------------------------------------

Total - 7 2 3 1 1

+ Legend: - = not applicable, 0 = solver not compatible or not installed, 1 = solver installed.

> You can solve LP problems using: 'gurobi' - 'glpk' - 'matlab' - 'dqqMinos' - 'pdco' - 'quadMinos' - 'lp_solve'

> You can solve MILP problems using: 'gurobi' - 'glpk'

> You can solve QP problems using: 'gurobi' - 'pdco' - 'qpng'

> You can solve MIQP problems using: 'gurobi'

> You can solve NLP problems using: 'matlab'

> Checking for available updates ...

> There are 6429 new commit(s) on <master> and 0 new commit(s) on <develop> [c72a7b @ master]

> You can update The COBRA Toolbox by running updateCobraToolbox() (from within MATLAB).

Elapsed time is 143.438534 seconds.

-----------------------------------------------------------------------------------------------------

MATLAB Version: 9.4.0.813654 (R2018a)

MATLAB License Number: 310290

Operating System: Mac OS X Version: 10.13.6 Build: 17G65

Java Version: Java 1.8.0_144-b01 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode

-----------------------------------------------------------------------------------------------------

MATLAB Version 9.4 (R2018a)

Simulink Version 9.1 (R2018a)

Aerospace Toolbox Version 2.21 (R2018a)

Bioinformatics Toolbox Version 4.10 (R2018a)

Communications System Toolbox Version 6.6 (R2018a)

Computer Vision System Toolbox Version 8.1 (R2018a)

Control System Toolbox Version 10.4 (R2018a)

Curve Fitting Toolbox Version 3.5.7 (R2018a)

DSP System Toolbox Version 9.6 (R2018a)

Database Toolbox Version 8.1 (R2018a)

Embedded Coder Version 7.0 (R2018a)

Financial Toolbox Version 5.11 (R2018a)

Fixed-Point Designer Version 6.1 (R2018a)

Fuzzy Logic Toolbox Version 2.3.1 (R2018a)

Global Optimization Toolbox Version 3.4.4 (R2018a)

Image Acquisition Toolbox Version 5.4 (R2018a)

Image Processing Toolbox Version 10.2 (R2018a)

Instrument Control Toolbox Version 3.13 (R2018a)

MATLAB Coder Version 4.0 (R2018a)

MATLAB Compiler Version 6.6 (R2018a)

MATLAB Compiler SDK Version 6.5 (R2018a)

Mapping Toolbox Version 4.6 (R2018a)

Model Predictive Control Toolbox Version 6.1 (R2018a)

Neural Network Toolbox Version 11.1 (R2018a)

Optimization Toolbox Version 8.1 (R2018a)

Parallel Computing Toolbox Version 6.12 (R2018a)

Partial Differential Equation Toolbox Version 3.0 (R2018a)

RF Toolbox Version 3.4 (R2018a)

Robust Control Toolbox Version 6.4.1 (R2018a)

Signal Processing Toolbox Version 8.0 (R2018a)

Simscape Version 4.4 (R2018a)

Simscape Multibody Version 5.2 (R2018a)

Simscape Power Systems Version 6.9 (R2018a)

Simulink Coder Version 8.14 (R2018a)

Simulink Control Design Version 5.1 (R2018a)

Simulink Desktop Real-Time Version 5.6 (R2018a)

Statistics and Machine Learning Toolbox Version 11.3 (R2018a)

Symbolic Math Toolbox Version 8.1 (R2018a)

System Identification Toolbox Version 9.8 (R2018a)

Wavelet Toolbox Version 5.0 (R2018a)

> Default shell : /bin/bash

> Version of shell : GNU bash, version 3.2.57(1)-release (x86_64-apple-darwin17)

> Architecture : MACI64

> MATLAB folder : /Applications/MATLAB_R2018a.app

> COBRA Toolbox root : /Users/mlavoie/cobratoolbox

> git version : git version 2.15.2 (Apple Git-101.1)

> curl version : curl 7.54.0 (x86_64-apple-darwin17.0)

> CBT_LP_SOLVER : gurobi

> CBT_MILP_SOLVER : gurobi

> CBT_QP_SOLVER : gurobi

> CBT_MIQP_SOLVER : gurobi

> CBT_NLP_SOLVER : matlab

> GUROBI_PATH : /Library/gurobi752/mac64/matlab

> ILOG_CPLEX_PATH :

> TOMLAB_PATH :

> MOSEK_PATH :

> ----------------------------------- END OF CONFIGURATION REPORT -----------------------------------

> Please send the report located in

/Users/mlavoie/cobratoolbox/COBRAconfigReport.log

to the developers or post it in the forum:

Michel_Lavoie

unread,

Sep 12, 2018, 1:10:46 PM9/12/18

to COBRA Toolbox

Hi all, I forgot to tell you that I checked for matter production when all exchange reactions are set to 0 and when the objective function is the ATP consumption (ATP + H2O -> ADP + Pi + H) or an artificial NADPH consumption (NADPH -> NADP + H) reaction and I found no non-zero fluxes with a FBA minimizing the taxicab norm using LP (minNorm optional input = 'one' in the function OptimizeCbModel()). However, when I use the default 'inNorm' option (i.e., 0) in the OptimizeCbModel() function , I did get some non-zero fluxes. So, I am still wondering why the optimal FBA solution is always without N inputs (or a N uptake rate = 0 as specified below).

Thanks for your help,

Michel

Chintan Joshi

unread,

Sep 12, 2018, 2:45:55 PM9/12/18

to COBRA Toolbox

Hi Michael,

It seems that NH4 uptake is competing with the growth objective. Try running flux variability on both models (with and without NH4 constraints). It is possible that one or more of the biomass components are consumed amongst reactions responsible for NH4 assimilation in your model. Also in the model without NH4 constraints, what does the flux through exchange reactions look like, as in try to identify how the model is trying to meet its N needs in absence of NH4 uptake.

Thanks,

Chintan

Michel_Lavoie

unread,

Sep 12, 2018, 11:27:06 PM9/12/18

to COBRA Toolbox

Hi Chintan and all,

Without constraints on NH4+ exchange, the FBA optimal solution is a growth rate of 0.0095 h-1 through the biomass equation (see data and code below)

>> printConstraints(model1, 1000, -1000)

MinConstraints:

maxConstraints:

>> printConstraints(model1, -1000, 1000)

MinConstraints:

R_EX_co2_e -5.78704

R_ATPS3m 2

R_NGAM 1

maxConstraints:

R_EX_co2_e -5.78704

R_EX_photon_e -1000

R_ATPS3m 2

R_NGAM 1

>> printFluxVector(model1, FBAsolution.x, false, true)

R_EX_12ppd__S_e 0

R_EX_glc__D_e 0

R_EX_h2o_e 0

R_EX_h_e -5.273

R_EX_cl_e -0.00724

R_EX_pi_e -5.498

R_EX_nh4_e 0

R_EX_fe3_e -2.42e-05

R_EX_k_e -0.00543

R_EX_ca2_e -0.001524

R_EX_mg2_e -0.001898

R_EX_mn2_e -2.858e-06

R_EX_cobalt2_e -5.716e-07

R_EX_zn2_e -1.905e-06

R_EX_co2_e -5.787

R_EX_cu2_e -9.526e-07

R_EX_o2_e 7.836

R_EX_fe2_e 0

R_EX_h2s_e 0

R_EX_no3_e 0

R_EX_so4_e -0.004823

R_EX_photon_e -1000

R_EX_thm_e -2.267e-06

R_EX_btn_e -1.353e-05

EX_si[C_e] -0.007621

sink_octdp[C_c] 0

sink_dscl[C_c] 0

R_Biomass 0.009526

With constraints on NH4+ exchange, setting lower bound for NH4 exchange = -1000 and a upper bound = -1, the optimal FBA solution is a growth rate of 0.0029 h-1 (see data and codes below). So, the growth rate decreases when NH4+ uptake from the medium becomes significant. All exchanges reactions decreased by around 3.2 except oxygen exchange, which decreases by 30% and phosphorus exchange, which decreases by 50% and except the fixed CO2 and photon fluxes, which did not change with or without NH4+ uptake from the medium.

>> printConstraints(model1, -1000, 1000)

MinConstraints:

R_EX_co2_e -5.78704

R_ATPS3m 2

R_NGAM 1

maxConstraints:

R_EX_nh4_e -1

R_EX_co2_e -5.78704

R_EX_photon_e -1000

R_ATPS3m 2

R_NGAM 1

>> printFluxVector(model1, FBAsolution.x, false, true)

R_EX_12ppd__S_e 0

R_EX_glc__D_e 0

R_EX_h2o_e 0

R_EX_h_e -2.628

R_EX_cl_e -0.002236

R_EX_pi_e -3.698

R_EX_nh4_e -1

R_EX_fe3_e -7.472e-06

R_EX_k_e -0.001677

R_EX_ca2_e -0.0004707

R_EX_mg2_e -0.0005861

R_EX_mn2_e -8.825e-07

R_EX_cobalt2_e -1.765e-07

R_EX_zn2_e -5.883e-07

R_EX_co2_e -5.787

R_EX_cu2_e -2.942e-07

R_EX_o2_e 5.92

R_EX_fe2_e 0

R_EX_h2s_e 0

R_EX_no3_e 0

R_EX_so4_e -0.001489

R_EX_photon_e -1000

R_EX_thm_e -7.001e-07

R_EX_btn_e -4.177e-06

EX_si[C_e] -0.002353

sink_octdp[C_c] 0

sink_dscl[C_c] 0

R_Biomass 0.002942

I draw a metabolic map using the cobra toolbox paint4Net and the function ‘draw_by_met’.

[directionRxns, involvedMets, deadEnds] = draw_by_met(model1, {'nh4[C_c]'}, 'true', 1, 'struc', {''}, FBAsolution.x)

I found that, with NH4+ uptake, virtually all NH4+ comes from the external medium although threonine deaminase (degradation of the amino acid threonine to NH4+) also contribute slightly to intracellular NH4+ production. By contrast, without NH4+ uptake, virtually all NH4+ is produced by a serine deaminase (an enzyme degrading the amino acid serine to NH4+). So, as suggested by Chintan, NH4+ uptake appears to compete with biosynthesis of amino acids in the growth equation (i.e., it decreases the modeled growth rate). In the case where no NH4+ uptake occurs from the medium, the NH4+ produced inside the cell comes mostly from serine, which is synthesized from glycine and metabolites of the ADOMET cycle. The precursors of the ADOMET cycle are synthesized by the glyoxylate cycle (glycine, glyoxylate) and glycine in this cycle is produced from a coenzyme a (coa), which mostly comes from beta oxidation of fatty acids. So, I tried running the model without beta oxidation or glyoxylate cycle (i.e., glycolate cycle or photorespiration), and I found that the FBA cannot find an optimal growth rate (infeasible FBA). So, there are lots of pathways degrading amino acids (which forms the biomass equation) in my model, which can provide NH4+. There are also several reactions leading to beta oxidation of fatty acids, which produces coa and acetyl-Coa for the TCA cycle. No acetyl-Coa is produced from pyruvate dehydrogenase and thus glycolysis is not connected to the TCA cycle; beta oxidation is instead connected to the TCA cycle, which is strange for a biomass equation maximizing ATP production and synthesis of amino acids, fatty acids, pigments, sugars, vitamins (although my biomass equation does not include 'coa'). There is a pathway of fatty acid beta oxidation to propanoyl-Coa (ppcoa), which is converted to virtually all cellular coenzyme a (coa), which contains nitrogen and contributes to NH4+ production inside the cells from no nitrogen input from the medium. Indeed, ‘coa’ is not de novo synthesized in my molecular network and is not in the biomass equation, which means that ‘coa’ synthesis is not optimized to a given concentration in the cell. How could I tell the model that NH4+ uptake is essential and that ‘coa’ is just recycled in the cell and that the biosynthetic pathway of this ‘coa’ is neglected?

As requested by Chintan, I also ran Flux Variability Analysis and I found that NH4+ exchange vary between 0 and 2.12 mmol g DW-1 h-1 when there is no fixed NH4+ uptake (lb = -1000, ub = 0). When I fixed ub= -1, I find a minimum flux of -1.001 and a maximum flux of -1. Any idea on how to make the model take up NH4+ from the medium as a real cell would be very welcome!

Thanks for your help,

Michel

Thomas Pfau

unread,

Sep 13, 2018, 12:51:13 AM9/13/18

to cobra-...@googlegroups.com

My best guess would be that you have amino acid exchangers, or some mass balancing issues.

Did you try running checkMassChargeBalance and checking the results for inconsistencies?

Other than that, we would need the model to see what might be wrong.

Best

Thomas

On 09/12/2018 08:45 PM, Chintan Joshi wrote:

Unsubscribe

It appears that you have subscribed to commercial messages from this sender. To stop receiving such messages from this sender, please unsubscribey

--

---
You received this message because you are subscribed to the Google Groups "COBRA Toolbox" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cobra-toolbo...@googlegroups.com.
For more options, visit http://secure-web.cisco.com/1_YdD5RnfSkNvKIbElvwE69ridBGupv2b-764ZkQg25AbRYulV8nTIVsdyZu1xCxeutuKS-Lg-wOHxeJEv5B8gUmwswteXxuaNYAEPSRDBxCJOFFfUL6WaoJiBxqzL1PHg1SwrtWFt-EWJfs0hv3lnZxqqi1K4QLyOhBO7EN-kvX9Og3w-X78vUrvXMhlqGFwi76MrZ_VvjucXYGSGjnbEwFxrVhq5eLHomj7zXOJa7rkS5qhcb8R3vO32C9LMttN_ZA7TJQx2gXP2Aj5yav1lfKWCmFyAByII1Hr0O6vyEvuvJHHiDeNTVaQXg-6y7PFkDNalETuzfk9kZ6dS6jyOBhhJMJYSAdTRQGYzhCvz2vpSWBxghocLDKj62-DKrbChoLVSae7S-JXybcDMHgP8w/l34%3Ahttps%3A%2F%2Fgroups.google.com%2Fd%2Foptoute.

-- 
Université du Luxembourg
Faculté des Sciences, de la Technologie et de la Communication
Campus Belval, Biotech II 423
6 avenue du Swing
L-4367 Belvaux
Tel: (+352) 46 66 44 5309
Email: thoma...@uni.lu

Michel_Lavoie

unread,

Sep 13, 2018, 4:01:15 PM9/13/18

to COBRA Toolbox

Hi Thomas and all,

Since my model was almost extensively done directly from the Bigg database, I assumed that there were no problems with mass and charge balance. However, I just tried the 'checkMassChargeBalance' function and there are apparently several problems in mass and charge balance. They are probably mostly due to the fact that in the Bigg database, there might be two charges and/or two formulas for one metabolite in some cases. So, I will have to uniformize the number of formula and charge by selecting only one formula and one charge for each metabolite. We will see if there are errors remaining here. But, I am still wondering if a cofactor such as 'coenzyme a (coa)' , which is not de novo synthesized in the model, could provide nitrogen to the system? I suppose the toolbox consider the amount of N in each compound and takes that into account when performing FBA so that no matter is created from nothing when there are no mass or charge balance errors. Am I right?

I also tested the presence of inconsistencies in the model using the following line of code : 'FluxConsistent = verifyModel(model,'fluxConsistency', true)'. I found that there are several inconsistencies. I read that these inconsistencies can be explained by holes in the network (e.g., blocked reactions, closed exchange reactions). Is it the right and the only definition of 'network inconsistency'? I would like to make sure I understand that.

I will continue to look for explanations for growth in the absence of N inputs.

Thanks for your help,

Michel

To unsubscribe from this group and stop receiving emails from it, send an email to cobra-toolbox+unsubscribe@googlegroups.com.

For more options, visit http://secure-web.cisco.com/1_YdD5RnfSkNvKIbElvwE69ridBGupv2b-764ZkQg25AbRYulV8nTIVsdyZu1xCxeutuKS-Lg-wOHxeJEv5B8gUmwswteXxuaNYAEPSRDBxCJOFFfUL6WaoJiBxqzL1PHg1SwrtWFt-EWJfs0hv3lnZxqqi1K4QLyOhBO7EN-kvX9Og3w-X78vUrvXMhlqGFwi76MrZ_VvjucXYGSGjnbEwFxrVhq5eLHomj7zXOJa7rkS5qhcb8R3vO32C9LMttN_ZA7TJQx2gXP2Aj5yav1lfKWCmFyAByII1Hr0O6vyEvuvJHHiDeNTVaQXg-6y7PFkDNalETuzfk9kZ6dS6jyOBhhJMJYSAdTRQGYzhCvz2vpSWBxghocLDKj62-DKrbChoLVSae7S-JXybcDMHgP8w/l34%3Ahttps%3A%2F%2Fgroups.google.com%2Fd%2Foptoute.

Thomas Pfau

unread,

Sep 14, 2018, 1:18:16 AM9/14/18

to cobra-...@googlegroups.com

Hi Michel,

On 2018-09-13 22:01, 'Michel_Lavoie' via COBRA Toolbox wrote:

Since my model was almost extensively done directly from the Bigg database, I assumed that there were no problems with mass and charge balance. However, I just tried the 'checkMassChargeBalance' function and there are apparently several problems in mass and charge balance. They are probably mostly due to the fact that in the Bigg database, there might be two charges and/or two formulas for one metabolite in some cases. So, I will have to uniformize the number of formula and charge by selecting only one formula and one charge for each metabolite. We will see if there are errors remaining here.

This is a likely issue for balancing, in particular proton balancing, so yes, the chemical formulas as well as the reaction equations should be created with the same protonation state in mind.

But, I am still wondering if a cofactor such as 'coenzyme a (coa)' , which is not de novo synthesized in the model, could provide nitrogen to the system? I suppose the toolbox consider the amount of N in each compound and takes that into account when performing FBA so that no matter is created from nothing when there are no mass or charge balance errors. Am I right?

Without having the model at hand it is hard to tell.
The toolbox assumes, that a model is atom and charge balanced. And if it is not, any FBA results are by definition invalid. During FBA we do not check the balancing, as there should not be any imbalance due to the reactions, which is a requirement for a proper function of FBA.

I also tested the presence of inconsistencies in the model using the following line of code : 'FluxConsistent = verifyModel(model,'fluxConsistency', true)'. I found that there are several inconsistencies. I read that these inconsistencies can be explained by holes in the network (e.g., blocked reactions, closed exchange reactions). Is it the right and the only definition of 'network inconsistency'? I would like to make sure I understand that.

The fluxConsistency Check in verifyModel only tests, whether there are reactions in the model, which cannot carry flux (and would thus be inconsistent). In general, yes, these inconsistent reactions are likely due to gaps or potentially constraints in the model and in the end, any inconsistent reaction (defined in this way) could be made consistent by adding the respective exchange reactions for all its products and substrates.

Hi Thomas and all, --

---
You received this message because you are subscribed to the Google Groups "COBRA Toolbox" group.

To unsubscribe from this group and stop receiving emails from it, send an email to cobra-toolbo...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

-- 
Université du Luxembourg
Faculté des Sciences, de la Technologie et de la Communication
Campus Belval, Biotech II 115
6 avenue du Swing
L-4367 Belvaux
Tel: (+352) 46 66 44 5309
Email: thoma...@uni.lu

Ronan M.T. Fleming

unread,

Sep 14, 2018, 6:10:11 AM9/14/18

to COBRA Toolbox

https://opencobra.github.io/cobratoolbox/stable/tutorials/tutorialReconToFBAmodel.html

--

--
Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.
----------------------------------------------------------------------------
Assistant Professor,

Division of Systems Biomedicine and Pharmacology,

Leiden Academic Centre for Drug Research,

Faculty of Science,

Leiden University.

https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming

&

H2020 Project Coordinator

Systems Medicine of Mitochondrial Parkinson’s Disease

http://sysmedpd.eu

----------------------------------------------------------------------------
Mobile: +353 873 413 072
Skype: ronan.fleming
----------------------------------------------------------------------------
(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)

Michel_Lavoie

unread,

Sep 18, 2018, 12:05:55 PM9/18/18

to COBRA Toolbox

Hi all, I corrected all mass and charge balance problems and, now, I cannot find an optimal solution using FBA no matter what constraints are tried. The FBA solution is always unfeasible. I tried to add sink and demand reactions to all components of the biomass equation as explained in the tutorial (https://opencobra.github.io/cobratoolbox/stable/tutorials/tutorialReconstructionSOP.html). However, I find that all biomass components with negative coefficients cannot be removed from the model and that all biomass components with positive coefficients cannot be produced by the model. Is there a way to find out why the FBA solution is always infeasible? Would you have any idea why I get this result?

Thank you very much for your help,

Michel

https://opencobra.github.io/cobratoolbox/stable/tutorials/tutorialReconToFBAmodel.html

RF Toolbox Version 3.4<sp

Michel_Lavoie

unread,

Sep 18, 2018, 12:26:21 PM9/18/18

to COBRA Toolbox

Hi all again, sorry for my last email. I did a mistake since there was a constraints on C uptake when I performed the test for biomass component synthesis or removal. So, without constraints, except with a constraints on non-growth associated maintenance (NGAM), I can find an optimal growth rate equal to 0. All biomass components with negative coefficients can be removed from the model , but not all biomass components with positive coefficients could be produced by the model.

You may want to see a little piece of code and output below. So, What should I do to debug that? Thanks , Michel

>> for i = 1 : length(BiomassComponentsPos)

model1_NEW = changeObjective(model1_NEW, strcat('sink_', BiomassComponentsPos(i)));

FBAsolution = optimizeCbModel(model1_NEW, 'min');

BiomassComponentsValuePos(i, 1) = FBAsolution.f;

end

[BiomassComponentsPos num2cell(BiomassComponentsValuePos)]

ans =

5×2 cell array

{'adp[C_c]'} {[-1000]}

{'pi[C_c]' } {[-1000]}

{'h[C_c]' } {[-1000]}

{'ppi[C_c]'} {[-1000]}

{'gdp[C_c]'} {[-1000]}

>> for i = 1 : length(BiomassComponentsNeg)

model1_NEW = changeObjective(model1_NEW, strcat('DM_', BiomassComponentsNeg(i)));

FBAsolution = optimizeCbModel(model1_NEW, 'max');

BiomassComponentsValueNeg(i, 1) = FBAsolution.f;

end

[BiomassComponentsNeg num2cell(BiomassComponentsValueNeg)]

ans =

69×2 cell array

{'atp[C_c]' } {[ 1000]}

{'glu__L[C_c]' } {[ 0]}

{'h2o[C_c]' } {[ 1000]}

{'amp[C_c]' } {[ 1000]}

{'pg120[C_c]' } {[ 0]}

{'pg140[C_c]' } {[ 0]}

{'pg141[C_c]' } {[ 0]}

{'pg160[C_c]' } {[ 0]}

{'pg161[C_c]' } {[ 0]}

{'pg180[C_c]' } {[ 0]}

{'pg181[C_c]' } {[ 0]}

{'gtp[C_c]' } {[ 1000]}

{'leu__L[C_c]' } {[ 0]}

{'asp__L[C_c]' } {[ 0]}

{'ala__L[C_c]' } {[ 0]}

{'cl[C_c]' } {[ 1000]}

{'gln__L[C_c]' } {[ 0]}

{'cys__L[C_c]' } {[ 0]}

{'ser__L[C_c]' } {[ 0]}

{'thr__L[C_c]' } {[ 0]}

{'his__L[C_c]' } {[ 0]}

{'ump[C_c]' } {[ 0]}

{'dcmp[C_c]' } {[ 0]}

{'damp[C_c]' } {[264.0857]}

{'gmp[C_c]' } {[ 1000]}

{'cmp[C_c]' } {[ 0]}

{'gly[C_c]' } {[ 0]}

{'met__L[C_c]' } {[ 0]}

{'asn__L[C_c]' } {[ 0]}

{'pro__L[C_c]' } {[ 0]}

{'lys__L[C_c]' } {[ 0]}

{'fe2[C_c]' } {[ 1000]}

{'phe__L[C_c]' } {[ 0]}

{'so4[C_c]' } {[ 1000]}

{'k[C_c]' } {[ 1000]}

{'gthrd[C_c]' } {[ 0]}

{'ile__L[C_c]' } {[ 0]}

{'glyb[C_c]' } {[ 0]}

{'cu2[C_c]' } {[ 1000]}

{'ca2[C_c]' } {[ 1000]}

{'mg2[C_c]' } {[ 1000]}

{'mn2[C_c]' } {[ 1000]}

{'cobalt2[C_c]' } {[ 1000]}

{'zn2[C_c]' } {[ 1000]}

{'fe3[C_c]' } {[ 1000]}

{'dgmp[C_c]' } {[413.4405]}

{'dtmp[C_c]' } {[ 0]}

{'trp__L[C_c]' } {[ 0]}

{'val__L[C_c]' } {[ 0]}

{'tyr__L[C_c]' } {[ 0]}

{'glcur[C_c]' } {[ 0]}

{'pydxn[C_c]' } {[ 0]}

{'thm[C_c]' } {[ 1000]}

{'btn[C_c]' } {[ 1000]}

{'13glucan[C_c]' } {[ 0]}

{'si[C_c]' } {[ 1000]}

{'tag161[C_c]' } {[ 0]}

{'tag160[C_c]' } {[ 0]}

{'tag140[C_c]' } {[ 0]}

{'pc140[C_c]' } {[ 0]}

{'pc160[C_c]' } {[ 0]}

{'pc180[C_c]' } {[ 0]}

{'pc161[C_c]' } {[ 0]}

{'cholphya[C_c]' } {[ 0]}

{'cholphyc1[C_c]'} {[ 0]}

{'cholphyc2[C_c]'} {[ 0]}

{'caro[C_c]' } {[ 0]}

{'diadinx[C_c]' } {[ 0]}

{'fxanth[C_c]' } {[ 0]}

>> printConstraints(model1, -1000, 1000)

MinConstraints:

R_NGAM 1

maxConstraints:

R_NGAM 1

Michel_Lavoie

unread,

Sep 18, 2018, 2:03:16 PM9/18/18

to COBRA Toolbox

Hi all,

I did again a test of synthesis or removal of biomass components, but this time including glucose uptake for my algae. So, not only photosynthesis (light + CO2), but also glucose uptake can contribute to algal growth. Now , I find that all components of the biomass equation can be synthesized or removed (see code and results below). However, when I perform a FBA after that, the optimal growth rate is 0 when I use constraints on glucose uptake (lb = -1000, ub = 0). When I fixed the glucose uptake to -1000, I find an infeasible growth rate. So, I am still wondering what could be wrong here. Any help would be very welcome.

for i = 1 : length(BiomassComponentsPos)

model1_NEW = changeObjective(model1_NEW, strcat('sink_', BiomassComponentsPos(i)));

FBAsolution = optimizeCbModel(model1_NEW, 'min');

BiomassComponentsValuePos(i, 1) = FBAsolution.f;

end

[BiomassComponentsPos num2cell(BiomassComponentsValuePos)]

for i = 1 : length(BiomassComponentsNeg)

model1_NEW = changeObjective(model1_NEW, strcat('DM_', BiomassComponentsNeg(i)));

FBAsolution = optimizeCbModel(model1_NEW, 'max');

BiomassComponentsValueNeg(i, 1) = FBAsolution.f;

end

[BiomassComponentsNeg num2cell(BiomassComponentsValueNeg)]

BiomassComponentsPos =

5×1 cell array

{'adp[C_c]'}

{'pi[C_c]' }

{'h[C_c]' }

{'ppi[C_c]'}

{'gdp[C_c]'}

BiomassComponentsNeg =

69×1 cell array

{'atp[C_c]' }

{'glu__L[C_c]' }

{'h2o[C_c]' }

{'amp[C_c]' }

{'pg120[C_c]' }

{'pg140[C_c]' }

{'pg141[C_c]' }

{'pg160[C_c]' }

{'pg161[C_c]' }

{'pg180[C_c]' }

{'pg181[C_c]' }

{'gtp[C_c]' }

{'leu__L[C_c]' }

{'asp__L[C_c]' }

{'ala__L[C_c]' }

{'cl[C_c]' }

{'gln__L[C_c]' }

{'cys__L[C_c]' }

{'ser__L[C_c]' }

{'thr__L[C_c]' }

{'his__L[C_c]' }

{'ump[C_c]' }

{'dcmp[C_c]' }

{'damp[C_c]' }

{'gmp[C_c]' }

{'cmp[C_c]' }

{'gly[C_c]' }

{'met__L[C_c]' }

{'asn__L[C_c]' }

{'pro__L[C_c]' }

{'lys__L[C_c]' }

{'fe2[C_c]' }

{'phe__L[C_c]' }

{'so4[C_c]' }

{'k[C_c]' }

{'gthrd[C_c]' }

{'ile__L[C_c]' }

{'glyb[C_c]' }

{'cu2[C_c]' }

{'ca2[C_c]' }

{'mg2[C_c]' }

{'mn2[C_c]' }

{'cobalt2[C_c]' }

{'zn2[C_c]' }

{'fe3[C_c]' }

{'dgmp[C_c]' }

{'dtmp[C_c]' }

{'trp__L[C_c]' }

{'val__L[C_c]' }

{'tyr__L[C_c]' }

{'glcur[C_c]' }

{'pydxn[C_c]' }

{'thm[C_c]' }

{'btn[C_c]' }

{'13glucan[C_c]' }

{'si[C_c]' }

{'tag161[C_c]' }

{'tag160[C_c]' }

{'tag140[C_c]' }

{'pc140[C_c]' }

{'pc160[C_c]' }

{'pc180[C_c]' }

{'pc161[C_c]' }

{'cholphya[C_c]' }

{'cholphyc1[C_c]'}

{'cholphyc2[C_c]'}

{'caro[C_c]' }

{'diadinx[C_c]' }

{'fxanth[C_c]' }

sink_adp[C_c] adp[C_c] <=>

sink_pi[C_c] pi[C_c] <=>

sink_h[C_c] h[C_c] <=>

sink_ppi[C_c] ppi[C_c] <=>

sink_gdp[C_c] gdp[C_c] <=>

DM_atp[C_c] atp[C_c] ->

DM_glu__L[C_c] glu__L[C_c] ->

DM_h2o[C_c] h2o[C_c] ->

DM_amp[C_c] amp[C_c] ->

DM_pg120[C_c] pg120[C_c] ->

DM_pg140[C_c] pg140[C_c] ->

DM_pg141[C_c] pg141[C_c] ->

DM_pg160[C_c] pg160[C_c] ->

DM_pg161[C_c] pg161[C_c] ->

DM_pg180[C_c] pg180[C_c] ->

DM_pg181[C_c] pg181[C_c] ->

DM_gtp[C_c] gtp[C_c] ->

DM_leu__L[C_c] leu__L[C_c] ->

DM_asp__L[C_c] asp__L[C_c] ->

DM_ala__L[C_c] ala__L[C_c] ->

DM_cl[C_c] cl[C_c] ->

DM_gln__L[C_c] gln__L[C_c] ->

DM_cys__L[C_c] cys__L[C_c] ->

DM_ser__L[C_c] ser__L[C_c] ->

DM_thr__L[C_c] thr__L[C_c] ->

DM_his__L[C_c] his__L[C_c] ->

DM_ump[C_c] ump[C_c] ->

DM_dcmp[C_c] dcmp[C_c] ->

DM_damp[C_c] damp[C_c] ->

DM_gmp[C_c] gmp[C_c] ->

DM_cmp[C_c] cmp[C_c] ->

DM_gly[C_c] gly[C_c] ->

DM_met__L[C_c] met__L[C_c] ->

DM_asn__L[C_c] asn__L[C_c] ->

DM_pro__L[C_c] pro__L[C_c] ->

DM_lys__L[C_c] lys__L[C_c] ->

DM_fe2[C_c] fe2[C_c] ->

DM_phe__L[C_c] phe__L[C_c] ->

DM_so4[C_c] so4[C_c] ->

DM_k[C_c] k[C_c] ->

DM_gthrd[C_c] gthrd[C_c] ->

DM_ile__L[C_c] ile__L[C_c] ->

DM_glyb[C_c] glyb[C_c] ->

DM_cu2[C_c] cu2[C_c] ->

DM_ca2[C_c] ca2[C_c] ->

DM_mg2[C_c] mg2[C_c] ->

DM_mn2[C_c] mn2[C_c] ->

DM_cobalt2[C_c] cobalt2[C_c] ->

DM_zn2[C_c] zn2[C_c] ->

DM_fe3[C_c] fe3[C_c] ->

DM_dgmp[C_c] dgmp[C_c] ->

DM_dtmp[C_c] dtmp[C_c] ->

DM_trp__L[C_c] trp__L[C_c] ->

DM_val__L[C_c] val__L[C_c] ->

DM_tyr__L[C_c] tyr__L[C_c] ->

DM_glcur[C_c] glcur[C_c] ->

DM_pydxn[C_c] pydxn[C_c] ->

DM_thm[C_c] thm[C_c] ->

DM_btn[C_c] btn[C_c] ->

DM_13glucan[C_c] 13glucan[C_c] ->

DM_si[C_c] si[C_c] ->

DM_tag161[C_c] tag161[C_c] ->

DM_tag160[C_c] tag160[C_c] ->

DM_tag140[C_c] tag140[C_c] ->

DM_pc140[C_c] pc140[C_c] ->

DM_pc160[C_c] pc160[C_c] ->

DM_pc180[C_c] pc180[C_c] ->

DM_pc161[C_c] pc161[C_c] ->

DM_cholphya[C_c] cholphya[C_c] ->

DM_cholphyc1[C_c] cholphyc1[C_c] ->

DM_cholphyc2[C_c] cholphyc2[C_c] ->

DM_caro[C_c] caro[C_c] ->

DM_diadinx[C_c] diadinx[C_c] ->

DM_fxanth[C_c] fxanth[C_c] ->

ans =

5×2 cell array

{'adp[C_c]'} {[-1000]}

{'pi[C_c]' } {[-1000]}

{'h[C_c]' } {[-1000]}

{'ppi[C_c]'} {[-1000]}

{'gdp[C_c]'} {[-1000]}

ans =

69×2 cell array

{'atp[C_c]' } {[ 1000]}

{'glu__L[C_c]' } {[473.2143]}

{'h2o[C_c]' } {[ 1000]}

{'amp[C_c]' } {[ 1000]}

{'pg120[C_c]' } {[ 74.2630]}

{'pg140[C_c]' } {[ 67.1465]}

{'pg141[C_c]' } {[ 66.4835]}

{'pg160[C_c]' } {[ 59.2825]}

{'pg161[C_c]' } {[ 58.9932]}

{'pg180[C_c]' } {[ 53.0301]}

{'pg181[C_c]' } {[ 51.8855]}

{'gtp[C_c]' } {[ 1000]}

{'leu__L[C_c]' } {[258.0438]}

{'asp__L[C_c]' } {[ 1000]}

{'ala__L[C_c]' } {[ 1000]}

{'cl[C_c]' } {[ 1000]}

{'gln__L[C_c]' } {[352.0408]}

{'cys__L[C_c]' } {[370.4319]}

{'ser__L[C_c]' } {[ 1000]}

{'thr__L[C_c]' } {[718.2540]}

{'his__L[C_c]' } {[209.1157]}

{'ump[C_c]' } {[ 125]}

{'dcmp[C_c]' } {[122.2216]}

{'damp[C_c]' } {[606.0268]}

{'gmp[C_c]' } {[ 1000]}

{'cmp[C_c]' } {[116.0026]}

{'gly[C_c]' } {[603.3333]}

{'met__L[C_c]' } {[206.8646]}

{'asn__L[C_c]' } {[ 500]}

{'pro__L[C_c]' } {[498.1203]}

{'lys__L[C_c]' } {[ 500]}

{'fe2[C_c]' } {[ 1000]}

{'phe__L[C_c]' } {[337.4461]}

{'so4[C_c]' } {[ 1000]}

{'k[C_c]' } {[ 1000]}

{'gthrd[C_c]' } {[186.9409]}

{'ile__L[C_c]' } {[378.8820]}

{'glyb[C_c]' } {[198.9011]}

{'cu2[C_c]' } {[ 1000]}

{'ca2[C_c]' } {[ 1000]}

{'mg2[C_c]' } {[ 1000]}

{'mn2[C_c]' } {[ 1000]}

{'cobalt2[C_c]' } {[ 1000]}

{'zn2[C_c]' } {[ 1000]}

{'fe3[C_c]' } {[ 1000]}

{'dgmp[C_c]' } {[ 1000]}

{'dtmp[C_c]' } {[120.7864]}

{'trp__L[C_c]' } {[175.7635]}

{'val__L[C_c]' } {[527.8293]}

{'tyr__L[C_c]' } {[364.2499]}

{'glcur[C_c]' } {[597.6190]}

{'pydxn[C_c]' } {[361.1944]}

{'thm[C_c]' } {[ 1000]}

{'btn[C_c]' } {[ 1000]}

{'13glucan[C_c]' } {[443.1217]}

{'si[C_c]' } {[ 1000]}

{'tag161[C_c]' } {[ 40.7600]}

{'tag160[C_c]' } {[ 41.0522]}

{'tag140[C_c]' } {[ 44.5973]}

{'pc140[C_c]' } {[ 53.3848]}

{'pc160[C_c]' } {[ 48.3326]}

{'pc180[C_c]' } {[ 44.1540]}

{'pc161[C_c]' } {[ 48.1319]}

{'cholphya[C_c]' } {[ 19.8506]}

{'cholphyc1[C_c]'} {[ 21.3156]}

{'cholphyc2[C_c]'} {[ 21.2474]}

{'caro[C_c]' } {[ 38.6213]}

{'diadinx[C_c]' } {[ 20.8333]}

{'fxanth[C_c]' } {[ 17.8571]}

>> FBAsolution = optimizeCbModel(model1,'max', 'one')

Ronan M.T. Fleming

unread,

Sep 18, 2018, 2:14:36 PM9/18/18

to COBRA Toolbox

Hi Michael,

try relaxedFBA. See how here

https://opencobra.github.io/cobratoolbox/stable/tutorials/tutorialRelaxedFBABounds.html

and also

https://opencobra.github.io/cobratoolbox/stable/tutorials/tutorialRelaxedFBA.html

Regards,

Ronan

--

---
You received this message because you are subscribed to the Google Groups "COBRA Toolbox" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cobra-toolbo...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Michel_Lavoie

unread,

Sep 18, 2018, 11:07:49 PM9/18/18

to COBRA Toolbox

Hi Ronan and all,

Thanks for your your help and the web links. In the presence of glucose , i can now find a reasonable growth rate when using reasonable constraints. However, when I try to model the growth rate in the light in absence of glucose, so the only source of carbon is CO2, I cannot find a way to synthesize all components of the biomass equation and I find an optimal growth rate equals to 0. The FBA returns an optimal result (for most constraints I tried) rather than an infeasible results. So, I could not use a relaxed FBA in this case because the FBA were not infeasible although I did tried for a case where the FBA was infeasible and it did not improve the results. I also ran the function (biomassPrecursorCheck(model1)), it appears that several components of the biomass equation cannot be synthesized even though there are no deadend in my model and all equations are mass and charge balanced. However, when I run the function ' verifyModel(model1,'fluxConsistency', true) ' , I find that several reactions (around 250 out of 800) are inconsistents in my model. How can I find a way the model can properly synthesize all components of the biomass equations and thus calculate a non-zero growth rate?