|Home | About | Journals | Submit | Contact Us | Français|
Unlike many localized infections, the development and resolution of bacteremia involves physical and immunological interactions between many anatomic sites. In an effort to better understand these interactions, we developed a computational model of bacteremia as a dynamical system fashioned after multicompartmental pharmacodynamic models, incorporating bacterial proliferation and clearance in the blood, liver, spleen, and lungs, and the transport of pathogens between these sites. A system of four first-order homogeneous ODEs was developed. Blood and organ bacterial burdens were measured at various time points from 3 to 48 h postinoculation using an LD25 murine model of Staphylococcus epidermidis bacteremia. Using these empiric data, solutions to the mathematical model were recovered. A bootstrap resampling method was used to generate 95% confidence intervals around the solved parameters. The validity of the model was examined in parallel experiments using animals acutely immunocompromised with cyclophosphamide; the model captured abnormalities in bacterial partitioning previously described with this antineoplastic agent. Lastly, the approach was used to explore possible benefits to clinically observed hyperdynamic blood flow during sepsis: in simulation, normal mice, but not those treated with cyclophosphamide, enjoyed significantly more rapid bacterial clearance from the bloodstream under hyperdynamic conditions.
Staphylococcus epidermidis, a member of the coagulase-negative staphylococci, is among the most frequently identified causes of bloodstream infections in the United States, a problem related to increasing use of short- and long-term implantable devices (1, 2). In recent clinical series, this group of organisms causes between 11% and 33% of nosocomial bloodstream infections, with an incidence of 15.8 per 10,000 hospital admissions (3, 4), and is a disproportionately frequent cause of bloodstream infection in intensive care units (5). Perhaps the strongest risk factor for developing a bloodstream infection is neutropenia. Short-term crude mortality after a neutropenic bloodstream infection exceeds 30%, and attributable mortality has been estimated to be approximately 12% in patients with hematological malignancy (6, 7).
Surviving a bloodstream infection requires effective filtration of bacteria from the blood and elimination of pathogens from organs that have been seeded during episodes of bacteremia. Host defense in this context seems similar to the pharmacokinetic clearance of drugs, a phenomenon that has been successfully conceptualized as a series of compartments into which agents may enter, between which agents may be transferred, and within which they may be degraded. However, the analogy is not perfect; unlike drugs, bacterial burden can increase in either the early hours of infection when host defenses are coming online or in the event of unsuccessful host defense. Numerically analyzing complex experimental data taken from multiple compartments at multiple time points in a manner that imposes physiological and immunological relationships is a significant departure from standard statistical techniques, where the most common analyses conclude at the rejection of a null hypothesis.
Considering bacterial clearance from the bloodstream in a physiologically based pharmacokinetic (PBPK) modeling paradigm was first described in 1983 (8). In that report, a limited 2-compartment model was fit to measures of blood and mesenteric lymph node bacterial content in mice over time. The authors assumed that bacteria did not proliferate during the experiment and computationally had the capacity to fit only the mean of the measured data from multiple animals. Accordingly, the conclusions of the work were limited to the observation that bacteria are cleared from the bloodstream in a manner approximated more closely by functions of a higher order than linear or single-exponential forms.
Improved methods of quantitatively thinking about host-pathogen interactions during life-threatening bloodstream infection could provide a useful framework for both experimental and analytical work in this area. Therefore, in the current work, we substantially expanded the concept of bacteremia as a PBPK phenomenon. A four-compartment mathematical model that related bacterial growth, clearance, and intercompartmental transfer between the blood, lung, liver, and spleen was used to characterize a murine LD25 S. epidermidis bacteremia model. This model was used to study two questions. First, what are the effects of acute immunocompromise induced by the antineoplastic agent cyclophosphamide as they relate to the resolution of an episode of bacteremia? Second, how does hyperdynamic blood flow affect the conclusions drawn from this computational model? Although more technically demanding than traditional statistical methods, analyzing experimental results from an animal model of multiorgan infection with a PBPK mathematical approach provides an opportunity to pursue novel questions the answers to which would otherwise be inaccessible.
Staphylococcus epidermidis RP62A, a strain originally isolated from an infected intravascular catheter, was obtained from American Type Culture Collection (ATCC 35984) (9, 10). Klebsiella pneumoniae Xen 39 is an engineered strain with a chromosomally inserted constitutively expressed lux operon and therefore can be visualized in vivo using high-sensitivity cameras (Xenogen, Inc, Hopkinton, Mass). Bacteria were grown from cryopreserved stock in tryptic soy broth or on solid (TSA) media (MP Biomedicals, Solon, Ohio) and incubated at 37°C. To obtain midlog phase cells, overnight cultures from single colonies were subcultured in fresh tryptic soy broth by 1:100 dilution and grown for 2 to 3 h until optical density (OD600nm) readings reached 0.4. Cultures were washed, resuspended in sterile 0.9% NaCl, and quantified turbidometrically (11). Six- to 10-week-old outbred male ICR mice (Harlan Sprague-Dawley, Indianapolis, Ind) were used in all experiments. Animal protocols were approved by the University of Michigan’s Animal Use Committee. The experiments were perfomed in adherence to the National Institutes of Health Guidelines on the Use of Laboratory Animals.
Mice were anesthetized using 2.5% inhaled isoflurane (Vedco, Inc, St Joseph, Mo) and injected into the deep dorsal vein of the penis with 200 μL of the bacterial suspension containing 109 colony-forming units (CFUs) of S. epidermidis. The magnitudes of bacterial inocula were confirmed by plating serial dilutions. In preliminary experiments, this dose was determined to produce a 7-day LD25 infection. At specified time points, animals were euthanized by exsanguination and cervical dislocation under deep isoflurane anesthesia. The entire liver, spleen, and lungs, and 1 mL of blood were harvested in a sterile fashion. Portions of each organ were weighed in tared tubes, homogenized in 3 mL of 0.9% NaCl, and then quantitatively cultured on TSA. For survival analysis, a separate cohort of animals was injected with the same dose of bacteria and monitored daily for 9 days.
Leukopenia was induced in some animals via intraperitoneal injection of cyclophosphamide (200 mg/kg; Fluka Biochemika 29875, Switzerland) 72 h before the study (12). Neutropenia was confirmed by automated complete blood cell count and by peripheral blood smear analysis. Neutropenia was maintained by giving a second injection of cyclophosphamide at the time of bacterial inoculation.
For analysis of very early postinoculation uptake, we used luminescent bacteria injected i.v. into isoflurane-anesthetized mice positioned in a dedicated imaging system (IVIS Imaging System, Xenogen/Caliper Life Sciences, Hopkinton, Mass). Image series were recorded using default instrument settings and 30-s exposures. Five animals were studied under identical experimental conditions. Unfortunately, luminescent S. epidermidis was not available at the time of this writing. Therefore, we alternatively chose K. pneumoniae Xen 39.
To establish the upper bound for the bacterial growth rate in whole blood, midlog phase bacteria were diluted to various concentrations, added to 200 μL whole mouse blood (1:10 by volume), and incubated at 37°C. At specified time points, 20-μL samples were taken from each tube and used for quantitative cultures on TSA. Growth constants (as h−1) were determined assuming first-order growth kinetics.
Bacteria within infected mice were assumed to exist as a well-mixed suspension of CFUs with diameter = 0 (i.e., no explicit relationship between bacterial and capillary dimensions was included) and to reside only within four compartments: the blood, liver, spleen, and lung (Fig. 1). Within each compartment, bacteria proliferated at rate a (a ≥ 0) and were eliminated at rate d (d ≥ 0), both as per hour (h−1), in proportion to the number of bacteria in the compartment. As there was no experimental means to distinguish a from d in vivo, the two terms were treated as a net bacterial growth rate (a–d) in the model. The upper bounds of a, and thus a–d, was established experimentally to be 0.92 h−1 in the whole blood bacterial growth assay described previously. Transfer between blood and other compartments was described as a function of the rate of blood flow, q, to each compartment, the volume, v, of each compartment, and a partitioning coefficient, p, where p is greater than 1 when bacteria were being filtered from the blood, and 0 is less than p is less than 1 when bacteria were being shed into the bloodstream from the compartment. Like proliferation or elimination, partitioning from one compartment into another was proportional to the number of bacteria present in the donor compartment.
These features were assembled into a system of four autonomous ordinary differential equations (ODEs) describing the instantaneous rate of change of bacterial burden in each compartment,
where l, s, h, and b are the bacterial concentrations, as CFU/mL of tissue, in the lung, spleen, liver (h, hepatic), and blood, respectively. Whereas circulating blood could transfer bacteria into the spleen, the spleen transferred bacteria only into the liver via the portal circulation. Hepatic circulation was divided into two components: splenic portal blood flow and the combination of hepatic arterial flow and all nonsplenic portal flow. Values for organ volume, v, were taken from organ mass measurements in the current experiments with an assumption that tissue density in each case was 1 g/mL. Values for organ blood flow, q, were linearly scaled by body mass from previously published microsphere-based flow measurements in rats with bacterial peritonitis (13). Modeling was performed using both normodynamic and hyperdynamic blood flow assumptions (Table 1).
Software for model fitting and all analysis was written in Matlab R2006b (The Mathworks, Inc, Natick, Mass). To generate median and 95% confidence interval (CI) estimates of each fitted parameter, a bootstrapping method was used as has been discussed elsewhere (14). For each bootstrap replicate, the liver, spleen, lung, and blood bacterial content of one mouse was selected randomly, with replacement, from each time point to generate a time series of multicompartment bacterial burden. For this time course of measurements, the model equations were solved using the appropriate ODE solver in Matlab. The parameters were estimated by fitting the solution of the differential equations to the bootstrap data using least squares optimization. In this way, each bootstrap replicate produced one estimate of all net bacterial growth rates and intercompartmental partitioning coefficients. One thousand replicates of this procedure were performed. An identical process was undertaken for neutropenic animals. Importantly, the model’s initial conditions were taken to be the first postinoculation time point (i.e., 3 h postinfection). This decision was made to avoid early intercompartmental redistribution that presumably occurred at rates too fast to be reliably measured by our intermittent culture technique.
To determine whether parameter estimates produced a stable solution (i.e., one which would ultimately lead to a fixed point, in this case, zero bacteria in all compartments) or an unstable one (i.e., in which the bacterial content increased in a nonsurvivable fashion), stability analysis was conducted using the Routh-Hurwitz criteria as has been described previously (15). Details of this procedure are presented in the appendix.
To establish the very early dynamics of bacterial distribution and clearance in this model, we used a luminescent strain of K. pneumoniae to conduct real-time imaging studies. Strain differences notwithstanding, sequential whole-body imaging revealed that within 30 s, bacteria began concentrating in the region of the liver. Within 4 min, the lung, liver, and spleen all demonstrated intense uptake of blood-borne bacteria (Fig. 2). Based on these data, the blood and these three organs were the targets of our subsequent modeling. In addition, given the very rapid distribution of bacteria in the opening minutes of a bacteremia episode, we concentrated our experimental measures and modeling on the slower process of bacterial removal.
When later time points (3–48 h) were studied by quantitative culture rather than bacterial luminescence, pathogen burden was found to decrease over time in all compartments, although at different rates (Fig. 3). In particular, liver bacterial burden increased over the first 24 h, whereas it decreased in other compartments. Clearance dynamics were captured by estimated growth rates and intercompartmental partitioning coefficients; bootstrap estimates of these parameter values and their 95% CIs are shown in Figure 4 and in Table 2. Of note, the model suggests that the lung, and not the liver, had the highest bacterial killing rate. A bacterial growth rate of −4.70 h−1 for the lung indicates that 90% of the accumulated burden would be expected to be cleared every 30 min in the lung, whereas the liver rate of −0.24 h−1 suggests that 90% of accumulated bacteria would be cleared in roughly 9 h. However, the bacterial partitioning in the liver was significantly more efficient (×8-fold), supporting the claim that despite differences in organ perfusion, the liver is the main site of bacterial clearance during a bacteremic episode, a conclusion in line with extensive previous experimental observations (16–18).
As expected, cyclophosphamide administration 72 h in advance of inoculation produced significant reductions in circulating leukocytes, particularly neutrophils (Figs. 5, A–C). Seventy percent (14/20) of normal mice survived 7 days after bacterial inoculation, whereas no bacteremic animals pre-treated with cyclophosphamide survived 24 h (0/22, Fig. 5D). Changes in compartmental bacterial burden and model parameter estimates for immunocompromised animals are shown alongside their normal counterparts in Figures 3 and and44 and in Table 2. Interestingly, the most prominent differences between normal and cyclophosphamide-treated animals were among intercompartmental partitioning coefficients, rather than bacterial clearance rates. This finding supports previous reports that neutrophils participate in bacterial clearance from the bloodstream by more than bactericide alone, and that the immune lesion produced by cyclophosphamide extends to the Kupffer cells and vascular endothelium (16, 18, 19).
An important feature of systems of ODEs is stability, the existence of a range of model parameters that, over time, will drive the state of the system toward a fixed point. In the current context, stable solutions were those that would ultimately result in complete eradication of bacterial burden and thus were a prerequisite for host survival. Conversely, unstable solutions were those that would over time produce bacterial burdens approaching infinity and necessarily would be fatal to the host. In the current model, the boundaries between stable and unstable solutions were a family of five membranes that existed in the 6-dimensional space whose axes were the three net growth rates (ah − dh, as − ds, and al − dl ) and three partition coefficients (ph, ps, and pl). As described in the appendix, the system of equations can be rewritten in matrix-vector form, with parameter coefficient matrix M. Based on our analysis, one of the membranes, along which the determinant of the matrix M was equal to zero (det[M] = 0), was found to separate the solutions of normal and cyclophosphamide-treated mice, and therefore defined the cyclophosphamide-treated parameters as necessarily unstable and not survivable. In Figure 6, ,2-dimensional2-dimensional sections of the 6-dimensional parameter space are shown with density contour plots of the 1,000 numerical solutions from normal and cyclophosphamide-treated mice and the calculated boundary between stable and unstable solutions. The solutions for the two groups lay on opposite sides of the stability boundary det(M) = 0, indicating that the solutions lead in the case of the normal animals to complete elimination and, in the case of chemotherapy-treated animals, to uncontrolled bacterial proliferation.
In Figure 7, the median parameters for normal and cyclophosphamide-treated mice were applied with either normodynamic or hyperdynamic blood flow values q (Table 1). For both experimental conditions, the model was provided with mean compartmental bacterial burdens measured at 3 h as its initial conditions. The bacterial burden over the ensuing 15 h (18 h postinoculation) was then forecast. Hyperdynamic blood flow offered a significant advantage to normal mice, accelerating the elimination of bacteria from the bloodstream. However, given the impaired partitioning described by the mathematical model, there was little quantitative difference when these blood flow values were applied to the cyclophosphamide-treated group.
In the current work, a physiologically based pharmacokinetic model was used to characterize S. epidermidis bacteremia in mice. The fitted models agreed well with experimental data, and the use of bootstrap resampling of empiric measurements allowed statistical comparisons to be made between normal animals and those with a clinically important test condition, cyclophosphamide-induced leukopenia. Formal stability analysis of the differential equations used to describe host-pathogen interactions demonstrated that model parameters in cyclophosphamide-treated mice were not only quantitatively different from those of normal mice, but lay on the far side of an analytical boundary beyond which bacterial proliferation would progress unchecked. Surprisingly, the parameters most altered by cyclophosphamide were not the rates of bacterial clearing, but of blood:organ partitioning, indicating defects in the shuttling of bacteria out of the bloodstream and not just in bacterial killing.
That bacterial disappearance from the bloodstream follows a pattern exponential decay has been recognized since the 1950s after classic studies using radiolabeled heat-killed Escherichia coli and Staphylococcus aureus (20). In those studies, serial measures of blood radioactivity were tracked and found to follow first-order kinetics. Because of the nature of those experiments, no multicompartmental modeling was possible. Mathematical limitations also existed at the time-the first description of a 2-compartment model for any pharmacological application (intercompartmental movement of inulin) had only been reported in 1949, and until the 1960s, pharmacokinetics was largely a nascent field (21, 22). In 1983, Cheewatrakoolpong and colleagues (8) used a 2-compartment model to explain the disappearance of K. pneumoniae from the bloodstream and from mesenteric lymph nodes in mice. They did not treat these anatomic sites as two interrelated compartments, but modeled bacterial counts from each site independently using equations of the form:
where B1 (analogous to our a – d) was the rate constant for the measured compartment (the blood or lymph nodes), and B2 was the rate constant for a hidden unmeasured second compartment. This approach allowed better fitting of experimental data than did a single exponential form but could not provide mechanistic insight into the underlying dynamics. In addition, because of the destructive nature of sampling (an issue in the current work as well), those investigators chose to fit only the mean values of replicate mice killed at each time point. As a result, their numerical approach provided estimates of A and B in the previous equation but no variance around these estimates.
The modeling strategy in the current work sought to overcome these shortcomings by requiring that all numerical solutions to fit simultaneous measurements of bacterial burden in each of the modeled compartments and by generating a large number of bacterial time courses using a resampling method such that parameter estimates, their variances, and joint probabilities could all be generated. The result is a model that can shed light on the process of bacterial clearance during an episode of bacteremia.
An important test of validity is the ability of a mathematical model to demonstrate behavior for which it has not explicitly been designed. In the case of this multicompartmental model, solutions obtained for cyclophosphamide-treated mice support the validity of our approach; pretreatment with this agent seemed to affect the shuttling of bacteria out of the bloodstream, not just to reduce the rate at which bacteria were killed. This was most notable in the liver and lung, where the difference was approximately an order of magnitude. Our experimental data, and our necessarily simplified model of intercompartmental transfer, do not reveal details about the blood-organ interface. However, in the liver, cooperativity between resident Kupffer cells and recruited circulating neutrophils against bacteria in the portal circulation is well described (23). Failure to compartmentalize circulating bacteria because of a lack of available neutrophils for this partnership is consistent with our results. To our knowledge, a similar situation in the lung intravascular compartment, perhaps between marginated monocytes and neutrophils taking place on an ad hoc basis, has not been described. Particularly intriguing are recent descriptions of neutrophil extracellular traps, tangles of neutrophil DNA in which bacteria may become embedded. Recent evidence indicates that after LPS stimulation, neutrophils that marginate in the lung and liver may deploy extracellular DNA as a means of increasing the efficiency of filtration of blood-borne bacteria (24). Our model predictions seem consistent with this experimental observation.
Another benefit of the current model is its ability to draw a clear distinction between normal and cyclophosphamide-treated mice, in an analytical, rather than purely statistical, way. We were able to take advantage of stability analysis to dichotomize model behavior into two outcomes after acute infection: progressive elimination, where total bacterial burden over time approaches zero, and progressive overwhelming proliferation, where the burden approaches infinity. These two outcomes are associated with two distinct sets of model parameters that lay on opposite sides of a multidimensional boundary defined by the Routh-Hurwitz criterion for system stability. The existence of distinctly survivable and lethal trajectories in a model system allows one to pose experimentally testable questions such as what therapies might shift a host’s model parameters from the lethal side of this boundary to the surviving side.
Lastly, by including blood flow in our model, we were able to consider the contribution of hyperdynamic blood flow during severe infection. Hyperdynamic sepsis, in which cardiac output and regional blood flow are increased over uninfected levels, is a common clinical occurrence. Using scaled blood flow values from rats with bacterial peritonitis, we explored the impact of increased organ blood flow on the multicompartmental behavior of bacteria during systemic infection. As might be predicted, elevated blood flow rates accelerated the disappearance of bacteria in normal animals. Because partitioning parameters in cyclophosphamide-treated mice were so much lower than in normal mice, increased blood flow in this context did not significantly benefit the host.
All models of illness, whether mathematical or experimental, fall short of fully reproducing clinical human experience, and our model of disseminated bacterial infection has several noteworthy limitations. The first is the nature and intensity of our bacterial challenge. We chose an i.v. inoculation system for this work because we wanted to capture the behavior of a host exposed to a single episode of bacterial showering into the bloodstream. This model is admittedly artificial but allowed us to develop a mathematical approach that did not require the inclusion of additional “black box” terms that could account for ongoing bacteria seeding that could not be measured experimentally.
A second limitation is that the mathematical model is autonomous or time invariant. It does not have a means of capturing host or pathogen behaviors that change over time; a bacterium entering the system 12 h after infection experiences the same host as one present in the initial inoculum. In reality, infected hosts recruit humoral and cellular defenses that require time to come on line, and hemodynamic changes, beneficial or harmful, also occur with some delay. We have shown previously in a murine model of pneumonia that the rate of bacterial elimination within the lung increases in the opening hours of infection as neutrophils are recruited into the airspaces (25). A reasonable next step in developing the theoretical model will be to incorporate a reliably measurable phenomenon, such as organ neutrophil content, that presumably will track temporally changing bactericidal or partitioning events.
Another unavoidable limitation of the current strategy is the destructive nature of measuring bacterial burden in four compartments. Linear dynamical modeling attempts to estimate a trajectory of future behavior based on the current state of a system. In our experimental system, organ and blood cultures were obtained by killing of the experimental animal, and as such, it was not possible for any mouse to contribute data at two time points. It was therefore impossible to know exactly an individual’s trajectory. We chose to overcome this problem by generating in essence a simulated repeated measures design by bootstrap resampling the available experimental data. This strategy provided the additional benefit of producing estimates of parameter variance, a crucial issue for future work that aims to test more complicated hypotheses or to make sample size calculations for work in involving advanced mathematical techniques. The validity of this bootstrapping is well documented, but the larger statistical implications of this choice of numerical method for data sets such as the one used here deserve, in our view, considerable additional study and development (14).
Lastly, we would point out that the modeling in this article, both experimental and computational, seeks to address only one aspect of clinical bacteremia. In humans, the downstream inflammatory response to a bacterial challenge is frequently more devastating than the bacterial exposure itself. Considerable work is underway to better characterize all aspects of this response in an analytical way (26–28). The current work is one additional step toward that goal.
The authors thank Adam P. Thompson for his assistance in the laboratory and the staff of the Center for Advanced Computing at the University of Michigan. The code used to generate Figure 6 was modified from a routine written by Peter Perkins, Mathworks, Inc, and is available at the Mathworks File Exchange, www.mathworks.com/matlabcentral/fileexchange. This work was supported by National Institute of General Medical Sciences grant R01 GM069438 and an associated supplement for computational biology. Additional funds were provided by a pilot grant from the University of Michigan Center for Computational Medicine and Biology. Ms. Chung was sponsored by a Society for Academic Emergency Medicine/Emergency Medicine Foundation student grant.
To compare the predicted fates of normal and cyclophosphamide-treated mice, we determined the conditions for stability of the model using standard analytical methods (15, 29). Our 4-compartment system can be expressed in matrix form as:
or in abbreviated form as:
where the 4 × 4 matrix M is the system’s coefficient matrix (e.g., m11 = [a−d]l−ql/vlpl). When this model is fit to experimental data, the output of the fitting procedure are quantities for each of the variables m11 through m44. To establish stability, the characteristic equation, p(λ), of the coefficient matrix was determined to be:
The Routh array for this equation was generated using the Symbolics Toolbox in Matlab. The first column of the Routh array was found to be:
where α1 and α2 are more complicated algebraic terms. For a solution to be stable, all five of these elements must be greater than zero. Numerical estimates of the second through fifth elements for normal and chemotherapy-treated mice, and for normal and chemotherapy-treated mice with hypothetical hyperdynamic blood flow, are shown in Table A1.
These results suggested that det(M) = 0 was the dominant stability boundary in our system. We therefore solved det(M) = 0 for various model parameters to locate the boundary between stability and instability in relationship to our fitted parameters. As 1,000 bootstrap solutions from both the normal and cyclophosphamide-treated animals were available, the boundary was calculated for all 2,000 sets of parameters. The median and 95% confidence bounds were used for the plots in Figure 5 in the main text.