## Abstract

NK cells are innate lymphocytes that mediate early host defense against viruses, such as cytomegalovirus. IL-15 is upregulated during viral infections and drives the expansion of NK cells. However, the influence of IL-15 on murine NK cell division and death rates has not been quantitatively studied. Therefore, we developed a series of two-compartment (representing quiescent and dividing NK cell subpopulations) mathematical models, incorporating different assumptions about the kinetic parameters regulating NK cell expansion. Using experimentally derived division and death rates, we tested each model’s assumptions by comparing predictions of NK cell numbers with independent experimental results and demonstrated that the kinetic parameters are distinct for nondividing and dividing NK cell subpopulations. IL-15 influenced NK cell expansion by modulating recruitment and division rates to a greater extent than death rates. The observed time delay to first division could be accounted for by differences in the kinetic parameters of nondividing and dividing subsets of NK cells. Although the duration of the time delay to first division was not significantly influenced by IL-15, the recruitment of nondividing NK cells into the replicating subpopulation increased with greater IL-15 concentrations. Our model quantitatively predicted changes in NK cell accumulation when IL-15 stimulation was reduced, demonstrating that NK cell divisional commitment was interrupted when cytokine stimulation was removed. In summary, this quantitative analysis reveals novel insights into the in vitro regulation of NK cell proliferation and provides a foundation for modeling in vivo NK cell responses to viral infections.

Natural killer cells are innate lymphocytes that recognize transformed or infected cells rapidly without prior sensitization. They are particularly important in early host resistance against viral pathogens while the adaptive immune system is being marshaled (1–3). NK cells are stimulated by cytokines (e.g., IL-15 and IL-12), as well as by direct recognition of infected or transformed cells through activation receptors, resulting in cytolysis of the abnormal cells and the production of immunomodulatory cytokines (e.g., IFN-γ).

In addition to these effector functions, NK cells rapidly proliferate during viral infections, resulting in increased NK cell numbers (4–7). This NK cell expansion is driven by pro-proliferative cytokines, such as IL-15, which are upregulated during viral infections (6, 8, 9). Exogenous IL-15 or overexpression of IL-15 in transgenic mice also stimulates NK cell proliferation, resulting in increased NK cell numbers (6, 10). The NK cell proliferative response during viral infections can be divided into three distinct phases: 1) early, nonselective, cytokine-driven NK cell proliferation; 2) preferential proliferation of a subset of NK cells expressing activation receptors that recognize infected cells; and 3) resolution of proliferation and contraction of the expanded NK cell numbers (5, 7, 9). In addition to driving viral-induced NK cell proliferation, IL-15 plays an important role in NK cell development, survival, and homeostasis (11–14). Dysregulation of NK cell proliferation and homeostasis can result in abnormally low levels of NK cells (13–15) or in lymphoproliferative disorders, such as chronic NK cell lymphocytosis (16–18).

The fundamental understanding garnered from qualitative models describing NK cell proliferation based on experimental studies (e.g., Refs. 5–7) can be substantially enriched and extended by quantitative insights from mathematical models. For example, Gett and colleagues established a mathematical framework of how T cells integrate signaling from cytokine, costimulatory, and TCRs (19), and three subsequent models built on this work to provide insight into the kinetic parameters controlling T cell proliferation (20–22). Furthermore, a recent study used mathematical modeling to characterize the integration of activation and inhibitory receptor signaling by Vav1 in the regulation of NK cell cytotoxicity (23).

Although mathematical models have been used to study T cell proliferation in vitro, no work has been done to determine how NK cell proliferation responses are quantitatively regulated. We had previously developed an in vitro CFSE-based assay of NK cell proliferation that was able to replicate the first two phases of viral-induced NK cell proliferation (7), and we were interested in using experimentally determined rate constants from this system in mathematical modeling to develop a more quantitative understanding of the kinetic parameters controlling NK cell proliferation. In this study, we present a novel two-compartment mathematical model to investigate the influence of IL-15 on the kinetic parameters governing NK cell proliferation and expansion. We demonstrate the model’s ability to predict changes in NK cell expansion that were verified with independent experiments. Observations from this model provide novel insights into the factors that regulate cytokine-driven NK cell proliferation.

## Materials and Methods

### Mice

Female C57BL/6 (B6) mice were obtained from the National Cancer Institute (Charles River, Wilmington, MA). They were maintained under specific pathogen-free conditions and used between 8 and 16 wk of age. All experiments were conducted in accordance with institutional guidelines for animal care and use.

### Cell preparation and culture

Single-cell suspensions were prepared from whole spleens using standard techniques (7). Splenic NK cells were enriched from ∼2.5% to 45–60% of lymphocytes after negative selection with magnetic beads (Miltenyi, Auburn, CA). Enriched NK cells were labeled with 0.25 μM CFSE (Molecular Probes, Eugene, OR) for 7 min at 37°C. The reaction was quenched by adding RPMI 1640 medium supplemented with 10% fetal calf serum, and the cells were subsequently incubated for 30 min at 37°C. The CFSE-labeled splenocytes were plated in 96-well plates (2.5 × 10^{4} NK cells/well) and cultured for 4 days with various concentrations of murine IL-15 (PeproTech, Rocky Hill, NJ).

### Quantifying the NK cell population

^{+} lymphocytes and enumerated with calibration beads (Spherotech, Lake Forest, IL) in conjunction with 7-aminoactinomycin (7-AAD; BD Pharmingen, CA). Flow cytometry analysis was performed with a FACSCalibur flow cytometer (BD Pharmingen, San Diego, CA), and the data were analyzed with FlowJo (Tree Star, Ashland, OR), Origin 7.5 (OriginLab, Northampton, MA), and Mathematica 7.0 (Wolfram, Champaign, IL).

### Estimating cell numbers and division progression

Precursor cohort analysis was conducted, as previously described (19–22), using OriginPro 7.5 software. In brief, a Levenberg–Marquardt algorithm was used to determine distinct Gaussian distributions representing cell cohorts with the same division history in the CFSE intensity histogram. The ratios of the area under each Gaussian distribution to the total area under the CFSE intensity histogram were used in calculating the distribution of NK cells in each distinct division. Precursor numbers were then calculated by dividing the number of cells present in division *i* by a factor of 2^{i}.

### Cytokine reduction

After NK cells were stimulated with 25 ng/ml IL-15 for 48 h, the majority of the IL-15 containing media inside each well was extracted. Approximately 30 μl media was left covering the NK cells at the bottom of each well. A total of 220 μl media without IL-15 was slowly added to prevent disruption of the NK cell cluster. This resulted in a final calculated concentration of 3 ng/ml IL-15. The NK cells were then cultured for an additional 48 h, and cell numbers were assessed periodically. In mock cytokine reduction experiments, the majority of the media was removed and 220 μl media containing 25 ng/ml IL-15 was slowly added, maintaining an IL-15 concentration of 25 ng/ml.

### Quality of prediction calculation

The quality of simulations and predictions from each mathematical model was assessed by computing the normalized root mean square deviation (NRMSD) between model computations and observed experimental values. The sum of squared deviations (SSD) at different time points were computed, and NRMSD was calculated as a percentage value by the following formula:An NRMSD value of 0% indicated absolute agreement between model computations and experimental results. Greater differences between model results and experimental data were reflected in higher NRMSD values.

### Models for NK cell proliferation

We hypothesized that quiescent NK cells had kinetic parameters distinct from those of actively dividing NK cells. Therefore, we devised a model that separated NK cells into two subsets or compartments: undivided (compartment U) and divided (compartment D) cells. Division and death rates were used to characterize the changes in the NK cell population in each compartment. Because division by compartment U cells represented recruitment into compartment D, the division rate of compartment U cells was equivalent to the recruitment rate. Three models were developed using systems of ordinary differential equations to solve for the number of cells in both compartments and to test the validity of the two-compartment hypothesis (Fig. 1).

Model 1

Model 2

Model 3

In these equations, *r* was the rate with which compartment U cells were recruited to compartment D, *d _{U}* was the death rate of compartment U cells,

*τ*was the time delay to first division, and

*k*and

*d*were the division and death rates of compartment D cells. (Mathematical symbols are defined in Table I.) The estimates

_{D}*r*,

*d*,

_{U}*τ*, and

*k*were obtained through precursor cohort analysis (19, 20, 22). The estimate for

*d*was derived using

_{D}*r*,

*d*, and

_{U}*k*.

*C*and

_{U}*C*were the total number of cells in compartments U and D, respectively. Initially,

_{D}*C*(0) was equal to the number of live NK cells at the start of the experiment, and

_{U}*C*(0) was equal to zero.

_{D}### CFSE histogram simulations

Simulated CFSE histogram profiles were generated from the differential equation predictions in our model using Matlab (Mathworks, MA). The experimental profile of NK cells after CFSE labeling was used as the initial cell distribution. CFSE data are organized by the flow cytometer into 1024 bins, corresponding to an intensity range of 1–10,000. The natural CFSE decay over the initial 24 h after CFSE labeling was accounted for in the simulation and resulted in the histograms shifting to a mean CFSE intensity of 1000–2000. During this initial observed decay, no dilution caused by cell division occurred. Corresponding to the limits of flow cytometer detection, up to 8 cycles of cell division were experimentally analyzed (20). Model differential equations were used with experimentally derived parameters to estimate cell distributions at specific time intervals. The number of NK cells that were recruited from the undivided compartment to the divided compartment, the number of NK cells that continued to divide, and the number of dead cells were calculated. These changes in the NK cell population were applied to the corresponding CFSE intensity bins using Matlab to define the simulated CFSE histograms at specific time points.

## Results

### IL-15 stimulation directly augments NK cell proliferation kinetics

We used [^{3}H]thymidine incorporation experiments to determine the appropriate range of IL-15 concentrations with which to study in vitro NK cell proliferation. [^{3}H]thymidine incorporation began to increase at 3 ng/ml IL-15 and approached maximum response levels at 200 ng/ml (Supplemental Fig. 1).

Based on these results, soluble IL-15 was added to CFSE-labeled NK cell cultures at concentrations of 3, 9, 25, and 75 ng/ml. NK cell numbers were assessed with flow cytometry at defined time points. Representative NK cell CFSE intensity histograms are shown for 30, 48, 61, and 78 h in Fig. 2. During the first 30 h of cell culture, no cell divisions were observed at any of the four IL-15 concentrations. Even after 78 h, negligible cell division was seen in NK cells incubated in 3 ng/ml IL-15, and the CFSE intensity distribution consisted of a single Gaussian peak (Fig. 2A). After stimulation with 9 ng IL-15/ml, a subset of NK cells divided up to four times; however, the majority of the NK cell population remained quiescent. Therefore, although the recruitment rate from the nondividing subpopulation of NK cells into the dividing population was greater after stimulation with 9 ng/ml than with 3 ng/ml, it remained relatively low. As IL-15 concentrations increased to 25 and 75 ng/ml, NK cells manifested progressively more divisions (Fig. 2C, 2D), reflected in increased CFSE dilution with multiple Gaussian peaks. After stimulation with greater IL-15 concentrations, a larger fraction of NK cells were found in later cohorts. This observation suggests that the division rate increased with increasing IL-15 concentration, but it is not possible to empirically distinguish the influence of division and death rates on the observed expansion of NK cells. To address this issue, we used a mathematical modeling approach incorporating two compartments representing quiescent and dividing subsets of NK cells.

### Two-compartment models

Three models were proposed to interpret and analyze the quantitative experimental data. Model 1 was designed to test the validity of the two-compartment hypothesis. In Model 1, the division rate of compartment D cells was equal to the recruitment rate of compartment U cells, and the death rates in both compartments were identical (Fig. 1A), representing essentially a single compartment system (20, 22). In contrast, Model 2 included a distinct division rate for compartment D cells and an explicit time delay to first division (Fig. 1B). The time delay to first division accounted for the observed lag in transition of NK cells from compartment U to D. Model 3 incorporated distinct death rates in both compartments, as well as distinct recruitment rate in compartment U and proliferation rate in compartment D. Model 3 tested the hypothesis that an explicit time delay was not necessary to account for the observed delay to cell recruitment in the context of distinct proliferation and death rates in compartments U and D (Fig. 1C).

### Model parameter estimation

Model parameter estimates were derived from CFSE experimental results using precursor cohort analysis. Cell division resulted in multiple CFSE “peaks” in CFSE histograms secondary to diluted CFSE in each generation of daughter cells. Each peak represented a cohort of cells that had divided the same number of times. These cohorts could be mathematically represented as a series of Gaussian distributions of CFSE intensities. *N _{i}*(

*t*) was defined as the number of cells in the

*i*

^{th}cohort, and the values of

*N*(

_{i}*t*) were determined as described in

*Materials and Methods*.

*C*was equal to the number of cells in the initial undivided cohort

_{U}*N*

_{0}.

*C*(the total number of cells in compartment D) was equal to the sum of cells in all subsequent cohorts.To determine the number of precursor cells (

_{D}*n*) from which the divided cohorts originated, we divided the cell numbers in each cohort by 2

_{i}^{i}.At the beginning of each experiment when all cells were in compartment U and no cells had divided, the precursor number

*n*

_{0}was equal to the number of cells in compartment U.

The starting NK cell population consisted of a finite number of precursors. After accounting for the decrease in precursors because of serial doubling, the only change to the precursor cell number was through cell death. Therefore, the rate of decrease in total precursor numbers over time was used to calculate an estimate of the death rate in compartment U (*d _{U}*):The term

*n*(

_{total}*t*) represented the total number of precursors in the NK cell population at time

*t*. Therefore, the slope of the linear regression of the ln[

*n*(

_{total}*t*)] versus

*t*plot provided an estimate of

*d*(Fig. 3A).

_{U}The recruitment rate of compartment U cells was calculated by tracking the changes in the number of undivided cells over time, because the rate of change in the number of undivided cells was the sum of recruitment and death rates.The recruitment rate of undivided cells (*r*) was therefore equal to the difference between the slope of the linear regression of ln[*C _{U}*(

*t*)] versus

*t*and

*d*(Fig. 3B).

_{U}The division rate of compartment D cells (*k*) was estimated using cell numbers from each cohort and the statistical concept of expected value (20). The value of *N*_{0} decreased with recruitment and death rates, as noted in the following equation:*N*_{0} represented the number of undivided cells, whereas *N _{i}* represented the number of cells that had divided

*i*times. Therefore, the average number of divisions that cells in compartment D had completed at time

*t*, [

*M*(

*t*)], could be calculated as follows:

*M*(

*t*) increased with time, and the slope of

*M*(

*t*) versus

*t*was equal to the rate of division of cells in compartment D (Fig. 3C).

Analysis of CFSE histograms in multiple experiments (Fig. 2) demonstrated minimal NK cell division at 30 h regardless of the IL-15 concentration. The time delay (*τ*) characterizing this prolonged time to first division was calculated by extrapolating the time at which the mean division number of precursor cells excluding undivided precursors, *m*(*t*), was equal to 1 (i.e., the average time to first division). The mean division number of precursor cells was calculated using the frequency distribution of precursor cells, *f _{i}*(

*t*) (i.e., the precursors that have divided

*i*times as a fraction of total precursor numbers). The concept of

*f*(

_{i}*t*) was based on normalized cohort cell numbers as described by De Boer and colleagues (20, 22).

*m*(*t*) was the sum of the product of precursor frequency distribution and precursor division number, excluding nondividing precursors using the factor 1 − *f*_{0}(*t*).The time delay (*τ*) was estimated as the time at which *m*(*t*) = 1 in a plot of *m*(*t*) versus *t* (Fig. 3D).

The values of *τ* ± SD calculated at 9, 25, 50, and 75 ng/ml IL-15 are shown in Fig. 3E. The average time delay to first division based on these experimental results was 32.5 h and appeared relatively independent of IL-15 concentration. However, the values calculated at 9 ng/ml IL-15 demonstrated substantial variability, suggesting that the time delay may be prolonged at lower IL-15 concentrations.

In contrast with the time delay, the fractional recruitment of NK cells from compartment U to D increased with increasing IL-15 concentrations (Fig. 3F). The fractional recruitment was defined as the ratio of precursors whose progeny were in compartment D to total precursors and was calculated by accounting for the fraction of precursors that remained in compartment U.Fractional recruitment was dependent on IL-15 concentration with an increasing fraction of the precursor NK cells in compartment U recruited to compartment D as IL-15 levels increased (Fig. 3F); however, fractional recruitment appeared to plateau as IL-15 concentrations increased to >50 ng/ml.

Although cell death was not directly measured in the CSFE experiments, the death rate of compartment D cells (*d _{D}*) was estimated by comparing the theoretical [

*N*

_{1}(

*t*)*] and observed [

*N*

_{1}(

*t*)] number of cells that had completed one cell division. The theoretical value

*N*

_{1}(

*t*)*, describing cells recruited from compartment U to D that had completed one division at time

*t*, was calculated using the recruitment rate

*r*and division rate

*k*, assuming no cell death in compartment D. Given a time delay of

*τ*, the number of cells that had not been recruited and remained in compartment U wasSimilarly, after time

*t*, the number of cells that had not been recruited wasThe nonrecruited cells consisted of quiescent and dead cells. If no recruitment occurred, the number of nonrecruited cells would remain constant over time. Therefore, the change in the number of nonrecruited cells was equal to the number of recruited cells.

*C*(0)

_{U}*e*

^{−}

*was used because recruitment was appreciable only after the apparent time delay period. Therefore, the number of cells recruited into compartment D at time*

^{rτ}*t*wasThe theoretical number of cells in first division at time

*t*must account for cells in compartment D that may have divided more than once by subtractingThus, for

*t*>

*τ*, the theoretical number of cells in division 1 [

*N*

_{1}(

*t*)*] was represented asThe observed value

*N*

_{1}(

*t*) was equal to the number of cells in the first division at time

*t*. Therefore, the death rate in compartment D (

*d*

_{D}) could be approximated from the slope of the linear regression of versus time. The values of

*d*at different IL-15 concentrations are listed in Table II.

_{D}The parameter estimates derived from CFSE experiments with 3, 25, and 75 ng/ml IL-15 are shown in Table II. Because of minimal cell division and transition to compartment D, only *d _{U}* and

*r*could be calculated for experiments using 3 ng/ml IL-15. At greater IL-15 concentrations, the proliferation rate of dividing NK cells (

*k*) was 3-fold greater than the recruitment rate (

*r*) of quiescent NK cells, suggesting that NK cells divide at a faster rate once they start to proliferate (Table II). Moreover, dividing NK cells also had greater death rates than quiescent NK cells after stimulation with either 25 or 75 ng/ml IL-15. Dead NK cells could also be directly identified by the uptake of 7-AAD (a fluorescent molecule that intercalates in double-stranded DNA but does readily pass through cell membranes). Estimates of death rates based on enumeration of 7-AAD

^{+}NK cells over time were consistent with our conclusion that NK cells stimulated with 75 ng/ml IL-15 had greater death rates than NK cells stimulated with 25 ng/ml IL-15 (Supplemental Fig. 2).

### Model simulation of experimental data

These experimentally derived proliferation and death rates were used in models to evaluate hypotheses about cytokine-driven NK cell proliferation. The model simulations of NK cell accumulation over 100 h were compared with experimental data from three independent experiments at 3, 25, and 75 ng/ml IL-15 (Fig. 4). Consistent with our experimental results (Fig. 2A), no NK cell expansion was predicted to occur after stimulation with 3 ng/ml IL-15 stimulation (data not shown). Model 1 used identical proliferation and death rates in compartments U and D (i.e., *k* = *r* and *d _{D}* =

*d*). The significant differences between the Model 1 simulation and experimental results (Fig. 4A, 4B) support the hypothesis that quiescent and dividing NK cells have distinct kinetic parameters. In contrast, the simulation results in Models 2 and 3 were much closer to experimental values at both IL-15 concentrations (Fig. 4C–F), as quantified by the NRMSD values. The similarity between the quality of the predictions from Models 2 and 3 when compared with experimental data suggests that differences in the recruitment rate, division rate, and distinct death rates between quiescent and dividing NK cells (Model 3) are sufficient to account for the observed initial lag in proliferation without the incorporation of an explicit time delay (Model 2).

_{U}### Model predictions of NK cell expansion under different IL-15 stimulation conditions

The ability to simulate experimental results is critical for a mathematical model, but the true test of this approach is in accurately predicting results under experimental conditions not used in deriving the parameter values. Therefore, we used parameter estimates determined from CFSE experiments with 3, 25, and 75 ng/ml IL-15 to generate logarithmic equations that defined the kinetic parameters as functions of IL-15 concentration (Fig. 5). These functions were used to interpolate parameter values at different IL-15 concentrations, including 9 and 50 ng/ml (Table II).

The predicted kinetic parameters at 9 and 50 ng/ml IL-15 were then used to calculate NK cell expansion in all three models, and the predictions were compared with results obtained from independent experiments (three experiments with 9 ng/ml IL-15 and two with 50 ng/ml IL-15; Fig. 6). Model 1 was unable to accurately predict NK cell expansion at either concentration of IL-15 (Fig. 6A, 6B), affirming the hypothesis that quiescent and dividing NK cells have distinct proliferation and death kinetics. In contrast, the predictions of Models 2 and 3 were comparable with much lower NRMSD values. Given that Model 3 could predict NK cell expansion, as well as Model 2 at IL-15 concentrations of 9, 25, 50, and 75 ng/ml, we concluded that the presence of distinct proliferation and death rates in the two-compartment NK cell proliferation model is sufficient to eliminate the necessity of incorporating an explicit time delay. Supporting this conclusion, Model 3 was able to accurately predict mean division numbers of NK cells at specific time points over the course of the experiment (Supplemental Fig. 3A). Therefore, we used Model 3 in the following work.

### Robust model predictions under dynamic cytokine stimulation conditions

Our two-compartment model was able to accurately predict the expansion of NK cells in response to a variety of steady-state IL-15 concentrations (Figs. 4, 6). However, in vivo IL-15 concentrations during stimuli such as viral infections are not static (6, 9), but rather are elevated and then return to baseline homeostatic levels. Therefore, we tested the robustness of our model in predicting NK cell numbers under conditions in which IL-15 stimulation was reduced to basal levels (3 ng/ml) after 48 h of stimulation with 25 ng/ml IL-15.

We hypothesized that the reduction of IL-15 levels would gradually modulate the kinetic parameters controlling NK cell expansion to values consistent with the new cytokine concentration. We estimated the length of time taken for cellular behavior to transition between the two IL-15 concentrations to be ∼12 h based on the following assumptions: 1) equilibration of IL-15 levels between the new media and the fluid remaining around the cell pellet would occur over a relatively short period, given a moderate amount of mixing when the new fluid was added; 2) ∼4 h would be required to reach steady-state levels of receptor occupancy at the new equilibrium IL-15 concentrations based on receptor-level modeling of soluble IL-15 binding to IL-15Rs on cell surfaces (following the approach outlined in Ref. 24; data not shown); and 3) cells already in the process of dividing would be minimally affected by changes in cytokine stimulation. Proliferating NK cells require at least 6 h to complete a division cycle (see Fig. 2C where NK cells stimulated with 25 ng/ml IL-15 undergo up to 5 rounds of division in ∼30 h after the they are recruited into the proliferating pool). We made the simplifying assumption that the temporal impact of the changing cytokine levels on cells already committed to dividing could be approximated by the time required for the cells in compartment D to undergo, on average, 0.5 division [*M*(*t*) = 0.5]. Given that *M*(*t*) = *kt* and *k* at 25 ng/ml IL-15 was 0.061 (Table II), the time for *M*(*t*) to equal 0.5 would be 8 h. Taking these three observations into account, we incorporated into the model the assumption that the IL-15–dependent kinetic parameters would vary as the effective concentration of IL-15 decreased linearly from the initial IL-15 concentration to the new concentration over a 12-h period.

With this modification, the model was able to accurately account for the expansion of NK cells under mock cytokine exchange and IL-15–reducing conditions when compared with independent experimental data (Fig. 7). The ability of the model to accurately predict NK cell numbers during mock exchange experiments (25 ng/ml IL-15 media replaced at 48 h with 25 ng/ml IL-15 media; Fig. 7A) demonstrated that the removal and replacement of media did not alter NK cell proliferation kinetics. When the IL-15 concentration was reduced to 3 ng/ml after 48 h of stimulation with 25 ng/ml IL-15, model predictions accurately accounted for the slowing and subsequent plateauing of NK cell expansion (Fig. 7B). The model was also able to approximate the decreased mean division numbers of NK cells in the cytokine reduction group compared with the mock exchange group observed during the experiments (Supplemental Fig. 3B). In addition, the model could predict the distribution of undivided and divided NK cells into distinct division cohorts. Indeed, model-simulated histograms of NK cell division cohorts compared favorably with CFSE histograms from mock exchange (Fig. 7C, 7D) and cytokine exchange (Fig. 7C, 7E) experiments. Consistent with experimental results, the simulated histograms of NK cell division cohorts representing the cytokine reduction experiments had fewer cells in cohorts that had undergone more than four divisions in addition to fewer NK cells overall (Fig. 7E). In summary, our two-compartment model was capable of robust predictions of both overall NK cell numbers and population distributions under changing IL-15 conditions.

## Discussion

To characterize NK cell proliferation and expansion, we developed a two-compartment mathematical model distinguishing quiescent and dividing populations of NK cells and used it to quantitatively investigate the proliferative response of NK cells to IL-15 stimulation in vitro. To our knowledge, this is the first work to quantitatively study the impact of IL-15 on the kinetic parameters (e.g., recruitment rate) governing cytokine-driven NK cell proliferation. In contrast with earlier modeling work on T cell proliferation that relied primarily on parameter estimates based on multiparameter, nonlinear curve fitting of experimental data (20–22), we directly derived rate constants describing recruitment, proliferation, and death rates within each compartment from experimental data. The predictions from our model were tested and verified in independent experiments. The analysis of the model provided a number of unique insights into the factors that control cytokine-driven NK cell proliferation: 1) IL-15 modulated recruitment and division rates to a greater extent than death rates, 2) increasing IL-15 concentrations paradoxically led to higher death rates in dividing NK cells, 3) the observed time delay to first division was not significantly influenced by IL-15 concentrations and could be directly accounted for by differences in the kinetic parameters of nondividing and dividing subsets of NK cells, 4) the recruitment of nondividing NK cells into the replicating pool increased with greater IL-15 concentrations, and 5) NK cell divisional commitment was interrupted when cytokine stimulation was removed.

IL-15 drives NK cell expansion by modulating proliferation rates (including recruitment rate from compartment U) and death rates. Both the proliferation and death rate of dividing cells increased with increasing IL-15 concentrations. However, increasing IL-15 had a greater impact on the proliferation rate than on the death rate (Table II), resulting in a net increase in NK cells as IL-15 concentrations increased. Although IL-15 is known to support NK cell survival (11, 12), our model parameter derivation found higher death rates as the IL-15 concentration increased (Table II), suggesting that rapidly proliferating NK cells undergo higher rates of death than quiescent cells. This conclusion was supported by estimates of NK cell death rates based on direct enumeration of 7-AAD^{+} NK cells (Supplemental Fig. 2). It is possible that this observation reflects physical constraints and/or increased cell–cell contact within the expanding cell mass of dividing cells rather than an inherent feature of IL-15–driven NK cell proliferation; however, experiments varying the NK cell density support the conclusion that replicating NK cells have a higher death rate than quiescent NK cells (Supplemental Fig. 4). Although unexpected, this observation is consistent with previous reports that the rate of T cell death is dependent on cell division number (22), and that rapidly dividing T cells are more susceptible to apoptosis (25–27).

Our model demonstrates that several characteristics of NK cell proliferation are similar to those of T cell proliferation. For example, our work found that the observed time delay to first division (*τ*) is relatively insensitive to IL-15 concentration (although more variability was seen at low IL-15 concentrations). This is consistent with the observation that varying IL-2 stimulation did not alter the average time to first division for CD3-stimulated T cells (21, 22). In contrast, the fractional recruitment (i.e., the proportion of cells entering first division) of both NK and T cells increased with increasing pro-proliferative cytokines (IL-15 for NK cells and IL-2 for T cells) (21, 22).

However, in a number of other ways, our model demonstrates that the proliferative behavior of NK cells is distinct from T cells. Our work confirms results from previous NK cell studies that pro-proliferative cytokine stimulation is sufficient to drive significant NK cell expansion (6, 7, 9) in contrast with T cells. Although T cells, to a certain extent, do undergo “bystander proliferation” (cytokine-driven division in the absence of signaling through their TCR) (28, 29), previous in vitro studies have found minimal proliferation with IL-2 alone (19, 21). In addition, our experiments demonstrate that NK cell divisional commitment is dynamic and can be interrupted when the cytokine milieu changes. This observation contrasts with previous reports that once CD8 T cells are stimulated, they enter into an apparent program of division and divide at least seven times without further antigen stimulation (30, 31). However, one caveat in this comparison is that the CD8 T cell studies were done with T cells stimulated through the TCR, whereas our studies of NK cells used stimulation solely by pro-proliferative cytokines and not through their activation receptors.

In summary, the combination of mathematical modeling and experimental observations allows the systematic interrogation of the influence of IL-15 on NK cell proliferation. Our quantitative study of the kinetic parameters controlling cytokine-driven NK cell proliferation and expansion demonstrated several similarities between NK and T cell proliferative behavior, but it also identified a number of differences between NK and T cell proliferation. In addition to the insights garnered about IL-15–driven NK cell proliferation, our model provides the basis for future studies incorporating signaling through NK cell activation receptors. A more complete understanding of the factors that govern NK cell proliferation and homeostasis may lead to novel therapies to modulate NK cells in clinically relevant settings such as intractable viral infections or NK cell lymphoproliferative disorders.

## Disclosures

The authors have no financial conflicts of interest.

## Footnotes

This work was supported by National Institute of Allergy and Infectious Diseases Grants R01 AI078994 and AI073552 and a Basil O’Connor award from the March of Dimes.

The online version of this article contains supplemental material.

Abbreviations used in this article:

- 7-AAD
- 7-aminoactinomycin
- NRMSD
- normalized root mean square deviation.

- Received October 18, 2011.
- Accepted January 18, 2012.

- Copyright © 2012 by The American Association of Immunologists, Inc.