**DOI:**10.1128/JVI.00266-10

## ABSTRACT

Seasonal and pandemic influenza A virus (IAV) continues to be a public health threat. However, we lack a detailed and quantitative understanding of the immune response kinetics to IAV infection and which biological parameters most strongly influence infection outcomes. To address these issues, we use modeling approaches combined with experimental data to quantitatively investigate the innate and adaptive immune responses to primary IAV infection. Mathematical models were developed to describe the dynamic interactions between target (epithelial) cells, influenza virus, cytotoxic T lymphocytes (CTLs), and virus-specific IgG and IgM. IAV and immune kinetic parameters were estimated by fitting models to a large data set obtained from primary H3N2 IAV infection of 340 mice. Prior to a detectable virus-specific immune response (before day 5), the estimated half-life of infected epithelial cells is ∼1.2 days, and the half-life of free infectious IAV is ∼4 h. During the adaptive immune response (after day 5), the average half-life of infected epithelial cells is ∼0.5 days, and the average half-life of free infectious virus is ∼1.8 min. During the adaptive phase, model fitting confirms that CD8^{+} CTLs are crucial for limiting infected cells, while virus-specific IgM regulates free IAV levels. This may imply that CD4 T cells and class-switched IgG antibodies are more relevant for generating IAV-specific memory and preventing future infection via a more rapid secondary immune response. Also, simulation studies were performed to understand the relative contributions of biological parameters to IAV clearance. This study provides a basis to better understand and predict influenza virus immunity.

Current strategies for preventing or decreasing the severity of influenza infection focus on increasing virus-neutralizing antibody titers through vaccination, as experience indicates that this is the best way to prevent morbidity and mortality. Influenza A virus (IAV) undergoes mutations of the genes encoding the hemagglutinin (HA) and neuraminidase (NA) proteins that the neutralizing antibodies are directed against. When the variation is low (antigenic drift), prior vaccination often confers substantial heterologous immunity against a new seasonal IAV strain. In contrast, major genetic changes (antigenic shift) can result in pandemic IAV strains, since for novel strains, the humoral immune response is a primary response, and heterologous immunity is lacking. The emergence of such pandemic strains and the fact that young children are more vulnerable to influenza diseases highlight the need to better understand which viral and immune parameters determine the outcome of infection with viruses novel to the individual.

Conventional experimental methods to measure influenza virus immunity have been limited to animal models and studies of adult human peripheral blood leukocytes. The advantages of using animal models include the ability to intensively sample multiple tissues and to utilize genetic and other interventions, such as blocking or depleting antibodies, to dissect the contribution of individual arms of the immune system. However, it is easy to question the relevance of these experiments to humans because of the many important biological differences between human and murine immune systems (29). In both the animal and human systems, we are limited to measuring those parameters and variables for which assays are available, most of them being *ex vivo*. Parameters such as cell-to-cell spread of the virus *in vivo*, trafficking of immune cells to the lung, and the *in vivo* interactions in an intact immune system are much more difficult or impossible to measure with contemporary techniques, particularly in humans. Computational approaches have the potential to offset some of these limitations and provide additional insight into the kinetics of the IAV infection and the associated immune response.

Animal models of influenza virus infection in which different arms of the immune system have been suppressed suggest that some components of the adaptive immune system are required for complete viral clearance, often termed a sterilizing immune response. For example, abrogation of the CD4 T-cell response by cytotoxic antibody therapy or through knockout of major histocompatibility complex (MHC) class II slightly delays viral clearance but has little overall effect on the ability to control the infection (21, 54, 55). Elimination of the CD8 T-cell response typically results in delayed viral clearance (12, 20, 47), although animals with intact CD4 T-cell and B-cell compartments are able to control the infection in the absence of CD8 T cells. Presumably, this occurs through antibody-mediated mechanisms (54). Most animals depleted of both CD8 T cells and B cells are not able to clear the virus, which results in death (14, 32, 53). CD4^{+} T cells certainly contribute to the control of IAV infection, although cytotoxic CD4 T cells are not frequently observed unless cultured *in vitro* (8, 22, 45). Thus, it is generally accepted that CD8 T cells and/or antibodies are sufficient for timely and complete IAV clearance. Studies that strictly separate the relative roles of CD8 T cells and virus-specific antibodies are less satisfying. Animals depleted of both CD4 and CD8 T cells generally do not control the infection, despite substantial production of anti-IAV IgM antibodies (4, 23, 33, 34). However, adoptive transfer of IAV-specific IgM or IgG antibodies is protective (40, 51), suggesting that the timing and magnitude of the antibody response, i.e., the affinity, avidity, and antibody isotype, are important protective factors.

While murine gene knockout or lymphocyte depletion studies are highly informative, they also have a number of limitations. Most importantly, the near-complete ablation of one component of the adaptive immune system often causes profound and unpredictable effects on many other immune components. Although the reported experimental measurements are highly quantitative, they often focus only on a limited number of time points or measurements and do not capture the complexity of the altered, or intact, immune response. In contrast, high-frequency experimental sampling, coupled with mathematical modeling techniques and new statistical approaches, can give insights into the complex biology of IAV infection and test the assumptions inherent in the model. We have learned from other systems, particularly HIV (19, 35, 37, 38, 56), that quantitative analysis of the biology can reveal important factors that are not intuitively obvious. For example, our current estimates for the rates of HIV production and the life span of productively infected cells *in vivo* were obtained via mathematical modeling (35).

Mathematical models have long been used to investigate viral dynamics and immune responses to viral infections, including influenza A virus (3, 5, 7, 15, 16, 31, 36, 48). We recently described a complex differential equation model to simulate and predict the adaptive immune response to IAV infection (24). This model involves 15 equations and 48 parameters, and because of its complexity, many of the parameter values that could not be directly measured were unidentifiable. Thus, it is difficult to estimate all model parameters by fitting experimental data directly to this complex model, although the model can be used to perform simulation predictions (25). This issue can, however, be addressed by reducing the model into smaller submodels with smaller but identifiable sets of parameters, which can be estimated from experimental data. In this paper, we describe such an approach which focuses on IAV infection and the immune response solely within the lung.

In the present report, we have fitted a model of primary murine influenza virus infection to data. In naïve subjects, our data suggested that there is no adaptive immune response (e.g., IAV-specific CD8^{+} T cells or antibodies) detectable in the spleen, lymph nodes, or lung until approximately 5 days after infection; therefore, we have divided the analysis into the following two phases: the initial preadaptive (innate) phase and the later adaptive phase. We use direct experimental data from infection of mice with the H3N2 influenza virus A/X31 strain (2, 24) to obtain key kinetic parameters. The model fitting has revealed that the duration of the infection depends on a small set of immune components, and even large fluctuations in other arms of the immune system or IAV behavior have surprisingly little impact on the outcome of the infection.

## MATERIALS AND METHODS

Murine kinetic experiments.C57BL/6 mice (Jackson Laboratory, ME) from 6 to 16 weeks of age (*n* = 340) were anesthetized with 2,2,2-tribromoethanol and intranasally inoculated with 0.03 ml of 0.96 × 10^{5} 50% egg infective doses (EID_{50}) of H3N2 A/Hong Kong/X31 influenza A virus. Baseline mice (day 0) were inoculated immediately prior to organ harvesting. A total of 14 of 340 mice were excluded from the study because of sampling issues (death or poor sample quality). Time points used for viral titers were 0, 0.125 (3 h), 0.25 (6 h), 0.5 (12 h), 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 8, 9, 10, 12, and 14 days postinfection. Tissue for determination of CD8 T-cell counts and serum for antibody measurements were collected at 0, 1, 2, 3, 4, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 14, 16, 20, 28, 35, 62, and 84 days postinfection. At each time point, data were collected from 3 to 15 mice (9 to 12 mice for most time points for flow cytometry and enzyme-linked immunospot [ELISPOT] assays and 3 to 6 mice per time point for viral titer measurements).

Tissue collection.On the day of organ harvest, mice were euthanized using a lethal dose of 2,2,2-tribromoethanol. Lung, mediastinal lymph node, and spleen tissues (whole) were isolated from individual mice. Single-cell suspensions were prepared by Dounce homogenization or mechanical disruption with fine mesh and suspended in minimal essential medium (MEM) with 5% fetal bovine serum. Red blood cells were lysed using buffered ammonium chloride solution (Gey's solution), and cells were resuspended in complete MEM (cMEM) for analysis. Total lymphocyte numbers were obtained by manual counting using a hemacytometer and confirmed in many samples by flow cytometry measurements.

Virus titers and serum antibody.Lungs processed for viral titer measurements were separated from those used for cell samples. Lung homogenates were serially diluted and inoculated into embryonated hen eggs to determine the 50% egg infective dose (EID_{50})/ml, as previously described (2, 24). Blood collected by cardiac puncture was allowed to clot, and the serum was collected and frozen for determination of IAV-specific antibody levels by enzyme-linked immunosorbent assay (ELISA). Briefly, ELISA plates were coated with 200 μl of 1 μg/ml of sucrose gradient-purified, detergent-disrupted X31 virus or anti-murine Ig(H+L) capture antibody, blocked with 0.5% bovine serum albumin (BSA) in phosphate-buffered saline (PBS), and then incubated with 5 μl of serum per 200 μl well at 37°C for 30 min. Plates were then washed, and the appropriate well was incubated with either anti-IgM or anti-IgG horseradish peroxidase (HRP)-conjugated secondary antibodies at 37°C for 30 min, washed once more in PBS with 0.1% Tween, and then developed as previously described (18). Data were collected on a 96-well plate reader. Absolute antibody concentrations were calculated from a standard curve generated with known concentrations of murine IgG and IgM (Sigma, St. Louis, NJ).

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 naïve B6.SJL (CD45.1^{+}) mice were used as antigen-presenting cells and infected with influenza virus (multiplicity of infection [MOI] = 1) in 1 ml serum-free media for 60 min. Infected cells were then washed and resuspended in cMEM. A total of 1 × 10^{6} antigen-presenting cells were added to 1 × 10^{6} responders for a total volume of 100 μl. GolgiPlug (BD) was then diluted to 1 μl/ml in cMEM, and 100 μl was added to each well. Cells were incubated for 5 to 6 h at 37°C. Samples were then surface stained for CD3, CD4, and CD8 and then washed and resuspended in Cytofix/Cytoperm (BD) at 100 μl/well for 15 min. After one Perm/Wash (BD), anti-gamma interferon (IFN-γ), anti-tumor necrosis factor alpha (TNF-α), and anti-interleukin-2 (IL-2) antibodies were added to Perm/Wash, and cells were incubated for 30 min on ice in the dark. Samples were resuspended in PBS-BSA for fluorescence-activated cell sorting (FACS). Samples were analyzed on a BD LSRII cytometer, and the results were processed using FlowJo software (Tree Star, Inc.). Background values in duplicate wells lacking virus were subtracted. The number of antigen-specific CD8 cells per lung was calculated as follows: (CD8 cytokine subset cell count [flow]/lymphocyte count [flow]) × (total number of cells in pooled lung/number of pooled lungs). 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.

Mathematical models and statistical methods.To describe the dynamic interactions among the components (epithelial cells, virus, CD8 CTLs, and antibody) in the lung, the following differential equation model was developed based on our previous work (24) to quantify immune responses in primary infection,
(1) where the definitions of the variables and parameters are listed in Table 1, and a schematic illustration of model component interactions is given in Fig. 1. In the first equation of model 1, the infectible target (epithelial) cells *E _{p}* have a growth rate ρ

_{E}, and viral infection occurs with a rate β

_{α}. In the second equation, infected epithelial cells are generated due to viral infection with the term β

_{α}

*E*and eliminated by CD8

_{P}V^{+}T effectors

*T*with a rate

_{E}*k*or by cell death (virus induced or apoptotic) with a constant rate δ

_{E}_{E}

_{*}. In the third equation, free infectious influenza virions

*V*are produced by infected epithelial cells with a rate π

_{α}and cleared nonspecifically with a rate

*c*or neutralized by virus-specific IgG (

_{V}*A*) and IgM (

_{G}*A*) antibodies, with rates

_{M}*k*and

_{VG}*k*, respectively.

_{VM}We measured the antigen-specific CD8^{+} T-cell counts (CTLs) in the lung, IAV-specific serum IgG and IgM antibody concentrations, and viral titers [denoted by *T _{E}*(

*t*),

*A*(

_{G}*t*),

*A*(

_{M}*t*), and

*V*(

*t*), respectively] in mice at various time points after primary IAV challenge with the X31 virus. IgM antibodies to influenza virus are first detectable during the immune response to a primary IAV infection on day 6 (Fig. 2 D), with a peak antibody concentration on day 8 postinfection. Virus-specific IgG antibodies appear 12 to 24 h later (Fig. 2C). In this study, the smoothed curves of the CD8

^{+}T-cell counts in the lung and the IAV-specific serum IgG and IgM antibody concentrations were substituted into model 1 and taken as known covariates. The viral titer data were used to fit the model (see “Data processing for model fitting”) and estimate the remaining unknown parameters.

According to the data (Fig. 2), the adaptive immune response, as measured by virus-specific CD8^{+} CTL cells, IgG and IgM antibodies, is negligible within the first 5 days after infection. We therefore assume that the viral dynamic model for preadaptive immune responses can be simplified to the following:
(2) In this study, model 2 is primarily used to fit the viral titer data from days 0 to 5 (preadaptive phase), and model 1 is used to fit the viral titer data from days 5 to 14 (the adaptive phase) as well as the combined two-phase data (from days 0 to 14).

Before fitting the models to data, it is critical to verify whether all parameter values can be uniquely determined from the model and data. We performed structural identifiability analyses (58) for models 1 and 2. All the kinetic parameters in models 1 and 2 were found structurally identifiable except the parameter π_{α}. However, we can fortunately prove that arbitrarily fixing π_{α} will not change the estimates of all other parameters, except for *E _{p}*(0). This was previously shown to be true for a similar HIV model (49). Our model fitting results also confirm this identifiability analysis result. For details about identifiability analysis, see reference 30.

When fitting models 1 and 2 to the experimental data, we used the following measurement model, log_{10}*Y*(*t _{i}*) =

*V*(

*t*) + ε(

_{i}*t*), where

_{i}*i*equals 1, 2, …

*n*;

*V*(

*t*) denotes the true viral titer value at time point

_{i}*t*;

_{i}*Y*denotes the measured value of

_{i}*V*(

*t*); and ε(

_{i}*t*) denotes the independently and identically distributed (

_{i}*i.i.d.*) measurement errors. Note that log transformation is used in the measurement model to stabilize the computational algorithms. The standard nonlinear least-squares (NLS) approach is employed to obtain parameter estimates. To reliably locate the global minimum of the residual sum of squares (RSS), the hybrid optimization method was used (25, 44). Additionally, parameters are scaled to the same order to improve the reliability of the estimation results since they are orders of magnitude different. The 95% confidence intervals of parameter estimates are calculated using the bootstrap method (28). For more details, see reference 59.

We explored various model assumptions to obtain reliable parameter estimates. Since the model structure changes according to different assumptions, we report the Akaikes Information Criterion corrected (AICc) scores (1, 9) for each model for the purpose of model comparison and selection:
where *n* is the number of observations, and *k* is the number of unknown parameters. A lower AICc score indicates a better model.

Data processing for model fitting.The viral titer data used for model fitting are plotted in Fig. 2A. To incorporate the information in the antigen-specific CD8^{+} T-cell (CTL) data (Fig. 2B) and the antibody data (Fig. 2C and D), we used a nonparametric approach (R routine “smooth.spline” function) to smooth the CD8^{+} T-cell and the antibody data first (the degrees of freedom used in smoothing are 16, 11, and 12 for CD8^{+} T cells, IgG, and IgM, respectively; also, all smoothed curves are subject to the nonnegative constraint). We then substitute the smoothed curves into our mathematical models as time-varying covariates without measurement errors, which is a justified method, as discussed in our previous study (26).

The viral load data below detection, as explicitly labeled in the raw data file, are imputed as log_{10}(0.5) for the model fitting purpose. Furthermore, we assume that the initial number of infected epithelial cells is zero at the beginning of infection; that is, (0) equals 0 cells per lung. The initial viral load can be determined from the data directly. A linear model was fitted to the log-transformed data collected before 1.5 days, excluding data at time zero. The reason that the data point at time zero was excluded is that it takes a certain amount of time for viral particles to spread out and infect target cells, such that the observed drop of viral titer from hour 0 to hour 3 is artificial. The estimated initial condition of viral load was as follows: *V̂*(0) ≈ 10^{3.1682} ≈ 1473 EID_{50}/ml.

Model simulation software.Our parameter estimation and simulations based on the proposed differential equation models were performed using a software tool, DEDiscover, developed by our Center for Biodefense Immune Modeling, which is publicly available on its website (https://cbim.urmc.rochester.edu/software). The models described in this article are also available as built-in modules in the software.

## RESULTS

Primary data for model fitting.Data obtained from the experimental murine infection are shown in Fig. 2. Influenza A/X31 viral titers, measured from lung tissue, are shown in Fig. 2A and peak between days 2 and 3, with complete viral clearance by day 11. Virus-specific CD8^{+} T cells are not detectable before day 5, peak between days 9.5 and 11, and return to baseline approximately 20 days after infection. Absolute amounts of virus-specific serum IgG and IgM levels are shown in Fig. 2C to D, with a peak IgM level at 8 to 10 days postinfection, and return to baseline levels after 20 days. In contrast, IgG levels peaked at day 25, with a slow decay over the next 60 days.

Estimation results for preadaptive immune response kinetics.In a primary influenza infection, virus-specific T- and B-cell responses are difficult to detect until 5 days after viral inoculation (39, 46). Assuming that the contributions of the adaptive antibody and cellular responses of the immune system are negligible during this period, we used model 2 to fit the viral titer data during the first 5 days of infection. All kinetic parameters could be estimated with statistical significance, and confidence intervals are provided with each estimate in Table 2. The net growth rate (the proliferation rate minus the death rate) of uninfected epithelial cells is estimated as 6.2 × 10^{−8} day^{−1}, which is statistically significant but essentially negligible. The estimated infection rate of epithelial cells is 2.4 × 10^{−6} ml·(EID_{50})^{−1}·day^{−1}. The estimated death rate of infected epithelial cells is 0.60 day^{−1}, corresponding to a half-life of 1.2 days. The estimated clearance rate of virus particles is 4.2 day^{−1}, corresponding to a half-life of 3.9 h. In addition, the predicted number of infectible target (epithelial) cells drops rapidly to nearly zero around day 2 (Fig. 3 A). This does not mean that the lung is devoid of epithelial cells. Rather, it implies that the number of permissive target cells for the virus is very low due to mechanisms not explicitly defined in the model, such as the antiviral effects of type I IFN.

Given that the net growth rate of uninfected epithelial cells (ρ_{E}) during the preadaptive immune response period (days 0 to 5) appeared to be negligible, we next asked whether ρ_{E} is a significant contributor to the outcome of the infection. To address this issue, we compared the model containing parameter ρ_{E} to the model without ρ_{E} (or equivalently, ρ_{E} ≠ 0 and ρ_{E} = 0, respectively) using the model selection criteria, AICc scores. We found that the AICc scores are distinct for the models with and without ρ_{E} (−15 for the model with ρ_{E} and −18 for the model without ρ_{E}), indicating that the model without ρ_{E} was a better model fit to the data. In addition, the estimates of all other parameters remained nearly the same when ρ_{E} was removed, strongly suggesting that the net growth rate of uninfected epithelial cells is negligible during the first 5 days of IAV infection.

The production rate of virus per infected epithelial cell (π_{α}) is difficult to measure directly, and there are no reliable estimates *in vivo*. Model estimation results suggest that fixing the value of π_{α} between 1 and 10^{4} EID_{50}/(ml·day) per infected cell does not change estimates of other parameters, such as the infection rate. This observation is consistent with the conclusion of the identifiability analysis. However, careful examination of the time trajectories of the number of uninfected and infected epithelial cells (Fig. 3) revealed that π_{α} has a strong negative correlation with the initial number of uninfected epithelial cells, *E _{p}*(

*t*= 0), and therefore with the number of uninfected (

*E*) and IAV-infected () epithelial cells during infection. That is, the effects of π

_{p}_{α}are compensated by the corresponding changes in

*E*and , such that the magnitude of

_{p}*E*and changes for different values of π

_{p}_{α}. Because of this balance, the estimates of all other parameters are essentially unaffected by π

_{α}. Finally, it should be noted that π

_{α}does have a significant effect on the viral titer [

*V̂*(

*t*), e.g., the peak] if the compensation by

*E*is removed; that is, if

_{p}*E*(

_{p}*t*= 0) is fixed instead of being an unknown parameter. Although the total number of epithelial cells in the mouse lung can be estimated, because not all such cells are targets of infection, then

*E*(

_{p}*t*= 0) would still not be known. In summary, the virus production rate (π

_{α}) is very important when predicting the preadaptive and adaptive IAV immune response given its significant effect on viral titer

*V̂*(

*t*), although it has no effect on other model parameter estimates when

*E*(0) is not fixed.

_{p}Estimation of adaptive immune response kinetics.From day 5 onward during primary influenza infection, IAV-specific CD8^{+} T cells become detectable in the lung, and IAV-specific antibody-secreting B cells can be detected in the lymph nodes, spleen, and peripheral blood, signaling the effective start of the adaptive immune response. We next analyzed the model during the adaptive phase by fixing the initial conditions at day 5 of uninfected *E _{p}*(

*t*= 5) and infected (

*t*= 5) epithelial cells and virus

*V*(

*t*= 5) at the values estimated from the preadaptive phase model (Table 2). We also varied the parameter π

_{α}from 1 to 10

^{4}EID

_{50}/(ml·day) per infected cell, and again, it did not significantly affect the estimates of all other parameters.

After model fitting, we found all parameter estimates are statistically significantly different from zero. The parameter estimates and the fitted curves for different π_{α} values are shown in Table 2 and Fig. 3B, respectively. The net growth rate of uninfected epithelial cells is 0.34 day^{−1} (doubling time of 2.0 days), the infection rate is 0.78 ml·(EID_{50})^{−1}·day^{−1}, the CTL killing rate is 6.4 × 10^{−5} cell^{−1} day^{−1}, the death rate of infected epithelial cells due to factors other than CTL killing is 0.15 day^{−1}, the estimated virus neutralization rate due to IgG is 3.1 × 10^{−5} ml/(pg·day), the estimated virus neutralization rate due to IgM is 4.4 ml/(pg·day), and the estimated virus clearance rate (due to all other factors) is 4.0 × 10^{−4} day^{−1}. These estimates along with the observation that the IgM and IgG serum concentrations differ only by 4-fold or less suggest that in primary IAV infection, the contribution of IgM to clearance of infectious virus is much greater than that of IgG. Similarly, the model fitting results (Fig. 4 C) suggest that elimination of infected cells is dominated by the CTL response compared to other mechanisms.

Estimation results obtained from the combined two-phase fitting.For completeness and comparison, model 1 was also fitted to the viral titer data from days 0 to 14 (the combined two phases). In this case, we assume that all kinetic parameters are constant for days 0 to 14 during the two phases. The parameter estimates and the curves fitted for different π_{α} values for the combined two-phase data are shown in Table 2 and Fig. 3C, respectively. The net growth rate of uninfected epithelial cells is 8.6 × 10^{−4} day^{−1}, which is again essentially negligible. The infection rate is 5.1 × 10^{−6} ml·(EID_{50})^{−1}·day^{−1}, the CTL killing rate is 1.4 × 10^{−5} cell^{−1} day^{−1}, the death rate of infected epithelial cells due to factors other than CTL killing is 1.2 day^{−1}, the estimated virus neutralization rate due to IgG is 4.5 × 10^{−8} ml/(pg·day), the estimated virus neutralization rate due to IgM is 0.078 ml/(pg·day), and the estimated virus clearance rate (due to all other factors) is 5.5 × 10^{−6} day^{−1}. The RSS (49.4) of fitting both phases is approximately 1% lower than the summation (49.9) of the RSS of fitting two separate phases. Since the conclusions from combined two-phase fitting are essentially similar to those from fitting the two phases separately, the following discussion continues to separately consider the two phases.

Relative contributions of immune response kinetics.Certain parameters, such as those characterizing the infection of target cells and virus replication, are shared by both the early (preadaptive) and the late (adaptive) phases of the infection. It is therefore of interest to compare the early and adaptive immune responses. For this purpose, we plotted the estimates of parameters common in models 1 and 2, and we also compared the contributions of different components in model 1, as shown in Fig. 4. Both the net epithelial cell growth rate (ρ_{E}) and the infection rate (β_{α}) during the adaptive immune response are more than 5 orders of magnitude higher than those rates during the early response. However, these increases may be artificial since the number of target cells (*E _{p}*) rapidly drops to nearly zero by day 2 (see Fig. 3), such that the two terms ρ

_{E}

*E*and β

_{p}_{α}

*E*in the first equation of model 1 become negligible during the adaptive immune response phase, and therefore, the estimates of ρ

_{P}V_{E}and β

_{α}cannot be reliably determined from the data.

To identify the main components contributing to the clearance of infected epithelial cells and of free virus particles during the adaptive phase of the immune response, we developed 10 submodels and compared them to the full model (model 1) using the model selection criteria, AICc scores. Due to the extensive computational time required for model fitting, we fixed π_{α} as 100 EID_{50}/(ml·cell·day), which is close to published virus production rates (200 virions per cell per day) (10, 17, 50). The results are reported in Table 3, from which we found that the model without the contributions of CTLs and antibody has the highest AICc score (2.69), suggesting that the effects of CTLs and antibodies are not dispensable. The model with only the effects of CTLs and IgM has the lowest AICc value (−24.1). This implies that CTLs and early IgM antibodies are the two main contributors to the clearance of infected epithelial cells and free virus particles, respectively, in a primary IAV infection, with the effect of IgG production being negligible even 1 day later.

For comparison, we also evaluated 9 submodels for the combined two-phase fitting. From Table 3, we found that the model containing only the effects of virus infection, infected cell death due to factors other than CTL killing, and IgM has the lowest AICc value (−40.1). This result confirms the significant contribution of virus-specific IgM antibody to virus clearance. Furthermore, this result also suggests, although the CTL effect dominates the killing of infected epithelial cells during the adaptive phase, the overall contribution of infected cell death due to factors other than CTL killing is also important (Table 3 and Fig. 4).

Effect of the net growth of uninfected epithelial cells.As suggested by the model fitting and model selection results, the effect of the new growth of uninfected epithelial cells was found to be negligible during days 0 to 14 after infection. One reason for this is that the replacement rate of epithelial cells is low relative to the rate of virus spread and production. For example, the average half-life of the ciliated epithelial cells in the lungs of uninfected mice was reported to be as long as 17 months (41), corresponding to a death rate of about 0.001 day^{−1}. Therefore, in steady state, the proliferation rate of the ciliated epithelial cells must be also low to balance the slow death of epithelial cells. In reality, the replacement of epithelial cells is accelerated when the epithelial cell layer is damaged, though it may still be slow in comparison to virus dynamics. In addition, our model does not explicitly consider other target cells such as macrophages and dendritic cells (24), which could contribute to virus production. All cells permissive to productive IAV are subsumed under the term “epithelial cells.” More rapid proliferation or infiltration of new target cells into the lung could account for disparities between this model and real-world situations, resulting in observations of higher titers or prolonged replication.

Experimental observations showed that the doubling time of epithelial cells could be as short as 8 h (6, 57), corresponding to a proliferation rate of 2.1 day^{−1}, which is 8 orders of magnitude higher than the estimated proliferation rate (6.2 × 10^{−8} day^{−1}) during the first 5 days. It is therefore of great interest to understand how the peak viral titer will change if the doubling time of uninfected epithelial cells is shortened, e.g., from 17 months to 8 h. For this purpose, we increased the proliferation rate in model 2 from 10^{−8} day^{−1} to 10^{0} day^{−1}, holding all other parameters constant. The peak viral titer doubled (plot not shown). Although doubling of the viral load is significant, it is low compared to the change of 8 orders of magnitude simulated for the proliferation rate. This suggests that epithelial proliferation is a minor contribution to viral load, even at supraphysiological rates.

For the adaptive immune response, if we assume that CTLs are the only elements that actively kill infected epithelial cells and that virus particle elimination is mediated by IgM, how are viral load kinetics and other parameters affected in the absence of other factors? To address this question, we fit model 1 to data from days 5 to 14 with the assumption that epithelial cell proliferation, infection, and death by viral cytopathic effects are all at zero (e.g., ρ_{E} *= 0*, β_{α} = 0, and δ_{E}_{*} = 0). Compared with the full model, the killing rate of infected cells by CTLs (*k _{E}*) increased from 6.4 × 10

^{−5}to 7.1 × 10

^{−5}cell

^{−1}·day

^{−1}, and the rate of IgM-mediated elimination of infectious virus particles (

*k*) increased from 4.4 to 5.2 ml/(pg·day). These changes in parameter estimates compensated such that the fitted curve of viral titer changed little. In other words, contributions of the infection and proliferation of uninfected epithelial cells and the death of infected epithelial cells by viral cytopathic effects do not significantly influence the viral kinetics during the adaptive immune response phase.

_{VM}We next asked if an accelerated IgM response would improve viral clearance in a primary IAV infection, assuming IgM secretion begins earlier. The data suggest that IgM reaches its peak around day 8. We shifted its values 3 days earlier to address this question. Fitting model 1 to the shifted IgM data and viral titer data from days 5 to 14 (Fig. 5) suggests that the initial drop at around day 5 for the shifted IgM data is much more dramatic than that for the original IgM data due to the sudden increase of the total amount of IgM antibody at day 5 after shifting. We also found that the relative importance of CTLs and IgM does not change.

Mortality associated with different viral strains.Uncontrolled virus replication can lead to high titers, and even a modest viral load cannot be maintained for more than a few days before death will occur (61). It is of interest to understand how the kinetics of viral load change with increasing infection rate β_{α}. For this purpose, we focused on the early, preadaptive phase of the infection. If we fix all other parameters in model 2 and increase only β_{α} from 1.0 × 10^{−6} to 1.0 ml·EID_{50}^{−1}·day^{−1}, the magnitude of viral titer peak essentially does not change. Instead, the peak time moves closer to day 0 (Fig. 6 A). This simulation result suggests that the infection rate has little effect on the magnitude of the peak viral load. As indicated in Fig. 6B, if we fix all other parameters in model 2 [including the initial number of target cells *E _{p}*(

*t*= 0)] and increase only the viral replication rate π

_{α}from 10 to 10

^{3}EID

_{50}·ml

^{−1}·cell

^{−1}·day

^{−1}, the peak viral titer increases from 10

^{6}to 10

^{8}EID

_{50}. Further increase of π

_{α}from 10

^{3}to 10

^{4}EID

_{50}·ml

^{−1}·cell

^{−1}·day

^{−1}does not augment the peak any more due to target cell limitation. In addition, if π

_{α}is as low as 1, the viral titer monotonically decreases to zero. This simulation study suggests that a viral strain, such as the recent H1N1 A/California strain, which has a higher replication rate but not a higher infection rate, might pose a more serious threat to public health.

Effect of the number of available target cells.We next used model 2 to investigate the influence of the available target cells on viral titer. It would be extremely difficult to control the initial number of uninfected cells under *in vivo* experimental conditions. However, the mathematical models can be used to explore interesting questions that cannot be directly answered by experiments. As indicated in Fig. 7, if we increase the initial number of target cells *E _{p}*(

*t*= 0) from 10

^{5}to 10

^{8}cells while all other parameters are held constant, the peak viral titer increases from 10

^{6}to 10

^{9}EID

_{50}/ml. This suggests that the availability of target cells can seriously limit or boost the increase of viral titer. However, if the initial number of target cells

*E*(

_{p}*t*= 0) is as low as 10

^{4}or 10

^{3}cells, viral titer monotonically decreases, and no peak can be observed.

## DISCUSSION

In this work, we have developed a differential equation model of influenza A virus infection and fitted this model to a comprehensive data set from primary IAV infection in the mouse. The data set is unique in that it represents the results from simultaneous high frequency (every 12 to 24 h) sampling of lymphocytes, antibody, and influenza virus in a large number of animals (*n* = 340), with multiple replicates at each data point. The large number of frequent measurements gives much greater power in fitting the model to data and led to several key findings, which we summarize here.

Our model divides the immune response to IAV infection into the following two phases: the preadaptive (innate) and adaptive immune response phases. During the innate response phase (days 0 to 5), adaptive immune response elements (IAV-specific Ig and lung-infiltrating CTLs) are negligible. Although clonal expansion of IAV-specific CTLs and B cells may be starting in the lymph nodes and spleen prior to day 5, measurements of the effectors in the lung are below detection limits during this period. It is during the innate phase that the viral infection is spread among susceptible respiratory epithelial target cells, and each infected cell becomes a viral production source, thus contributing to peak viral load. The total number of epithelial cells includes those that are susceptible to viral infection (targets), those that are infected, those that are immune to infection, and those that are regenerating cells (either target or immune cells). *A priori*, key parameters that determine the viral load during the innate phase are virus infectivity and replication rates, host parameters such as the number of infected, target, and regenerating epithelial cells, and the rate of epithelial cell death. The adaptive phase of the immune response limits viral load by active viral clearance and by the killing of IAV-infected cells. We use our model to investigate the sensitivity of the mathematical models to changes in specific parameters and to estimate kinetic parameters that are difficult to measure directly.

A key set of model predictions, regarding the factors that affect the peak viral load, are as follows: (i) virus infectivity changes the timing but not the magnitude of the peak viral load, (ii) peak viral load is target cell limited, and (iii) peak viral load is very sensitive to the viral replication rate but not to infectivity. The magnitude of peak viral load is relatively insensitive to infectivity, which may seem counterintuitive. However, within physiologic parameter ranges, the viral replication rate is much more significant, and moderate to high rates of replication during this period of low viral clearance are quantitatively the primary determinant of peak viral load. These results imply that viruses with higher production would be more pathogenic than those with high infectivity. Although not investigated here, infectivity may be more relevant to spread of the virus from person to person rather than a primary determinant of virus levels within an individual.

Consistent with these findings regarding the significance of virus production versus spread in the host is the prediction that peak viral load is limited by the available target cells. Thus, once established in an epithelial cell, the pace of virus production and the amount of virus produced have a greater influence than cell-to-cell spread, particularly as the infection progresses and the number of available targets decrease. During the adaptive phase, the actions of antibody and CTLs outpace virus replication and spread, causing titers to decline. Thus, we predict that therapies that impede or slow virus replication would have the most significant impact during the innate phase of infection. Similarly, robust therapies that limit target cell availability may also have a strong effect. In fact, the type I interferon system is an example of how biology has evolved a mechanism to achieve exactly this effect.

Using several techniques, our modeling predicts that epithelial cell proliferation has a negligible impact on determining the kinetics of IAV infection. Parameter ρ_{E}, the net growth rate of uninfected epithelial cells per day, is included in the absolute number of target cells available for the virus to infect.
The reliability of the estimate of ρ_{E} is greatly affected by the magnitude of the entire term, reflecting the net growth of uninfected epithelial cells ρ_{E}*E _{p}*. Once the infection is established, estimating the term for epithelial proliferation becomes unreliable, primarily because the number of infectible target cells,

*E*, is nearly zero. Even if the proliferation rate of the epithelial cells increases, the value of the entire term (ρ

_{p}_{E}

*E*) is still low once most of the target cells are infected (e.g.,

_{p}*E*≈ 0) and has no effect on the adaptive model predictions. Similar analysis holds for estimation of the number of viral infections per unit time (β

_{p}_{α}). Essentially, a very small number of target cells (

*E*) and a large amount of virus (

_{p}*V*) cannot be used to reliably estimate the epithelial proliferation (ρ

_{E}) or infection rates (β

_{α}) because most of the target cells have become infected ().

In our models, infected epithelial cell death occurs through the following two modes: killing by CTLs and death by other factors that are a composite of immune (complement, ADCC, and NK cells, etc.) and natural (necrosis and apoptosis) factors. Considering the overall infection (days 0 to 14), death by factors other than CTL killing contributed the most to loss of infected cells. However, during the adaptive phase (days 5 to 14), model selection analysis suggests that the contribution of CTL killing is more significant than the contribution of other factors. Specifically, based on the adaptive phase-fitting results, CTLs are the main contributors to the clearance of infected epithelial cells, with a killing rate (*k _{E}*) of 6.4 × 10

^{−5}per CTL per day. The antigen-specific CD8 T-cell number peaks at >5 × 10

^{4}cells, yielding an effective peak killing rate (

*k*·

_{E}*T*) of >3.2 day

_{E}^{−1}. Comparing the innate to the adaptive phase of the immune response, we found that the half-life of infected epithelial cells decreases from an average of 1.16 to 0.52 days, primarily due to the effects of CTLs. The death rate (δ

_{E}

_{*}) of infected epithelial cells due to factors other than CTLs is 0.16 per day during the adaptive phase, while δ

_{E*}equals 0.6 per day during the first 5 days. Whether this is a real effect or an artifact caused by the much stronger effect of CTL killing is uncertain. Nonetheless, since most of the infected epithelial cells are killed by CTLs more efficiently, the other causes of infected cell death do not alter the infected epithelial cell half-life significantly. To further evaluate the impact of the CTL killing on viral clearance, we also performed computer simulations to compare the outputs of model 2 with and without CTL killing (equivalently,

*k*equaling 6.4 × 10

_{E}^{−5}per cell per day and

*k*equaling 0 per cell per day). We find that the virus can be completely cleared about 10 days after infection when CTL killing is included in the model. When we omit the CTL killing effect, however, our model predicts that the virus cannot be cleared (instead, the virus titer reaches a nadir of about 3.5 log

_{E}_{10}EID

_{50}on day 9 and then starts to rebound). This observation suggests that the CTL effect has a big impact on influenza virus clearance.

Our estimates of killing rates for murine CTLs are consistent with estimates of the killing rates for CD8 T cells for lymphocytic choriomeningitis virus (LCMV)-infected cells (4.5 × 10^{−5} per CTL per day), even with differences in the target cells and virus (13). Others have also estimated CD8 CTL killing rates in murine LCMV infection by fitting a closed-form model to experimental data, with the estimated half-life of infected cells being about 10 min (6.9 × 10^{−3} days) (42, 43, 60). There are several factors that could explain the difference in target cell half-life between the two systems. LCMV infection was measured in the spleen, where infected target cells and effector CTLs are in close contact with each other, with greater freedom to migrate. In the IAV-infected epithelium, target cells are arranged in an essentially two-dimensional surface, and CTL migration to infected target cells in the airways may be more limited. Therefore, the effective half-life of the target cells is longer in the influenza model because of limited target accessibility, even though the per cell killing rates are similar.

Our results suggest that IgM is crucial for IAV clearance in a primary infection, whereas the contribution of IgG is less important, despite the eventual presence of higher and more sustained levels of IAV-specific IgG. During the adaptive phase of infection (>5 days), parameter estimates show IgM-mediated clearance of free virus (*k _{VM}*) at a rate of 4.4 ml/(pg·day). In contrast, IgG antibody had a relatively low contribution to viral particle clearance, with a rate (

*k*) of 1.0 × 10

_{VG}^{−5}ml/(pg·day). During the adaptive phase, our estimate of the clearance rate of virus particles due to factors other than antibody neutralization (

*c*) is <0.01 per day, which is much lower than that during the first 5 days of infection. This translates into a shortened half-life of infectious viral particles from 4 h (days 0 to 5) to 1.8 min (days >5), on average, predominantly due to IAV-specific IgM. Even if we fix the value of

_{V}*c*at the innate level when estimating the rates of virus inactivation by IgG and IgM (

_{V}*k*and

_{VG}*k*, respectively), the estimates nearly do not change, indicating that the overall model is insensitive to changes in

_{VM}*c*. This is a particularly interesting finding in that early IgM responses would likely be without significant affinity maturation and thus have lower affinity for the viral antigens than later IgG. In considering these observations, keep in mind that virus titers drop 3 or more orders of magnitude from days 5 to 8. Experimental observations show that class-switched IgG antibody, which is generally dependent on virus-specific CD4

_{V}^{+}T helper cells, does not become detectable in significant quantities in our data until day 7, a time when virus titers are relatively low. Therefore, the majority of free virus elimination occurs when IgM, plus the uptake of virus by target cells, phagocytosis, and mucociliary processes, numerically dominate. Nevertheless, we realize that the measurement of total IgG and IgM antibodies may overestimate the amount (and type) of antibodies ultimately involved in clearance of viral particles. Future experiments focused on measuring virus-antibody complexes may yield better insights.

It is worth discussing some of the unexpected ways in which these predictions and fitting results work together. If the viral load is dependent primarily on the rate of virus production, but not spread, then it becomes intuitive that the elimination of the virus “factories” by CTLs has the largest impact on reducing virus titers. This prediction is borne out by both the biological observations (2, 11, 12, 52) and the model fitting of the adaptive phase to the experimental data. In the same way, the effect of antibody-mediated inhibition of virus spread by IgG would have relatively little influence since it effectively diminishes infectivity. Infectivity of the virus was not found to significantly predict influenza viral titers.

Our results are also consistent with the other previously published modeling studies of the IAV immune response in human subjects (3). Despite many biological and physical differences between mice and human subjects, many of the parameter estimates are remarkably close. In our study, the estimated overall infection rate of epithelial cells is 5.1 × 10^{−6} ml·EID_{50}^{−1}·day^{−1}, while the overall death rate of infected epithelial cells is 1.2 per day. Baccam et al. (3) investigated the kinetics of influenza A virus infection in human subjects using a target cell-limited model fitted to data on nasal IAV shedding obtained from six human subjects. The infection rates were estimated to be from 2.7 × 10^{−6} to 1.1 × 10^{−3} (TCID_{50}/ml)^{−1} day^{−1}, which are in the same range of our estimates. The death rates of infected target cells were estimated to be between 2.1 and 13.5 day^{−1}, higher than our estimate of 1.2 day^{−1}. The difference in the epithelial cell death rates could be due to the relatively high virus production value assigned in our study, which was chosen based on direct measurements from *in vitro* systems using cultured airway epithelial cells (6, 27, 57). Based on these estimates, we fixed virus production at 100 EID_{50}·ml^{−1}·cell^{−1}·day^{−1} and estimated innate virus clearance (*c _{V}*) as 4.2 day

^{−1}. Baccam et al. (3) estimated the viral production rates per cell to be 3.0 × 10

^{−3}to 5.9 × 10

^{−1}(TCID

_{50}/ml) day

^{−1}and the viral clearance rate to be 2.1 to 13.5 day

^{−1}. In addition, Handel et al. (16) investigated the development and spread of neuraminidase inhibitor resistance in humans infected by influenza virus. The estimated death rate of infected cells shown in the work of Handel et al. (16) is 0.5 to 1.3 day

^{−1}, which is very close to our estimates listed in Table 2 (i.e., 0.60 day

^{−1}for the early phase, 0.15 day

^{−1}for the adaptive phase, and 1.2 day

^{−1}for the combined two phases). The estimated virus death rate shown in the work of Handel et al. (16) is 0.081 to 3.0 day

^{−1}, which is lower than our estimated innate virus clearance rate of 4.2 day

^{−1}.

In conclusion, we have constructed and validated a robust and quantitative model of primary influenza infection. Our data suggest that an effective adaptive immune response appears to begin around day 5 after infection, when viral titers are high. Our model indicates that the viral replication rate is more important than viral infectivity, with respect to viral titers, for a primary infection. In addition, virus-specific IgM and CD8 cytotoxic T lymphocytes are the most important factors in the adaptive immune response to clear virus in a primary infection. These findings suggest that future work might focus on increasing the early IgM response and vaccines which boost virus-specific CD8 T-cell production and memory.

## ACKNOWLEDGMENTS

We acknowledge the excellent and dedicated technical support of Tina Pellegrin for the murine experiments. We also thank Amber Smith for very useful comments and suggestions.

This work was funded by the National Institute of Allergy and Infectious Diseases contracts N01-AI-50020 and R01 AI069351 (to M.S.Z.). Portions of this work were done under the auspices of the U.S. Department of Energy under contract DE-AC52-679 06NA25396 (to A.S.P.).

## FOOTNOTES

- Received 4 February 2010.
- Accepted 8 April 2010.
- *Corresponding author. Mailing address for David Topham (immunology questions): 601 Elmwood Avenue, Box 609, Rochester, NY 14642. Phone: (585) 273-1403. Fax: (585) 273-2452. E-mail: david_topham{at}urmc.rochester.edu. Mailing address for Hulin Wu (modeling questions): 601 Elmwood Avenue, Box 630, Rochester, NY 14642. Phone: (585) 241-0705. Fax: (585) 256-2541. E-mail: hwu{at}bst.rochester.edu
↵▿ Published ahead of print on 21 April 2010.

## REFERENCES

- American Society for Microbiology