## Abstract

The biological parameters that determine the distribution of virus-specific CD8^{+} T cells during influenza infection are not all directly measurable by experimental techniques but can be inferred through mathematical modeling. Mechanistic and semimechanistic ordinary differential equations were developed to describe the expansion, trafficking, and disappearance of activated virus-specific CD8^{+} T cells in lymph nodes, spleens, and lungs of mice during primary influenza A infection. An intensive sampling of virus-specific CD8^{+} T cells from these three compartments was used to inform the models. Rigorous statistical fitting of the models to the experimental data allowed estimation of important biological parameters. Although the draining lymph node is the first tissue in which Ag-specific CD8^{+} T cells are detected, it was found that the spleen contributes the greatest number of effector CD8^{+} T cells to the lung, with rates of expansion and migration that exceeded those of the draining lymph node. In addition, models that were based on the number and kinetics of professional APCs fit the data better than those based on viral load, suggesting that the immune response is limited by Ag presentation rather than the amount of virus. Modeling also suggests that loss of effector T cells from the lung is significant and time dependent, increasing toward the end of the acute response. Together, these efforts provide a better understanding of the primary CD8^{+} T cell response to influenza infection, changing the view that the spleen plays a minor role in the primary immune response.

Current strategies for preventing or decreasing the severity of influenza infection focus on increasing virus-neutralizing Ab titers through vaccination, as experience indicates that this is the best way to prevent morbidity and mortality. It is generally accepted that T cells provide a substantial degree of protection from disease, and cytotoxic CD8^{+} T cells control the infection by eliminating infected cells when neutralizing Abs are absent (1–3). Because influenza viruses are unique in their ability to undergo antigenic shift, resulting in a wholesale evasion of neutralizing Abs, T cells that cross-react with conserved antigenic regions in the viruses are perhaps the most important immune component against emerging novel strains of the virus. Understanding the factors that regulate the T cell response is key to the development of new immunization strategies designed to promote cross-reactive immunity.

With the exception of certain highly pathogenic strains of influenza (4–6), the virus infection is usually restricted to the airways and lung tissues. This phenomenon is believed to be due to the limited expression of host enzymes required to cleave the viral hemagglutinin protein into its active conformation (7, 8). Because of this, it takes time for Ag (and APCs) to reach the draining lymph node (9, 10) and also for virus-specific effector CD8^{+} T cells to migrate through the bloodstream back to the site of infection. The dynamic nature of this process makes it very difficult to study directly. Conventional approaches rely on tissue sampling at defined time points to assemble a set of “snapshots” that can be linked together to understand the system as a whole. Unfortunately, this understanding still depends on a number of assumptions. For example, a given tissue sampling tells us where the cells are at that point in time but does not reveal where they came from. In addition, the rates of cell proliferation and trafficking into and out of a given tissue can, at best, be inferred. In the end, the number of cells at any single tissue site will depend on the combined, simultaneous processes of proliferation, trafficking, retention, and cell loss due to death, phagocytosis, or sloughing.

Mathematical modeling and computational approaches to the problem offer one of the best means to estimate the key biological parameters that cannot be directly and unambiguously measured. Mathematical models have long been used to investigate viral dynamics and immune responses to viral infections, including influenza A virus (IAV) (11–18). In particular, we recently developed complex differential equation models to quantify key kinetic parameters and predict early and adaptive immune responses against IAV infection (19, 20). However, limited by available experimental information, our previous studies mainly focused on the lung. Better understanding of immune cell kinetics among multiple compartments is thus a natural and necessary next step. Toward this goal, in this article we propose new mathematical models of multicompartment cell trafficking and fit them to experimental data from mice infected with the H3N2 influenza A/X31 strain. Based on model fitting results, we compare and verify several mechanistic hypotheses that cannot be directly inferred without mathematical models.

## Materials and Methods

### Murine kinetic experiments

Female C57BL/6 mice (The Jackson Laboratory, Bar Harbor, ME) from 6 to 13 wk of age (*n* = 78) were anesthetized with 2,2,2-tribromoethanol and intranasally inoculated with 0.03 ml of 1 × 10^{5} EID_{50} H3N2 A/Hong Kong/X31 IAV. Tissues for determination of CD8^{+} T cell counts were collected at days 0, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, and 28 postinfection. At each time point, data were collected from six mice (for most of the time points) and from 12 mice on day 12 for flow cytometry. Data from days 5 to 14 were used in modeling fitting to understand the cell trafficking during the adaptive phase. All experiments involving animals were reviewed and approved by the University Committee for Animal Resources (institutional animal care and use committee).

### Tissue collection

On the day of organ harvest, mice were euthanized using a lethal dose of 2,2,2-tribromoethanol. Whole lung, mediastinal lymph nodes, and spleen tissues were isolated from individual mice. Single-cell suspensions were prepared by Dounce homogenization or mechanical disruption with a fine mesh and suspended in MEM with 5% FBS. RBCs were lysed using buffered ammonium chloride solution (Gey’s solution), and cells were resuspended in complete minimum essential medium (cMEM) for analysis. Total lymphocyte numbers were obtained by manual counting using a hemocytometer and confirmed in many samples by flow cytometry measurements.

### Determination of virus-specific CD8^{+} T cell immune responses

Cell phenotypes were determined by flow cytometry after surface and intracellular staining. To determine the number of virus-specific CD8^{+} T cells, spleen cells from naive B6.SJL (CD45.1^{+}) mice were used as APCs and infected with influenza (multiplicity of infection = 1) in 1 ml serum-free media for 60 min. Infected cells were then washed and resuspended in cMEM. APCs (1 × 10^{6}) were added to 1 × 10^{6} responders for a total volume of 100 μl. GolgiPlug (BD Pharmingen) was then diluted to 1 μl/ml in cMEM and 100 μl added to each well. Cells were incubated for 5∼6 h at 37°C. Samples were then surface stained for CD3, CD4, and CD8 and then washed and resuspended in Cytofix/Cytoperm (BD Pharmingen) 100 μl/well for 15 min. After one Perm/Wash (BD Pharmingen), anti–IFN-γ, anti–TNF-α, and anti–IL-2 Abs were added in Perm/Wash and cells incubated for 30 min on ice in the dark. Samples were resuspended in PBS/BSA for FACS. Samples were analyzed on a BD Biosciences LSRII cytometer, and results were processed using FlowJo software (Tree Star). The number of Ag-specific CD8^{+} cells per lung was calculated as follows:

For this analysis, the total number of CD3^{+} CD8^{+} cells expressing one or more of the three cytokines TNF-α, IFN-γ, and IL-2 (TNF^{−}IFN^{+}IL-2^{+}, TNF^{−}IFN^{+}IL-2^{−}, TNF^{+}IFN^{−}IL-2^{+}, TNF^{+}IFN^{−}IL-2^{−}, TNF-IFN^{−}IL-2^{+}, TNF^{+}IFN^{+}IL-2^{+}, and TNF^{+}IFN^{+}IL-2^{−}) was used as the number of virus-specific cells. See Fig. 1 for the observed CD8^{+} T cell counts in lymph node, spleen, and lung.

### Detection of Ag presentation by lymph node dendritic cells

APCs were isolated essentially as described (10, 21). C57BL/6 mice for influenza experiments were obtained from The Walter and Eliza Hall Institute of Medical Research. Mice were maintained under specific pathogen-free conditions. Experiments with mice began when they were between 6 and 10 wk of age. Mice were anesthetized with methoxyflurane and then intranasally infected with a nonlethal challenge of H3N2 A/Hong Kong/X31 IAV diluted in 25 μl sterile PBS. The regional (mediastinal) draining lymph nodes were digested for 20 min at room temperature with collagenase/DNase and then treated for 5 min with EDTA to disrupt T cell–dendritic cell (DC) complexes. Previously, it has been shown that depletion of cells expressing CD11c removed Ag-specific stimulatory capacity indicating that Ag presentation was limited to CD11c^{+} DCs (9, 22). The LacZ-inducible hybridomas specific for D^{b} nucleoprotein (NP)_{366–374} (BWZ-IFA.NP4) (9) and D^{b} acid polymerase (PA)_{224–233} (BWZ-IFA.PA1) were maintained in DMEM containing 10% FCS, 50 μM 2-mercaptoethanol, 2 mM l-glutamine, 100 U/ml penicillin, and 100 μg/ml streptomycin (hybridoma medium). These clones were used to analyze Ag presentation by lymph node cells, as previously described (23).

### Mathematical models and assumptions

First, we consider a mechanistic model involving lymph node, spleen, and lung to describe CD8^{+} T cell response to primary IAV infection. For the mediastinal lymph node (MLN), we assume that there exist proliferation, migration to spleen, migration to lung, migration to other compartments, and apoptosis (death) or loss due to other mechanisms, and all these dynamic rates are proportional to the number of activated, effector CD8^{+} T cells, , in the MLN. Furthermore, we assume that the influx of CD8^{+} T cells from other compartments to the MLN can be ignored compared with the other kinetic components mentioned above. For spleen, the number of CD8^{+} T cells, , is assumed to be mainly affected by proliferation in place, influx from MLN, migration to lung, and loss due to death or other mechanisms. For CD8^{+} T cells in lung, , direct influx from both lymph node and spleen is considered, and the loss is again due to death or other reasons. We further make the initial assumptions that the activation of CD8^{+} T cells that leads to expansion of the virus-specific population in both MLN and spleen is stimulated mainly by Ag-bearing DCs, with negligible expansion due to proliferation in lung. Based on these tenets, we formulated the following mechanistic model:(1)where denotes the number of mature Ag-bearing DCs in MLN; is the number of mature DCs in spleen; *τ* is the time delay of the effects of DCs on CD8^{+} T cell activation; and are the proliferation rates of CD8^{+} T cells stimulated by DCs in MLN and spleen, respectively; , , and are the loss rates in MLN, spleen, and lung, respectively; is the migration rate from MLN to spleen, the migration rate from MLN to lung, and the migration rate from spleen to lung. Finally, data for are not available in this study or the literature to the best knowledge of the authors. is thus used to replacein all calculations, with the assumptions that the APCs in the MLN and spleen share similar temporal patterns but can differ in total number proportional to the size of the two organs.

We make explicit and strong assumptions in the model above. However, we recognize that numerous factors may affect CD8^{+} T cell trafficking after influenza virus infection, and their mechanisms are not fully understood; therefore, the mechanistic model above may be oversimplified and thus not able accurately to interpret the experimental observations. We hence explore three more hypotheses: first, that the effects of other undefined factors and mechanisms that affect CD8^{+} T cell growth and disappearance are not negligible; second, that the rise and fall of CD8^{+} T cell numbers in the infection is determined by clonal expansion followed by increased programmed cell death when the pathogen is cleared, versus a constant rate of loss over the acute phase; third, that an influx of T cells from other tissue compartments such as the liver (24) could explain late increases in splenic CD8^{+} T cell numbers not due to proliferation (Fig. 1).

To study the first hypothesis, we consider the following model:(2)where two unknown time-varying parameters and are introduced to replace and and should be interpreted as the net growth rates of CD8^{+} T cells in MLN and spleen, respectively, not distinguishing proliferation from disappearance/death. To study the second hypothesis, we introduce a time-varying loss rate of CD8^{+} T cell in lung, removing the two previous unknown time-varying parameters introduced in model 2, so the third model is

To study the third hypothesis, we consider a time-varying flux of CD8^{+} T cells from other additional compartments to spleen at a non-parametric rate , and we also keep the time-varying loss rate of CD8^{+} T cell in lung. This alternative model can be written as

In models 2 to 4, because the time-varying parameters do not depend on explicit assumptions about how CD8^{+} T cells are stimulated or disappear in spleen, we call such models “semimechanistic.”

All these alternative hypotheses are explored and tested from our experimental data using the statistical methods introduced below. We summarize our mathematical notations and parameter definitions in Table I. Also, the schematic illustration of the CD8^{+} T cell migration pathways is depicted in Fig. 2.

### Statistical methods

We performed structural and practical identifiability analysis for the four ordinary differential equation (ODE) models presented earlier before we applied statistical methods to estimate the parameters in these ODE models (25, 26). From the structural identifiability analysis, all the parameters were found theoretically identifiable. Further, we performed practical identifiability analysis to confirm whether all the parameters are estimable from experimental data with noise using Monte Carlo simulations (25, 27, 28). Our results demonstrate that almost all the parameters are practically identifiable except for , , and .

We estimated all the kinetic parameters including both constant and time-varying parameters as well as the delay time and the initial conditions for the ODE models from the time course data of activated CD8^{+} T cell numbers in MLN, spleen, and lung using the nonlinear least squares method (29–32). We used a hybrid global optimization algorithm DESQP in Ref. 30 to minimize the residual sum of squares. The log-transformation and standardization of the data were used to stabilize the estimation algorithm. The time-varying parameters in models 2 to 4 were approximated by a sum of third-order B-splines of df 3 (33), for example, and , so that the time-varying parameter ODE models can be transformed into standard ODE models with constant parameters (32). The number of knots and basis functions can be determined using the model selection criterion such as the Akaike Information Criterion (AIC) (34, 35). We used Fisher Information Matrix to calculate 95% confidence interval of the parameters.

We explore and evaluate various alternative models and model assumptions using the AIC score defined as (34, 35):,where RSS is the residual sum of squares, *n* is the number of observations, and *k* is the number of unknown parameters. The AIC is a score that can be used to evaluate different models by compromising the fitting residuals and model complexity (measured by the number of parameters in the model). A smaller AIC score indicates a better model. In this study, the number of observations (*n*) is 180 in total (from day 5 to day 14), and the number of parameters (*k*) is 9 in model 1, 12 in model 2, 11 in model 3, and 14 in model 4.

Some of our parameter estimations and simulations for ODE models were performed using DEDiscover, a publicly available tool developed by the University of Rochester Center for Biodefense Immune Modeling (https://cbim.urmc.rochester.edu/software). DEDiscover is a user-friendly modeling environment for both mathematical modelers and immunologists that provides a workflow-oriented interface and a comprehensive set of computing algorithms. In particular, to help modelers and immunologists to compare and select models built upon alternative assumptions, DEDiscover also reports the AIC score after model fitting.

## Results

### Model fitting to determine key cellular parameters

Our objective using the modeling approach is to estimate biological parameters that are difficult or impossible to measure directly. These include such processes as cell trafficking from one anatomical compartment to another and growth of the virus-specific CD8^{+} T cell population. Initially, we fitted model 1 to the experimental data from mice for virus-specific CD8^{+} T cells in the MLN, spleen, and lung. The virus-specific T cells were identified and measured using an intracellular cytokine assay with live-virus restimulation and staining for IL-2, IFN-γ, and TNF-α as previously described in our original global influenza model (19). CD8^{+} T cells that stained positive for any combination of one to three cytokines were counted as positive. Prior to day 5, the difference between virus-restimulated and control cultures was close to zero and inconsistent, making it impossible accurately to fit the model to the data. In contrast, virus-specific CD8^{+} T cells become easily detectable from day 5 onward by a variety of assays, including the intracellular cytokine approach (19, 24, 36–39). Therefore, data from days 5 to 14 were used for the model fitting. The experimental data and the fitted curves are plotted in Figs. 1 and 3, with the estimated parameters reported in Table II.

### Growth/loss rates from mechanistic model (model 1)

In the mechanistic model (model 1; see *Materials and Methods*), we make the assumption that the expansion of CD8^{+} T cells in both MLN and spleen is mainly a consequence of Ag presentation analogous to mature Ag-bearing DCs that can activate naive T cells. We also allow a time delay for stimulation of T cell expansion, the length of which is unknown but can be estimated from data. The concentrations of DCs in MLN and spleen are based on the two different parameters and in our model. In more explicit terms, our assumptions are that: 1) the stimulation strength for activation and growth of the virus-specific CD8^{+} T cell population is proportional to the number of DCs; and 2) the temporal patterns of stimulation in MLN and spleen are similar, but the numbers can be very different based on the size of the organs. Based on these assumptions, we fit this model, and the results are summarized in Table II. Both the residual sum of squares (15.77) and AIC score (−420.29) of this model are the largest compared with the other three models, which suggests that this mechanistic model does not fit the data as well, possibly due to an oversimplification of the model assumptions or omission of key biological parameters included in the alternative models.

The estimated DC/APC-dependent effector CD8^{+} T cell population expansion rates for MLN and spleen, and , respectively have peak values of 1.3 d^{−1} and 3.5 d^{−1} at around day 7, corresponding with a doubling time of 12 h and 4.8 h. This suggests that the growth of the virus-specific CD8^{+} T cell population in spleen is ∼2-fold greater than that in the MLN. Similar conclusions were drawn from the results derived from models 3 and 4 described later (see Table II and Fig. 4).

The time delay for stimulation of T cell expansion by DCs is calculated as 3.08 d, matching well with direct measurements in the literature (40–43). The estimates of the migration rates of effector CD8^{+} T cells from MLN to spleen and from spleen to lung, *γ _{ms}* and

*γ*, are 0.16 d

_{sl}^{−1}and 0.5 d

^{−1}, respectively, which suggests that migration from spleen to lung is substantial. The constant loss rate of effector CD8

^{+}T cells from the lung,

*δ*, could be reliably (

_{l}*p*< 0.001) estimated at ∼4 d

^{−1}, corresponding with a half-life of 4.2 h. However, it was not possible reliably to determine the values for the constant loss rates (

*δ*and

_{m}*δ*) of CD8

_{s}^{+}effectors from the MLN and spleen or the migration rate of activated CD8

^{+}T cells from MLN to lung (

*γ*). This is due to the fact that the estimated parameter values are low and close to zero and because of this may be (statistically) obscured by noise in the data set.

_{ml}### Estimation of the growth rate of effector CD8^{+} T cells

Because it is difficult directly to measure or estimate explicit values for cell loss other than migration in the MLN and spleen, in model 2 we considered an alternative model in which we introduce combined terms for both cell proliferation and loss as the time-varying pure growth rates, and . The estimates of these two parameters are plotted in Fig. 4. Model 2 makes no assumptions on the patterns and magnitudes of the two time-varying parameters and . That is, we did not make an explicit assumption of what factors regulate the expansion of the CD8 effector population. We believe this approach has its advantages especially when some mechanisms are not completely clear. The patterns of and are similar to one another but, as with the population growth rates estimated in the first model, the value of is more than 2-fold greater than , which suggests a noticeably faster net generation of effector CD8^{+} T cells in the spleen than in MLN. For MLN, decreases from 2.1 d^{−1} at day 5 to 0.9 d^{−1} at around day 10; for spleen, decreases from 6.0 d^{−1} at day 5 to 3.6 d^{−1} at around day 12. The corresponding maximum doubling times of the virus-specific CD8^{+} T cell populations in MLN and spleen are 8 h and 2.8 h, respectively, which are more than 30% more rapid than the estimated doubling times in model 1 (at 12 h and 4.8 h in MLN and spleen, respectively). Therefore, it is likely that factors important to the expansion of the effector CD8^{+} T cell population in MLN and spleen are missing from model 1. The smaller AIC score for model 2 than that for model 1 (−428.27 versus −420.29) reinforced this interpretation indicating that model 2 is a better fit to the data. Also note that the short doubling time in the spleen suggests additional factors at work, such as the contribution of other compartments to the growth rate. This is considered explicitly later.

However, the estimated migration rates from MLN to spleen and from spleen to lung as well as the disappearance rate in lung (1.1 d^{−1}, 4.2 d^{−1}, and 32 d^{−1}, respectively) are all an order of magnitude higher in model 2, raising questions about the accuracy of these estimates. The high values could be due to the strong correlations between the proliferation parameters and the migration and disappearance rates.

### Estimates of effector CD8^{+} T cell migration and disappearance

Next, the parameters for cell trafficking and constant (models 1 and 2) or time-varying death or loss (model 3) were estimated. These constant loss rates represent the disappearance of cells from the system due to death or physical removal (sloughing or phagocytosis), which cannot be easily resolved from one another and are treated as a single term. As mentioned earlier, in comparison with other parameters, the estimates for the constant pure loss rates for MLN and spleen, *δ _{m}* and

*δ*, were essentially zero and could not be reliably determined. Similarly, the migration rate of virus-specific CD8

_{s}^{+}T cells from MLN to lung,

*γ*, was also estimated as close to zero. A low AIC score for models in which these parameters are excluded, indicating a good fit, reinforces the conclusion that these rates have a minimal impact on the fidelity of the model and the behavior of the system.

_{ml}Notably, when we change the constant disappearance rate *δ _{l}* in lung in model 1 to time-varying as in model 3, the fit of the model to the data improved further, and AIC score decreased from −420.29 to −473.40, which suggests the pattern of cell loss in lung cannot be simplified as a constant over time. Instead, the estimated time trajectory of is plotted in Fig. 4

*B*, which shows that the disappearance rate of Ag-specific CD8

^{+}T cells in lung monotonically increases from 0 at day 5 to 9.7 d

^{−1}at day 14. The shortest half-life of CD8

^{+}T cells therefore occurs at day 14 with a value ∼2 h, demonstrating an increase in the rate of effector cell loss at the end of the immune response.

Collectively, these results suggest increased effector cell clearance toward the end of the acute response, particularly in the lung. In addition, there is a strong indication that the spleen plays a much more substantial role in the overall CD8^{+} T cell response to the infection than may have previously been appreciated.

### Modeling the contribution of other tissues

In contrast to the growth of the virus-specific CD8^{+} T cell population, T cell trafficking into peripheral tissues is an Ag-independent process (44, 45). Activated T cells that enter the circulation have the ability to traffic to many tissues other than the site of infection in the lung (46, 47). Because T cell migration through peripheral tissues is dynamic and continuous, and assuming some T cells pass through the tissues and re-enter the circulation via the lymphatics, we investigated the contribution of influx of activated CD8^{+} T cells from compartments other than the three we considered. To accomplish this, model 4 was developed, which included an additional compartment corresponding with all other tissue sites combined (see model 4 in *Materials and Methods*). Also, as suggested by model 3, we retain the time-varying disappearing parameter in lung. The fitting results suggest that a time-varying influx *η _{s}*(

*t*) from other tissues to spleen is significant, contributing a maximum of 8.2 × 10

^{4}cells per day. Also, the estimates of all other constant parameters are very close to the estimates from models 1 and 3, suggesting that unlike model 2, the introduction of the two time-varying parameters in model 4 do not affect the estimates of other constant parameters, and the results are therefore reliable. The smallest AIC score (−477.59) also suggests that model 4 is superior to all three other models. Such results suggest that it is very likely that there exists a significant, time-varying influx of CD8

^{+}T cells from other tissues to spleen; also, the rate of disappearance of effector CD8

^{+}T cells in lung should be monotonically increasing (see Fig. 4). For convenience, we plotted the product of and in Fig. 5. We see that the number of cells lost per day increases from 0 at day 5 to 2.7 × 10

^{5}cells per day at day 9 and then decreases down to 8 × 10

^{4}cells per day at day 14. These results suggest that sites other than the MLN and spleen are making substantial contributions to the effector CD8

^{+}T cell pool and moderate the obvious contributions of the spleen as a sole organ. Instead, these results are consistent with the view that, in addition to in situ expansion, the spleen may collect effector T cells in the circulation that have come from other tissues. These findings do not diminish the conclusion that cells from the spleen substantially contribute to the response in the lung.

## Discussion

We have developed ODEs and time-delay ODEs to model the population growth, trafficking, and loss of virus-specific CD8^{+} T cells in MLN, spleen, and lung in mice with primary influenza A infection. These models allowed us to challenge various assumptions and test hypotheses about how the system works. In addition, the model fitting process generated estimates of several key biological parameters that are technically difficult or impossible to measure empirically. Several assumptions that were held prior to the start of this work proved difficult to justify in the models. In particular, the role of the spleen in the primary CD8^{+} T cell response to influenza infection may be more significant than previously thought.

There is a widely held view that during a primary influenza infection, the development of the virus-specific CD8^{+} T cell response occurs in the following sequential fashion: DCs bearing viral Ags migrate from the lung to the draining MLN where they activate naive virus-specific T cells. The T cells proliferate in the MLN and then leave the lymph node to enter the blood, appearing next in the lung after extravasation (24, 48–51). This whole process takes 4–5 d before virus-specific CD8^{+} T cells can be detected at any one site. Very careful measurements made with sensitive techniques such as MHC class I tetramer staining (24, 36) during the very early phase of the acute response (days 0–6) support this general view (24, 51–53). However, these approaches only tell us where the T cells first become detectable, not where they came from. Nor do they consider the relative sizes of the T cell populations in each site, or the population growth and trafficking rates that would be required for the MLN to be the major source of effectors. Nevertheless, when looking at such data, the earliest tetramer^{+} CD8^{+} T cells can be detected in the MLN at day 4–5, then appear simultaneously in the lung, spleen, and organs such as the liver at 5–6 d (24, 52, 53). After this point, the numbers of cells at each site rises very rapidly, presumably driven by a combination of proliferation and migration, with little initial loss. Unfortunately, these major contributing biological parameters and sites become hard experimentally to deconvolute through empirical observations. The various rates of activation, proliferation, migration, and death are difficult to measure directly and are often confounded by the difficulty to distinguish experimentally whether, for example, a T cell that has upregulated Ki67, diluted its CFSE label, or incorporated BrdU did so in the compartment that was sampled or in another site the T cell just came from (54, 55). However, by comparing alternative models through a process of fitting, we can resolve some of the factors.

We developed four models and compared them using mainly the AIC score. It is possible for two models to have close scores (e.g., model 3 versus model 4); however, as suggested by Burnham and Anderson (35), the difference instead of the absolute value of AIC score matters. Even a small decrease in AIC scores will still suggest superiority in model structures. In this study, the model that produced both the best fit to the experimental observations and the best AIC score is model 4. In this model, the expansion of the CD8^{+} effector T cell population is dependent on the number and kinetics of the professional, Ag-bearing APCs rather than the viral load. The effector CD8^{+} T cell populations expand in both the draining MLN and the spleen, but the spleen is estimated to make a much larger contribution to the pool of cells that migrate to the lung. Because cell loss in these sites is estimated to be very low, the size of the T cell population in the spleen, in turn, is likely dependent upon a combination of factors that include in situ proliferation and an influx of T cells from other tissue compartments. These other compartments would include cells returning to the circulation from peripheral nonlymphoid tissues not infected with the virus, but may also include lymphoid organs such as the cervical lymph nodes and bone marrow that have the capacity to contribute to the virus-specific T cells response (53, 56–58) but were not explicitly considered in our modeling or the experimental data set. Alternatively, the growth of the CD8^{+} T cell numbers in the spleen could also include those cells that were activated to proliferate outside of the spleen but continued a program of cell division (59). Other factors such as the local cytokine environment to promote T cell expansion were also not explicitly modeled, though such effects should still be proportional to the DC kinetics as suggested though our model fitting.

In addition to a time-varying term for the contribution of other sources of T cells to the spleen, the best fit to the data was provided by the models with a time-varying loss rate for the effector T cells, with values increasing toward the end of the acute response and highest in the lung. Models with constant rates of cell loss did not reliably reflect the experimental data. Regardless, effector CD8^{+} T cell loss rates from the MLN and spleen were estimated to be very low and close to zero compared with other parameters. These results suggest that T cell survival factors that include Ag or cytokines must wane toward the end of the response. It is also possible that the features unique to influenza infection of the lung, such as the epithelial site of infection and effector cell trafficking into the mucosa, drive some of these effects.

A notable advantage of time-varying parameters is that uncertainties in mechanisms can be well accounted for. Although this approach usually requires frequent time-course data, it could be the only right choice for complex dynamic problems. For example, although the introduction of one time-varying parameter will result in three more constant parameters in our models, the overall performance measured by AIC scores of time-varying parameter models still outperformed model 1, which has only constant parameters. In more than one of the models, we found that the growth rate of the virus-specific, effector CD8^{+} T cell population in spleen is twice as large as that in MLN, and both the experimental data and fitted-model simulations indicate that the MLN precedes the spleen by about 1 d. Furthermore, the numbers of cells and rate of accumulation in the spleen and the lung cannot be solely explained by proliferation in and/or migration from the MLN, suggesting that the MLN is not the sole source of antiviral effectors. Direct migration from the MLN to the lung and death rates in the MLN and spleen were all estimated to be very small relative to other contributors. Combined, these estimates support that the system operates differently than we originally assumed.

The finding of a minimal contribution of effector CD8^{+} T cell migration from MLN to the growth of this population in the lung deserves some qualification. These results do not imply that CD8^{+} T cells do not traffic from the MLN to the lung, nor do they change the view that the adaptive response is initiated in the draining lymph node. Instead, these estimates should be viewed in consideration of the entire systemic response, the relative sizes of the virus-specific T cell populations at each site over time, and in comparison with the rates of cell migration and accumulation in sites other than the lung. For example, the numbers of virus-specific effector CD8^{+} T cells at day 5, , , and in model 3, are estimated as 4.4 × 10^{3}, 3.5 × 10^{4}, and 1.3 × 10^{3} cells per MLN, spleen, and lung, respectively; the spleen being an order of magnitude above the other sites. Shortly after virus-specific CD8^{+} T cells become detectable in the system, the spleen becomes the major reservoir of these cells, and the distribution cannot be mathematically explained by growth and trafficking of the virus-specific effector population from the MLN alone. In contrast, the migration rates of activated CD8^{+} T cells from MLN to spleen and from spleen to lung in models 1, 3, and 4 are all larger than from MLN to lung and close to each other at around 0.2 d^{−1} and 0.5 d^{−1}, respectively.

This means that although common sense suggests that some of the cells must come directly from the MLN, in fact the contribution of this source is very small and primarily observable very early in the response when cell numbers are low. Later in the response, the only way to compensate would be to impose extraordinarily high cell growth and migration rates for the MLN.

The picture that emerges fits with observations made by Masopust, Reinhardt, and Marshall in 2001 showing that activated effector T cells traffic all over the body and not just to the sites of infection (46, 47, 53). Experiments in which the MLN is carefully digested with enzymes to liberate T cells complexed to APCs (60) do not change the values enough to alter this conclusion. This view is less controversial when considering that the spleen can function as a site of primary immunity under conditions that disrupt trafficking of naive T cells into the lymph nodes (56–58, 61). The original experiments, performed more than a decade ago, demonstrated a more than 10-fold decrease in virus-specific CTLs (and Th precursors) in the lung in mice that had been splenectomized (56). The conclusions about the role of the spleen stand in contrast to early reports of experiments with splenectomized mice that concluded the spleen made little contribution to the infection (62). A more likely interpretation is that splenectomy does not affect control of mild influenza infection because there are more than enough CTLs generated to control the virus. This is consistent with the observation that the absence of a spleen reduces the total number of CD8^{+} T cells generated.

One of the valuable features of modeling and simulation is the ability explicitly to explore alternative hypotheses and model structures. For example, we compared two models in which the T cell response was driven by either the amount of virus in the system (results not shown) or the number of professional, Ag-bearing APCs (DCs, or *D* as a variable in this article). The DC-driven model provided a substantially, and statistically better, fit to the actual experimental data. This suggests that the kinetics and magnitude of the CD8^{+} T cell response is limited more by the number of APCs available than the amount of Ag.

Several other key parameters were estimated from the models. These include the growth and disappearance of T cells in different compartments. In model 4, the growth rates [ in MLN and in spleen] increased rapidly between days 5 and 7 and then dropped with the decrease in DCs to reach a minimum around day 13 (Fig. 3). The maximum population doubling time in the MLN was ∼20 h, and the spleen followed a similar pattern but was estimated to have a maximum doubling time of 6.9 h, or ∼3.5 doublings per day. This is significantly slower than other cell-based estimates that range from as high as 5.3 to 12 cell divisions per day (40), but lower than others that estimate 17–19 h per division (63). These discrepancies may be accounted for by differences in population growth rates versus cell division rates.

There is also recent evidence that proliferation of CD8^{+} T cells occurs in the lung itself (54, 55), but this was not initially considered in our models. However, if we explicitly explore the hypothesis that the T cells proliferate in the lung by adding the term , then the death and proliferation rates in lung become so strongly correlated such that we cannot distinguish them from each other using this model and data set. Even when the proliferation rate in the lung is fixed to some value, only the estimate of the death rate in lung was affected. That is, the proliferation rate in lung cannot be determined separately based on current models and data. Therefore, we did not explicitly include the term for proliferation in the lung in our model. Instead, the parameter is, in effect, the net of the proliferation and disappearance rates for CD8^{+} cells in lung. The possibility of cell proliferation in the lung does not negate the conclusion that the spleen and other sites are a major contributor to the overall CD8^{+} T cell response. Migration of T cells that have proliferated in the lung back to the MLN or spleen is likely to be very small or negligible. Furthermore, although evidence of CD8^{+} T cell proliferation in the lung exists, it makes a substantial but not complete contribution to the total number of CD8^{+} T cells in the lung and seems to contribute the most at early time points in the infection (54).

Notably, the lung was the only site in which cell loss, presumably due to death and other factors, was estimated to be time varying with a maximum value of 9.7 per day, or a half-life of 1.7 h. The number of cells lost per day reached its maximum around day 9 (2.7 × 10^{5} cells per day, see Fig. 5). This is not out of line with the notion that continued encounter with Ag-bearing target cells and migration out into the infected airways (from which there is no return) account for a significant disappearance of the effector cells. If the cells in the lung were also proliferating as has been postulated, then the actual loss rate would have to be higher still. Keep in mind, however, that we treated cell loss as a time-varying parameter instead of a constant one. We believe our modeling methods and techniques can provide a more powerful platform for immunologists to understand viral and immune system dynamics.

Collectively, these efforts change the way the immune response to influenza infection is viewed. The spleen can no longer be considered as a location used merely to contain the “spillover” of excess virus-specific CD8^{+} T cells in the blood that were generated in the MLN and not taken up by the lung or other tissues. Instead, the spleen is a major contributor of effector CD8^{+} T cells against the virus infection. This view also implies that the spleen contains Ag-bearing APCs that can drive T cell proliferation and is an environment where T cell expansion can occur. As such, it may also be an important source of the memory T cells that remain at the end of the acute immune response.

## Disclosures

The authors have no financial conflicts of interest.

## Acknowledgments

We acknowledge the excellent and dedicated technical support of Tina Pellegrin and Joe Hollenbaugh during the murine experiments.

## Footnotes

This work was supported by National Institute of Allergy and Infectious Diseases Contract HHSN272201000055C and Grant R01 AI069351 (to M.S.Z.). Portions of this work were performed under the auspices of the U.S. Department of Energy under Contract DE-AC52-679 06NA25396 (to A.S.P.).

Abbreviations used in this article:

- AIC
- Akaike Information Criterion
- cMEM
- complete minimum essential medium
- DC
- dendritic cell
- IAV
- influenza A virus
- MLN
- mediastinal lymph node
- NP
- nucleoprotein
- ODE
- ordinary differential equation
- PA
- acid polymerase
- RSS
- residual sum of squares.

- Received May 16, 2011.
- Accepted August 23, 2011.

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