Model incapable of using metabolite generated by other model

[rcolpo]

Hi all!

I'm a beginner with MICOM, and I am trying to understand how it works. To better grasp its functionalities, I built two toy "metabolic models" using COBRApy: toy1 and toy2. In COBRApy, toy1 can grow when executing toy1.optimize(). However, toy2 cannot grow because it's unable to generate the metabolites "C_c" or "C_e". Conversely, toy1 can produce "C_e". On MICOM, I anticipated that toy1 would generate "C_e" and toy2 would intake it, enabling both to grow.

For toy1 to produce "C_e", this must be a constraint in the community model, as generating "C_e" would reduce toy1's biomass. My question is: How can I enforce a solution where toy2 grows at a specific ratio relative to toy1? For instance, can toy2 have 50% of the toy1 biomass?

Below is the code I used on my example.

First, I created the toy models using cobrapy:

import cobra


#crete model toy1. is able to grow. It can also generate C_e.

toy1 = cobra.Model('toy1')

my_reaction = cobra.Reaction("EX_A")
toy1.add_reactions([my_reaction])
my_reaction.reaction = " -> A_e"
        
my_reaction = cobra.Reaction("A_transport")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_e <-> A_c"

my_reaction = cobra.Reaction("C_transport")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "C_c <-> C_e"

my_reaction = cobra.Reaction("A_C")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> C_c"

my_reaction = cobra.Reaction("biomasToy1")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> D_c"

my_reaction = cobra.Reaction("DM_A")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> "

my_reaction = cobra.Reaction("DM_Cc")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "C_c -> "

my_reaction = cobra.Reaction("DM_Ce")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "C_e -> "

my_reaction = cobra.Reaction("DM_D")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "D_c -> "

toy1.objective = "biomasToy1"

for met in toy1.metabolites:
    if met.id.endswith('_e'): met.compartment = 'e'
    else: met.compartment = 'c'
    
cobra.io.write_sbml_model(toy1, "C:/UFZ/Jens/Micom/models/toy1.xml")


#crete model toy2. Should not be able to grow, because thre is no input of C_e

toy2 = cobra.Model('toy2')

my_reaction = cobra.Reaction("EX_A")
toy2.add_reactions([my_reaction])
my_reaction.reaction = " -> A_e"

my_reaction = cobra.Reaction("EX_C")
toy2.add_reactions([my_reaction])
my_reaction.reaction = " -> C_e"
        
my_reaction = cobra.Reaction("A_transport")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "A_e <-> A_c"

my_reaction = cobra.Reaction("C_transport")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "C_e <-> C_c"

my_reaction = cobra.Reaction("biomasToy2")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "A_c + C_c -> B_c"

my_reaction = cobra.Reaction("DM_B")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "B_c -> "

my_reaction = cobra.Reaction("DM_A")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> "

toy2.objective = "biomasToy2"

for met in toy2.metabolites:
    if met.id.endswith('_e'): met.compartment = 'e'
    else: met.compartment = 'c'
    
cobra.io.write_sbml_model(toy2, "C:/UFZ/Jens/Micom/models/toy2.xml")

Then, I used MICOM to generate the community model.

## generating community model using MICOM
import pandas as pd
import micom

pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)

tax = pd.DataFrame({"id": ["toy1", "toy2"], "file": ["C:/UFZ/Jens/Micom/models/toy1.xml", "C:/UFZ/Jens/Micom/models/toy2.xml"], 'abundance':[1,0.5]})
com = micom.Community(tax)

sol = com.optimize(fluxes=True,pfba=True)

Below are the fluxes I got:

>>> print(sol.fluxes)
    
reaction     A_C  A_transport  C_transport  DM_A  DM_B  DM_Cc  DM_Ce       DM_D       EX_A     EX_A_m  biomasToy1  biomasToy2
compartment                                                                                                                  
medium       NaN          NaN          NaN   NaN   NaN    NaN    NaN        NaN        NaN -66.666599         NaN         NaN
toy1         0.0    99.999899          0.0   0.0   NaN    0.0    0.0  99.999899  99.999899        NaN   99.999899         NaN
toy2         NaN     0.000000          0.0  -0.0  -0.0    NaN    NaN        NaN   0.000000        NaN         NaN         0.0

Thanks,
Rodrigo.

[cdiener]

Hi, when I run your code both organisms grow.

In [14]: sol.members
Out[14]: 
              abundance  growth_rate  reactions  metabolites
compartments                                                
medium              NaN          NaN          2            2
toy1           0.666667    99.999899          9            5
toy2           0.333333    99.999899          7            5

Maybe you set another medium? By default, the growth solutions won't be unique though. Using cooperative tradeoff will push the solution towards growth for both taxa (see here for some more background). So running sol = com.cooperative_tradeoff(fraction=1.0, fluxes=True) should get you what you want. If you specifically want to constrain that both growth rates fulfill a certain ratio, that is possible too by adding another constraint. Here is an example where BM(toy) = 2*BM(toy2):

bm1 = com.reactions.biomasToy1__toy1
bm2 = com.reactions.biomasToy2__toy2
const = com.problem.Constraint(
    bm1.flux_expression - 2 * bm2.flux_expression, 
    name="biomass_constraint",
    lb=0, ub=0
)
with com:
    com.add_cons_vars([const])
    const_sol = com.optimize()
const_sol.members

Which will give:

              abundance  growth_rate  reactions  metabolites
compartments                                                
medium              NaN          NaN          2            2
toy1           0.666667        100.0          9            5
toy2           0.333333         50.0          7            5

[rcolpo]

Dear Christian,

Thank you for taking the time to answer my question.
The reason why I got a different result is that I forgot to comment out the my_reaction = cobra.Reaction("EX_C") reaction from toy2 in my original question. This changes the media I used, as you suggested.

With "EX_C" removed from toy2, I included a new constraint to force flux on "biomasToy2", but got the error "WARNING solver encountered an error infeasible". Below is how I included the constraint:

bm2 = com.reactions.biomasToy2__toy2
const = com.problem.Constraint(
	bm2.flux_expression,
	name="bm2_constraint",
	lb=1, ub=1000
)

with com:
	com.add_cons_vars([const])
	const_sol = com.optimize()

toy2 can not grow because it does not have "C_e" available (because now it was removed from the model)... but toy1 should be able to provide it. This behavior is just like it does not recognize the "C_e" coming from toy1, as the same "C_e" that toy2 is trying to uptake. Maybe metabolites can not be exchanged if they are not initially present in the media?

P.S.:
To double-check whether toy1 could really produce "C_e", I forced flux on com.reactions.DM_Ce__toy1, and it worked. But it failed when I tried to force the transport of "C_e" inside "toy2" (com.reactions.C_transport__toy2).

Below you find how I tried to force "C_e" out of toy1 and into toy2.

C_transport_toy1 = com.reactions.C_transport__toy1
const_C_transport_toy1 = com.problem.Constraint(
	C_transport_toy1.flux_expression,
	name="DM_Ce_constraint",
	lb=10, ub=1000
)

C_transport_toy2 = com.reactions.C_transport__toy2
const_C_transport_toy2 = com.problem.Constraint(
	C_transport_toy2.flux_expression,
	name="C_transport_constraint",
	lb=10, ub=1000
)

with com:
	com.add_cons_vars([const_C_transport_toy1])
	com.add_cons_vars([const_C_transport_toy2])
	const_sol = com.optimize(fluxes=True)

[cdiener]

Okay, understood. So if you remove EX_C from model toy2 that model will not be able to import/take up C so that's why it will never import it. Similarly toy needs an exchange reaction so C can be secreted into the environment. The way to force this interaction is to leave in the import reaction but set the exchange flux in the community medium to zero (not present in the environment). That way toy2 can take up C but it has to be produced by other taxa.

So your models would be built like that (note the <-> arrows to allow export of C in model 1):

import cobra


#crete model toy1. is able to grow. It can also generate C_e.

toy1 = cobra.Model('toy1')

my_reaction = cobra.Reaction("EX_A")
toy1.add_reactions([my_reaction])
my_reaction.reaction = " <-> A_e"
 
my_reaction = cobra.Reaction("EX_C")
toy1.add_reactions([my_reaction])
my_reaction.reaction = " <-> C_e"
       
my_reaction = cobra.Reaction("A_transport")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_e <-> A_c"

my_reaction = cobra.Reaction("C_transport")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "C_c <-> C_e"

my_reaction = cobra.Reaction("A_C")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> C_c"

my_reaction = cobra.Reaction("biomasToy1")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> D_c"

my_reaction = cobra.Reaction("DM_A")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> "

my_reaction = cobra.Reaction("DM_Cc")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "C_c -> "

my_reaction = cobra.Reaction("DM_Ce")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "C_e -> "

my_reaction = cobra.Reaction("DM_D")
toy1.add_reactions([my_reaction])
my_reaction.reaction = "D_c -> "

toy1.objective = "biomasToy1"

for met in toy1.metabolites:
    if met.id.endswith('_e'): met.compartment = 'e'
    else: met.compartment = 'c'
    
cobra.io.write_sbml_model(toy1, "C:/UFZ/Jens/Micom/models/toy1.xml")


#crete model toy2. Should not be able to grow, because thre is no input of C_e

toy2 = cobra.Model('toy2')

my_reaction = cobra.Reaction("EX_A")
toy2.add_reactions([my_reaction])
my_reaction.reaction = " <-> A_e"

my_reaction = cobra.Reaction("EX_C")
toy2.add_reactions([my_reaction])
my_reaction.reaction = " <-> C_e"
        
my_reaction = cobra.Reaction("A_transport")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "A_e <-> A_c"

my_reaction = cobra.Reaction("C_transport")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "C_e <-> C_c"

my_reaction = cobra.Reaction("biomasToy2")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "A_c + C_c -> B_c"

my_reaction = cobra.Reaction("DM_B")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "B_c -> "

my_reaction = cobra.Reaction("DM_A")
toy2.add_reactions([my_reaction])
my_reaction.reaction = "A_c -> "

toy2.objective = "biomasToy2"

for met in toy2.metabolites:
    if met.id.endswith('_e'): met.compartment = 'e'
    else: met.compartment = 'c'
    
cobra.io.write_sbml_model(toy2, "C:/UFZ/Jens/Micom/models/toy2.xml")

That way you get the trophic interaction even if there is no C in the environment:

import pandas as pd
import micom
 
pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)
 
tax = pd.DataFrame({"id": ["toy1", "toy2"], "file": ["toy1.xml", "toy2.xml"], 'abundance':[1,0.5]})
com = micom.Community(tax)

com.reactions.EX_C_m.bounds = 0, 0  # <- remove from env
sol = com.cooperative_tradeoff(fraction=1.0, fluxes=True)

In [15]: sol.members
Out[15]: 
              abundance  growth_rate  reactions  metabolites
compartments                                                
medium              NaN          NaN          2            2
toy1           0.666667         80.0         10            5
toy2           0.333333         40.0          7            5

In [16]: sol.fluxes[["EX_C", "EX_C_m"]]
Out[16]: 
reaction     EX_C  EX_C_m
compartment              
medium        NaN     0.0
toy1        -20.0     NaN
toy2         40.0     NaN

[rcolpo]

Dear Christian,

Thank you for your detailed answer.
The sol.fluxes.EX_C['toy2'] is equal to 40.0. This means that "C_e" is available in the media to toy2, right? But I don't see signs that the C_e is actually coming from toy1.

What I want to test is if toy2 can grow with a metabolite produced and excreted by toy1. Therefore, toy1 can produce C_c in reaction "A_C", and excrete it in reaction "C_transport". Regarding the production and transport of "C_c", the reactions I would expect to see active in the models are:

In "toy1":
" <-> A_e"
"A_e <-> A_c"
"A_c -> C_c"
"C_c <-> C_e"

In "toy2":

"C_e <-> C_c"
"A_c + C_c -> B_c"

The " -> C_e" reaction should not be active in toy1 or toy2.

As you mentioned, we can simulate the absence of "C_e" in the toy1's and toy2's media by changing the boundaries of the "EX_C" reaction. But, would it be better to make instead com.reactions.EX_C__toy1.bounds = -1000, 0 and com.reactions.EX_C__toy2.bounds = -1000, 0? However, when I use these commands, toy2 can not grow.

[cdiener]

Individual taxa are compartmentalized in MICOM. So for C to be used in a trophic interaction it has to be produced by toy1, be secreted into the medium (compartment "m") via EX_C__toy1, be imported into toy2 by EX_C__toy2 and then be used to produce biomass in toy2. All of this is happening in my example above. Note that the compartment "e" is not the shared environment for both taxa, that is the compartment "m". "e" is the close extracellular space for each individual taxon that is not shared with other organisms. So if you have a transport <-> C_e MICOM will convert this to C_m <-> C_e and add a global <-> C_m as well.

Finally, in the fluxes you see that toy1 is producing C and toy2 is consuming it. Those fluxes are normalized to 1gDW of toy1 and toy2 but because their abundances are in a 2:1 ratio the net fluxes is actually balanced and all the secretion flux of C from toy1 is consumed by toy2.