BACKGROUND
metagenomic sequencing facilitates the study of a large number of microorganisms in environmental samples  <cit> .
this technique has been used to study the composition of gut microbiota  <cit> , its role in human metabolism  <cit>  and its relation to diseases such as atherosclerosis  <cit> , obesity  <cit>  and crohn’s disease  <cit> .
in functional metagenomic studies, it is common to use kegg orthology   <cit>  to annotate gene functions  <cit> .
ko can be used to predict the composition ratio of microbial gene families and pathways from the human microbiome project  <cit> .
the functional annotation for a large number of sequenced bacteria, nearly  <dig>  bacterial genomes, is stored in the integrated microbial genomes  database, and the genomes are mapped to kegg pathway images  <cit> .
this tremendous new genomic information adds a new opportunity for systems biology, as it enables use of information about genome content for prediction of metabolic phenotypes of species in the gut  <cit> , or to develop community systems  <cit>  or supra-model organisms  <cit> .
therefore, it is relevant to reconstruct accurate genome scale metabolic models  from ko annotated by metagenomic analysis.
several methods have been developed to reconstruct genome scale models from gems of closely related organisms  <cit> , kegg  <cit> , and the model seed  <cit> .
the raven toolbox  <cit>  has been used to generate gems for the eukaryotic microorganisms pichia stipitis and pichia pastoris using iin <dig>  a gem of saccharomyces cerevisiae <cit> .
however, this method requires a gem of a closely related organism.
the raven toolbox has another function to solve this problem by assigning gene to ko using muscle  <cit>  and hmmer  <cit> .
then it generates the draft model by mapping ko to kegg reactions.
the web-based methods, fame  <cit>  and microbesflux  <cit> , are able to produce draft models for ~ <dig> and ~ <dig>  kegg genomes, respectively.
the disadvantage of both fame and microbesflux is that they are limited to organisms already annotated in kegg.
the model seed can generate a draft model for a desired organism based on rast annotation of genes  <cit> .
even though some of these methods have computational gap filling methods, there is still a need for manual curation to obtain a functional model.
manual curation is generally cumbersome and time-consuming.
the lack of visualization, such as organization and readability of reactions and genes names into the model excel file and kegg maps, can hamper manual curation of generated draft models.
to facilitate generation of a draft model and manual gap filling, we assembled an organized reference reaction data set consisting of common microbial reactions, where every reaction is assigned with kos.
the reactions were collected from high quality gems and from rhea, a manually annotated database of chemical reactions  <cit> , but not from kegg reactions.
in spite of the accuracy of kegg reactions, some reactions need to substantially manual curation for substrate and co-factor usage, and the reactions in reconstructed gems generally have to be well annotated in terms of substrate and co-factor usage and elemental balancing.
our reaction data set can, in principle, be used to generate draft models for all  <dig>  bacteria in the img database, kegg organisms, or other user-defined organisms annotated by ko.
here, we used this reaction data set to generate high quality gems for two bacterial genomes from the img database: b. adolescentis l2- <dig> and f. prausnitzii a2- <dig>  bifidobacterium is a dominating genus in the phyla actinobacteria present in the human gut microbiota and faecalibacterium is the most abundant genus among the firmicutes.
firmicutes, bacteroidetes and actinobacteria are the most highly abundant phyla in the human gut microbiota  <cit> .
both bifidobacterium and faecalibacterium interact with bacteroidetes <cit> .
moreover, bifidobacterium produces acetate to protect the host from infection  <cit> , and faecalibacterium has a relation with crohn disease  <cit> .
furthermore, the production of butyrate by faecalibacterium, among others, has been associated with a healthy state  <cit> .
finally, we simulated and compared the interactions between the two organisms using two approaches for community modeling: flux balance analysis  <cit>  and optcom  <cit> .
methods
figure  <dig> shows a summary of the methods developed and employed in this study.
reconstruction was based on the assignment of genes to kos and this was used to generate a draft model from the reaction data set.
a kegg map viewer is helpful in identifying gaps in the reconstruction.
model reconstructions
step  <dig>
we assembled a reference reaction data set containing reactions assigned with ko for each kegg map .
for organisms escherichia coli k- <dig> mg <dig>  staphylococcus aureus n <dig> and saccharomyces cerevisiae, we mapped the genes from each kegg map to the corresponding reactions contained in the respective gems iaf <dig>  <cit> , isb <dig>  <cit>  and ito <dig>  <cit> .
the gem reactions, together with the associated kos, were then added to the reference reaction data set.
in the case the organism had no genes for the enzyme or all gems lacked the reaction, we downloaded it from rhea by the matlab function get_reaction_from_rhea.
this function downloads rhea reactions as xml format and prints them.
to avoid the mismatched metabolite names between rhea reactions and the reference reaction set, we evaluated the metabolite names in the rhea reaction using its corresponding kegg reaction.
we retained rhea metabolite names if they did not exist in the reference reaction data set, and the rhea reaction and its corresponding ko were subsequently inserted into the reference reaction data set.
iaf <dig> reactions not present in any of the kegg maps, especially transporters and reactions occurring in the cell wall, were inserted into the reference reaction data set with the corresponding ko as the corresponding gene.
moreover, we added the capsular polysaccharide and teichoic acid biosynthesis reactions, which are required for cell wall biosynthesis, from the gem of lactobacillus plantarum wcfs  <cit> , because cell wall could be a significant fraction of gram dry weight of gram-positive bacteria.
also we added the methane metabolism from the gem of methanosarcina barkeri <cit>  and added the siderophore group biosynthesis from the gem of mycobacterium tuberculosis <cit> .
finally, the reactions were organized, ordered and made readable to facilitate the manual gap filling process.
step <dig>
we downloaded the gene annotation for each studied organism from the img database  <cit> .
we extracted the set gk=  for each organism from the downloaded img file using the function get_gene_ko_from_img.
for organisms available in kegg, the set gk can be obtained directly using the function get_gene_ko_from_kegg_org_id, otherwise the users can build the set gk themselves.
the set gk was passed to the function builddraftmodel to extract reactions using ko identifiers from the reaction data set.
the draft model was exported to an excel file by the function savedraftmodel.
finally, we removed the exchange reaction for metabolites that were not participating in any cytosolic reaction.
all the described functions are provided in additional file  <dig> and can be used in the raven toolbox.
step  <dig>
the gaps in each model were filled manually by mapping the model to kegg maps and inserting the required reactions to ensure full connectivity in the model.
to find genes for the filled reactions or metabolic genes, we extended the search to other available gene annotation in the img database.
both studied organisms have genes annotated by pfam  <cit> , tigrfams  <cit> , tc families  <cit>  and metacyc  <cit> .
moreover, b. adolescentis l2- <dig> has gene annotations by seed.
annotation with tc families has no specific gene assignment, so we ran a bidirectional blast between the tc protein sequence and each organism sequence.
also, we searched for the residual filled reactions by their enzyme name or ec number in the pfam database to get the corresponding pfam identifiers, and the corresponding genes were searched in img file.
additional file 3: table s <dig> contains the results of this analysis.
flux balance analysis
flux balance analysis   <cit>  was used to reconstruct and validate both models using the raven toolbox and mosek  as a linear programming solver under the matlab programming environment .
equation  describes the main formulation of fba, where s is a stoichiometric matrix, v is a vector of flux values for all reactions and c is a weight vector for each flux in vector v. typically all values in c are zero except the flux of the biomass reaction  which is fixed to one.
ub is the upper bound for the flux, and lb is the lower bound for the flux.
in fba, we assumed that the model grows in a growth medium  which represents a set of metabolites d, such as a carbon source, ammonia, phosphate, sulfur.
certain exchange reactions, which carry flux from the medium to the model , were fixed to the model uptake rate g
d
, for the metabolite d. for example 1 mmol/gdw/h of glucose uptake rate is fixed to flux  which carries glucose from the medium to the model.
maxz=ctvs.t.s.v=0v≤ub−v≤−lbvd=gd,∀d∈gm
both b. adolescentis l2- <dig> and f. prausnitzii a2- <dig> models grow anaerobically in rich media.
we assumed that each model consumed ammonia as a source of nitrogen, phosphate, h2s or cysteine as a source of sulfur, nicotinate and all amino acids having transporter reactions.
in addition, the f. prausnitzii a2- <dig> model consumed folic acid from the medium.
xanthine, uracil and urea transporters were closed.
the compositions of protein, rna and dna in the biomass were estimated from neidhardt et al.
<cit> .
the compositions of peptidoglycan and capsular polysaccharide are the same as in the gem of lactobacillus plantarum wcfs  <cit> .
we used flux variability analysis  <cit>  to evaluate the predicted fluxes by the cobra toolbox function fluxvariability <cit>  .
to determine if a studied organism is able to grow on different carbon sources such as galactose, xylose or fructose, the transporter flux was fixed to 1 mmol/gdw/h and the biomass formation was optimized.
community modeling
in community studies each model was allowed to take up as much glucose as possible to maximize growth.
in fba simulations, g
d
splits into two kinds of reactions with the first being a distribution of metabolite d to the different microorganisms in the community and the second being different transport reactions, where each reaction represents transport of glucose to one organism  <cit> .
in optcom, g
d
became two variables: uptake variable for growth medium metabolite d, uval
d
k
, or/and export variable for growth media, eval
d
k
, where k = bif or fap   <cit>  and where bif is b. adolescentis l2- <dig> and fap is f. prausnitzii a2- <dig>  equation  gives the general optcom problem formulation for this community.
this problem has nonlinear constraint, so it cannot be solved using mosek.
we therefore wrote a function generateoptcommodel to convert raven matlab models to gams language  and used the baron solver hosted on the neos servers  <cit> , a free optimization server, to solve the optimization problem.
maxz=vbiomassbif+vbiomassfaps.t.maxzbif=vbiomassbif,maxzfap=vbiomassfaps.t.s.t.sbifvbif=0sfapvfap=0vbif≤ubbifvfap≤ubfap−vbif≤−lbbif−vfap≤−lbfapvglcbif≤uvalglcbifvglcfap≤uvalglcfapvacetatebif=evalacetatebifvacetatefap≤uvalacetatefapuvalglcbif+uvalglcfap=total_glcuvalacetatefap=evalacetatebif
description of community using xml
we described a community structure without details of each model as xml format , because libsbml fails to read an sbml containing a user defined attribute or xml tag for community features  <cit> .
both ibif <dig> and ifap <dig> competed for glucose while ifap <dig> consumed acetate produced by ibif <dig>  the functions generateoptcommodel and generatecommodel used xml files to generate optcom and fba models.
RESULTS
reconstruction of reference reaction data set
gems elucidate how organisms consume nutrients, carbon source, ammonia, phosphate and autotrophic metabolites to build their biomass precursors and produce chemical byproducts  <cit> .
the biochemical reactions included in a gem are based on experimental or predicted function of enzymes contained by the studied organism  <cit> .
reconstructed gems share many components: exchange flux, transport, central metabolism, nucleotide, amino acids, cofactor biosynthesis, cell wall and lipid.
most gems use kegg maps and literature to illustrate the content of each component, where it contains one or more kegg maps.
additional file 1: figure s <dig>  shows these components and the kegg map name for each component and how flux distributes to each component and builds the necessary biomass precursors.
to cover the reactions in additional file 1: figure s <dig>  we built a reference reaction data set from published gems and a manually curated reaction database rhea.
this reaction data set contained reactions for the central carbon metabolism , amino acids and nucleotide biosynthesis, cell wall  and cofactors .
we adopted the fatty acids biosynthesis and glycerophospholipd metabolism in iaf <dig>  then we included reactions that connect other carbon resources, such as galactose and maltose to the main network.
the reactions were organized to facilitate manual revision and editing of newly reconstructed gems.
each reaction was assigned with kos obtained from kegg maps.
the reference reaction data set comprises  <dig>  reactions out of which  <dig> came from rhea,  <dig> unique metabolite and  <dig> kos .
the reference reaction data set was used to generate a draft model with an input file containing a gene and its ko.
additional file 7: figure s <dig> shows how the reaction set covers kegg maps.
gems description
table  <dig> shows statistics for the reconstructed gems and comparison with draft models generated by model seed.
the ibif <dig> model comprises  <dig> reactions,  <dig> unique metabolites and  <dig> genes constituting  <dig> % of the total number of genes.
it contains  <dig> genes from tigrfams and one gene from pfam.
it needed  <dig> reactions to become a connected model.
model ifap <dig> comprises  <dig> reactions,  <dig> unique metabolites and  <dig> genes, constituting  <dig> % of the total number of genes.
it contains  <dig> genes from tigrfams and two genes from pfam.
it needed  <dig> reactions to become a connected model.
additional file  <dig> contains both models as excel and sbml format.
both ibif <dig> and ifap <dig> have reactions for the central carbon metabolism and can utilize other sole carbon sources than glucose, as reported in in vivo studies.
ibif <dig> can utilize galactose, fructose and maltose, which is consistent with in vivo studies  <cit> .
the ifap <dig> can also utilize galactose and maltose but cannot utilize xylose, which is also consistent with in vivo studies  <cit> .
ibif <dig> features the bifid shunt pathway or the f6ppk pathway representing a special bifidobacteria pathway converting glucose to pyruvate   <cit> .
the f6ppk pathway includes fructose-6-phosphate phosphoketolase converting d-fructose 6-phosphate to acetyl phosphate and d-erythrose 4-phosphate, compared to the common part of glycolysis with 6-phosphofructokinase and fructose- <dig> -bisphosphate aldolase.
ifap <dig> has a faecalibacterium prausnitzii butyrate producing pathway   <cit> .
neither ibif <dig> nor ifap <dig> produces anything when the glucose uptake rate is 0 mmol/gdw/h and the objective function is biomass or atp non-growth association maintenance, so the models did not generate energy or matter from nothing.
in spite of the two models were validated with fba in the following two sections, the comprehensive validation of the gems needs extensive experiments.
since these bacteria are not yet well-studied, we think the two models based on recently sequenced, may assist a lot to overcome some important questions, computing phenotypic states and describing the genotype-phenotype relationships.
bifidobacterium adolescentis adolescentis l2- <dig> validation
bifidobacterium has predicted genes for biosynthesis of all  <dig> amino acids, purines and pyrimidines  <cit> .
however, bifidobacterium only grows in complex media, probably because some of the genes in the biosynthetic pathway for amino acids are non-functional  <cit> .
we assumed that ibif <dig> grows in media containing  <dig> amino acids for which it has transporter reactions.
the ibif <dig> model did not produce lactate when biomass or atp production was optimized, while bifidobacterium produces acetate, lactate, formate and ethanol in vivo.
under glucose limitation, bifidobacterium does not produce lactate because it tries to maximize energy production by cleaving pyruvate to acetyl phosphate and formate  <cit> .
furthermore, the specific rate of sugar consumption affects the amount of lactate production.
for example, the organism produces a large amount of lactate when it has a rapid sugar consumption, but produces a small amount of lactate when it consumes a less preferred sugar like oligofructose  <cit> .
to study the ability of the model ibif <dig> to produce lactate, we maximized atp production for non-growth association maintenance, i.e., the reaction atp + h2o = > adp + phosphate + h. the model produced 3 mmol of atp per 1 mmol of glucose and produced only acetate, formate and ethanol.
when the model was constrained to produce 1 mmol of lactate per 1 mmol of glucose, it produced  <dig>  mmol of atp and  <dig>  mmol of acetate per 1 mmol of glucose.
figure  <dig> shows the effect of lactate production on the ibif <dig> model when the biomass is an objective function.
in figure 2a, the model achieved maximum biomass when there was no production of lactate.
in figure 2b, with increased production of lactate, the production rate of acetate, ethanol and formate decreased.
the model was still able to produce acetate without formate because bifidobacterium has the fructose-6-phosphate phosphoketolase enzyme, which converts fructose 6-phosphate into erythrose 4-phosphate and acetyl phosphate and the latter can yield an atp when metabolized to acetate.
flux variability analysis showed the differences between maximum and minimum fluxes for acetate, formate and ethanol were  <dig> ,  <dig>  and  <dig>  mmol per 1 mmol of glucose.
the last results showed that the model aimed to generate atp by converting pyruvate to acetate through acetyl-coa and acetyl-phosphate and it therefore has to regenerate nad + by forming ethanol, as this is only way this co-factor can be balanced when there is formation of acetate .
although the model predicts a flux distribution for the theoretical ratio between acetate and lactate in bifidobacterium, it fails to predict the amount of lactate just like previous gems of lactic acid bacteria  <cit> .
to overcome this problem, oliveira et al.
constrained the pyruvate formate lyase reaction to an interval to deal with lactate production in a gem of l. lactic <cit> .
bas teusink et al.
fixed the measured flux in the gem of l. plantarum wcfs <dig>  <cit> .
milan et al.
added new enzyme turnover parameter to avoid metabolism overflow in gem of l.lactis <cit>  based on flux balance analysis with molecular crowding  <cit> .
finally, bas teusink et al.
showed that l. plantarum optimizes its yield when it grows with glycerol to support the prediction of gem in lactic acid bacteria  <cit> .
in the present work, we constrained the lactate flux in the model.
when we constrain lactate production with a yield of  <dig> , the model produces acetate with a yield of  <dig>  , as listed in table  <dig>  flux variability analysis shows acetate yield having a  <dig>   difference between maximum and minimum fluxes.
when the atp production for non-growth association maintenance was used as an objective function, it produced acetate with a yield of  <dig>  .
flux variability analysis showed that the acetate yield was not different between the maximum and minimum fluxes.
finally, when the model was constrained to produce  <dig>  mol of lactate per mol of glucose, it had a growth yield of  <dig>  gdw per mol of glucose, which was very close to the in vivo.
bifidobacterium growth yield of  <dig>  gdw per mol of glucose  <cit> .
faecalibacterium prausnitzii a2- <dig> validation
to study the effect of external acetate on butyrate production in faecalibacterium prausnitzii a2- <dig>  biomass production was used as an objective function.
the model produces butyrate with a yield of  <dig>   and co-consumes  <dig>  mmol of acetate per mmol of glucose.
the ratio of acetate uptake to butyrate production was  <dig> %, which is close to the 85-90% observed in in vivo studies of f. prausnitzii <cit> .
flux variability analysis shows that the acetate, butyrate, and formate have a difference of  <dig> ,  <dig> , and  <dig>  mmol/gdw/h between maximum and minimum fluxes, respectively.
table  <dig> shows the comparison between these values with in vivo studies  <cit> , where f. prausnitzii consumed 10 mm of glucose and  <dig>  ±  <dig>  mm of acetate to produce  <dig>  ±  <dig>  mm of butyrate.
figure  <dig> shows a sensitivity analysis of the effect of acetate and glucose on ifap <dig> with biomass as an objective function.
the acetate and glucose uptake rate varied from  <dig> to 1 mmol/gdw/h.
ifap <dig> cannot grow without glucose and grows poorly without acetate.
this result is similar to in vivo studies of f. prausnitzii a2- <dig>  that shows poor growth without acetate in the medium  <cit>  and the f. prausnitzii strains l2- <dig> and atcc  <dig> cannot grow without acetate  <cit> .
the model predicted increase in growth rate with acetate supplied was  <dig> -fold, close to the in vivo value, where the increase in growth rate is  <dig> -fold  <cit> .
community simulation
both optcom and fba methods were applied to ibif <dig> and ifap <dig> to simulate how b. adolescentis l2- <dig> and f. prausnitzii a2- <dig> co-culture together.
both organisms compete for 1 mmol/gdw/h of glucose to maximize their growth.
the model ibif <dig> generates acetate and ifap <dig> consumes acetate to produce butyrate, which plays a critical role in colonic homeostasis and cancer prevention  <cit> .
figure  <dig> depicts the prediction of fluxes calculated using both methods: optcom and fba, where biomass was an objective function.
ifap <dig> consumed  <dig>  mmol/gdw/h of glucose and all the acetate produced by ibif <dig> to produce  <dig>  mmol/gdw/h of butyrate.
with flux variability analysis, the differences between maximum and minimum yield of acetate and butyrate were  <dig>  and  <dig>   respectively.
there was no significant difference between optcom and fba simulations.
however, for this simple system, fba provided faster simulation than optcom.
finally, fba was used to simulate how a consortium of ibif <dig> and ifap <dig> interact.
the composition of the consortia was varied from 0% to 100% ifap <dig>  while keeping the total biomass constant, and the objective function was minimization of glucose uptake rate in both ibif <dig> and ifap <dig>  the biomass production rate was fixed at × % in ifap <dig> and  % in ibif <dig> of total biomass, where the total biomass growth rate was  <dig>  .
figure  <dig> shows the increase in butyrate production as ifap <dig> composition increases.
CONCLUSIONS
we assembled a reaction set from published gems, where each reaction is assigned with ko.
this reaction set was used to generate draft gems for each non-kegg organisms.
it represented a simple method to generate bacterial draft models from kegg ko, instead of generating it from kgml  <cit> .
the description of a community as a xml format can be used together with the two community simulation methods optcom and fba.
this saves time and effort when performing community modeling.
community simulations of an acetate producer b. adolescentis adolescentis l2- <dig> and an acetate consumer f. prausnitzii a2- <dig> provided insights into metabolic cross talk between these two members of the gut microbiota.
it shows the importance of acetate supply to butyrate production, since the growth and production of faecalibacterium prausnitzii is severely hampered by limited acetate supply.
this is an initial attempt to approach the very complex ecosystem and metabolic organ that the gut microbiota constitutes.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
iee-l developed the methods, reconstructed the two genome-scale metabolic models, performed the simulations and drafted the manuscript.
fk and ss assisted in the process of model reconstruction and edited the manuscript.
in, ths and jn supervised the work and edited the manuscript.
jn conceived and designed the project.
all authors read and approved the final manuscript.
supplementary material
additional file 1: figure s1
mapping the main component of gem to kegg pathway maps .
click here for file
additional file 2
matlab funtions used in the paper.
click here for file
additional file 3: table s1
genes with annotations: tigrfams, kos, seed, tc number.
click here for file
additional file 4: table s2
flux variability analysis for prediction fluxes in the study.
click here for file
additional file 5
the structure of community as xml format.
click here for file
additional file 6: table s3
reference reaction data set.
click here for file
additional file 7: figure s2
mapping the assembled reaction set to kegg metabolic pathways maps.
click here for file
additional file 8
the two models in sbml and ms excel format.
the two models are also available through the human metabolic atlass http://www.metabolicatlas.com/downloads/micro.
click here for file
additional file 9: figure s3
main carbon metabolism in bifidobacterium adolescentis adolescentis l2- <dig>
click here for file
additional file 10: figure s4
main carbon metabolism in faecalibacterium prausnitzii a2- <dig>
click here for file
acknowledgements
we thank ali zomorrodi and prof. costas maranas for providing us with optcom implementation.
ibrahim e. el-semman appreciates the egyptian cultural affairs and missions sector for their financial support.
we acknowledge funding from the knut and alice wallenberg foundation.
we thank chalmers library for funding the open access charge.

 metagenomic sequencing facilitates the study of a large number of microorganisms in environmental samples. this technique has been used to study the composition of gut microbiota, its role in human metabolism and its relation to diseases such as atherosclerosis, obesity and crohn’s disease. in functional metagenomic studies, it is common to use kegg orthology to annotate gene functions. ko can be used to predict the composition ratio of microbial gene families and pathways from the human microbiome project. the functional annotation for a large number of sequenced bacteria, nearly bacterial genomes, is stored in the integrated microbial genomes database, and the genomes are mapped to kegg pathway images. this tremendous new genomic information adds a new opportunity for systems biology, as it enables use of information about genome content for prediction of metabolic phenotypes of species in the gut, or to develop community systems or supra-model organisms. therefore, it is relevant to reconstruct accurate genome scale metabolic models from ko annotated by metagenomic analysis. several methods have been developed to reconstruct genome scale models from gems of closely related organisms, kegg, and the model seed. the raven toolbox has been used to generate gems for the eukaryotic microorganisms pichia stipitis and pichia pastoris using iin a gem of saccharomyces cerevisiae. however, this method requires a gem of a closely related organism. the raven toolbox has another function to solve this problem by assigning gene to ko using muscle and hmmer. then it generates the draft model by mapping ko to kegg reactions. the web-based methods, fame and microbesflux, are able to produce draft models for ~ and ~ kegg genomes, respectively. the disadvantage of both fame and microbesflux is that they are limited to organisms already annotated in kegg. the model seed can generate a draft model for a desired organism based on rast annotation of genes. even though some of these methods have computational gap filling methods, there is still a need for manual curation to obtain a functional model. manual curation is generally cumbersome and time-consuming. the lack of visualization, such as organization and readability of reactions and genes names into the model excel file and kegg maps, can hamper manual curation of generated draft models. to facilitate generation of a draft model and manual gap filling, we assembled an organized reference reaction data set consisting of common microbial reactions, where every reaction is assigned with kos. the reactions were collected from high quality gems and from rhea, a manually annotated database of chemical reactions, but not from kegg reactions. in spite of the accuracy of kegg reactions, some reactions need to substantially manual curation for substrate and co-factor usage, and the reactions in reconstructed gems generally have to be well annotated in terms of substrate and co-factor usage and elemental balancing. our reaction data set can, in principle, be used to generate draft models for all bacteria in the img database, kegg organisms, or other user-defined organisms annotated by ko. here, we used this reaction data set to generate high quality gems for two bacterial genomes from the img database: b. adolescentis l2- and f. prausnitzii a2- bifidobacterium is a dominating genus in the phyla actinobacteria present in the human gut microbiota and faecalibacterium is the most abundant genus among the firmicutes. firmicutes, bacteroidetes and actinobacteria are the most highly abundant phyla in the human gut microbiota. both bifidobacterium and faecalibacterium interact with bacteroidetes. moreover, bifidobacterium produces acetate to protect the host from infection, and faecalibacterium has a relation with crohn disease. furthermore, the production of butyrate by faecalibacterium, among others, has been associated with a healthy state. finally, we simulated and compared the interactions between the two organisms using two approaches for community modeling: flux balance analysis and optcom. methods figure  shows a summary of the methods developed and employed in this study. reconstruction was based on the assignment of genes to kos and this was used to generate a draft model from the reaction data set. a kegg map viewer is helpful in identifying gaps in the reconstruction. model reconstructions step we assembled a reference reaction data set containing reactions assigned with ko for each kegg map. for organisms escherichia coli k- mg staphylococcus aureus n and saccharomyces cerevisiae, we mapped the genes from each kegg map to the corresponding reactions contained in the respective gems iaf, isb and ito. the gem reactions, together with the associated kos, were then added to the reference reaction data set. in the case the organism had no genes for the enzyme or all gems lacked the reaction, we downloaded it from rhea by the matlab function get_reaction_from_rhea. this function downloads rhea reactions as xml format and prints them. to avoid the mismatched metabolite names between rhea reactions and the reference reaction set, we evaluated the metabolite names in the rhea reaction using its corresponding kegg reaction. we retained rhea metabolite names if they did not exist in the reference reaction data set, and the rhea reaction and its corresponding ko were subsequently inserted into the reference reaction data set. iaf reactions not present in any of the kegg maps, especially transporters and reactions occurring in the cell wall, were inserted into the reference reaction data set with the corresponding ko as the corresponding gene. moreover, we added the capsular polysaccharide and teichoic acid biosynthesis reactions, which are required for cell wall biosynthesis, from the gem of lactobacillus plantarum wcfs, because cell wall could be a significant fraction of gram dry weight of gram-positive bacteria. also we added the methane metabolism from the gem of methanosarcina barkeri and added the siderophore group biosynthesis from the gem of mycobacterium tuberculosis. finally, the reactions were organized, ordered and made readable to facilitate the manual gap filling process. step we downloaded the gene annotation for each studied organism from the img database. we extracted the set gk= for each organism from the downloaded img file using the function get_gene_ko_from_img. for organisms available in kegg, the set gk can be obtained directly using the function get_gene_ko_from_kegg_org_id, otherwise the users can build the set gk themselves. the set gk was passed to the function builddraftmodel to extract reactions using ko identifiers from the reaction data set. the draft model was exported to an excel file by the function savedraftmodel. finally, we removed the exchange reaction for metabolites that were not participating in any cytosolic reaction. all the described functions are provided in additional file and can be used in the raven toolbox. step the gaps in each model were filled manually by mapping the model to kegg maps and inserting the required reactions to ensure full connectivity in the model. to find genes for the filled reactions or metabolic genes, we extended the search to other available gene annotation in the img database. both studied organisms have genes annotated by pfam, tigrfams, tc families and metacyc. moreover, b. adolescentis l2- has gene annotations by seed. annotation with tc families has no specific gene assignment, so we ran a bidirectional blast between the tc protein sequence and each organism sequence. also, we searched for the residual filled reactions by their enzyme name or ec number in the pfam database to get the corresponding pfam identifiers, and the corresponding genes were searched in img file. additional file 3: table s contains the results of this analysis. flux balance analysis flux balance analysis was used to reconstruct and validate both models using the raven toolbox and mosek as a linear programming solver under the matlab programming environment. equation describes the main formulation of fba, where s is a stoichiometric matrix, v is a vector of flux values for all reactions and c is a weight vector for each flux in vector v. typically all values in c are zero except the flux of the biomass reaction which is fixed to one. ub is the upper bound for the flux, and lb is the lower bound for the flux. in fba, we assumed that the model grows in a growth medium which represents a set of metabolites d, such as a carbon source, ammonia, phosphate, sulfur. certain exchange reactions, which carry flux from the medium to the model, were fixed to the model uptake rate g d, for the metabolite d. for example 1 mmol/gdw/h of glucose uptake rate is fixed to flux which carries glucose from the medium to the model. maxz=ctvs.t.s.v=0v≤ub−v≤−lbvd=gd,∀d∈gm both b. adolescentis l2- and f. prausnitzii a2- models grow anaerobically in rich media. we assumed that each model consumed ammonia as a source of nitrogen, phosphate, h2s or cysteine as a source of sulfur, nicotinate and all amino acids having transporter reactions. in addition, the f. prausnitzii a2- model consumed folic acid from the medium. xanthine, uracil and urea transporters were closed. the compositions of protein, rna and dna in the biomass were estimated from neidhardt et al.. the compositions of peptidoglycan and capsular polysaccharide are the same as in the gem of lactobacillus plantarum wcfs. we used flux variability analysis to evaluate the predicted fluxes by the cobra toolbox function fluxvariability. to determine if a studied organism is able to grow on different carbon sources such as galactose, xylose or fructose, the transporter flux was fixed to 1 mmol/gdw/h and the biomass formation was optimized. community modeling in community studies each model was allowed to take up as much glucose as possible to maximize growth. in fba simulations, g d splits into two kinds of reactions with the first being a distribution of metabolite d to the different microorganisms in the community and the second being different transport reactions, where each reaction represents transport of glucose to one organism. in optcom, g d became two variables: uptake variable for growth medium metabolite d, uval d k, or/and export variable for growth media, eval d k, where k = bif or fap  and where bif is b. adolescentis l2- and fap is f. prausnitzii a2- equation gives the general optcom problem formulation for this community. this problem has nonlinear constraint, so it cannot be solved using mosek. we therefore wrote a function generateoptcommodel to convert raven matlab models to gams language and used the baron solver hosted on the neos servers, a free optimization server, to solve the optimization problem. maxz=vbiomassbif+vbiomassfaps.t.maxzbif=vbiomassbif,maxzfap=vbiomassfaps.t.s.t.sbifvbif=0sfapvfap=0vbif≤ubbifvfap≤ubfap−vbif≤−lbbif−vfap≤−lbfapvglcbif≤uvalglcbifvglcfap≤uvalglcfapvacetatebif=evalacetatebifvacetatefap≤uvalacetatefapuvalglcbif+uvalglcfap=total_glcuvalacetatefap=evalacetatebif description of community using xml we described a community structure without details of each model as xml format, because libsbml fails to read an sbml containing a user defined attribute or xml tag for community features. both ibif and ifap competed for glucose while ifap consumed acetate produced by ibif the functions generateoptcommodel and generatecommodel used xml files to generate optcom and fba models.  reconstruction of reference reaction data set gems elucidate how organisms consume nutrients, carbon source, ammonia, phosphate and autotrophic metabolites to build their biomass precursors and produce chemical byproducts. the biochemical reactions included in a gem are based on experimental or predicted function of enzymes contained by the studied organism. reconstructed gems share many components: exchange flux, transport, central metabolism, nucleotide, amino acids, cofactor biosynthesis, cell wall and lipid. most gems use kegg maps and literature to illustrate the content of each component, where it contains one or more kegg maps. additional file 1: figure s shows these components and the kegg map name for each component and how flux distributes to each component and builds the necessary biomass precursors. to cover the reactions in additional file 1: figure s we built a reference reaction data set from published gems and a manually curated reaction database rhea. this reaction data set contained reactions for the central carbon metabolism, amino acids and nucleotide biosynthesis, cell wall and cofactors. we adopted the fatty acids biosynthesis and glycerophospholipd metabolism in iaf then we included reactions that connect other carbon resources, such as galactose and maltose to the main network. the reactions were organized to facilitate manual revision and editing of newly reconstructed gems. each reaction was assigned with kos obtained from kegg maps. the reference reaction data set comprises reactions out of which came from rhea, unique metabolite and kos. the reference reaction data set was used to generate a draft model with an input file containing a gene and its ko. additional file 7: figure s shows how the reaction set covers kegg maps. gems description table  shows statistics for the reconstructed gems and comparison with draft models generated by model seed. the ibif model comprises reactions, unique metabolites and genes constituting % of the total number of genes. it contains genes from tigrfams and one gene from pfam. it needed reactions to become a connected model. model ifap comprises reactions, unique metabolites and genes, constituting % of the total number of genes. it contains genes from tigrfams and two genes from pfam. it needed reactions to become a connected model. additional file contains both models as excel and sbml format. both ibif and ifap have reactions for the central carbon metabolism and can utilize other sole carbon sources than glucose, as reported in in vivo studies. ibif can utilize galactose, fructose and maltose, which is consistent with in vivo studies. the ifap can also utilize galactose and maltose but cannot utilize xylose, which is also consistent with in vivo studies. ibif features the bifid shunt pathway or the f6ppk pathway representing a special bifidobacteria pathway converting glucose to pyruvate. the f6ppk pathway includes fructose-6-phosphate phosphoketolase converting d-fructose 6-phosphate to acetyl phosphate and d-erythrose 4-phosphate, compared to the common part of glycolysis with 6-phosphofructokinase and fructose- -bisphosphate aldolase. ifap has a faecalibacterium prausnitzii butyrate producing pathway. neither ibif nor ifap produces anything when the glucose uptake rate is 0 mmol/gdw/h and the objective function is biomass or atp non-growth association maintenance, so the models did not generate energy or matter from nothing. in spite of the two models were validated with fba in the following two sections, the comprehensive validation of the gems needs extensive experiments. since these bacteria are not yet well-studied, we think the two models based on recently sequenced, may assist a lot to overcome some important questions, computing phenotypic states and describing the genotype-phenotype relationships. bifidobacterium adolescentis adolescentis l2- validation bifidobacterium has predicted genes for biosynthesis of all amino acids, purines and pyrimidines. however, bifidobacterium only grows in complex media, probably because some of the genes in the biosynthetic pathway for amino acids are non-functional. we assumed that ibif grows in media containing amino acids for which it has transporter reactions. the ibif model did not produce lactate when biomass or atp production was optimized, while bifidobacterium produces acetate, lactate, formate and ethanol in vivo. under glucose limitation, bifidobacterium does not produce lactate because it tries to maximize energy production by cleaving pyruvate to acetyl phosphate and formate. furthermore, the specific rate of sugar consumption affects the amount of lactate production. for example, the organism produces a large amount of lactate when it has a rapid sugar consumption, but produces a small amount of lactate when it consumes a less preferred sugar like oligofructose. to study the ability of the model ibif to produce lactate, we maximized atp production for non-growth association maintenance, i.e., the reaction atp + h2o = > adp + phosphate + h. the model produced 3 mmol of atp per 1 mmol of glucose and produced only acetate, formate and ethanol. when the model was constrained to produce 1 mmol of lactate per 1 mmol of glucose, it produced  mmol of atp and  mmol of acetate per 1 mmol of glucose. figure  shows the effect of lactate production on the ibif model when the biomass is an objective function. in figure 2a, the model achieved maximum biomass when there was no production of lactate. in figure 2b, with increased production of lactate, the production rate of acetate, ethanol and formate decreased. the model was still able to produce acetate without formate because bifidobacterium has the fructose-6-phosphate phosphoketolase enzyme, which converts fructose 6-phosphate into erythrose 4-phosphate and acetyl phosphate and the latter can yield an atp when metabolized to acetate. flux variability analysis showed the differences between maximum and minimum fluxes for acetate, formate and ethanol were, and  mmol per 1 mmol of glucose. the last results showed that the model aimed to generate atp by converting pyruvate to acetate through acetyl-coa and acetyl-phosphate and it therefore has to regenerate nad + by forming ethanol, as this is only way this co-factor can be balanced when there is formation of acetate. although the model predicts a flux distribution for the theoretical ratio between acetate and lactate in bifidobacterium, it fails to predict the amount of lactate just like previous gems of lactic acid bacteria. to overcome this problem, oliveira et al. constrained the pyruvate formate lyase reaction to an interval to deal with lactate production in a gem of l. lactic. bas teusink et al. fixed the measured flux in the gem of l. plantarum wcfs. milan et al. added new enzyme turnover parameter to avoid metabolism overflow in gem of l.lactis based on flux balance analysis with molecular crowding. finally, bas teusink et al. showed that l. plantarum optimizes its yield when it grows with glycerol to support the prediction of gem in lactic acid bacteria. in the present work, we constrained the lactate flux in the model. when we constrain lactate production with a yield of, the model produces acetate with a yield of, as listed in table  flux variability analysis shows acetate yield having a difference between maximum and minimum fluxes. when the atp production for non-growth association maintenance was used as an objective function, it produced acetate with a yield of. flux variability analysis showed that the acetate yield was not different between the maximum and minimum fluxes. finally, when the model was constrained to produce  mol of lactate per mol of glucose, it had a growth yield of gdw per mol of glucose, which was very close to the in vivo. bifidobacterium growth yield of gdw per mol of glucose. faecalibacterium prausnitzii a2- validation to study the effect of external acetate on butyrate production in faecalibacterium prausnitzii a2- biomass production was used as an objective function. the model produces butyrate with a yield of and co-consumes  mmol of acetate per mmol of glucose. the ratio of acetate uptake to butyrate production was %, which is close to the 85-90% observed in in vivo studies of f. prausnitzii. flux variability analysis shows that the acetate, butyrate, and formate have a difference of,, and  mmol/gdw/h between maximum and minimum fluxes, respectively. table  shows the comparison between these values with in vivo studies, where f. prausnitzii consumed 10 mm of glucose and  ±   mm of acetate to produce  ±   mm of butyrate. figure  shows a sensitivity analysis of the effect of acetate and glucose on ifap with biomass as an objective function. the acetate and glucose uptake rate varied from to 1 mmol/gdw/h. ifap cannot grow without glucose and grows poorly without acetate. this result is similar to in vivo studies of f. prausnitzii a2- that shows poor growth without acetate in the medium and the f. prausnitzii strains l2- and atcc cannot grow without acetate. the model predicted increase in growth rate with acetate supplied was -fold, close to the in vivo value, where the increase in growth rate is -fold. community simulation both optcom and fba methods were applied to ibif and ifap to simulate how b. adolescentis l2- and f. prausnitzii a2- co-culture together. both organisms compete for 1 mmol/gdw/h of glucose to maximize their growth. the model ibif generates acetate and ifap consumes acetate to produce butyrate, which plays a critical role in colonic homeostasis and cancer prevention. figure  depicts the prediction of fluxes calculated using both methods: optcom and fba, where biomass was an objective function. ifap consumed  mmol/gdw/h of glucose and all the acetate produced by ibif to produce  mmol/gdw/h of butyrate. with flux variability analysis, the differences between maximum and minimum yield of acetate and butyrate were and respectively. there was no significant difference between optcom and fba simulations. however, for this simple system, fba provided faster simulation than optcom. finally, fba was used to simulate how a consortium of ibif and ifap interact. the composition of the consortia was varied from 0% to 100% ifap while keeping the total biomass constant, and the objective function was minimization of glucose uptake rate in both ibif and ifap the biomass production rate was fixed at × % in ifap and % in ibif of total biomass, where the total biomass growth rate was. figure  shows the increase in butyrate production as ifap composition increases.  we assembled a reaction set from published gems, where each reaction is assigned with ko. this reaction set was used to generate draft gems for each non-kegg organisms. it represented a simple method to generate bacterial draft models from kegg ko, instead of generating it from kgml. the description of a community as a xml format can be used together with the two community simulation methods optcom and fba. this saves time and effort when performing community modeling. community simulations of an acetate producer b. adolescentis adolescentis l2- and an acetate consumer f. prausnitzii a2- provided insights into metabolic cross talk between these two members of the gut microbiota. it shows the importance of acetate supply to butyrate production, since the growth and production of faecalibacterium prausnitzii is severely hampered by limited acetate supply. this is an initial attempt to approach the very complex ecosystem and metabolic organ that the gut microbiota constitutes. competing interests the authors declare that they have no competing interests. authors’ contributions iee-l developed the methods, reconstructed the two genome-scale metabolic models, performed the simulations and drafted the manuscript. fk and ss assisted in the process of model reconstruction and edited the manuscript. in, ths and jn supervised the work and edited the manuscript. jn conceived and designed the project. all authors read and approved the final manuscript. supplementary material additional file 1: figure s1 mapping the main component of gem to kegg pathway maps. click here for file additional file 2 matlab funtions used in the paper. click here for file additional file 3: table s1 genes with annotations: tigrfams, kos, seed, tc number. click here for file additional file 4: table s2 flux variability analysis for prediction fluxes in the study. click here for file additional file 5 the structure of community as xml format. click here for file additional file 6: table s3 reference reaction data set. click here for file additional file 7: figure s2 mapping the assembled reaction set to kegg metabolic pathways maps. click here for file additional file 8 the two models in sbml and ms excel format. the two models are also available through the human metabolic atlass http://www.metabolicatlas.com/downloads/micro. click here for file additional file 9: figure s3 main carbon metabolism in bifidobacterium adolescentis adolescentis l2- click here for file additional file 10: figure s4 main carbon metabolism in faecalibacterium prausnitzii a2- click here for file acknowledgements we thank ali zomorrodi and prof. costas maranas for providing us with optcom implementation. ibrahim e. el-semman appreciates the egyptian cultural affairs and missions sector for their financial support. we acknowledge funding from the knut and alice wallenberg foundation. we thank chalmers library for funding the open access charge.
the assembled reaction set is a useful tool to generate bacterial draft models from kegg orthology. here we assembled a reaction data set containing reactions obtained from existing genome-scale metabolic models, where each reaction is assigned with kegg orthology. the reaction data set was then used to reconstruct two genome scale metabolic models for gut microorganisms available in the img database bifidobacterium adolescentis l2- which produces acetate during fermentation, and faecalibacterium prausnitzii a2- which consumes acetate and produces butyrate. this tremendous new genomic information adds new opportunities for systems biology to reconstruct accurate genome scale metabolic models. the integrated microbial genomes database contains functional annotation for nearly bacterial genomes. the gut microbiota plays an important role in human health and disease by acting as a metabolic organ. f. prausnitzii is less abundant in patients with crohn’s disease and has been suggested to play an anti-inflammatory role in the gut ecosystem. the f. prausnitzii model, ifap comprises reactions and unique metabolites. the consortium of ibif and ifap was applied to predict f. prausnitzii’s demand for acetate and production of butyrate which plays an essential role in colonic homeostasis and cancer prevention. we used optcom and flux balance analysis to simulate how both organisms interact. metagenomic sequencing has shown how dysbiosis in the gut microbiota is associated with human metabolic diseases such as obesity and diabetes. fast and accurate reconstruction of metabolic models for members of the gut microbiota, as well as methods to simulate a community of microorganisms, are therefore needed. bifidobacterium adolescentis l2-32faecalibacterium prausnitzii a2-165genome-scale metabolic modelmetabolic modeling of gut microbiota the b. adolescentis model, ibif comprises reactions and unique metabolites.
 the gut microbiota plays an important role in human health and disease by acting as a metabolic organ. metagenomic sequencing has shown how dysbiosis in the gut microbiota is associated with human metabolic diseases such as obesity and diabetes. modeling may assist to gain insight into the metabolic implication of an altered microbiota. fast and accurate reconstruction of metabolic models for members of the gut microbiota, as well as methods to simulate a community of microorganisms, are therefore needed. the integrated microbial genomes database contains functional annotation for nearly bacterial genomes. this tremendous new genomic information adds new opportunities for systems biology to reconstruct accurate genome scale metabolic models.  here we assembled a reaction data set containing reactions obtained from existing genome-scale metabolic models, where each reaction is assigned with kegg orthology. the reaction data set was then used to reconstruct two genome scale metabolic models for gut microorganisms available in the img database bifidobacterium adolescentis l2- which produces acetate during fermentation, and faecalibacterium prausnitzii a2- which consumes acetate and produces butyrate. f. prausnitzii is less abundant in patients with crohn’s disease and has been suggested to play an anti-inflammatory role in the gut ecosystem. the b. adolescentis model, ibif comprises reactions and unique metabolites. the f. prausnitzii model, ifap comprises reactions and unique metabolites. each model was validated with in vivo data. we used optcom and flux balance analysis to simulate how both organisms interact.  the consortium of ibif and ifap was applied to predict f. prausnitzii’s demand for acetate and production of butyrate which plays an essential role in colonic homeostasis and cancer prevention. the assembled reaction set is a useful tool to generate bacterial draft models from kegg orthology. bifidobacterium adolescentis l2-32faecalibacterium prausnitzii a2-165genome-scale metabolic modelmetabolic modeling of gut microbiota 