Introduction

Natural Killer (NK) cells are lymphocytes of the innate immune system and provide important protection against viral infections and tumors (Caligiuri, 2008; Vivier et al., 2008). NK cells develop in the bone marrow and express a diverse set of germline-encoded activating and inhibitory receptors (Abel et al., 2018; Biassoni, 2008; Huntington et al., 2007). Mouse cytomegalovirus (MCMV) infection induces antigen-specific expansion of naïve Ly49H+ NK cells in mice, followed by a contraction phase, and then the formation of long-lived adaptive NK cells displaying recall responses to the same antigen (Cerwenka and Lanier, 2016; Robbins et al., 2004; Schlub et al., 2011; Sun et al., 2009). The activating Ly49H receptor, encoded by Klra8, recognizes the MCMV-encoded m157 glycoprotein and leads to activation of Ly49H+ NK cells (Orr et al., 2009; Voigt et al., 2003). In humans, a similar activation of NK cells by human cytomegalovirus is observed with the preferential expansion of NKG2C+ NK cells (Guma et al., 2004; Hendricks et al., 2014; Ishiyama et al., 2022; Kim et al., 2019; Lopez-Verges et al., 2011). The expansion, contraction, and memory formation of Ly49H+ NK cells during MCMV infection bears similarity to antigen-specific responses of adaptive CD8+ T cells that generate memory CD8+ T cells (Adams et al., 2020; Buchholz et al., 2013; Kaech et al., 2002; Schlub et al., 2011).

A range of experiments investigated epigenetic changes, RNA transcription, and surface protein expression in single and bulk Ly49H+ NK cells undergoing expansion and contraction in mice responding to MCMV infection (Flommersfeld et al., 2021; Grassmann et al., 2019; Lau et al., 2018; Martinet et al., 2015; Min-Oo et al., 2014; Nabekura et al., 2014; Potempa et al., 2022; Riggan et al., 2022). These observations show that the NK cell compartment undergoes phenotypic changes, resulting in changes in the abundances of varying subpopulations of cells. The development of these subpopulations depends on the three main signal inputs, antigen simulation of a specific receptor (e.g., Ly49H), co-stimulation, and cytokine signaling (Ni et al., 2012; Peng and Tian, 2017; Romee et al., 2016). These subpopulations express distinct epigenetic, transcriptomic, and protein markers imprinted by signaling and gene regulatory processes. Observations of clonal expansion of barcoded single Ly49H+ NK cells adoptively transferred into Rag2^−/−Il2rg^−/− mice revealed another layer of heterogeneity where sizes and compositions of NK cell clones originating from single naïve Ly49H+ NK cells vary widely (Flommersfeld et al., 2021; Grassmann et al., 2019). This heterogeneity in the NK cell population is regulated by the stochastic kinetics of cell proliferation, cell death, and differentiation involving the subpopulations that emerge as the NK cells respond to the infection. Analysis of similar heterogeneities in clonal CD8+ T cell populations have shown generation of a highly proliferative short-lived effector phenotype and intrinsic stochastic fluctuations in the differentiation and proliferation processes, which account for heterogeneity of the CD8+ T cell clones (Buchholz et al., 2013). These effector CD8+ T cells are CD62L- and differentiate during the expansion phase from upstream memory precursors that are CD62L+(Schlub et al., 2010; Yang et al., 2011).

Despite the qualitative similarities between the antigen-specific responses of NK cells and CD8+ T cells there are several marked differences. Single-cell NK clones typically undergo a relatively smaller 1000× increase in number compared to the 40,000× increase in CD8+ T cell populations (Adams et al., 2020). Additionally, mature CD27-NK cells grow more slowly than their immature CD27+ counterparts (Buchholz et al., 2013; Hayakawa and Smyth, 2006; Paul et al., 2016), whereas effector CD62L-CD8+ T cells proliferate faster than precursor CD62L+ subsets (Adams et al., 2020; Kretschmer et al., 2020). These data suggest there may be mechanistic differences in antigen-specific NK and CD8+ T cell responses.

Here, we addressed whether distinct mechanisms underlie the expansion of NK cells and CD8+ T cells during antigen-specific responses by performing stochastic population dynamic modeling on published bulk and single-cell NK cell data collected during the expansion phase following MCMV infection in C57BL/6 mice. Our analysis and modeling show that the mechanistic rules that describe the expansion of CD8+ T cells are inadequate for describing the development of MCMV-specific NK cells. We found that mature CD27-Ly6C-NK cells undergo a greater rate of cell proliferation compared to immature CD27+ NK cells and mature CD27-Ly6C+ NK cells, and mature CD27-NK cells undergo increased death relative to the CD27+ immature cells. The variations in the timescales of proliferation and death processes coupled with the intrinsic stochastic fluctuations within these processes among the NK cell subsets regulate the heterogeneity in the sizes of the NK cell clones originated from individual immature NK cells during MCMV infection. Quantification of the population kinetics of Ly49H+ NK cells in MCMV infection in the expansion phase shows a 10-fold increase in the cell growth rate compared to that in homeostatic conditions, whereas the effect of the infection on the maturation rate appears to be negligible. We also find that growth rates of adoptive NK cells in transfer experiments during MCMV infection in the expansion phase are higher than those of endogenous NK cells.

Results

A progressive two-state immature to mature differentiation model cannot capture the heterogeneity of the NK cell clones in the expansion phase

We began our investigations by studying whether a minimal two-state progressive differentiation model (Fig. 1a) could describe the size and composition of NK cell clones that originated from single immature Ly49H+ CD27+ NK cells at day 8 post-infection with MCMV, a time typically at the peak of expansion. The differentiation of immature CD27+ to mature CD27-NK cells during MCMV infection and in homeostasis has been widely documented and established in previous experiments (Chiossone et al., 2009; Grassmann et al., 2019; Min-Oo et al., 2014; Mitrovic et al., 2012). We used data collected by Flommersfeld et al.(Flommersfeld et al., 2021), in which they adoptively transferred single genetic marker barcoded Ly49H+ NK cells to a host that is subsequently infected with MCMV. We characterized the NK cell clones by the means and variances of the populations of CD27+ or CD27-NK cells and the co-variance between the two cell types in the clones. The NK cell clones reported by Flommersfeld et al. contained mixtures of CD27+ and CD27-NK cells. We evaluated the percentage of CD27+ NK cells in each clone and computed the correlation (C_size-CD27+) of the size of the clone with the percentage of CD27+ NK cells in the clones. The larger size NK cell clones contained a greater proportion of the CD27-NK cells and thus produced a negative correlation between the clone size and the percentage of CD27+ NK cells in individual clones (Fig. 1b). This finding is in line with that found in a similar experiment by Grassmann et al. (Grassmann et al., 2019). These authors also transferred congenically-marked CD27- and CD27+ Ly49H+ NK cells in 1:1 ratio at day 0 and then assayed the expanded CD27+/-NK cells at day 8 post infection, showing a higher proportion of mature CD27-NK cells compared to the immature CD27+ NK cells (Figure S2E in Grassmann et al. (Grassmann et al., 2019)). This result is supported by other experiments (Hayakawa and Smyth, 2006; Paul et al., 2016), and we used this observation to constrain our model parameter space. Heterogeneity of clones of CD8+ T cells generated by antigenic SIINFEKL stimulation showed a similar pattern where the largest clones in the expansion phase (day 7 post-stimulation) are composed of higher fractions of the mature CD62L-CD8+ T cells (Buchholz et al., 2013). Because mature CD62L-CD8+ T cells proliferate faster than the immature CD62L+ CD8+ T cells (Buchholz et al., 2013; Kretschmer et al., 2020), it is intuitive that a larger size of a clone can be attributed to the larger proportion of the mature CD8+ T cells in the clone (Fig. 1c). In contrast, NK cells proliferate slower when they differentiate to the mature CD27-phenotype. Thus, it is unclear how the similar relationship between the clone size and the composition of the clones can be understood. In order to probe this question further we used a stochastic model (Fig. 1a) that describes the proliferation and differentiation processes in single Ly49H+ NK cells in response to MCMV infection. We do not model the MCMV infection explicitly and assume that the parameters of the model capture the effects of antigen, co-receptor, and cytokine stimulation and signaling and epigenetic modifications on the rates of differentiation and the proliferation of the cell types. In this model, immature CD27+ NK cells differentiate with a rate r to mature CD27-NK cells, and each cell type proliferates with rates, k_I and k_M, respectively. The proliferation and differentiation processes are modeled as stochastic events occurring with the above rates. We start our simulations with a single CD27+ NK cell and evaluated the number of CD27+ and CD27-NK cells at 8 days, which represent the size and the composition of a single NK cell clone. Repetitions of the stochastic simulations with the same initial configuration and rates will produce clones of different sizes and compositions due to intrinsic noise fluctuations in the stochastic simulations. We computed the correlation (C_size-CD27+) between the clone sizes and the percentages of the CD27+ NK cells in each clone from our simulations. Our parameter scan of randomly sampled values of k_I, k_M, and r revealed that C_size-CD27+ is negative only when k_I < k_M (Fig. 1d). However, the observations that CD27+ cells should grow faster than CD27-cells imply k_I > k_M for the NK cells. Therefore, the two-state model cannot capture the heterogeneity of the NK cell clones.

Figure 1:

Next, we investigated if inclusion of death of NK cells can make the model consistent with the negative correlation between the clone size and the percentage of CD27+ NK cells in each clone. We updated our proposed model to include both proliferation and death of each cell population (Fig. 1e). In such a model, each cell has a probability of proliferating according to birth rate b and probability of dying indicated by the death rate d. Cells can go through many proliferation events before undergoing death, which would indicate a positive growth rate. The net growth of the CD27+ and CD27-NK cell populations are determined by k_I=b_I-d_I and k_M=b_M-d_M, respectively, where b_I and d_I (or b_M and d_M) are the proliferation and death rates of the immature CD27+ (or mature CD27-) NK cells. Repeating a parameter scan by randomly sampling the 5 parameters (b_I, d_I, b_M, d_M, and r) from uniform distributions and evaluating the correlation C_size-CD27+ shows that some parameter configurations can generate the negative values at the parameter values, k_I > k_M (Fig. 1f). Further inspection shows that parameter values in the regime k_I > k_M that produce the negative correlation share a common trait: mature CD27-cells proliferate and die faster than immature CD27+ cells. To gain further insight into this behavior we investigated the dependence of C_size-CD27+ on the parameter values and stochastic fluctuations analytically and found that larger proliferation and death rates of the CD27-NK cells produce large stochastic fluctuations in the sizes of the CD27-cell populations. Thus, mature CD27-NK cells in a clone could undergo many proliferation events before encountering death, creating larger clones with higher proportion of mature CD27-NK cells. Conversely, the CD27-cells could go through many death events before a proliferation event causing the clone size to shrink (details in Supplementary Text 1).

We investigated whether this model could accurately fit the first (means) and second (variances and covariance) moments of the measured CD27+ and CD27-NK cell populations at 8 days post-infection. We fit the model by minimizing the sum of square errors (SSE) subject to the constraints C_size-CD27+ < −0.2 and k_I > k_M. This parameter configuration cannot accurately capture the moments of the clonal burst data (Fig. 1g). Thus, although the model that incorporates both NK cell proliferation and death is able produce the negative correlation between the clone size and the percentage of CD27+ NK cells in each clone, it poorly describes the distributions of the immature and mature NK cells in the clones.

To further investigate the role of NK cell death in explaining the clone size distribution we compared the above two-state models without or with NK cell death based on the theory of model selection using modified Akkaike Information Criterion (AIC_c) (Burnham et al., 2011). We computed the probability distribution function of the number of immature CD27+ (N_I) and mature CD27-(N_M) cells at any time obtained from our stochastic model by solving the corresponding Master equation numerically and computed the Likelihood function for the available clone size distribution data (details in the Materials and Methods section), and then evaluated modified-AIC (AIC_c) for the two models. To keep the computation in a feasible range we used a part (38 out of 56 clone sizes) of the available data where the maximum clone size is ~ 1000. The subset of clone size data used here also shows a negative value of C_size-CD27+. The AICc value for the two-state model with death is over 50 points lower than that for the two state model without death, which clearly show that the data favor the model with NK cell death.

Although this approach only models a subset of the data, it also finds that birth and death rates for mature cells are greater than that of immature cells, and that the growth rate of immature cells is greater than that of mature cells. Further details regarding the analysis are provided in the Materials and Methods section and the Supplementary Material. Because of the computational intensive nature of this approach and its agreement with the simpler method of fitting moments with a constraint, we continue with the later approach for subsequent models.

A progressive three-stage differentiation model can quantitatively recreate clonal burst dynamics

We modified our two-state model to include an intermediate state (Fig. 2a) in the model. This reflects many investigations probing NK cell maturation, which show NK cell development goes through several maturation states and some of those states can potentially represent the proposed intermediate state in the model. A well-established NK cell maturation scheme is CD27+CD11b-→CD27+CD11b+ →CD27-CD11b+ (Chiossone et al., 2009). However, the clonal burst data show a negligible proportion of CD11b-NK cells at 8 days post-infection. Another potential maturation scheme could be CD27+KLRG1-→CD27+KLRG1+ →CD27-KLRG1+ (Kamimura and Lanier, 2015). In addition, recent single cell RNA-seq experiments show the presence of CD27-Ly6C+ and CD27-Ly6C-cells at 8 days post-MCMV infection (Riggan et al., 2022), thus the intermediate state could represent a CD27-Ly6C-state as well. Here we present results for the model where the intermediate state represents a mature CD27-state, the other variations of the model are shown in the Supplementary Material (Fig. S3-S5). In order to determine parameter ranges that give rise to higher growth of immature NK cells compared to their immature counterparts, we computed the average size of a clone that originates from either a single CD27+ NK cell or single CD27-NK cell by solving ODEs describing the kinetics of the average NK cell populations. We constructed an effective growth difference metric Δ_k, which describes the difference between the average NK cell populations at day 8 that originate from a single CD27+ or CD27-NK cell at t=0. If Δ_k is positive, then CD27+ cells grow faster than CD27-cells (Materials and Methods). We performed a parameter scan on this three-stage model and find that some parameter configurations with Δ_k > 0 can capture the negative values of C_size-CD27+ (Fig. 2b). We then evaluated whether this three-stage model can accurately describe the distributions of the immature and mature NK cells in clones. We fitted the means, variances, and the covariance of the NK subpopulations in the clones by varying the model parameters with constraints Δ_k > 0 and C_size-CD27+ < −0.2 and found that the best-fit three-stage model shows an excellent fit to the data (Fig. 2c, Table 1). Thus, adding an intermediate CD27-stage to the model can capture all the key features of the heterogeneity in the NK cell clonal populations.

Table 1:

Figure 2:

The three-stage model above makes no assumptions about the identities of the differing phenotypes of the CD27-cell subsets. However, the best fit parameters indicate that the first CD27-subset has a high proliferation rate, and the final CD27-subset has a high death rate. Riggan et al. (Riggan et al., 2022) previously determined that, in the contraction phase, Ly6C+ NK cells differentiate from Ly6C-cells, and that Ly6C+ cells die more rapidly than Ly6C-cells. Our model shows that similar activity may be observed during the expansion phase as well, as we show higher death of Ly6C+ cells and a unidirectional maturation from Ly6C-→ Ly6C+ cells. Because our analysis pertains to the expansion phase of the NK cell response to MCMV we thus aimed to validate if the mature states proposed in the model can be identified with CD27-Ly6C- and CD27-Ly6C+ NK cells. We found that our best fit parameters correctly recreate the experimentally observed strong negative correlation between the %CD27+ and %Ly6C+ NK cells across clones (Grassmann et al., 2019), and correctly estimate that most clones have <20% Ly6C+ cells at day 8 post-infection (Fig. 2d-e). This independent validation gives evidence that our parameter estimates can be meaningfully interpreted for a model where NK cells differentiate along a trajectory from CD27+Ly6C-→ CD27-Ly6C-→ CD27-Ly6C+.

We sought to understand the mechanism by which the three-stage model can generate negative values for the correlation C_size-CD27+ as well as quantitatively describe the heterogeneity in the NK cell populations in the clones. The two-stage model predictions for the first and the second moments for the NK cell populations are substantially less than the actual values, indicating that our best-fit parameters are not producing enough NK cells (Fig. 1g). The presence of the intermediate state in the three-stage model alleviates this issue, which can rapidly grow and increase the number of CD27-Ly6C+ NK cells while the rapidly dying CD27-Ly6C-NK cells can help meet the constraint that a population of CD27-cells grows slower than one of CD27+ cells. This variation increases stochastic fluctuations in the observed numbers of CD27-cells, as some clones will have many CD27-cells from high growth of the intermediate Ly6C-subset, while other clones will have few CD27-cells from low growth of the Ly6C+ subset.

Furthermore, the presence of the CD27-Ly6C+ NK cells expand the range of the model parameters that can give rise to negative values of C_size-CD27+.

Comparison of kinetics between endogenous NK cells, adoptively transferred cells, and homeostatic cells identifies potential identical maturation rates in all conditions

Here we compare the kinetics of NK cell maturation occurring in response to MCMV infection and in homeostasis. In particular, we investigated the changes in the rates of cell growth and differentiation induced by MCMV infection. First, we quantified the growth and the differentiation rates using time-stamped flow cytometry data on homeostatic NK cells with a one-time treatment of tamoxifen in fluorescent NK cell reporter mice reported by Adams et al.(Adams et al., 2021). Treating the mice with tamoxifen determines whether observed NK cells are newly populated from bone marrow or existed in peripheral blood at the onset of the infection. We fitted a two-state CD27+ → CD27-maturation model described by ordinary differential equations (see Materials and Methods) to the data collected at days 3, 7, 14, 23, 29, and 35 to estimate the growth rates of the immature CD27+ (k_CD27+) and the mature CD27-(k_CD27-) NK cells, the rate of differentiation (r) of the CD27+ to CD27-NK cells, and the rate (λ) with which CD27+ NK cells are populated from bone marrow. We estimated these rates both for Ly49H+ and Ly49H-NK cells. The estimated rates for Ly49H+ cells (Table 2, Fig. 3a) show that a CD27+ cell population would double every 6.9 days if there were no differentiation, whereas CD27-cells have a negative growth rate with a half-life of 5.6 days. This implies that the death rate of a population of CD27-cells is greater than their proliferation rate. The differentiation rate from CD27+ → CD27-phenotypes is similar to that of the death rate of Ly49H+ CD27-cells. Approximately 1.37% of the Ly49H+ NK pool is replenished from bone marrow every day. Comparing these rates to those of the Ly49H-NK pool reveals that Ly49H-NK cells are repopulated from bone marrow at almost twice the rate (~2.4%/day). All other estimated rates for Ly49H-NK cells are within the confidence bounds of the Ly49H+ estimated rates (Table 2, Fig. 3b). Given the constraint governing λ which dictates

Table 2:

Figure 3:

(Materials and Methods) the difference in estimated λ values can be attributed to a higher percentage of CD27-cells. Next, we compared these rates with corresponding rates estimated in our three-state model using the clonal burst data. The estimate for CD27+ → CD27-differentiation rate is equivalent between Ly49H+ NK cells in homeostasis and in MCMV infection. Our estimations for the clonal burst data (Table 1) show that varying subsets of CD27-cells have varying growth rates ranging −0.7 day⁻¹ to 1.3 day⁻¹in MCMV infection, a wide range not in disagreement with the estimate for homeostatic data. Ly49H+ CD27+ cells grow ~10× faster in MCMV infection than in homeostatic conditions.

We investigated whether the kinetics of NK cell maturation in the expansion phase during MCMV infection in single-cell adoptive transfer experiments is quantitatively similar to that of an endogenous population of NK cells. To this end, we collected time-stamped mass cytometry data on C57BL/6 mice splenocytes at days 0, 4, and 7 post-MCMV infection (Potempa et al., 2022). We then gated Ly49H+ NK cells from each mouse and determined the mean relative abundances across mice of CD27+Ly6C-, CD27-Ly6C-, and CD27-Ly6C+ subsets among the NK cell population (Fig. S6). We multiplied these relative abundances by mean absolute numbers of Ly49H+ endogenous NK cells during MCMV infection from Robbins et al.(Robbins et al., 2004) to obtain time-stamped population sizes of NK cell subsets. We then modeled these data with the ordinary differential equation (ODE) solution in Equations S17-S25 using the parameters from Table 1 and an initial condition of the mean uninfected in vivo abundances of endogenous NK cells. The ODE result is multiple orders of magnitude above the measured abundances (Fig. S7a). This indicates that the kinetics of endogenous NK cells differ significantly from that of NK cells in adoptive transfer experiments.

To further quantify the differences in the kinetic rates between the adoptively transferred and endogenous NK cells, we then fit the 3-stage model with 2 mature stages to the endogenous cell population. We used the gated Ly49H+ populations of CD27+/Ly6C-, CD27-/Ly6C-, and CD27-/Ly6C+ NK cells and fit to a system of ODEs that describe the net growth rates of each population and the differentiation rates between them (Materials and Methods). The best fit to this data has unrealistic estimates for the death rate of CD27-/Ly6C+ NK cells and the differentiation rate from CD27+ → CD27-cells (Table S5, Fig. S7b). We then attempted to constrain the fitting procedure by setting the CD27+ → CD27-differentiation rate to those of the previous estimations for the clonal burst kinetics and homeostatic conditions (Tables 1 and 2).

Under this constraint, the only rate estimated to be within the confidence interval for the equivalent rate in adoptive transfer experiments is the differentiation rate from Ly6C-→ Ly6C+ (Table 3, Fig. 3c). The growth rate of CD27+Ly6C-cells is ~4× faster for adoptively transferred cells, the growth rate of CD27-Ly6C-cells is 2-3× faster, and the death rate of mature cells is ~1.5× faster. However, the resulting best fit to the endogenous NK cell population does show that the intermediate CD27-/Ly6C-cell population has the fastest growth rate, and the mature CD27-/Ly6C+ cells have a large death rate, confirming the qualitative results from the adoptive transfer experiments. Thus, although the growth kinetics are different between endogenous and adoptive transfer experiments, both systems have the intermediate population proliferating rapidly and the terminally mature population dying rapidly.

Table 3:

Both the adoptively transferred cells and endogenous population are estimated to have a large death rate for CD27-cells. To test this result, we analyzed NK cell viability in mice that were infected with MCMV at different times post-infection (Fig. 4a). Surprisingly, mature (CD27-) cells indeed had a greater fraction of dead cells than immature (CD27+) cells during expansion (days 4 and 7 post-infection) (Fig. 4b). To determine if the higher rates of cell death in the CD27-NK cell population in the three-state model give rise to a higher fraction of dead NK cells we extended our three-state model to evaluate the number of dead mature and immature NK cells (Fig. S8c) which showed that the best fit parameters for the adoptive transfer experiments are consistent with the results reported above (Fig. 4c). Because we only have best fit growth rates of the endogenous NK cells and not separate proliferation and death rates, we tested a variety of birth and death rates consistent with the best fit parameters such that b-d=k for each cell subset.

Figure 4:

Many of the possible configurations agree with the higher percentage of live CD27+ cells. This gives additional confirmation that CD27-cells may in fact be dying at faster rates than CD27+ cells.

Discussion

We employed stochastic minimal mechanistic in silico models to investigate mechanisms that underlie development of MCMV-specific NK cells at early times (0 – 8 days post-infection). We found Ly49H+ NK cells undergoing antigen-specific expansion follow a maturation scheme from (immature) CD27+Ly6C-→ (mature I) CD27-Ly6C-→ (mature II) CD27-Ly6C+, characterized by a highly proliferative mature CD27-Ly6C-phenotype and a more mature CD27-Ly6C+ phenotype with higher rates of cell death (Fig. 5). The differences in the rates of proliferation and cell death are crucial to give rise to the heterogeneous clones that originate from single immature CD27+ NK cells, in particular, the presence of larger sized NK cell clones composed of mature CD27-NK cells at 8 days post-infection. Similar heterogeneity in the clones of CD8+ T cells generated in the expansion phase (<8 days post-infection) from an acute antigen-specific expansion of single immature CD62L+ cells in C57BL/6 mice has been observed (Buchholz et al., 2013); however, the larger CD8+ T cell clones composed of effector CD62L-arise from the higher rate of proliferation of the CD62L-CD8+ T cells. This is in contrast with the NK cell maturation program where the larger sized clones of mature NK cell populations arise due to the interplay of uneven increase in the proliferation rates and increase in the rates of cell death as the NK cells mature, and intrinsic stochastic fluctuations occurring in these processes. Innate heterogeneity in adoptively transferred single NK cell phenotypes may also contribute to observed stochasticity. In CD8+ T cells, precursor CD62L+ cells are deemed memory precursors, matching their lowly proliferative phenotype (Kretschmer et al., 2020). An equivalent population of memory precursor NK cells has not been identified. CD27-Ly6C+ cells match previously defined phenotypes for memory NK cells and could potentially be memory precursors. However, our model shows that these cells experience rapid cell death. Therefore, for these cells to live longer to become memory cells additional processes such as a change in environmental conditions or an epigenetic change that promotes survival in these cells at later times during the contraction phase could become relevant. This differentiation trajectory is also in contrast with CD8+ T cells, as NK cell memory precursors would be formed downstream of effector (Ly6C-) NK cells, whereas CD8+ effectors (CD62L-) are formed downstream of CD8+ T cell memory precursors (CD62L+) in acute infection. Alternately, another mature NK subset such as Ly6C-KLRG1-that is lowly proliferative and long-lived could give rise to the adaptive NK cell phenotype. Inclusion of these processes as well as the later time kinetics could be a future extension of our approach. We note a similarity of our three state NK cell maturation model with maturation of T cells in persistent (13-140 days) T. gondi infection in genetically resistant mice reported by Robey and colleagues (Chu et al., 2016) where the infection gave rise to a continuous production of CXCR3-KLRG1+ terminally differentiated effector T cells from a quiescent memory CXCR3+KLRG1-T cells via an intermediate CXCR3+KLRG1+ state. The intermediate T cell state showed higher proliferation and increased death compared to the memory and effector T cells. However, these T cells were uniformly CD62L-in contrast to the memory and effector T cells that are produced in response to acute antigenic SIINFEKL infection.

Figure 5:

NK cell clones in MCMV infection go through a smaller (~1000×) fold increase in the expansion phase compared that (~40,000×) of antigen-specific expansion of CD8+ T cells. This could be because CD8+ T cells can generate autocrine cytokines such as IL-2 that help the cells to proliferate (Kalia and Sarkar, 2018), whereas NK cells rely on other cells to make the necessary cytokines for proliferation (such as IL-2 and IL-15) and survival (such IL-21) (Fauriat et al., 2010). Furthermore, NK cells are more dependent on cytokines such as IL-15 compared to the CD8+ T cells (DeGottardi et al., 2016). Therefore, cytokine competition for IL-15 and IL-2 might also explain increased death of mature NK cells in the expansion phase. A test of this hypothesis could be investigation of MCMV-induced NK cell maturation in Il15+/- mice.

We quantified the population kinetics of antigen-specific Ly49H+ NK cells during the expansion phase of MCMV infection (<8 days post-infection) using both endogenous and adoptively transferred NK cells. This antigen-specific response was further compared with the kinetics of NK cell maturation in the absence of MCMV infection. Adoptively transferred Ly49H+ NK cells grow at an increased rate compared to the endogenous Ly49H+ NK cells (~3× faster for Ly6C-cells). A possible explanation of this phenomenon is the stronger competition for resources (e.g. antigens and/or cytokines) that would occur in endogenous NK cells could prevent them from growing as quickly as the less dense population of adoptively transferred NK cells. Interestingly, our estimate of the differentiation rate of the immature CD27+ to mature CD27-NK cells showed negligible effect of MCMV infection on this process. This shows that the signaling processes governing NK cell antigen-specific expansion and cytotoxicity are not impacting transcription and translation of CD27, but rather preferential expansion of certain subsets is creating altered phenotype abundances observed in MCMV infection.

Our stochastic modeling of heterogeneous NK cell clones showed that larger stochastic variations in the mature NK cell populations arising due to higher intrinsic noise fluctuations generated by the interplay between the higher proliferation and cell death processes are important for producing large clones of mature NK cells at 8 days post-infection. This is in contrast with stochastic kinetics of CD8+ T cells where the larger clones of effector CD62L-CD8+ T cells can arise due to a higher rate of proliferation of these cells even in the absence of any cell death (Buchholz et al., 2013). Our effort to use different models to quantitatively describe the heterogeneity of NK cell clones demonstrate the critical role of the correlation between the clone size and the proportion of the mature NK cells in the clones, C_size-CD27+, in screening out mechanisms that are unable to give rise to the observed heterogeneity.

Materials and methods

NK cell development model

We modeled the NK cell proliferation, death, and maturation as a continuous time Markov processes where a proliferation or a death event in a NK cell population increases or decreases the population number, respectively, by one cell. A differentiation event changing a NK cell from immature (I) to mature (M) decreases the NK cell population in state I and increases the population in state M simultaneously by a single NK cell. The kinetics of the system are given by the Master equation describing the time evolution of the conditional probability for the system to be in a particular state containing different numbers of immature and mature NK cells at any time t given a specific state at t=0. The ODEs describing the kinetics of first and second moments can be derived from the Master equation. The solutions are shown in the Supplementary Materials. The first and second moments for the clone size distributions for all our models were calculated by solving the corresponding ODEs.

Stochastic simulations

We use the Gillespie algorithm(Gillespie, 1977) for all our stochastic simulations. For two-stage and three-stage models, we initialized the simulation with a single CD27+ NK cell. We evolved the system until 8 days and numbers of CD27+ and CD27-cells in the NK cell clones were recorded. We calculated the correlation C_size-CD27+ for the clones by averaging over 10⁴ stochastic trajectories.

Fitting model to moments

The first and second moments for the clone size distributions obtained for two- and three-stage models were computed by solving ODEs (Eqs. S2-S6 in Supplementary Text 1 and the ODEs in Supplementary Text 2). The parameter vector represents the rates for proliferation, differentiation, and death in the models. To compare the predicted moments from the model with the actual moments y computed from the data, we used the cost function below:

where, σ_i is the standard deviation corresponding to the data used for computing y_i. A penalty function, is introduced to implement the constrained imposed by sign of the correlation C_size-CD27+. We set to a high value (= 10⁴) when C_size-CD27+ > −0.2 at any value of Δ_k or C_size-CD27+ < 0 and Δ_k < 0, otherwise .This constraint describes the experimental observation that immature CD27-NK cells grow faster than mature CD27+ NK cells in the expansion phase. σ_i is calculated by taking the standard deviation of each moment across 10,000 bootstrapped samples of the clones. Minimization of the cost function was performed with the lmfit python package, using the ‘least_squares’ Trust Region Reflective algorithm. Parameter search space was confined from 0.0 to 2.0 day⁻¹ for all parameters.

Effective Growth Rate for 3-stage system

We developed a metric to quantify the difference in the growth rate of immature CD27- and mature CD27+ NK cells for the three-stage model. We computed mean populations of the NK cell subsets at day 8 by solving Eqs. S17-S25 in Supplementary Text 2 for two different initial conditions. In one case we initialized with the immature CD27+Ly6C-cells, and in the other case we initialized with equal proportions of mature CD27-Ly6C- and mature CD27-Ly6C+ cells. We then computed the “effective growth difference” (Δ_k) by evaluating the difference in the mean populations of the immature CD27+ and the mature CD27-NK cells at day 8. We performed an asinh(Δ_k) transformation to help visualization of the data in the graphs shown in the manuscript.

Parameter scans

For parameter scans showing the relationship between growth rates and the correlation between clone size and %CD27+ cells, each parameter is randomly selected from a uniform distribution from 0 to 1.6 day⁻¹. For each parameter set, 10⁴ Gillespie simulations were carried out to compute the correlation C_size-CD27+ for the NK cell clones at day 8. We performed scans for 10⁵ different parameter sets.

Ordinary differential equation fits to homeostatic data

The growth and differentiation kinetics of populations of CD27+Tom+ (n₁), CD27-Tom+ (n₂), CD27+Tom- (n₃), and CD27-Tom- (n₄) cells in the periphery in homeostasis are given by the ODEs below:

The CD27+ and CD27-cells grow with rates k_CD27+ and k_CD27-, respectively, and the differentiation CD27+ → CD27-occurs at a rate, r. The influx of immature CD27-cells from the bone marrow to the periphery occurs at a rate α. At homeostasis it is reasonable to assume that the sum of all Ly49H+ and the sum of all Ly49H- NK cells were assumed to be constant through time, i.e.,

The kinetics of the percentages (N_i = 100×(n_i/n_TOT)) of each cell type follow the ODEs below:

We define , which describes the influx rate of either the total Ly49H+ or total Ly49H-population as a percentage of the given population. This makes the ODEs which we fit:

We assume the total number of NK cells and the relative abundances of varying subsets do not change with time in homeostasis, therefore,

These constraints lead to equalities between parameters,

Thus, by estimating only k_CD27+ and k_CD27-, and using the mean relative abundances of CD27+ and CD27-cells across mice and time, we can solve all four rate parameters in the homeostatic state.

Ordinary differential equation fits to mass cytometry data

We gated Ly49H+ NK cells in the CyTOF studies into three subsets: CD27+Ly6C-, CD27-Ly6C-, and CD27-Ly6C+. Details of the gating are described in the Supplementary Material. We describe the growth and differentiation kinetics of the populations the immature CD27+Ly6C-(N_I), mature CD27-Ly6C- (N_INT), and mature CD27-Ly6C+ (N_M) by the ODEs below,

Where, k_I, k_INT, and k_M denote the rates of growth for N_I, N_INT, and N_M, respectively. The rates of differentiation for CD27+Ly6C-→ CD27-Ly6C- and CD27-Ly6C-→ CD27-Ly6C+ are given by r₁ and r₂, respectively. Since the above system of ODEs is autonomous and linear, we can obtain the time dependent solutions analytically and exactly. The analytical solution was numerically evaluated for different values of the rate parameters represented by .The analytical solutions were evaluated numerically for specific rate constants and initial conditions. We used the estimated population sizes (details in the main text) for the above NK cell subsets in the CyTOF data at t=0 as the initial conditions and computed the population sizes at day 4 and day 7 by solving the above ODEs for a parameter set .We then estimated the parameters by by fitting the population sizes at day 4 and day 7 using the cost function below.

Bootstrapping confidence intervals for parameter estimates

To assess the confidence intervals for parameter estimates, we randomly sampled the clones with replacement with equal probability, calculated each moment and its standard deviation, and fit the moments as described above. We repeated the process 1,000 times, yielding 1,000 estimates for the parameters θ. We then sorted the parameter estimates by magnitude and took the 2.5^th and 97.5^th percentile for each. These estimates made up the 95% confidence intervals reported here.

For estimations of the homeostatic and endogenous NK populations, we sampled a population of mice before sampling errors in measurements for each bootstrap estimate.

Animals and MCMV infection

C57BL/6 (B6) B6 mice were purchased from the National Cancer Institute at 6 weeks of age and housed in the specific pathogen-free animal facility at the University of California San Francisco (UCSF) in accordance with the guidelines of the Institutional Animal Care and Use Committee. Experiments used mice of both genders, aged 8-12 weeks.

B6 mice were infected with 1000 plaque-forming units (PFU) of Smith strain MCMV via intraperitoneal injection. Spleens were acquired different time point following MCMV infection and examined by flow cytometry.

Flow Cytometry

Splenocytes from uninfected or infected mice were processed into single cell suspensions, RBC lysed, and then stained for flow cytometry. Surface staining was performed using pre-conjugated antibodies in the presence of the viability marker ZombieRed (BioLegend) and were then analyzed on an LSR II cytometer (BD Biosciences). The antibodies used for flow cytometry were: CD3χ (145-2C11), CD19 (6D5), CD49b (DX5), NK1.1 (PK136), KLRG1 (2F1), Ly6C (HK1.4), CD27 (LG.3A10), and Ly49H (3D10). All antibodies were purchased from BioLegend, BD Biosciences, or ThermoFisher Scientific. An unconjugated anti-CD16/CD32 antibody (clone 2.4G2, UCSF Antibody Core) was used to block non-specific binding.

Mass cytometry

Samples were prepared for mass cytometry as previously described (Potempa et al., 2022). Briefly, splenocytes from each mouse were obtained by mechanical dissociation through a 40-μm filter, treated with ACK lysis buffer, and resuspended in PBS + 5 mM EDTA (PBS/EDTA). Cells were stained for viability by mixing 1:1 with PBS/EDTA + 100 mM cisplatin (Enzo Life Sciences) for 1 minute before quenching with PBS/EDTA + 0.5% BSA (Staining Media). One million splenocytes from each mouse was stained with metal-conjugated antibodies against CD16 (Nd144) and CD32 (Eu153) for 15 minutes at room temperature in Staining Media before supplementation with an Fc blocking Ab (clone 2.4G2). After washing in Staining Media, cells were stained with a master mix of 35 additional metal-conjugated antibodies (Table S6) for 30 minutes at room temperature. At the conclusion of the staining period, cells were washed in PBS/EDTA before fixation with PBS/EDTA + 3.2% paraformaldehyde and stored at 4º C for 3-5 days. One day prior to data acquisition, samples were barcoded using a 20-plex Pd barcoding kit (Fluidigm) per manufacturer protocol, pooled, and fixed again, overnight at 4º C, in a PBS/EDTA + 3.2% paraformaldehyde solution supplemented with a 1:1000 dilution of a 191/193 DNA intercalator (Fluidigm). Just prior to acquisition, barcoded samples were pelleted, washed in double-deionized water, and resuspended in double-deionized water with a 1:20 dilution of EQ four Element Calibration Beads (Fluidigm). Data were acquired using a CyTOF2-Helios mass cytometer (Fluidigm), which were then normalized and debarcoded prior to analysis.

Acknowledgements

This work was supported by NIH grant R01AI146581 to J.D. N.M.A. was supported by a Medical Scientist Training Program grant from the National Institute of General Medical Sciences (T32GM007739 to the Weill Cornell/Rockefeller/Sloan Kettering Tri-Institutional MD-PhD Program) and by an F30 Predoctoral Fellowship from the National Institute of Allergy and Infectious Diseases (F30 AI136239). We would like to thank the Ohio Supercomputer Center for computational resources.

Additional files

Supplement