Trouble Translating Solution-Derived Constraints to Generic Model

[Rohak72]

Hi! Recently, I used MEWpy to run the GIMME algorithm, which returned me a condition-specific solution object with a shortened flux vector (comprised of just 3120 reactions as opposed to the 20677 contained in my GEM). To account for these changes, I created a copy of the original model, made a list of all model reactions NOT included in the GIMME model, and attempted to use remove_reactions() with remove_orphans = True to generate an updated, GIMME-specific model.

This updated version successfully includes all the desired reactions, but when I go to run model.optimize(), I see the same objective as the original model, even though the GIMME output (and solution) specified a different one. I'm a bit stumped as to why this is; I tried looking at the constraints of the GIMME solution and applied those to the model update process, but that yielded no change. Is there a method or inbuilt cobra functionality I'm missing that, given a solution object and dataframe/series of reaction fluxes, builds a specialized model? If not, how can I ensure that the change in the objective is reflected in my new model?

Thanks in advance, looking forward to hearing from you!

Best,
Rohak

[gregmedlock]

Can you share your code? Most likely, either you are missing a constraint in the new model or the objective is not the same as specified in GIMME.

It may also be helpful to know what your goal is. The GIMME solution is often used to create a pruned model like this to take an inventory on active reactions from a specific context and compare them with another, but what you're asking seems to be generating a pruned model and then re-applying the GIMME constraints/objective and optimizing (in which case the solution should be identical to that in the first model). What is your purpose in doing this?

[Rohak72]

Hi @gregmedlock,

Thanks for your response! To clarify, my goal in doing this is simply to convert the GIMME solution to a tissue-specific model. Currently, the GIMME output does return the pruned solution I'm looking for, but this output is not rendered in the form of a cobra model suitable for downstream pipelines (e.g. gMCS, fast-SL, gene KO simulations). The goal of my project is to screen for robust synthetic lethal pairs in dengue virus, so having a complete model would aid tremendously in the later analysis I plan to conduct. Does that make sense?

Below, I've included my code:

from cobra.io.sbml import read_sbml_model
from copy import deepcopy
from mewpy.omics import ExpressionSet
from mewpy.simulation import get_simulator
from mewpy.omics import eFlux, GIMME, iMAT

import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings("ignore")

# 1: reading in cobra model and expression data
model = read_sbml_model("Inputs/THG.xml")
raw = pd.read_csv('Inputs/ExpressionData.csv')[["EnsemblID", "Value1", "Value2"]]
raw.columns = ["EnsemblID", "CNTL", "DENV"]
raw = raw.groupby("EnsemblID").mean().reset_index()

# 2: preparing inputs to mewpy's gimme function (expr, sim)
genes = list(raw["EnsemblID"].values)
conditions = ["CNTL", "DENV"]
matrix = raw[["CNTL", "DENV"]].to_numpy()
expr = ExpressionSet(genes, conditions, matrix)
sim = get_simulator(model)

#3: running GIMME
gimme_res_UT = GIMME(sim, expr, biomass = "MAR13082", condition = 0, growth_frac = 0.9)

#4: updating model
gimme_UT = model.copy()
rxns_to_remove = list(set(gimme_UT.reactions.list_attr("id")) ^ set(gimme_res_UT.dataframe.index.values))
gimme_UT.remove_reactions(rxns_to_remove, True)

At the end of this sequence, both models retain the same objective value of 187.354. The objective function is also uniform across both models: 1.0*MAR13082 - 1.0*MAR13082_reverse_11d67.

Hope this provides some clarity, please let me know what you think.

Best,
Rohak

[cdiener]

Hi, sorry for the confusion. Right now remove_reactions will not touch the objective. So if you remove reactions that are in the objective they will remain there. However, they are now disconnected from the model and thus fixed to zero. I opened #1396 to discuss this. So you would have to update your objective yourself.

[Rohak72]

Ah, I see -- thanks for the clarification, @cdiener! Sorry for the basic question, but how exactly would I update the objective, considering that my intended objective expression would be the same across both models? Is there any way I could manually configure it? I tried temporarily changing the objective and then reverting back to MAR13082, but that didn't seem to make a difference.

[Rohak72]

Hello, just to follow-up: is it possible to update the objective in accordance to the reaction removal? Are there any inbuilt methods in place and if not, what steps could I take to move forward? I find that when I do model.summary() and pass in the solution values (fluxes, objective value, etc.), everything works properly. Could I potentially rehash this summary into a model of its own? Please let me know what you think, thanks!

[Rohak72]

Thanks for your response, @cdiener! So, my issue doesn't quite involve updating the objective expression -- it is supposed to remain the same in both the GIMME solution and the corresponding COBRA model. Yet, the main problem I'm facing is that the objective value of the solution compared to the objective value of the model (both containing the same constraints and reactions) are not the same.

As such, is there a way for the value of the objective to change in response to remove_reactions? I know you mentioned earlier that removed reactions will still remain in the model despite being disconnected, but I'm wondering if there's a way I could circumvent that. If that's not possible, I completely understand.

[cdiener]

Oh sorry, I definitely misunderstood then. So if two models have indeed the same reactions, constraints and objective the objective value has to be the same because it is unique (fluxes themselvels are not, but the objective is). So it would be like @gregmedlock said, either the models are different or a different objective was used. Do you know what the objective returned by Mewpy is? GIMME usually minimizes and inconsistency score, maybe that is the objective Mewpy returned?

[Rohak72]

No problem! MEWpy returned an objective of around 177.9 and its objective expression matched that of my original model: 1.0*MAR13082 - 1.0*MAR13082_reverse_11d67. I don't think -- at least from my knowledge -- that a different objective was used for the GIMME computation. What's really interesting, though, is when I iteratively run the knock_out() function on every reaction I want to remove (instead of just remove_reactions(list, True)), the objectives on the GIMME solution and context-specific model match up exactly. However, those knocked out reactions haven't been deleted from the model and if I try to do so, the problem re-appears.

That's sort of where I'm at right now. What do you think?

[cdiener]

Hmm, mysterious. Maybe related to the remove_orphans option removing metabolites that are in the biomass function. If you run it with remove_reactions(*, False) what objective value do you get?

[gregmedlock]

Can you share what the reaction MAR13082 is? Your objective for a model with GIMME constraints and objective will be minimization of each reaction's flux multiplied by the reaction's associated gene expression difference from a threshold. From the Becker 2008 paper:

It's possible MAR13082 is constructed to represent this summation, but pretty unlikely.

[Rohak72]

A mismatched objective turned out to be the issue! Thanks so much for letting me know, Greg -- I've fixed the objective issue internally and everything works great now. Thanks to the both of you (@cdiener, @gregmedlock) for your help, I really appreciate it!