biosustain / p-thermo

Genome-scale metabolic model of P. thermoglucosidasius NCIMB 11955
Apache License 2.0
5 stars 0 forks source link

Nucleotide triphosphate cycling #106

Closed surtfire closed 4 years ago

surtfire commented 4 years ago

The last change I want to implement at the moment deals with a few nucleotide bi/tri phosphate cycles involving 2.7.1.40 and 2.7.4.6 (see diagram), and presumably the redundant ones that are not currently selected too.

Ideally I want to set all 2.7.1.40 to only be able to go from ATP -> ADP, yet the bacillus model suggests that both sets of these reactions should be reversible. But I'm not sure how to break this cycle.

What sounds a good way of solving this?

image

Values below have O2 : baseline 30, 20, 10, 8, 6, 4, 3, 2, 0

image

image

vh-mol commented 4 years ago

yet the bacillus model suggests that both sets of these reactions should be reversible.

As far as I'm aware, the Bacillus model doesn't have proper reversibility set for a lot of the reactions. I've come across it before... Normally, reversibility is set based on thermodynamics, or specific information that is available (only for some reactions this is available).

The E. coli model (iML1515) has this checked better, and as thermodynamics isn't really organism specific we can also consider the reversibility from that model. For the 2.7.1.40 reaction, I just checked and (in both our model and the iML1515) the reaction is ADP + Phosphoenolpyruvate --> ATP + Pyruvate. This is also the final step of glycolysis, critical to the production of ATP, and so I think we should leave this as it is.

Based on thermodynamics, the 2.7.4.6 reaction should be reversible. I think the main thing that is weird is why we need to produce so much GTP when we decrease the amount of oxygen supplied? I think this is the main cause of the cycle: maybe we can look into why that is happening and then it will resolve the problem? As I think these two reactions themselves seem fine in themselves.

On another note, I'm a bit confused as to why your 'baseline' doesn't carry any flux through PYK? When I run the FBA, I get flux through PYK as 'i would expect (being a part of glycolysis)?

surtfire commented 4 years ago

'Ran it all again with the most up to date version of the model and looked at every use of GTP/GDP, CTP/CDP and more generally at all the nucletoide triphosphate variants. Once again the ATP PYK does not have any flux. Ultimately I think it all comes back to the nucleotide phosphate preference issue I mentioned when we discussed the NADK variants. Which all narrows down to the following: When any cap is applied to oxygen, the 2.7.1.40/2.7.4.6 cycle of R01858/7 shifts to R00430.

image

What’s driving me crazy is that I can’t find anything that produces enough GDP for the large negative flux of R00430. There does not appear to be any reaction which can spend the GTP at the rate it is being accumulated either, with R00330 showing no flux change,

image

Other than ATP, only the mono-phosphate nucleotides are used in the biomass equation, as they were what DNA and RNA were substituted for. As GTP/GDP are not used in biomass equation it appears their main current use is to cycle in this reaction couple. No other reactions which use GTP or CTP have flux above 0.5. Indeed, we'd assume from this that the model is actually accumulating GTP from R00430.

Thus given all of this, I think what you said the E.coli model is doing (Setting PYKs to operate only in the –TP -> -DP direction) is the best way of addressing this. Generally commentary on PYK suggests that “This reaction, although appearing reversible, is essentially irreversible under physiological conditions, thus helping control the metabolic flux in glycolysis”, so it is likely something we’ve missed that should be set unidirectional anyways! If this sounds ok, I’ll implement it.

vh-mol commented 4 years ago

Once again the ATP PYK does not have any flux. Ultimately I think it all comes back to the nucleotide phosphate preference issue I mentioned when we discussed the NADK variants.

From the data you show here this does seem to be the case. But again interestingly, when I try to run the optimization from my side (in the baseline) I get quite a decent flux through PYK. I also get a little flux through two other variants, but this is so much less that I don't think it is an issue. See below: image Also, when inspecting the baseline model, it only consumes 29.85 mmol/gcdw/h of oxygen. So in my case, capping at 30 mmol/gcdw/h doesn't change anything... So before i can really look and try to help I think we should figure out where this discrepancy is coming from? What command are you using for your optimization? Also, could you confirm what your medium is here? I have the following (which is also set as the model default): image

Also, when I decrease the oxygen supply (setting the lower bound of the EX_o2_ereaction) anywhere from -30 to -0.5 mmol/gcdw/h, the ratio between the R00430 and R00330 doesn't change: the flux predominantly goes through the ATP variant, and only slight amount through the GTP or CTP versions...

@BenjaSanchez any idea where this may be coming from?

Setting PYKs to operate only in the –TP -> -DP direction

And in the e. coli model they do not operate in this direction you mention here. The PYK reaction is generally used to produce ATP from PEP. In our model however, all the PYK reactions are already unidirectional, resulting in the formation of NTP and pyruvate from NDP and PEP. So this is already correct and doesn't need to be changed. Note that the 'reaction' column in your table doesn't reflect the reversibility of the reactions actually encoded in the model. (I'm not sure where they then come from, I don't know if this is matlab way of showing the reactions?)

So maybe if we figure out the reason why your analysis is giving such different values, we can see where to look further (or maybe the issue is then even solved, as I don't see it from my side?) In my analysis there is no GTP cycle, and the amount of GTP produced and consumed is quite small: mainly just to satisfy the carbohydrate demand in the biomass recation.: (in default settings) image Decreasing the oxygen supply doesn't increase the GTP demand either, it actually decreases it further...

surtfire commented 4 years ago

This is indeed very strange! Once again this morning after re-pulling everything I have PYK4 at -5.83426094832241 and PYK at 0

You're right about the PYK directionality, I think I've just been staring too long at the same problem.

The media is exactly the same too.

Currently a complete running for the baseline would be achieved with: initCobraToolbox setenv('ILOG_CPLEX_PATH','C:\Program Files\IBM\ILOG\CPLEX_Studio128\cplex\matlab\x64_win64'); Model = readCbModel('g-thermo.xml'); Model = changeObjective(Model,'biomass'); FBA = optimizeCbModel (Model,'max');

The only other thing I can think of that is different is that I use IBM ILOG CPLEX Studio v12.8 as the solver for the problem, but I've not noticed discrepancies this large in the past.

Either way, if there (kinda frustratingly) isn't a cycle here after all, then there isn't a problem with the model, and it's presumably something to do with my setup. I'll focus instead for the time being on picking out useful areas of our transcriptomics data for creating the aerobic and anaerobic models.

vh-mol commented 4 years ago

This is so weird.. hopefully Ben may have some insights here? Let's wait for his reply...

Then if there may not be a cycle in the end, I will re-run the GAM fitting and merge that into the model. Then I'll start trying to practice fitting the C13 data to the model. Did you manage to find Charlotte Wards data files on that?

BenjaSanchez commented 4 years ago

@surtfire I don't have much experience with the CPLEX solver, maybe Gurobi is more stable in this regard? In any case, if you are using optimizeCbModel I would always recommend using the minNorm option set as one, to avoid multiple optimal solutions. Perhaps switch to that setting and let us know if you see the cycle pop up again.

surtfire commented 4 years ago

I'd avoided that optional input in the past but setting it to one does solve the problem! Do you know why/know of any good papers that explain why taxicab geometry is more beneficial in these types of model?

Additionally, gurobi and CPLEX give the same solutions which is also good news.

BenjaSanchez commented 4 years ago

solving FBA is a lineal problem that is largely undetermined: in the g-thermo case, 1176 variables (reactions) - 890 equations (metabolites) = 286 degrees of freedom, so it's quite common to have coupled reactions that can have different combinations of fluxes. Adding a flux minimization step (can be taxicab or others) will always aid in this, as the most simple solution will be chosen and will get rid of most of this variability. You can learn more about alternative optima at https://doi.org/10.1016/j.ymben.2003.09.002

surtfire commented 4 years ago

Thanks Ben! I'll read that.