Predicting Experimental Sepsis Survival with a Mathematical Model of Acute Inflammation

Sepsis is characterized by an overactive, dysregulated inflammatory response that drives organ dysfunction and often results in death. Mathematical modeling has emerged as an essential tool for understanding the underlying complex biological processes. A system of four ordinary differential equations (ODEs) was developed to simulate the dynamics of bacteria, the pro- and anti-inflammatory responses, and tissue damage (whose molecular correlate is damage-associated molecular pattern [DAMP] molecules and which integrates inputs from the other variables, feeds back to drive further inflammation, and serves as a proxy for whole-organism health status). The ODE model was calibrated to experimental data from E. coli infection in genetically identical rats and was validated with mortality data for these animals. The model demonstrated recovery, aseptic death, or septic death outcomes for a simulated infection while varying the initial inoculum, pathogen growth rate, strength of the local immune response, and activation of the pro-inflammatory response in the system. In general, more septic outcomes were encountered when the initial inoculum of bacteria was increased, the pathogen growth rate was increased, or the host immune response was decreased. The model demonstrated that small changes in parameter values, such as those governing the pathogen or the immune response, could explain the experimentally observed variability in mortality rates among septic rats. A local sensitivity analysis was conducted to understand the magnitude of such parameter effects on system dynamics. Despite successful predictions of mortality, simulated trajectories of bacteria, inflammatory responses, and damage were closely clustered during the initial stages of infection, suggesting that uncertainty in initial conditions could lead to difficulty in predicting outcomes of sepsis by using inflammation biomarker levels.


INTRODUCTION
Sepsis is defined as life-threatening organ dysfunction caused by an overwhelming immune response to an infection (Singer et al., 2016). Sepsis and septic shock often lead to organ failure and are leading causes of morbidity and mortality (Singer et al., 2016;Lelubre and Vincent, 2018). In the course of a typical infection, pathogens (e.g., Gram-negative bacteria) infect the host, triggering a pro-inflammatory innate immune response targeted at eliminating the pathogen from the system. The body also mounts an anti-inflammatory response that works to maintain homeostasis and prevent overwhelming inflammation to the system. The acute inflammatory response in sepsis becomes detrimental when it can no longer be contained at the locus of infection (where innate immune mechanisms can help eliminate bacteria) and becomes systemic, leading to collateral damage/ dysfunction in the surrounding healthy tissue. Thus, a major goal in sepsis research is to determine why the host mounts an overwhelming inflammatory response to infections which leads to sepsis in some cases but not in others (Horiguchi et al., 2018;Grondman et al., 2020).
Acute inflammatory responses to infection are further complicated by additional variations among sepsis patients (Chen et al., 2011), the complex interplay between pro-and anti-inflammatory mediators and the innate and adaptive immune responses (Delano and Ward, 2016), and a triphasic distribution of patient deaths that occur days, weeks, and years after initial infection (Delano and Ward, 2016). These variations have highlighted the need for patient-specific treatments and novel approaches to sepsis drug design and clinical trials Buchman et al., 2016;Day et al., 2018).
Recent decades have brought about improved sepsis treatment practices, including the rapid administration of antibiotics, infection source control, appropriate choice of fluid (crystalloids) for fluid resuscitation, and administration of vasopressors (norepinephrine) (Schorr and Dellinger, 2014). These evolutionary improvements in patient care have reduced in-hospital mortality but have shifted the mortality distribution so that most sepsis-related deaths occur months after treatment (Delano and Ward, 2016). While improved care often allows patients to overcome initial septic episodes, it has been speculated that underlying physiologic/biochemical aberrations combined with sepsis-initiated immune dysfunction place patients at longterm risk for sepsis mortality (Yende et al., 2008). Secondary infections occurring during late stages (>15 days) were shown to correlate with higher levels of opportunistic fungi and bacteria compared to earlier secondary infections (<6 days) (Otto et al., 2011), and impairments in immune function and cytokine secretion were observed in individuals that survived a septic episode (Arens et al., 2016). Despite evidence of late-stage complications or mortality from sepsis (Nedeva, 2021), mortality is still common within a week of sepsis onset; interventions in those cases rely on early identification of septic trajectories, which is often difficult.
The current study implemented a mechanistic model of the immune response to an Escherichia coli (E. coli) infection and aimed to determine whether early data on both the pathogen and the host response are sufficient to predict a health or disease outcome. A parsimonious system of four ordinary differential equations (ODEs) was used to track changes in bacteria, tissue damage, a pro-inflammatory response, and an anti-inflammatory response following the administration of escalating inoculums of E. coli encapsulated in a fibrin clot (an experimental model of Gram-negative peritonitis (Ahrenholz and Simmons, 1980)) in rats. Ultimately, the mathematical model defined in this study was used to evaluate and predict the time dynamics of the bacterial infection, which were observed to vary depending on the initial bacterial dose. The model demonstrated that small changes in the pathogen growth rate or characteristics of the immune response could lead to significantly different mortality times despite identical initial bacterial loads.

Experimental Method
Preparation and Dose Estimation of E. coli-Impregnated Fibrin Clot All animal studies were carried out following approval by the University of Pittsburgh Institutional Animal Use and Care Committee (IACUC approval #0807947) and complied with the NIH Guide for the Care and Use of Laboratory Animals. In this study, experimental data were obtained from rats with peritonitis induced by an E. coli-impregnated fibrin clot, similar to the methods established by Ahrenholz and Simmons (1980) and modified by Namas et al. (2012). In brief, the animals were subjected to a varying dose of E. coli (strain ATCC 25922; American Type Culture Collection, Manassas, VA, United States) inoculum in a fibrin clot introduced into the peritoneum via laparotomy. The E. coli colonies were grown to a specified optical density using a spectrophotometer (DU 530 UV/VIS; Beckman Coulter, Brea, CA, United States) equivalent to a concentration ranging between 1 × 10 8 to 5 × 10 8 colonyforming units [CFUs]/clot on the day of bacterial fibrin clot implantation. After the addition of fibrinogen (1%) and Thrombin (15u), the clot was placed in the peritoneum via laparotomy, as described previously (Namas et al., 2012). Since quantification of bacteria at the time of implantation of the fibrin clot via optical density is imprecise, bacteria were quantified after each implantation by limiting dilution plating to obtain a count of implanted E. coli in a given rat (Namas et al., 2012).
We determined that rats that received a fibrin clot inoculum of 1 × 10 8 -2.0 × 10 8 CFUs/clot had a mortality rate of ∼40-45% after 48 h of clot implantation while rats that received inocula higher than 2.0 × 10 8 CFUs/clot had a mortality rate of 80% during the first 24 h of implantation (data not shown). All rats in these experiments exhibited signs of septicemia in the form of lethargy, hypothermia, reluctance to feed, tachycardia, and tachypnea after 24 h of E. coli fibrin clot impregnation into the peritoneum. After closing the abdomen, a topical anesthetic was applied over the surgical wound, and the rats were returned to their cages and allowed food and water ad libitum.

Model Calibration and Validation
The experimental studies detailed above generated one experimental data set that was used to calibrate the model and consisted of measurements of bacterial levels at different timepoints in 31 rats injected with increasing doses of bacteria (1.28, 2.48 and 5.05 × 10 8 bacteria with standard error of the mean of 1.5 × 10 7 , 2.8 × 10 7 , and 1.2 × 10 8 , respectively). A second (separate) experimental data set was used to validate the model and consisted of time to death (mortality; 24, 48, 72, or 96 h) for 27 rats injected with varying doses of bacteria (1-5 × 10 8 bacteria). We detail the process of model calibration and validation below.
Bacterial levels in the peritoneal cavity measured in individual septic rats, each euthanized at a distinct time point from 24 to 96 h after being impregnated with ≈1 × 10 8 -5 × 10 8 colonyforming units of E. coli, were used to calibrate the model ("calibration data set", see Figure 1A; Supplementary Table S1). Observed times to mortality for rats injected with ≈1.4 × 10 8 -4.8 × 10 8 colonyforming units of E. coli were used to evaluate the predictive capabilities of the ODE model ("validation data set," see Figure 1B; Supplementary Table S2). These two data sets were derived from two different rat populations.
The calibration data set ( Figure 1A; Supplementary Table S1) was produced using several experiments involving multiple rats in each experiment and three different bacterial levels with mean values of 1.28 × 10 8 , 2.48 × 10 8 , and 5.05 × 10 8 bacteria. A peritoneal lavage was used to assess the number of bacteria in the peritoneal cavity in each individual animal euthanized at one of the time points indicated in Supplementary Table S1, as detailed in (Namas et al., 2012). This produced an estimate for the expected bacterial levels over time in a rat subject to one of the targeted loads, as seen in Figure 1A. The time points were chosen to be closer together for rats injected with higher bacterial loads, since these animals were more likely to die earlier than animals receiving a lower bacterial inoculum.
In the validation data set ( Figure 1B; Supplementary Table  S2), the rats were monitored continuously until 96 h, at which time surviving rats were euthanized. If a rat was observed to be dead at a 24-h checkpoint (24, 48, 72, or 96 h), that time point was noted as the "observed mortality time" (OMT) for that rat and the rat was removed from the study. Rats labeled with an OMT 24, 48, or 72 h either died within the 24 h before their OMT or were euthanized at their OMT (in accordance with IACUC and NIH guidelines) due to being severely unhealthy. Rats labeled with a 96 OMT either died between 72 and 96 h, were severely unhealthy at 96 h, or were relatively healthy at 96 h. Given these distinctions, it is possible that some of the 96 h OMT rats would have survived longer. The calibration data set is depicted using closed circles (C) for each measured bacterial level. Bacterial measurements were obtained for the following three initial bacterial loads: 1.28 × 10 8 bacteria (cyan), 2.48 × 10 8 bacteria (blue), and 5.05 × 10 8 bacteria (red). Lines connect the geometric means of the measurements at each time point. (B) The validation data set includes the observed mortality time for twentyseven rats injected with varying levels of bacteria via an E. coli-impregnated fibrin clot (each measurement is represented with a C). The bottoms and tops of the boxes in the box-and-whiskers plot indicate the edges of the first and third quartiles, respectively, and the whiskers extend to the smallest and largest values in each group that are less than 2.7 within-group standard deviations from the mean. The data points for Rat 11 and Rat 20 (Comparison of model simulations and experimental data suggest a high mortality) are denoted by a red and cyan C, respectively.
Frontiers in Systems Biology | www.frontiersin.org November 2021 | Volume 1 | Article 755913 Figure 1B depicts the mortality data set (closed circles, •) collected for 27 rats with the observed mortality time on the x-axis and the initial bacterial load on the y-axis. A Bartlett test revealed that the variances of the four groups (OMT 24,48,72,96 h) were not significantly different (p-value 0.689) while ten out of ten normality tests (Öner and Deveci Kocakoç, 2017) on the deviations in the data from their respective group norms were passed, suggesting that the data in each group are approximately normally distributed. Analysis of variance (ANOVA), which assumes normally distributed data with equal variances across groups, was then performed to show that the initial bacteria levels for the three OMT 48-, 72-, and 96-h groups with OMT >24 h were not significantly different from each other (p-value 0.251). More specific delineation of the health and disease status of the rats in the 96-h group may have led to more variance between the OMT 48-, 72-, and 96-h groups, but those data were not available at the time of this study. ANOVA on all four groups and use of Tukey's procedure with 95% confidence intervals suggested that the OMT 24-h group had significantly different levels of initial bacterial load compared to the OMT 48-, 72-, and 96-h groups (p-value < 0.0051). The analysis was performed with Matlab (normalitytest (Öner and Deveci Kocakoç, 2017) and Matlab's vartestn, anova1, and multcompare; Matlab code available via Github, see Supplementary Information), and the resulting conclusions regarding the differences between the groups were consistent with the relationships suggested by the box and whiskers plot shown in Figure 1B. The analysis suggested that there were effectively only two different experimental groups: rats with OMT 24 h and rats with OMT >24 h. These two categories were therefore used when assessing the predictive capabilities of the model.

Model Variables and Interactions
The mathematical model developed in this sepsis study was based on previous models of the immune response Reynolds et al., 2006;Arciero et al., 2010;Arciero et al., 2013a;Barber et al., 2013). It was used to predict the dependence of health or disease outcomes on pathogen properties and the immune response. Following the administration of E. coli into the peritoneum, the model was used to predict changes in the number of bacteria in the peritoneal cavity of the animal (B), the level of tissue damage (ε), a pro-inflammatory response (M), and an anti-inflammatory response (A). Tissue damage serves as a proxy for whole-organism health status. In vivo, the molecular correlates to tissue damage are damage-associated molecular pattern [DAMP] molecules such as high-mobility group box 1 (HMGB1) (Wang et al., 2004;Lotze and Tracey, 2005). Due to limited calibrating data and a desire for model simplicity, the numerous cells and cytokines that generate the pro-and antiinflammatory responses were grouped into two general populations, as in previous models Arciero et al., 2010;Arciero et al., 2013b). Specifically, tissuespecific effects of bacterial infection were assumed to be driven by both resident tissue macrophages and inflammatory cells such as neutrophils and macrophages that infiltrate the tissue from the blood (Cavaillon and Annane, 2006;Remick, 2007;Qiu et al., 2019;Kumar, 2020). Nominally pro-inflammatory cytokines (e.g., TNF, IL-1b, IFN-γ, IL-17A) and anti-inflammatory cytokines (e.g., TGF-β1, IL-10, IL-4, IL-13) were not measured in these experiments, and therefore these cytokines (along with a multitude of chemokines and DAMPs) were grouped into general populations whose effects were assumed independent of inflammatory mediator type. Figure 2 provides a schematic of the assumed interactions (Remick, 2007) among model populations (schematic adapted from ). A time-dependent dosing function, D(t), defines the release of bacteria from the clot into the surrounding tissue. In response, pro-inflammatory cells activate and begin to destroy the bacteria while also causing collateral tissue damage and promoting self-recruitment. Both pro-inflammatory cells and tissue damage trigger an antiinflammatory response, which in turn inhibits the growth of the pro-inflammatory response and damage levels. These interactions are labeled in Figure 2 with arrow heads indicating upregulation and blunted ends indicating downregulation.
Although it is nearly impossible to state that no other model could account for the variables and interactions depicted in Figure 2 (Ganusov, 2016), the approach utilized in this study adapted previous models Arciero et al., 2010) of bacterial infection that were successful at reproducing multiple pathophysiological behaviors. The current model was based on clear assumptions, including ones that stemmed from the experimental technique for bacterial administration into the rats. More specifically, Ahrenholz and Simmons (Ahrenholz and Simmons, 1980) originally used the experimental method implemented in this study to consider the theory that increased fibrin levels decrease the severity of infections during E. coli peritonitis. Namas et al. (2012) used the same procedure to explore the effectiveness of hemadsorption techniques on treating septic peritonitis. They performed several initial experiments in which they injected rats with fibrin clots containing varying levels of bacteria. The rats usually died or resolved in a relatively short amount of time, so only the acute immune response was assumed to play a significant role in infection dynamics. Additionally, the experimental technique for bacterial administration introduced bacteria more gradually into the system than instantaneous bacterial injection into the blood. Both observations played a significant role in the modeling choices employed here.

Model Equations
The rate of change of the number of bacteria in an animal is governed by Eq. 1: The first term, D(t), defines the release of bacteria from the fibrin clot as an exponentially decaying dosing function given by D(t) D 1 exp(-k D t). The second term corresponds to logistic growth of bacteria with growth rate k 1 . The third term corresponds to a non-specific, local, innate immune response that is assumed to eliminate a small amount of pathogen without activating a full systemic inflammatory response . The details of the parameter value derivation for this term are given in . The final term defines the elimination of bacteria by the systemic pro-inflammatory response, which is inhibited by systemic anti-inflammation. As a first-order approximation, all bacteria initially in the clot were assumed to empty into the surrounding tissue without reproducing or dying, which is equivalent to where B source is the initial amount of bacteria in the fibrin clot (B source 128 × 10 6 , 248 × 10 6 , or 505 × 10 6 bacteria in the calibration data set). To satisfy this condition for the conservation of clot bacteria, D 1 must equal k D B source . B was also scaled by a factor B s (Parameter Values) to account for the difference in bacterial levels between the rat experiments described in this study and the bacteria levels considered in the previous (human) models . This factor also converts number of bacteria to concentration of bacteria (bacteria/cm 3 ) and allows previous parameter values Arciero et al., 2010) to be used here.
We note that multiple dosing functions were considered before choosing the exponential decay model. For example, an additional variable tracking the bacterial levels inside the clot was introduced and allowed to experience growth and/or transfer between the clot and surrounding tissue. Despite providing three additional degrees of freedom, this alternate approach did not yield model fits that were significantly better than using a simple exponential function to define the rate at which bacteria emptied into the surrounding tissue. We also explored making parameter k d a saturating sigmoidal function of B source , since experimental observations suggested that a clot with a small number of bacteria had a slower release rate than a clot containing a large inoculum of bacteria. This assumption, however, also did not lead to improved model fits. Thus, an exponential function with constant rate of decay, k d , was implemented in the model and yielded reasonable fits to the data.
Equation 2 defines the rate of change of the systemic proinflammatory response (M): dM dt The pro-inflammatory response is increased by pro-inflammatory cells, bacteria, and damage and is inhibited by the systemic anti-inflammatory response. Natural decay of the response occurs at rate µ M . Recruitment of inactive pro-inflammatory cells occurs at a rate proportional to ] 1 , as given in Reynolds et al., 2006). The possible effect of LPS inflicting toxicity (without any bacteria present) has been subsumed into the parameter values for the onset and persistence of inflammation as a consequence of bacterial infection. T helper one cells (Th1) and T helper 17 cells (Th17) are also assumed to be included in model population M.
The rate of change of the anti-inflammatory response (A) is defined in Eq. 3: The first term accounts for the nonzero background level (s A ) of anti-inflammatory mediators that resides in the body. In the second term, the pro-inflammatory response and damage activate the anti-inflammatory response. The anti-inflammatory response also self-regulates, which is modeled using the (1 + k A A) term in the denominator. Natural decay of the response occurs at rate µ A . The final model equation (Eq. 4) gives the rate of change of tissue damage levels (ε): In the first term, damage is repaired at rate 1/τ. In the second term, similar to (Arciero et al., 2010), it was assumed that the system can tolerate a certain level of inflammation before proinflammatory cells begin causing damage. This assumption is represented using a threshold function: T gives the threshold above which immune cells begin to cause damage in the nearby tissue, and f gives the rate at which the cells cause damage.

Initial Conditions
At the onset of each simulation, the rats were assumed to be healthy, corresponding to initial conditions of B(0) 0, M(0) 0, A(0) s A /µ A , and ε(0) 0. Since animals exhibit background levels of anti-inflammatory mediators (e.g., transforming growth factor-β1 (Morikawa et al., 2016)), nonzero initial background levels of anti-inflammatory response mediators were assumed. Parameter B source corresponds to the number of bacteria inserted into the clot for each experiment and controls the magnitude of the dosing function, D(t). Thus, B source was changed in each simulation (i.e., for each rat) to correspond to the experiment being modeled. For a sufficiently high B source , it was observed that the clots were highly saturated to the point that bacteria left the clot nearly instantaneously. To model this, B(0) was set to B source -B c , where B c is a fitted parameter thought of as the "clot capacity." This assumes the number of bacteria that exceed B c instantaneously enter the surrounding tissue. With only B c bacteria assumed to remain in the clot after the initial injection, the dosing function is adjusted to depend only on those bacteria (B c ) remaining in the clot: Parameter Values 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 Arciero et al., 2010;Barber et al., 2013). To be consistent with experimental measures, the units for bacteria were given as 10 6 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 antiinflammatory 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 , ] 4 , B c , and k 1 ) 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 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 6 bacteria experienced healthy recovery outcomes, simulated rats injected with 248 × 10 6 bacteria experienced aseptic outcomes, and simulated rats injected with 505 × 10 6 bacteria experienced septic outcomes. Such an approach is supported by the following assumptions: 1) the relatively good fit to the data (R 2 0.70); 2) relatively low levels of bacteria [O(10 6 )]in the 248 × 10 6 initial bacterial load calibration experiments after 72 h (i.e., probably not septic); 3) persisting high levels of bacteria [O(10 7 )] in the 505 × 10 6 initial bacterial load calibration experiments after 48 h (i.e., probably septic); and 4) high mortality rates for loads above 200 × 10 6 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 logtransformed.
In the optimization procedure, the following squared residuals of the log-transformed data were minimized: 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 modelpredicted 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 2 and ν 1 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 2 and ν 1 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 1 was fit during the optimization procedure to yield its default value of k 1 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.

Sensitivity Calculation
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 FIGURE 3 | Validity of fitting the mathematical model to log-transformed data. (A) Calibrating data set (C) plotted with model-predicted bacterial dynamics over the experimental times considered (solid curves). The following three initial bacterial loads were provided in the data set: 128 × 10 6 bacteria (cyan), 248 × 10 6 bacteria (blue), and 505 × 10 6 bacteria (red). (B) Plot of the standardized residuals for the given fit to log-transformed data. These plots demonstrate the validity of using logtransformation and fitting via least squares.
Frontiers in Systems Biology | www.frontiersin.org November 2021 | Volume 1 | Article 755913 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 ].
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: 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: 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 2 and s l are not identifiable from the data because increasing k 2 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 2 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 . 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.

Bacterial Infection Model Dynamics Reproduce Key Sepsis Outcomes
The bacterial infections simulated in this study led to three possible outcomes: healthy recovery, aseptic death, or septic death . 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.
Frontiers in Systems Biology | www.frontiersin.org November 2021 | Volume 1 | Article 755913 In Figure 4, all three outcomes were predicted to occur given the same level of initial bacterial dose (B source 128 × 10 6 bacteria) but for three different levels of pathogen growth rate, k 1 . For a low pathogen growth rate (k 1 1.2/h), the proinflammatory 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 1.3/h), the proinflammatory 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 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).

Sensitivity Analysis Supports Choice of Initial Model Parameters
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 1 , k 2 , μ 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 2 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 2 , μ l and k 3 ) and second-most sensitive to pathogen and dosing related parameters (i.e., k 1 , k D , B s , B inf ). The other parameters affecting the inflammatory response and damage displayed lower sensitivities.

Comparison of Model Simulations and Experimental Data Suggest a High Mortality Prediction Capacity
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 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, FIGURE 4 | Predicted dynamics for healthy recovery, septic, and aseptic outcomes as pathogen growth rate is varied. Model predictions of bacteria (A), pro-inflammatory (B), antiinflammatory (C), and damage (D) levels as pathogen growth rate (k 1 ) was varied. Three different outcomes were predicted for an initial bacterial infection of B source 128 × 10 6 bacteria: healthy recovery (k 1 1.2/h, cyan), asepsis (k 1 1.3/h, blue), and sepsis (k 1 1.4/h, red). November 2021 | Volume 1 | Article 755913 10 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 1 ) 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 6 bacteria in Rat 11 (red) and 172 × 10 6 bacteria in Rat 20 (cyan). Specifically, k 1 1.27/h for Rat 11, and k 1 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 1 ), strength of the local immune response (k 2 ), and proinflammatory activation rate (] 1 ). As k 1 was increased or as k 2 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 Blue curves correspond to rats that had a predicted observed mortality time of 24 h, and red curves correspond to rats that were predicted to survive for more than 24 h. The clustering of many of the trajectories suggests that small uncertainties in the initial state of the system could lead to significantly different systemic predictions.

Outcome Dependence on Parameter Values
Frontiers in Systems Biology | www.frontiersin.org November 2021 | Volume 1 | Article 755913 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 6 bacteria. This suggests that, regardless of the value of k 1 or k 2 , an aseptic outcome is not possible for a small enough initial bacterial load. Parameter ] 1 is the maximum activation rate of the proinflammatory response. As shown in Figure 8C, there was a very narrow region (199 × 10 6 bacteria < B source < 201 × 10 6 bacteria) in which the model predicted that increasing ] 1 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 6 bacteria), a moderate activation rate (] 1 0.08 M-units/h) divides the septic region (] 1 < 0.08 M-units/h) from the aseptic region (ν 1 > 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 6 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 ] 1 ). 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 ] 1 had only a minor impact. In addition, if B source > 258 × 10 6 bacteria, bacterial levels approached a high steady state level, regardless of the size of the interaction term. So, if ] 1 < 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 ] 1 > 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 6 bacteria.

DISCUSSION
In this study, a mathematical model of the immune response to a bacterial infection was developed from previous inflammation models Arciero et al., 2010;Barber et al., 2013) and applied to an experimental study in which variable doses of E. coli were administered to rats. The model was calibrated to bacterial levels in the rats at multiple time points following implantation of the E. coli fibrin clot. The model successfully predicted outcomes of observed mortality times in 27 rats given several different initial bacterial loads. The differences in model-predicted levels of bacteria, proinflammatory response, anti-inflammatory response, and damage for rats with OMT 24 h and rats with OMT >24 h, however, were not as large as expected. Thus, the model accurately predicts observed outcomes when initial conditions are well-known, but predicting final outcomes based on data at early time points is challenging when large uncertainty surrounds the appropriate initial conditions to use for model or statisticsbased predictions. The model also showed that variability in experimental outcomes can result from variability in the pathogen growth rate, strength of the local immune response, or maximum activation rate of the pro-inflammatory response, since small changes in these parameter values generated large changes in system outcomes.
An important goal for mathematical models of disease is to yield accurate patient-specific outcomes. Ideally, given early time course data on bacterial levels and cytokine levels in a patient, a model could be used to predict whether that patient will follow a healthy recovery (resolving), septic, or aseptic outcome. Then, depending on the forecasted outcome, appropriate countermeasures could be prescribed, allowing for an efficient and accurate use of resources. In many studies, pure statistics are used to predict the likely status or outcome of a given patient. Applying pure statistics to the validation data set in this study (e.g., using a predictor that predicts death within 24 h if the initial bacterial load exceeds a certain amount) would produce similar successful predictions but would not be able to predict septic vs aseptic outcomes, explain failed predictions, or reveal the underlying system dynamics. While more complex statistics using more data per patient can help to improve the accuracy of statistical methods, ultimately such methods still do not identify the mechanisms leading to observed outcomes as mechanistic mathematical models do. An optimal theoretical approach for understanding disease should include a combination of mathematical modeling and statistical methods (Albers et al., 2018;Baker et al., 2018). The study presented here provides an example of such a combined mathematical and statistical approach.

Mathematical and Statistical Predictions
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 antiinflammatory 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 .
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.

Comparison to Prior Mathematical Modeling Studies of Sepsis
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)  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 agentbased 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.

Limitations
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.

Concluding Remarks
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.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author https://github.com/ jarobarb/code_for_sepsis.git.

ETHICS STATEMENT
The animal study was reviewed and approved by the University of Pittsburgh Institutional Animal Use and Care Committee.

AUTHOR CONTRIBUTIONS
JA and YV both equally contributed as senior authors, with JA supervising the modeling components and YV supervising the experiments. RN performed experiments and collected and compiled the experimental data set. JA, JB, TB, AC, and AT developed the model equations and analyzed model behavior. JB finalized all model calibrations to experimental data and statistically analyzed model fits and subsequent predictions. JB, JA, YV, and RN wrote the article. All authors contributed to the article and approved the submitted version.

FUNDING
YV, National Institutes of Health (NIH; www.nih.gov), NIH Grant U01EB021960. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the article. JA, NIH Grant 1R01EY030851. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the article. JA-National Science Foundation (NSF; www.nsf.gov) NSF DMS-1654019. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the article. JA and JB, NSF DMS-1852146. The grant provided financial support for the REU projects performed by AC and AT that eventually led to this publication. JB, NSF DMS-1951531. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the article.