Adding metabolites to model clears model.exchanges

[matthiaskoenig]

This is a strange bug when adding exchange reactions to an existing model. See discussion here: https://groups.google.com/d/msgid/cobra-pie/e5ba650b-d377-4b4a-8e6c-f134db83a214n%40googlegroups.com?utm_medium=email&utm_source=footer I could confirm the strange behavior with the latest develop version. Not sure what is going on here, but definitely the expected behavior.

import cobra

# load model
model = cobra.io.read_sbml_model('model copy.xml')

# check amount of reactions and metabolites
print('before:')
print('reactions:', len(model.reactions))
print('metabolites:', len(model.metabolites))
print('genes:', len(model.genes))
print('exchange reactions:', len(model.exchanges))
print('medium = ', model.medium)

# define metabolites (THIS REMOVES THE model.exchanges !!!)
M_3php_e = cobra.Metabolite(id="3php_e", compartment="e", name="3phosphohydroxypyruvate")
M_3php_p = cobra.Metabolite(id="3php_p", compartment="p", name="3phosphohydroxypyruvate")
model.add_metabolites([M_3php_e, M_3php_p])
print('exchange reactions after metabolite:', len(model.exchanges))

# define and add transport reaction 
R_3PHPtex = cobra.Reaction(id="3PHPtex", name="3phosphohydroxypyruvate transport via diffusion extracellular to periplasm", 
                           lower_bound=-1000, upper_bound=1000)
R_3PHPtex.add_metabolites({M_3php_e: -1.0, M_3php_p: 1.0})

model.add_reactions([R_3PHPtex]) # metabolites are added automatically
print('exchange reactions after reaction:', len(model.exchanges))

# define and add exchange reaction
model.add_boundary(model.metabolites.get_by_id('3php_e'), type='exchange')

# now to check if the reactions have the names you expect
print('\n reactions added:\n', model.reactions.get_by_id('3PHPtex')) # '\n' is a line break in python
print(model.reactions.get_by_id('EX_3php_e'))
print('exchange reactions:', len(model.exchanges))

# add 3php_e to the medium
[model copy.zip](https://github.com/opencobra/cobrapy/files/6024239/model.copy.zip)

medium = model.medium
print(medium)
medium['EX_3php_e'] = 10
model.medium = medium

# print info again
print('\n after:')
print('reactions :', len(model.reactions))
print('metabolites :', len(model.metabolites))
print('genes :', len(model.genes))
print('exchange reactions:', len(model.exchanges))
print('medium = ', model.medium)

Results in the following output (note exchanges becoming 0):

before:
reactions: 2894
metabolites: 1802
genes: 1861
exchange reactions: 338
medium =  {'EX_glc__D_e': 10.0, 'EX_h2o_e': 10.0, 'EX_h_e': 10.0, 'EX_cl_e': 10.0, 'EX_pi_e': 10.0, 'EX_nh4_e': 10.0, 'EX_fe3_e': 10.0, 'EX_k_e': 10.0, 'EX_ca2_e': 10.0, 'EX_mg2_e': 10.0, 'EX_mn2_e': 10.0, 'EX_cobalt2_e': 10.0, 'EX_zn2_e': 10.0, 'EX_cu2_e': 10.0, 'EX_o2_e': 10.0, 'EX_fe2_e': 10.0, 'EX_mobd_e': 10.0, 'EX_so4_e': 10.0}
exchange reactions after metabolite: 0
exchange reactions after reaction: 0

 reactions added:
 3PHPtex: 3php_e <=> 3php_p
EX_3php_e: 3php_e <=> 
exchange reactions: 1
{'EX_3php_e': 1000.0}

 after:
reactions : 2896
metabolites : 1804
genes : 1861
exchange reactions: 1
medium =  {'EX_3php_e': 10}

Model attached.

[cdiener]

I think your attachment did not work.

Looks like it does not identify the external compartment anymore. Maybe the initial model was lacking the "e" compartment? What is the output of running:

from cobra.medium import find_external_compartment

# initial loading
model.comaprtments
find_external_compartment(model)

# adding the reactions
model.compartments
find_external_compartment(model)

[cdiener]

Found the model. The problem are incorrect names for the compartment:

In [4]: model.compartments
Out[4]: {'C_c': 'cytosol', 'C_p': 'periplasm', 'C_e': 'extracellular space'}

So this would need to be corrected or the new metabolites have to be added with the correct names:

# check amount of reactions and metabolites
print('before:')
print('reactions:', len(model.reactions))
print('metabolites:', len(model.metabolites))
print('genes:', len(model.genes))
print('exchange reactions:', len(model.exchanges))
print('medium = ', model.medium)

# define metabolites (THIS REMOVES THE model.exchanges !!!)
M_3php_e = cobra.Metabolite(id="3php_e", compartment="C_e", name="3phosphohydroxypyruvate")
M_3php_p = cobra.Metabolite(id="3php_p", compartment="C_p", name="3phosphohydroxypyruvate")
model.add_metabolites([M_3php_e, M_3php_p])
print('exchange reactions after metabolite:', len(model.exchanges))

# define and add transport reaction 
R_3PHPtex = cobra.Reaction(id="3PHPtex", name="3phosphohydroxypyruvate transport via diffusion extracellular to periplasm", 
                           lower_bound=-1000, upper_bound=1000)
R_3PHPtex.add_metabolites({M_3php_e: -1.0, M_3php_p: 1.0})

model.add_reactions([R_3PHPtex]) # metabolites are added automatically
print('exchange reactions after reaction:', len(model.exchanges))

# define and add exchange reaction
model.add_boundary(model.metabolites.get_by_id('3php_e'), type='exchange')

# now to check if the reactions have the names you expect
print('\n reactions added:\n', model.reactions.get_by_id('3PHPtex')) # '\n' is a line break in python
print(model.reactions.get_by_id('EX_3php_e'))
print('exchange reactions:', len(model.exchanges))

medium = model.medium
print(medium)
medium['EX_3php_e'] = 10
model.medium = medium

# print info again
print('\n after:')
print('reactions :', len(model.reactions))
print('metabolites :', len(model.metabolites))
print('genes :', len(model.genes))
print('exchange reactions:', len(model.exchanges))
print('medium = ', model.medium)

which gives:

before:
reactions: 2894
metabolites: 1802
genes: 1861
exchange reactions: 338
medium =  {'EX_glc__D_e': 10.0, 'EX_h2o_e': 10.0, 'EX_h_e': 10.0, 'EX_cl_e': 10.0, 'EX_pi_e': 10.0, 'EX_nh4_e': 10.0, 'EX_fe3_e': 10.0, 'EX_k_e': 10.0, 'EX_ca2_e': 10.0, 'EX_mg2_e': 10.0, 'EX_mn2_e': 10.0, 'EX_cobalt2_e': 10.0, 'EX_zn2_e': 10.0, 'EX_cu2_e': 10.0, 'EX_o2_e': 10.0, 'EX_fe2_e': 10.0, 'EX_mobd_e': 10.0, 'EX_so4_e': 10.0}
exchange reactions after metabolite: 338
exchange reactions after reaction: 338

 reactions added:
 3PHPtex: 3php_e <=> 3php_p
EX_3php_e: 3php_e <=> 
exchange reactions: 339
{'EX_glc__D_e': 10.0, 'EX_h2o_e': 10.0, 'EX_h_e': 10.0, 'EX_cl_e': 10.0, 'EX_pi_e': 10.0, 'EX_nh4_e': 10.0, 'EX_fe3_e': 10.0, 'EX_k_e': 10.0, 'EX_ca2_e': 10.0, 'EX_mg2_e': 10.0, 'EX_mn2_e': 10.0, 'EX_cobalt2_e': 10.0, 'EX_zn2_e': 10.0, 'EX_cu2_e': 10.0, 'EX_o2_e': 10.0, 'EX_fe2_e': 10.0, 'EX_mobd_e': 10.0, 'EX_so4_e': 10.0, 'EX_3php_e': 1000.0}

 after:
reactions : 2896
metabolites : 1804
genes : 1861
exchange reactions: 339
medium =  {'EX_glc__D_e': 10.0, 'EX_h2o_e': 10.0, 'EX_h_e': 10.0, 'EX_cl_e': 10.0, 'EX_pi_e': 10.0, 'EX_nh4_e': 10.0, 'EX_fe3_e': 10.0, 'EX_k_e': 10.0, 'EX_ca2_e': 10.0, 'EX_mg2_e': 10.0, 'EX_mn2_e': 10.0, 'EX_cobalt2_e': 10.0, 'EX_zn2_e': 10.0, 'EX_cu2_e': 10.0, 'EX_o2_e': 10.0, 'EX_fe2_e': 10.0, 'EX_mobd_e': 10.0, 'EX_so4_e': 10.0, 'EX_3php_e': 10}

If reactions are not annotated with SBO terms, cobrapy has a list of common names for the external compartment it looks through. "e" is in that list but "C_e" is not, so it will use the new compartment here.

One can work around that though and specify the external compartment, see https://github.com/opencobra/cobrapy/issues/747 for in depth discussion on that.

[matthiaskoenig]

@cdiener Thanks for the quick solution.