BACKGROUND
clostridium difficile is a gram-positive, spore-forming, anaerobic bacterium, which infects or colonizes various animal species.
clinical manifestations in humans range from asymptomatic colonization to mild diarrhea, pseudomembranous colitis, and death  <cit> .
infection by this bacterium is associated not only with significant patient morbidity and mortality, but also with a large economic burden for healthcare systems  <cit> .
the primary risk factor for development of c. difficile infection among hospitalized patients is antibiotic use, which promotes toxicogenic c. difficile strains to proliferate, produce toxins, and induce disease  <cit> .
infection by this bacterium is most commonly associated with antibiotics such as clindamycin and amoxicillin  <cit> .
current recommendations for treatment of c. difficile infection  call for other antibiotics, such as metronidazole for mild infection cases and vancomycin for more severe cases  <cit> .
the emergence of hypervirulent and antibiotic-resistant strains of this bacterium has motivated the search for novel methods of treating cdi.
one method involves searching the bacterial central metabolic pathways for drug targets to create the next generation of antibiotics  <cit> .
the quest to better understand this bacterium and identify novel drug targets against it can benefit vastly from a model of the genotype-phenotype relationship of its metabolism.
methods to model the genotype-phenotype relationship range from stochastic kinetic models  <cit>  to statistical bayesian networks  <cit> .
kinetic models are limited as extensive experimental data is required to determine the rate laws and kinetic parameters of biochemical reactions.
an alternative to kinetic models is metabolic modeling, which has been used to depict a range of cell types without the need for difficult-to-measure kinetic parameters  <cit> .
metabolic models have been able to predict cellular functions, such as cellular growth capabilities on various substrates, effect of gene knockouts at genome scale  <cit> , and adaptation of bacteria to changes in their environment  <cit> .
metabolic models require a well-curated genome-scale metabolic network of the cell.
such networks contain all the known metabolic reactions in an organism, along with the genes that encode each enzyme involved in a reaction.
the networks are constructed based on genome annotations, biochemical characterizations, and published literature on the target organism.
the different scopes of such networks include metabolism, regulation, signaling, and other cellular processes  <cit> .
despite the success of metabolic modeling in capturing large-scale biochemical networks, the approach is limited as it describes cellular phenotype simply in terms of biochemical reaction rates and is thereby disconnected from other biological processes that impact phenotype.
moreover, metabolic models cannot account for changes in the metabolism of the bacterium in response to different environmental conditions.
recent advances in the omic technologies, such as genomics , transcriptomics , and proteomics , have enabled quantitative monitoring of the abundance of biological molecules at various levels in a high-throughput manner.
integration of transcriptomic data has been shown to be effective in improving metabolic model predictions of cellular behavior in different environmental conditions  <cit> .
here we present an integrated model of the metabolism of c.difficile strain  <dig>  we expanded the network imltc806cdf  <cit>  with regards to various pathways, such as fatty acid, glycerolipid, and glycerophospholipid metabolism.
fatty acids are not only important components of bacterial cell membranes but they are also important intermediate metabolites in the production of vitamins, lipid a, and quorum sensing molecules  <cit> .
the metabolism of phospholipids is also of interest as these compounds have been found to be closely tied to the growth phase in bacteria such as bacillus subtilis  <cit> .
to bridge the gap between gene expression data and protein abundance, we accounted for the codon usage bias of the bacterium.
during translation of a mrna to a protein, the information contained in the form of nucleotide triplets  in the rna is decoded to derive the amino acid sequence of the resulting protein.
most amino acids are coded by two to six synonymous codons.
these codons, which code for the same amino acid, are surprisingly used differentially in protein-encoding sequences  <cit> .
the codon usage has been found to alter the translation time and the abundance of the resulting protein  <cit> .
to our knowledge, this is the first time codon usage has been quantified and incorporated into a genome-scale metabolic reconstruction.
we used the modified network and flux balance analysis  <cit>  to simulate the steady-state metabolism of the bacterium.
to understand the behavior of the bacterium in different environments, we integrated gene expression data.
we incorporated the codon usage of the bacterium to bridge the gap between gene expression levels and protein abundance in the model.
we then validated our metabolic models against the literature on the bacterium.
following this validation process, we used our models to identify potential drug targets.
essential genes have been previously proposed as potential therapeutic targets.
<cit> .
we propose an additional method of predicting therapeutically-relevant genes through metabolic pathway sensitivity analysis and calculation of flux control coefficients.
the choice of gene to target can be further refined by eliminating genes with a human homolog to reduce the off-target effects of the selected drug  <cit> .
methods
construction and validation of the metabolic model
icdf834: an expansion of the imltc806cdf network
in modifying the imltc806cdf network  <cit> , we consulted kegg  <cit>  and incorporated some of the output from the review and curation of the metacyc  <cit>  database for c. difficile, which was released on march  <dig>   <dig>  during curation, we manually considered the directionality and gene-reaction associations of each reaction in the existing network.
we also manually expanded the existing network according to the procedure specified by thiele et al.
in  <cit> .
we supported additions to the network with published literature on the bacterium.
for example, the fatty acid profile found in clostridium difficile is mostly dominated by c16: <dig>  c16: <dig>  c18: <dig>  and c18: <dig>  <cit> .
the major phospholipid types in this bacterium are phosphatidylglycerol analogs, with pg, pg, pg, pg constituting the majority of these species  <cit> .
our modified network icdf <dig> modifies and expands pathways concerning lipid metabolism in the existing network, such as those where compounds and reactions involved had been grouped together.
by expanding the metabolism of the bacterium, we can also account for the wide array of fatty acids c. difficile can metabolize from its environment.
this can provide important insights as many gram-positive bacteria have been found to be able to incorporate and metabolize extracellular fatty acids  <cit> .
when defining metabolic pathways in the expanded network, we used kegg pathway identifiers so to remain consistent with the conventions employed in imltc806cdf  <cit> .
the lipid component of the biomass equation of imltc806cdf had been obtained from the metabolic network of staphylococcus aureus  <cit> , where lipid compounds had been lumped together.
there is a paucity of analyses on the chemical content of c. difficile’s biomass.
therefore, upon increasing the granularity of the network, we assumed coefficients from the biomass equation of the ibsu <dig> metabolic network developed for bacillus subtilis, where these lumped lipid and teichoic acid species have been replaced by explicit species.
constraint-based reconstruction and modeling approach
one constraint-based method for simulating the metabolic steady-state of a cell is flux-balance analysis , which can be used to analyze the metabolic network solely on the basis of systemic mass-balance and reaction capacity constraints.
fba simulations have been able to capture microorganism growth, nutritional resource consumption, and waste-product secretion rates of various cell types  <cit> .
the first step of fba involves representing the metabolic network in the form of a numerical matrix s of size .
this matrix contains the stoichiometric coefficients of each of the m metabolites in the n different reactions.
in the matrix, each row represents one unique metabolite and each column represents one reaction.
the stoichiometric matrix helps enforce a mass balance constraint on the system.
the mass balance on the cell for i= <dig> …,m metabolites and j= <dig> …,n reactions constrains the metabolite concentrations x
i, as shown in eq.
<dig>  where v
j is the flux through reaction j.
<dig> dxidt=∑j=1nsijvj,i= <dig> …,m.
under the steady state assumption dxidt= <dig> ∀i, the total amount of any compound being produced equals the total amount being consumed:
<dig> ∑j=1nsijvj= <dig> i= <dig> …,m.
in most metabolic models, there are more reactions than there are compounds  <cit> .
because there are more unknown variables than equations , any v that satisfies eq.
<dig> is considered to be in the null space of s.
fba can be used to find and determine points within the solution space that are most representative of the biological system using linear programming methods.
studies have revealed that metabolic fluxes in microorganisms are best predicted by maximizing the cellular objectives of growth  <cit> .
to determine the point corresponding to the maximum growth rate within the constrained space, the objective function shown in eq.
<dig> was maximized:
<dig> f=ctv,
where c is a vector of weights and indicates how much each reaction flux contributes to the biomass objective function.
the maximum growth rate can be achieved by determining the flux distribution v that results in maximal biomass flux.
additional constraints can be added through the upper bound vju and the lower bound vjl for the flux v
j. these bounds mandate the minimum and maximum fluxes allowed for a certain reaction and further decrease the space of allowable flux distributions for the relevant system.
the mathematical representation of the metabolic reactions, the objective function, and the capacity constraints define a linear system as shown in eq.
<dig>
<dig> maxctvsubject tosv=0vjl≤vj≤vju,j= <dig> …,n.
the model fluxes are usually given units of mmol/gdw·h, where gdw is the dry weight of cell mass in grams and h is the reaction time in hours.
the bounds enforce thermodynamic constraints by dictating whether reactions are reversible or irreversible.
the lower and upper flux bounds were arbitrarily chosen to be - <dig> mmol/gdw·h and  <dig> mmol/gdw·h for reversible reactions.
for irreversible reactions, vjl was chosen to be  <dig> mmol/gdw·h and vju was set to  <dig> mmol/gdw·h.
for our analysis, we used the cobra toolbox  <dig>   <cit>  in matlab .
multi-objective optimization in metabolic models
one limitation of using only biomass as the objective is that goals in metabolism are often different and simultaneously competing so the scalar notion of “optimality” does not hold; examples of such trade-offs include maximizing energy production while minimizing protein costs  <cit> .
moreover, the biomass objective vector is usually perpendicular to one of the surfaces of the solution space of the fba problem.
consequently, biomass maximizing flux states are usually degenerate; there exist multiple flux distributions that yield the same maximal biomass value  <cit> .
to choose between the various flux distributions, additional criteria must be considered.
for these reasons, we modeled metabolism as a multiobjective phenomenon.
by modeling the metabolism of bacterium as a multi-objective problem, we address a conflict problem whereby maximizing one objective  might involve a trade-off in the other objective ; cells are thought to face a trade-off that is described by the set of pareto-optimal solutions.
we used a multi-objective optimization approach to address the z objectives, as shown in eq.
<dig>
<dig> maxf=,f <dig> …,fz)subject tosv=0vjl≤vj≤vju,j= <dig> …,n.
note that, without loss of generality, we assumed that all the functions have to be maximized since minimizing a function f is equivalent to maximizing −f.
various works have attempted to systematically evaluate the ability of different objectives functions to reliably predict intracellular flux  <cit> .
according to their findings, bacterial metabolism can be better described by the objective of maximization of biomass or atp production paired with the objective of minimization of intracellular flux  <cit> .
introducing the minimization of intracellular flux as a secondary objective allows for economic allocation of resources by the bacterium by selecting for metabolic routes that contain the fewest number of steps  <cit> .
thus, for our analyses we used maximization of biomass, along with minimization of intracellular flux as our objectives.
in a maximization multi-objective problem, a vector that is part of the feasible space is considered to be pareto-optimal if all other vectors have the same or a lower value for at least one of the objective functions.
therefore, a pareto-optimal solution is found when there exists no other feasible solution which would increase one objective without decreasing another objective.
the set of pareto-optimal solutions constitutes the pareto-optimal front  <cit> .
in the absence of additional information, no one pareto-optimal solution can be said to be better than the other; higher-level information is required to choose one of the solutions  <cit> .
as proposed by costanza et al.
<cit> , to solve this multi-objective optimization problem one can use bilevel linear programming coupled with evolutionary algorithms, namely stochastic optimization methods that simulate the process of natural evolution.
evolutionary algorithms are well suited to multi-objective problems because they can generate multiple pareto-optimal solutions after one run and can use recombination to make use of the similarities of solutions  <cit> .
the input to the evolutionary algorithm is a set of arrays, also called individuals, representing potential solutions to the problem.
these arrays are then ranked based on the values of their objective functions.
potential optimal solutions are generated by retaining the best individuals and by generating new individuals through the use of variation.
this process is continued until no further improvements are detected on the pareto front.
the population size and the number of populations used with this algorithm were  <dig> and  <dig>  respectively.
to solve the linear programs, we used the gurobi solver   <cit> .
to validate our choice of objectives, we conducted a genetic analysis using multi-objective optimization.
in this analysis, binary “knockout” vectors were created, with each containing a  <dig> in the location of a gene set to be off  <cit> .
this analysis allowed us to determine how the growth of the organism changes in different environments, when genes may be turned on or off.
robustness analysis
a facet of living organisms is their homeostasis, otherwise known as their ability to remain robust to external and internal perturbations within a certain range.
external perturbations include changes in temperature or food supply while internal perturbations include spontaneous mutations.
the robustness of biological systems is partly due to the presence of parallel metabolic pathways.
robustness represents the insensitivity of a system to changes in system parameters.
global robustness  analysis can be used to survey the parameter space to determine the region where the cell exhibits specific features.
more specifically, we perturbed the flux bounds of the metabolic model and observed the resulting effects on biomass production.
the perturbation function γ where γ applies noise σ, assumed to be gaussian, to the system ψ for the trial τ. as proposed in  <cit> , a robust trial is associated with a ρ of 1:
<dig> ρ= <dig> if|ϕ−ϕ|≤ε <dig> otherwise
where ε is the robustness threshold.
the gr was defined as the percentage of trials determined to be robust.
we arbitrarily defined ε to be 1% of the metric ϕ and we arbitrarily limited the noise to 1%.
incorporating transcriptomic and codon usage data in genome-scale models
to increase the reliability of the model, gene expression data was added to the fba framework .
to relate this gene expression data to protein abundance, codon usage bias data was also incorporated.
the translation rate of a codon is determined in part by the speed of diffusion of a translationally-competent trna to the ribosome.
because trnas are differentially abundant in the cell, codons pairing to high-abundance trnas are translated faster than those pairing to low-abundance trnas.
although synonymous codons produce the same amino acid sequence, they can alter the translation speed and the protein expression levels depending on the abundance of their associated trna  <cit> .
studies have revealed that a large codon bias generally resulted in higher protein expression levels  <cit> .
therefore, the inclusion of codon bias can help improve the metabolic model predictions by helping link gene expression levels to protein levels.
fig.
<dig> framework for modeling the metabolism of c.difficile.
the updated metabolic network of the bacterium was used to create a metabolic model that was assessed using sensitivity and robustness analyses.
integrating gene expression and codon usage data yielded context-specific metabolic models that were evaluated against biological rationale and found fit for clinical applications.
the augmented metabolic models were then used to identify potential therapeutic targets using gene essentiality analysis, posa, and flux control coefficient calculations
the codon usage table for c. difficile was obtained from the kazusa codon usage database  <cit> , which lists the frequency of different codons in the genome.
the weights for synonymous codons was determined as the ratio between the observed frequency of the codon k and the frequency of the most preferred synonymous codon for that amino acid:
<dig> wk=fkmax,wherek,m∈.
we obtained the mrna sequence associated with the  <dig> genes of c. difficile from uniprot  <cit> .
the counts of different codons were determined for each mrna sequence.
to obtain a measure of the codon bias, we calculated the codon adaptation index  for each gene.
the cai represents the relative adaptiveness of the codon usage of the relevant gene to the codon usage of highly expressed genes  <cit> .
the cai ranges from  <dig> to  <dig>  with a value of  <dig> indicating high expression and, by correlation, high abundance of the associated protein.
the cai represents the geometric mean of the weights corresponding to the codons in the sequence:
<dig> cai=e1l∑l=1lln),
where l is the number of codons in the genes and w
k is the weight associated with codon type k for lth codon along the length l of the gene.
because a large codon bias has been shown to result in higher protein expression levels, the gene expression data g
t for each gene t was scaled by cai such that genes with the low codon bias had lower expression gt′:
<dig> gt′=gt·.
each of the reactions in the metabolic model depends on a gene set, which is represented through the use of and/or operators.
in this formulation, if a gene set is composed of two genes and an and operator, both genes are required to carry out the corresponding reaction.
on the other hand, if two genes connected by or, one gene is sufficient in carrying out the reaction.
this formulation can be transformed to derive the gene set expression gse
j for gene set j of reaction j from the expression of individual genes gt′, which in our case has been scaled by their respective codon usage.
when two genes are connected through an and operator, the gene set expression for reaction i, g
i, is the minimum of the scaled expression of the individual genes t making up the gene set.
the gene set expression for two genes connected by an or operator is the sum of the scaled expression of the individual genes.
in each reaction of the model, to map the gene set expression into a specific condition of the model, we used the piecewise muliplicative function h and the associated h
j was adopted as a multiplicative factor for the flux bounds  <cit>  :
vjlh≤vj≤vjuh,  where
<dig> h=|)gsej−1|gsej−1|ifgsej∈ℝ+∖{1}1ifgsej= <dig>
the function h was chosen because at high mrna abundance, an increase in mrna abundance has been found to produce a relatively small increase in the protein synthesis rate.
on the other hand, at low mrna abundance, an increase in mrna abundance has been found to produce a large increase in the protein synthesis rate  <cit> .
finally, we validated our context-specific metabolic models by incorporating codon usage and differential gene expression data into our model.
we then compared trends in our models’ biomass predictions to literature on the bacterium.
prediction of therapeutic targets
essential gene analysis
for each gene in the model, essential gene analysis involved removing reactions catalyzed by the gene or by a complex involving that gene and then using fba  <cit>  to predict growth.
genes were considered essential if following their removal, the predicted maximum growth rate was zero.
the c. difficile r <dig> , for which gene essentiality data was available for comparison with our in silico results, had been grown on tryptone-glucose-yeast extract  broth.
to approximate this medium, we used the complex medium defined by larocque et al.
during essential gene analysis of imltc806cdf  <cit> .
pathway-oriented sensitivity analysis
the growing research attention on metabolic pathways, rather than on specific reactions, is motivated by novel methods that allow for a better understanding of the functionality of complex webs of metabolic reactions.
to date, much of the study of metabolic pathways, their crosstalks, and their role in the overall metabotype has been carried out with statistical and model-based approaches  <cit> .
sensitivity analysis is used to identify model inputs that have a large influence on the model outputs.
to find the metabolic pathways that have the largest effect on the outputs of imltc806cdf and icdf <dig>  we used pathway-oriented sensitivity analysis   <cit> .
posa involves genetically manipulating the metabolic model to find the sensitive pathways, which make a large impact on model outputs.
in other words, we perturbed pathways by mutating the genes that govern their biochemical reactions and analyzed the result on the outputs.
in the knockout vector y={b
<dig> b
<dig> …,b
s,…,b
p}, b
s represents the perturbations on the genes governing the metabolic pathway s, where |b
s|=w
s .
because the gene knockouts are represented through the use of binary variables, we perform combinatorial perturbations, namely the bits in b
s are switched from  <dig> to  <dig> or from  <dig> to 0; note that if a gene in b
s is set to  <dig>  this gene is knocked-out in the model.
according to  <cit> , the pathway elementary effect  for the genetic perturbation b
s can be defined as follows:
<dig> pees=∥f−f∥Δs,
where bs~ represents the genetic manipulation of the input b
s; y~ is the mutation carried out on the knockout vector y; f is the vector v of fluxes as produced by the model; finally, Δ
s is a scale factor defined as:
<dig> Δs=1ws∑i=1wsbs~,s= <dig> …,p.
next, the sensitivity indices μ and σ are determined by calculating the mean and the standard deviation of the distribution of the pee for each input.
pathways with a large μ have a large influence on the output.
a large σ indicates an input whose influence highly depends on the value of other inputs.
by perturbing the genes through the use of knockouts and comparing the outputs of the model with and without the genetic manipulations, we detected the most sensitive pathways of the metabolic models.
calculation of flux control coefficients
posa provides valuable information on sensitive pathways that can be targeted by therapies, but often more specific drug target predictions are desired.
to understand how a metabolic pathway is controlled and can be altered, its control structure has to be determined.
the flux control coefficient  <cit>  is the flux v
ydh through a particular reaction, catalyzed by enzyme ydh, of the metabolic pathway with respect to the concentration x
xase of an enzyme xase:
<dig> cxxasevydh=∂vydh∂xxase·xxasevydh=∂lnvydh∂lnxxase
in our calculations, the enzyme concentration was assumed to be equal to the gene expression level adjusted by cai.
when calculating the flux control coefficients, we considered a 1% perturbation in the enzyme concentration x
xase.
flux control coefficients provide a quantitative measure of the degree of control an enzyme exerts on a metabolite flux and can quantitatively substitute for the qualitative concept of essential gene  <cit> .
thus, they can be used to identify steps that should be modified to achieve a successful alteration of the flux in outputs of clinical  relevance.
analysis of cdna microarrays
we used microarray analysis to determine the combination of genes which were up-regulated or down-regulated in different environmental conditions.
we used limma  <cit> , a package in bioconductor  <dig> , for statistical analysis of gene expression.
we preprocessed the data through background correction, within-array normalization, and between-array normalization.
after normalization, we used filtering to remove probes that did not appear to be expressed in any of the experimental conditions.
next, we used linear models to analyze the microarray data.
to conduct statistical analysis and assess differential expression, we used an empirical bayes method to modulate the standard errors of the log-fold changes.
to test for the comparisons of interest, we used an analysis of variance  model.
RESULTS
expansion and modification of imlt806cdf to icdf834
the genome of c. difficile strain  <dig> is composed of a circular chromosome of  <dig> , <dig> bp coding for  <dig> open reading frames , along with a plasmid containing  <dig> bp coding for  <dig> orfs  <cit> .
the modified metabolic network draft contains 21% of the orfs present in the chromosomal genome of the bacteria with  <dig> orfs, a modest improvement upon imltc806cdf, which contains  <dig> orfs, as shown in table  <dig>  our expanded metabolic network also consists of  <dig> metabolites and  <dig> reactions.
the final version of the network is available as an sbml file and as an excel file that indicate the reactions, metabolites, genes, and compartments involved in the metabolic network, along with references to literature that support additions or modifications to the existing network.
the new network has two additional dead-end metabolites as compared with those found in imltc806cdf.
the excel and sbml file, along with the the justification for keeping the dead-end metabolites in the model, have been uploaded to http://github.com/ssahebkashaf/peptoclostridiumdifficile <dig>  the code for all of the analyses employed in our work is also freely available on this repository.
genomic information of c. difficile
reconstructed models
imltc806cdf
icdf834
we repeated analyses previously conducted by larocque et al.
to validate imltc806cdf .
namely, we compared the ability of icdf <dig> and imltc806cdf to identify essential amino acids and metabolizable carbon sources.
the removal of amino acids that were not found to be essential or to affect growth, did not affect model-predicted biomass production in both models.
moreover, no biomass was produced in the absence of essential amino acids   <cit>  in both models.
therefore, similar to the previous network, our network is able to account for the essentiality of various amino acids on the growth of c.difficile.
with regards to carbon sources, both models were able to correctly predict a range of carbon sources that are utilized by the bacterium.
moreover, the bacterium was able to generate biomass in the absence of other carbon substrates, such as fructose, mannose, mannitol, and sorbitol.
this finding is consistent with literature, which maintains that c. difficile is not restricted to metabolizing sugars and can ferment other compounds, even amino acids, to obtain both its carbon and energy  <cit> .
validation of metabolic models
genetic analysis using multi-objective optimization
our modeling approach is intended to simulate the conflicting objectives faced by the bacterium, where optimal performance in one objective coincides with sub-optimal performance in another objective.
we used a knockout parameter space to find the genetic designs that would optimize the two objectives.
in fig.
<dig>  we show the areas of objective space discovered by the genetic algorithm during the genetic analysis from the first generation to generation  <dig>  the optimization algorithm adaptively moves to regions that maximize biomass while minimizing the total intracellular flux, as evident in the curvature of the plot in fig.
<dig>  after conducting the genetic analysis, the pareto front, shown in black in the inset of fig.
<dig>  was determined.
the pareto front is the set of nondominated solutions that represents the range of phenotypes resulting from different trade-offs between the two objectives.
the presence of a pareto front, as opposed to a singular dominated solution, aligned with our a priori expectations regarding the metabolic plasticity inherent to the bacterium  <cit> .
our findings, along with previous literature on the choice of objectives, supported our choice of objectives to model c. difficile’s metabolism.
fig.
<dig> genetic analysis using multi-objective optimization.
regions of objective space explored by the optimization algorithm for the objectives of maximization of biomass and minimization of total intracellular flux.
solutions are represented by progressively warmer colors depending on the time step of the algorithm in which they had been adaptively generated from the initial point.
the pareto front is shown in black in the inset
robustness analysis
we gauged the robustness of our model by determining the change in the maximal biomass flux in response to different perturbations.
global robustness  analysis revealed that the biomass production was fully robust to perturbations for a flux bound perturbation  and a tolerance  of 1% .
the gr falls when σ is increased or ε is decreased.
this facet of the bacterium’s metabolism was biologically relevant as bacteria such as c. difficile are able to grow despite small fluctuations in their physical environment.
robustness analysis illustrated that the global behavior of our metabolic model matches our expectations from biological rationale and supported the use of our models to predict the behavior of the bacterium in different environments.
changes in c. difficile’s growth in different conditions
we obtained the relevant microarray datasets from the gene expression omnibus  database  <cit>  under the accession numbers gse <dig> and from the arrayexpress database  <cit>  under the accession numbers e-geod- <dig> and e-bugs- <dig>  context-specific models for c. difficile were generated by incorporating gene expression data obtained for the bacterium in different environmental conditions.
to improve the reliability of the model, we also integrated codon usage data.
model predictions of these context-specific models were compared to expectations about the organism’s behavior from literature.
previous work suggests that sub-mic concentrations of amoxicillin, metronidazole, and clindamycin slowed growth of toxigenic c. difficile as compared with the controls  <cit> .
to test these findings in silico, we incorporated gene expression levels of c. difficile in response to sub-mic levels of different antibiotics into our model.
as compared with the c. difficile grown on bhi broth, toxigenic strains of c. difficile grown on sub-inhibitory concentrations of antibiotics exhibited reductions in their biomass, with those grown on amoxicillin showing the smallest growth .
this finding is supported by literature  <cit>  that has shown that in vitro, amoxicillin is effective against c. difficile.
these findings have lead to speculations that in vivo, this antibiotic is effective aginst vegetative forms of the bacterium but not against c. difficile spores  <cit> .
another potential explanation is that this broad-spectrum antibiotic may impair the intestinal microflora in a way that supports proliferation of c. difficile.
the microarray data for each condition was obtained from the geo or arrayexpress databases, using the specified accession numbers.
the differential gene expression levels obtained from analysis of this microarray data was used to make a metabolic model for each condition.
these context-specific metabolic models were used to predict change in biomass production for each condition compared with the control of each microarray dataset
additionally, the decline in biomass production following heat shock from  <dig> to  <dig> °c shown in table  <dig> could be due to the general stress response employed by the bacterium.
the heat shock response of c. difficile has been found to be involve gene clusters homologous to e. coli heat-shock operons  <cit> .
the heat shock response in e. coli has been found to be associated with a decrease in central carbon metabolism and a decline in cellular growth  <cit> .
literature on related bacteria is thereby in agreement with the model’s prediction of a significant reduction in growth in c. difficile following the heat shock.
additionally, according to the work of dubois et al., the supplementation of 10mm cysteine to the medium did not affect c. difficile’s growth  <cit> .
after integrating the microarray data from their work, we found that our in silico findings agreed with their experimental results.
validation of the findings of our context-specific metabolic models against the literature on the bacterium showed that metabolic models allow for an enriched view of omic data and may be valuable tools for better understanding the behavior of c. difficile in different conditions.
prediction of therapeutic targets
gene essentiality analysis
essential genes have been cited as promising targets for development of new antimicrobials due to their importance for bacterial survival  <cit> .
using fba, we performed an in silico gene deletion study to predict potential essential genes that may lead to the identification of new drug targets.
this analysis had already been conducted for imltc806cdf based on a 5% threshold, and gene essentiality results had been compared to genes deemed essential for b. subtilis, for which this data had been available  <cit> .
we performed gene essentiality analysis for both imltc806cdf and icdf <dig> and validated our results using recently available literature on the essential genes of the c. difficile r <dig>   <cit> .
while imltc806cdf predicted  <dig> essential genes and had a  <dig> % accuracy in predicting gene essentiality, icdf <dig> predicted  <dig> essential genes and had a  <dig> % accuracy .
pathway-oriented sensitivity analysis and flux control coefficients
for our posa analysis, we chose the gene expression profile of the bacterium when grown on bhi broth.
each pathway was assessed through random perturbations of its reactions, and the average perturbation μ and the standard deviation σ were computed as a result.
we performed the pathway-based sensitivity analysis and identified sensitive pathways before and after modifying the metabolic model as shown in fig.
<dig>  the pathway with the largest μ, and thereby the greatest control on biomass production or growth in both imltc806cdf and icdf <dig> is the valine, leucine, and isoleucine metabolism pathway.
these three amino acids are essential to the bacterium and their metabolism was also expected to be essential.
the second most sensitive pathway is alanine, aspartate, and glutmate metabolism in imltc806cdf and glycolysis/gluconeogenesis in icdf <dig>  additional sensitive pathways in icdf <dig> include pyrimidine metabolism and pyruvate metabolism.
model findings suggest that therapies against infection may likely be more effective if they target key enzymes in these sensitive pathways.
fig.
<dig> posa was used to compare the most sensitive pathways of imltc806cdf and icdf <dig>  the imltc806cdf model is composed of  <dig> metabolic pathways and the icdf <dig> model is composed of  <dig> metabolic pathways.
biomass production is most sensitive to pathways with higher calculated μ
to find more specific therapeutic targets, flux control coefficients for enzymes on biomass production in the metabolic model were determined and compared for bhi broth , cysteine supplementation , and heat shock  gene expression data.
the four enzymes with largest flux control coefficients in each condition are shown in fig.
<dig>  while the complete list of flux control coefficients in different conditions has been uploaded to the public repository.
these flux control coefficients were interestingly involved in pathways deemed sensitive during posa.
these enzymes varied amongst the four conditions, suggesting that access to the in vivo gene expression profile of c. difficile can be used to predict better drug targets for patients.
therapies aimed at reducing growth of c. difficile should target enzymes with high flux coefficients as, according to our model, their activity is most closely tied to biomass production.
fig.
<dig> genes encoding the enzymes with the largest flux control coefficients for biomass production in different conditions .
table of metabolic pathway hosting the genes and of gene descriptions  <cit>  .
a flux control coefficient of  <dig> implies full control of the metabolite flux by the associated enzyme
CONCLUSIONS
in this study, we expanded the existing metabolic network for c. difficile and used it to create context-specific metabolic models of its metabolism that allow us to understand how the bacterium alters its metabolism depending on its environment.
to predict the bacterium’s behavior in different environmental conditions, the model was integrated with transcriptomic and codon usage data to generate reliable and context-specific metabolic flux distributions.
we validated the model by conducting robustness and sensitivity analyses.
we further assessed its predictive potential by comparing model predictions with published experimental data to gauge the consistency of model findings with the current knowledge of c. difficile’s metabolism.
through this literature-based validation, we found that the model is a valuable tool for qualitatively understanding the behavior of the bacterium in different settings.
the model can also be used to find potential therapeutic targets by allowing for determination of essential genes and context-specific sensitive pathways and flux control coefficients.
context-specific metabolic models can allow for a better understanding different medically-relevant conditions  and can be continuously refined by integrating novel information regarding c. difficile’s metabolism.
our model can be used by biomedical researchers to study the bacterium and devise targeted treatments.
our approach can also be scaled-up to simulate the interactions between the gut microbiota and the host using a bottom-up modeling approach.
by accounting for the gut microbiota-host interactions, we can construct a whole gut model response to infections and other inflammatory events, paving the path towards more informed and effective treatments.
availability and requirements
project name:
icdf834
project home page:
http://github.com/ssahebkashaf/pep-toclostridiumdifficile630
operating system: platform independent
programming language: matlab
other requirements: gurobi
license: university of cambridge
any restrictions to use by non-academics: license needed
additional files
additional file  <dig> carbon source usage, carboxylic acids secreted, and essential amino acids in icdf <dig> as compared with imltc806cdf and literature.
additional file  <dig> gr of biomass production in icdf <dig> to internal and external perturbations.
additional file  <dig> gene essentiality analysis of icdf <dig> and imltc806cdf.
abbreviations
caicodon adaptation index
fbaflux balance analysis
grglobal robustness
orfopen reading frame
peepathway elementary effect
posapathway-oriented sensitivity analysis
tgytryptone-glucose-yeast extract

 clostridium difficile is a gram-positive, spore-forming, anaerobic bacterium, which infects or colonizes various animal species. clinical manifestations in humans range from asymptomatic colonization to mild diarrhea, pseudomembranous colitis, and death. infection by this bacterium is associated not only with significant patient morbidity and mortality, but also with a large economic burden for healthcare systems. the primary risk factor for development of c. difficile infection among hospitalized patients is antibiotic use, which promotes toxicogenic c. difficile strains to proliferate, produce toxins, and induce disease. infection by this bacterium is most commonly associated with antibiotics such as clindamycin and amoxicillin. current recommendations for treatment of c. difficile infection call for other antibiotics, such as metronidazole for mild infection cases and vancomycin for more severe cases. the emergence of hypervirulent and antibiotic-resistant strains of this bacterium has motivated the search for novel methods of treating cdi. one method involves searching the bacterial central metabolic pathways for drug targets to create the next generation of antibiotics. the quest to better understand this bacterium and identify novel drug targets against it can benefit vastly from a model of the genotype-phenotype relationship of its metabolism. methods to model the genotype-phenotype relationship range from stochastic kinetic models to statistical bayesian networks. kinetic models are limited as extensive experimental data is required to determine the rate laws and kinetic parameters of biochemical reactions. an alternative to kinetic models is metabolic modeling, which has been used to depict a range of cell types without the need for difficult-to-measure kinetic parameters. metabolic models have been able to predict cellular functions, such as cellular growth capabilities on various substrates, effect of gene knockouts at genome scale, and adaptation of bacteria to changes in their environment. metabolic models require a well-curated genome-scale metabolic network of the cell. such networks contain all the known metabolic reactions in an organism, along with the genes that encode each enzyme involved in a reaction. the networks are constructed based on genome annotations, biochemical characterizations, and published literature on the target organism. the different scopes of such networks include metabolism, regulation, signaling, and other cellular processes. despite the success of metabolic modeling in capturing large-scale biochemical networks, the approach is limited as it describes cellular phenotype simply in terms of biochemical reaction rates and is thereby disconnected from other biological processes that impact phenotype. moreover, metabolic models cannot account for changes in the metabolism of the bacterium in response to different environmental conditions. recent advances in the omic technologies, such as genomics, transcriptomics, and proteomics, have enabled quantitative monitoring of the abundance of biological molecules at various levels in a high-throughput manner. integration of transcriptomic data has been shown to be effective in improving metabolic model predictions of cellular behavior in different environmental conditions. here we present an integrated model of the metabolism of c.difficile strain we expanded the network imltc806cdf with regards to various pathways, such as fatty acid, glycerolipid, and glycerophospholipid metabolism. fatty acids are not only important components of bacterial cell membranes but they are also important intermediate metabolites in the production of vitamins, lipid a, and quorum sensing molecules. the metabolism of phospholipids is also of interest as these compounds have been found to be closely tied to the growth phase in bacteria such as bacillus subtilis. to bridge the gap between gene expression data and protein abundance, we accounted for the codon usage bias of the bacterium. during translation of a mrna to a protein, the information contained in the form of nucleotide triplets in the rna is decoded to derive the amino acid sequence of the resulting protein. most amino acids are coded by two to six synonymous codons. these codons, which code for the same amino acid, are surprisingly used differentially in protein-encoding sequences. the codon usage has been found to alter the translation time and the abundance of the resulting protein. to our knowledge, this is the first time codon usage has been quantified and incorporated into a genome-scale metabolic reconstruction. we used the modified network and flux balance analysis to simulate the steady-state metabolism of the bacterium. to understand the behavior of the bacterium in different environments, we integrated gene expression data. we incorporated the codon usage of the bacterium to bridge the gap between gene expression levels and protein abundance in the model. we then validated our metabolic models against the literature on the bacterium. following this validation process, we used our models to identify potential drug targets. essential genes have been previously proposed as potential therapeutic targets.. we propose an additional method of predicting therapeutically-relevant genes through metabolic pathway sensitivity analysis and calculation of flux control coefficients. the choice of gene to target can be further refined by eliminating genes with a human homolog to reduce the off-target effects of the selected drug. methods construction and validation of the metabolic model icdf834: an expansion of the imltc806cdf network in modifying the imltc806cdf network, we consulted kegg and incorporated some of the output from the review and curation of the metacyc database for c. difficile, which was released on march during curation, we manually considered the directionality and gene-reaction associations of each reaction in the existing network. we also manually expanded the existing network according to the procedure specified by thiele et al. in. we supported additions to the network with published literature on the bacterium. for example, the fatty acid profile found in clostridium difficile is mostly dominated by c16: c16: c18: and c18:. the major phospholipid types in this bacterium are phosphatidylglycerol analogs, with pg, pg, pg, pg constituting the majority of these species. our modified network icdf modifies and expands pathways concerning lipid metabolism in the existing network, such as those where compounds and reactions involved had been grouped together. by expanding the metabolism of the bacterium, we can also account for the wide array of fatty acids c. difficile can metabolize from its environment. this can provide important insights as many gram-positive bacteria have been found to be able to incorporate and metabolize extracellular fatty acids. when defining metabolic pathways in the expanded network, we used kegg pathway identifiers so to remain consistent with the conventions employed in imltc806cdf. the lipid component of the biomass equation of imltc806cdf had been obtained from the metabolic network of staphylococcus aureus, where lipid compounds had been lumped together. there is a paucity of analyses on the chemical content of c. difficile’s biomass. therefore, upon increasing the granularity of the network, we assumed coefficients from the biomass equation of the ibsu metabolic network developed for bacillus subtilis, where these lumped lipid and teichoic acid species have been replaced by explicit species. constraint-based reconstruction and modeling approach one constraint-based method for simulating the metabolic steady-state of a cell is flux-balance analysis, which can be used to analyze the metabolic network solely on the basis of systemic mass-balance and reaction capacity constraints. fba simulations have been able to capture microorganism growth, nutritional resource consumption, and waste-product secretion rates of various cell types. the first step of fba involves representing the metabolic network in the form of a numerical matrix s of size. this matrix contains the stoichiometric coefficients of each of the m metabolites in the n different reactions. in the matrix, each row represents one unique metabolite and each column represents one reaction. the stoichiometric matrix helps enforce a mass balance constraint on the system. the mass balance on the cell for i= …,m metabolites and j= …,n reactions constrains the metabolite concentrations x i, as shown in eq. where v j is the flux through reaction j. dxidt=∑j=1nsijvj,i= …,m. under the steady state assumption dxidt= ∀i, the total amount of any compound being produced equals the total amount being consumed: ∑j=1nsijvj= i= …,m. in most metabolic models, there are more reactions than there are compounds. because there are more unknown variables than equations, any v that satisfies eq. is considered to be in the null space of s. fba can be used to find and determine points within the solution space that are most representative of the biological system using linear programming methods. studies have revealed that metabolic fluxes in microorganisms are best predicted by maximizing the cellular objectives of growth. to determine the point corresponding to the maximum growth rate within the constrained space, the objective function shown in eq. was maximized: f=ctv, where c is a vector of weights and indicates how much each reaction flux contributes to the biomass objective function. the maximum growth rate can be achieved by determining the flux distribution v that results in maximal biomass flux. additional constraints can be added through the upper bound vju and the lower bound vjl for the flux v j. these bounds mandate the minimum and maximum fluxes allowed for a certain reaction and further decrease the space of allowable flux distributions for the relevant system. the mathematical representation of the metabolic reactions, the objective function, and the capacity constraints define a linear system as shown in eq. maxctvsubject tosv=0vjl≤vj≤vju,j= …,n. the model fluxes are usually given units of mmol/gdw·h, where gdw is the dry weight of cell mass in grams and h is the reaction time in hours. the bounds enforce thermodynamic constraints by dictating whether reactions are reversible or irreversible. the lower and upper flux bounds were arbitrarily chosen to be - mmol/gdw·h and mmol/gdw·h for reversible reactions. for irreversible reactions, vjl was chosen to be mmol/gdw·h and vju was set to mmol/gdw·h. for our analysis, we used the cobra toolbox in matlab. multi-objective optimization in metabolic models one limitation of using only biomass as the objective is that goals in metabolism are often different and simultaneously competing so the scalar notion of “optimality” does not hold; examples of such trade-offs include maximizing energy production while minimizing protein costs. moreover, the biomass objective vector is usually perpendicular to one of the surfaces of the solution space of the fba problem. consequently, biomass maximizing flux states are usually degenerate; there exist multiple flux distributions that yield the same maximal biomass value. to choose between the various flux distributions, additional criteria must be considered. for these reasons, we modeled metabolism as a multiobjective phenomenon. by modeling the metabolism of bacterium as a multi-objective problem, we address a conflict problem whereby maximizing one objective might involve a trade-off in the other objective; cells are thought to face a trade-off that is described by the set of pareto-optimal solutions. we used a multi-objective optimization approach to address the z objectives, as shown in eq. maxf=,f …,fz)subject tosv=0vjl≤vj≤vju,j= …,n. note that, without loss of generality, we assumed that all the functions have to be maximized since minimizing a function f is equivalent to maximizing −f. various works have attempted to systematically evaluate the ability of different objectives functions to reliably predict intracellular flux. according to their findings, bacterial metabolism can be better described by the objective of maximization of biomass or atp production paired with the objective of minimization of intracellular flux. introducing the minimization of intracellular flux as a secondary objective allows for economic allocation of resources by the bacterium by selecting for metabolic routes that contain the fewest number of steps. thus, for our analyses we used maximization of biomass, along with minimization of intracellular flux as our objectives. in a maximization multi-objective problem, a vector that is part of the feasible space is considered to be pareto-optimal if all other vectors have the same or a lower value for at least one of the objective functions. therefore, a pareto-optimal solution is found when there exists no other feasible solution which would increase one objective without decreasing another objective. the set of pareto-optimal solutions constitutes the pareto-optimal front. in the absence of additional information, no one pareto-optimal solution can be said to be better than the other; higher-level information is required to choose one of the solutions. as proposed by costanza et al., to solve this multi-objective optimization problem one can use bilevel linear programming coupled with evolutionary algorithms, namely stochastic optimization methods that simulate the process of natural evolution. evolutionary algorithms are well suited to multi-objective problems because they can generate multiple pareto-optimal solutions after one run and can use recombination to make use of the similarities of solutions. the input to the evolutionary algorithm is a set of arrays, also called individuals, representing potential solutions to the problem. these arrays are then ranked based on the values of their objective functions. potential optimal solutions are generated by retaining the best individuals and by generating new individuals through the use of variation. this process is continued until no further improvements are detected on the pareto front. the population size and the number of populations used with this algorithm were and respectively. to solve the linear programs, we used the gurobi solver. to validate our choice of objectives, we conducted a genetic analysis using multi-objective optimization. in this analysis, binary “knockout” vectors were created, with each containing a in the location of a gene set to be off. this analysis allowed us to determine how the growth of the organism changes in different environments, when genes may be turned on or off. robustness analysis a facet of living organisms is their homeostasis, otherwise known as their ability to remain robust to external and internal perturbations within a certain range. external perturbations include changes in temperature or food supply while internal perturbations include spontaneous mutations. the robustness of biological systems is partly due to the presence of parallel metabolic pathways. robustness represents the insensitivity of a system to changes in system parameters. global robustness analysis can be used to survey the parameter space to determine the region where the cell exhibits specific features. more specifically, we perturbed the flux bounds of the metabolic model and observed the resulting effects on biomass production. the perturbation function γ where γ applies noise σ, assumed to be gaussian, to the system ψ for the trial τ. as proposed in, a robust trial is associated with a ρ of 1: ρ= if|ϕ−ϕ|≤ε otherwise where ε is the robustness threshold. the gr was defined as the percentage of trials determined to be robust. we arbitrarily defined ε to be 1% of the metric ϕ and we arbitrarily limited the noise to 1%. incorporating transcriptomic and codon usage data in genome-scale models to increase the reliability of the model, gene expression data was added to the fba framework. to relate this gene expression data to protein abundance, codon usage bias data was also incorporated. the translation rate of a codon is determined in part by the speed of diffusion of a translationally-competent trna to the ribosome. because trnas are differentially abundant in the cell, codons pairing to high-abundance trnas are translated faster than those pairing to low-abundance trnas. although synonymous codons produce the same amino acid sequence, they can alter the translation speed and the protein expression levels depending on the abundance of their associated trna. studies have revealed that a large codon bias generally resulted in higher protein expression levels. therefore, the inclusion of codon bias can help improve the metabolic model predictions by helping link gene expression levels to protein levels. fig. framework for modeling the metabolism of c.difficile. the updated metabolic network of the bacterium was used to create a metabolic model that was assessed using sensitivity and robustness analyses. integrating gene expression and codon usage data yielded context-specific metabolic models that were evaluated against biological rationale and found fit for clinical applications. the augmented metabolic models were then used to identify potential therapeutic targets using gene essentiality analysis, posa, and flux control coefficient calculations the codon usage table for c. difficile was obtained from the kazusa codon usage database, which lists the frequency of different codons in the genome. the weights for synonymous codons was determined as the ratio between the observed frequency of the codon k and the frequency of the most preferred synonymous codon for that amino acid: wk=fkmax,wherek,m∈. we obtained the mrna sequence associated with the genes of c. difficile from uniprot. the counts of different codons were determined for each mrna sequence. to obtain a measure of the codon bias, we calculated the codon adaptation index for each gene. the cai represents the relative adaptiveness of the codon usage of the relevant gene to the codon usage of highly expressed genes. the cai ranges from to with a value of indicating high expression and, by correlation, high abundance of the associated protein. the cai represents the geometric mean of the weights corresponding to the codons in the sequence: cai=e1l∑l=1lln), where l is the number of codons in the genes and w k is the weight associated with codon type k for lth codon along the length l of the gene. because a large codon bias has been shown to result in higher protein expression levels, the gene expression data g t for each gene t was scaled by cai such that genes with the low codon bias had lower expression gt′: gt′=gt·. each of the reactions in the metabolic model depends on a gene set, which is represented through the use of and/or operators. in this formulation, if a gene set is composed of two genes and an and operator, both genes are required to carry out the corresponding reaction. on the other hand, if two genes connected by or, one gene is sufficient in carrying out the reaction. this formulation can be transformed to derive the gene set expression gse j for gene set j of reaction j from the expression of individual genes gt′, which in our case has been scaled by their respective codon usage. when two genes are connected through an and operator, the gene set expression for reaction i, g i, is the minimum of the scaled expression of the individual genes t making up the gene set. the gene set expression for two genes connected by an or operator is the sum of the scaled expression of the individual genes. in each reaction of the model, to map the gene set expression into a specific condition of the model, we used the piecewise muliplicative function h and the associated h j was adopted as a multiplicative factor for the flux bounds: vjlh≤vj≤vjuh, where h=|)gsej−1|gsej−1|ifgsej∈ℝ+∖{1}1ifgsej= the function h was chosen because at high mrna abundance, an increase in mrna abundance has been found to produce a relatively small increase in the protein synthesis rate. on the other hand, at low mrna abundance, an increase in mrna abundance has been found to produce a large increase in the protein synthesis rate. finally, we validated our context-specific metabolic models by incorporating codon usage and differential gene expression data into our model. we then compared trends in our models’ biomass predictions to literature on the bacterium. prediction of therapeutic targets essential gene analysis for each gene in the model, essential gene analysis involved removing reactions catalyzed by the gene or by a complex involving that gene and then using fba to predict growth. genes were considered essential if following their removal, the predicted maximum growth rate was zero. the c. difficile r, for which gene essentiality data was available for comparison with our in silico results, had been grown on tryptone-glucose-yeast extract broth. to approximate this medium, we used the complex medium defined by larocque et al. during essential gene analysis of imltc806cdf. pathway-oriented sensitivity analysis the growing research attention on metabolic pathways, rather than on specific reactions, is motivated by novel methods that allow for a better understanding of the functionality of complex webs of metabolic reactions. to date, much of the study of metabolic pathways, their crosstalks, and their role in the overall metabotype has been carried out with statistical and model-based approaches. sensitivity analysis is used to identify model inputs that have a large influence on the model outputs. to find the metabolic pathways that have the largest effect on the outputs of imltc806cdf and icdf we used pathway-oriented sensitivity analysis. posa involves genetically manipulating the metabolic model to find the sensitive pathways, which make a large impact on model outputs. in other words, we perturbed pathways by mutating the genes that govern their biochemical reactions and analyzed the result on the outputs. in the knockout vector y={b b …,b s,…,b p}, b s represents the perturbations on the genes governing the metabolic pathway s, where |b s|=w s. because the gene knockouts are represented through the use of binary variables, we perform combinatorial perturbations, namely the bits in b s are switched from to or from to 0; note that if a gene in b s is set to this gene is knocked-out in the model. according to, the pathway elementary effect for the genetic perturbation b s can be defined as follows: pees=∥f−f∥Δs, where bs~ represents the genetic manipulation of the input b s; y~ is the mutation carried out on the knockout vector y; f is the vector v of fluxes as produced by the model; finally, Δ s is a scale factor defined as: Δs=1ws∑i=1wsbs~,s= …,p. next, the sensitivity indices μ and σ are determined by calculating the mean and the standard deviation of the distribution of the pee for each input. pathways with a large μ have a large influence on the output. a large σ indicates an input whose influence highly depends on the value of other inputs. by perturbing the genes through the use of knockouts and comparing the outputs of the model with and without the genetic manipulations, we detected the most sensitive pathways of the metabolic models. calculation of flux control coefficients posa provides valuable information on sensitive pathways that can be targeted by therapies, but often more specific drug target predictions are desired. to understand how a metabolic pathway is controlled and can be altered, its control structure has to be determined. the flux control coefficient is the flux v ydh through a particular reaction, catalyzed by enzyme ydh, of the metabolic pathway with respect to the concentration x xase of an enzyme xase: cxxasevydh=∂vydh∂xxase·xxasevydh=∂lnvydh∂lnxxase in our calculations, the enzyme concentration was assumed to be equal to the gene expression level adjusted by cai. when calculating the flux control coefficients, we considered a 1% perturbation in the enzyme concentration x xase. flux control coefficients provide a quantitative measure of the degree of control an enzyme exerts on a metabolite flux and can quantitatively substitute for the qualitative concept of essential gene. thus, they can be used to identify steps that should be modified to achieve a successful alteration of the flux in outputs of clinical relevance. analysis of cdna microarrays we used microarray analysis to determine the combination of genes which were up-regulated or down-regulated in different environmental conditions. we used limma, a package in bioconductor, for statistical analysis of gene expression. we preprocessed the data through background correction, within-array normalization, and between-array normalization. after normalization, we used filtering to remove probes that did not appear to be expressed in any of the experimental conditions. next, we used linear models to analyze the microarray data. to conduct statistical analysis and assess differential expression, we used an empirical bayes method to modulate the standard errors of the log-fold changes. to test for the comparisons of interest, we used an analysis of variance model.  expansion and modification of imlt806cdf to icdf834 the genome of c. difficile strain is composed of a circular chromosome of, bp coding for open reading frames, along with a plasmid containing bp coding for orfs. the modified metabolic network draft contains 21% of the orfs present in the chromosomal genome of the bacteria with orfs, a modest improvement upon imltc806cdf, which contains orfs, as shown in table our expanded metabolic network also consists of metabolites and reactions. the final version of the network is available as an sbml file and as an excel file that indicate the reactions, metabolites, genes, and compartments involved in the metabolic network, along with references to literature that support additions or modifications to the existing network. the new network has two additional dead-end metabolites as compared with those found in imltc806cdf. the excel and sbml file, along with the the justification for keeping the dead-end metabolites in the model, have been uploaded to http://github.com/ssahebkashaf/peptoclostridiumdifficile the code for all of the analyses employed in our work is also freely available on this repository. genomic information of c. difficile reconstructed models imltc806cdf icdf834 we repeated analyses previously conducted by larocque et al. to validate imltc806cdf. namely, we compared the ability of icdf and imltc806cdf to identify essential amino acids and metabolizable carbon sources. the removal of amino acids that were not found to be essential or to affect growth, did not affect model-predicted biomass production in both models. moreover, no biomass was produced in the absence of essential amino acids in both models. therefore, similar to the previous network, our network is able to account for the essentiality of various amino acids on the growth of c.difficile. with regards to carbon sources, both models were able to correctly predict a range of carbon sources that are utilized by the bacterium. moreover, the bacterium was able to generate biomass in the absence of other carbon substrates, such as fructose, mannose, mannitol, and sorbitol. this finding is consistent with literature, which maintains that c. difficile is not restricted to metabolizing sugars and can ferment other compounds, even amino acids, to obtain both its carbon and energy. validation of metabolic models genetic analysis using multi-objective optimization our modeling approach is intended to simulate the conflicting objectives faced by the bacterium, where optimal performance in one objective coincides with sub-optimal performance in another objective. we used a knockout parameter space to find the genetic designs that would optimize the two objectives. in fig. we show the areas of objective space discovered by the genetic algorithm during the genetic analysis from the first generation to generation the optimization algorithm adaptively moves to regions that maximize biomass while minimizing the total intracellular flux, as evident in the curvature of the plot in fig. after conducting the genetic analysis, the pareto front, shown in black in the inset of fig. was determined. the pareto front is the set of nondominated solutions that represents the range of phenotypes resulting from different trade-offs between the two objectives. the presence of a pareto front, as opposed to a singular dominated solution, aligned with our a priori expectations regarding the metabolic plasticity inherent to the bacterium. our findings, along with previous literature on the choice of objectives, supported our choice of objectives to model c. difficile’s metabolism. fig. genetic analysis using multi-objective optimization. regions of objective space explored by the optimization algorithm for the objectives of maximization of biomass and minimization of total intracellular flux. solutions are represented by progressively warmer colors depending on the time step of the algorithm in which they had been adaptively generated from the initial point. the pareto front is shown in black in the inset robustness analysis we gauged the robustness of our model by determining the change in the maximal biomass flux in response to different perturbations. global robustness analysis revealed that the biomass production was fully robust to perturbations for a flux bound perturbation and a tolerance of 1%. the gr falls when σ is increased or ε is decreased. this facet of the bacterium’s metabolism was biologically relevant as bacteria such as c. difficile are able to grow despite small fluctuations in their physical environment. robustness analysis illustrated that the global behavior of our metabolic model matches our expectations from biological rationale and supported the use of our models to predict the behavior of the bacterium in different environments. changes in c. difficile’s growth in different conditions we obtained the relevant microarray datasets from the gene expression omnibus database under the accession numbers gse and from the arrayexpress database under the accession numbers e-geod- and e-bugs- context-specific models for c. difficile were generated by incorporating gene expression data obtained for the bacterium in different environmental conditions. to improve the reliability of the model, we also integrated codon usage data. model predictions of these context-specific models were compared to expectations about the organism’s behavior from literature. previous work suggests that sub-mic concentrations of amoxicillin, metronidazole, and clindamycin slowed growth of toxigenic c. difficile as compared with the controls. to test these findings in silico, we incorporated gene expression levels of c. difficile in response to sub-mic levels of different antibiotics into our model. as compared with the c. difficile grown on bhi broth, toxigenic strains of c. difficile grown on sub-inhibitory concentrations of antibiotics exhibited reductions in their biomass, with those grown on amoxicillin showing the smallest growth. this finding is supported by literature that has shown that in vitro, amoxicillin is effective against c. difficile. these findings have lead to speculations that in vivo, this antibiotic is effective aginst vegetative forms of the bacterium but not against c. difficile spores. another potential explanation is that this broad-spectrum antibiotic may impair the intestinal microflora in a way that supports proliferation of c. difficile. the microarray data for each condition was obtained from the geo or arrayexpress databases, using the specified accession numbers. the differential gene expression levels obtained from analysis of this microarray data was used to make a metabolic model for each condition. these context-specific metabolic models were used to predict change in biomass production for each condition compared with the control of each microarray dataset additionally, the decline in biomass production following heat shock from to °c shown in table could be due to the general stress response employed by the bacterium. the heat shock response of c. difficile has been found to be involve gene clusters homologous to e. coli heat-shock operons. the heat shock response in e. coli has been found to be associated with a decrease in central carbon metabolism and a decline in cellular growth. literature on related bacteria is thereby in agreement with the model’s prediction of a significant reduction in growth in c. difficile following the heat shock. additionally, according to the work of dubois et al., the supplementation of 10mm cysteine to the medium did not affect c. difficile’s growth. after integrating the microarray data from their work, we found that our in silico findings agreed with their experimental results. validation of the findings of our context-specific metabolic models against the literature on the bacterium showed that metabolic models allow for an enriched view of omic data and may be valuable tools for better understanding the behavior of c. difficile in different conditions. prediction of therapeutic targets gene essentiality analysis essential genes have been cited as promising targets for development of new antimicrobials due to their importance for bacterial survival. using fba, we performed an in silico gene deletion study to predict potential essential genes that may lead to the identification of new drug targets. this analysis had already been conducted for imltc806cdf based on a 5% threshold, and gene essentiality results had been compared to genes deemed essential for b. subtilis, for which this data had been available. we performed gene essentiality analysis for both imltc806cdf and icdf and validated our results using recently available literature on the essential genes of the c. difficile r. while imltc806cdf predicted essential genes and had a % accuracy in predicting gene essentiality, icdf predicted essential genes and had a % accuracy. pathway-oriented sensitivity analysis and flux control coefficients for our posa analysis, we chose the gene expression profile of the bacterium when grown on bhi broth. each pathway was assessed through random perturbations of its reactions, and the average perturbation μ and the standard deviation σ were computed as a result. we performed the pathway-based sensitivity analysis and identified sensitive pathways before and after modifying the metabolic model as shown in fig. the pathway with the largest μ, and thereby the greatest control on biomass production or growth in both imltc806cdf and icdf is the valine, leucine, and isoleucine metabolism pathway. these three amino acids are essential to the bacterium and their metabolism was also expected to be essential. the second most sensitive pathway is alanine, aspartate, and glutmate metabolism in imltc806cdf and glycolysis/gluconeogenesis in icdf additional sensitive pathways in icdf include pyrimidine metabolism and pyruvate metabolism. model findings suggest that therapies against infection may likely be more effective if they target key enzymes in these sensitive pathways. fig. posa was used to compare the most sensitive pathways of imltc806cdf and icdf the imltc806cdf model is composed of metabolic pathways and the icdf model is composed of metabolic pathways. biomass production is most sensitive to pathways with higher calculated μ to find more specific therapeutic targets, flux control coefficients for enzymes on biomass production in the metabolic model were determined and compared for bhi broth, cysteine supplementation, and heat shock gene expression data. the four enzymes with largest flux control coefficients in each condition are shown in fig. while the complete list of flux control coefficients in different conditions has been uploaded to the public repository. these flux control coefficients were interestingly involved in pathways deemed sensitive during posa. these enzymes varied amongst the four conditions, suggesting that access to the in vivo gene expression profile of c. difficile can be used to predict better drug targets for patients. therapies aimed at reducing growth of c. difficile should target enzymes with high flux coefficients as, according to our model, their activity is most closely tied to biomass production. fig. genes encoding the enzymes with the largest flux control coefficients for biomass production in different conditions. table of metabolic pathway hosting the genes and of gene descriptions. a flux control coefficient of implies full control of the metabolite flux by the associated enzyme  in this study, we expanded the existing metabolic network for c. difficile and used it to create context-specific metabolic models of its metabolism that allow us to understand how the bacterium alters its metabolism depending on its environment. to predict the bacterium’s behavior in different environmental conditions, the model was integrated with transcriptomic and codon usage data to generate reliable and context-specific metabolic flux distributions. we validated the model by conducting robustness and sensitivity analyses. we further assessed its predictive potential by comparing model predictions with published experimental data to gauge the consistency of model findings with the current knowledge of c. difficile’s metabolism. through this literature-based validation, we found that the model is a valuable tool for qualitatively understanding the behavior of the bacterium in different settings. the model can also be used to find potential therapeutic targets by allowing for determination of essential genes and context-specific sensitive pathways and flux control coefficients. context-specific metabolic models can allow for a better understanding different medically-relevant conditions and can be continuously refined by integrating novel information regarding c. difficile’s metabolism. our model can be used by biomedical researchers to study the bacterium and devise targeted treatments. our approach can also be scaled-up to simulate the interactions between the gut microbiota and the host using a bottom-up modeling approach. by accounting for the gut microbiota-host interactions, we can construct a whole gut model response to infections and other inflammatory events, paving the path towards more informed and effective treatments. availability and requirements project name: icdf834 project home page: http://github.com/ssahebkashaf/pep-toclostridiumdifficile630 operating system: platform independent programming language: matlab other requirements: gurobi license: university of cambridge any restrictions to use by non-academics: license needed additional files additional file carbon source usage, carboxylic acids secreted, and essential amino acids in icdf as compared with imltc806cdf and literature. additional file gr of biomass production in icdf to internal and external perturbations. additional file gene essentiality analysis of icdf and imltc806cdf. abbreviations caicodon adaptation index fbaflux balance analysis grglobal robustness orfopen reading frame peepathway elementary effect posapathway-oriented sensitivity analysis tgytryptone-glucose-yeast extract 
the standard metabolic model cannot account for changes in the bacterial metabolism in response to different environmental conditions. to our knowledge, this is the first time codon usage has been quantified and integrated into a metabolic model. to determine potential therapeutic targets using the model, we conducted gene essentiality and metabolic pathway sensitivity analyses and calculated flux control coefficients. after an extensive validation process, we used icdf to obtain state-of-the-art predictions of therapeutic targets for c. difficile. we show how context-specific metabolic models augmented with codon usage information can be a beneficial resource for better understanding c. difficile and for identifying novel therapeutic targets. importantly, to bridge the gap between gene expression levels and protein abundance, we accounted for the synonymous codon usage bias of the bacterium in the model. to account for this limitation, we also integrated transcriptomic data, which details the gene expression of the bacterium in a wide array of environments. we validated our context-specific metabolic models using sensitivity and robustness analyses and compared model predictions with literature on c. difficile. clostridium difficile is a bacterium which can infect various animal species, including humans. clostridium difficilemetabolic networksmetabolic pathwaysmetabolic modelinggenome scale modelingflux balance analysissensitivity analysisantibiotic resistanceissue-copyright-statement© the author 2017 we present icdf an updated metabolic network of c. difficile that builds on imltc806cdf and features reactions, genes, and metabolites. we obtained % accuracy in predicting gene essentiality when compared to experimental data for c. difficile r homologs. genome-scale metabolic models contain all known biochemical reactions of a microorganism and can be used to investigate this relationship. we used this metabolic network to reconstruct the metabolic landscape of this bacterium. the model predicts interesting facets of the bacterium’s metabolism, such as changes in the bacterium’s growth in response to different environmental conditions. the metabolic fluxes were defined as a function of protein abundance.
 clostridium difficile is a bacterium which can infect various animal species, including humans. infection with this bacterium is a leading healthcare-associated illness. a better understanding of this organism and the relationship between its genotype and phenotype is essential to the search for an effective treatment. genome-scale metabolic models contain all known biochemical reactions of a microorganism and can be used to investigate this relationship.  we present icdf an updated metabolic network of c. difficile that builds on imltc806cdf and features reactions, genes, and metabolites. we used this metabolic network to reconstruct the metabolic landscape of this bacterium. the standard metabolic model cannot account for changes in the bacterial metabolism in response to different environmental conditions. to account for this limitation, we also integrated transcriptomic data, which details the gene expression of the bacterium in a wide array of environments. importantly, to bridge the gap between gene expression levels and protein abundance, we accounted for the synonymous codon usage bias of the bacterium in the model. to our knowledge, this is the first time codon usage has been quantified and integrated into a metabolic model. the metabolic fluxes were defined as a function of protein abundance. to determine potential therapeutic targets using the model, we conducted gene essentiality and metabolic pathway sensitivity analyses and calculated flux control coefficients. we obtained % accuracy in predicting gene essentiality when compared to experimental data for c. difficile r homologs. we validated our context-specific metabolic models using sensitivity and robustness analyses and compared model predictions with literature on c. difficile. the model predicts interesting facets of the bacterium’s metabolism, such as changes in the bacterium’s growth in response to different environmental conditions.  after an extensive validation process, we used icdf to obtain state-of-the-art predictions of therapeutic targets for c. difficile. we show how context-specific metabolic models augmented with codon usage information can be a beneficial resource for better understanding c. difficile and for identifying novel therapeutic targets. we remark that our approach can be applied to investigate and treat against other pathogens. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords clostridium difficilemetabolic networksmetabolic pathwaysmetabolic modelinggenome scale modelingflux balance analysissensitivity analysisantibiotic resistanceissue-copyright-statement© the author 2017 