FVA objective cap

[kamoors]

Hi, me again..

We are using iMAT++'s integrated FVA function, while also using the CS models separately. While doing FVA on those models using COBRA we were puzzled as to how iMAT++ does FVA.

In the FluxVariability() function from COBRA one can set the % of the objective (e.g. biomass), but I don't see anything like that in the parameters.

Could you please clarify?

Thank you!

Karlis

[XuhangLi]

Hi Karlis,

Thank you very much for your interest in our FVA function. I am excited to see in-depth user like you for our algorithm!

Yes, it is a great question – I should clarify about the CS models. As you may see in the code description, iMAT++ produced a Mixed-Integer Linear Problem Model (MILP) instead of regular LP. This is because iMAT++ is designed to work with integer constraints such as the gene on/off status. By saying that, it means the MILP model cannot be analyzed by regular FVA algorithms designed for LP, for example, the FluxVariability function in CORBA. I am not sure if you will see an error if you run on FluxVariability().

Our FVA_MILP() function is a very simple wrapper function: it just runs the maximization and minimization of each query reaction as objective function given a pre-constrained MILP. We had to make this function because MILP is not compatible with regular FVA functions in COBRA.

For your purpose, it looks like you wanted to set an additional objective with its constraint when doing FVA. This can be easily achieved with some adjustments before using the FVA_MILP() function in our package.

Say if you want to constrain the reaction ‘RXN0001’ in your model. You have the original model as ‘model’ and the iMAT++ MILP as ‘MILP’. You wanted to constrain the RXN0001 to 90% of its maximum. Just do this:

[myObj_lb, myObj_ub] = FVA_MILP(MILP, model,{ ‘RXN0001’}); % obtain the max of the objective function RXN0001 MILP.lb(strcmp(model.rxns, ‘RXN0001’)) = myObj_ub * 0.9; % constrain the minimal flux of the objective to 90% of its maximum [lb, ub] = FVA_MILP(MILP, model, model.rxns); % do FVA for all reactions

Unfortunately, I don’t have any simple idea to make the MILP model compatible with regular COBRA functionalities that are for LP models. Hope the above workaround can help with your research!

Best, Hang

From: kamoors @.> Date: Friday, April 28, 2023 at 7:09 AM To: WalhoutLab/MERGE @.> Cc: Subscribed @.***> Subject: [WalhoutLab/MERGE] FVA objective cap (Issue #7)

Hi, me again..

We are using iMAT++'s integrated FVA function, while also using the CS models separately. While doing FVA on those models using COBRA we were puzzled as to how iMAT++ does FVA.

In the FluxVariability() function from COBRA one can set the % of the objective (e.g. biomass), but I don't see anything like that in the parameters.

Could you please clarify?

Thank you!

Karlis

— Reply to this email directly, view it on GitHubhttps://github.com/WalhoutLab/MERGE/issues/7, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AJ6DQPJ33QKGLJPPTRCDY3TXDOQO3ANCNFSM6AAAAAAXPC2AIM. You are receiving this because you are subscribed to this thread.Message ID: @.***>

[kamoors]

Hi!

Thank you for the quick response.

Why I am asking is that when I run FVA with your algorithm on the MILP output, or if I obtain the active reactions and then overlay that with the base model, I get different values for the same reaction, as well as different ratios between samples.

So, e.g. with MILP ratio of max flux for model_sample1 / model_sample2 = 1.18, whereas with the active rxns overlayed I get model_sample1 / model_sample2 = 2

Since we want to obtain a COBRA model at the end of this pipeline instead of a MILP problem it is important to not have a discrepancy... Do you know why this could be happening?

Also, this is how I obtain the ''active reactions'' (ie. 1 and 0 for each reactions). Please let me know if you think it's a problem:

is_active = levels_f + levels_r; is_active(is_active > -2) = 1; is_active(is_active == -2) = 0;

So, basically if the reaction can be active in OFD or AFD then I set it to 1, of it's -1 in levels_f and levels_r then I set it to 0.

Then I take all reactions in the base iCEL1314 model that == 0 and set lb and ub = 0.

Karlis

[XuhangLi]

Hi Karlis,

Yes, you are right – you will get different values by MILP FVA as compared with the overlayed base model. It is because the MILP model is more sophisticatedly constrained than the overlayed model you mentioned.

You did the overlayed model very correctly: it is similar to the tissue-specific models we used as the FPA input in the paper (FPA is a simple LP algorithm). The advantage is as what you are looking for – it is still a COBRA model; the disadvantage is that by doing such active/inactive constraining, we lose the flexibility in the MILP.

In MILP, we control things with integer constraints that force flux on ‘on’ genes and block flux on ‘off’ genes. The biggest difference is that the MILP forces the flux through ‘on’ genes; but when you convert it into overlayed LP, the ‘on’ gene control is lost. There are other more detailed differences in MILP, but let’s simplify it for now.

I think you are totally fine to make a context-specific COBRA model by simply blocking the inactive reactions, like what you said. If I could recommend, I would further recommend doing so with reaction directionality: if the levels_f = -1 and levels_r = 1, I would set the upper boundary of the reaction to 0 but leave the lower boundary free. (in your proposal, you leave both upper and lower boundaries free). This may give you more context in your resulting model. If you do so, it is identical to the tissue-specific models we derived in our paper for downstream FPA analysis. I believe there are some example codes to do that in the repo.

Hope it helps! Thanks

Hang

From: kamoors @.> Date: Friday, April 28, 2023 at 11:25 AM To: WalhoutLab/MERGE @.> Cc: Li, Xuhang @.>, Comment @.> Subject: Re: [WalhoutLab/MERGE] FVA objective cap (Issue #7)

Hi!

Thank you for the quick response.

Why I am asking is that when I run FVA with your algorithm on the MILP output, or if I obtain the active reactions and then overlay that with the base model, I get different values for the same reaction, as well as different ratios between samples.

So, e.g. with MILP ratio of max flux for model_sample1 / model_sample2 = 1.18, whereas with the active rxns overlayed I get model_sample1 / model_sample2 = 2

Since we want to obtain a COBRA model at the end of this pipeline instead of a MILP problem it is important to not have a discrepancy... Do you know why this could be happening?

Also, this is how I obtain the ''active reactions'' (ie. 1 and 0 for each reactions). Please let me know if you think it's a problem:

is_active = levels_f + levels_r; is_active(is_active > -2) = 1; is_active(is_active == -2) = 0;

So, basically if the reaction can be active in OFD or AFD then I set it to 1, of it's -1 in levels_f and levels_r then I set it to 0.

Then I take all reactions in the base iCEL1314 model that == 0 and set lb and ub = 0.

Karlis

— Reply to this email directly, view it on GitHubhttps://github.com/WalhoutLab/MERGE/issues/7#issuecomment-1527731711, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AJ6DQPJVHZYYKT2UJ6UUTBDXDPONBANCNFSM6AAAAAAXPC2AIM. You are receiving this because you commented.Message ID: @.***>