Building metabolic models with signaling compound production or response
Once users are familiar with the add_signal and add_multitoxin functions, the first step to running simulations with signaling compounds is building metabolic models that can accommodate the production of or response to these compounds. We will detail several model-building processes, including the approach taken in our paper.
Adding new signals to existing models
Many users will be most interested in adding the signaling capability to existing genome-scale metabolic models. We will walk through this process by using an existing model as an example.
E_cobra = cobra.io.load_model('textbook') # basic IJO1366 model
As with any COMETS simulation, cobra models need to be converted to COMETS models prior to usage. This is a straightforward process in the case of existing models.
import cometspy as c
E = c.model(E_cobra)
E.open_exchanges() # ensures that models can take up any nutrients provided in the COMETS media
E.obj_style = "MAX_OBJECTIVE_MIN_TOTAL" # type of flux solution
Setting up a COMETS simulation also involves specifying the initial populations, creating a layout, and seeding the media with the relevant components. More information about this process is available in the COMETS documentation, as many of the specific details will vary significantly based on the desires of the user.
E.initial_pop = [0, 0, 1.e-7] # the first two values are x, y coordinates, third is gDW for starting biomass
layout = c.layout(E)
E_media = {'h_e': 1000.0,
'h2o_e': 1000.0,
'nh4_e': 1000.0,
'o2_e': 1000.0,
'pi_e': 1000.0}
for key, value in E_media.items():
layout.set_specific_metabolite(key, value)
layout.set_specific_metabolite("glc__D_e", 5e-4)
Once users have decided on many of the core parameters for their simulations, they may want to add a signaling compound. In this example, we will consider a nondegradable, bactericidal toxin whose purpose is to kill E. coli in a linear way. This toxin will exist in the media from the outset of simulations.
To make a toxin, we will assign a signal to the E. coli model that alters a reaction in that model based upon the concentration of an environmental metabolite. We will make a new metabolite for this purpose, although this isn't strictly necessary, as we will see later.
When setting a new metabolite as a toxin or signaling compound, the affected model will need to be altered to have an exchange reaction for that metabolite. We will use functions from cobra to do that.
from cobra import Metabolite, Reaction
toxin_e = Metabolite(id = "toxin_e", compartment = "e") # create a toxin metabolite in the extracellular environment
EX_toxin_e = Reaction(id = "EX_toxin_e",
lower_bound = -1000.,
upper_bound = 1000.) # create an exchange reaction for the toxin
EX_toxin_e.add_metabolites({toxin_e: -1}) # add toxin metabolite and stoichiometric coefficient, which in this case is -1
E_cobra.add_reactions([EX_toxin_e]) # add the exchange reaction to E coli
Now, the E. coli metabolic model can sense the toxic metabolite in the environment, which is essential for its sensitivity to the compound. However, we have yet to define how the compound will impact the growth rate or survival of E. coli, which we will do with the add_signal function. Here, we are interested in defining a toxin that will kill E. coli ('death') without changing the external concentration of the metabolite. These changes need to be made to the COMETS version of the model.
response = E.reactions.loc[E.reactions.REACTION_NAMES == "Biomass_Ecoli_core",:].ID.item() # reaction altered by the toxin
signal_exch = E.reactions.loc[E.reactions.REACTION_NAMES == "EX_toxin_e",:].EXCH_IND.item() # exch_id of the toxin
E.add_signal('death', signal_exch, 'met_unchanged', 'linear', parms = [0, 0.1])
The death rate of E. coli will now be increased based on the extracellular concentration of the toxin based on the linear parameter values 0 and 0.1.
We are now prepared to run a simulation with an E. coli model that will be sensitive to a toxin. In this case, we will need to add that toxin into the extracellular environment as an additional media component prior to running our simulation.
toxin_mmol = 0.6 # set the concentration of toxin in the environment
layout.set_specific_metabolite('toxin_e', toxin_mmol) # add toxin into the environment
Turning existing metabolites into signaling compounds
In some cases, instead of creating an entirely new metabolite that serves as a signaling compound, users will want models to respond to existing metabolites. We detail that process here, based on the textbook E. coli model.
E_cobra = cobra.io.load_model('textbook')
We can consider a situation where the extracellular build-up of the existing metabolite acetaldehyde (EX_acald_e) reduces the function of phosphofructokinase (PFK). The basics of the model set up will not change much from our previous example.
# create COMETS model
E = c.model(E_cobra)
E.open_exchanges()
E.initial_pop = [0, 0, 0.01]
# add signaling relationship without adding a new metabolite
PFK_num = E.reactions.loc[E.reactions.REACTION_NAMES == "PFK", "ID"].values[0]
acald_exch_id = E.reactions.loc[E.reactions.REACTION_NAMES == "EX_acald_e", "EXCH_IND"].values[0]
In this case, we will need to be sensitive to the maximum flux through the PFK reaction in order to appropriately reduce it through the build up of acetaldehyde. To find the maximum flux, we can solve using the cobra model:
sol = E_cobra.optimize()
max_ub = sol.fluxes["PFK"]
print(max_ub)
intercept = max_ub
This value will be useful for creating our final signal. Because we are reducing the function of a metabolic reaction rather than growth rate, we will also change the linear relationship; the flux through PFK can be negative, so we will use the linear function_relationship rather than bounded_linear:
slope = -1000 # this is a fairly arbitrary relationship for the slope between acetaldehyde concentration and flux through PFK
E.add_signal(PFK_num, acald_exch_id, "ub", "linear", [slope, max_ub])
Here, we have defined a signal such that the maximum flux through PFK max_ub will reduce linearly with a slope of -1000 as the concentration of acetaldehyde increases. We can run our simulations in a couple of ways. First, as before, we will observe the consequences of a constant concentration of extracellular signaling compound.
layout = c.layout(E)
layout.set_specific_metabolite("glc__D_e", 5.)
layout.set_specific_metabolite("nh4_e", 1000.)
layout.set_specific_metabolite("pi_e", 1000.)
layout.set_specific_metabolite("acald_e", 0.25)
Using this approach, our signaling compound is present in the media at an initial concentration that will be reduced over simulation time as it is taken up by E. coli.
We may also be interested in a signaling compound that is added dynamically, whether to replicate pulses by a producer or another scenario. The simplest way to do this is by using several available tools in COMETS, one of which is the ability to "refresh" compounds, adding a certain mmol of compound per spatial box per hour.
layout.set_specific_metabolite("acald_e", 0.) # reset
layout. set_specific_refresh("acald_e", 0.25)
Making toy producer models using cobra
Finally, users may want to make toy models, such as the ones used in our publication. Functionally, these models will not significantly differ from existing models. The basic premise of their construction is identical to the previous examples. Their benefit is that they are more easily manipulated, since every exchange and reaction is chosen and specifically defined by the user.
The most straightforward way to generate toy metabolic models is to create a function that defines the desired models using a variety of functions and object types from cobrapy. So far, we have only seen how to add signaling compounds that exist in the extracellular environment and impact a strain growing in isolation. This function provides a way to create a metabolic model that will produce the signaling compound. This same basic process can be applied to existing models, if users would like to add toxin production to existing models.
def make_RPS_cobra_models_biomass_cost(growth_rate = 1., toxin_cost = 0.02, resistance_cost = 0.0, toxin_prod = 1.):
# create metabolites for carbon and toxin
carbon_e = Metabolite(id = "carbon_e",
compartment = "e")
carbon_c = Metabolite(id = "carbon_c",
compartment = "c")
toxin_c = Metabolite(id = "toxin_c", compartment = "c")
toxin_e = Metabolite(id = "toxin_e", compartment = "e")
# create the reactions - biomass, carbon exchange and transport, and toxin exchange and transport
#exchange reactions
EX_carbon_e = Reaction(id = "EX_carbon_e",
lower_bound = -growth_rate, # growth rate
upper_bound = 1000.)
EX_carbon_e.add_metabolites({carbon_e: -1})
EX_toxin_e = Reaction(id = "EX_toxin_e",
lower_bound = -1000.,
upper_bound = 1000.)
EX_toxin_e.add_metabolites({toxin_e: -1})
# transport reactions
carbon_transfer = Reaction(id = "carbon_transfer",
lower_bound = 0.,
upper_bound = 1000.)
carbon_transfer.add_metabolites({carbon_e: -1,
carbon_c: 1})
toxin_transfer = Reaction(id = "toxin_transfer",
lower_bound = -1000.,
upper_bound = 0.)
toxin_transfer.add_metabolites({toxin_e: -1,
toxin_c: 1})
# biomass reactions
Biomass = Reaction(id = "Biomass",
lower_bound = 0.,
upper_bound = 1000.)
Biomass.add_metabolites({carbon_c: -1.})
Biomass_producer = Reaction(id = "Biomass",
lower_bound = 0.,
upper_bound = 1000.)
Biomass_producer.add_metabolites({carbon_c: -(1. + toxin_cost + resistance_cost),
toxin_c: toxin_prod * (1. + toxin_cost + resistance_cost)})
Biomass_resistant = Reaction(id = "Biomass",
lower_bound = 0.,
upper_bound = 1000.)
Biomass_resistant.add_metabolites({carbon_c: -(1. + resistance_cost)})
# generate the 3 model types with the appropriate reactions and objectives
producer = Model("producer")
producer.add_reactions([EX_carbon_e, carbon_transfer, Biomass_producer,
EX_toxin_e, toxin_transfer])
producer.objective = Biomass_producer
resistant = Model("resistant")
resistant.add_reactions([EX_carbon_e, carbon_transfer, Biomass_resistant])
resistant.objective = Biomass_resistant
susceptible = Model("susceptible")
susceptible.add_reactions([EX_carbon_e, carbon_transfer, Biomass,
EX_toxin_e])
susceptible.objective = Biomass
return((producer, resistant, susceptible))
This function takes four arguments (growth rate, the cost of toxin production, the cost of toxin resistance, and the amount of toxin produced per gram of biomass produced) to create three model genotypes (susceptibles, toxin producers, and resistant cheaters). Susceptible models have four metabolic reactions (carbon exchange, carbon transport, biomass production, and extracellular to intracellular toxin exchange), producers have five metabolic reactions (carbon exchange, carbon transport, biomass production, intracellular to extracellular toxin exchange, and toxin transport), and resistant cheaters have three metabolic reactions (carbon exchange, carbon transfer, and biomass production).
When manipulating this function to generate toy models, particular attention should be paid to the construction of the biomass reactions, which will allow users to change the relative growth rates of each model (the default here is that they all have the same growth rate, aside from minimal penalities for the cost of toxin production or resistance) along with a variety of other assumptions. Here, we have tied toxin production to producer biomass and carbon uptake, which differs from previous examples, where signaling compounds already existed in the media:
Biomass_producer.add_metabolites({carbon_c: -(1. + toxin_cost + resistance_cost),
toxin_c: toxin_prod * (1. + toxin_cost + resistance_cost)})
This setup defines toxin production in such a way that toxins are only created during growth and that the costs associated with toxin production (toxin_cost) or resistance to toxin (resistance_cost) vary directly with growth rate. It also establishes that toxin production is directly proportional to carbon uptake, regardless of growth rate, although the amount of toxin created per unit of carbon can be changed (toxin_prod). If desired, users can alter these constraints to fit the appropriate biological details of the signaling compounds being modeled.
Notably, the function defined above also penalizes the biomass production of resistant cheaters to impose a cost of toxin resistance (resistance_cost), seen here:
Biomass_resistant.add_metabolites({carbon_c: -(1. + resistance_cost)})
Many of these modeling assumptions can be adjusted by changing the biomass reactions for any of the three models, or by altering the four arguments in the function call:
growth_rate = 0.75 # growth rate for all three genotypes
toxin_cost = 0.02 # 2% reduction in growth rate as a result of producing energetically expensive toxin
resistance_cost = 0 # no reduction in growth rate as a result of being resistant to toxin
toxin_prod = 15 # 15 mmol of toxin produced per gram of producer biomass
producer, resistant, susceptible = make_RPS_cobra_models_biomass_cost(growth_rate = growth_rate, toxin_cost = toxin_cost, resistance_cost = resistance_ost, toxin_prod = toxin_prod)
Making COMETS models from cobra toy models
As we have already shown, once cobrapy toy models have been constructed, they need to be converted into models that are recognizable by cometspy. Exactly how this is done will vary slightly for each type of model and its associated biological details, and is somewhat more complicated when using toy models than pre-existing ones.
P = c.model(producer) # create the comets model from the cobra producer
R = c.model(resistant)
S = c.model(susceptible)
P.obj_style = "MAX_OBJECTIVE_MIN_TOTAL" # set flux style as FBA (pFBA is also possible)
R.obj_style = "MAX_OBJECTIVE_MIN_TOTAL"
S.obj_style = "MAX_OBJECTIVE_MIN_TOTAL"
Diverging from curated models, when using a toy model, during the initial conversion, COMETS will erroneously assume that the biomass equation is an exchange for the susceptible and resistant models, because the reactions are unbalanced (unlike the producer, for which carbon uptake is balanced by toxin production). That will need to be fixed first.
# tell COMETS that the biomass reaction is not an exchange
R.reactions.loc[R.reactions.REACTION_NAMES == "Biomass", "EXCH"] = False # set false
index = R.reactions.loc[R.reactions.REACTION_NAMES == "Biomass", "EXCH_IND"].values[0]
R.reactions.loc[R.reactions.REACTION_NAMES == "Biomass", "EXCH_IND"] = 0
R.reactions.loc[R.reactions.EXCH_IND > index, "EXCH_IND"] -= 1
S.reactions.loc[S.reactions.REACTION_NAMES == "Biomass", "EXCH"] = False
index = S.reactions.loc[S.reactions.REACTION_NAMES == "Biomass", "EXCH_IND"].values[0]
S.reactions.loc[S.reactions.REACTION_NAMES == "Biomass", "EXCH_IND"] = 0
S.reactions.loc[S.reactions.EXCH_IND > index, "EXCH_IND"] -= 1
Next, we will need to add the signaling capability to the susceptible model, so that the uptake of toxin results in a proportionate reduction in growth rate.
gr_absense_of_toxin = growth_rate # baseline growth rate
toxin_conc_where_effect_starts = 0. # threshold concentration where effect starts
slope_of_toxin_effect = -growth_rate # slope of the relationship betwen toxin concentration and growth rate
toxin_conc_where_effect_saturates = 1. # threshold concentration where the effect ends
add_signal_parameter = 'bounded_linear'
S.add_signal(biomass_id, toxin_exch_id, 'ub', add_signal_parameter,
parms = [gr_absense_of_toxin,
toxin_conc_where_effect_starts,
slope_of_toxin_effect,
toxin_conc_where_effect_saturates])
Using these arguments will create a signal such that the growth rate of the susceptible model will reduce linearly in proportion to the concentration of toxin in the environment, with the effect on growth rate occurring whenever toxin is present and saturating when the toxin concentration reaches 1.