Blocked reactions with High Confidence Score

[migp11]

Hi,

I'm using CORDA to generate context-specific model based Recon2.2.

After generating a my CS-model, I found that its contains blocked reaction with High Confidence Score
which are not blocked in Recon2.2. This looks a bit buggi to me. I'm reading the original paper on corda to get a proper insight on the algorithm. Then I'll give a look at the code.

Hope to help
Best
Miguel

[cdiener]

After optimization do those appear in the *.impossible list of your CORDA object. In order to guarantee numerical stability I currently require that any high confidence reeaction can at least carry a flux of 1 when all reaction bounds are raised to (-1e-6, 1e6). Since the default bounds of Recon 2 are in the (-1000, 1000) this basically means it will consider everything as blocked that can not carry a flux of at least 1e-3. Maybe it makes sense to rather use a scaling of the bounds than those absolute values... For now you can control this using the *.tflux attribute of the CORDA object. for instance with:

opt = CORDA(...)
opt.tflux = 1e-3
opt.build()

Recon 2.2 has a lot of blocked reactions (>2150) even when opening all inputs (which is weird). How did you evaluate which reactions you consider blocked in Recon2?

[migp11]

Hi @cdiener

Thanks for your quick answers (both of them) and for this great implementation. I don't know why most researcher keep developing novel cobra algorithm for Matlab :-/ is a pity.

Respect to the blocked reactions issue I've done the following.

1- I detect blocked reactions in Recon2.2, opening every exchange flux (in both directions) and performing an FVA (in fact, I did that just by calling the function find_blocked_reactions). As you already mentioned, Recon2.2 has many blocked reactions. In the 2.2 version (https://link.springer.com/article/10.1007/s11306-016-1051-4) I've found exactly 1921, which is around ~24% (*) and I'll refer to this set as meta-blocked.

2- I ran CORDA using Recon2.2 and passing reactions confidence computed from gene expression confidence for a set of cancer cell lines. The resulting CS-models exhibit blocked reaction, most of them, belonging to the set of meta-blocked (those already blocked in Recon2.2) but some of them (76) which are not blocked in Recon2.2. Of course, to find the blocked reactions I used the same procedure than in 1.

If you think it will help, I can send you the model (Recon2.2), a table with the reaction confidence for one of the models, as well a notebook, so you can reproduce this thing).

Again, thank you for your feedback and for this contribution!

Best,
Miguel

(*) I think is a reasonable number considering that, as a references, E. coli iJO1366 exhibit ~16% of blocked reactions and mammalian cells metabolism are in principle much complex.

[cdiener]

Yes if you want to you can send me the confidences. However if you run it yourself and just paste the output of the "impossible" list that would be enough. I am pretty sure it's related to the target flux cutoff and am thinking about a solution that is transparent but also numerically consistent.

[migp11]

Hi,

Sorry for the delay in my answer.. busy weeks :-/
I've checked again my previous calculation and I've found that the CSM produced by corda include 77 blocked reactions, which are not blocked in the Recon2.2

I did another experiment, I've remove the metablocked retacions from Recon2.2 an re run corda to create the CSM. the produced model also contains the same 77 blocked reactions.

Attached you will find the jupyter notebook as well the data I've used to do the experiment.
Hope it helps.

Best
Miguel

corda_test.zip

[cdiener]

Hi, sorry for the delay. So the short answer is numerical accuracy and here comes the longer one.

The conditions to test for blocked reactions that corda uses are pretty different from what find_blocked_reactions uses. For corda all bounds are raised to (-1e6, 1e6) (respecting reversibility of course) and that also applies for exchange reactions. After that corda considers all reactions with a flux larger 1e-6 as non-blocked. The problem is that in some rare cases (depends on numerical stability of the universal model) you actually need fluxes smaller than 1e-6 in some reactions to give
flux to the target reaction. However that can never be evaluated in a numerically stable manner so corda will no include those reactions. This can lead to some reactions that are blocked, however those reactions depend on flux values of other reaction in an instable manner (those other reactions might strictly violate bounds or equality constraints).
For the future I think it makes sense to switch to a scaling based on your initial boundaries. That should make those things less common since the used conditions in the model construction will be closer to the initial bounds. I hope I can get to that soon...

[migp11]

Hi,

Thank you for your replay. Indeed, I also thought it was a problem of numerical stability. I think the main problem relies on the fact the the biomass equation contains stoichiometric coefficients that span across (at least) six orders of magnitude and this causes a flux distribution producing biomass to have fluxes spanning across this same range. This range is very close to the numerical precision of most of the solvers and thus is quite difficult to define a general cut-off. In the case of CORDA, I think the main point is that reactions included in the high confidence core (confidence=3) should be active in the CSM unless these reactions were already blocked in the global (or reference) model. One possible walk around could be the following:

1. Check blocked reactions in the global model.
2. Prune blocked reaction independently of its confidence (may throw some warning)
3. Create the CSM and check that there are no blocked reactions.

Hope it helps
Best,
Miguel

[cdiener]

Yeah, might be good to mention in the docs that the input model should not contain blocked reactions. For newer models such as Recon3 this is probably not necessary since it already comes with a version that has no blocked reactions...

[cdiener]

I added a new rescaling strategy in #14 that is more representative for the original model. I hope that should give better coherence between the reconstruction and the base model in terms of blocked reactions. Still need to check whether that is true but feel free to try it out.