Clonal stochasticity in early NK cell response to mouse cytomegalovirus is generated by mature subsets of varying proliferative ability
- Steve and Cindy Rasmussen Institute for Genomic Medicine, Nationwide Children’s Hospital, United States
- Biomedical Sciences Graduate Program, The Ohio State University, United States
- Department of Microbiology and Immunology and the Parker Institute for Cancer Immunotherapy, University of California, San Francisco, United States
- Department of Physics, The Ohio State University, United States
- Immunology Program, Memorial Sloan Kettering Cancer Center, United States
- GIG Statistical Consulting, LLC, United States
- Department of Pathology, New York University Grossman School of Medicine, United States
- Department of Pediatrics, The Ohio State University, United States
- Pelotonia Institute for Immuno-Oncology, College of Medicine, The Ohio State University, United States
Figures
Figure 1 with 4 supplements
Two-stage linear models cannot capture desired properties of experimental clonal burst data.
(a) Schematic representation of the two-stage model. The model describes immature (CD27+) cells and mature (CD27-) cells. Each cell type has a distinct growth rate (kI for immature, kM for mature), and immature cells differentiate to mature cells at rate r. (b) Experimental data from Flommersfeld et al. display negative value for the correlation Csize-CD27+ between clonal burst size and %CD27+ cells. (c) Schematic description of antigen-specific clonal proliferation and differentiation in CD8+ T cells. Because mature CD62L- cells proliferate more rapidly compared to the immature CD62L+ cells, initial clones made primarily of CD62L- cells can grow to generate the larger size clones. (d) A parameter scan of all possible combinations of kI, kM, and r shows that negative correlations (Csize-CD27+<0) are not realized when kI > kM. Each point represents a value of Csize-CD27+ obtained from the model at a unique value of set of model parameters. The parameter configurations that populate the bottom-right quadrant satisfy our two constraints imposed by experimental observations regarding the negative values of Csize-CD27+ and higher growth rate of immature NK cells compared to their mature counterparts. (e) Schematic representation of the three-stage model including cell death. The model has separate parameters for birth of immature (CD27+) cells (bI), death of immature cells (dI), birth of mature (CD27-) cells (bM), death of mature cells (dM), and differentiation of immature cells to mature (r). In this model, the growth rates are defined as kI = bI dI and kM = bM dM. These net rates relate to the rates of growth or loss of clonal populations. (f) A parameter scan for this model shows that some parameter configurations can meet the two constraints where kI > kM and the parameters result in negative values for Csize-CD27+. Each point represents a unique parameter configuration. Red points denote parameter configurations where bM > bI and dM > dI. (g) Shows the best fit of the model in (e) to the moments with the constraints kI > kM and Csize-CD27+<–0.2. Circles represent the fitted values of the first and second moments for the distribution of the NK cell clones (n = 55) obtained from the model. Blue and red circles represent mean populations of CD27+, and CD27- NK cells, respectively. The variances of the CD27+ and CD27- NK cells are shown in cyan and pink circles, and the covariance between the CD27+ and CD27- NK cells is shown in purple. Horizontal lines around the circles show the standard deviations of bootstrapped moments (10,000 bootstraps). The solid line is the line y=x, which demonstrates where values would lie if we had a perfect fit.
Figure 1—figure supplement 1
Analysis of the stochastic kinetics using the Master equation.
The Master equation and its solution are described in Appendix 1. (a) The approximation for the correlation Csize-CD27+ maintains the sign of the actual correlation. The approximation for Csize-CD27+ (y-axis) was compared with the empirical solution (x-axis) for several parameter sets. Each parameter set is represented by a red star. The black dashed line represents the line y = x. For each parameter set, empirical Csize-CD27+ was calculated by simulating 1000 clones. (b) Histogram of counts of CD27+ and CD27- NK cells shown for the 38 data points used in estimating the model parameters. (c) (left) Comparison of the numerical solution of the Master equation at t = 8 days with upper bounds in NI for a pure birth process (NI → NI +1, rate bI) (dashed red) with the exact solution (solid black) with no upper bound in NI. The system starts at t = 0 with NI = 1; bI is taken as the largest value taken in our parameter range. (Right) shows the absolute difference in the probability distribution function of the exact solution and numerical solution. The largest difference is of the order of 10–9. (d) shows changes in the -log-likelihood with parameter iterations for the model without death (top) or with death (bottom). The minimum log-likelihood values correspond to –545.4 (top) and –514.6 (bottom). The optimal parameter values are bI = 0.97, bM = 1.07, r = 0.05 for two states without death, and bI = 0.97, bM = 1.7, dI = 0.15, dM = 0.92 and r = 0.05 for two states with death.
Figure 1—figure supplement 2
Stochasticity in the duration of NK cell activation helps contribute to negative correlations.
(a) The model in (a) differs from the original two-stage model in that new CD27+ cells are introduced at varying times during the course of infection. A in the model clone is stimulated for a random length of time between 0–8 days drawn from a lognormal distribution. (b) Scanning different parameter values in the model shows a range of parameters (shown within the box) can produce negative values in Csize-CD27size-CD27+ for rates kI > kM. The correlation was calculated for 104 clones and could recreate negative correlations as well. (c) shows the variation of the clone sizes with the percentage of CD27+ cells in each clone for the rates taken from Appendix 3—table 1 but with growth rates in place of separate birth and death processes, where kI = 0.533, kM = 0.477, and r=0.135. The correlation Csize-CD27size-CD27+ was calculated for 56 clones (same number of clones as the experimental data), which shows a statistically insignificant (p=0.16) negative value. (d) shows additional clones for the same parameter set as in (c) at different times to demonstrate how stochasticity in total activation time contributes to negative correlations. Color represents the length of activation, so clones that underwent division and differentiation for only one day are shown in blue.
Figure 1—figure supplement 3
Linearly increasing differentiation rate with number of divisions results in negative correlations for a two-stage model.
(a) A schematic representation of our model. The clonal population is initiated by a single founder cell. The daughter cells differ from the mother cell but inherit some attributes. In our case, at each cell division the rate of differentiation in the daughter cells increases linearly from that of the mother cell. In this representation, the rate of differentiation after the first division, r0=0.0675 and the rate of differentiation in the daughter cells after ndivisions divisions is given by, rcell=r0×ndivisions. (b) Variation of the clone sizes with the percentage of CD27+ cells in the clones for 56 clones generated in our simulation of the model is shown (a). Parameters are taken from Appendix 3—table 1 but with growth rates in place of separate birth and death processes, where kI = 0.533, kM = 0.477, and r0=0.0675. The growth rate of CD27+ cells was higher than that of CD27- cells, yet the correlation is –0.532.
Figure 1—figure supplement 4
Modeling asymmetric division of CD27+ cells in a two-stage model.
(a) Schematic representation of the two-stage model with asymmetric division. It is altered from the model in Figure 1e. So that it has two different terms for division of immature cells. Immature (CD27+) cells can divide symmetrically as before according to rate bsym, and additionally can divide asymmetrically into a single immature and single mature (CD27-) cell according to rate basym. In this model, the growth rates are defined as kI = bsym + basym– dI and kM = bM – dM. (b) A parameter scan for this model where dI = dM = 0 shows that adding asymmetric division can also account for negative correlations. Each spot represents a unique parameter configuration, with color indicating the value of basym. (c) Best fit of this model to the moments with constraints Δk > 0 and Csize-CD27+ < –0.2.
Figure 2 with 3 supplements
A three-stage linear model can capture experimental observations of NK cell clones if intermediate double-positive NK cells are the fastest growing subset.
(a) Schematic representation of the three-stage model. In this model there are three stages of NK cell maturation: an immature (CD27+X-), intermediate (CD27-X-), and terminally mature (CD27-X+) subset, where the terminally mature state is marked by an expression of a yet to be determined marker protein X. Each subset has a distinct proliferation and death rate, and cells progress from immature to intermediate maturity according to rate r1 and from intermediate to terminal maturity according to rate r2. (b) A parameter scan for the model in (a) shows that adding a third stage of maturation can account for the negative values in Csize-CD27+. Each point represents a unique parameter configuration. Positive values of Δk indicate that immature CD27+ cells grow faster than mature CD27- cells and vice versa. See ‘Materials and methods’ for further details. (c) Best fit of this model to the moments with the constraints Δkk0 and Csize-CD27+ < –0.2. Circles represent the fitted values of the moments to the model. Blue and red circles represent mean populations of CD27+, and CD27- NK cells, respectively. The variances of the CD27+ and CD27- NK cells are shown in cyan and pink circles, and the covariance between the CD27+ and CD27- NK cells is shown in purple. Horizontal lines around the circles show the standard deviations of bootstrapped moments. The solid line is the line y=x, which demonstrates where values would lie if we had a perfect fit. (d) Clones stochastically simulated with Gillespie’s algorithm (‘Materials and methods’) from the best fit parameters and the resulting clone sizes, CD27+ percentages, and correlation. (e) Observed clonal compositions of Ly6C+ and CD27+ cells from stochastically simulated clones. (f) Simulated clones from the three-stage model in (a) at the best-fit parameter values shown in (c) when the three stages are taken as immature CD27+Ly6C-, mature CD27-LyC- and the terminally mature CD27-Ly6C+, respectively. Each point is a simulated clone resulting from a single NK cell.
Figure 2—figure supplement 1
Three-stage models without death are able to capture negative correlations under the constraint Δk > 0.
(a) Schematic representation of a three-stage model where the first two stages are CD27+ and the last stage is CD27-. (b) A parameter scan for the model shown in (a). Color represents the value bI-bDPWe randomly initialized these simulations with an immature or double-positive cell according to the relative abundances of CD27+CD11b- and CD27+CD11b+ NK cells in the mass cytometry data. Thus, some clones are initialized with an immature cell, and others are initialized with an intermediate double-positive cell. (c) Schematic representation of a three-stage model where only the first stage is CD27+ and the last two stages are CD27-. (d) A parameter scan for the model shown in (c). Color represents the value bI-bINT.
Figure 2—figure supplement 2
A three-stage model with two CD27+ stages quantitatively captures clonal burst moments.
(a) Schematic representation of a three-stage model where the first two stages are CD27+ and the last stage is CD27-. Each subset undergoes both proliferation and cell death. (b) Best fit to the clonal burst data from the model in (a) given the constraints Csize-CD27+ < –0.2 and Δk > 0. We used an initial condition for moment solutions and Δk according to the relative abundances of CD27+CD11b- and CD27+CD11b+NK cells in the mass cytometry data. Thus, some clones are initialized with an immature cell, and others are initialized with an intermediate double-positive cell.
Figure 2—figure supplement 3
Mechanistic explanation for how a three-stage model with two CD27+ stages can create negative correlations under constraint Δk > 0.
(a) Some parameter combinations from the parameter scan in Figure 2—figure supplement 1b have bI≅r1, meaning that the change in the first subset is net zero on average. Thus, we investigated a model where this population is held constant to see if negative correlations are still possible. This model is represented here, with a constant inflow a into a CD27+ phenotype. This model can be interpreted as having a ‘stem cell-like progenitor’. (b) A parameter scan for the model shown in (a). Color represents the value of a. (c) Best fit to the clonal burst data from the model in (a) given the constraints Csize-CD27+ < –0.2 and Δk > 0.
Figure 3 with 2 supplements
Kinetics of homeostatic NK cells and endogenous NK cells responding to MCMV infection.
(a, b) Best fit to (a) Ly49H+ and (b) Ly49H- homeostatic NK cells. Y axis refers to percentage of the total Ly49H+ or Ly49H- population occupied by that NK subset. Circles and Xs represent the mean values of tamoxifen-induced td-Tomato positive and negative NK cells observed in the data, respectively. Solid and dashed error bars represent standard deviations around the means of tamoxifen-induced positive and negative NK abundances respectively (n = 7–9). Smooth curves represent model fits for tamoxifen-induced positive (solid) and negative (dashed) NK cells. Shaded regions refer to model 95% confidence bands, as determined by bootstrapping model fits. (c) Best fit to profiled endogenous NK cells responding to MCMV infection. Points represent mean cell abundances, and error bars represent standard deviations of bootstrapped means. Smooth curves represent model fit.
Figure 3—figure supplement 1
Gating strategies for mass cytometry data.
All gating strategies are representative from a single mouse’s NK cells. (a) Gating strategy for CD27 and CD11b from Ly49H+ NK cells. (b) Gating strategy for CD27 and Ly6C from Ly49H+ NK cells. (c) Gating strategy for Ly49H.
Figure 3—figure supplement 2
Supplemental fits to mass cytometry data.
(a) Simulated trajectories given the best-fit parameters for the clonal burst data in Table 1. (b) Best fit to mass cytometry with constraints that each parameter should be less than 2 and greater than –2 (or greater than 0 for differentiation rates r1 and r2).
Figure 4 with 1 supplement
Mature CD27- NK cells undergo rapid death during the expansion phase.
(a) Schematic of experimental setup. (b) Relative abundances of dead cells for CD27+ (blue) and CD27- (red) ex vivo NK cell subsets as measured by flow cytometry experiment. Data are representative of experiments with n=4–5 mice per timepoint. Graph shows means ± SEM; *p<0.033, ***p<0.001 represent statistically significant difference between CD27+ and CD27– NK cells as determined by two-way ANOVA. (c) Parameter scans for determining percentage of live cells for CD27+and CD27- subsets using the model shown in Figure 4—figure supplement 1c. Blue points represent a parameter scan where birth and death rates were varied by random sampling such that the net growth rates kI = bI - dI, kINT = bINT - dINT, and kM = bM - dM are equivalent to the kinetic estimates for endogenous NK cells shown in Table 3. Red points represent randomly sampled dead cell clearance rates with other rates set to those describing adoptive transfer kinetics in Table 1. Black dotted line is y=x, and points below this line show a higher percentage of live cells for immature CD27+ cells than for mature CD27- cells.
Figure 4—figure supplement 1
Models to measure percentage of live cells for each subset demonstrate that percentage of live cells is a good proxy for death rates.
(a) Schematic representation of a three-stage model with two CD27+ subsets, where all dead cells are cleared at a rate C. Dead cell abundances are measured, and %live cells can be calculated for CD27+ and CD27- cells by dividing the abundance of live cells by the sum of live and dead cells for each subset. (b) A parameter scan given the initial percentages of live cells observed in the experiment in Figure 3b shows that many parameter configurations generate a higher percentage of living CD27+ cells than that of CD27- cells. Of those parameter configurations, the difference in death rates is shown in this histogram. Because all CD27- death rates are higher, this indicates that a higher %live cells for CD27+ cells requires a higher death rate of CD27- cells. (c) Schematic representation of a three-stage model with two CD27- subsets, where all dead cells are cleared at a rate C. (d) A parameter scan given the initial percentages of live cells observed in the experiment in Figure 3b shows that many parameter configurations generate a higher percentage of living CD27+ cells than that of CD27- cells. Of those parameter configurations, the difference in death rates is shown in this histogram. Because all CD27- death rates are higher, this indicates that a higher %live cells for CD27+ cells requires a higher death rate of CD27- cells.
Figure 5
Proposed mechanism of NK clonal expansion in response to MCMV infection.
Three distinct subsets of NK cells are necessary to generate full NK response to infection. The first is CD27+, Ly6C-, and moderately proliferative, with a doubling time of ~0.3 days and negligible cell death. The second is CD27-, Ly6C-, and highly proliferative, with a doubling time of ~0.2 days. The last is CD27-, Ly6C-, and prone to rapid cell death with a half-life of ~0.45 days. This picture paints CD27-Ly6C- NK cells as having a highly proliferative and highly effective phenotype, similar to effector CD8+ T cells, and CD27+Ly6C- cells as long-lived, similar to memory precursor CD8+T cells.
Appendix 1—scheme 1
Schematic diagram showing the values of the numbers of immature (NI) and mature (NM) for which the Master equation was solved numerically.
The boundary layers ( and ) at the upper bounds of the numbers used to make the computation feasible are marked in grey. The values where NI and NM assume negative values are shown in blue.
Tables
Table 1
Best-fit parameters to three-stage model describing NK clonal bursts given in Figure 2a.
Confidence intervals are determined by bootstrapping.
| Parameter | Estimate (confidence interval) |
|---|---|
| bI | 1.064 day–1 (1.051–1.215) |
| bINT | 1.613 day–1 (1.561–1.996) |
| bM | 0.700 day–1 (0.001–0.727) |
| dI | 0.054 day–1 (0.052–0.232) |
| dINT | 0.285 day–1 (0.261–0.433) |
| dM | 1.370 day–1 (1.330–2.000) |
| r1 | 0.126 day–1 (0.088–0.172) |
| r2 | 0.164 day–1 (0.151–0.330) |
Table 2
Best-fit parameters to homeostatic Ly49H+ and Ly49H- NK cells.
Confidence intervals are determined by bootstrapping.
| Parameter | Ly49H+ estimate (confidence interval) | Ly49H- estimate (confidence interval) |
|---|---|---|
| λ | 1.366%day (1.117–1.594) | 2.412%/day (2.028–2.915) |
| r | 0.125 day–1 (0.095–0.159) | 0.147 day–1 (0.118–0.173) |
| kCD27+ | 0.097 day–1 (0.064–0.135) | 0.088 day–1 (0.053–0.119) |
| kCD27- | –0.120 day–1 (−0.152 to –0.091) | –0.107 day–1 (−0.126 to –0.086) |
Table 3
Best fit parameters to endogenous NK cells responding to MCMV infection with constrained r1.
Confidence intervals are determined by bootstrapping.
| Parameter | Estimate (confidence interval) |
|---|---|
| kI | 0.236 day–1 (0.153–0.350) |
| kINT | 0.537 day–1 (0.065–0.874) |
| kM | –0.484 day–1 (−1.194–0.371) |
| r1 | 0.126 day–1 (fixed) |
| r2 | 0.447 day–1 (0.000–0.874) |
Appendix 3—table 1
Best-fit parameters to the two-stage model with birth and death describing NK cell clonal bursts given in Figure 1e.
| Parameter | Estimate |
|---|---|
| bI | 0.662 day–1 |
| bM | 1.441 day–1 |
| dI | 0.129 day–1 |
| dM | 0.964 day–1 |
| r | 0.135 day–1 |
Appendix 3—table 2
Best-fit parameters to the asymmetric division model describing NK cell clonal bursts given in Figure 1—figure supplement 4.
| Parameter | Estimate (confidence interval) |
|---|---|
| bsym | 0.716 day–1 (0.704–0.811) |
| basym | 0.708 day–1 (0.663–0.725) |
| bM | 1.406 day–1 (1.362–1.442) |
| dI | 0.039 day–1 (0.026–0.097) |
| dM | 0.432 day–1 (0.401–0.481) |
| r | 0.017 day–1 (0.004–0.025) |
Appendix 3—table 3
Best-fit parameters to the three-stage model with two CD27+ stages describing NK cell clonal bursts given in Figure 2—figure supplement 1a.
| Parameter | Estimate (confidence interval) |
|---|---|
| bI | 0.899 day–1 (0.803–0.963) |
| bDP | 1.450 day–1 (1.247–1.801) |
| bM | 0.995 day–1 (0.992–2.000) |
| dI | 0.000 day–1 (0.000–0.007) |
| dDP | 0.004 day–1 (0.000–0.271) |
| dM | 0.076 day–1 (0.067–1.055) |
| r1 | 0.017 day–1 (0.000–0.018) |
| r2 | 0.591 day–1 (0.418–0.962) |
Appendix 3—table 4
Best-fit parameters to the two-stage model with constant inflow describing NK cell clonal bursts given in Figure 2—figure supplement 3a.
| Parameter | Estimate |
|---|---|
| bI | 1.464 day–1 |
| bM | 1.596 day–1 |
| dI | 0.332 day–1 |
| dM | 0.521 day–1 |
| r | 0.263 day–1 |
| a | 0.886 cell/day |
Appendix 3—table 5
Best-fit parameters to endogenous NK cells responding to MCMV infection.
| Parameter | Estimate (confidence interval) |
|---|---|
| kI | 0.108 day–1 (0.034–0.239) |
| kINT | 1.775 day–1 (0.018–2.000) |
| kM | –2.000 day–1 (−2.000–0.367) |
| r1 | 3.85×10–18 day–1 (0.000–0.106) |
| r2 | 1.482 day–1 (0.000–1.720) |
Appendix 3—table 6
Best-fit parameters to the three-stage model with two CD27- stages describing NK cell clonal bursts given in Figure 2a, but while fitting to comparative growth rates of CD27+ vs. CD27- cells rather than constraining to Δk > 0.
Because r2 ≃ 0 and Δk < 0, we report the constrained fit in the main text.
| Parameter | Estimate |
|---|---|
| bI | 1.375 day–1 |
| bINT | 2.000 day–1 |
| bM | 0.263 day–1 |
| dI | 0.350 day–1 |
| dINT | 0.866 day–1 |
| dM | 2.000 day–1 |
| r1 | 0.015 day–1 |
| r2 | 0.000 day–1 |
Appendix 3—table 7
Mass cytometry antibody panel.
| Metal | Channel | Target |
|---|---|---|
| In | 113 | CD45 |
| In | 115 | – |
| La | 139 | CD62L |
| Pr | 141 | Ly6C |
| Nd | 142 | CD132 |
| Nd | 143 | NKR-P1B |
| Nd | 144 | CD16 |
| Nd | 145 | – |
| Nd | 146 | CD96 |
| Sm | 147 | 2B4 |
| Nd | 148 | IFNAR |
| Sm | 149 | DNAM-1 |
| Nd | 150 | CD127 |
| Eu | 151 | Ly49D |
| Sm | 152 | CD49b |
| Eu | 153 | CD32 |
| Sm | 154 | CD69 |
| Gd | 155 | CD8 |
| Gd | 156 | CEACAM-1 |
| Gd | 157 | CD3 |
| Gd | 158 | NKG2ACE |
| Tb | 159 | Ly49H |
| Gd | 160 | NK1.1 |
| Dy | 161 | TIM-3 |
| Dy | 162 | – |
| Dy | 163 | TIGIT |
| Dy | 164 | CD2 |
| Ho | 165 | CD11b |
| Er | 166 | CD137 |
| Er | 167 | CD26 |
| Er | 168 | Ly49CI |
| Tm | 169 | NKp46 |
| Er | 170 | NKG2D |
| Yb | 171 | LAG-3 |
| Yb | 172 | KLRG1 |
| Yb | 173 | CD19 |
| Yb | 174 | Ly49G2 |
| Lu | 175 | IL18Ra |
| Yb | 176 | CD90 |
| Bi | 209 | CD27 |
Additional files
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Clonal stochasticity in early NK cell response to mouse cytomegalovirus is generated by mature subsets of varying proliferative ability
eLife 14:RP104951.
https://doi.org/10.7554/eLife.104951.4
