Looping through Model Bounds for Different Conditions

[AmalReji]

I am looping through bounds for different conditions and finding the solution for each:

input:

model = cobra.io.sbml.read_sbml_model(“file.xml”)
conditions =[“condition_1”, “condition_2”]
for condition in conditions:
      if condition == “condition_1”:
            model.Reactions.get_by_id(“R1”).bounds=0,1000
            model.Reactions.get_by_id(“R2”).bounds=-1000,1000

      if condition == “condition_2”:
            model.Reactions.get_by_id(“R1”).bounds=0,500
            model.Reactions.get_by_id(“R2”).bounds=-500,500

      with model:
            solution= model.optimize()
      print(“H2 produced when condition = ”, condition, “:”, solution.fluxes[“SK_H2[c]”])

output:

H2 produced when condition = condition_1: 3.08
H2 produced when condition = condition_2: 2.32

However I have found that optimising for a second solution (say, with another objective of interest) in the loop somehow changes the output for the following conditions:

input:

model = cobra.io.sbml.read_sbml_model(“file.xml”)
conditions =[“condition_1”, “condition_2”]
for condition in conditions:
      if condition == “condition_1”:
            model.Reactions.get_by_id(“R1”).bounds=0,1000
            model.Reactions.get_by_id(“R2”).bounds=-1000,1000

      if condition == “condition_2”:
            model.Reactions.get_by_id(“R1”).bounds=0,500
            model.Reactions.get_by_id(“R2”).bounds=-500,500

      with model:
            solution= model.optimize()
      print(“H2 produced when condition = ”, condition, “:”, solution.fluxes[“SK_H2[c]”])

      with model:
            model.objective = new_objective
            new_solution= model.optimize()
      print(“H2 produced when condition = ”, condition, “and objective = new_objective:”, new_solution.fluxes[“SK_H2[c]”])

output:

H2 produced when condition = condition_1: 3.08
H2 produced when condition = condition_1 and objective = new_objective: 0.00

H2 produced when condition = condition_2: 0.00
H2 produced when condition = condition_2 and objective = new_objective: 0.00

I took this result to mean that the model is somehow being modified by the second model.optimize(). I have noticed this same result when performing

with model:
      flux_variability_analysis(model, model.reactions.get_by_id('SK_H2_c[c]'))

after a model.optimize().

Is the only resolution to redefine the model for each condition?:

for condition in conditions:
      if condition == “condition_1”:
            model = cobra.io.sbml.read_sbml_model(“file.xml”)
            model.Reactions.get_by_id(“R1”).bounds=0,1000
            model.Reactions.get_by_id(“R2”).bounds=-1000,1000

      if condition == “condition_2”:
            model = cobra.io.sbml.read_sbml_model(“file.xml”)
            model.Reactions.get_by_id(“R1”).bounds=0,500
            model.Reactions.get_by_id(“R2”).bounds=-500,500

I was reluctant to do this as I have an extensive 'tuning' process in between defining the model and the for loop.

Any help will be appreciated! Thanks in advance.

[cdiener]

You need to move the bounds setting into the with: blocks as well. That guarantees they are automatically reverted when you exit the with block.

[AmalReji]

Thanks for your response. I have tried incorporating the bounds setting and optimization steps within the same with: block, as suggested. Here is the modified code:

def setBoundsFromCondition(model, condition):
      if condition == “condition_1”:
            model.Reactions.get_by_id(“R1”).bounds=0,1000
            model.Reactions.get_by_id(“R2”).bounds=-1000,1000

      if condition == “condition_2”:
            model.Reactions.get_by_id(“R1”).bounds=0,500
            model.Reactions.get_by_id(“R2”).bounds=-500,500

clean_model = model
for condition in conditions:
      # biomass objective
      with clean_model as model:
            setBoundsFromCondition(model, condition)
            solution= model.optimize()
            print(“H2 produced when condition = ”, condition, “:”, solution.fluxes[“SK_H2[c]”])
      # new objective
      with clean_model as model:
            setBoundsFromCondition(model, condition)
            model.objective = new_objective
            new_solution= model.optimize()
            print(“H2 produced when condition = ”, condition, “and objective = new_objective:”, new_solution.fluxes[“SK_H2[c]”])

Summary of what I did to the code:

• Bounds changes and optimizations are now within all the same with:
• Added clean_model to store the model before any possible changes in the for loop
• clean_model is used as the model for every new conditon and every new objective
  - so changes from the previous condition couldn't possibly reflect on subsequent conditions
• Clearly stated the objectives for both reactions

But I still have the problem of H2 flux = 0 after the initial condition.

I suspected one of the model.optimize() and model.objective functions may be responsible, as they call on the model, but the error is now still occuring despite using clean_model for each condition/objective.

I was planning to make a separate post about this, but it may be helpful for the current problem:

• I am getting different fluxes for all reactions (except for biomass) when I change the order of the conditions list at the start
• Whichever condition is first will have H2 flux > 0, but following conditions will have
• Does cobrapy use a stochastic optimisation process for this to occur?

Thanks again!

EDITS: Tried using a clean_model for each condition and added extra information about the problem

[cdiener]

It's not stochastic but you are trying to reproduce non-unique solutions and that will be hard. So basically your objective value is guaranteed to be unique but your fluxes are not since there may be many different flux solutions yielding an optimal value. The solver just picks whatever is easiest for the solver. When you optimize several times the previous solution is used to seed the solver which speeds it up but also can lead to the behavior you are seeing.

The assignment clean_model = model unfortunately does nothing. They both still point to the same object (Python assigns by reference). If you want to reinitialize your model in every iteration you can do something like for model.copy() as mod: ... which will reset the model but will also be pretty slow.

Though, I think the take home message here is that the H2 flux is not relevant to your objective, since it can either be zero or some larger value and still give you optimality. You could also run an FVA on each step for only the H2 flux to get the largest and smallest H2 flux under that condtion.

[AmalReji]

Thank you for your reply and recommendations. Yes, I ended up using

 model = copy.deepcopy(original_model)

which slowed the process down but gives me the results I am looking for. Thanks for explaining the result, too!