Plasticity of the Systemic Inflammatory Response to Acute Infection during Critical Illness: Development of the Riboleukogram

Background Diagnosis of acute infection in the critically ill remains a challenge. We hypothesized that circulating leukocyte transcriptional profiles can be used to monitor the host response to and recovery from infection complicating critical illness. Methodology/Principal Findings A translational research approach was employed. Fifteen mice underwent intratracheal injections of live P. aeruginosa, P. aeruginosa endotoxin, live S. pneumoniae, or normal saline. At 24 hours after injury, GeneChip microarray analysis of circulating buffy coat RNA identified 219 genes that distinguished between the pulmonary insults and differences in 7-day mortality. Similarly, buffy coat microarray expression profiles were generated from 27 mechanically ventilated patients every two days for up to three weeks. Significant heterogeneity of VAP microarray profiles was observed secondary to patient ethnicity, age, and gender, yet 85 genes were identified with consistent changes in abundance during the seven days bracketing the diagnosis of VAP. Principal components analysis of these 85 genes appeared to differentiate between the responses of subjects who did versus those who did not develop VAP, as defined by a general trajectory (riboleukogram) for the onset and resolution of VAP. As patients recovered from critical illness complicated by acute infection, the riboleukograms converged, consistent with an immune attractor. Conclusions/Significance Here we present the culmination of a mouse pneumonia study, demonstrating for the first time that disease trajectories derived from microarray expression profiles can be used to quantitatively track the clinical course of acute disease and identify a state of immune recovery. These data suggest that the onset of an infection-specific transcriptional program may precede the clinical diagnosis of pneumonia in patients. Moreover, riboleukograms may help explain variance in the host response due to differences in ethnic background, gender, and pathogen. Prospective clinical trials are indicated to validate our results and test the clinical utility of riboleukograms.


Introduction
Critical illness is marked by organ dysfunction, the need for vital support, and a high risk of death, occurring against a backdrop of systemic immune activation. This immune activation may begin as an adaptive response to the initial injury, however, as the disease progresses, the immune response may become maladaptive or paralyzed [1,2]. Critical illness-associated immune dysregulation has been described as the interplay between pro-and antiinflammatory responses [3], although recent evidence suggests a mixed inflammatory state is common [4,5]. While this process has been qualitatively described, there are no quantitative diagnostic or prognostic tools that have been validated clinically to assess immune status in the critically ill [6]. Consequently, infectious complications are not only common in intensive care units but also difficult to diagnose [7]. This has contributed to inappropriate use of broad-spectrum antibiotics and the emergence of multi-drug resistant organisms [8,9].
A few years ago, studies employing cultured human cells suggested that instead of a single molecule (e.g., IL-6), a constellation of molecules could be used to monitor the complexities of the inflammatory response, serving as markers of infection [10,11]. Miniaturized, multiplexed assays provide a rapid method for the unbiased screening of thousands of molecular species in a single assay [12]. These technological advances provided the potential for investigators to leverage high-throughput assays to better study the host response to and recovery from critical illness and injury [13,14]. Improved molecular diagnostics and prognostics, a better understanding of the complexity of systemic inflammation, and new therapeutic targets are expected deliverables, as reviewed recently [4,12].
Based upon our ability to diagnose abdominal sepsis in pilot mouse studies [15], we hypothesized that the host response to infection could not only differentiate between infected and noninfected states, but could also be used clinically to differentiate between the host response to infectious agents and to model the host response to and recovery from infectious perturbations. Pneumonia was chosen as an immune system perturbation, given its relative frequency and considerable cost in terms of patient morbidity and health care expense [7,16]. A bench-to-bedside, translational approach was employed to study the host response to pneumonia in critically ill subjects, comparing the informational content of standard clinical parameters and plasma cytokines to changes in the RNA abundance in circulating leukocytes.

Mice, experimental procedures, and samples
Care and use of laboratory animals were conducted in accordance with a protocol approved by the Washington University Animal Studies Committee, in compliance with guidelines (N01-RR-2-2135) prepared by the Committee on Care and Use of Laboratory Animals, Division of Research Resources, National Institutes of Health. Seven to nine week-old, male C57BL/6 mice were purchased (Harlan, Inc. Madison WI) and allowed to acclimatize for at least one week in a temperature-and light-controlled, pathogen-free barrier facility. Treated animals and concurrently studied controls were observed at 24 hour intervals for survival over eight days. In additional cohorts, whole blood was collected at 24 hours after injury.
In three additional cohorts of animals, blood was collected into an EDTA-coated syringe from the inferior vena cava being careful to avoid contamination of the needle with other tissues. Blood was diluted 1:1 with normal saline, pooled for the 8 animals in each treatment group, and separated into cells and plasma. Plasma was stored at 280uC until use. Erythrocytes were lysed hypotonically and RNA from peripheral leukocytes was harvested using RLT (Qiagen) and stored at 280uC until use. The 24 hour time point after injury was chosen as a time before appreciable mortality develops in animals with significant lung injury.

Target cRNA and gene expression signal
Each RNA sample was run on one GeneChip (a total of 15 mouse blood GeneChips from 120 animals). Total RNA from mouse blood was extracted as previously described [19]. cRNA target for GeneChip hybridization was prepared from total RNA (Affymetrix, Santa Clara, CA). Both total RNA and cRNA were electrophoretically assessed for quality (Agilent Bioanalyzer). The mouse blood cRNA samples were hybridized with the U74Av2 GeneChip (approximately 12,400 probe sets). Fluorescent hybridization signal was detected using a GeneChip Scanner 3000 (Affymetrix). These mouse microarray data (and those for patients, see below) have been deposited in NCBI's Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) and are accessible through GEO Series accession number GSE6377.

Data analysis and statistical tests for differential expression
Expression values were calculated from GeneChip .cel files using Robust Multichip Average (RMA) software [20]. Differentially expressed genes were identified using a mixed-model analysis of variance (ANOVA) and linear contrasts (PartekH Infer TM software) as previously reported [15]. Leave-one out crossvalidation (k-nearest neighbors, k = 2) was used to determine the reproducibility within this experimental set. Principal components analysis (PCA) was used to visually explore global effects for genome-wide trends, unexpected effects, and outliers in the expression data (PartekH Pro TM software, www.partek.com).
Patient studies. After obtaining informed consent, venous blood (7 ml) was collected from mechanically ventilated, nonseptic patients according to a protocol approved by the Washington University Institutional Review Board (#2004-0294). Patients were candidates for enrollment if they were on mechanical ventilation in the surgical ICU medical ICU, neurological ICU, or cardiothoracic ICU (CTICU) for 48 hours, were expected by the attending ICU physician to need mechanical ventilation for at least another 48 hours at the time of enrollment, and could provide written informed consent (from the patient or legal surrogate). VAP was diagnosed by the ICU attending physician, consistent with recently reported recommendations [16], without input from the investigators. Clinical data were entered into a VAP database, including gender, ethnic background, age, admitting diagnosis, type of ICU, APACHE II score, airway sampling technique and culture results, initial antibiotic therapy, and maximal clinical pulmonary infection score (CPIS) calculated based upon available data (several patients lacked daily arterial blood gas measurements to calculate P a O 2 / F i O 2 ratios) [21]. Patient blood was processed as described previously to minimize red blood cell RNA artifact [13]; briefly, samples were centrifuged (4006g 109 RT) to form a buffy coat and to separate plasma from cells. Plasma was withdrawn and stored at 280uC until use. Blood cells were diluted into EL buffer (90 ml) (Qiagen) and incubated on ice for 159. Leukocytes were pelleted by centrifugation (4006g 109 4uC), washed with EL buffer (30 ml) and lysed into RLT buffer (Qiagen) containing 1% bmercaptoethanol. Genomic DNA was sheared by repeated passage through an 18 gauge needle and the resultant material was stored at 280uC.
Patient blood leukocyte mRNA profiling. RNA was extracted, amplified and assessed for quality as described for murine samples. cRNA was hybridized against the HG-FOCUS array (Affymetrix, .8700 probe sets encoding ,8400 genes) and imaged as described for murine samples. Orthologs of murine genes were identified by comparison of the GeneChip Identifiers using the NetAffx Toolkit (Affymetrix). Consistent with recently published consensus statements [16,21], clinical data were judged to determine when (if ever) each patient developed ventilatorassociated pneumonia, with each patient acting as her/his own control. The timeline for each patient was defined such that the day of VAP diagnosis by the ICU attending physician was defined as day 0. A seven day time window from the gene expression time series was chosen as days 23, 21, 0, 1, and 3, with day 0 being the day that a patient was diagnosed as having VAP by the attending physician. Because blood samples were collected every other day, patients' samples were collected either on days 22, 0, +2 or on days 23, 21, +1, +3 relative to the VAP day of diagnosis described as time 0. For the purpose of analyzing the data from those patients who had samples collected on ''odd'' days, the time 0 data for these patients were interpolated. Those mRNA species whose abundance changed concordantly among the patients during the 7-day window surrounding the date of diagnosis were identified using extraction of Differential Gene Expression (EDGE) software [23]. Online databases were used to determine gene annotation and functional categorization (DAVID 2.0 accessed 16 November 2006) [24].
Clustering algorithm. Patient genes identified by EDGE were clustered as described previously [25]. Briefly, a gene coexpression network was constructed by connecting every gene to the top d genes (d = 5 in this study) to which its expression profile is most similar. The network then was partitioned into a set of communities, i.e., relatively densely connected sub-graphs, by a spectral graph algorithm [26]. The genes within each community formed a cluster. The number of clusters was determined automatically by the algorithm in order to maximize a modularity score [25,26]. Gene expression data were normalized prior to clustering such that the expression levels of each gene for each patient had a mean of zero and a standard deviation of one. Similarity between two gene expression profiles was measured by Pearson correlation coefficient.
Karhunen-Loeve decomposition of microarray gene expression data. To determine the dynamics of the host response to pneumonia, we constructed first a raw gene expression matrix corresponding to the ith pneumonia patient after RMA normalization [20] to be X i = [x i (1) … x i (N T )] where x i (k)MR N , (k = 1,…,N T , i = 1,…11, and N = 8793 genes) is the gene expression vector at k th time point. k = 1,…,N T where N T = 11 for most patients, correspond to sample collecting days 1, 3, 5, …, 21. Note that for few patients, only a portion of the time series (i.e., less than 11 points) was available.
For the EDGE-selected genes, the data were projected onto a smaller dimensional space using the series expansion method similar to principal components analysis, Karhunen-Loeve Decomposition (KLD) [27]. In order to preserve the alignment of the time series with respect to day 0 of VAP, we first obtained the average expression vectors x k ð Þ, k = 1,2,…9, by averaging the expression values at days corresponding to 23,21,…,13 in all patients (corresponding to the nine time points that most patients had samples collected). The KLD method looks for a basis y 1 ,y 2 ,…y N so that one can expand x k ð Þ as where a i k ð Þ~Sx k ð Þ,y i T and AE?,?ae stands for the standard inner product.
The orthonormal basis y 1 ,y 2 ,…y N can be selected as the eigenvectors of the correlation matrix C 1 MR N6N , obtained as The first principal mode y 1 corresponds to a constant bias term.
Hence the most important variation is captured by y 2 ,y 3 ,…y N respectively in the decreasing order. Once the orthonormal basis is obtained, each patient data can be projected onto this basis as a j i k ð Þ~Sx j k ð Þ,y i T for patients j = 1, 2, … . The discrete derivative of the coefficients a j i k ð Þ at time k is obtained as Validation of microarray results. Select genes were subjected to real-time quantitative PCR (RTq-PCR) for independent confirmation of relative expression levels. cDNA was generated and 100 ng was subjected to routine SybrGreen RT-PCR as per manufacturer's instruction (Applied Biosystems, Foster City, CA). In addition, 85 genes were selected at random from the total number on the GeneChip. The gene expression signal from these genes was analyzed in a manner identical to that described above. This procedure was repeated 100 times to estimate the informational value of randomly selected genes.

Validation of the riboleukogram
After enrollment of the first 20 patients, a second cohort of 7 patients was analyzed to validate the informational content of the leukocyte RNA species (genes) that changed in abundance in response to critical illness complicated by VAP. The blood handling, processing, and GeneChip analysis protocols were identical to those described above.

Murine model
Intratracheal (i.t.) installation of saline caused no deaths over 7days whereas i.t. introduction of P. aeruginosa endotoxin resulted in death of 90% of the animals studied within 96 hours ( Figure 1A). The 7-day mortality caused by live P. aeruginosa was adjusted by varying the size of the inoculum and injuries causing 50% and 90% 7-day mortality were achieved. A dose of S. pneumoniae was given i.t. that resulted in 85% 7-day mortality. Once these injuries were established, three separate cohorts of mice were used for each experimental group in the subsequent studies.
Peripheral blood was collected and pooled from groups of 4-7 mice 24 h after surgery, prior to appreciable mortality in any group ( Figure 1A). All murine RNA were of good quality based on the peak profiles of 18S and 28S ribosomal RNA. cRNA generated from these samples had a uniform size distribution. All hybridizations were of good quality; both the number of features present (35-40%) and the signals on each array fell within acceptable ranges [28].
Analysis of normalized gene expression data identified 219 probe sets (40 unannotated ESTs, 10 redundant probe sets, 169 annotated genes (Table S1) whose expression levels differentiated between the five groups. Leave one out cross-validation using knearest neighbors (k = 2) resulted in a 93% classification accuracy. The single misclassified sample was from the low-dose P. aeruginosa infection and was classified as ''saline''. These 219 genes differentiated between the host responses to Gram-negative bacteria, Gram-negative bacterial toxin (LPS) and Gram-positive bacteria. The probe sets clustered into six groups ( Figure 1B) and these groups defined the gene expression cartography of the murine response to pneumonia ( Figure 1C). Genes that fell within clusters 2 and 3 exhibited increased RNA abundance in animals responding to high lethality insults. Genes that fell within clusters 4-6 were transcriptionally suppressed during high lethality insults. Cluster 1 bridged the two groups. Gene ontology assignments identified enriched molecular functions in the three distinct groups. The bridging cluster (cluster 1) was enriched for genes involved in the immune response (N = 10, P = 4610 27 ) and genes with NTPase activity (N = 6, P = 5610 24 ). Clusters 2 and 3 were enriched for genes involved in intracellular signaling pathways Figure 1. (A) Eight-day survival curves of mice challenged with intra-tracheal injection of one of 5 solutions, each dosed to produce the observed mortality. Significant differences (p,0.05) were observed between the 80-90% mortality (Pseudomonas bacteria, Streptococcus bacteria, and Pseudomonas LPS) versus 50% mortality (Pseudomonas bacteria) versus 0% mortality groups. Blood samples were obtained at 24 hours (red arrow), a time prior to appreciable mortality. (B) Clustering of the 219 probe sets that differentiated the five treatment groups separated the probe sets into six clusters. (C) Co-expression network analysis of these six mouse gene clusters were used to explore the gene expression cartography of the leukocyte response to pneumonia. Genes in common with the human coexpression network ( Figure 3B) are circled. (D) Principal components (PC) analysis of an algorithm-selected subset of the 219 probe sets whose microarray-measured RNA abundance in leukocytes isolated 24 h after the onset of pneumonia. PC2 appeared to explain in part expression signal variance due to mortality rates, while PC3 explained in part the variance due to type of insult. doi:10.1371/journal.pone.0001564.g001 (N = 18, P = 6610 27 ). Clusters 4-6 were enriched for genes that encode nuclear proteins (N = 24, P = 9610 24 ). Principal components analysis of the microarray-derived transcript abundances of the genes selected based on cross-validation clearly differentiated the five experimental groups ( Figure 1D) based on the 7-day mortality (principal component 1) and the type of agent used (Gram positive, sterile, Gram negative, principal component 2).

Clinical study-training cohort
After being mechanically ventilated for .48 h, 27 patients were enrolled into the study. The first 20 patients enrolled were assigned arbitrarily to a training cohort; the other 7 were assigned to a validation cohort. Blood samples were taken at ,48-hour intervals during the study period and then separated into plasma and leukocytes (see below). Of the 20 patients in the training cohort, eight patients either were extubated without developing VAP or withdrew from the study. Of the 12 patients in this cohort who developed microbiologically-confirmed VAP (Table 1), 11 met our analysis criteria of having samples before and after the attending physician made a diagnosis of VAP (that is, one patient, #9, was excluded from analysis for developing VAP on the study entry day). Clinical pulmonary infection scores (CPIS) increased in all 11 patients coincident with the diagnosis of VAP. Three of the patients were culture positive for a Gram positive agent (S. aureus) and the remaining eight patients were culture positive for one or more Gram negative agents. In every case, initial intravenous antibiotic therapy was appropriate for the cultured organism responsible for VAP, based upon cultured organism susceptibilities. Nine of the 11 patients developed VAP 3-6 days after enrollment in the study (''early VAP'') while two patients developed VAP after prolonged mechanical ventilation (''late VAP'', Table 1). All of the patients survived and were discharged from the ICU. Patient-specific timelines were aligned for analysis by assigning ''day 0'' to the day that the attending physician diagnosed VAP.
RNA isolated from patient samples was of high quality and hybridizations met standard performance criteria (vide supra). To assess whether the genes identified in the murine model conveyed information in the patient study, the microarray abundance of the human orthologs of the 219 genes that distinguished the murine pneumonias were numerically analyzed. Principal components analysis of the average RMA-normalized expression levels of these 109 ortholog genes resulted in gene expression trajectories that described the cohort of patients as they developed, were treated for, and recovered from VAP ( Figure 2A). Trajectory translocation along the X-coordinate (principal component 2) appeared to be informative with regard to the onset of VAP -beginning immediately before the diagnosis and ending approximately six days after appropriate antibiotic therapy was initiated. Principal components analysis of plasma cytokine abundance in these patients showed a qualitatively similar trajectory, but with large error bars ( Figure 2B). Nevertheless, translocation along the Xcoordinate (principal component 2) again appeared to coincide with the onset of VAP.
Independent analysis of patient microarray data resulted in the identification of 85 probe sets whose abundance changed significantly during the course of VAP (Table 2). Of the 109 human orthologs that were used to calculate the trajectories shown in Figure 2A, 5 probe sets (4.6%) were present in the list of human probe sets (lactotransferrin, cathelicidin antimicrobial peptide, phospholipid scramblase 1, inhibin beta A, and hydroxyprostaglandin dehydrogenase 15-(NAD)). Network analysis found that the expression behavior of these 85 genes segregated into four clusters ( Figure 3A). Transcript abundance in clusters 1 and 2 generally increased and transcript abundance in clusters 3 and 4 generally decreased around the time of VAP diagnosis.
Principal components analysis of the microarray expression profiles of these 85 genes defined a common response to pulmonary infection complicating critical illness ( Figure 3C). Importantly, trajectory translocation in the X-coordinate (principal component 2) occurred days prior to the clinical diagnosis of VAP. In addition, the informational content of 85 genes chosen at random was determined iteratively 100 times for the first 11 VAP patients. As shown in Figure 3D, only the 85 genes identified by EDGE as significant (FDR#0.10) produced a trajectory; the other sets of genes were scattered randomly around the origin of the graph.
Three patients exhibited contrary gene expression profiles within two of these four clusters ( Figure 3A). Patients 1, 7 and 11 showed decreased expression of genes in cluster 2 and patients 7 and 11 showed increased expression of genes in clusters 3 and 4. Based on patient demographics (Table 1 and data not shown), the only clear difference between patients 7 and 11 and the other patients in the study is that these two patients developed VAP later in their ICU course (study day 18 and 14 respectively) whereas the remaining patients developed VAP between study days 3 and 6. This was also evident in PCA analysis of the 11 individual trajectories (data not shown). The changes in transcript abundance of selected genes were validated by quantitative RT-PCR ( Figure  S1). Finally, we tested whether host ethnicity, host gender, host age, or the cell wall phenotype of the infectious agent had an effect on the number of informational genes. The most informational of these demographic variables was host ethnicity (Figure 4).
We observed in Figure 3C that the aggregate riboleukogram variance in principal component 2 decreased as patients recovered from acute infection. This finding suggested that principal components analysis of microarray-measured gene expression described an attractor, as gene expression time series can be described in terms of dynamical system theory as trajectories in the phase space defined by the main principal components. By plotting the change in PC2 against PC2 over time, we found indeed that the mapped gene expression information appeared to converge toward a common point, suggesting that PC2 represents the expression of the infection-inducible genes ( Figure 5). Consistent with differences in patient age, gender, ethnicity, preexisting co-morbidities, and nature of injury insult, each patient's individual trajectory started at a different point and described a patient-specific arc (data not shown).

Validation cohort
A second cohort of 7 patients was analyzed to evaluate the informational value of the 85 genes that were identified in the first 11 patients with VAP. Two of these 7 additional patients were diagnosed with ''late'' VAP by the attending physician, while another 2 cases were described by the attending as ''possible'' VAP ( Table 1). The individual riboleukograms for these 7 patients demonstrate the existence of immune recovery (basins of attraction) as well as the heterogeneity of the host response. In general, the individual riboleukograms follow a path moving from left to right, that is, from critical illness to recovery (Figure 6, green and red shaded areas, respectively). The development of an infectious complication is typically associated with a change in riboleukogram trajectory. For example, the paths of patients 13,14,15,16, and 17 change directions abruptly coincident with changes in clinical status and concern for VAP or sepsis (see   Figure 6A, but have a similar shape and slope. Both patients 13 and 18 had intracranial hemorrhage. Patient 13 was not diagnosed with VAP (no infiltrate on CXR) but was treated with antibiotics for a fever of 39.4uC and WBC of 31,300 (day 5), tracheal secretions that subsequently grew out Acinetobacter and Stenotrophomonas (CPIS 6), and concern for catheter-related sepsis.
In Figure 6B, the aggregate riboleukogram for the training cohort of VAP patients (aligned for day of VAP diagnosis, see also Figure 3C) is compared to the aggregate riboleukogram for all patients aligned for day of study entry (that is, both training and validation cohorts, N = 11+7, irrespective of VAP day). Again noted are the PCA domains of critical illness and recovery. The aggregate VAP riboleukogram diverges from the aggregate critical illness riboleukogram and rejoins the latter at the point of recovery.

Discussion
Using a bench-to-bedside approach, we have implemented a mouse model of pneumonia and found that RNA abundance profiles obtained from blood samples taken prior to appreciable mortality were able to distinguish between the two variables tested in the assay: lethality of the insult and type of infectious agent. These data extend our previous observations in a mouse model of abdominal sepsis, wherein microarray-measured expression profiles from circulating leukocytes distinguished between infectious and non-infectious etiologies of systemic inflammation in a deidentified cohort [15]. Thus, the mouse circulating leukocyte transcriptional response to infection can not only distinguish between infectious and non-infectious inflammatory insults, but also the type of infectious agent and its associated mortality. Network analysis suggested that the pneumonia-induced transcriptional changes reprioritized mouse leukocytes for the initiation of an immune response, the transcriptional regulation of intracellular signaling cascades, and the induction of numerous transcription factors and other nuclear genes.
Results from the mouse model suggested that the transcriptional activity of buffy coat-isolated cells may be used to monitor the onset of and recovery from acute infection. We tested this hypothesis by calculating principal components using microarray expression profiles of RNA isolated from mechanically ventilated patients at risk for pneumonia. Initially, we examined the behavior of the human orthologs to the genes identified in the mouse pneumonia study. The onset of acute infection corresponded with translation along PC 2 in Figure 2. Translation along PC2 ceased 5-6 days after appropriate antibiotic therapy in the patients was started, consistent with recovery. These data show that there are specific transcriptional programs instituted by circulating immune cells during acute infection which have diagnostic potential in the setting of critical illness. Principal components analysis of plasma cytokine abundance generated a qualitatively similar trajectory; however, in line with previous reports, plasma cytokine abundance (including procalcitonin) was insufficient to diagnose acute infection in this small cohort of critically ill patients. As with other tissues, changes in RNA abundance observed in circulating leukocyte do not necessarily reflect changes in protein abundance, and vice versa. While principal components analysis of informational murine genes in authentic human disease showed there was a conserved and informational peripheral leukocyte transcriptional response to localized infection, the information contained in those genes would not appear to be more useful than current clinical  criteria (that is, the translation in PC2 did not begin until the day the attending physician made the diagnosis). However, by explicitly accounting for variance over time, a set of 85 genes were identified subsequently in our first 11 VAP patients whose microarray expression levels changed consistently before the clinical diagnosis of VAP. These 85 genes clustered into four groups with their abundance either increasing or decreasing throughout the 7-day window bracketing the onset of infection ( Figure 3A). There was a significant association between genes known to play key roles in defense against bacterial pathogens and those genes that increased in apparent abundance (cluster 1 and 2 probes sets) coincident with the diagnosis of VAP. Consistent with activation of the host response to pneumonia, all of the genes with the ''defense against bacteria'' ontology were found in cluster 2 and encode primarily granulocytic, antimicrobial proteins, and adhesion molecules. In contrast, the genes that decreased in apparent abundance (clusters 3-4) showed different behavior depending on whether the patient developed VAP early or late in the study. Although no consistent biological theme emerged from this list of the 59 transcripts, this finding provides some insight into the transcriptional basis of differences in the critically ill host's response to early-versus late-onset VAP [29]. Of interest, the RNA abundance of five genes were altered by pneumonia in both the mouse and human gene sets (Figures 1C and  3B): lactotransferrin (LTF), cathelicidin antimicrobial peptide (CAMP), phospholipid scramblase 1 (PLSCR1), inhibin beta A (INHBA), and hydroxyprostaglandin dehydrogenase 15-(NAD) (HPGD). Pathway analysis indicates that these five genes connect  in a network rich with interactions between important mediators of inflammation, as seen in Figure S2 and Table S2, (Ingenuity Pathway Analysis, Ingenuity Systems, Redwood, CA). Principal components analysis of the microarray expression data for the 85 transcripts generated a curve with qualitative similarities to curves derived from either cytokines or the 109 human orthologs; however, the patient-to-patient variance in the measured abundance of these genes was significantly smaller and translation along the second principal component preceded the diagnosis of VAP by 24 to 72 hours ( Figure 3C). In contrast, equivalent analysis of randomly selected sets of genes provided no information about the host response to critical illness complicated by pneumonia.
Macroscopically, pneumonia is diagnosed when there are symptomatic changes in clinical status, as manifested typically by increased CPIS scores. However, microscopically pneumonia occurs at the transition from colonization to infection concomitant with failure of multiple host barriers protecting the bronchoalveolar epithelium. These events lead to tissue exposure to bacteria and bacterial products, and in many cases to toxicosis and bacteremia. Our data suggest that during VAP, circulating leukocytes are exposed to bacterial products 24-72 hours prior to diagnosis by an attending physician. Initiation of antibiotic therapy earlier during this time window is expected to significantly improve outcome [16]. These findings were further evaluated in a The time-dependent behavior of these genes was classified into four clusters. Shown is the normalized abundance of the average cluster member for each of the five microarrays that bracket the attending physician's diagnosis of VAP for each patient. (B) Co-expression network analysis of these four clusters was used to generate the gene expression cartography of the human blood response to acute infection superimposed on critical illness. Clusters 1 and 2 are tightly associated with one another, as are clusters 3 and 4. The red circles identify the 5 genes in common with the mouse coexpression network ( Figure 1C). (C) Principal components analysis of the abundance of these 85 genes in the training data set (11 patients with VAP). PC1 (not shown) represents a constant bias term, PC2 and PC3 are shown. The arrow indicates where the attending physician's diagnosed VAP. The green and red circles indicate the points where the patients entered and exited the study, respectively. (D) Principal components analysis of microarray data generated by 100 iterations of randomly chosen sets of 85 genes, all plotted on the same two axes. The randomly chosen sets are not informational for healing from critical illness complicated by VAP. The only set that describes a discernable path is the list of 85 genes derived from EDGE analysis (inset magnified view, same riboleukogram data as in panel C). doi:10.1371/journal.pone.0001564.g003 small validation patient cohort, confirming the value of the 85gene riboleukogram to quantitate the host response to and recovery from critical illness. Individual riboleukograms changed direction coincident with clinical findings of pneumonia or sepsis (frequently days before the clinical diagnosis and maximal CPIS scores). Nevertheless, marked heterogeneity was observed in some patient responses that could not be linked to monitored variables (e.g., patients 13 and 18, Figure 6A). In these patients, perhaps the well-described influence of the underlying acute illnesses (pulmonary contusion and intracranial hemorrhage) on host immune responses could provide part of the explanation [30,31]. As patients healed from critical illness, the riboleukograms converged, a finding consistent with the existence of an immunological ''attractor'' state. Comparison of the aggregate VAP versus critical illness riboleukograms indicated that a large portion of the signal for the 85 genes identified is a reflection of recovery from critical illness (PC2 in Figure 6B). Nevertheless, VAP signal became evident among the ''noise'' of critical illness once the individual riboleukograms were aligned for day of VAP diagnosis, evident as a deflection of the VAP riboleukogram upward in PC3 ( Figure 6B). Thus, we submit that riboleukograms are a molecular analytical tool with substantial potential to improve diagnostics, prognostics, and our understanding of the host response to critical illness complicated by acute infection.
An important consequence of our observation that patientspecific riboleukograms converged is that the variance in leukocyte gene expression for these 85 genes decreased significantly over the time course in patients with VAP. Studies of physical systems suggest that the stable states of networks can be represented as attractors, a set of points in the phase space to which the genetic network evolves over time [32,33]. In particular, every trajectory initiated within the bounds of an attractor terminates inside the attractor. Recently this hypothesis has been confirmed experimentally in vitro [11,34]. Our findings suggest that the immune system of a critically ill patient who recovers from disease returns back to a stable state and that the immune response trajectory can be a considered formally as a dynamic system. Borrowing from concepts in the physical sciences [35], we hypothesized that an individual's baseline immune system is within a basin of attraction prior to injury, nosocomial infection perturbs that system, and after the infection cleared the system would return to it's initial attractor. In phase space analysis, we found that the patientspecific trajectories appeared to converge and that the onset of pneumonia 'pushes' the gene expression data away from that basin-a perturbation that is ameliorated by appropriate antibiotic therapy. This result provides what we believe is the first evidence in patients for a basin of attraction for the genetic network associated with immunological health.
Because of the heterogeneity of this patient population (for example, age, gender, pre-existing health status, type and severity of critical illness), it is not surprising that the error bars on the PC analysis are initially quite large. The end-points of the trajectories, however, appear to converge in a smaller space, suggesting that the patients' immune systems are returning to a health attractor (Figures 5 and 6). As observed in the validation cohort, patientspecific trajectories are not smooth curves, but appear to have inflection points that may correlate with hospital intervention or the onset of infection. Our data further suggest that differences in host ethnic backgrounds and gender are more important than differences in infecting organism as determinants of the host leukocyte transcriptional response. These observations are consistent with a recent study comparing GeneChip signal from cells lines derived from Asians and Caucasians, indicating differences in the gene expression levels of 25% of the 4000 genes studied in these two ethnic groups [36]. Moreover, as the rate of sepsis, the youngest age of sepsis onset, and the sepsis mortality rate are highest for African-American males [37], our data suggest that further study of riboleukograms is indicated to gain insight into these health disparities. These profiles also suggest that there is a transcriptional component that is conserved across that diversity of the human population.
In summary, our analysis demonstrates the plasticity of the blood leukocyte response to bacteria in vivo, extending previous in vitro studies [10,11] into the clinical realm. For the first time, we provide evidence that riboleukogram gene expression analysis can be applied to a heterogeneous clinical population to monitor the host response to and recovery from critical illness complicated by acute infection. Moreover, we identify new, conserved gene targets that appear to be informational for recovery at the transcriptional level, many of which are involved with granulocyte maturation and chemokine (not cytokine) responsiveness. This conclusion is supported by the relative weak diagnostic information provided by the plasma cytokine data, consistent with conclusions reached at recent a consensus conference [16]. In particular, plasma procalcitonin levels were not informational. There are, however, important limitations to our study, most notably, the preliminary nature of the work, the small number of patients, and the well-known difficulty of ruling out false-positive and false-negative diagnoses of VAP based on clinical parameters. This uncertainty is both the motivation for (and challenge of) developing novel molecular diagnostics for VAP. In addition, there is no consensus on how to leverage the dynamic nature of the clinical, RNA, and protein data collected to build hybrid models that improve diagnostics and prognostics. Thus, we conclude that our data demonstrate the technical feasibility and clinical potential of the riboleukogram approach, but proof of clinical utility will require further study.
Thus, as graphs of myocardial electrical information (electrocardiograms) were tapped over a century ago to provide an objective means to aid heart diagnostics, we submit that riboleukograms will aid in the diagnosis and prognosis of acute infectious and inflammatory disease. The diagnostic potential of riboleukograms is supported by two very recent independent reports that corroborate our mouse data ( Figure 1D), indicating  Table 1). Note that the largest number of probe sets were associated with differences in ethnic background (African-American compared to Caucasian). Of note, as opposed to the mouse response, there we no genes that differentiated between the human response to bacterial cell wall products (Gram negative versus Gram-positive bacteria), suggesting that signal variance due to ethnic background, gender, and age is greater than that due to infecting organism. doi:10.1371/journal.pone.0001564.g004 that circulating leukocyte RNA signatures in patients differentiate between the host responses to sterile versus infectious causes of systemic inflammation and between the host response to Gramnegative versus Gram-positive pathogens [38,39]. Prospective clinical trials are indicated to validate our results and determine the value of this new technology; to optimize gene selection methods that account for differences in patient ethnicity, gender, and age; and to develop computational approaches that integrate clinical and molecular data to improve diagnostics [40]. Figure S1 The changes in transcript abundance of selected genes from Table 2 Figure 3C). The other 7 curves are the individual riboleukograms of the patients in the validation cohort. The inset magnifies the trajectories of patients 13-17 (see Table 1) and demonstrates abrupt changes in riboleukogram course typically coincident with an increase in CPIS score (first occurrence of maximal CPIS value is indicated by the arrows). The paths of patients 13, 18 and 19, are atypical (see text for additional details). (B) The aggregate 11 patient VAP riboleukogram (black curve, same as panel A) is compared to the aggregate riboleukogram of all patients aligned by study day (that is, training and validation cohorts, irrespective of VAP day of diagnosis, dotted blue curve). Note that the VAP riboleukogram deviates from the ''critical illness'' riboleukogram (black arrows) prior to VAP diagnosis (lighting bolt), but after treatment, the VAP riboleukogram converges with the critical illness riboleukogram at the point of recovery. The green and red circles indicate where the patients entered and exited the study, respectively. doi:10.1371/journal.pone.0001564.g006