accessing the corda solution

[resendislab]

accessing the corda solution

[cdiener]

So more or less I would say. CORDA runs several iterations of that minimization. The lowest level access to a single iteration is through the CORDA.associated function that identifies the required reactions to activate a target. For instance using the example model:

In [1]: from corda import test_model
   ...:
   ...: mod = test_model()
   ...:

In [2]: conf = {}
   ...: for r in mod.reactions: conf[r.id] = -1
   ...: conf["r60"] = 3
   ...:

In [3]: from corda import CORDA
   ...:
   ...: opt = CORDA(mod, conf)
   ...:

In [4]: opt.associated(["r60"])
Out[4]:
array(['r10_reverse_a2a1f', 'r12', 'r14', 'r15', 'r20',
       'r21_reverse_3489c', 'r22_reverse_ad87b', 'r23_reverse_efabc',
       'r24_reverse_c4377', 'r25_reverse_049e4', 'r26_reverse_798ae',
       'r29_reverse_924fc', 'r30_reverse_cbe40', 'r31',
       'r35_reverse_a9257', 'r36_reverse_eb2c2', 'r38', 'r39', 'r4', 'r40',
       'r41', 'r42', 'r43', 'r46', 'r47', 'r50', 'r51_reverse_ce7c7',
       'r54', 'r59_reverse_b74cc', 'r6', 'r7', 'r8'],
      dtype='<U17')

If you use GLPK you can cheat and get the solution from that intermediate result since GLPK does not invalidate the primal values when the objective is reset (cplex and gurobi do so this would raise an error).

In [5]: from cobra.core import get_solution

In [6]: sol = get_solution(opt.model)

In [7]: sol
Out[7]: <Solution 2884.676 at 0x7f3ee5897080>

Here 2884 is the associated cost. However, there is no way to get the actual optimization criterion. So if that feature is interesting to you it might be better to implement it... Do you have a particular use case in mind? That would help in selecting the correct data structure to log those... 😄

[jotech]

Hi @cdiener sorry for my late reply! Thank you very much for your answer!

The reason for having a look at the flux distribution of corda is that I'm wondering if it couldn't be used as useful flux distribution which is close to the expression data. So I'm interested in the flux of reactions shown under minimization of the cost-consuming reaction.

I was playing a little bit with your examples and put my code in a notebook.
While it seems to work with the ecoli core model, I encountered problems when trying my luck with the recon2 model. Corda can deliver a reduced model but the corresponding fluxes are empty?!
Here I'm somehow confused. Doesn't corda including the active reactions in the reduced model?
Do you have a clue why it isn't working?

I'm also wondering how corda is handling the exchange reactions. Are the bounds of the reactions relaxed or could the stay constrained? (e.g. if no oxygen is available, then the optimization will start with a lower bound for oxygen of 0)

[cdiener]

Hi, thanks for including a notebook. That helped understanding what you are trying to do. So the CORDA algorithm has three stages and the results you are seeing are the ones from the last stage (add free reactions). The important thing here is that CORDA iterates by reaction. So it tries to find the cheapest way to activate a single reaction. In the Recon case you are probably seeing a random result for an import flux which requires no other reactions than itself. That output is probably not what you want. If you would want fluxes associated with a single reaction activation this would have to be done with a fresh CORDA object and using the associated function like shown above.

However, I would say the best flux distribution for your data is the pFBA result of the full reconstruction. Which you can extract with opt.cobra_model().

from cobra.flux_analysis import pFBA
from corda import CORDA

opt = CORDA(model, conf)
opt.build()
reconstruction = opt.cobra_model()
sol = pFBA(reconstruction)

However, that may still give you a solution with many zero fluxes since not all are required for growth.

In general, it will be hard to use CORDA in the same way as method that get fluxes directly from gene expression measurements since the philosophy is different. CORDA does not assume that all reactions with active gene-reaction rules have to be able to carry flux at the same time, however all of those must be able to carry flux. This is based on experience with bacterial cells that usually express many more enzymes than is actually required for maximum growth to provide some bet hedging. If you want to change this assumption you would have to change you objective aver the reconstruction. For instance by maximizing the sum of absolute fluxes or something like that which would give you a flux value for each reaction.

[cdiener]

Regarding the bound: the upper bound of the irreversible reactions are not respected due to numerical issues. Those will be raised to a very high value in order to disallow zero-flux reactions to be included in the solution. However, they will bes respected if they equal zero (reversibility or knock outs are respected). The lower bounds of the irreversible reactions will be respected. So forcing a flux through a reaction will work. That behavior is different from the Matlab version which does not respect bounds in that way.

[jotech]

Thanks for the clarification! Didn't know that corda iterates per reaction!

Concerning the associated function, do I get only associated reaction or also the flux of these reactions?
Using pFBA is a good idea, I have tried this. But in case of the human model it still maximize biomass which is not a justified objective. I was thinking about having cost associated fluxes for each reaction and try to use them instead. Or minimizing the distance to these fluxes using moma. Perhaps the later might work even in case of corda's iterating per reaction?

[cdiener]

Yeah that is currently difficult due to the redundancy detection since CORDA also checks for alternative pathways with the same cost. I could log the fluxes for each pathway and each reaction, however those fluxes are still obtained with different upper/lower bounds since for numerical stability raises the outer bounds to 1e6 and then tries to obtain a minimum flux of 1 through the target reaction. It would not be too hard to write a version of the associated function that does what you want. Maybe raising the solver tolerance as well... How would you deal with reversible reactions? For instance if a reversible reactions has a confidence of 3 would you want to test both directions?

[cdiener]

Hi, will be closing this for now. Feel free to reopen if you don't agree.