Article tools
  • View PDF
  • Cite this article
  • E-mail this article
  • Send a letter
  • Information on reuse
  • Standard abbreviations
  • Article usage
Author information
Need help?

Research Article

Evidence for bistable bacteria-neutrophil interaction and its clinical implications

Roy Malka1, Baruch Wolach2,3, Ronit Gavrieli2, Eliezer Shochat4 and Vered Rom-Kedar1

1Department of Computer Science and Applied Mathematics, Weizmann Institute of Science, Rehovot, Israel. 2Laboratory for Leukocyte Function and 3Department of Pediatrics, Meir General Hospital, Kfar-Sava, Israel. 4Pharma Research and Early Development, Hoffmann-La Roche, Basel, Switzerland.

Address correspondence to: Vered Rom-Kedar, Department of Computer Science and Applied Mathematics, Weizmann Institute of Science, Rehovot 76100, Israel. Phone: 972.8934.3170; Fax: 972.8934.4122; E-mail:

First published July 23, 2012
Submitted: February 3, 2012; Accepted: June 7, 2012.

Neutropenia, which may develop as a consequence of chemotherapy, increases the risk of bacterial infection. Similarly, increased risk of bacterial infection appears in disorders of phagocytic functions, such as the genetic disorder chronic granulomatous disease. To elucidate the organizing principles behind these distinct immunodeficiency conditions, we investigated the interaction between in vitro bacteria and human neutrophils by experiments and mathematical modeling. The model and the experiments showed that the in vitro bacterial dynamics exhibit bistability for a certain range of neutrophil concentration and function. Thus, there is a critical bacterial concentration above which infection develops, and below which neutrophils defeat the bacteria. Whereas with normal neutrophil concentration and function, an infection may develop when the initial bacterial concentration is very high, under neutropenic conditions or when there is neutrophil dysfunction, the critical bacterial concentration can be lower, within the clinically relevant range. We conclude that critical bacterial concentration has clinically relevant implications. The individual maximum bearable bacterial concentration depended on neutrophil concentration, phagocytic activity, and patient barrier integrity; thus, the resulting maximal bearable bacterial concentration may vary by orders of magnitude between patients. Understanding the interplay between neutrophils and bacteria may enhance the development of new therapeutic approaches to bacterial infections.

See the related Commentary beginning on page 2776.


Neutropenia is common in cancer patients treated with intensive chemotherapy regimens. When the peripheral neutrophil blood count is below 500 × 103/ml, there is a dramatic increase in the risk of severe pyogenic infections (1), often resulting in organ damage and death (2). The reasons for this well-documented critical threshold are not fully understood. In addition, conditions associated with neutrophil dysfunction, even with normal neutrophil counts, are known to lead to inability to fight infections (35).

To explain these clinical findings, in vitro bactericidal experiments have been conducted. Initially, these in vitro studies suggested that the occurrence of infection depends on the ratio between bacterial concentration (B) and neutrophil concentration (N) (68); we refer to this herein as the ratio-dependent model (Figure 1A). Further investigations of Staphylococcus epidermidis–neutrophil interaction in vitro (9, 10), in conjunction with a linear mathematical model, suggested that there exists a single critical neutrophil concentration (CNC), below which neutrophils are completely unable to control bacterial growth, and above which they can control any bacterial concentration (neutrophil-threshold model; Figure 1B). Most recently, mathematical models that take into account nonlinear effects such as saturation in bacterial growth and in neutrophil killing rates were proposed and analyzed (11, 12), resulting in the monostable (Figure 1C) and bistable (Figure 1D) models.

Bacteria versus neutrophils: 4 plausible behaviors.Figure 1

Bacteria versus neutrophils: 4 plausible behaviors. Red regions correspond to experimental B and N concentrations that lead to an increase in the bacterial population (B locally wins); blue regions correspond to concentrations that lead to a decrease in the bacterial population (N locally wins). (A) Ratio-dependent models propose that the initial B/N ratio determines the result (68). (B) The neutrophil-threshold model proposes that a single neutrophil threshold exists, below which the bacterial population grows, and above which the bacterial population decreases (9, 10). (C) Monostable models propose that when the bacterial natural growth conditions are mild, the presence of neutrophils leads to a gradual decrease in the bacterial maximum capacity, down to extinction (11). (D) Bistable models propose that when the bacterial natural growth conditions are more favorable, the presence of neutrophils leads to bistable behavior (see Figure 1 in ref. 11). Here, we focused on the experimentally relevant region. The BNC (dashed and solid black curves) separates regions of bacterial growth from regions of bacterial decrease. Only in the bistable case (D) does the BNC have a positive slope, not conserving the B/N ratio.

Generally, monostability refers to a system with a single stable steady state, whereas bistability refers to a system that has 2 stable states for some range of parameter values. Here, N serves as a controlling parameter, and the 2 possible pronounced (stable) states are bacterial population capacity and no-bacteria or extinction state. In the bistable regime, both states are stable and are separated by an additional intermediate unstable steady state. B values above this unstable state grow to the population capacity concentration, whereas concentrations below it decay to the extinction state. In the bistable case in particular, for N above some threshold, the ability of neutrophils to overcome bacterial infection depends on the initial B (Figure 1D). It was shown previously that the neutrophil bactericidal experiments of Li et al. using S. epidermidis (10) are consistent with such bistable behavior (11).

These behaviors may be represented mathematically by ordinary differential equations governing the rate of change in B, where N is held fixed. The rate of change in B may be divided into the natural growth function of B (which depends on B and on the control experimental conditions) and to the bacteria-neutrophil interaction term. The interaction term depends on B, N, and neutrophil function in the given experimental condition. These 4 behaviors may be realized by considering different classes of natural growth terms and interaction terms (see ref. 12 and Methods). Such models were extended to explore the bistable bacteria-neutrophil dynamics under the influence of chemotherapy and G-CSF intervention by allowing N to vary (see refs. 12, 13, and Methods).

Here, we designed and performed an experiment to examine whether in vitro Staphylococcus aureus bacteria growth under attack from human neutrophils exhibits bistability. Our most important finding was a strong indication of bistability in all data examined, independent of the explicit forms of the mathematical models. We further used these data to parameterize the bistable bacteria-neutrophil equation via a specific previously proposed model (12), and the parameterized equation indeed fit the in vitro data. The parameterized model was then used to study the bacterial dynamics under neutropenic conditions. Under mild prolonged neutropenic conditions, the bacterial threshold for infection (i.e., maximum bearable bacterial influx) was found to be sensitive to the observed variability in the bactericidal activity (BA) of the neutrophils.


The S. aureus–human neutrophil interaction is bistable. We performed a series of bactericidal experiments in which S. aureus bacteria were incubated for 60 minutes in 24-well plates in suspensions containing 10% autologous sera and fixed N values. After 60 minutes, the neutrophils were lysed to stop the interaction, and the bacteria were assayed by colony culturing (see Methods). The ranges of experimental B and N values were chosen to detect the bistable regime. This series of experiments was repeated independently 4 times, each time using neutrophils and sera from a different healthy volunteer; each of the 4 experiments contained around 20 data points. Our experimental findings are summarized in Table 1.

Table 1

Bactericidal experiment

Next, we examined which of the 4 plausible behaviors shown in Figure 1 was consistent with our data. In all cases, there was a critical bacteria-neutrophil curve (BNC) that separated regions of bacterial growth from regions of reduction. On this curve, B is fixed: bacterial growth is exactly balanced by neutrophil elimination. In the ratio-dependent behavior, the BNC is a straight line through the origin, along which the B/N ratio is constant (Figure 1A). In the neutrophil-threshold behavior, the BNC is a vertical line (cutting the B axis at the CNC described in refs. 9, 10; Figure 1B). In the monostable case, the BNC is a decreasing function of N (Figure 1C). In the bistable case, for the relevant range of B, the BNC is an increasing function of N: bacterial load increases if the initial B exceeds the BNC, but decreases if initial B falls below the BNC (Figure 1D).

Here, we estimated the experimental BNC and established that its properties were in best agreement with the bistable behavior. Note that the bistable case is the only one of the plausible behaviors in which this curve has a positive slope and a nonconstant B/N ratio. For the ratio-dependent case, the B/N ratio along the curve is constant; for monostable and neutrophil-threshold behaviors, the curve has negative and vertical slopes, respectively.

By using the sign test (see Methods), we demonstrated that the experimentally derived BNC had a positive slope (Figure 2). This strongly indicated that the neutrophil-threshold and monostable behaviors are inconsistent with the data. Additionally, the B/N ratio along the experimentally derived curve was not constant (see Methods and Figure 3). Hence, the ratio-dependent behavior is inconsistent with our data and with the data of Li et al. (10). Therefore, we concluded that of the 4 plausible bacteria-neutrophil behaviors, only the bistable behavior is consistent with our data using S. aureus in suspension and with the previously reported data using S. epidermidis in gel (10).

Sign diagrams.Figure 2

Sign diagrams. At each initial experimental data point [B(0),N], we indicate whether the final bacterial population in each of the 2 duplicates after 60 minutes has a larger or smaller value than the initial one (+ or O, respectively). As each initial data point corresponds to 2 experiments (duplicates; see Methods), 2 symbols are plotted for each point. Points with 2 different signs appear near the BNC, where small experimental deviations between the duplicates may lead to opposite outcomes. (A) A slanted line (here in log scale) divided the symbols better than a vertical line (n = 78). (B) S. epidermidis data (n = 34), adopted from ref. 10. (C and D) Addition of the BNC (black curve), determined from the fitted mathematical model of Equation 3 and Table 2 (the dashed BNC curve of Figure 1D in a logarithmic plot). The entire data set (n = 86) is presented in C, and the single subject P1 (n = 22) in D, with the extent of increase/decrease (i.e., magnitude of Y = log[B(60 minutes),N/B(0)]) shown by color. As expected, near the BNC, the rates were close to 0 (greenish symbols).

PBK diagram.Figure 3

PBK diagram. PBK is plotted as a function of N, with lines connecting data points having the same experimental B. The points at which these PBK lines crossed the N axis were designated N*[B(0)]. The B/N ratios at the N*[B(0)] points approximated the B/N ratios on the BNC; thus, (B(0),N*[B(0)]) belongs to the experimentally extracted BNC. (A) Growth rate lines for S. aureus with experimental B of 104, 106, and 107 CFU/ml for the collective data set (n = 86). Data represent mean ± SD. Note that for experimental B of 104 CFU/ml, the intersection occurred at N* ≈ 105, in contrast to N* ≈ 2 × 106 for 107 CFU/ml; that is, the B/N ratio changed by at least 1 order of magnitude along the estimated BNC. (B) Growth rate lines for the S. epidermidis data (n = 34) adopted from ref. 10 with the indicated experimental B values. For experimental B of 2 × 103 CFU/ml, the intersection occurred at N* ≈ 2 × 105, in contrast to N* ≈ 2 × 106 for 1.6 × 108 CFU/ml; that is, the B/N ratio changed by at least 3 orders of magnitude along the BNC.

The above findings were qualitative and independent of the explicit mathematical form of the relevant models. We next used the data to parameterize a specific bistable model of the bacteria-neutrophil interaction (11). Table 2 presents the parameter fit for the S. aureus data for each of the 4 subjects, in which neutrophil functionality is measured with the autologous sera of each subject; for the pooled data set; and for the S. epidermidis data extracted graphically from the study of Li et al. (10). A large clinical study is needed to establish the optimal functional form of a bistable model and the typical average and variability of the model’s parameters (see Discussion). Nonetheless, there was good correspondence between the parameterized model and the fitted data sets (see Methods). Hence, we concluded that the parameterized model is adequate for describing the available data sets.

Table 2

Personal and pooled parameter estimates for the model in Equation 3

For each of these data sets, we also determined the predicted Nc, the point at which the fitted BNC intersects the B = 0 axis. Below this CNC, a small initial B would grow exponentially, whereas above it, a small initial B would be cleared. Whereas our in vitro experiments were conducted at B values of 104 CFU/ml and above, this prediction provides insights regarding the dynamics at much lower B. In particular, when N is very close to the predicted Nc, the small blood B of septicemic patients (which may reach 102–103 CFU/ml; refs. 1417) may be above the BNC, and thus develop to an acute infection. In our sample of 4 healthy individuals, the average Nc for the 4 subjects was 557 × 103 neutrophils/ml, with a coefficient of variation of 18.5%. Notably, the standard classification for grade IV/severe neutropenia is an absolute neutrophil count (ANC) less than 500 × 103 neutrophils/ml (2).

The chemotherapy-induced severe neutropenia scenario. To better illustrate the possible clinical implications of the variability in the model parameters, we considered the following common clinical situation. The neutrophil levels of a patient drop for a period of time and then recover. While this neutropenic cycle occurs, the patient is exposed to a constant bacterial influx. Such a scenario often arises after chemotherapy, as the neutrophil counts are temporarily reduced due to bone marrow damage, and bacterial influx may be increased due to degradation of the integrity of the patient’s barriers. We examine whether the variability in neutrophil function (i.e., in the model parameters shown in Table 2) may be critical in such clinical scenarios.

To examine such issues, we incorporated our experimental data set into the bacteria-neutrophil equation of our previously proposed in vivo mathematical model (see Methods and ref. 12). This model was designed specifically to study scenarios of infection when the neutrophil supply is limited, as in neutropenia. We used it here to isolate the influence of neutrophil functionality and barrier integrity from all other factors. First, we chose the parameters governing blood N of a hypothetical patient that undergoes chemotherapy and receives G-CSF support (see refs. 12, 13, and Methods). The chosen parameters provided a neutropenic time profile [N(t)] of a patient with prolonged severe neutropenia. Without G-CSF support, the patient’s N would be below the critical threshold of 500 × 103 neutrophils/ml for 10 days, whereas with 10 daily G-CSF injections, N would be mostly above the threshold of severe neutropenia, yet would remain neutropenic (i.e., N less than 1,000 × 103 neutrophils/ml) for about 7 days (Figure 4A). Note that for this hypothetical patient, standard daily shots suffice to abrogate the severe neutropenia, and he belongs to the lower G1 class previously reported (13), with acute marrow capacity of 0.15.

Evolution of infection in hypothetical patients during a 21-day cycle of chFigure 4

Evolution of infection in hypothetical patients during a 21-day cycle of chemotherapy-induced neutropenia. (A) Hypothetical N with no G-CSF support (untreated) and with 10 daily injections (treated). Thresholds for neutropenia (1,000 × 103 neutrophils/ml) and severe neutropenia (500 × 103 neutrophils/ml) are shown by dashed and solid red lines, respectively. (BD) Bacterial dynamics for 3 hypothetical patients with increasing bacterial influx (s). Hypothetical patients P1-treated and P4-treated received G-CSF support and thus had mild neutropenia. Their neutrophil functionality parameters are taken to be those of volunteers P1 and P4, respectively (see Table 2). Hypothetical patient P1-untreated did not receive G-CSF support and thus had severe neutropenia. His neutrophil functionality parameters are taken to be those of volunteer P1. (B) P1-treated could bear the 300-CFU/ml/h load, whereas for P4-treated and P1-untreated, this small s led to the development of acute infection. (C) With a more than 350-fold increase in s, P1-treated still succeeded in overcoming the infection. (D) With s increased 400-fold, all hypothetical patients developed an infection, regardless of G-CSF support.

We then examined the resulting bacterial growth for 3 hypothetical patients. The first hypothetical patient received G-CSF support, and his neutrophils’ functionality parameters were set to be the same as those of volunteer P1; therefore, we refer to this patient hereafter as P1-treated (Figure 4A). P1-treated overcame both small and quite large bacterial influx contaminations; in fact, a 400-fold increase of the small bacterial influx was needed to cause this patient to get an acute infection (Figure 4, B–D). The second hypothetical patient, P4-treated, had neutrophil functionality parameters set to those of volunteer P4. P4-treated also received G-CSF support and had a neutropenic profile identical to that of P1-treated, so his N was mostly above 500 × 103 neutrophils/ml. Given their similarities, the change in response caused by different neutrophil functionality was dramatic: P4-treated could not overcome even the smallest bacterial influx (Figure 4B). Finally, the third hypothetical patient, P1-untreated, had the high neutrophil BA plus the autologous serum of volunteer P1 in the absence of G-CSF support. For P1-untreated, neutrophil counts were below the threshold value for almost 10 days (Figure 4A), as expected, which indicated that the hypothetical patient could not overcome even the smallest bacterial influx.

Because we have no information about the real responses to chemotherapy of healthy volunteers P1 and P4, we examined a case in which P1-treated and P4-treated have identical changes in N in response to chemotherapy and G-CSF support. We further assumed that the neutrophil functionality of each volunteer was unchanged by the chemotherapy. The mathematical model allowed us to examine such hypothetical medical scenarios and isolate the influences of specific factors. The simulations demonstrated that interpatient variability in neutrophil functionality with autologous sera alone led to vast differences in the response to different levels of bacterial influx. Even with G-CSF support, P4-treated could not overcome a bacterial influx of 300 CFU/ml/h (Figure 4B); importantly, as described above, some septicemic patients demonstrate blood B of 102–103 CFU/ml (1417). Conversely, P1-treated overcame a 350-fold increase of this influx (Figure 4C). Put differently, for a given neutropenic profile, the patient’s neutrophil functionality determined the critical bacterial influx above which severe infections developed. It was previously shown that in general, this critical bacterial influx depends sensitively on N, neutrophil functionality, and bacterial growth rates (12). Here, we demonstrated that with the recorded interpatient variability, this critical value changed by at least 2 orders of magnitude. We concluded that this phenomenon is dramatic and is expected to be clinically observable.

The above results applied to neutropenic patients with neutrophil counts close to the medical threshold value of 500 × 103 neutrophils/ml (Figure 4A). Simulations of the same model with noncritical neutropenic profiles do not exhibit such parameter sensitivity (12). Indeed, when the neutrophil nadir counts were far above (270% increase) or far below (28% decrease) the threshold value, all patients eventually recovered or rapidly developed acute infections, respectively.


Here, the bistable model was shown to provide an organizing structure relating neutropenia-associated conditions, neutrophil dysfunctions, and bacterial influx. It fully explained the current knowledge of in vitro bacteria-neutrophil dynamics and was corroborated by our experimental data. Using a mathematical model, we examined the clinical implications of these findings on infection development in neutropenic patients and in patients with malfunctioning neutrophils. We propose that all major clinical observations regarding the development of infections in neutropenic patients are consistent with our model.

First, the severity and duration of the neutropenic nadir are known to be correlated with increased risk of infection (1, 18). The existence of the BNC indicates that the severity and duration of the nadir directly influence the likelihood of developing a severe infection, thus explaining the observed correlation.

Second, at some ranges of nadir levels and nadir durations, the risk of infection greatly increases (1, 18). For these near-critical ranges, there was strong interpatient variability in infection development. Our model provided a mechanism for handling such variability: the interpatient variability in neutrophil function led to variability in the individual critical Nc and thus to strong variability in infection development among patients with similar neutropenic profiles.

Third, sterility and isolation are known to reduce the risk of infection in some neutropenic patients (19, 20). In our in vivo model, the barriers’ integrity and the sterility conditions influenced one combined parameter: bacterial influx. Occasional contamination appears as a local sudden change in B. The existence of the BNC may explain the sensitivity of neutropenic patients to occasional contamination, hence the need for isolation. Aseptic conditions help to keep the bacterial load below this critical curve. Occasional large bacterial contamination may push the bacterial load above the BNC and lead to a severe infection. The shift of this curve to higher N values when bacterial influx is incorporated (12) explains the importance of sterile conditions for neutropenic patients. This shift means that the same occasional contamination that may be kept under control under sterile conditions would lead to severe infection when the body fights a constant bacterial influx from a nonsterile environment. It also explains the increased risk of infection in chemotherapy-induced neutropenia compared with other neutropenic conditions; indeed, mucosal barrier injury is a known side effect of radiotherapy and chemotherapy treatments (21). Finally, our model showed that in the case of severe neutropenia (N < Nc), even a small bacterial influx (e.g., from the intestines) may cause infection; thus, in this case, isolation cannot prevent the development of infections (22).

Fourth, the bistable behavior may explain specific puzzling clinical scenarios, such as the following 2 cases. The first patient, an immunocompromised neonate (1 month old), also suffers from severe congenital neutropenia (SCN), with ANC less than 200 × 103 neutrophils/ml. Supportive therapy with G-CSF significantly improved neutrophil count and function (23). Together with antibiotic therapy, the treatment controlled life-threatening infections. The second patient, an adult with acute lymphoblastic leukemia who received chemotherapy, developed a high fever with severe neutropenia (ANC of 380 × 103 cell/ml). He received specific antibiotics and G-CSF. Despite the recovery of neutrophil counts, the patient died. As shown in Figure 4, our model demonstrated that the disparity in these scenarios may be explained both by the degradation of the adult’s barriers by chemotherapy and by the possible relatively weak functionality of the adult’s neutrophils, caused by either normal interpatient variability or the effect of chemotherapy on neutrophil function.

In addition to neutrophil count and functionality, other factors affect the efficiency of bacterial clearance in cases of septicemia or bacterial growth in the blood: these include bacterial clearance in liver and spleen, blood monocyte concentration, and blood opsonin concentration. Nevertheless, neutrophil BA is known to be essential and amenable to medical intervention. Indeed, the common medical practice for avoiding microbial spreading in neutropenic conditions is early restitution of neutrophil levels by G-CSF treatment as well as administration of antibiotics, as described above. Clinical observations indicate that G-CSF increases both neutrophil counts and neutrophil function (2325); however, other studies found variability in the effect of G-CSF on in vitro neutrophil BA (26). Finally, with respect to the neutrophil functionality of the volunteers studied here, we measured the functionality of their neutrophils with their autologous sera, which was extracted at the time the blood sample was taken. In future studies, the effect of sera on neutrophil functionality may be explored.

Clinical insights regarding CGD patients. Chronic granulomatous disease (CGD) is an immunodeficiency syndrome caused by defects in NADPH oxidase, the enzyme that generates ROS in phagocytizing leukocytes. The lack of ROS generation is caused by mutations in any one of the genes encoding structural components of NADPH oxidase. CGD patients have an impairment in NADPH oxidase processing, which leads to life-threatening infections and consequently to granuloma formation in different tissues. This granuloma formation is now considered to be a hyperinflammatory reaction to various stimulants in CGD patients, causing granulomatous lesions that eventually result in organ dysfunction (27).

The diagnosis of CGD was established by neutrophil functional and genetic studies (3, 28). In the standard functionality test, S. aureus at a final B of 4.4 × 106 CFU/ml is incubated with 4.4 × 106 neutrophils/ml for 60 minutes, after which the neutrophils are lysed and B is assayed. Since most CGD patients lack or have very low ROS production, their BA (calculated as –log[B(t)/B(0)]) is low, at less than 0.5. Equivalently, these patients have a percentage of bacteria killed (PBK; calculated as [1 – B(t)/B(0)] × 100; e.g., ref. 7) that is smaller than normal, at 68%.

Recently, a 5-year-old patient suffering from repeated infections was tested (28). His BA of 0.62 was not considered compatible with CGD. Inspired by the bistability findings, the test was repeated with a larger initial B of 2 × 107 CFU/ml. The BA of the patient decreased to 0.05, whereas the BA of the control decreased to 0.16. Given that on the BNC, both measures of neutrophil functionality (BA and PBK) vanished (Figure 3), the higher initial B was closer to the BNC for the patient and further below it for the healthy control; that is, the patient’s BNC was lower than that of the control. This difference may explain the patient’s repeated infections. We suspect that the low BA of most CGD patients indicates that their BNCs are considerably lower than normal in both autologous and homologous sera.

Clinical perspective. Of the 998 patients referred for immune evaluation at the Laboratory for Leukocyte Function in Meir General Hospital due to repeated infections, the etiology of the immunodeficiency could be established in only 35% of the cases (29). We propose that some of the patients’ undefined disorders are related to a combination of several mild impairments, each one by itself lying within the lower end of the normal population distribution. Designing a protocol for testing immunodeficiency may help in identifying such patients.

Potentially new diagnostic uses. The bistability paradigm suggests that the standard tests for immune evaluation could be optimized by using assays with multiple B and N values. The following protocol may be considered: (a) take a single blood sample; (b) use the multiconcentration bacteria-neutrophil assay (see Methods) using both autologous and homologous sera; (c) find the experimental BNC of the patient for each serum (see Methods); (d) estimate the individual parameters of a bistable model for each serum (see Methods); and (e) compare the individual autologous BNC and the individual parameters with those of the homologous, a control, and a corresponding control population.

We expect that the proposed protocol would be useful to 2 groups of patients. First, for patients with repeated infections associated with deficiencies in neutrophil functionality, we expect that the autologous BNC will be lower than the control’s BNC and that some of the model parameters may be abnormal. Such an assessment may help to determine the proper therapy, for example, by providing a quantitative prediction for the needed G-CSF and antibiotic support. Moreover, comparison of the fitted parameters with the various controls may possibly help to diagnose the deficiency.

Second, for patients needing to undergo chemotherapy treatments that bear significant neutropenic risk, the individually parameterized model and the resulting BNC may help in an a priori assessment of a patient’s critical zone. We propose that such an assessment may help to identify those individuals at high risk of infection due to their intrinsic immune vulnerability. Moreover, such an assessment may help to determine the needed support therapy, including the isolation protocol (accounting for both patient neutrophil count and BNC). Finally, we propose that by repeating the above-described 5-step protocol at several time points after chemotherapy (see ref. 30), separately fitting at step (e) the parameters at each time point, one may examine whether a significant change in neutrophil function has occurred.

Implementation of these ideas in the clinic requires additional empirical studies. To begin with, a large population of healthy subjects should be assayed, and different common experimental settings should be explored; for example, different laboratories use different percentages of autologous serum (9) and different types of bacteria. Then, mathematical forms of bistable models that provide the best fit to the gathered data should be selected — the previously developed model (11) used herein presents only one of many other options (see ref. 12). These experimental and theoretical explorations will generate estimated population parameters and provide their between-subject variability. With these estimations at hand, we envision the suggested protocol will provide a diagnostic tool based on a standard objective measure. We hope that the emerging principle — that the combined effect of mild impairments may lead to significant clinical consequences — can contribute to the current practice for treating patients that suffer from repeated infections or undergo chemotherapy treatments.


Neutrophil preparation

Heparinized blood was obtained by venipuncture from healthy donors. Medicated volunteers or those that suffered from an infection-related condition in the 2 weeks prior to sample were excluded. Isolated neutrophils were prepared as described previously (31). Neutrophil purity was greater than 95%, as determined using the Advia 120-Automated Hematology Analyzer (Bayer). Purified neutrophils were resuspended to a final N of 1 × 107 neutrophils/ml in PBS supplemented with 0.2% glucose and 0.1% bovine serum albumin (Sigma-Aldrich) (PBS-G-BSA). Cells were kept at room temperature and used within 1 hour.

Serum preparation

Autologous human serum was prepared from 3 ml whole blood (without anticoagulant) from the same donor whose neutrophils were used in the experiment. After 20 minutes at room temperature, serum was separated from the clot by centrifugation at 1,500 g for 5 minutes and then removed to a clean tube. We use the 10% serum, which has proven to be the optimal concentration for S. aureus opsonization in our experience and that of others (32).

Bacteria preparation

S. aureus (502A-ATCC 27217) were freshly grown before each experiment and allowed to enter an early stationary phase (4 hours, with continuous shaking, at 37°C), as previously described (6). Bacteria were suspended with PBS-G-BSA and kept at 4°C until use. The final B was determined spectrophotometrically at 590 nm, by reference to a previously determined growth curve relating the CFU of S. aureus to A590.


BA (killing) was assayed as previously described (6, 32), with modifications. In brief, 100 μl bacteria (final B, 104–107 CFU/ml) with 10% autologous serum (vol/vol) was placed in a 24-well plate and incubated with continuous shaking at 37°C for various durations, according to the specific lag time of each experimental B [see below; baseline value of B measured after the concentration-dependent lag time is referred to throughout as B(0)]. When the bacteria entered the growth phase, 100 μl neutrophils (final N, 2.5 × 105 to 44 × 105 neutrophils/ml) were added and incubated for another 60 minutes under the same conditions. The concentration of viable bacteria was assayed. Neutrophils were lysed with distilled water at 11.0 pH; the suspensions were diluted and plated in triplicate on mannitol salt agar plates (Novamed Diagnostics) for 18 hours at 37°C, after which the colonies were counted.

Experimental setup

Our experiments followed a previously described protocol (32); the novelty in the present study is the judicious, mathematically motivated choice of initial data points. Here, we briefly discuss the general setting and elaborate on specific relevant issues. Several preliminary bactericidal experiments were performed to identify the range of N and corresponding B at which bistable behavior could be obtained. These preliminary experiments also suggested that 2 important ingredients must be calibrated first. Recall that the experiment is designed to demonstrate the competition between rapid natural bacteria growth and neutrophil phagocytosis at varying concentrations. Hence, we must verify that the bacteria population is indeed growing when no phagocytes are present and that the methodology for stopping the experiment is not concentration dependent. We first studied the natural growth curves of S. aureus under the above experimental setup and determined the concentration-dependent lag time of each experimental B: 80 minutes for 104 and 105 CFU/ml, 60 minutes for 106 CFU/ml, and 40 minutes for 107 CFU/ml. As expected, the lag time of S. aureus decreased with increasing experimental B. In addition, we chose a lysis method that did not exhibit concentration dependence (lysis with distilled water at 11.0 pH; ref. 32 and data not shown).

We then performed 4 sets of experiments. In each of the experiments, S. aureus and fresh neutrophils from a healthy individual were inserted into wells in which the bacteria were multiplying (i.e., after the corresponding lag time). B(0) at the time of neutrophil insertion was measured (in 2 wells), 2 control wells with no neutrophils were set, and an array of duplicate wells with varying initial N was set. The experiment was stopped after 60 minutes by pH lysis, and CFUs were then counted. Table 1 shows all the reliable experimental results. Only 1 case was eliminated due to bacteria overgrowth, which prohibited reliable CFU counts.


Each of the 4 different qualitative behaviors in Figure 1 may be represented by many mathematical models. Our goals were to (a) examine which of these behaviors is consistent with the data and (b) examine the suitability of the parameterized bistable model of Table 2 to describe the data set.

We used 2 different sets of statistical and mathematical tools. First, a probabilistic nonparametric sign-test model was constructed and used to establish that it was highly probable that the critical curve (i.e., the BNC) that splits the data between increasing and decreasing concentrations has a positive slope. By examining PBK curves, we then identified the estimated range of the experimental BNC and found it did not conserve the B/N ratio. Second, parameters were fitted to 3 variants of our nonlinear model. The predictive ability of the models was studied using nonparametric hypothesis testing (Kolmogorov-Smirnov, refs. 33, 34; Mann-Whitney-Wilcoxon, ref. 35; and Supplemental Figure 1; supplemental material available online with this article; doi: 10.1172/JCI59832DS1). The best model of the 3 variants was shown to provide a good fit to the data sets (Figure 5).

Parameter fit.Figure 5

Parameter fit. Using the P1 parameters of Table 2, we integrated the fitted model of Equation 3 and plotted the bacterial dynamics for experimental N of 2.5 × 105 (A), 5 × 105 (B), 7.5 × 105 (C), and 10 × 105 (D) neutrophils/ml. Solid curves are monotonically increasing and correspond to the bacterial growth regime; dashed curves correspond to the bacterial decrease regime. The data points of P1 (see Table 1) were well described by these curves. Plotted are [t*,B(0)] — the time t* at which the simulation crosses a B(0) value — and the corresponding value of B(60 minutes), at time t* + 60. As expected, B(60 minutes) values were very close to the fitted bacterial-dynamics curves.

Sign test. In each of the individual subject data sets, we found 2 or 3 N values for which the bacteria population increased at high experimental B and decreased at low experimental B (Table 1). This finding was consistent with the bistable and ratio-dependent behaviors and inconsistent with the neutrophil-threshold and monostable behaviors. To provide a robust statistical perspective regarding this observation, we used the nonparametric sign test. First, we assigned a binary signifier denoting growth or decrease to each data point of the S. aureus experiment (+ or O, respectively). Figure 2 depicts this coarse version of the data presented in Table 1. Figure 2, A and C, present all 4 experiments together, and Figure 2D presents the data set for volunteer P1 (see Table 1). Here, the BNC divided the [B(0),N] plane into a region of growth and a region of decrease. In statistical terms, curves like the BNC are called classifiers.

In this part of the analysis, we used only N at 2.5 × 105–106 neutrophils/ml, the concentrations for which there was a difference among the 4 possible behaviors; the data points at N of 0 and 4.4 × 106 neutrophils/ml corresponded to a behavior common to all 4 cases, and were thus excluded.

The sign test examines whether it is possible that by simple coin tosses, one could obtain the same success rate in distributing binary signifiers on either side of a given classifier. This is the null hypothesis (that fair coin tosses may explain the data). If this null hypothesis cannot be rejected with confidence (usually P values of 0.05 or lower), then the model is possibly no more predictive than a fair coin toss and is thus falsified. We examined whether this null hypothesis was acceptable for lines with positive slopes, nearly vertical lines, and lines with negative slopes. Since the data points are sparse, many classifiers have the same score; thus, classifiers of the simplest possible form must be chosen. Here, we selected the 2-parameter family of slanted lines (each defined by its slope and N intercept). One may view this choice as a rough approximation to the nonlinear BNC.

Briefly, we found that there were classifiers with positive slopes for which the null hypothesis could be rejected with confidence; the best of these had P values on the order of 10–7. Such lines were excellent classifiers in the sense that it is highly unlikely that their good division of increasing versus decreasing regions was accidental. The results were statistically more significant for the individual data sets than for the pooled data, possibly due to the variability in neutrophil activity among the different subjects (see Table 2). On the other hand, we found that for all the lines we took with nearly vertical or negative slope, the null hypothesis could not be rejected; the best line in this category had a P value of 0.43. In all cases, there was no statistical evidence that vertical lines (as in the neutrophil-threshold model; Figure 1B) or lines with a negative slope (as in the monostable model; Figure 1C) divided the data set into growing and decreasing regions with more predictive power than a fair coin. See Supplemental Methods for details.

Experimentally extracted BNC. As described above, on the BNC, bacterial growth was exactly balanced by the the elimination term; hence, on this curve, PBK exactly vanished. For a given experimental B, values of N to the left of the BNC led to bacterial growth and thus negative PBK, whereas values to the right corresponded to bacterial decrease and thus positive PBK. In other words, interpolation of the experimental data provides a rough estimate of the experimental BNC. Figure 3 depicts the graphical realization of such estimates. For each fixed experimental B, we used a line to connect the PBK results with increasing N values. The point at which this PBK line crossed the N axis was designated N*[B(0)]. We concluded that (B(0),N*[B(0)]) belongs to the experimentally extracted BNC.

In analyzing the PBK diagram for the S. aureus experiment (Figure 3A), we found that the B/N ratio along the BNC changed by at least 1 order of magnitude, for both individual and pooled data sets. This ratio changed by at least 3 orders of magnitude for the S. epidermidis data extracted from the study of Li et al. (ref. 10 and Figure 3B). In contrast, for any ratio-dependent model, the BNC must maintain a constant slope. We concluded that both data sets are inconsistent with ratio-dependent models. Moreover, in both experiments, the sequence of intersection points of the experimental B–PBK curves with the N axis for increasing initial B monotonically increased with N (see the perfect ordering of these curves’ intersection with the x axis in Figure 3). This behavior was consistent with bistable models. Combining the results of the sign test and the experimental BNC estimates, we concluded that the bistable behavior is the most suitable for describing the 2 data sets of bacteria-neutrophil dynamics.

Bacteria-neutrophil model and parameter fitting. We further used the data sets to parameterize a specific form of a bistable model. Then, beyond simple fitting, we explored the model predictions for certain in vivo scenarios.

The 6-parameter model shown in Equation 1 is consistent with reasonable biological assumptions regarding the in vitro behavior of bacteria under the attack of neutrophils (11).

(Equation 1)

All parameters (r, δ, β, α, γ, η, and s) were assumed to be non-negative. The first 3 parameters govern the natural dynamics of the bacteria: r and δ control the natural linear growth and death rates of the bacteria, respectively, and β controls the natural saturation of the bacterial growth rate at high concentrations. The other 3 parameters control bacteria killing by neutrophils: α is the neutrophils’ bacterial killing rate at low concentrations, and γ and η control the saturation and interference in the killing rate, respectively, as B and N increase. Finally, s corresponds to the bacterial influx from the environment; thus, it is set to 0 in the current setting of the in vitro experiment. For small s, this model is bistable in the case of Equation 2, and monostable when the inequality is reversed.

(Equation 2)

See ref. 12 for further details on this and other related models.

Below, we describe the process of fitting the parameters of this model to the data. We had 2 types of data sets: one corresponding to natural growth curves, and the other corresponding to bacteria-neutrophil dynamics. We first estimated the parameters associated with the natural bacterial dynamics, then used the bactericidal experiments with non-0 N values to find the parameters of the killing term.

Using the first data sets, we found that the bacteria exhibited roughly exponential growth (β = 0), with an estimated growth rate ρaureus (i.e., rd) of 0.027 (see below; for data from ref. 10, we used the authors’ estimate for the exponential growth rate ρepidermis of 0.013). Using these values in Equation 1, we arrived at the 3-parameter model shown in Equation 3, in which α > ρη (note that this natural condition emerged from the fitting and was not posed as a constraint).

(Equation 3)

For γ greater than 0, this model exhibits the bistable behavior, and its BNC is a ray with slope (α – ρη)/ργ, which emanates from Nc of ρ/(α – ρη). We fit the data from the bactericidal experiments to 3 different variants of the model shown in Equation 3. The first variant was achieved by setting γ and η at 0 and fitting the single parameter α. This linear variant was just the neutrophil-threshold model dB/dt = ρB – αNB, with the vertical BNC. The second variant was achieved by setting η at 0 and fitting the 2 parameters α and γ — a model admitting a BNC with a positive slope that does not incorporate the neutrophil interference effect. Finally, in the third variant, all 3 parameters, α, γ, and η, were fitted to the data.

In Supplemental Methods and Supplemental Figure 1, we compared in detail the goodness of fit of these 3 variants and their prediction error distributions. The model in which the kill term has saturation in both bacteria and neutrophil populations was superior to the others. We rejected the null hypothesis of equal distributions of the errors with P values of 0.05 using Kolmogorov-Smirnov and Mann-Whitney-Wilcoxon nonparametric tests (33). Thus, we concluded that the inclusion of neutrophil saturation in the kill term (i.e., interference) is indeed significant.

Figure 5 presents an example of the quality of the fit. The parameterized model provided a good fit to the experimental data set for patient P1. Similar fits were achieved in all other cases.

Figure 2, C and D, show the BNCs of the parameterized model embedded in the corresponding sign diagrams. The BNCs of the fitted models divided the increasing and decreasing data points with success rates of 80%–90%, much better than the nearly random success rates of the vertical and negative slope line classifiers. In fact, these were similar to the success rates of the best classifiers that were found by optimizing the sign test. Moreover, all misclassified points (i.e., + to the right/below or O to the left/above the BNC) had trend magnitudes close to 0. Measurement errors could explain these mismatches. Notably, this separation was an emerging property of the model, achieved with no optimization with regard to classifying the data points nor by specific constraints on the parameters.

Although other explicit mathematical forms might provide a better fit to the experimental data and to the classification problem, Figure 5 and the emerging classification property demonstrated that the fitted model (Equation 3) was adequate for describing the available data sets.

Mathematical model for the severe chemotherapy–induced neutropenia scenario. To produce Figure 4, the model in Equation 4 (see refs. 12, 13) was used to calculate N(t) and G-CSF concentration [G(t)] in blood, and thus the bacterial load [B(t)], for hypothetical patients with severe prolonged neutropenia.

(Equation 4)

Note that the equations for G(t) and N(t) do not depend on B(t): it was assumed that the neutropenia is so severe that neutrophil and G-CSF dynamics are externally controlled. Neutrophil influx depends on both chemotherapy and G, the latter of which was mainly controlled here by the pharmacokinetics of G-CSF injections. See Equation 2 and Supplemental Table 1 for the functional form of IN,G(t) and the parameters used in Equation 4.

To produce Figure 4, B–D, bacterial dynamics parameters were taken from Table 2 with constant s of 5, 1,966, and 2,000 CFU/ml/min, respectively (note that rates in Equation 4 are per minute). This model allows N in the tissue to immediately equilibrate with the blood levels and assumes that the neutrophils’ killing efficiency persists. In reality, adhesion, chemotaxis, and toxicity will further reduce the neutrophils’ ability to fight the bacteria and slow the transition rate. Hence, we propose that in neutropenic patients, this model will be useful to identify scenarios in which acute infection is likely to develop. Borderline cases and cases in which neutrophil dynamics are coupled to bacterial load need more detailed modeling. See ref. 12 for more details regarding the properties, predictions, and limitations of this model.

Study approval

The study protocol was approved by the local Committee on Human Volunteers at Meir Medical Center, in accordance with the Declaration of Helsinki. Blood donors provided written informed consent.

Supplemental data


The authors thank Yakar Kannai and Elchanan Mossel for helpful discussions and Dirk Roos and John Higgins for comments on this manuscript. This work was partially supported by the Israel Science Foundation (grant 273/07) and the Minerva Foundation. V. Rom-Kedar is the Estrin Family Chair of Computer Science and Applied Mathematics at Weizmann Institute of Science.


Conflict of interest: Eliezer Shochat is an employee of Hoffmann-La Roche.

Citation for this article:J Clin Invest. 2012;122(8):3002–3011. doi:10.1172/JCI59832.

Roy Malka’s present address is: Center for Systems Biology and Department of Pathology, Massachusetts General Hospital, and Department of Systems Biology, Harvard Medical School, Boston, Massachusetts, USA.

See the related Commentary beginning on page 2776.


  1. Crawford J, Dale DC, Lyman GH. Chemo­therapy–induced neutropenia: risks, consequences, and new directions for its management. Cancer. 2004;100(2):228–237.
    View this article via: PubMed CrossRef
  2. WHOICD-10: International Statistical Classification of Diseases and Related Health Problems. 10th ed. Vol. 1. Geneva, Switzerland: World Health Organization; 1992.
  3. Wolach B, et al. Chronic granulomatous disease in israel: clinical, functional and molecular studies of 38 patients. Clin Immunol. 2008;129(1):103–114.
    View this article via: PubMed CrossRef
  4. van den Berg JM, et al. Chronic granulomatous disease: The european experience. PLoS One. 2009;4(4):e5234.
    View this article via: PubMed CrossRef
  5. Kuijpers TW, et al. Natural history and early diagnosis of lad-1/variant syndrome. Blood. 2007;109(8):3529–3537.
    View this article via: PubMed CrossRef
  6. Clawson CC, Repine JE. Quantitation of maximal bactericidal capability in human neutrophils. J Lab Clin Med. 1976;88(2):316–327.
    View this article via: PubMed
  7. Leijh PCJ, van den Barselaar MT, van Furth R. Kinetics of phagocytosis of Staphylococcus aureus and Escherichia coli by human granulocytes. Immunology. 1979;37(2):453–465.
    View this article via: PubMed
  8. Hammer MC, Baltch AL, Sutphen NT, Smith RP, Conroy JV. Pseudomonas aeruginosa: quantitation of maximum phagocytic and bactericidal capabilities of normal human granulocytes. J Lab Clin Med. 1981;98(6):938–948.
    View this article via: PubMed
  9. Li Y, Karlin A, Loike DJ, Silverstein CS. A critical concenration of neutrophils is required for effective bacterial killing in suspension. Proc Natl Acad Sci U S A. 2002;99(12):8289–8294.
    View this article via: PubMed CrossRef
  10. Li Y, Karlin A, Loike DJ, Silverstein CS. Determination of the critical concentration of neutrophils required to block bacterial growth in tissues. J Exp Med. 2004;200(5):613–622.
    View this article via: PubMed CrossRef
  11. Malka R, Shochat E, Rom-Kedar V. Bistability and bacterial infections. PLoS One. 2010;5(5):e10010.
    View this article via: PubMed CrossRef
  12. Malka R, Rom-Kedar V. Bacteria-phagocytes dynamics, axiomatic modelling and mass-action kinetics. Math Biosci Eng. 2011;8(2):475–502.
    View this article via: PubMed CrossRef
  13. Shochat E, Rom-Kedar V. Novel strategies for g-csf treatment of high- risk severe neutropenia suggested by mathematical modeling. Clin Cancer Res. 2008;14(20):6354–6363.
    View this article via: PubMed CrossRef
  14. Werner AS, Cobbs CG, Kaye D, Hook EW. Studies on the bacteremia of bacterial endocarditis. JAMA. 1967;202(3):199–203.
    View this article via: PubMed CrossRef
  15. La Scolea LJ Jr, Dryja D. Quantitation of bacteria in cerebrospinal fluid and blood of children with meningitis and its diagnostic significance. J Clin Microbiol. 1984;19(2):187–190.
    View this article via: PubMed
  16. Wain J, et al. Day quantitation of bacteria in bone marrow from patients with typhoid fever: relationship between counts and clinical features. J Clin Microbiol. 2001;39(4):1571–1576.
    View this article via: PubMed CrossRef
  17. Peters RP, et al. Quantitative detection of Staphylococcus aureus and Enterococcus faecalis DNA in blood to diagnose bacteremia in patients in the intensive care unit. J Clin Microbiol. 2007;45(11):3641–3646.
    View this article via: PubMed CrossRef
  18. Rahman Z, Esparza-Guerra L, Yap HY, Fraschini G, Bodey G, Hortobagyi G. Chemotherapy-induced neutropenia and fever in patients with metastatic breast carcinoma receiving salvage chemotherapy. Cancer. 1997;79(6):1150–1157.
    View this article via: PubMed CrossRef
  19. Spiers AS, Dias SF, Lopez JA. Infection prevention in patients with cancer: microbiological evaluation of portable laminar air flow isolation, topical chlorhexidine, and oral non-absorbable antibiotics. J Hyg (Lond). 1980;84(3):457–465.
    View this article via: PubMed CrossRef
  20. Schimpff SC, et al. Infection prevention in acute nonlymphocytic leukemia. Laminar air flow room reverse isolation with oral, nonabsorbable antibiotic prophylaxis. Ann Intern Med. 1975;82(3):351–358.
    View this article via: PubMed
  21. van Vliet MJ, Harmsen HJ, de Bont ES, Tissing WJ. The role of intestinal microbiota in the development and severity of chemotherapy-induced mucositis. PLoS Pathog. 2010;6(5):e1000879.
    View this article via: PubMed CrossRef
  22. Fenelon LE. Protective isolation: who needs it? J Hosp Infect. 1995;30(suppl):218–222.
    View this article via: PubMed CrossRef
  23. Wolach B, Gavrieli R, Pomeranz A. Effect of granulocyte and granulocyte macrophage colony stimulating factors (g-csf and gm-csf) on neonatal neutrophil functions. Pediatr Res. 2000;48(3):369–373.
    View this article via: PubMed CrossRef
  24. Spiekerman K, Roesler J, Emmendoerffer A, Elsner J, Welte K. Functional features of neutrophils induced by G-CSF and GM-CSF treatment: differential effects and clinical implications. Leukemia. 1997;11(4):466–478.
    View this article via: PubMed CrossRef
  25. Roilides E, Walsh TJ, Pizzo PA, Rubin M. Granulocyte colony-stimulating factor enhances the phagocytic and bactericidal activity of normal and defective human neutrophils. J Infect Dis. 1991;163(3):579–583.
    View this article via: PubMed CrossRef
  26. Yang K, et al. Del- claux. Bactericidal activity response of blood neutrophils from critically ill patients to in vitro granulocyte colony-stimulating factor stimulation. Intensive Care Med. 2003;29(3):396–402.
    View this article via: PubMed
  27. Rieber N, Hector A, Kuijpers T, Roos D, Hartl D. Current concepts of hyperinflammation in chronic granulomatous disease. Clin Dev Immunol. 2012;2012:252460.
    View this article via: PubMed
  28. Wolach B, et al. Two X-linked chronic granulomatous disease patients with unusual NADPH-oxidase properties. J Clin Immunol. 2011;31(4):560–566.
    View this article via: PubMed CrossRef
  29. Wolach B, Gavrieli R, Roos D, Berger-Achituv S. Lessons learned from phagocytic function studies in a large cohort of patients with recurrent infections. J Clin Immunol. 2012;32(3):454–466.
    View this article via: PubMed CrossRef
  30. Hubel K, et al. Suppressed neutrophil function as a risk factor for severe infection after cytotoxic chemotherapy in patients with acute nonlym- phocytic leukemia. Ann Hematol. 1999;78(2):73–77.
    View this article via: PubMed CrossRef
  31. Böyum A. Isolation of mononuclear cells and granulocytes from human blood. Isolation of monuclear cells by one centrifugation, and of granulocytes by combining centrifugation and sedimentation at 1 g. Scand J Clin Lab Invest Suppl. 1968;97:77–89.
    View this article via: PubMed
  32. Green JN, Winterbourn CC, Hampton MB. Analysis of neutrophil bactericidal activity. Methods Mol Biol. 2007;412(21):319–332.
    View this article via: PubMed
  33. von Mises R. Mathematical Theory Of Probability And Statistics. New York, New York, USA: Academic Press; 1964.
  34. Press W, Teukolsky S, Vetterling W, Flannery B. Numerical Recipes in C. 2nd ed. Cambridge, United Kingdom: Cambridge University Press; 1992.
  35. Hogg RV, Craig AT. Introduction to Mathematical Statistics. 5th ed. Upper Saddle River, New Jersey, USA: Prentice-Hall; 1995.