Metabolic Modeling of Cystic Fibrosis Airway Communities Predicts Mechanisms of Pathogen Dominance

Few taxonomic groups dominate the CF airway community samples.Principal-component analysis (PCA) was performed on the normalized read data of the 75 samples to evaluate sample heterogeneity. The first three principal components (PCs) captured 77.8% of the data variance, with the first PC capturing 57.3% of variance and most heavily weighting the most abundant genera Pseudomonas, Streptococcus, and Prevotella as expected (see Table S1 in the supplemental material). A considerable degree of heterogeneity was evident from a plot of the 75 samples in the coordinates defined by the first three PCs (Fig. 2A). Most striking were the outlier samples from three patients infected with Enterobacteriaceae (samples 25 to 27), Burkholderia (samples 19 to 21), or Achromobacter (samples 31 and 32) compared to the patients lacking these three organisms (i.e., the remaining 67 samples).

• Open in new tab
• Download powerpoint

FIG 2

PCA performed on the normalized read data. (A) PCA performed for all 75 samples with the normalized reads for each taxonomic group plotted using the first three principal components (PCs) that explained 57.3%, 12.3%, and 8.2%, respectively, of the data variance. Samples containing Enterobacteriaceae (samples 25 to 27), Burkholderia (samples 19 to 21), and Achromobacter (samples 31 and 32) appeared as outliers. (B) PCA performed for 67 samples when the 8 samples containing Enterobacteriaceae, Burkholderia, and Achromobacter were removed. The normalized reads for each taxonomic group were plotted using the first two PCs, which explained 72.6% and 11.7%, respectively, of the data variance.

Because each pathogen infected only a single patient among the 46 included patients, we generated a smaller data set of 67 samples by removing these 8 samples. When PCA was performed on this reduced data set, the first three PCs explained 92.6% of the data variance (Table S2), suggesting substantially reduced heterogeneity compared to the full data set. These three PCs heavily weighted only the four taxonomic groups Pseudomonas, Streptococcus, Prevotella, and Haemophilus, with the first PC representing high Pseudomonas and low Streptococcus, the second PC component representing high Streptococcus and moderate Pseudomonas, and the third PC representing high Haemophilus, low Pseudomonas, and low Streptococcus. Considerable heterogeneity was evident when the 67 samples were plotted using the first two PCs accounting for 84.2% of the variance (Fig. 2B). Here the first PC represented high Pseudomonas, low Streptococcus, moderate Prevotella, and moderate Haemophilus, and the second PC represented low Pseudomonas, high Streptococcus, low Prevotella, and low Haemophilus.

Based on these results, we focused our community modeling efforts on predicting the infrequent dominance of the pathogens Enterobacteriaceae, Burkholderia, and Achromobacter and the heterogeneity in the abundances of Pseudomonas, Streptococcus, Prevotella, and Haemophilus across the remaining samples. Pseudomonas, Streptococcus, and Prevotella have been found by directly sampling the lung of CF patients via bronchoalveolar lavage (40), while Haemophilus is a widely accepted CF pathogen (7). The other 10 genera (Table 1) were maintained in the model to simulate competition/cooperation with the more dominant species and to determine if the relatively low abundances of these genera could be predicted.

The community model can reproduce dominance of CF pathogens.We simulated the growth of each species individually to compare their monoculture growth rates with the nominal community nutrient uptake rates (Table S3). Interestingly, the three highest growth rates belonged to the rare pathogens Escherichia, Burkholderia, and Achromobacter, while the next three highest growth rates belonged to the common pathogens Pseudomonas, Streptococcus, and Staphylococcus (Fig. 3A; species numbered as in Table 1). These predictions were consistent with our modeling results for the gut microbiome (41), where opportunistic pathogens consistently had higher growth rates than commensal species. The other two species, Prevotella and Haemophilus, commonly observed in the 75 patient samples were predicted to have much lower in silico growth rates. The three species representing Fusobacterium, Granulicatella, and Porphyromonas did not grow individually due to their inability to meet the defined ATP maintenance demand, although they could grow when strategically combined with other modeled species. For example, Fusobacterium, Granulicatella, and Porphyromonas were predicted to grow in coculture with Ralstonia, Prevotella, and Actinomyces, respectively. The species abundances predicted for a specified nutrient condition depended on both the monoculture growth rates and the ability of each species to efficiently utilize secreted metabolites to enhance its growth rate. These emergent cross-feeding relationships allowed otherwise slower-growing species to coexist with species that exhibited high monoculture growth rates.

• Open in new tab
• Download powerpoint

FIG 3

Single-species and community simulations performed with the nominal nutrient uptake rates in Table S3. (A) Single-species growth rates with the species numbered according to Table 1. (B) Comparison of predicted species abundances to the average of the normalized reads for the single patient infected with Enterobacteriaceae/Escherichia (samples 25 to 27). (C) Comparison of predicted species abundances to the average of the normalized reads for the single patient infected with Burkholderia (samples 19 to 21). (D) Comparison of predicted species abundances to the average of the normalized reads for the single patient infected with Achromobacter (samples 31 and 32). (E) Comparison of predicted species abundances to the average of the normalized reads for the 43 patients not infected with Enterobacteriaceae/Escherichia, Burkholderia, or Achromobacter (samples 1 to 18, 22 to 24, 28 to 30, and 33 to 75).

We conducted simulations using the nominal nutrient uptake rates (Table S3) to determine if the community model could capture dominance of each rare pathogen in the absence of the other two rare pathogens. Each simulation was performed by constraining the abundances of the other two pathogens to zero, effectively producing reduced communities of 15 species. The predicted abundances from each simulation were compared to the normalized reads averaged over the patient samples which contained the associated pathogen: Enterobacteriaceae/Escherichia (samples 25 to 27) (Fig. 3B), Burkholderia (samples 19 to 21) (Fig. 3C), or Achromobacter (samples 31 and 32) (Fig. 3D). For each simulated case, the model correctly predicted dominance of the associated pathogen. For the Burkholderia- and Achromobacter-infected patients, the abundances of the dominant pathogen as well as less prevalent species were well predicted.

We performed simulations for the remaining 43 patients by reducing the community to 14 species by constraining the abundances of all three rare pathogens to zero. The model-predicted abundances were compared to the normalized reads averaged over the 67 samples remaining when the 8 rare pathogen-containing samples were removed (Fig. 3E). The model correctly predicted that Pseudomonas, Streptococcus, and Prevotella would dominate the community, although the Prevotella abundance was overpredicted at the expense of Streptococcus as well as several less abundant genera. The only other genus present in the simulated community was Staphylococcus, while the averaged reads showed a greater amount of diversity. Compared to the averaged data, individual samples showed less diversity, which is more consistent with model predictions as discussed below.

The community model can reproduce pathogen heterogeneity across airway samples.The CF airway communities exhibited a substantial degree of sample-to-sample heterogeneity when rare pathogens were present (Fig. 2A) or absent (Fig. 2B). We performed simulations to assess the extent to which sample-to-sample differences in taxonomic group reads could be explained by heterogeneity in the metabolic environment of the CF lung. More specifically, we randomized the community nutrient uptake rates around their nominal values (Materials and Methods; also see Table S3) to mimic heterogeneous lung environments shown to occur across CF patients (42, 43) and in longitudinal samples from a single patient (44). An objective of our future research will be to model sample-by-sample variability in individual patients as a function of disease state (e.g., clinically stable, pulmonary exacerbation, and antibiotic treatment). In this study, each simulation with a set of randomized uptake rates was termed a “simulated sample,” and we tested the hypothesis that the experimental samples could be interpreted as having been drawn from the much larger set of simulated samples we generated. Due to the relatively small number of Enterobacteriaceae/Escherichia-, Burkholderia-, and Achromobacter-containing samples, we performed only 100 randomized community simulations for each of these pathogens. In contrast, 1,000 randomized simulations were performed for communities without these three rare pathogens since the associated patient sample size was comparatively large. The single model simulation that best represented a particular patient sample was determined by the minimum least-squares error between the normalized measured reads and the predicted abundances across all simulations. For the 8 rare pathogen-containing samples, we plotted the measured reads and predicted abundances of the best-fit models for the five most common genera (Pseudomonas, Streptococcus, Prevotella, Haemophilus, and Staphylococcus) and the pathogen of interest (Fig. 4; Table S4). For the remaining 67 samples, we plotted the measured reads and predicted abundances of the best-fit models for the five most common genera plus the next most abundant genus according to measured reads (Fig. 5; Table S5).

• Open in new tab
• Download powerpoint

FIG 4

Taxonomic reads for patient samples containing rare pathogens compared to species abundances predicted from community models with randomized nutrient uptake rates. The genera Pseudomonas, Streptococcus, Prevotella, Haemophilus, and Staphylococcus and the indicated rare pathogen (Enterobacteriaceae/Escherichia, Burkholderia, or Achromobacter) are shown for each case. (A) Individual models that best fit the 3 Enterobacteriaceae/Escherichia-containing samples, 25 to 27, selected from an ensemble of 100 15-species models without Burkholderia or Achromobacter. (B) Individual models that best fit the 3 Burkholderia-containing samples, 19 to 21, selected from an ensemble of 100 15-species models without Enterobacteriaceae/Escherichia or Achromobacter. (C) Individual models that best fit the 2 Achromobacter-containing samples, 31 and 32, selected from an ensemble of 100 15-species models without Enterobacteriaceae/Escherichia or Burkholderia. Each abundance for a patient sample is shown in the first, dark-colored bar, and each abundance predicted by the corresponding model is shown in the second, light-colored bar.

• Open in new tab
• Download powerpoint

FIG 5

Taxonomic reads for patient samples without rare pathogens compared to species abundances predicted from community models with randomized nutrient uptake rates. The genera Pseudomonas, Streptococcus, Prevotella, Haemophilus, and Staphylococcus and the next most abundant genera are shown for each case. Individual models that best fit the 67 patient samples were selected from an ensemble of 1,000 14-species models without Enterobacteriaceae/Escherichia, Burkholderia, or Achromobacter. (A) Three representative samples for which the least-squares error measures were within the smallest third of all samples. (B) Three representative samples for which the least-squares error measures were within the middle third of all samples. (C) Three representative samples for which the least-squares error measures were within the largest third of all samples. Each abundance for a patient sample is shown in the first, dark-colored bar, and each abundance predicted by the corresponding model is shown in the second, light-colored bar.

Randomized nutrient simulations were able to generate model predictions that reproduced the major features of the 3 Enterobacteriaceae/Escherichia-containing samples (Fig. 4A), including the high-Enterobacteriaceae/Escherichia reads and the presence of the other main community members (Pseudomonas, Streptococcus, and Prevotella). The Streptococcus reads were predicted relatively accurately, while Pseudomonas reads were underpredicted and Prevotella reads were overpredicted. As measured by the least-squares error, improved predictions were obtained for the 3 Burkholderia-containing samples (Fig. 4B). The Burkholderia reads were accurately reproduced, and Streptococcus was correctly predicted to be the second most abundant genus, suggesting a synergism between these two genera. This prediction has experimental support from in vitro experiments showing that mucin-degrading anaerobes such as streptococci promote the growth of CF pathogens such as Burkholderia cenocepacia when mucins are provided as the sole carbon source (45). The two Achromobacter-containing samples were well predicted in terms of Achromobacter reads and Pseudomonas being the other dominant genus (Fig. 4C). These predictions are consistent with an in vitro study showing that Achromobacter sp. enhanced the ability of multiple P. aeruginosa strains to form biofilms (46). Furthermore, a clinical study with 53 patients having positive cultures for A. xylosoxidans showed that all 6 patients who were chronically infected by A. xylosoxidans were coinfected with P. aeruginosa (47). Complete comparisons of the normalized measured reads and model predicted abundances for the 8 samples with the rare pathogens are presented in Table S4, which shows that the model generally produced less diverse communities as measured by the richness (number of species with abundances exceeding 1%) and the equitability (the inverse Simpson metric [48]).

The lack of patient samples containing Enterobacteriaceae/Escherichia, Burkholderia, and Achromobacter limited our ability to analyze heterogeneity of communities with these pathogens. In contrast, the 67 samples remaining when the 8 samples containing these three pathogens were removed offered a much larger data set for heterogeneity analysis. Each of these 67 samples was matched to one of the 1,000 randomized model simulations according to the smallest least-squares error between the normalized reads of the sample and the predicted abundances of the model (Table S5). Representative results are shown for patient samples with relatively small (0.004 to 0.087, Fig. 5A), moderate (0.089 to 0.116, Fig. 5B), and large (0.123 to 0.574, Fig. 5C) error values. Samples which were most accurately reproduced generally contained high Pseudomonas reads (84% ± 15%) with the remainder of the community consisting of Streptococcus and Prevotella (Fig. 5A). These 22 samples were best matched by 11 distinct models, suggesting that patient samples dominated by Pseudomonas contained a higher degree of heterogeneity than the simulated samples.

The 22 samples which produced moderate prediction errors were characterized by lower and more variable Pseudomonas reads (48% ± 28%) as well as more variable distributions of Streptococcus and Prevotella reads (Fig. 5B). The ensemble of randomized models could capture the relative amounts of these three genera but often predicted the presence of Staphylococcus not observed in the patient samples. This discrepancy could be attributable to the unmodeled ability of Pseudomonas to secrete diffusible toxins which inhibit Staphylococcus respiration and render Staphylococcus less metabolically competitive in partially aerobic environments (49) such as the CF lung. Interestingly, the model ensemble could reproduce the relatively high Ralstonia reads in sample 1 while also predicting no Ralstonia in samples 15 and 69. The 23 samples which produced the largest prediction errors were characterized by much lower Pseudomonas reads (13%), higher reads of Streptococcus and Prevotella (34% and 19%, respectively; e.g., samples 26 and 74 in Fig. 5C), and higher representation of less common genera. These samples also produced higher Haemophilus reads, primarily due to two Haemophilus-dominated samples (e.g., sample 39 in Fig. 5C). While the model ensemble generally was able to reproduce the observed Streptococcus and Prevotella reads in these samples, the models tended to overpredict Pseudomonas and Staphylococcus at the expense of the less common genera. In particular, the ensemble underpredicted the abundances of Rothia, Fusobacterium, and Gemella while the average reads of these three genera across the 23 samples summed to 16% This discrepancy could suggest that these 23 samples were obtained from patients with less advanced CF lung disease, which correlates to higher diversity communities in vivo (30, 50).

To gain further insights into the ability of the community model to mimic sample-to-sample heterogeneity in the absence of rare pathogens, we compared read data and abundance predictions in the PC space calculated from the 67 patient samples. Each of the 1,000 model simulations was mapped into the two-dimensional space defined by the first two PCs (Fig. 2B), which explained 84.2% of normalized read data variance (Table S2). The model ensemble was able to reproduce most of the observed variability as reflected by the cloud of model simulations overlapping 56 of the 67 patient samples (Fig. 6A). The patient and simulated samples covered the same range of the first PC, which was heavily weighted by Pseudomonas, Streptococcus, and Prevotella (Table S2). Importantly, this consistency shows that heterogeneity across these three dominant genera could be predicted from variations in the CF lung metabolic environment, as we hypothesized.

• Open in new tab
• Download powerpoint

FIG 6

Principal-component analysis (PCA) of taxonomic reads for patient samples without rare pathogens and species abundances predicted from 14-species community models with randomized nutrient uptake rates. (A) Representation of the 67 patient samples (blue crosses labeled with sample number) in the two-dimensional space defined by the first two principal components (PCs) obtained when PCA is performed on the normalized reads of these patient samples. Predicted species abundances (red circles) from an ensemble of 1,000 models transformed into the PC space of the normalized read data. (B) Enlarged view of the lower left portion of the PCA plot in panel A. (C) Average genus reads obtained for 12 samples (samples 5, 6, 10, 39, 42, 43, 49, 57, 61, 68, 70, and 74) in panel B with elevated Prevotella representation compared to the average abundances predicted from the best-fit models for these samples with the species number as in Table 1.

The model ensemble also could reproduce variations in the second PC, which was heavily weighted by the three dominant genera and Haemophilus, for sufficiently large values of the first PC, which corresponded to relatively high Pseudomonas and low Streptococcus and Prevotella. In contrast, the model ensemble did not cover the patient samples in the lower left quadrant of the PC plot (Fig. 6B). These samples were characterized by unusual combinations of relatively high Prevotella, Haemophilus, Rothia, and/or Fusobacterium that the model could not reproduce in its present form. Of these 12 poorly modeled samples, Prevotella was highly represented in 8 samples. When the normalized reads of these 8 samples and their associated best-fit abundances were averaged, the models overpredicted Pseudomonas, Streptococcus, and Staphylococcus at the expense of the less common genera (Fig. 6C).

The community model predicts that pathogen dominance is driven by metabolite cross-feeding.To investigate putative metabolic mechanisms by which pathogens may establish dominance in the CF lung, we used model predictions to quantify rates of metabolite cross-feeding between species. For each rare pathogen (Escherichia, Burkholderia, and Achromobacter), 100 simulations performed with randomized community uptake rates were used to calculate average exchange rates of the five most significantly cross-fed metabolites between Pseudomonas, Streptococcus, and the pathogen of interest. The overall metabolite exchange rate from one species to another species was calculated by determining the minimum uptake or secretion rate for each exchanged metabolite and then summing these minimum rates over all exchanged metabolites.

Escherichia was predicted to consume the organic acids acetate, formate, and l-lactate produced by Streptococcus, while Streptococcus benefitted from the amino acids serine and threonine secreted by Escherichia (Fig. 7A and D). Due to the existence of alternative optima with respect to the secretion products (51), l-lactate secretion was not predicted in Streptococcus monoculture even through the metabolic reconstruction supported l-lactate production (21) (www.vmh.life). While Streptococcus strains are well known to product l-lactate as the primary product via homolactic fermentation (52, 53), we chose not to manually curate the metabolic reconstruction since in silico l-lactate synthesis was induced by the presence of other community members such as Escherichia. Pseudomonas was minimally involved in metabolite exchange due to its low average abundance (∼1%) across the 100 simulations. Hence, our model suggested that organic acid cross-feeding could play a role in Enterobacteriaceae propagation in the CF lung.

• Open in new tab
• Download powerpoint

FIG 7

Predicted metabolite cross-feeding relationships for 15-species communities containing Escherichia, Burkholderia, or Achromobacter. Negative rates denote metabolite uptake, and positive rates denote metabolite secretion. The overall metabolite exchange rate from one species to another species was calculated by determining the minimum uptake or secretion rate for each exchanged metabolite and then summing these minimum rates over all exchanged metabolites. The arrow thickness is proportional to the overall metabolite exchange rate between the two species. (A) Average exchange rates of the five highest cross-fed metabolites between the three most abundant species for 100 model ensemble simulations containing Escherichia. (B) Average exchange rates of the five highest cross-fed metabolites between the three most abundant species for 100 model ensemble simulations containing Burkholderia. (C) Average exchange rates of the five highest cross-fed metabolites between the three most abundant species for 100 model ensemble simulations containing Achromobacter. (D) Schematic representation of overall metabolite exchange rates for Escherichia-containing communities corresponding to panel A. Pseudomonas was omitted due to its low exchange rates compared to the other two species. (E) Schematic representation of overall metabolite exchange rates for Burkholderia-containing communities corresponding to panel B. (F) Schematic representation of overall metabolite exchange rates for Achromobacter-containing communities corresponding to panel C.

More complex cross-feeding relationships were predicted for Burkholderia-containing communities that supported average Pseudomonas and Streptococcus abundances both exceeding 10%. The highest exchange rates were predicted for formate and acetate produced by Streptococcus and consumed by Burkholderia (Fig. 7B and E). The two species also exchanged amino acids, with Streptococcus providing alanine to Burkholderia and Burkholderia producing aspartate and serine for Streptococcus. Burkholderia provided the same two amino acids to Pseudomonas while receiving a small exchange of acetate in return. Pseudomonas also consumed formate secreted by Streptococcus. These model predictions suggested that acetate, formate, and alanine produced by Streptococcus via heterolactic fermentation (52) could promote Burkholderia growth in vivo. Indeed, in vitro experiments have shown that mucin-degrading anaerobes such as streptococci may promote the growth of CF pathogens such as B. cenocepacia by secreting acetate (45).

Compared to the other two pathogens, Achromobacter was predicted to be less efficient at cross-feeding, having only low uptake rates of alanine, l-lactate, and threonine secreted by the other two species. In contrast, Pseudomonas was predicted to benefit from relatively high uptake rates of formate produced by Streptococcus and succinate produced by Achromobacter. Collectively, these model predictions could help explain the enhanced ability of Burkholderia to dominate the simulated CF airway communities compared to Achromobacter (Fig. 4) despite the single-species growth rates of the two species being similar (Fig. 3).

Similar cross-feeding analyses were performed for 1,000 simulations with randomized nutrient uptake rates in 14-species communities lacking Escherichia, Burkholderia, and Achromobacter. To investigate the possibility of differential cross-feeding patterns, the simulations were split into 500 cases with the highest Pseudomonas abundances and 500 cases with the lowest Pseudomonas abundances (Fig. 8A). For each set of 500 simulations, the average exchange rates of the five most significantly cross-fed metabolites between the four most abundant species (Pseudomonas, Streptococcus, Prevotella, and Staphylococcus) were calculated. The overall metabolite exchange rate between any two species were calculated from the individual metabolite uptake and secretion rates as before.

• Open in new tab
• Download powerpoint

FIG 8

Predicted metabolite cross-feeding relationships for 14-species communities without Escherichia, Burkholderia, and Achromobacter. One thousand model ensemble simulations were performed and split into 500 cases with relatively high Pseudomonas abundances and 500 cases with relatively low Pseudomonas abundances. (A) Average abundances of the five most highly represented species for the high- and low-Pseudomonas-abundance cases. (B) Average exchange rates of the five highest cross-fed metabolites between the four most abundant species for high-Pseudomonas-abundance cases. (C) Average exchange rates of the five highest cross-fed metabolites between the four most abundant species for low-Pseudomonas-abundance cases. (D) Schematic representation of overall metabolite exchange rates for high-Pseudomonas-abundance cases corresponding to panel B. (E) Schematic representation of overall metabolite exchange rates for low-Pseudomonas-abundance cases corresponding to panel C.

When Pseudomonas abundances were predicted to be relatively high (average of 61%), community interactions were dominated by Pseudomonas consumption of formate, ethanol, acetate, and aspartate secreted by the other three species (Fig. 8B). Formate cross-feeding was predicted to be particularly important, which was consistent with an in vitro study showing that expression of the P. aeruginosa fdnH gene (encoding a formate dehydrogenase) was elevated in synthetic sputum medium compared to glucose minimal medium (54). Similarly, the expression of P. aeruginosa adhA (encoding an alcohol dehydrogenase) was elevated in patient-derived CF sputum compared to in vitro rich medium (55). Since P. aeruginosa strains have the capability to take up both formate and ethanol (56, 57), these in vitro studies suggest that this cross-feeding mechanism could occur in CF airway communities. Staphylococcus was the major source of exchanged formate and ethanol (Fig. 8D), a prediction consistent with studies showing that P. aeruginosa benefits from the presence of S. aureus (49, 58). Both alanine and aspartate have been shown to serve as preferred carbon sources for P. aeruginosa in a minimal medium supplemented with lyophilized CF sputum (54). However, the ensemble model did not predict exchange of l-lactate between P. aeruginosa and S. aureus, which differs from coculture experiments that mimic the CF lung environment (49). Strong interactions between P. aeruginosa and various streptococci also have been reported (30), although the importance of metabolite cross-feeding in mediating these interactions remains incompletely understood (59). Finally, in the model Pseudomonas supplied small amounts of d-lactate for Prevotella and Staphylococcus consumption, a prediction consistent with an in vitro study showing P. aeruginosa anaerobic production of the LldA enzyme catalyzing d-lactate synthesis (60).

When Pseudomonas abundances were predicted to be relatively low (average of 32%), metabolite cross-feeding remained dominated by Pseudomonas consumption of secreted by-products and amino acids (Fig. 8C). Pseudomonas was predicted to have high consumption rates of formate produced by all three other species and l-lactate synthesized only by Streptococcus, consistent with the ability of Streptococcus salivarius (61) and P. aeruginosa (49) to synthesize and consume l-lactate, respectively. Higher exchange rates between Streptococcus and Staphylococcus were predicted when Pseudomonas abundances were relatively low (Fig. 8E). The two species cross-fed alanine and l-lactate produced by Streptococcus and aspartate and ethanol secreted by Staphylococcus. Our predicted cross-feeding relationships in Pseudomonas- and Streptococcus-dominated communities could provide insights into CF disease progression, as high abundances of Streptococcus relative to Pseudomonas have been shown to correlate with higher-diversity airway communities and improved CF clinical stability (30). Younger CF patients also are known to have more diverse airway communities (62), so such interpretations would need to be made with care.