Trauma in silico: Individual-specific mathematical models and virtual clinical populations

See allHide authors and affiliations

Science Translational Medicine  29 Apr 2015:
Vol. 7, Issue 285, pp. 285ra61
DOI: 10.1126/scitranslmed.aaa3636


Trauma-induced critical illness is driven by acute inflammation, and elevated systemic interleukin-6 (IL-6) after trauma is a biomarker of adverse outcomes. We constructed a multicompartment, ordinary differential equation model that represents a virtual trauma patient. Individual-specific variants of this model reproduced both systemic inflammation and outcomes of 33 blunt trauma survivors, from which a cohort of 10,000 virtual trauma patients was generated. Model-predicted length of stay in the intensive care unit, degree of multiple organ dysfunction, and IL-6 area under the curve as a function of injury severity were in concordance with the results from a validation cohort of 147 blunt trauma patients. In a subcohort of 98 trauma patients, those with high–IL-6 single-nucleotide polymorphisms (SNPs) exhibited higher plasma IL-6 levels than those with low IL-6 SNPs, matching model predictions. Although IL-6 could drive mortality in individual virtual patients, simulated outcomes in the overall cohort were independent of the propensity to produce IL-6, a prediction verified in the 98-patient subcohort. In silico randomized clinical trials suggested a small survival benefit of IL-6 inhibition, little benefit of IL-1β inhibition, and worse survival after tumor necrosis factor–α inhibition. This study demonstrates the limitations of extrapolating from reductionist mechanisms to outcomes in individuals and populations and demonstrates the use of mechanistic simulation in complex diseases.


Traumatic injury represents the most common cause of death for young people and is a significant source of morbidity and mortality for all ages (1, 2). Hemorrhage and trauma, like sepsis, induce an acute inflammatory response involving a coordinated mobilization of numerous inflammatory cells and circulating mediators, with repercussions for all organ systems (36). Acute inflammation after injury is complex: A robust, early tumor necrosis factor–α (TNF-α) response appears to be crucial for the survival of blunt trauma patients and experimental animals, although in these same patients elevated interleukin-6 (IL-6) levels are associated with higher morbidity and mortality (7, 8). Further, there is a broad-based elevation of multiple innate immune and T cell inflammatory pathways after blunt trauma, which is in turn associated with complications such as multiple organ dysfunction syndrome (MODS) and nosocomial infection (9, 10). These and other studies suggest that an appropriately adaptive acute inflammatory response is beneficial and crucial for tissue recovery after trauma; however, if inappropriately exaggerated or sustained, the inflammatory response can also compromise healthy tissue, further exacerbating inflammation through positive feedback loops and causing further impairment to tissue integrity (11).

The outcomes landscape in blunt trauma has shifted from high mortality (now only ~5%) to MODS, nosocomial infection/sepsis, prolonged length of stay (LOS), and other complications (10, 12). Therefore, most patients have an indeterminate outcome upon admission but are highly likely to survive. Thus, the current challenge lies in stratifying likely morbidity to better implement surveillance, prophylaxis, and resource use. We have used both quasi-mechanistic, data-driven modeling approaches (9, 10, 13, 14) as well as biologically motivated mechanistic mathematical models based on ordinary differential equations (ODEs) (1518) to gain insights into the whole-organism response to traumatic injury. We have explored the use of multicompartment ODE models for explaining individual-specific trajectories of inflammation and organ dysfunction in outbred swine subjected to endotoxemia (19). We and others have also used mechanistic mathematical models as a tool for carrying out in silico clinical trials for sepsis (2022) and wound healing (23, 24). We calibrated these previous studies in sepsis qualitatively to reflect major clinical outcomes and reproduce key features of previous clinical trials (20, 21). We have also shown that mechanistic computational models can recapitulate the dynamic inflammatory trajectories and ultimate outcomes in individual human subjects subjected to localized tissue trauma (24). To date, however, there have been no reported studies in which computational models were either calibrated or validated against prospectively obtained data in patients. In addition, the impact of trauma on inflammation and organ dysfunction has not been explored in simulated clinical cohorts. More fundamentally, there has been much interest in, but little clinically oriented work on, the question of whether it is possible to extrapolate from reductionist mechanism to the responses of individual patients or patient cohorts (25).

Herein, we extended our previously published, pig-specific two-compartment ODE model (19) to include an additional tissue compartment in which trauma occurs (Fig. 1). Using this augmented model, we tested its use for predicting outcomes for large cohorts of human blunt trauma patients.

Fig. 1. Schematic of three-compartment mathematical model of human acute inflammation.

(A) Human trauma model structure. The insult of trauma and pathogen/endotoxin induce systemic inflammation across the three compartments. (B) Progression of inflammation and positive feedback loops in the model after traumatic injury insult. Clinical outputs are shown in red.


Demographics of initial cohort of blunt trauma survivors

A central goal of the current study was to generate patient-specific mathematical models of inflammation, organ dysfunction, and outcomes. To do so, blood samples from 33 trauma survivors (19 males and 14 females) of motor vehicle accidents (26), falls (5), or motorcycle accidents (2) were studied. The overall demographics of the patients were as follows: age, 44 ± 3 years; Injury Severity Score (ISS), 24 ± 3; intensive care unit (ICU) LOS, 7 ± 1; total LOS, 14 ± 2; and number of days on a ventilator (vent days), 4 ± 1. Plasma was sampled within the first 6 hours after trauma, to a maximum of 30 days after trauma.

Patient-specific trauma models

We previously reported a two-compartment (“blood” and “lung”), equation-based model of inflammation, organ (dys)function (including lung physiology and mean arterial pressure), and whole-animal health status, which was calibrated with data on endotoxemic swine (19). We obtained similar data from the 33 blunt trauma patients as described above, including plasma IL-1β, IL-6, IL-10, IL-12p40, TNF-α, and NO2/NO3 from blunt trauma patients. We augmented the pig-specific, two-compartment mathematical model (19) by adding a third (“tissue”) compartment in which injury takes place, with resultant effects on the blood, lungs, and whole-animal health status (Fig. 1; see Materials and Methods and Supplementary Materials for details). We calibrated the model initially using data obtained from the individual survivors of blunt trauma described above. Figure 2 (A to C) shows the results of three patient-specific mathematical models, in which a small subset of model parameters was modified to match trajectories of IL-1β, IL-6, IL-10, IL-12p40, TNF-α, and NO2/NO3 in individual trauma patients. The inflammation and subsequent organ dysfunction were initiated in each virtual patient by a simulated initial injury. This initial simulated injury was patient-specific, scaled in magnitude to the actual ISS assigned to that particular (real) patient.

Fig. 2. Patient-specific model fits.

(A) Model time-course prediction and data for trauma patient with low insult (ISS = 5). (B) Model prediction and data for trauma patient with intermediate/high insult (ISS = 30). (C) Model prediction and data for trauma patient with severe insult (ISS = 50).

Additionally, model calibration revealed that, in certain cases, modulation of a key inflammatory mediator such as IL-6 (8) could determine the difference between survival and nonsurvival of a trauma patient (see parameter sensitivity analysis below and in the Supplementary Materials). This was a result of a bifurcation in the general “damage” predictor, whereby damage reached a self-sustaining level that did not return to its original baseline state. The specific parameters used to obtain these outcomes, as well as their descriptions, can be found in table S1. As an example, fig. S1A shows two virtual patients who are identical in all ways, except for the IL-6 production phenotype. In this case, the two patients have roughly the same predicted outcome for the same traumatic insult. Conversely, fig. S1B shows two analogous virtual patients who have completely different predicted outcomes because of differences in the IL-6 production phenotype.

Virtual population outputs

We next used these patient-specific mathematical models to create a large population of virtual blunt trauma patients. Two of the most important clinical outcomes after traumatic injury are the overall degree of injury in a given patient (which can be approximated by the ISS and which correlates to susceptibility to experience multiple organ dysfunction) and the amount of time taken to resolve damage and return to a baseline state, which correlates to ICU LOS. We hypothesized that these two clinical outcomes could both be associated with the general damage term defined in the mathematical model. Specifically, area under the damage curve was equated with global susceptibility to MODS, and time until the damage curve resolves to 15% of peak level was mapped to recovery time/ICU LOS.

Results of the population-level simulation indicate a wide range of clinically relevant outputs that are not strictly dependent on ISS, but nonetheless show the general trend that more severe injury leads to worse clinical outcomes. Figure 3 (A to B) shows the distributions of the two outcomes for the three levels of injury severity used for these analyses. Notably, most virtual trauma patients resolve injury somewhere between 2 and 10 days after injury, consistent with observed ICU LOS in the original cohort of 33 patients (7 ± 1 days). Additionally, the ability of the model to reproduce a wide range of outcomes—reliant on subsequent, variable, simulated trajectories of circulating inflammatory mediators as well as the initial traumatic injury—is consistent with clinical observation (26). Figure 3C shows the predicted mean total LOS as a function of ISS.

Fig. 3. Model-predicted and observed multiple organ dysfunction and ICU LOS.

(A) Distributions of “Total damage” or susceptibility to multiple organ dysfunction for various levels of traumatic injury severity. (B) Distributions of “Damage recovery time” or ICU LOS for various levels of traumatic injury severity. (C) Predicted mean total LOS as a function of ISS. (D) Verification of model predictions shown in (C), in a separate cohort of 147 blunt trauma patients.

To validate this prediction, we used a separate cohort of 147 blunt trauma patients that did not include the original 33 patients used to generate the 10,000 virtual patient cohort. This validation cohort, all survivors (99 males and 48 females), included victims of motor vehicle accidents (100), falls (25), and motorcycle accidents (23). The overall demographics of the validation cohort patients were as follows: age, 41.1 ± 1.2 years; ISS, 20.8 ± 0.9; ICU LOS, 6.6 ± 0.6; total LOS, 12.8 ± 0.8; and vent days, 3.5 ± 0.5. Indeed, distributions of ISS, total LOS, and ICU LOS from the 147-patient validation cohort largely matched those predicted by the ODE model (Fig. 3D).

Predicted damage bifurcation and trauma-induced mortality

Having shown that the predicted damage trajectory of a virtual trauma patient could vary from that of survival to that of death (fig. S1), we next sought to determine whether the simulated cohort of virtual trauma patients could exhibit this same potential for survival or death after injury at any given severity. Although the mathematical model was calibrated with only data from survivors of blunt trauma, simulation results of the 10,000 virtual trial population included patients who never resolved the trauma-induced damage. This dichotomy is a result of a bifurcation in the damage equation that was not programmed into the mathematical model, the result of which was that output continues to increase or becomes stable at a nonzero value rather than return to baseline. We equate this latter state with trauma-induced mortality. Figure S2 (A to C) includes damage trajectories for each of the 10,000 virtual patients for the three levels of traumatic insult. Table 1 quantifies the number of predicted nonsurviving patients for each injury severity level. In the absence of explicit training of the model to this target, the predicted 30-day survival rate of 96.5% (that is, mortality rate of 3.5%) for all virtual human trauma patients was similar to clinical observations for cases that do not require additional surgery [that is, mortality rate of 4.5% (12)].

Table 1. Summary of simulated 30-day survival rates after simulated blunt trauma.
View this table:

Responsiveness analysis and verification in a large cohort of blunt trauma patients

We next determined the parameters to which the mathematical model was most sensitive, hypothesizing that these parameters would represent biologically meaningful aspects of the response to trauma. High relative responsiveness indicates model sensitivity to parameter perturbations, indirectly indicating the principal drivers of trauma-induced inflammation for various ISS in the context of the abstractions inherent in this mathematical model. We focused on the same clinically relevant model outputs as in Fig. 3 [“Damage recovery time” (analogous to ICU LOS) and “Total damage” (analogous to severity of MODS)]. The relative responsiveness effects on these model outputs after alterations in cytokine production rates are summarized in Fig. 4 and suggest that IL-1β and IL-6 are key drivers of outcome for the virtual trauma patients.

Fig. 4. Model-predicted impact of specific inflammatory mediators on simulated ICU LOS and multiple organ dysfunction.

(A) Predicted responsiveness of “Damage recovery time” (that is, ICU LOS) to simulated degree of cytokine production. (B) Predicted responsiveness of “Total damage” (that is, total LOS) to simulated degree of cytokine production).

This responsiveness analysis further suggested that the proinflammatory mediators IL-1β and TNF-α are more important in driving outcomes at low injury severity (Fig. 4). In contrast, the mediators that can exert at least some anti-inflammatory effect, namely, IL-6 and IL-10, were more important in driving outcomes at moderate and high ISS. This difference indicates that the predominant determinants of inflammation are not static, even within an individual patient, but vary with the traumatic insult to which the patient is subjected.

A complex role for IL-6 in trauma outcomes inferred from in silico trials and verified in a large cohort

Having established that IL-6–related parameters in our ODE model are a key determinant of mortality and survival in silico in the highest-severity patients, we examined this observation in greater detail to verify key aspects in separate cohorts of blunt trauma patients. We first examined the predicted distribution of IL-6 concentrations as a function of ISS. Figure 5A shows the predicted dynamics of IL-6 over the first 6 days after injury, and Fig. 5B shows the actual dynamics of IL-6 as function of ISS in the 147-patient validation cohort. Figure 5C shows a comparison of predicted IL-6 area under the curve (AUC) as a function of ISS, and Fig. 5D shows the qualitative concordance of IL-6 AUC as a function of ISS for the 147-patient validation cohort.

Fig. 5. Model-predicted and observed IL-6 dynamics.

(A) Predicted IL-6 AUC as a function of injury severity. (B) Actual IL-6 AUC in a validation cohort of trauma patients. (C) Predicted IL-6 AUC as a function of ISS. (D) Qualitative concordance of IL-6 AUC as a function of ISS for the 147-patient validation cohort.

The predicted influence of IL-6 levels on individual outcomes in silico (Fig. 4) led to the seemingly logical hypothesis that the outcomes of virtual patients in the cohort of 10,000 that produce higher levels of IL-6 in response to trauma would be worse than the outcomes of the virtual patients that produce low levels of IL-6. In silico, IL-6 AUC was predicted to vary both as a function of ISS and of the propensity to produce IL-6 in the virtual patients (Fig. 6, A to C).

Fig. 6. Model-predicted and observed impact of IL-6 production on clinical outcomes.

(A to C) Predicted circulating IL-6 concentrations, as a function of propensity to produce IL-6, in virtual blunt trauma patients subjected to low/intermediate (A), intermediate/high (B), or severe (C) injury severity. (D) In a 98-patient validation cohort, plasma IL-6 levels are shown as a function of IL-6 SNP. Plasma IL-6 levels in healthy, 11 noninjured volunteers of similar age and gender as the validation cohort. (E) Propensity to produce IL-6 is predicted by the mathematical model to have no statistically significant association with clinical outcomes. (F) In a validation cohort, blunt trauma patients with high IL-6 SNPs exhibit no statistically significant differences in clinical outcomes as compared to patients with low IL-6 SNPs.

To test the hypothesis that patients predisposed to make higher IL-6 levels in response to traumatic injury would have worse outcomes than patients predisposed to make low levels of IL-6, we examined a subcohort of 98 blunt trauma patients derived from the 147-patient cohort that contained subjects with both high–IL-6 single-nucleotide polymorphism (SNP) (G/G, G/C) and low–IL-6 SNP (C/C). There was no statistically significant difference in age or ISS between the high– and low–IL-6 SNP groups. Analysis of the dynamics of plasma IL-6 confirmed that circulating IL-6 varied as a function of IL-6 SNPs (Fig. 6D) (reference values of healthy, uninjured human volunteers are given for comparison and as a baseline value). However, as shown in Fig. 6E, model simulations predicted no statistically significant differences in damage recovery time (a proxy for ICU LOS as described above) or severity within the moderate/high ISS virtual cohort as a function of IL-6 production. These results were validated in the 98-patient subcohort: as seen in Fig. 6F, there were no statistically significant differences in total LOS, ICU LOS, or requirement for mechanical ventilation between the high– and the low–IL-6 SNP groups of patients with moderate to severe injury. Thus, a seemingly intuitive hypothesis derived from the simulations of individual patients and supported by numerous prior studies was both predicted and demonstrated to be incorrect at the level of the population.

In silico randomized clinical trials in blunt trauma

The striking discrepancy between reductionist mechanisms and the behavior of patient cohorts led us to study this population behavior further in the context of simulated clinical trials. To date, there is no approved therapy that modulates the immuno-inflammatory response in trauma. We (2124) and others (20) have shown the potential use of in silico clinical trials for suggesting the characteristics of patients that might benefit, or be harmed by, therapies for inflammatory diseases. Model sensitivity analysis suggested differential, injury- and patient-dependent sensitivity of outcomes based on the propensity to produce IL-1β, IL-6, and TNF-α. Accordingly, we sought to leverage the three-compartment mathematical model described above to carry out in silico clinical trials of selective inhibition of each of these three cytokines as a single therapy administered in a randomized fashion. Using methods similar to those we have described previously (21) (see Materials and Methods), we compared the effects of neutralizing each of these inflammatory mediators to the effects of a placebo. This approach allowed us to predict either benefit or harm of a given intervention (Table 2) and suggested that inhibition of IL-1β would have little effect, that inhibition of IL-6 could provide a modest benefit, and—seemingly paradoxically but in agreement with our previous studies in blunt trauma (7)—that inhibition of TNF-α would worsen mortality. These results again highlight the often non-intuitive behavior of populations, both in silico and in reality.

Table 2. Summary of in silico randomized clinical trials in severe blunt trauma.

Simulated clinical trials were carried out as described in the Materials and Methods. Percent survival was evaluated at 30 days of simulated time.

View this table:


Trauma is a leading cause of death worldwide, often due to the immune dysregulation that accompanies complications such as nosocomial infection and MODS (2729). Properly regulated inflammation allows for timely recognition and effective reaction to injury. However, inflammation can become self-sustaining if the insult is of too great a magnitude or duration, or if the genetic/epigenetic makeup of the individual predisposes toward maintaining elevated inflammation (11, 30, 31). For the individual patient, the envisioned role of personalized medicine is the prediction of outcomes for the individual, and systems biology is generally described as the main tool required for accomplishing this goal (32). Nowhere is this goal more challenging than for rapidly evolving inflammation and organ dysfunction after injury (33), although to date, systems biology approaches have not been used clinically in this setting. For populations of trauma victims, the reductionist and fragmented process of the current paradigm for clinical trial design has failed to yield any approved therapy for inflammation and MODS. Thus, there is a pressing need to address the complexity and interconnection between inflammation and physiology, complexity that likely manifests as interindividual heterogeneity (3437). In short, the whole is greater than—or at least different from—the sum of the parts. Herein, we sought to reconcile how cellular- and molecular-level insights about the interplay of inflammation and organ dysfunction in a cohort of individuals (the parts) affect the ultimate behavior of a cohort of patients (the whole).

To accomplish this goal, we turned to dynamic computational modeling. This class of mechanistic models has been used to explain the non-intuitive behavior that often emerges from reductionist biological interactions (3842), to suggest interactions not yet described by experimental data (43), and to address experimental controversies (44). We previously constructed mechanistic computational models based on data from mice (1518, 45, 46), rats (47), swine (19), and humans (24, 48, 49). With this approach, we have shown the potential for predicting individual trajectories of inflammation and tissue damage for localized tissue trauma in humans (24). We have also demonstrated the potential for in silico clinical studies using such computational models (2124).

Mechanistic models are distinct in their use from statistical models. Data-driven and statistical models can help define key regulators of inflammation in hemorrhaged mice (13) and biomarkers of outcome in trauma patients (9, 14). However, they have only a limited ability to predict the outcomes of trauma patients based on circulating inflammatory mediators and even then generally are only accurate in large patient cohorts, not in individual patients (26, 5054).

In response to calls in the literature to gain greater insight from relatively small patient cohorts with computational modeling (25), we extended prior work in large animals (19) to the responses of individual trauma victims. These patient-specific models captured, to varying degrees, individual trajectories of inflammation and organ dysfunction. The initial conditions for each of these patient-specific simulations were scaled relative to the severity of their injury, as quantified by the ISS. The ISS involves a mathematical manipulation of various factors associated with body region–specific injuries; these mathematical transformations were designed in an attempt to use ISS for stratification of risk of subsequent morbidity or mortality. However, although the ISS is a well-accepted method for patient stratification based on the degree of injury, this retrospective scoring system is not used to predict prospectively the outcomes of individual patients except at the extremes of injury severity (55, 56). Herein, we used the ISS as a means for quantifying the magnitude of the initial stimulus for acute inflammation in our simulations; outcome prediction was accomplished by running a given simulation and assessing both specific aspects of pathophysiology (hypotension and lung dysfunction) as well as the overarching “damage” variable.

Thus, although ISS is in essence the initial trigger for inflammation and subsequent dysfunction in our mathematical model, the ultimate prediction of outcome is based on the nonlinear way in which multiple factors interact within a given patient’s simulation. One of the factors that interact with the magnitude of injury is the manifestation of acute inflammation in a given patient, which in turn is driven in part by SNPs in the coding regions affecting the expression of inflammatory mediators. We deployed advanced statistics and data mining techniques on the in silico data set generated in our simulations to identify IL-6 and IL-10 as two key drivers of inflammation after severe injury, whereas IL-1β–related parameters appeared to drive simulated outcomes following mild injury. Thus, our simulations suggest that SNPs that control the levels of these inflammatory mediators may ultimately affect patient outcomes, although likely in complex and indirect ways.

We focused initially on the role of IL-6 after injury. This cytokine has been long recognized as a biomarker of adverse outcome in numerous types of acute inflammation (8), most notably in studies from the large-scale Inflammation and the Host Response to Injury consortium, which showed that the predominant biomarker of morbidity and mortality in males after blunt trauma was IL-6 (57). We first confirmed that trajectories of circulating IL-6 after injury varied as a function of IL-6 SNP, with patients with SNPs coding for high IL-6 production exhibiting higher plasma IL-6 than patients with SNPs coding for low IL-6 production. This result verified an association between genotype and dynamic levels of IL-6.

Our simulations also suggested a sensitivity of damage to IL-6 in patient-specific models: Elevated virtual production of IL-6 could convert a survivor of severe injury to a nonsurvivor on an individual basis. It should be noted that the mathematical model was trained on data of trauma survivors only; thus, it was unexpected (and strikingly different from statistically based models) that the simulation captured such a bifurcation of ultimate outcomes without training on data from nonsurvivors.

Despite reports suggesting improved outcomes in IL-6–deficient mice or in animals subjected to various IL-6 neutralization or retargeting strategies for trauma/hemorrhage (5862), and our own experience with beneficial outcomes of neutralizing IL-6 at the time of resuscitation in hemorrhaged mice (63), implementation of IL-6–directed therapy has proven challenging due to the number of receptors and other binding proteins for IL-6, as well as complex effects of different doses of neutralizing anti–IL-6 antibodies (64). This difficulty may be a result in part of the beneficial, procoagulant role of IL-6 in trauma/hemorrhage (65) or to its dual (pro- and anti-inflammatory) effects (8, 64, 66, 67). Our simulations suggest that, although individual outcomes may be altered by changing the production of or the response to IL-6, the clinical outcomes of trauma patient populations do not vary in a simple way as a function of genetic predisposition to produce IL-6. This observation confirms that individualized treatment must consider both the genetics of the individual as well as the severity of the injury, age, comorbidities, and contribution of multiple inflammatory mediators acting in concert.

The interrelations among inflammatory mediators, genetic predisposition, and organ dysfunction are a key aspect of our modeling approach. The past two decades have seen little success in improving clinical outcomes in sepsis by modulating acute inflammation (68). This clinical experience has led to the suggestion that inflammation may not be a key driver of outcomes in sepsis and hence is not a viable target (69). Our models of dynamic networks in porcine sepsis (70) also suggest that inflammation is connected to organ dysfunction in a positive feedback loop. This interconnection is present in in silico clinical trials of acute inflammation in sepsis, which replicated the failed outcomes of modulating single mediators and predicted the failure of combinations of mediators (20, 21). As discussed above, we carried out an in silico clinical trial of IL-6 neutralization in severe blunt trauma, and our results predict only a slight beneficial effect of this strategy. In contrast, and in support of our previous studies suggesting a beneficial role of early increased TNF-α production after injury (7), simulated neutralization of TNF-α was predicted to result in worse outcomes as compared to placebo.

These results do not mean that inflammation as a whole, or IL-6 or TNF-α specifically, are not viable therapeutic targets. Instead, rather than administering such therapies in a randomized fashion, these therapies are likely to require careful, biomarker-driven stratification of patient subgroups (33). In such an approach, models such as the one described herein could be used to define patients with specific early inflammatory mediator profiles or other related characteristics that indicate they are likely to benefit from immunomodulation. To aid in this patient selection process, virtual populations could be better defined to match observations across a larger array of clinical measurements (for example, white blood cell counts or various physiological or genetic measurements). Ultimately, such models can be combined with patient-specific, adaptive approaches for inflammation control (71) to effect truly personalized solutions for acute inflammatory diseases.

The present study should be evaluated in the context of several possible limitations. First and foremost, like all models, the three-compartment model presented herein is an abstraction of reality, and thus omits multiple mechanisms associated with the response to trauma as well as organs besides the lung that are known to become dysfunctional after injury. This abstraction also means that the number of inflammatory mediators involved, the cells that produce them, and the compartments in which they localize are likewise absent. A second limitation concerns the patient cohort used for the construction and validation of the model. This cohort represents seriously injured individuals, but not so seriously injured that they are unable to survive at least 24 hours after injury. These patients are the ones likely to require ICU care and extensive hospitalization, which was the focus of this study. Also, although we have fit the model to the data and have used the model to predict features outside of the training data set, it is always formally possible that the parameter set used represents a local rather than a global optimum, and thus may make incorrect predictions upon further testing.

Despite these limitations, our findings raise a note of caution with regard to the extrapolation from reductionist analysis of disease pathways to the ultimate clinical outcomes. Biomedical research has been in many cases unable to effectively translate basic mechanistic knowledge into therapeutics, particularly in attempts to understand and modulate “systems” processes/disorders such as critical illness induced by trauma and sepsis (42). There have been numerous advances in defining novel molecules, signaling and synthetic pathways, and gene regulatory networks contributing to these diseases (7274). Nevertheless, these advances were produced and remain in scientific silos, and therefore have not been integrated into a systems-level understanding. An unfortunate consequence is the dearth of available therapeutics for these deadly diseases. Inflammation is a highly robust communication framework; thus, it should not be surprising that single-target approaches have failed, especially given the highly dynamic nature of the response to traumatic injury (10). Accordingly, we suggest that computational modeling can help to improve the design, targeting, dosage, timing, and patient selection for therapies that modulate inflammation in acute illness.


Study design

Sample size for validating model-predicted circulating levels of inflammatory mediators as a function of ISS was determined based on a power calculation based on the expected differences in circulating levels of IL-6 across the three ISS patient subcohorts (low/intermediate, intermediate/high, and severe), as well as the expected variance in these levels within a single patient cohort. This analysis suggested that a sample size of 30 to 50 would be sufficient to obtain significant differences. A power calculation also suggested that a similar sample size would be sufficient for validation of model-predicted differences in circulating IL-6 as a function of IL-6 SNP. Patient samples were obtained from a larger cohort of blunt trauma patients as part of an observational cohort described recently (10), and sample collection was stopped upon reaching 500 total patients as per the protocol approved by the University of Pittsburgh Institutional Review Board (see below). All data obtained from blunt trauma patients meeting the inclusion and exclusion criteria (see below) were used for this study. All inflammatory mediator data were used; any samples whose values were found to be outside of the detection range of the respective assays (see below) were re-assayed, and thus, outliers were not excluded.

Clinical trial design and summary

All human sampling was done after approval by the University of Pittsburgh Institutional Review Board and informed consent was obtained from each patient or next of kin as per Institutional Review Board regulations, and the overall characteristics of the large patient cohort that was used for the present studies have been detailed previously (10). Patients eligible for enrollment in the study were at least 18 years of age, admitted to the ICU after being resuscitated, and, per treating physician, were expected to live more than 24 hours. Reasons for ineligibility were isolated head injury, pregnant women, and penetrating trauma. Laboratory results and other basic demographic data were recorded in the database via direct interface with the electronic medical record. Three plasma samples, starting with the initial blood draw upon arrival, were assayed within the first 24 hours after trauma and then from days 1 to 7 after injury. To establish baseline levels of inflammatory mediators, a cohort of 11 healthy, non-injured human volunteers (7 males, 4 females; age, 41.9 ± 2.1 years) was recruited, and blood was withdrawn once. All blood samples were centrifuged, and plasma aliquots were stored in cryoprecipitate tubes at −80°C for subsequent analysis of inflammatory mediators. Relevant clinical data for all patients are provided in table S3.

Analysis of inflammatory mediators

IL-1β, IL-6, IL-10, IL-12p40, and TNF-α were measured using human-specific Luminex bead sets (BioSource-Invitrogen) using a Luminex 100 IS apparatus (MiraiBio). NO2/NO3 was measured by the nitrate reductase method using a commercially available kit (Cayman Chemical). All inflammatory mediator data are provided in the Supplementary Materials (table S3).

DNA isolation and genomic analyses

Genomic DNA was isolated from a subcohort of 98 blunt trauma patient’s blood (all survivors; see below) derived from the 147-patient validation cohort described above using the Genomic DNA Purification Kit (Qiagen) according to the manufacturer’s instructions. Isolated genomic DNA was then amplified for the IL-6 promoter SNP region [−174C 175 base pairs (bp), and −175G 175 bp] using a Sequence-Specific Primed polymerase chain reaction (PCR-SSP) mapped on a cytokine genotyping tray (One Lambda). After serial PCR amplification, the correct size of the amplified DNAs was confirmed by 2% agarose gel electrophoresis (Invitrogen). Expression-associated IL-6 genotypes [G/G (high IL-6), G/C (high IL-6), and C/C (low IL-6)] were determined on the basis of the presence or absence of a specific amplified DNA fragment according to the worksheet provided in the kit (One Lambda).

The 98-patient blunt trauma subcohort was categorized as having high–IL-6 SNP (G/G combined with G/C; n = 78) versus low–IL-6 SNP (C/C; n = 20) producers. The overall demographics of the high–IL-6 SNP patients were as follows: male = 44, female = 34; age 41.9 ± 1.8 years; ISS, 19.2 ± 1; ICU LOS, 5.3 ± 0.7; total LOS, 11.2 ± 0.9; and vent days, 2.3 ± 0.5. The mechanisms of injury were as follows: motor vehicle accidents (53), falls (14), or motorcycle accidents (11). Patients with low–IL-6 SNP demographics were as follows: male = 19, female = 1; age, 38.9 ± 3.15 years; ISS, 14.5 ± 2.2; ICU LOS, 4.7 ± 1; total LOS, 8.8 ± 1.8; and vent days, 2.15 ± 0.8. The mechanisms of injury were as follows: motor vehicle accidents (12), falls (5), and motorcycle accidents (3). All clinical, SNP, and IL-6 data for this patient cohort are provided in the Supplementary Materials (table S3).

Statistical analyses

All values are presented as means ± SEM. Determination of statistically significant differences was made using one-way analysis of variance (ANOVA) on ranks, with P < 0.05 considered significant. In addition, the output of the mechanistic mathematical model (see below) was tested to determine if there was a significant difference in predicted ICU LOS between different IL-6 production phenotypes. The relative IL-6 production level of an individual virtual patient is directly described by a single model parameter that is allowed to vary among patients. To determine if IL-6 production phenotype could be a possible predictor of patient outcome, the predicted ICU LOS of the top 10% IL-6 producers was compared against the bottom 10% IL-6 producers (as determined by the model parameter value). A determination of statistically significant difference was made using Student’s t test, with P < 0.05 considered significant.

Mathematical model development and analysis

A previously developed two-compartment mathematical model of porcine endotoxemia (19) was augmented to create a three-compartment model of human trauma. The previous model consisted of blood and lung epithelium compartments as an abstraction aimed at describing quantitatively lung dysfunction and hemodynamic collapse after injection of endotoxin in individual swine (19). The goal of the present study was to gain insights into the interplay between inflammation and organ dysfunction induced by blunt trauma in multiple tissues of individual human patients. This shift in subject, insult, and emphasized endpoints required the following tasks: (i) extension of the porcine endotoxemia model to reproduce blunt trauma in tissue; (ii) retuning of porcine inflammation parameters to generic human endotoxemia model; (iii) calibration of new trauma parameters based on data from individual human blunt trauma patients, with the ultimate goal of creating a large in silico population based on the data from this cohort of patients.

These steps are detailed in the Supplementary Materials. Schematically, the model structure is represented in Fig. 1A, and the modeled progression of inflammation is represented in Fig. 1B.

Reproduction of patient-specific data trajectories

Using the general model of human trauma described above, a subset of parameters (~10%) was varied between patients to allow the model to accurately reproduce the time course data at the individual level. This subset of parameters was chosen on the basis of maximum interpatient variability and parameter sensitivity, as well as biologically reasonable sources of variability (presumed to be related to interpatient genetic variability, especially in SNPs affecting the production of inflammatory cytokines). This list includes parameters related to the following: monocyte half-life, neutrophil half-life, endothelial cell (EC) half-life, effect of damage on EC inflammation, effect of TNF-α on EC inflammation, TNF-α production, IL-1β production, IL-6 production, IL-10 production, IL-12 production, NO production, and damage recovery rate. Clinical outputs and biomarkers were fit simultaneously for each individual using a genetic algorithm.

Generation and simulation of virtual patient population

A large population (n = 10,000) of virtual human trauma patients was generated by sampling from a uniform distribution of parameter values, as defined by the observed ranges in table S2. This virtual population was simulated under various levels of traumatic insult and key clinical outputs (total damage level and time to resolve damage) were recorded. Although not mapping directly to the current ISS classification (75), we simulated a low/intermediate ISS range (5 to 20), an intermediate/high ISS range (20 to 35), and a very severe ISS range (35 to 50). The data used for this process were based on patients with an actual ISS range of 4 to 50, and so our simulations in essence queried the behavior of the virtual patients outside of the primary data obtained from these patients. The ISS for each individual of the simulated population was assigned randomly based on these ranges. This large virtual population of human trauma patients gave a number sufficient to mine data properly and fully understand the mechanistic relationship of individual components of the system.

In silico randomized clinical trials

Model sensitivity analysis (see Supplementary Materials) suggested relationships between clinical outcomes and the propensity to produce IL-1β, IL-6, and TNF-α. Therefore, hypothetical treatments meant to inhibit these particular mediators selectively in the circulation as well as in tissues were tested in silico in the setting of severe blunt trauma. In the absence of a full, physiologically based pharmacokinetic model, the prototype for each treatment involves increasing the decay rate of the target mediator by a factor of 1000, essentially eliminating this mediator from the system at the time of therapy administration, mimicking the effects of treatment with a neutralizing antibody (21). To carry out this simulation, 10,000 virtual patients were simulated with severe traumatic injuries under four treatment scenarios: placebo, anti–TNF-α, anti–IL-6, and anti–IL-1β. All simulation projections were stored in a relational database. Monte Carlo techniques were used to sample realistic clinical trial-sized populations (n = 500 per therapy) repeatedly and observe their predicted clinical results over 5000 virtual trials. Although the Monte Carlo technique is not necessary to observe the mean of the virtual population, this technique has the advantage of being able to fully characterize the distribution of clinical trial response over many virtual trials. This negates the need for assumptions about normality and allows for an explicit definition of the margins of the distribution.

Extension of porcine endotoxemia model. Our previously published porcine endotoxemia model (19) provided the foundation of cytokine-mediated intra- and intercellular interactions, the role of inflammation in hemodynamic stability, and gas exchange across lung epithelium. Additionally, the overall health of the animal was summarized with a generic “damage” term that quantified the detrimental effects of hypoperfusion and hypotension. To extend the model to human trauma, we added a third compartment to represent all nonpulmonary epithelial tissue. Blunt trauma is applied to this “general tissue,” which releases damage-associated molecular patterns (DAMPs) in the tissue and activates tissue epithelial cells to initiate systemic inflammation (76, 77), much like the original insult of endotoxin to the blood in the two-compartment mathematical model (19). Additionally, the general “damage” term was updated to include the direct effects of this initial traumatic insult. The model is detailed in the Supplementary Materials. Schematically, the model structure is represented in Fig. 1A, and the modeled progression of inflammation is represented in Fig. 1B.

Parameter estimation to a generic human endotoxemia model. Because of the shift in species from swine to humans, it was necessary to retune parameters associated with both inflammatory mediator production and the sensitivity of the modeled immune cells to these individual mediators. Using a genetic algorithm, these parameters were estimated such that the model’s simulated trajectories of systemic inflammatory mediators matched observed biomarker trajectories collected from human endotoxemia experiments at the population level (78). Having incorporated endotoxemia data at the individual level, production parameters of the various inflammatory mediators included in the model were then varied to understand their acceptable ranges that corresponded to individual variability observed in the data. These ranges were then useful when fitting to individual human trauma patients as well as when producing virtual populations.

Calibration of new trauma parameters. A key model calibration step involved tuning the model with regard to sensitivity of epithelial cells to DAMPs. For constraint purposes, we assumed that model damage predictions should return to baseline quickly at low ISS regardless of the inflammatory response (as defined by IL-6 production rate). However, at higher ISS values, we assumed that these same small changes in inflammation response should have a greater impact on clinical outcome and, hence, the damage prediction. Figure S1 (A and B) demonstrates this phenomenon in a single patient following the final estimation of trauma related parameters. Table S1 shows the specific parameters used to obtain patient-specific outcomes. A formal parameter sensitivity analysis (model responsiveness analysis) was carried out to assess the validity of model parameters.

Generation and simulation of virtual patient population. With inflammatory mediator and clinical outcome trajectories matched across a wide range of patients and traumatic injury, the one-dimensional range of the varied parameters was gathered. For the immediate purpose of understanding the relative role of specific inflammatory mediators in trauma-induced inflammation and for comparison to the independently conducted principal components analysis, only the ranges of parameters related to inflammatory mediator production, as well as the parameter related to damage recovery rate (kd), were gathered (table S2).

Model responsiveness analysis. To define quantitatively the effect of the free parameters on model output, the sensitivity of these parameters was explored. One general method by which to quantify such influences is to consider the sensitivity of the model to the parameters as defined by the partial derivative of each model output (at prescribed times) with respect to each of the free parameters. Intuitively, a sensitive parameter will exhibit a large change in output when perturbed slightly, whereas insensitive parameters will show little change in output even when perturbed significantly. However, sensitivity is a very local property of the system; the result computed for any one parameter depends on the values assumed for all other free parameters. Additionally, the sensitivity of a given output with respect to a given parameter will also vary with time. In the present setting, our interest, however, was to characterize the general importance of these parameters over the range of values they can assume, without imposing assumptions about the specific values of the estimated parameters or limiting the analysis to a single time point.

In order to assess the relative influences of the parameters in a more global fashion, a “responsiveness” metric was defined over many sample virtual patients from the virtual population. For each patient, one-dimensional perturbations were applied to each parameter independently, and the model was simulated to compute the sensitivities of each parameter for each patient. For each parameter being studied in a patient, a vector of the relative changes observed for each (output, time) pair of interest was stored. A matrix was then constructed where each column corresponds to a parameter; specifically, each column was the concatenation of all of the vectors obtained at every patient when varying that parameter. Each row of the matrix therefore corresponded to the relative change observed at a unique (patient, analyte, time) tuple.

This matrix has a very natural interpretation: The norm of each column gives a scalar quantity that represents the global influence exerted by the parameter to which it corresponds. By simply sorting the columns according to their norms, a straightforward ranking of the importance of the parameters in question is obtained.


Fig. S1. Sensitivity of damage to IL-6 varies with ISS.

Fig. S2. Model-predicted survival trajectories.

Table S1. Parameters varied to reproduce patient-level trajectories.

Table S2. Observed parameter ranges from individual-level fits.

Table S3. All clinical, demographic, and inflammatory mediator data described in this publication (Excel file).

Appendix 1. Model equations and parameter values.


  1. Funding: This work was supported by NIH grants P50-GM-53789 to Y.V., T.R.B., and S.C.C., and R33-HL-089082 to Y.V., S.C.C., and G.N. Author contributions: D.B. performed patient data analysis and modeling, interpreted data, and contributed in writing the manuscript. R.A.N., K.A., and A.Z. performed patient data analysis. J.S. performed patient data analysis and modeling. D.A.B. and J.Y. performed the Luminex and NO2/NO3 measurements and collected samples for analysis. A.G. performed patient data analysis. A.A. performed patient data analysis and contributed to writing the manuscript. G.C. performed patient data analysis and modeling. G.N. participated in experimental design and contributed to writing the manuscript. R.Z. analyzed clinical data, performed statistical analysis, and contributed to writing the manuscript. S.C.C. participated in experimental design and contributed to writing the manuscript. T.R.B. participated in experimental design, interpreted data, and contributed to writing the manuscript. Y.V. participated in experimental design, interpreted data, contributed to writing the manuscript, and coordinated and supervised the overall research effort. Competing interests: Y.V. and T.R.B. are co-founders of and stakeholders in Immunetrics Inc., a company that builds and validates biologically accurate, predictive models. The other authors declare that they have no competing interests. Data and materials availability: All data, including model computations, are included in the paper or in the Supplementary Materials.
View Abstract

Stay Connected to Science Translational Medicine

Navigate This Article