model.optimize() seems to disregard model.medium, even when bounds are manually set

[Daniel Farkas]

Dear cobrapy community,

I have ran into a problem with a community model. I used gapseq to generate and gapfill the models and the merge function on kbase to create a compartmentalised community model, after adding exchange reactions for metabolites of interest for individual models.

The FBA solution returns a flux distribution in which many metabolites are taken up, but which were not included in the model.medium dictionary.

It seems even if I manually set lower bounds of exchange reactions of interest to -10, and the rest of the exchanges I constrain to 0, the solution is still the flux distribution with the original model.medium.

I really hope I am making a rookie mistake because it has been an unsolved problem for days.

my code:

model = cobra.io.read_sbml_model("model.xml")

model.medium = {"EX_cpd25328_e0": 10, #triacontane
                "EX_cpd25825_e0": 10, #tetracosane
                "EX_cpd25628_e0":10,#n-decane
                "EX_cpd01034_e0":10,#toluene
                "EX_cpd04136_e0":10,#p-xylene
                "EX_cpd04446_e0":10,#o-xylene
                "EX_cpd01007_e0":10,#benzene
                "EX_cpd00618_e0":10,#naphthalene
                'EX_cpd00013_e0': 100, #NH3
                'EX_cpd00099_e0': 100, #Cl-
                'EX_cpd00254_e0': 100, #Mg
                'EX_cpd00048_e0': 100, #Sulfate
                'EX_cpd10516_e0': 100, #Fe+3
                'EX_cpd00063_e0': 100, #Ca2+
                'EX_cpd00009_e0': 100, #Phosphate
                'EX_cpd00205_e0': 100, #K+
                'EX_cpd00030_e0': 100, #Mn2+
                'EX_cpd00149_e0': 100, #Co2+
                'EX_cpd00058_e0': 100, #Cu2+
                'EX_cpd00034_e0': 100, #Zn2+
                'EX_cpd00067_e0': 0.1, #H+
                'EX_cpd00001_e0': 1000} #H2O

model.objective = model.reactions.get_by_id("bio1_biomass")
solution = model.optimize()
solutions = pd.DataFrame(solution.fluxes)

#check uptake fluxes

solutions.to_csv("fluxes.csv")
uptakes = pd.read_csv("uptakes.csv", sep=",")
uptakes['Metabolite'] = uptakes['Metabolite'].str.strip()
uptakes["met_name"] = "nan"
for i in range(len(uptakes)):
    uptakes["met_name"][i] = model.metabolites.get_by_id(uptakes["Metabolite"][i]).name
print(uptakes)

>>> The output is a flux distribution which does not take into account model.medium.

>>> Manually setting the exchanges:

exchanges = ["cpd25328_e0", #triacontane
                "cpd25825_e0", #tetracosane
                "cpd25628_e0",#n-decane
                "cpd01034_e0",#toluene
                "cpd04136_e0",#p-xylene
                "cpd04446_e0",#o-xylene
                "cpd01007_e0",#benzene
                "cpd00618_e0",#naphthalene
                'cpd00013_e0', #NH3
                'cpd00099_e0', #Cl-
                'cpd00254_e0', #Mg
                'cpd00048_e0', #Sulfate
                'cpd10516_e0', #Fe+3
                'cpd00063_e0', #Ca2+
                'cpd00009_e0', #Phosphate
                'cpd00205_e0', #K+
                'cpd00030_e0', #Mn2+
                'cpd00149_e0', #Co2+
                'cpd00058_e0', #Cu2+
                'cpd00034_e0', #Zn2+
                'cpd00067_e0', #H+
                'cpd00001_e0']

# close all exchanges

for i in range(len(model.exchanges)):
    model.exchanges[i].lower_bound = 0

# reopen correct exchanges
for e in range(len(exchanges)):

    model.reactions.get_by_id("EX_"+str(exchanges[e])).lower_bound = -10

solution = model.optimize()

>>> again, same fluxes uptaking glycerol, etc, that I did not allow.

The step when I add exchange reactions of metabolites of interest (not reproducible code as it draws from modelSEED database etc).

nested_mets_dict has the information of the metabolites, including id, name, formula.

for m  in range(len(models)):
    print(models[m])
    ### Load Model
    model = cobra.io.read_sbml_model(path
                                     +str(models[m]))
    print(len(model.exchanges))
   
    for met in nested_mets_dict:                                                
        print("adding "+nested_mets_dict[met]["id"]+" for c0...")
        model.add_metabolites(Metabolite(str(nested_mets_dict[met]["id"]+"_c0"),
                                            formula=str(nested_mets_dict[met]["formula"]),
                                            name=str(nested_mets_dict[met]["name"]),
                                            compartment="c0"))
       
       
        #this is needed to fix a bug with the assignment of compartment above
        model.metabolites.get_by_id(str(nested_mets_dict[met]["id"]+"_c0")).compartment = "c0"
        print("added")
       
   
        print("adding "+nested_mets_dict[met]["id"]+" for e0...")
        model.add_metabolites(Metabolite(str(nested_mets_dict[met]["id"]+"_e0"),
                                         formula=str(nested_mets_dict[met]["formula"]),
                                         name=str(nested_mets_dict[met]["name"]),
                                         compartment="e0"))
       
        #this is needed to fix a bug with the assignment of compartment above
        model.metabolites.get_by_id(str(nested_mets_dict[met]["id"]+"_e0")).compartment = "e0"
       
       
        print("added")
       
        #create exchange reaction
        print("creating exchange for met "+nested_mets_dict[met]["id"]+"...")
        try:
            model.add_boundary(model.metabolites.get_by_id(str(nested_mets_dict[met]["id"]+"_e0")), type="exchange")
            print("created")
        except ValueError:
            print("ValueError")
            continue

        print("defining transport reaction...")
        #define  Transport  reaction
        reaction = Reaction(str(met+"_transport"))
        reaction.name = str(nested_mets_dict[met]["name"]+"_transport")
        reaction.subsystem = "hydrcarbon_transport"
        reaction.lower_bound = -1000.
        reaction.upper_bound = 1000.

        print("done")
       
        # This might be breaking the model in regards to exchanges - duplicated of adding mets above?
        #define intracellular  metabolite
        #print("defining intracellular metabolite...")
        #globals()[str(nested_mets_dict[met]["id"]+"_c")] = Metabolite(str(nested_mets_dict[met]["id"]+"_c0"),
        #                                                                        formula=str(nested_mets_dict[met]["formula"]),
        #                                                                       name=str(nested_mets_dict[met]["name"]),
        #                                                                        compartment="c0")

        #print("done")
        #print("defining extracellular metabolite...")
        #define extracellular  metabolite
        #globals()[str(nested_mets_dict[met]["id"]+"_e")] = Metabolite(str(nested_mets_dict[met]["id"]+"_e0"),
        #                                                                        formula=str(nested_mets_dict[met]["formula"]),
        #                                                                        name=str(nested_mets_dict[met]["name"]),
        #                                                                       compartment="e0")
        #print("done")
       
        print("adding metabolites to rxn...")
        reaction.add_metabolites({model.metabolites.get_by_id(str(nested_mets_dict[met]["id"]+"_c0")):1.0,
                                  model.metabolites.get_by_id(str(nested_mets_dict[met]["id"]+"_c0")):-1.0})

        print("adding reaction to model...")
        model.add_reactions([reaction])
        print("done")
   
   
    print("writing model...............")
    cobra.io.write_sbml_model(model, str(models[m].split(".fasta")[0]+"_exchanges_added.xml"))

If anyone can spot anything, or heard of a similar problem, please let me know!!!!!

Thank you very much in advance!

Daniel

[ch.d...@gmail.com]

Hard to debug without having access to the model. I would first check that model.exchanges looks okay and do something like model.exchanges[0].reaction to see if the default direction is what you would expect. Also, how large are your impossible uptakes? Deviations within the solver tolerances (<~1e-6) are to be expected.