JOURNAL OF VIROLOGY, July 2010, p. 6687–6698 0022-538X/10/$12.00 doi:10.1128/JVI.00266-10 Copyright © 2010, American Society for Microbiology. All Rights Reserved.
Vol. 84, No. 13
Quantifying the Early Immune Response and Adaptive Immune Response Kinetics in Mice Infected with Influenza A Virus䌤 Hongyu Miao,1 Joseph A. Hollenbaugh,2 Martin S. Zand,3 Jeanne Holden-Wiltse,1 Tim R. Mosmann,2 Alan S. Perelson,4 Hulin Wu,1* and David J. Topham2* Department of Biostatistics and Computational Biology, University of Rochester, Rochester, New York 146421; David H. Smith Center for Vaccine Biology & Immunology, Department of Microbiology and Immunology, University of Rochester, Rochester, New York 146422; Department of Medicine, Division of Nephrology, University of Rochester, Rochester, New York 146423; and Theoretical Biology and Biophysics, Los Alamos National Laboratory, Los Alamos, New Mexico 875454 Received 4 February 2010/Accepted 8 April 2010
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. 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
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
* 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:
[email protected]. 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:
[email protected]. 䌤 Published ahead of print on 21 April 2010. 6687
6688
MIAO ET AL.
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 virusspecific 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
J. VIROL.
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 ⫻ 105 50% egg infective doses (EID50) 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 (EID50)/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 ⫻ 106 antigen-presenting cells were added to 1 ⫻ 106 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.
VOL. 84, 2010
QUANTIFYING IMMUNE RESPONSE TO INFLUENZA INFECTION
6689
TABLE 1. Parameter and variable definitions for preadaptive and adaptive immune response models Parameter or variable
V(0) Ep(0) E*p(0) Ep E*p V TE AG AM ␣ E ␣ ⫽ ␣E kE ␦E* ␣ ⫽ V/␣ cV kVG kVM
Definition (unit)
Value
Initial infective viral titer (EID50/ml) per lung Initial no. of uninfected and infectible epithelial cells per lung Initial no. of infected epithelial cells per lung Uninfected and infectible epithelial cells per lung Infected epithelial cells per lung Infective viral titer (EID50/ml) Effector CD8⫹ T cells per lung Serum IgG antibody (pg/ml) Serum IgM antibody (pg/ml) A constant for unit standardization Net growth rate of uninfected epithelial cells (day⫺1) ␣ ⫻ the infection rate of uninfected epithelial cells 关ml 䡠 (EID50 䡠 day)⫺1兴 Depletion (killing) rate of infected epithelial cell by CD8⫹ effector cells (day⫺1) Death rate of infected epithelial cells (day⫺1) Production rate of infectious virus per infected epithelial cell divided by unknown constant ␣ (EID50 䡠 cell⫺1 䡠 day⫺1) Virus clearance rate (reduction of infectious viral titer) (day⫺1) Depletion (neutralization) rate of virus (infective viral titer) by IgG antibody 关ml/(pg 䡠 day)兴 Depletion (neutralization) rate of virus (infective viral titer) by IgM antibody 关ml/(pg 䡠 day)兴
Estimated Estimated Fixed (equals 0) Calculated Calculated Measured Measured Measured Measured Unknown Estimated Estimated Estimated Estimated Fixed
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,
冦
d E ⫽ EE p ⫺  ␣E PV dt p d E* ⫽ ␣EPV ⫺ kEE*PTE共t兲 ⫺ ␦EⴱE*p dt p d V ⫽ ␣E*P ⫺ cVV ⫺ kVGVAG共t兲 ⫺ kVMVAM共t兲, dt
with initial conditions
再
Ep共0兲 ⫽ Ep0 * Ep*共0兲 ⫽ EP0 V共0兲 ⫽ V0
(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 Ep have a growth rate E, and viral infection occurs with a rate ␣. In the second equation, infected epithelial cells Ep* are generated due to viral infection with the term ␣EPV and eliminated by CD8⫹ T effectors TE with a rate kE or by cell death (virus induced or apoptotic) with a constant rate ␦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 cV or neutralized by virus-specific IgG (AG) and IgM (AM) antibodies, with rates kVG and kVM, respectively. We measured the antigen-specific CD8⫹ T-cell counts (CTLs) in the lung,
FIG. 1. Schematic diagram of the differential equation models for influenza A virus infection in the lung.
Estimated Estimated Estimated
IAV-specific serum IgG and IgM antibody concentrations, and viral titers [denoted by TE(t), AG(t), AM(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. 2D), 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:
冦
d E ⫽ EE p ⫺  ␣E PV dt p d E* ⫽ ␣EPV ⫺ ␦EⴱE*P, with initial conditions dt p d V ⫽ ␣E*P ⫺ cVV dt
再
Ep共0兲 ⫽ EP0 * Ep*共0兲 ⫽ EP0 . V共0兲 ⫽ V0
(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 Ep(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, log10Y(ti) ⫽ V(ti) ⫹ ε(ti), where i equals 1, 2, . . . n; V(ti) denotes the true viral titer value at time point ti; Yi denotes the measured value of V(ti); and ε(ti) denotes the independently and identically distributed (i.i.d.) measurement errors. Note that log transformation is used in the measurement model to stabilize the computational algorithms. The standard nonlinear leastsquares (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
6690
MIAO ET AL.
J. VIROL.
FIG. 2. Experimental data. (A) Viral titer; (B) antigen-specific CD8⫹ T-cell data and the smoothed curve; (C) IgG data and the smoothed curve; (D) IgM data and the smoothed curve. 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:
冉 冊
RSS AICc ⫽ n ln n
2nk ⫹ n⫺k⫺1
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 log10(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, Ep*(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 ˆ (0) ⬇ 103.1682 ⬇ 1473 EID50/ml. condition of viral load was as follows: V 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 virusspecific 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
VOL. 84, 2010
QUANTIFYING IMMUNE RESPONSE TO INFLUENZA INFECTION Response
Early
Adaptive
Overall
6.2E⫺08 (3.9E⫺09, 9.6E⫺08) 3.4E⫺01 (7.6E⫺02, 9.5E⫺01) 8.6E⫺04 (4.8E⫺06, 1.00E⫹00)
E (day⫺1)
2.4E⫺06 (3.8E⫺07, 1.1E⫺05) 7.8E⫺01 (5.8E⫺02, 9.4E⫺01) 5.1E⫺06 (9.1E⫺07, 1.7E⫺05)
␣ (ml 䡠 EID50⫺1 䡠 day⫺1) kE (cell⫺1 day⫺1)
6.0E⫺01 (4.5E⫺02, 2.7E⫹01) 1.5E⫺01 (4.5E⫺05, 1.2E⫹00) 1.2E⫹00 (7.4E⫺01, 1.6E⫹00)
␦E* (day⫺1)
4.2E⫹00 (4.3E⫺01, 1.2E⫹02) 4.0E⫺04 (1.1E⫺04, 8.1E⫺02) 5.5E⫺06 (2.2E⫺06, 5.3E⫺01)
cV (day⫺1)
3.1E⫺05 (5.3E⫺07, 8.6E⫺05) 4.5E⫺08 (4.7E⫺09, 8.1E⫺06)
kVG 关ml/(pg 䡠 day)兴
4.4E⫹00 (1.2E⫹00, 10E⫹00) 7.8E⫺02 (5.1E⫺02, 8.1E⫹00)
kVM 关ml/(pg 䡠 day)兴
Estimates of kinetic parameters
6.4E⫺05 (2.0E⫺05, 7.6E⫺05) 1.4E⫺05 (1.1E⫺10, 5.0E⫺05)
TABLE 2. Estimates of kinetic parameters for the early immune response, the adaptive immune response, and the combined two phasesa
Ep(0) (cells per lung)
5.8E⫹05 (1.1E⫹05, 2.3E⫹07)
2.4E⫹05 (6.4E⫹04, 10E⫹06)
ml 䡠 (EID50)⫺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. 3A). 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 104 EID50/(ml 䡠 day) per infected cell does not change estimates of other parameters, such as the
a The early immune response occurs during the first 5 days (model 2), the adaptive immune response occurs during days 5 to 14 of infection (model 1), with ␣ equaling 100 EID50/(ml 䡠 day 䡠 cell), and the combined two phases occur during days 0 to 14 (model 1), with ␣ equaling 100 EID50/(ml 䡠 day 䡠 cell). The 95% confidence intervals are listed in parentheses after each estimate.
FIG. 3. Fitted viral titer curves for the innate and adaptive immune responses and corresponding predicted epithelial cell dynamics. (A) Fit of model 2 to the viral titer data from day 0 to day 5; (B) fit of model 1 to the viral titer data from day 5 to day 14; (C) fit of model 1 to the viral titer data from day 0 to day 14. Viral titer data are replotted from Fig. 2A. Best-fit parameter values used to generate the theoretical curves are given in Table 2.
6691
6692
MIAO ET AL.
J. VIROL.
FIG. 4. Comparison of parameter estimates in preadaptive (model 2 in blue), adaptive (model 1 in red), and overall (model 1 in black) immune responses and relative contributions of different components.
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, Ep(t ⫽ 0), and therefore with the number of uninfected (Ep) and IAV-infected (Ep*) epithelial cells during infection. That is, the effects of ␣ are compensated by the corresponding changes in Ep and Ep*, such that the magnitude of Ep and Ep*changes for different values of ␣. Because of this balance, the estimates of all other parameters are essentially unaffected by ␣. Finally, it should be noted that ␣ does have ˆ (t), e.g., the peak] if the a significant effect on the viral titer [V compensation by Ep is removed; that is, if Ep(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 Ep(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 ˆ (t), although it has no effect on other model on viral titer V parameter estimates when Ep(0) is not fixed. 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 Ep(t ⫽ 5) and infected Ep*(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 104 EID50/(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 䡠 (EID50)⫺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. 4C) 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 䡠 (EID50)⫺1 䡠 day⫺1, the CTL killing rate is 1.4 ⫻ 10⫺5
VOL. 84, 2010
QUANTIFYING IMMUNE RESPONSE TO INFLUENZA INFECTION
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 (Ep) rapidly drops to nearly zero by day 2 (see Fig. 3), such that the two terms EEp and ␣EPV in the first equation of model 1 become negligible during the adaptive immune response phase, and therefore, the estimates of 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 EID50/(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).
6693
TABLE 3. Model comparison and selection for the adaptive phase and the combined two phases Phase
Modela
Adaptive
Full model 1 kE ⫽ 0, kVG ⫽ 0, kVM ⫽ 0 E ⫽ 0, ␣ ⫽ 0 E ⫽ 0, ␣ ⫽ 0, cV ⫽ 0 E ⫽ 0, ␣ ⫽ 0, kE ⫽ 0 E ⫽ 0, ␣ ⫽ 0, ␦E* ⫽ 0 E ⫽ 0, ␣ ⫽ 0, kVG ⫽ 0 E ⫽ 0, ␣ ⫽ 0, kVM ⫽ 0 E ⫽ 0, ␣ ⫽ 0, kVG ⫽ 0, kVM ⫽ 0 E ⫽ 0, ␣ ⫽ 0, kVG ⫽ 0, cV ⫽ 0 E ⫽ 0, ␣ ⫽ 0, kVG ⫽ 0, cV ⫽ 0, ␦E* ⫽ 0
11.7 24.3 11.7 11.7 15.9 11.8 11.7 12.4 19.4 11.7 11.8
Overall
Full model 1 ␣ ⫽ 0 ␦E* ⫽ 0, kE ⫽ 0, kVM ⫽ 0 E ⫽ 0 E ⫽ 0, cV ⫽ 0 E ⫽ 0, kVG ⫽ 0 E ⫽ 0, cV ⫽ 0, kVG ⫽ 0 E ⫽ 0, cV ⫽ 0, kVG ⫽ 0, ␦E* ⫽ 0 E ⫽ 0, cV ⫽ 0, kVG ⫽ 0, kVM ⫽ 0 E ⫽ 0, cV ⫽ 0, kVG ⫽ 0, kE ⫽ 0
49.4 2738 651 49.4 49.4 49.4 49.4 885 486 50.0
a
RSS
AICc score
⫺9.24 2.69 ⫺16.2 ⫺19.2 ⫺9.99 ⫺19.1 ⫺19.2 ⫺17.6 ⫺7.15 ⴚ22.2 ⴚ24.1 ⫺33.2 315 185 ⫺34.2 ⫺36.2 ⫺36.2 ⴚ39.2 210 158 ⴚ40.1
The two best models for each phase are in boldface.
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 100 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 simu-
6694
MIAO ET AL.
FIG. 5. Comparison of viral titers between the shifted IgM data and the original IgM data from days 5 to 14 based on the prediction of model 2 (␣ ⫽ 100 EID50 䡠 ml⫺1 䡠 cell⫺1 䡠 day⫺1).
lated 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 (kE) 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 (kVM) 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. 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 䡠 EID50⫺1 䡠 day⫺1, the magnitude of viral titer peak essentially does not change. Instead, the peak time moves closer to day 0 (Fig. 6A). This simulation result suggests that the infec-
J. VIROL.
tion 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 Ep(t ⫽ 0)] and increase only the viral replication rate ␣ from 10 to 103 EID50 䡠 ml⫺1 䡠 cell⫺1 䡠 day⫺1, the peak viral titer increases from 106 to 108 EID50. Further increase of ␣ from 103 to 104 EID50 䡠 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 Ep(t ⫽ 0) from 105 to 108 cells while all other parameters are held constant, the peak viral titer increases from 106 to 109 EID50/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 Ep(t ⫽ 0) is as low as 104 or 103 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 (IAVspecific 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
VOL. 84, 2010
QUANTIFYING IMMUNE RESPONSE TO INFLUENZA INFECTION
6695
FIG. 6. Variation in infection patterns of epithelial cells and resulting viral titers predicted by fitting model 2 to the early viral titer data (days 0 to 5). Note that the data and simulation results are plotted through the adaptive phase (days 5 to 14) as well. (A) Simulation results for models fit to data and different influenza virus infection rates (␣) from 10⫺6 to 1.0 ml 䡠 (EID50 䡠 day)⫺1 (fixed ␣ ⫽ 100 EID50 䡠 ml⫺1 䡠 cell⫺1 䡠 day⫺1). There is little variation over a 7-log range of ␣ in uninfected epithelial cells (Ep) and infected epithelial cells (E*P). However, viral titer simulations correlate with measured data only during the first 5 days of infection, which is consistent with the fact that the model does not take into account the changes occurring during the adaptive phase (e.g., IgM clearance of virus and CTL killing of infected epithelial cells). Circles in the plot of viral titers represent actual data. (B) Simulation results for different viral production rates. Note that early viral titers vary considerably with different viral production rates. Please see the text for detailed discussions. Parameter values used to generate the theoretical curves are given in Table 2.
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
FIG. 7. Change of viral titer for different initial numbers of uninfected epithelial cells Ep(0) from 103 to 108 cells during the early immune response using model 2 (fixed ␣ ⫽ 100 EID50 䡠 ml⫺1 䡠 cell⫺1 䡠 day⫺1; other parameters are given in Table 2).
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 avail-
6696
MIAO ET AL.
ability 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. d E ⫽ EE p ⫺  ␣E PV dt p 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 EEp. Once the infection is established, estimating the term for epithelial proliferation becomes unreliable, primarily because the number of infectible target cells, Ep, is nearly zero. Even if the proliferation rate of the epithelial cells increases, the value of the entire term (EEp) is still low once most of the target cells are infected (e.g., Ep ⬇ 0) and has no effect on the adaptive model predictions. Similar analysis holds for estimation of the number of viral infections per unit time (␣). Essentially, a very small number of target cells (Ep) and a large amount of virus (V) cannot be used to reliably estimate the epithelial proliferation (E) or infection rates (␣) because most of the target cells have become infected (Ep*). 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 (kE) of 6.4 ⫻ 10⫺5 per CTL per day. The antigen-specific CD8 T-cell number peaks at ⬎5 ⫻ 104 cells, yielding an effective peak killing rate (kE 䡠 TE) of ⬎3.2 day⫺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, kE equaling 6.4 ⫻ 10⫺5 per cell per day and kE 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
J. VIROL.
cannot be cleared (instead, the virus titer reaches a nadir of about 3.5 log10 EID50 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 IgMmediated clearance of free virus (kVM) 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 (kVG) of 1.0 ⫻ 10⫺5 ml/(pg 䡠 day). During the adaptive phase, our estimate of the clearance rate of virus particles due to factors other than antibody neutralization (cV) 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 cV at the innate level when estimating the rates of virus inactivation by IgG and IgM (kVG and kVM, respectively), the estimates nearly do not change, indicating that the overall model is insensitive to changes in cV. 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⫹ 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
VOL. 84, 2010
QUANTIFYING IMMUNE RESPONSE TO INFLUENZA INFECTION
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 䡠 EID50⫺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 (TCID50/ 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 EID50 䡠 ml⫺1 䡠 cell⫺1 䡠 day⫺1 and estimated innate virus clearance (cV) 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 (TCID50/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.
6697
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.). REFERENCES 1. Akaike, H. 1973. Information theory and an extension of the maximum likehood principle, p. 267–281. In B. Petrov and F. Csa´ki (ed.), Second International Symposium on Information Theory. Akademiai Kiado, Budapest, Hungary. 2. Allan, W., Z. Tabi, A. Cleary, and P. C. Doherty. 1990. Cellular events in the lymph node and lung of mice with influenza. Consequences of depleting CD4⫹ T cells. J. Immunol. 144:3980–3986. 3. Baccam, P., C. Beauchemin, C. A. Macken, F. G. Hayden, and A. S. Perelson. 2006. Kinetics of influenza A virus infection in humans. J. Virol. 80:7590– 7599. 4. Baumgarth, N., O. C. Herman, G. C. Jager, L. E. Brown, L. A. Herzenberg, and J. Chen. 2000. B-1 and B-2 cell-derived immunoglobulin M antibodies are nonredundant components of the protective response to influenza virus infection. J. Exp. Med. 192:271–280. 5. Beauchemin, C., J. Samuel, and J. Tuszynski. 2005. A simple cellular automaton model for influenza A viral infections. J. Theor. Biol. 232:223–234. 6. Benali, R., J. M. Tournier, M. Chevillard, J. M. Zahm, J. M. Klossek, J. Hinnrasky, D. Gaillard, F. X. Maquart, and E. Puchelle. 1993. Tubule formation by human surface respiratory epithelial cells cultured in a threedimensional collagen lattice. Am. J. Physiol. 264:L183–L192. 7. Bocharov, G. A., and A. A. Romanyukha. 1994. Mathematical model of antiviral immune response. III. Influenza A virus infection. J. Theor. Biol. 167:323–360. 8. Brown, D. M., E. Roman, and S. L. Swain. 2004. CD4 T cell responses to influenza infection. Semin. Immunol. 16:171–177. 9. Burnham, K. P., and D. R. Anderson. 2004. Multimodel inference: understanding AIC and BIC in model selection. Sociol. Methods Res. 33:261–304. 10. Cairns, H. J. F., and M. Edney. 1952. Quantitative aspects of influenza virus multiplication. I. Production of “incomplete” virus. J. Immunol. 69:155–160. 11. Doherty, P. C., D. J. Topham, R. A. Tripp, R. D. Cardin, J. W. Brooks, and P. G. Stevenson. 1997. Effector CD4⫹ and CD8⫹ T-cell mechanisms in the control of respiratory virus infections. Immunol. Rev. 159:105–117. 12. Eichelberger, M., W. Allan, M. Zijlstra, R. Jaenisch, and P. C. Doherty. 1991. Clearance of influenza virus respiratory infection in mice lacking class I major histocompatibility complex-restricted CD8⫹ T cells. J. Exp. Med. 174:875–880. 13. Ganusov, V. V., and R. J. De Boer. 2008. Estimating in vivo death rates of targets due to CD8 T-cell-mediated killing. J. Virol. 82:11749–11757. 14. Graham, M. B., and T. J. Braciale. 1997. Resistance to and recovery from lethal influenza virus infection in B lymphocyte-deficient mice. J. Exp. Med. 186:2063–2068. 15. Hancioglu, B., D. Swigon, and G. Clermont. 2007. A dynamical model of human immune response to influenza A virus infection. J. Theor. Biol. 246:70–86. 16. Handel, A., I. M. Longini, Jr., and R. Antia. 2007. Neuraminidase inhibitor resistance in influenza: assessing the danger of its generation and spread. PLoS Comput. Biol. 3:e240. 17. Handel, A., A. Yates, S. S. Pilyugin, and R. Antia. 2009. Sharing the burden: antigen transport and firebreaks in immune responses. J. R. Soc. Interface 6:447–454. 18. Henn, A. D., J. Rebhahn, M. A. Brown, A. J. Murphy, M. N. Coca, O. Hyrien, T. Pellegrin, T. Mosmann, and M. S. Zand. 2009. Modulation of single-cell IgG secretion frequency and rates in human memory B cells by CpG DNA, CD40L, IL-21, and cell division. J. Immunol. 183:3177–3187. 19. Ho, D., A. Neumann, A. Perelson, W. Chen, J. Leonard, and M. Markowitz. 1995. Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection. Nature 373:123–126. 20. Hou, S., P. C. Doherty, M. Zijlstra, R. Jaenisch, and J. M. Katz. 1992. Delayed clearance of Sendai virus in mice lacking class I MHC-restricted CD8⫹ T cells. J. Immunol. 149:1319–1325. 21. Hou, S., X. Y. Mo, L. Hyland, and P. C. Doherty. 1995. Host response to Sendai virus in mice lacking class II major histocompatibility complex glycoproteins. J. Virol. 69:1429–1434. 22. Jelley-Gibbs, D. M., J. P. Dibble, S. Filipson, L. Haynes, R. A. Kemp, and S. L. Swain. 2005. Repeated stimulation of CD4 effector T cells can limit their protective function. J. Exp. Med. 201:1101–1112. 23. Lee, B. O., J. Rangel-Moreno, J. E. Moyron-Quiroz, L. Hartson, M. Makris, F. Sprague, F. E. Lund, and T. D. Randall. 2005. CD4 T cell-independent antibody response promotes resolution of primary influenza infection and helps to prevent reinfection. J. Immunol. 175:5827–5838. 24. Lee, H. Y., D. J. Topham, S. Y. Park, J. Hollenbaugh, J. Treanor, T. R. Mosmann, X. Jin, B. M. Ward, H. Miao, J. Holden-Wiltse, A. S. Perelson, M. Zand, and H. Wu. 2009. Simulation and prediction of the adaptive immune response to influenza A virus infection. J. Virol. 83:7151–7165.
6698
MIAO ET AL.
25. Liang, H., H. Miao, and H. Wu. 2010. Estimation of constant and timevarying dynamic parameters of HIV infection in a nonlinear differential equation model. Ann. App. Stat. 4:460–483. 26. Liang, H., and H. Wu. 2008. Parameter estimation for differential equation models using a framework of measurement error in regression model. J. Am. Stat. Assoc. 103:1570–1583. 27. Look, D. C., M. J. Walter, M. R. Williamson, L. Pang, Y. You, J. N. Sreshta, J. E. Johnson, D. S. Zander, and S. L. Brody. 2001. Effects of paramyxoviral infection on airway epithelial cell Foxj1 expression, ciliogenesis, and mucociliary function. Am. J. Pathol. 159:2055–2069. 28. Ma, S., and M. R. Kosorok. 2005. Robust semiparametric M-estimation and the weighted bootstrap. J. Multivar. Anal. 96:190–217. 29. Mestas, J., and C. C. W. Hughes. 2004. Of mice and not men: differences between mouse and human immunology. J. Immunol. 172:2731–2738. 30. Miao, H., X. Xia, A. S. Perelson, and H. Wu. On identifiability of nonlinear ODE models with applications in viral dynamics. SIAM Rev. Soc. Ind. Appl. Math., in press. 31. Mohler, L., D. Flockerzi, H. Sann, and U. Reichl. 2005. Mathematical model of influenza A virus production in large-scale microcarrier culture. Biotechnol. Bioeng. 90:46–58. 32. Mozdzanowska, K., M. Furchner, K. Maiese, and W. Gerhard. 1997. CD4⫹ T cells are ineffective in clearing a pulmonary infection with influenza type A virus in the absence of B cells. Virology 239:217–225. 33. Mozdzanowska, K., K. Maiese, and W. Gerhard. 2000. Th cell-deficient mice control influenza virus infection more effectively than Th- and B cell-deficient mice: evidence for a Th-independent contribution by B cells to virus clearance. J. Immunol. 164:2635–2643. 34. Palladino, G., K. Mozdzanowska, G. Washko, and W. Gerhard. 1995. Virusneutralizing antibodies of immunoglobulin G (IgG) but not of IgM or IgA isotypes can cure influenza virus pneumonia in SCID mice. J. Virol. 69:2075– 2081. 35. Perelson, A., A. Neumann, M. Markowitz, J. Leonard, and D. Ho. 1996. HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science 271:1582–1586. 36. Perelson, A. S. 2002. Modelling viral and immune system dynamics. Nat. Rev. Immunol. 2:28–36. 37. Perelson, A. S., P. Essunger, Y. Cao, M. Vesanen, A. Hurley, K. Saksela, M. Markowitz, and D. D. Ho. 1997. Decay characteristics of HIV-1-infected compartments during combination therapy. Nature 387:188–191. 38. Perelson, A. S., and P. Nelson. 1999. Mathematical analysis of HIV-1 dynamics in vivo. SIAM Rev. Soc. Ind. Appl. Math. 41:3–44. 39. Polakos, N. K., J. C. Cornejo, D. A. Murray, K. O. Wright, J. J. Treanor, I. N. Crispe, D. J. Topham, and R. H. Pierce. 2006. Kupffer cell-dependent hepatitis occurs during influenza infection. Am. J. Pathol. 168:1169–1178. 40. Ramphal, R., R. C. Cogliano, J. W. Shands, Jr., and P. A. Small, Jr. 1979. Serum antibody prevents lethal murine influenza pneumonitis but not tracheitis. Infect. Immun. 25:992–997. 41. Rawlins, E. L., and B. L. Hogan. 2008. Ciliated epithelial cell lifespan in the mouse trachea and lung. Am. J. Physiol. Lung Cell. Mol. Physiol. 295:L231– L234. 42. Regoes, R. R., D. L. Barber, R. Ahmed, and R. Antia. 2007. Estimation of the rate of killing by cytotoxic T lymphocytes in vivo. Proc. Natl. Acad. Sci. U. S. A. 104:1599–1603. 43. Regoes, R. R., A. Yates, and R. Antia. 2007. Mathematical models of cytotoxic T-lymphocyte killing. Immunol. Cell Biol. 85:274–279. 44. Rodriguez-Fernandez, M., J. A. Egea, and J. R. Banga. 2006. Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems. BMC Bioinformatics 7:483.
J. VIROL. 45. Roman, E., E. Miller, A. Harmsen, J. Wiley, U. H. Von Andrian, G. Huston, and S. L. Swain. 2002. CD4 effector T cell subsets in the response to influenza: heterogeneity, migration, and function. J. Exp. Med. 196:957–968. 46. Sangster, M. Y., J. M. Riberdy, M. Gonzalez, D. J. Topham, N. Baumgarth, and P. C. Doherty. 2003. An early CD4⫹ T cell-dependent immunoglobulin A response to influenza infection in the absence of key cognate T-B interactions. J. Exp. Med. 198:1011–1021. 47. Sarawar, S. R., M. Sangster, R. L. Coffman, and P. C. Doherty. 1994. Administration of anti-IFN-gamma antibody to beta 2-microglobulin-deficient mice delays influenza virus clearance but does not switch the response to a T helper cell 2 phenotype. J. Immunol. 153:1246–1253. 48. Smith, A., F. Adler, and A. Perelson. 2010. An accurate two-phase approximate solution to an acute viral infection model. J. Math. Biol. 60:711–726. 49. Stafford, M. A., L. Corey, Y. Cao, E. S. Daar, D. D. Ho, and A. S. Perelson. 2000. Modeling plasma virus concentration during primary HIV infection. J. Theor. Biol. 203:285–301. 50. Stray, S. J., and G. M. Air. 2001. Apoptosis by influenza viruses correlates with efficiency of viral mRNA synthesis. Virus Res. 77:3–17. 51. Subbarao, K., J. McAuliffe, L. Vogel, G. Fahle, S. Fischer, K. Tatti, M. Packard, W. J. Shieh, S. Zaki, and B. Murphy. 2004. Prior infection and passive transfer of neutralizing antibody prevent replication of severe acute respiratory syndrome coronavirus in the respiratory tract of mice. J. Virol. 78:3572–3577. 52. Topham, D. J., R. A. Tripp, and P. C. Doherty. 1997. CD8⫹ T cells clear influenza virus by perforin or Fas-dependent processes. J. Immunol. 159: 5197–5200. 53. Topham, D. J., R. A. Tripp, A. M. Hamilton-Easton, S. R. Sarawar, and P. C. Doherty. 1996. Quantitative analysis of the influenza virus-specific CD4⫹ T cell memory in the absence of B cells and Ig. J. Immunol. 157:2947–2952. 54. Topham, D. J., R. A. Tripp, S. R. Sarawar, M. Y. Sangster, and P. C. Doherty. 1996. Immune CD4⫹ T cells promote the clearance of influenza virus from major histocompatibility complex class II ⫺/⫺ respiratory epithelium. J. Virol. 70:1288–1291. 55. Tripp, R. A., S. R. Sarawar, and P. C. Doherty. 1995. Characteristics of the influenza virus-specific CD8⫹ T cell response in mice homozygous for disruption of the H-2lAb gene. J. Immunol. 155:2955–2959. 56. Wei, X., S. Ghosh, M. Taylor, V. Johnson, E. Emini, P. Deutsch, J. Lifson, S. Bonhoeffer, M. Nowak, and B. Hahn. 1995. Viral dynamics in human immunodeficiency virus type 1 infection. Nature 373:117–122. 57. Wikenheiser, K. A., D. K. Vorbroker, W. R. Rice, J. C. Clark, C. J. Bachurski, H. K. Oie, and J. A. Whitsett. 1993. Production of immortalized distal respiratory epithelial cell lines from surfactant protein C/simian virus 40 large tumor antigen transgenic mice. Proc. Natl. Acad. Sci. U. S. A. 90: 11029–11033. 58. Xia, X., and C. H. Moog. 2003. Identifiability of nonlinear systems with application to HIV/AIDS models. IEEE Trans. Automat. Contr. 48:330–336. 59. Xue, H., H. Miao, and H. Wu. Sieve estimation of constant and time-varying coefficients in nonlinear ordinary differential equation models by considering both numerical error and measurement error. Ann. Stat., in press. http:// www.imstat.org/aos/future_papers.html. 60. Yates, A., F. Graw, D. L. Barber, R. Ahmed, R. R. Regoes, and R. Antia. 2007. Revisiting estimates of CTL killing rates in vivo. PLoS One 2:e1301. 61. Zheng, B. J., K. W. Chan, Y. P. Lin, G. Y. Zhao, C. Chan, H. J. Zhang, H. L. Chen, S. S. Wong, S. K. Lau, P. C. Woo, K. H. Chan, D. Y. Jin, and K. Y. Yuen. 2008. Delayed antiviral plus immunomodulator treatment still reduces mortality in mice infected by high inoculum of influenza A/H5N1 virus. Proc. Natl. Acad. Sci. U. S. A. 105:8091–8096.
