Table 1 provides the values, units, and sources for all twenty-six model parameters. Because of the lumped nature of variables M, A, and ε, general units (“M-units”, “A-units”, and “ε-units”) were used to describe these populations, as treated previously (Reynolds et al., 2006; Arciero et al., 2010; Barber et al., 2013). To be consistent with experimental measures, the units for bacteria were given as 10⁶ bacteria. Several parameter estimation techniques were utilized to yield maximum parameter estimation efficiency and physiologically realistic results.

The unknown parameter space dimension was decreased by estimating several parameter values using previous studies. Following Reynolds et al. (Reynolds et al., 2006), k _A was chosen so that anti-inflammatory inhibition is at most 75% (Isler et al., 1999), and μ _A was chosen to be significantly smaller than typical cytokine decay rates (Bacon et al., 1973; Bocci, 1991; Fuchs et al., 1996; Huhn et al., 1997) because these cytokines have longer lasting downstream effects in the anti-inflammatory cascade. The quality of the parameter fit varied relatively little (<0.1%) when parameter T was allowed to vary. Choosing T too large made [fM ‒ T]₊ = 0 for all simulations while choosing T too small made [fM ‒ T]₊ ≈ fM, which violates the assumption that the system can tolerate a certain significant level of inflammation (Arciero et al., 2010). Setting T to one avoided both extremes without affecting fit quality.

Seven parameters remained that still had relatively high levels of uncertainty regarding their values. Those parameters (B _∞ , f, k _D , B _s , ν ₄, B _c , and k ₁) were estimated using a constrained least squares optimization on log-transformed data. The model fit and corresponding data is shown in Figure 3A. The fit was constrained based on qualitative observations regarding the general inflammation paradigm where typically one of three qualitative outcomes occur (Remick, 2007). In the first type of outcome, termed “healthy recovery,” individuals are affected by an infection, but their inflammatory response is quick and successfully eliminates the infection. In the second type of outcome, termed “aseptic death,” the inflammatory response eliminates many of the bacteria but overwhelms surrounding healthy tissue until death occurs. In the third type, termed “septic death,” large numbers of bacteria survive alongside a mounting inflammatory response and death results. Typical dynamics in other studies (Kumar at al., 2004; Clermont et al., 2004; Reynolds et al., 2006; Day et al., 2006) suggest that lower initial pathogen levels often produce healthy recovery outcomes, mid-level initial pathogen levels often produce aseptic outcomes, and high levels produce septic outcomes. To produce a model capable of similar behavior, the optimization process was constrained so that simulated rats injected with 128 × 10⁶ bacteria experienced healthy recovery outcomes, simulated rats injected with 248 × 10⁶ bacteria experienced aseptic outcomes, and simulated rats injected with 505 × 10⁶ bacteria experienced septic outcomes. Such an approach is supported by the following assumptions: 1) the relatively good fit to the data (R ² = 0.70); 2) relatively low levels of bacteria [O(10⁶)]in the 248 × 10⁶ initial bacterial load calibration experiments after 72 h (i.e., probably not septic); 3) persisting high levels of bacteria [O(10⁷)] in the 505 × 10⁶ initial bacterial load calibration experiments after 48 h (i.e., probably septic); and 4) high mortality rates for loads above 200 × 10⁶ bacteria (Namas et al., 2012). The large range of measured bacterial levels suggested that log-transforming the data was also reasonable. Importantly, the standardized residuals associated with the log-transformed measurements and simulated data (Figure 3B) passed ten out of ten normality tests (Öner and Deveci Kocakoç, 2017). Similar outcomes were not possible for data that were not log-transformed.

In the optimization procedure, the following squared residuals of the log-transformed data were minimized: ∑ i = 1 N m ( log ( y i ) − log ( y ^ i ) ) 2 (6) where y _i corresponds to all experimentally measured bacterial levels at the time points and initial bacterial loads listed in Supplementary Table S1, ŷ _i corresponds to the model-predicted bacterial levels at the same time points and initial bacterial loads, and N _m = 31 is the total number of measurements. Although a range of values for k ₂ and ν₁ are later used to assess their effects on system dynamics (see Outcome Dependence on Parameter Values), during the optimization, and by default unless otherwise stated, k ₂ and ν₁ were set to 0.6 1/M-units/h and 0.08 M-units/h, respectively, to be consistent with the values used in Reynolds et al. (2006). Parameter k ₁ was fit during the optimization procedure to yield its default value of k ₁ = 1.27/h; this value was used throughout the study except when assessing the impact of varying this parameter (Outcome Dependence on Parameter Values). Multiple standard optimization techniques were used during the process to ensure the finding of a reasonable optimum, including the Nelder-Mead simplex method, the steepest descent algorithm, and random restarts (Matlab code available via Github, see Supplementary Information). The resulting parameter estimates are given in Table 1.

Due to the high dimensionality of parameter space (26 parameters) and the relatively small amount of data (18 independent time points), only a small subset of the parameters should be optimized at any given time to avoid excessive computations and identifiability issues. Identifiability is discussed later in this section but briefly, the larger the number of fitted parameters, the more likely we will encounter an unidentifiable parameter set where multiple sets of parameter values can be used to produce the exact same optimal fit quality. Having no unique optimal parameter set can cause traditional optimization techniques to fail. Here, seven parameters were chosen to be fit based on their associated levels of uncertainty. It is possible that improved fits could be obtained by choosing a different parameter set [e.g., by fixing some of the fit parameters to other estimates found in the literature (Reynolds et al., 2006)]. In addition, while calibrating the model to a given data set will limit the general dynamics that a model may exhibit, it is possible that multiple optimal parameter sets exist whose dynamics significantly differ.

To explore 1) the suitability of the current model, 2) the set of fit parameters, and 3) the possible effects of varying other parameters, a local sensitivity analysis around the current best fit was performed. The relative sensitivities, s _ij , of the i ^th simulated measurement (ŷ _i) with respect to the j ^th parameter (p _j ) were estimated using the relative change in ŷ _i divided by the relative change in the parameter value: s ij = ( y ^ i ( p j + p j δ j ) − y ^ i ( p j ) y ^ i ( p j ) ) ( ( p j + p j δ j ) − p j p j ) = ( y ^ i ( p j + p j δ j ) − y ^ i ( p j ) δ j y ^ i ( p j ) ) (7) In Eq. 7, the magnitude of δ _j is chosen to be the square root of machine epsilon. For some parameters, differently signed δ _j ’s would alter the original simulation outcomes (healthy recovery, aseptic, septic). While such changes can happen in reality, this results in significant jumps in simulated measurement values and complicates the analysis. A complete analysis should consider behavior both when such jumps occur and when they do not and average the results. To include all parameters in the jump data analysis, δ _j would need to be altered so that all parameters experience such jumps. Since the jumps are typically of approximately the same size, an alternate strategy would need to be developed to compare these similarly sized jumps. Because of such complications, only sensitivities with no outcome behavior changes were considered here. At the same time, we note that larger calculated values of |s _ij| (calculated in the absence of outcome changes) correlate well with increased likelihood of these more significant outcome changes. To gauge the general sensitivity of all measurements to particular parameters and to understand how much a given parameter may affect outcomes, the root mean squared sensitivity (RMSS, S _i ) was calculated: S i = 1 N m ∑ i = 1 N m s ij 2 (8) where N _m = 31 gives the number of measurements made in the calibrating data set. This sensitivity value (S _i ) is provided in Supplementary Table S3 and provides an estimate for the relative importance of each parameter in the model fitting process with higher sensitivities corresponding to parameters that affect the fit and corresponding dynamics more. Better and quicker fits tend to result when parameters are more sensitive. For any given set of parameters, a sum of their RMSSs can give an estimate to how sensitive the measurements are to that specific parameter set.

Identifiability is another important consideration when assessing how parameter values affect system dynamics. As an example, model parameters k ₂ and s _l are not identifiable from the data because increasing k ₂ by 1% will produce the exact same change in the variable outputs as increasing s _l by 1%. This is because they only appear in the model as (k ₂ s _l ), never separately elsewhere. In such scenarios, it is impossible to identify the values of such parameters without additional outside information, as was the case in (Reynolds et al., 2006). Studying the identifiability of a large set of parameters with an associated nonlinear system can be difficult. While standard linear algebra methods can be used for linear systems (Cobelli and Distefano, 1980; Meshkat et al., 2009), nonlinear systems require more specialized and often more computationally and theoretically intensive approaches such as methods utilizing rational functions (Karlsson et al., 2012), methods focusing on local expansions or transformations (Meshkat et al., 2009), and methods centering on observability (Villaverde et al., 2016). Here we employed a local analysis, including linearization about the current model’s parameter values, to consider the issue of identifiability. In the local analysis, identifiability is investigated by considering the collinearity of the parameters, where the effects of changing one parameter can be reproduced by a linear combination of changes in the other parameters. In such cases, parameters can be redundant, with their variations capable of being reproduced by changing the other parameters. The idea of collinearity can be written in terms of the previously mentioned sensitivities: s _ik = ∑_j≠k α _j s _ij ∀i. A corresponding collinearity index can be defined as: CI_k = 1/σ_min where σ_min is the smallest singular value of the associated sensitivity matrix defined by components s _ij (Brun et al., 2001). Lower collinearity indices indicate more weakly correlated parameters that better describe the local parameter space and typically produce better overall fits and corresponding calibrated models. Such analysis identifies parameters sets that should not be used for fitting (high CI_k values) because of the possible presence of redundant parameters. The index can also be used to identify sets of parameters that may produce better fits.

The bacterial infections simulated in this study led to three possible outcomes: healthy recovery, aseptic death, or septic death (Reynolds et al., 2006). In this model, a healthy recovery outcome corresponds to a return to a steady state where pathogenic bacteria are fully eliminated (B = 0) and tissue damage levels return to baseline pre-infection levels (ε = 0). Aseptic death occurs when the bacteria are cleared (B = 0) but damage levels remain elevated (ε > 0). Septic death results when bacteria levels (B > 0) and damage (ε > 0) remain elevated at steady state.

In Figure 4, all three outcomes were predicted to occur given the same level of initial bacterial dose (B _source = 128 × 10⁶ bacteria) but for three different levels of pathogen growth rate, k ₁. For a low pathogen growth rate (k ₁ = 1.2/h), the pro-inflammatory response eliminated the bacteria, and healthy recovery (B/B _max, M/M _max, (A-A _healthy )/A _max, and ε/ε _max all less than 1%) was predicted to be restored by 37 h (cyan curve). For a slightly higher pathogen growth rate (k ₁ = 1.3/h), the pro-inflammatory response also eliminated the bacteria, but caused an elevated and sustained inflammatory response corresponding to aseptic death (blue curve) due to a forward-feedback loop of inflammation → damage → inflammation (Namas et al., 2012). For a high pathogen growth rate (k ₁ = 1.4/h), the pro-inflammatory response was incapable of eliminating the bacteria, so bacterial levels, damage, and the anti-inflammatory response all remained elevated at the steady state (corresponding to septic death, red curve).

Calculating local collinearities and sensitivities allowed for the identification of other parameter sets that are better (from a sensitivity and collinearity point of view) than the best fit identified in this paper. As an example, k _D , B _s , B _c , k ₁, k ₂, μ _l , and s _A have a sum of RMSSs of 137 and a collinearity index of 17 while the fit used in this paper has a sum of RMSSs of 62 and a collinearity of 55 (less sensitive and more collinear). Despite worse metrics, optimization of this alternative parameter set (originally suggested by our local sensitivity and collinearity analyses) did not significantly improve the overall parameter fit (R ² still approximately 0.70). This suggests that the parameters currently chosen for the model fitting are reasonable despite lower sensitivities and higher collinearities compared to other parameter sets.

The values of RMSS (Supplementary Table S3) showed that the model results are most sensitive to local immune response parameters (i.e., s _l , k ₂, μ _l and k ₃) and second-most sensitive to pathogen and dosing related parameters (i.e., k ₁, k _D , B _s , B _inf ). The other parameters affecting the inflammatory response and damage displayed lower sensitivities.

Prior studies utilizing reduced models of sepsis-induced inflammation were not docked to data, especially with regard to the level or trajectory of the damage variable that equated with death (Kumar et al., 2004; Day et al., 2006; Reynolds et al., 2006). To address this shortcoming, we compared the behavior of the damage variable in our model to mortality observations in an experimental model of Gram-negative bacterial sepsis.

In Figure 5A, the ability of the model to use damage levels to predict deaths occurring within 24 h (OMT = 24 h) and deaths occurring after 24 h (OMT >24 h) was quantified by comparing damage-based predictions of the mathematical model with the experimental observations in the validation data set. Using the mathematical model, the time of death (i.e., mortality time), t _c , for a simulated rat was assumed to be the model-predicted time at which the value of the damage variable, ε, reached a critical damage level (ε _crit ). If t _c ≤ 24 h, the rat was predicted to die within 24 h and to belong to the OMT = 24-h group. Similarly, if 24 < t _c ≤ 48, the model-predicted value for the OMT was 48 h and the rat was predicted to belong to the OMT >24-h group. The initial levels of bacteria given to each rat from the validation data set (27 total) were used as inputs to the calibrated model to yield model predictions for survival time (OMT = 24 h or OMT >24 h) for each rat in the study. Model accuracy was defined as the ratio of correct predictions to total number of predictions. Both the predictions and model accuracy are functions of ε _crit , as shown in Figure 5A. A relatively high level of accuracy (>80%) was found for a wide range of critical damage levels (49.2 ≤ ε _crit ≤ 123.4 ε-units). A maximum accuracy of 96% was attained for ε _crit values between 49.8 and 55 ε-units.

As an additional assessment tool of the predictive capability of the model, an area under the curve-receiver operating characteristics (AUC-ROC) curve was constructed (Figure 5B). In this figure, a positive outcome or identification corresponds to death occurring within 24 h. As such, ‘true positive’ corresponds to a model prediction agreeing with an experimental observation of a rat that dies within 24 h, while ‘false positive’ corresponds to a model prediction of a rat dying within 24 h when the experiment indicated that the rat did not die within 24 h. The recall/sensitivity/true positive rate = (number of true positives)/(number of correct predictions) and the false positive rate = (number of false positives)/(number of correct predictions). These rates are calculated for 0 ≤ ε _crit ≤ 200 ε-units to produce the AUC-ROC curve in Figure 5B. An area under this curve that is close to one indicates a model that can classify outcomes correctly. Here, the area under the curve was 0.98.

Since the experimental data in this study suggest significant differences between the OMT = 24-h and the OMT >24-h groups, simulations were conducted to determine whether the model could reveal additional differences that may exist between these two groups. In Figure 6, a simulation was run for each of the 27 different bacterial loads (B _source ) administered to the rats in the validation data set (Figure 1B) assuming that ε _crit = 53 ε-units (corresponding to the maximum accuracy determined in Figure 5A). Simulations that correspond to a predicted OMT = 24 h are shown in blue and those corresponding to an OMT >24 h are shown in red. As observed in Figure 6, the OMT 24-h group was predicted to have significantly higher bacteria, inflammatory, and damage levels. While there was a notable divide between model-predicted septic and aseptic cases in terms of bacterial levels (high vs near zero levels) after approximately 36 h, there was no particularly notable divide between the predicted OMT = 24-h (blue) and the septic members of the OMT >24-h (red) groups.

The model was also used to give possible explanations for the lack of statistical difference among the three OMT >24-h groups. As shown in Figure 1B and Supplementary Table S2, some rats with similar levels of initial bacterial loads exhibited significantly different outcomes. For instance, Rats 11 and 20 received similar initial bacterial loads, but Rat 11 died within 48 h while Rat 20 lived until 96 h. The model was used to investigate these large discrepancies by varying model parameters and evaluating their impact on system outcomes. For instance, these investigations showed that small changes in pathogen growth rate (k ₁) had a substantial effect on the predicted outcomes of the system. Figure 7 shows the model-predicted impact of decreasing this parameter slightly, given similar initial bacterial loads of 163 × 10⁶ bacteria in Rat 11 (red) and 172 × 10⁶ bacteria in Rat 20 (cyan). Specifically, k ₁ = 1.27/h for Rat 11, and k ₁ = 1.2/h for Rat 20. With these parameter values, the model predicted that Rat 11 died (i.e., its damage levels crossed the critical damage threshold of ε _crit = 53 ε-units) at approximately 34 h (corresponding to OMT = 48 h). Rat 20, however, was predicted to be healthy at 96 h (OMT = 96 h) with a low steady state value of bacteria and zero steady state value of damage (since pro-inflammatory levels never rise high enough to cause damage to the surrounding tissue). Similar changes in model predicted outcomes were observed when other parameters, such as those governing the immune response, were varied (Outcome Dependence on Parameter Values and Figure 8).

Figure 8 summarizes the dependence of stable outcomes (health is cyan, asepsis is blue, and sepsis is red) on the number of bacteria in the fibrin clot (B _source ), pathogen growth rate (k ₁ ), strength of the local immune response (k ₂), and pro-inflammatory activation rate (ν ₁). As k ₁ was increased or as k ₂ was decreased, the number of bacteria that the system could handle and still recover decreased (Figures 8A,B). The dashed lines in the figures indicate the three levels of B _source that were administered to the rats in the calibration data set. These multiple parameter bifurcation plots suggest that a healthy recovery outcome was not possible for many (high) levels of B _source . For lower levels of B _source , the model predicted that regions of healthy recovery become aseptic and then septic as the pathogen growth rate was increased (Figure 8A) or the strength of the local immune response was decreased (Figure 8B). As B _source was decreased, the width of the aseptic region decreased and completely disappeared when B _source = 6 × 10⁶ bacteria. This suggests that, regardless of the value of k ₁ or k ₂ , an aseptic outcome is not possible for a small enough initial bacterial load.

Parameter ν ₁ is the maximum activation rate of the pro-inflammatory response. As shown in Figure 8C, there was a very narrow region (199 × 10⁶ bacteria < B _source < 201 × 10⁶ bacteria) in which the model predicted that increasing ν ₁ causes a change in stable outcome from septic death to healthy recovery and then to aseptic death. This suggests that increasing the recruitment of pro-inflammatory cells and mediators first helps then overwhelms the system. For higher levels of B _source (B _source > 258 × 10⁶ bacteria), a moderate activation rate (ν ₁ = 0.08 M-units/h) divides the septic region (ν ₁ < 0.08 M-units/h) from the aseptic region (ν₁ > 0.08 M-units/h), suggesting that a highly responsive pro-inflammatory activation rate leads to aseptic outcomes.

The nearly equal split between septic and aseptic cases shown in Figure 8C for B _source > 258 × 10⁶ can be explained by the interactions between the bacteria and inflammatory response. The interaction term involving BM in the bacteria equation (Eq. 1) was initially very small since not enough M had been recruited. Thus, early dynamics for the bacteria equation were governed primarily by the remaining non-interaction terms (which are almost entirely unrelated to ν ₁). If B values increased above a (relatively low) threshold, the non-interaction terms became large (and positive) enough to outweigh any of the mounting contributions from the interaction term. This means that resulting high values of B depended on initial dynamics when M values were low and ν ₁ had only a minor impact. In addition, if B _source > 258 × 10⁶ bacteria, bacterial levels approached a high steady state level, regardless of the size of the interaction term. So, if ν ₁ < 0.08 M-units/h (right boundary of septic region in Figure 8C), the pro-inflammatory levels never grew large enough to outweigh the non-interaction terms in Eq. 1. If ν ₁ > 0.08 M-units/h, however, then the pro-inflammatory levels caused the interaction term in Eq. 1 eventually to outweigh the non-interaction terms, and the bacteria were eliminated. This elimination took place after the system attained high bacterial levels and was independent of B _source so long as B _source > 258 × 10⁶ bacteria.

Statistical analyses applied to the validation data set suggested no significant difference between the bacterial levels in rats with observed mortality times of 48, 72, and 96 h. Therefore, this study considered only two groups from the data set: rats that died before 24 h (OMT = 24 h) and rats that died after 24 h (OMT >24 h).

No clear separation among bacteria, inflammatory, or damage levels was predicted by the model for the OMT = 24-h and OMT >24-h groups. Interestingly, the model predicted very similar trajectories for the pro- and anti-inflammatory response. For example, there was a tight grouping of the pro- and anti-inflammatory response in the OMT = 24-h group. As a result, statistical models for this system would have difficulty predicting whether a septic or aseptic outcome should be expected, especially due to the uncertainty associated with detecting when the infection begins to take hold. These model predictions imply that pro-inflammatory, anti-inflammatory, or bacterial levels at early time points (e.g., 8 h) cannot be used to predict outcomes from a statistical perspective, thereby motivating the need for a combined statistical and mechanistic modeling approach.

Modeling helps to determine the relative importance (and validity) of proposed mechanisms and potential targets for successful interventions. For instance, the sensitivity analysis in this study showed that the local immune response is important. Obtaining patient specific parameters for the local immune response using well-designed measurements could improve parameter estimation, mitigate uncertainty in initial conditions, and enable more accurate patient-specific model-generated predictions regarding potential septic outcomes, corresponding dynamics, and proposed treatments. Such measurement planning and usage rely on well-chosen statistical techniques and could include both patient-specific measurements and/or general population measurements. In the case of the latter, maps of parameter spaces for general populations could be generated and then a few well-chosen measurements could be used to estimate the location of a specific patient in that space and their patient-specific corresponding parameter values. In addition, future models could make similar improvements in predictions by including mechanisms that involve specific (rather than general) cytokines (Chow et al., 2005) or those involving spatial organization of the system (Barber et al., 2013).

The validation data set showed that similar bacterial loads led (in some cases) to significantly different mortality times in rats. Just as previous models have demonstrated the impact of stochasticity on model outcomes (Cockrell and An, 2017; Mavroudis et al., 2019), the current model demonstrated that experimental or clinical variability in sepsis outcomes can be explained by very small differences in parameters governing bacteria or the host response. Additionally, the different outcomes in two rats with nearly identical bacterial loads could result from small differences in initial conditions due to underlying stressors at the beginning of the experiment. As a result, in order to enable patient-specific care, different parameter sets and/or initial conditions must be used for different individuals to capture variability among patients.

Several sepsis modeling studies have been described previously. Kumar et al. (2004) used a model similar to the model developed in the current study to predict that sepsis may have multiple negative outcomes (e.g., septic and aseptic death) that may require different treatment approaches. However, their study did not explicitly incorporate any experimental data. Yamanaka et al. (2019) created a model founded on clinical data focused solely on septic shock, an important subset of sepsis. McDaniel et al. (2019) introduced a “whole-body” sepsis model (via the BioGears Physiology Engine) to be used as an in silico septic simulator. While these and several other sepsis models (An, 2001; An, 2004; Chow et al., 2005; Day et al., 2006; Prince et al., 2006; Vodovotz et al., 2006; Zuev et al., 2006; Song et al., 2012; Shi et al., 2015; Cockrell and An, 2017) have offered many useful insights, the present study emphasized the acute inflammatory dynamics that take place during early sepsis development using a relatively simple and novel modeling framework that can be used to identify the mechanisms that underly vastly different outcomes despite nearly identical initial conditions. Notably, our results using experimental sepsis induced in genetically identical rats support those of Cockrell and An (2017), who used an agent-based model of human sepsis to suggest that mortality could occur under diverse conditions and influences which in turn would defy stratification based on inflammation biomarkers.

As is necessary in theoretical modeling, the model presented in this study applies simplifying assumptions to make predictions about the average health status in a system. For example, several immune system mediators are grouped together into two general populations for the pro- and anti-inflammatory responses. Also, the model considered only virulent bacteria; non-virulent (e.g., probiotic) bacteria (Arciero et al., 2010) were not included. Although additional model components could improve predictions, the simplicity of this model allows for very useful analysis of the impact of parameters and interactions on the health of a rat. The simplicity of the model, however, does limit the ability to capture more complicated underlying dynamics. For instance, a single variable, ϵ_crit, is used to predict outcomes based on system damage when outcomes likely depend on additional factors. Similarly, differences between species are subsumed into a single parameter (B _s ) when converting bacteria-associated quantities from the human system developed by Reynolds et al. (2006) to the rat system considered in this study.

The model was calibrated using both quantitative data (bacterial levels at multiple time points) and qualitative observations (coexistence of healthy recovery/septic/aseptic steady states for some parameter values). The data used in this study did not include measurements of cytokine levels, and thus using data sets with both bacterial and cytokine levels could help to improve model calibration and design. In addition, the validation data only included time points of 24, 48, 72, and 96 h. There is probably a more continuous change in mortality outcomes, and thus a more accurate threshold for death could be obtained given additional time point data. Finally, the model has been calibrated to peritoneal injections that cause sepsis in rats. More experimental data would allow for more general insights.

With complicated processes such as inflammation where hundreds of molecules, cells, and other factors play roles, mathematical models are essential for providing a mechanistic understanding of the system, as in the current sepsis study. This is especially true when very small differences in initial conditions or parameter values in complex systems can have a major impact on outcome, as illustrated in this study. Thus, model reduction is needed to facilitate analysis and interpretation. The parsimonious model of sepsis presented here, after calibration, reproduced experimental results, identified an inherent level of uncertainty associated with experimental data and associated predictions, predicted trends as bacterial load, pathogen growth rate, strength of the local immune response, and activation rate of the pro-inflammatory response were varied, and provided a simplifying paradigm that can be used to understand the life vs death outcomes for septic individuals. As recent studies in blunt trauma (Abboud et al., 2016a), traumatic brain injury (Abboud et al., 2016b), and pediatric acute liver failure (Zamora et al., 2016; Zamora et al., 2019) have shown, a dichotomy in patient outcomes was observed as soon as measurements were taken, suggesting that inflammation, regardless if it results from trauma, disease, or sepsis, exhibits the same features illustrated in this work. Using mathematical modeling provides an understanding of possible mechanisms that could explain such dichotomies in outcomes, in contrast to methods solely based on statistical assessments of clusters or outcomes.