With a functional Ensemble in hand, you're ready to perform simulations. In medusa, most simulations are performed by setting the model structure to represent an individual member, using cobrapy functions for the actual simulation, then repeating for all or many ensemble members.
Flux balance analysis (FBA) is one of the most widely used techniques in systems biology. See What is flux balance analysis? for an introduction to FBA, and the cobrapy documentation to see how FBA is performed with a single model.
When using medusa for FBA, the environmental conditions and objective function should be specified in ensemble.base_model, just as if it were a normal cobrapy Model:
In [1]:
import medusa
from medusa.test import create_test_ensemble
ensemble = create_test_ensemble("Staphylococcus aureus")
In [2]:
ensemble.base_model.objective.expression
Out[2]:
The current objective function is the biomass reaction (bio1)--to change this, just set the objective to another reaction. Let's change the objective to CO2 exchange, then change it back to biomass production:
In [3]:
ensemble.base_model.objective = 'EX_cpd00011_e'
print(ensemble.base_model.objective.expression)
ensemble.base_model.objective = 'bio1'
print(ensemble.base_model.objective.expression)
Similarly, you can manipulate the environmental conditions as in cobrapy. The base model for this example ensemble is from ModelSEED, so exchange reactions are specified with the 'EX_' prefix, followed by the metabolite id. Let's take a look at the exchange reactions that are currently open:
In [4]:
medium = ensemble.base_model.medium
medium
Out[4]:
That's a lot of open exchange reactions! Let's make them a bit more realistic for an in vitro situation. We'll load a file specifying the base composition of the media in biolog single C/N growth conditions, and set the media conditions to reflect that. The base composition is missing a carbon source, so we'll enable uptake of glucose. In the medium dictionary, the numbers for each exchange reaction are uptake rates. If you inspect the actual exchange reactions, you will find that the equivalent to an uptake rate of 1000 units is a lower bound of -1000, because our exchange reactions are defined with the boundary metabolite as the reactant, e.g. cpd00182_e -->.
In [5]:
import pandas as pd
biolog_base = pd.read_csv("../medusa/test/data/biolog_base_composition.csv", sep=",")
biolog_base
Out[5]:
In [6]:
# convert the biolog base to a dictionary, which we can use to set ensemble.base_model.medium directly.
biolog_base = {'EX_'+component:1000 for component in biolog_base['ID']}
# add glucose uptake to the new medium dictionary
biolog_base['EX_cpd00182_e'] = 10
# Set the medium on the base model
ensemble.base_model.medium = biolog_base
ensemble.base_model.medium
Out[6]:
With the medium set, we can now simulate growth in these conditions:
In [7]:
from medusa.flux_analysis import flux_balance
fluxes = flux_balance.optimize_ensemble(ensemble,return_flux='bio1')
In [16]:
# get fluxes for the first 10 members
fluxes.head(10)
Out[16]:
In [10]:
import matplotlib.pylab as plt
fig, ax = plt.subplots()
plt.hist(fluxes['bio1'])
ax.set_ylabel('# ensemble members')
ax.set_xlabel('Flux through biomass reaction')
plt.show()
You may want to perform simulations with only a subset of ensemble members. There are two options for this; either identifying the desired members for simulation with the specific_models parameter, or passing a number of random members to perform simulations with the num_models parameter.
In [14]:
# perform FBA with a random set of 10 members:
subflux = flux_balance.optimize_ensemble(ensemble, num_models = 10, return_flux = "bio1")
subflux
Out[14]:
In [15]:
submembers = [member.id for member in ensemble.members[0:10]]
print(submembers)
subflux = flux_balance.optimize_ensemble(ensemble, specific_models = submembers, return_flux = "bio1")
subflux
Out[15]:
In [ ]:
In [ ]: