Influenza Virus Infection Model With Density Dependence Supports Biphasic Viral Decay

Mathematical models that describe infection kinetics help elucidate the time scales, effectiveness, and mechanisms underlying viral growth and infection resolution. For influenza A virus (IAV) infections, the standard viral kinetic model has been used to investigate the effect of different IAV proteins, immune mechanisms, antiviral actions, and bacterial coinfection, among others. We sought to further define the kinetics of IAV infections by infecting mice with influenza A/PR8 and measuring viral loads with high frequency and precision over the course of infection. The data highlighted dynamics that were not previously noted, including viral titers that remain elevated for several days during mid-infection and a sharp 4–5 log₁₀ decline in virus within 1 day as the infection resolves. The standard viral kinetic model, which has been widely used within the field, could not capture these dynamics. Thus, we developed a new model that could simultaneously quantify the different phases of viral growth and decay with high accuracy. The model suggests that the slow and fast phases of virus decay are due to the infected cell clearance rate changing as the density of infected cells changes. To characterize this model, we fit the model to the viral load data, examined the parameter behavior, and connected the results and parameters to linear regression estimates. The resulting parameters and model dynamics revealed that the rate of viral clearance during resolution occurs 25 times faster than the clearance during mid-infection and that small decreases to this rate can significantly prolong the infection. This likely reflects the high efficiency of the adaptive immune response. The new model provides a well-characterized representation of IAV infection dynamics, is useful for analyzing and interpreting viral load dynamics in the absence of immunological data, and gives further insight into the regulation of viral control.

1. Introduction

Influenza A virus (IAV) is a leading cause of lower respiratory tract infections and causes a significant amount of morbidity and mortality (Simonsen et al., 2000; Taubenberger and Morens, 2008; Medina and García-Sastre, 2011), with over 15 million individuals infected and more than 200,000 hospitalizations each year in the U.S. (Thompson et al., 2004). Vaccination against influenza viruses remains the most effective measure to prevent infection, but the large number of antigenically distinct strains, the emergence of new strains, and the low efficacy of antivirals make combatting the disease challenging. New therapeutic strategies are thus necessary and may require modulation of different viral control mechanisms, which are not entirely understood for IAV infection. Thus, it is critical to gain a deeper understanding of the infection kinetics, including determining the time scales, magnitudes, contribution, and interrelatedness of different control processes throughout IAV infection.

Kinetic modeling of in vivo infection processes provides important insight into viral growth and decay, host immune responses, antiviral actions, and multi-pathogen interactions. Remarkably, as few as 3–4 equations for target cells, infected cells, and virus can accurately describe viral load dynamics for a variety of virus infections [e.g., IAV, HIV, HCV, Zika virus, and West Nile Virus (Perelson et al., 1996; Neumann et al., 1998; Baccam et al., 2006; Banerjee et al., 2016; Best et al., 2017)]. For IAV infections, numerous studies have used these simple models with great success to elucidate mechanisms during IAV infection and during IAV coinfection with bacterial pathogens (reviewed in Smith and Ribeiro, 2010; Beauchemin and Handel, 2011; Smith and Perelson, 2011; Smith and McCullers, 2014; Boianelli et al., 2015). However, investigating mechanisms of immune control is often inhibited by insufficient data, which limits effective model calibration and selection. Further, it can be difficult to distinguish between mechanisms because a viral kinetic model that excludes equations and terms for specific immune responses can fit viral load dynamics with ease.

To aid interpretation of model results and gain insight into the mechanisms of infection, previous studies have used linear regression and approximate solutions to the viral kinetic model (derived by Smith et al., 2010) to identify how different processes (e.g., virus infection, production, and clearance) contribute to viral load dynamics throughout the course of infection (Miao et al., 2010; Smith et al., 2010, 2011a; Holder et al., 2011a,b; Halloran et al., 2012; Li and Handel, 2014; Kakizoe et al., 2015; Pinky and Dobrovolny, 2016; Best et al., 2017; Palmer et al., 2017; Smith, 2017). In the initial hours of infection, virus quickly infects cells or is cleared. Following an eclipse phase, virus production begins and virus increases exponentially for ~2 d. This initial growth can be approximated by a linear function of the log₁₀ viral titers or by V(t) = e^λt, where λ is a combination of all infection processes and is equivalent to the log-linear slope (Smith et al., 2010, 2011a). After this growth phase, virus peaks and begins to decline until the infection is resolved. Virus decay is typically exponential in nature and can be approximated in a similar fashion as the growth phase. That is, V(t) = e ^{− δt}, where δ is the infected cell death rate and the sole process dictating the viral decay dynamics. Here, the log-linear slope is an estimate of the infected cell death rate (Smith et al., 2010, 2011a).

Although these dynamics and approximations have improved our knowledge of viral kinetics, some dynamical features, such as the plateauing of virus following the peak (reviewed in Smith and Perelson, 2011) cannot be explained by current kinetic models that exclude equations for immune factors. One model could reproduce the plateauing of virus through modeling interferon and an interferon-induced adaptive immune response (Pawelek et al., 2012). The study concluded that specific equations for the innate and adaptive responses were necessary. However, quantitative immunological data was not used to support model selection, parameterization, and conclusions. This type of data is scarce and has been a limiting factor of modeling studies. With viral loads as the most prevalent type of data, models that limit the number of parameters and equations remain desirable. However, even most viral load data is insufficiently quantitative to confidently detect features like a mid-infection plateau and build appropriate mathematical models.

Here, we first sought to increase the quality and quantity of viral load data in order to improve predictive power of mathematical models and gain a deeper insight into the kinetics of viral resolution. To do this, we measured viral loads daily from groups of BALB/cJ mice infected with influenza A/Puerto Rico/8/34 (H1N1) (PR8). In addition, we tightly controlled the experimental conditions and repeated the experiment numerous times to ensure reproducibility and identify data with meaningful biological heterogeneity (i.e., due to an underlying mechanism) vs. data with experimental heterogeneity (i.e., due to poor technique). The high resolution of these data defined important dynamical features, including a long plateau phase followed by a rapid decay phase. Because current viral kinetic models cannot reproduce these data, we developed a new model that incorporated a density-dependent decay of infected cells and could accurately describe the observed viral load dynamics. We used a rigorous fitting scheme to estimate the model parameters and infer important dynamics. Subsequent linear regression analysis and sensitivity analysis aided effective interpretation of the model results and direct comparison with published results. The data, model, and analyses provide a robust quantification of IAV infection kinetics and indicate that the rate of virus clearance changes with respect to the density of infected cells.

2. Materials and Methods

2.1. Use of Experimental Animals

All experimental procedures were approved by the Animal Care and Use Committee at SJCRH under relevant institutional and American Veterinary Medical Association guidelines and were performed in a Biosafety level 2 facility that is accredited by AALAAS.

2.2. Mice

Adult (6 week old) female BALB/cJ mice were obtained from Jackson Laboratories (Bar Harbor, ME). Mice were housed in groups of 5 mice in high-temperature 31.2 × 23.5 × 15.2 cm polycarbonate cages with isolator lids. Rooms used for housing mice were maintained on a 12:1 2-h light:dark cycle at 22 ± 2°C with 50% humidity in the biosafety level 2 facility at St. Jude Children's Research Hospital (Memphis, TN). Prior to inclusion in the experiments, mice were allowed at least 7 days to acclimate to the animal facility such that they were 7 weeks old at the time of infection. Laboratory Autoclavable Rodent Diet (PMI Nutrition International, St. Louis, MO) and autoclaved water were available ad libitum. All experiments were performed under an approved protocol and in accordance with the guidelines set forth by the Animal Care and Use Committee at St. Jude Children's Research Hospital.

2.3. Infectious Agents

All experiments were done using the mouse adapted influenza A/Puerto Rico/8/34 (H1N1) (PR8).

2.4. Infection Experiments

The viral infectious dose (TCID₅₀) was determined by interpolation using the method of Reed and Muench (Reed and Muench, 1938) using serial dilutions of virus on Madin-Darby canine kidney (MDCK) cells. Mice were intranasally inoculated with 75 TCID₅₀ PR8 in 100 μl. Mice were weighed at the onset of infection and each subsequent day for illness and mortality. Mice were euthanized if they became moribund or lost 30% of their starting body weight. Each experiment was repeated numerous times to ensure reproducibility. Two complete experiments (10 animals per time point) were used for these studies. The raw data is shown in Figure 1A and is available upon request.

FIGURE 1

Figure 1. Phases of viral kinetics and fit of the density-dependent viral kinetic model. (A) Fit of the density-dependent viral kinetic model (Equations 1–5) to viral titers from the lungs of individual mice (10 mice per time point) infected with 75 TCID₅₀ PR8. Each dot is an individual mouse and the solid black line is the optimal solution of the model. Phases I–V of the viral kinetics are illustrated, where virus (I) quickly infects cells, (II) increases exponentially, (III) peaks, (IV) decays slowly, then (V) decays rapidly and clears. (B) Optimal fit of the model (solid black line) shown with the model solutions using parameter sets within the 95% CIs (gray shading). Data are shown as the mean ± standard deviation. Parameters are given in Table 1.

2.5. Lung Titers

Mice were euthanized by CO₂ asphyxiation. Lungs were aseptically harvested, washed three times in PBS, and placed in 500 μl sterile PBS. Whole lungs were digested with collagenase (1 mg/ml, Sigma C0130), and physically homogenized by syringe plunger against a 40 μm cell strainer. Cell suspensions were centrifuged at 4°C, 500 × g for 7 min. The supernatants were used to determine the viral titers using serial dilutions on MDCK monolayers.

2.6. Mathematical Model

The standard viral kinetic model used to describe IAV infection kinetics tracks 4 populations: susceptible epithelial (“target”) cells (T), two classes of infected cells (I₁ and I₂), and virus (V) (Baccam et al., 2006):

$\begin{array}{ll}\frac{dT}{dt}=-\beta TV& \left(1\right)\end{array}$

$\begin{array}{ll}\frac{d{I}_1}{dt}=\beta TV-k{I}_1& \left(2\right)\end{array}$

$\begin{array}{ll}\frac{d{I}_2}{dt}=k{I}_1-\delta (I_2)I_2& \left(3\right)\end{array}$

$\begin{array}{ll}\frac{dV}{dt}=p{I}_2-cV& \left(4\right)\end{array}$

In this model, target cells become infected with virus at rate βV per cell. Once infected, these cells enter an eclipse phase (I₁) before transitioning at rate k per cell to produce virus (I₂). Virus production occurs at rate p per cell. Virus is cleared at rate c and virus-producing infected cells (I₂) are cleared according to the function δ(I₂). The standard viral kinetic model assumes that infected cells are cleared at a constant rate (δ(I₂) = δ_s) (Baccam et al., 2006). The subscript s is used to denote “standard.” This model could not recapitulate the data (see Table S1 and Figure S1) and a modification of the model was necessary. Given that the rate of infected cell clearance (δ(I₂)) drives the virus decay dynamics (Smith et al., 2010), we let the clearance rate vary with the number of infected cells such that

$\begin{array}{ll}\delta (I_2)=\frac{{\delta }_d}{K_{\delta }+I_2},& \left(5\right)\end{array}$

where δ_d/K_δ is the maximum rate of infected cell clearance and K_δ is the half-saturation constant. The subscript d is used to denote “density-dependent.” Modifications to other terms were examined, but none could replicate the data.

2.7. Parameter Estimation

Given a parameter set θ, the cost $C(\theta )={\displaystyle \sum _{v_i}}{(V(\theta ,t_i)-v_i)}^2$ was minimized across parameter ranges using an Adaptive Simulated Annealing (ASA) global optimization algorithm (details in the Supplementary Material) to compare experimental and predicted values of log₁₀TCID₅₀/lung. A sample search pattern is shown in Figure S2. Errors of the log₁₀ data were assumed to be normally distributed. To explore and visualize the regions of parameters consistent with the models, we fit Equations 1–5 to 1,000 bootstrap replicates of the data. For each bootstrap data set, the model was fit 10 times beginning from the best-fit parameters estimate θ^best that was found by fitting the model to the data then perturbing each parameter estimate uniformly within ±50% of its best-fit value. If the three best bootstrap fits were within χ² = 0.05 of the best-fit, then the bootstrap was considered successful (Smith et al., 2011a, 2013). For each best fit estimate, we provide 95% confidence interval (CI) computed from the bootstrap replicates. All calculations were performed in MATLAB.

Estimated parameters included the rates of virus infection (β), virus production (p), virus clearance (c), eclipse phase transition (k), infected cell clearance (δ_d), and the half saturation constant (K_δ). Bounds were placed on the parameters to constrain them to physically realistic values. Because biological estimates are not available for all parameters, ranges were set reasonably large based on preliminary results and previous estimates (Smith et al., 2011a). The rate of infection (β) was allowed to vary between $10^{-6}\text{ TCI}{\text{D}}_{50}^{-1}{\text{d}}^{-1}$ and $10^{-1}\text{ TCI}{\text{D}}_{50}^{-1}{\text{d}}^{-1}$, and the rate of viral production (p) between $10^{-1}\text{ TCI}{\text{D}}_{50}\text{ cel}{\text{l}}^{-1}{\text{d}}^{-1}$ and $10^3\text{ TCI}{\text{D}}_{50}\text{ cel}{\text{l}}^{-1}{\text{d}}^{-1}$. Bounds for the viral clearance rate (c) were 1 d⁻¹ (t_1/2 = 16.7 h) and 10³ d⁻¹ (t_1/2 = 1 min). Previous estimates of the eclipse phase rate (k) for IAV infection in mice resulted in estimates that fell outside the biologically feasible range of 4–6 h (Smith et al., 2011a). To insure biological feasibility, the lower and upper bounds for the eclipse phase rate (k) were 4 d⁻¹ and 6 d⁻¹. Limits for the half-saturation constant (K_δ) were 10² − 10⁶ cells, and limits for the infected cell clearance parameter (δ_d) were 1 × 10⁶ − 4 × 10⁶ cells/d.

The initial number of target cells (T₀) was set to 10⁷ cells (Smith et al., 2011a, 2013). Because the initial viral inoculum rapidly infects cells and/or is cleared within 4 h pi, as indicated by the undetectable viral titers at this time point (Figure 1), the initial number of infected cells I₁(0) was set to 75 cells to reflect an initial dose of 75 TCID₅₀. This is an alteration from previous studies, including our own, that either estimate the initial amount of virus (V₀) or set its value to the true viral inoculum. Fixing V(0) = 75 TCID₅₀ or estimating its value did not improve the fit and could not be statistically justified (see, for example, Table S1). Further, fixing V(0) = 75 TCID₅₀ yielded an unreasonably high estimate for the rate of virus clearance (c) due to the attempt to fit the sharp decay between 0–4 h pi. Estimating I₁(0) could also not be justified and did not improve the model fit (e.g., as in Table S1 and Figure S1). The initial number of productively infected cells (I₂(0)) and the initial free virus (V₀) were set to 0.

2.8. Linear Regression

We used the function polyfit in MATLAB to perform linear regression of the log₁₀ values of viral titer during the growth phase (4 h, 1 d pi) and the two decay phases (2–6 d pi and 7–8 d pi).

3. Results

3.1. Phases of Viral Load Kinetics

Mice infected with 75 TCID₅₀ PR8 have viral load kinetics that can be separated into five distinct phases (Figure 1A). This is in contrast to the three phases that we previously defined (Smith et al., 2010). In the first phase, virus quickly infects cells and is undetectable within 4 h pi. In the second and third phases, virus increases exponentially and peaks after ~2 d pi. Following the peak, the viral decline can be separated into two phases. In the first decay phase (2–6 d pi), virus decays slowly at a relatively constant rate. In the second decay phase (7–8 d pi), virus declines rapidly (4–5 log₁₀ TCID₅₀). Sixty percent of mice had no detectable virus by 8 d pi. The remaining mice resolved the infection by 9 d pi.

These data reduced the heterogeneity observed in a previous data set from infection with the same virus (Smith et al., 2011a). We discovered that the majority of heterogeneity in the previous data set could be attributed to inconsistent infections and, thus, inocula that varied. We further reduced heterogeneity by normalizing the viral titer to the total lung volume, rather than using units of TCID₅₀/ml lung homogenate. As expected, some heterogeneity remains at 1 d pi and at 8 d pi. These time points correspond to when virus is rapidly increasing and decreasing, respectively.

3.2. Kinetic Model With Density Dependent Viral Clearance

We first fit the standard viral kinetic model, which is given by Equations (1)–(4) and assumes only one mechanism of constant clearance (δ(I₂) = δ_s) (Baccam et al., 2006), to the viral load data (see Supplementary Material). This model was unable to capture the entire time course of viral load dynamics, but was able to fit the data from infection initiation to 7 d pi (Figure S1). To more accurately model IAV kinetics and simultaneously recapitulate the two phases of viral decline, we modified the rate of infected cell clearance (δ(I₂)) so that the rate changes with respect to the density of the infected cell population. That is, δ(I₂) = δ_d/(K_δ + I₂) (Equation 5), where δ_d/K_δ is the maximum rate of clearance and K_δ is the number of productively infected cells where the rate is half of its maximum.

Fitting this new model to the viral load data illustrated that the model can accurately reproduce the data and simultaneously capture both phases of viral decline while excluding specific immune responses. The resulting dynamics are shown in Figure 1, the parameter values and 95% confidence intervals (CIs) are given in Table 1, and the parameter ensembles are shown in Figure 2 and Figure S3. For this model, the basic reproduction number (R₀) is given by

$\begin{array}{ll}{R}_0=\frac{\beta p{T}_0{K}_{\delta }}{c{\delta }_d}& \left(6\right)\end{array}$

Given the parameters in Table 1, R₀ = 8.8.

TABLE 1

Table 1. Parameter values and 95% confidence intervals obtained from fitting the density-dependent viral kinetic model (Equations 1–5) to viral titers from mice infected with 75 TCID₅₀ PR8.

FIGURE 2

Figure 2. Parameter ensembles and histograms. Parameter ensembles and histograms resulting from fitting the density-dependent viral kinetic model (Equations 1–5) to viral titers from mice infected with 75 TCID₅₀ PR8. Two main correlations are evident between the rates of virus production (p) and clearance (c) and between the rates of infection (β) and infected cell clearance (δ_d/K_δ). The axes limits reflect imposed bounds (except k ∈ [4, 6]). All parameters except the eclipse phase rate (k) are well bounded (i.e., the 95% CIs do not reach the imposed bounds). Additional ensemble plots (e.g., for R₀) are in Figure S3.

To understand how the addition of δ(I₂) = δ_d/(K_δ + I₂) influences the other parameters during the fitting scheme, we plotted the resulting histograms and 2D parameter projections (Figure 2). As expected, strong correlations exist between the rates of virus production (p) and virus clearance (c) and between the rate of infection (β) and the infected cell death rate (δ_d/K_δ). The other correlations visible in Figure 2 were a consequence of these two relations. Of note, δ_d was not strongly correlated with any of the other model parameters (Figure S3). In addition, the confidence interval was small, particularly compared to the other parameters. Estimates for the other parameters (β, p, c, and K_δ) with the exception of the eclipse phase rate (k) were well bounded such that the 95% CIs fell within the upper and lower bounds imposed in the estimation scheme. Similar to previous studies (Baccam et al., 2006; Smith et al., 2011a), the eclipse phase rate (k) was restricted to biologically realistic values and was not well defined on the given interval. In support, the ASA algorithm search patterns show a longer search time for k compared to the other parameters (Figure S2).

To further determine how the addition of δ(I₂) = δ_d/(K_δ + I₂) influences the sensitivity of the model solution to changes in parameter values, we performed a one-at-a-time sensitivity analysis (Figure 3). The infected cell clearance parameter (δ_d) is the most sensitive parameter and largely dictates the viral decay. Decreasing δ_d significantly delays viral clearance while increasing δ_d leads to rapid viral resolution (Figure 3). In accordance with previous results (Smith et al., 2010), all other parameters are less sensitive and collectively affect the exponential growth phase and peak.

FIGURE 3

Figure 3. Sensitivity of the density-dependent viral kinetic model. Solutions of the density-dependent viral kinetic model (Equations 1–5) for the best-fit parameters (black line, Table 1) and with the indicated parameter increased (red) or decreased (blue) 50% from the best-fit value.

As illustrated in Figure 4A, the rate of infected cell clearance is rapid when these cells are in small numbers. Given the parameters in Table 1, the maximum clearance rate is $\delta (I_2)=12.7{\text{d}}^{-1}$, which corresponds to half-life t_1/2 = 1.3 h. The rate begins to slow when $I_2>10^4\text{ cells}$ and is minimal when I₂ is at its maximum (8 × 10⁶ cells). When I₂ is maximal, $\delta (I_2)=0.21{\text{d}}^{-1}$ and t_1/2 = 78 h. In our previous work, we discovered that linear regression analysis could be used to accurately estimate the exponential growth rate, which was a combination of all model parameters, and that the slope of the viral decay could provide an estimate of δ(I₂) (Smith et al., 2010, 2011a). To evaluate how these relations correlate to parameters in the model with density dependence, we performed a linear regression on the data during the growth phase (4 h–1 d pi) and the two decay phases (2–6 d pi and 7–8 d pi) (Figure 4B). The slope of the growth phase is 4.7 log₁₀TCID₅₀/d (red line in Figure 4B). In accordance with the previous studies, this slope is a good approximation to the model until shortly before the peak. The model deviates slightly from this estimate and suggests that the virus growth rate briefly increases prior to the peak and that the decay phase begins prior to 2 d pi. This nonlinearity in the growth can be attributed to the decreasing infected cell clearance rate as the number of infected cells increases. These results are in contrast to the standard viral kinetic model, which suggests that the virus growth rate strictly decreases following exponential growth (Smith et al., 2010, 2011a). In the first phase of decay, the slope is −0.2 log₁₀ TCID₅₀/d (green line in Figure 4B), which corresponds to $\delta (I_2)=0.4{\text{d}}^{-1}$ (green diamond in Figure 4A). In the second phase of decay, the slope is −3.8 log₁₀ TCID₅₀/d (blue line in Figure 4B), which corresponds to $\delta (I_2)=8.7{\text{d}}^{-1}$ (blue dot in Figure 4A).

FIGURE 4

Figure 4. Density-dependent infected cell clearance rate and correlation to linear regression. (A) The infected cell clearance rate (δ(I₂), Equation 5) is plotted for different values of infected cells (I₂). The green diamond and the blue dot indicate the corresponding infected cell clearance rates during the slow and fast phases of virus clearance, respectively. These correspond to linear regression estimates in (B). (B) Linear regression fits to the viral load data (white squares) during the growth phase (4 h–1 d pi, red line), the first phase of virus decay (2–6 d pi, green line), or the second phase of virus decay (7–8 d pi, blue line). The dashed black line is the fit of the density-dependent viral kinetic model (Equations 1–5) to the viral load data.

4. Discussion

Mathematical models have been widely used to investigate IAV dynamics (reviewed in Smith and Ribeiro, 2010; Beauchemin and Handel, 2011; Smith and Perelson, 2011; Smith and McCullers, 2014; Boianelli et al., 2015). The viral kinetic model given by Equations (1)–(4) with δ(I₂) = δ_s (Baccam et al., 2006) has been the standard in the field for over 10 years. We previously used this model together with data from murine infection to gain insight into IAV virulence factors (Smith et al., 2011a) and into coinfection with bacterial pathogens (Smith et al., 2013; Smith and Smith, 2016; Smith, 2017). Although some predictions made using this model have been experimentally tested and deemed accurate (Ghoneim et al., 2013; Smith and McCullers, 2014; Warnking et al., 2015; Smith and Smith, 2016), the data here suggested that some dynamical features could not be accounted for and thus a new model was necessary. The model we introduced here includes density-dependent infected cell clearance and better captures the entire course of IAV infection dynamics, including the two-phase viral decay following the peak (Figure 1). Importantly, the model added only a single parameter (the half-saturation constant, K_δ) while significantly improving the model fit to viral loads from IAV infection without including additional equations detailing immune responses

By sampling with high frequency and controlling for experimental heterogeneity, we were able to obtain more accurate data (i.e., smaller standard deviations and better reproducibility) that highlighted several important dynamics, some of which were not previously observed. Our data showed that viral loads are maintained at a high level between 2 d and 7 d pi (Figure 1). Sustained viral loads have been observed in several studies (Jao et al., 1970; Douglas et al., 1975; Larson et al., 1976; Reuman et al., 1989; Bjornson et al., 1991; Toapanta and Ross, 2009; Smith et al., 2011a). In some data sets, the peak appears more pronounced and is often followed by the plateau phase or a second, lower peak (Jao et al., 1970; Douglas et al., 1975; Larson et al., 1976; Bjornson et al., 1991; Bender and Small, 1993; Hayden et al., 1998; Baccam et al., 2006). Our murine data do not indicate a second peak, although there is a subtle increase in viral loads at 5 d pi that may be biologically significant. Previous influenza modeling studies suggest that these dynamics required equations/terms for the innate and adaptive immune responses (Baccam et al., 2006; Pawelek et al., 2012; Cao et al., 2016). However, HIV modeling studies have used similar density-dependent terms to achieve a two-phase viral load decay (Holte et al., 2006; Burg et al., 2009). Importantly, the model here provides a means for capturing the changes in viral load decay without complicating the model or inferring information about specific immune mechanisms, which are not well understood. However, the change in clearance rate could reflect the change from innate to adaptive immunity. If this is the case, our estimates would suggest that the adaptive response is 25 times more effective than the innate response (−0.2 log₁₀ TCID₅₀/d between 2–6 d pi vs. −3.8 log₁₀ TCID₅₀/d between 7–8 d pi; Figure 4B).

It is well accepted that the rapid decline in virus during the second decay phase is due to the infiltration of CD8⁺ T cells (reviewed in McMichael et al., 1983; Kim et al., 2011; Grant et al., 2016). These cells typically enter the infection site between 5–6 d pi and peak between 8–9 d pi (e.g., as in Toapanta and Ross, 2009). The rapid rate of viral decline between 7–8 d pi suggests that these cells are highly effective. However, the initial infiltration begins at least 1–2 d before a change in the rate of virus decay is visible. Thus, there may be a nonlinearity to this response or it may reflect a simultaneous increase in infections and killing of infected cells or a change in effectiveness proportional to the density of infected cells or to the density of CD8⁺ T cells. A handling-time effect, which can represent the time required for immune cells to kill each infected cell and/or the time for immune cells to become activated (e.g., as in Pilyugin et al., 1997; Graw and Regoes, 2009; Smith et al., 2011b; Gadhamsetty et al., 2014; Li and Handel, 2014; Le et al., 2015), can slow the per capita rate of clearance. Spatial constraints (e.g., crowding effect), where the number of immune cells within an area is limited, may also play a role. In contrast, numerous clearance mechanisms (e.g., interferons, macrophages, neutrophils, natural killer (NK) cells) are thought to be important during early- and mid-infection, but their contribution to the viral load kinetics is unclear. Using a model to distinguish between these mechanisms is challenging given the close fit of simple kinetic models to viral load data (Figure 1 and Figure S1). Further, neither the data nor the models can discriminate whether the maintenance of high viral loads is due to a lack of clearance of infected cells (i.e., long infected cell lifespan/ineffective clearance) or to the balance of new infections and clearance (i.e., short infected cell lifespan/rapid clearance coupled with rapid virus infection/production). Thus, new experimental designs and more diverse data are necessary.

Viral titers remain the most frequently used data to calibrate models and assess infection dynamics. This is because collecting immunological data is more laborious and expensive. Thus, we seek models that are simple yet accurate and that can be used in the absence of immunological data. The standard viral kinetic model includes the minimal number of parameters and equations needed to recapitulate viral load dynamics. However, viral load data is typically insufficient to uniquely define all 6 parameters (Miao et al., 2011; Smith et al., 2011a). Fortunately, this has not limited our ability to make robust predictions about the underlying biology or to estimate accurate parameter values even when correlations are present (Gutenkunst et al., 2007; Smith et al., 2013; Smith and Smith, 2016). Here, the resulting parameter ensembles were well-bounded (i.e., the 95% CIs did not include the imposed bounds) and correlations were observed in two sets of parameters (Figure 2). The correlation between the rates of virus production (p) and virus clearance (c) indicates the balance of these processes. This is expected because viral loads measure the amount of virus present and slow virus production/clearance would be indistinguishable from fast production/clearance. Similarly, the rates of infection (β) and infected cell clearance (δ_d/K_δ) were correlated, which indicates a balance of cells becoming infected and being cleared. This is visible in Figure 4B, where the log-linear fit to the data in the growth phase (red line) deviates from the model solution (black dashed line).

Analyzing infection kinetics with mathematical models provides a means to quantify different infection processes. By modeling viral load data, we can make meaningful predictions about the time scales, magnitudes, and rates of different processes even if we cannot directly define specific mechanisms. Further, having a well-characterized model allows us to design new experiments and to perform in silico experiments that evaluate situations where data is challenging to obtain. Here, our data, model, and analyses suggest that the clearance rate of infected cells is variable and depends on their density such that clearance slows when infected cells are numerous and is fast when they are in low numbers. Determining what processes give rise to this density dependence remains an open question. Understanding how and why the rate changes should facilitate a deeper understanding of other viral infections and of immunological data, as it becomes available. Further establishing how the virus and host components work together and how they can be manipulated will undoubtedly aid the development of therapies that prevent or treat IAV infections.

Author Contributions

AMS conceived the idea, performed the experiments, developed the model, ran the simulations, performed the analysis, and wrote the manuscript. APS and DM performed the experiments. VB implemented the code and ran the simulations.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past co-authorship with one of the authors AMS.

Acknowledgments

This work was supported by NIH grants AI100946 and AI125324, and by ALSAC. A portion of this work was completed while all authors were at St. Jude Children's Research Hospital. We thank Alan Perelson and Laura Liao for their helpful comments.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2018.01554/full#supplementary-material

References

Baccam, P., Beauchemin, C., Macken, C., Hayden, F., and Perelson, A. (2006). Kinetics of influenza A virus infection in humans. J. Virol. 80, 7590–7599. doi: 10.1128/JVI.01623-05

PubMed Abstract | CrossRef Full Text | Google Scholar

Banerjee, S., Guedj, J., Ribeiro, R., Moses, M., and Perelson, A. (2016). Estimating biologically relevant parameters under uncertainty for experimental within-host murine West Nile virus infection. J. R. Soc. Interface 13:20160130. doi: 10.1098/rsif.2016.0130

PubMed Abstract | CrossRef Full Text | Google Scholar

Beauchemin, C., and Handel, A. (2011). A review of mathematical models of influenza A infections within a host or cell culture: lessons learned and challenges ahead. BMC Public Health 11(Suppl 1):S7. doi: 10.1186/1471-2458-11-S1-S7

PubMed Abstract | CrossRef Full Text | Google Scholar

Bender, B., and Small, P. Jr. (1993). Heterotypic immune mice lose protection against influenza virus infection with senescence. J. Infect. Dis. 168, 873–880. doi: 10.1093/infdis/168.4.873

PubMed Abstract | CrossRef Full Text | Google Scholar

Best, K., Guedj, J., Madelain, V., de Lamballerie, X., Lim, S., Osuna, C., et al. (2017). Zika plasma viral dynamics in nonhuman primates provides insights into early infection and antiviral strategies. Proc. Natl. Acad. Sci. U.S.A. 114, 8847–8852. doi: 10.1073/pnas.1704011114

PubMed Abstract | CrossRef Full Text | Google Scholar

Bjornson, A., Mellencamp, M., and Schiff, G. (1991). Complement is activated in the upper respiratory tract during influenza virus infection. Am. Rev. Respir. Dis. 143(5 Pt 1):1062–1066.

PubMed Abstract | Google Scholar

Boianelli, A., Nguyen, V., Ebensen, T., Schulze, K., Wilk, E., Sharma, N., et al. (2015). Modeling influenza virus infection: a roadmap for influenza research. Viruses 7, 5274–5304. doi: 10.3390/v7102875

PubMed Abstract | CrossRef Full Text | Google Scholar

Burg, D., Rong, L., Neumann, A., and Dahari, H. (2009). Mathematical modeling of viral kinetics under immune control during primary HIV-1 infection. J. Theor. Biol. 259, 751–759. doi: 10.1016/j.jtbi.2009.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, P., Wang, Z., Yan, A., McVernon, J., Xu, J., Heffernan, J., et al. (2016). On the role of CD8+ T cells in determining recovery time from influenza virus infection. Front. Immunol. 7:611. doi: 10.3389/fimmu.2016.00611

PubMed Abstract | CrossRef Full Text | Google Scholar

Douglas, R. Jr., Betts, R., Simons, R., Hogan, P., and Roth, F. (1975). Evaluation of a topical interferon inducer in experimental influenza infection in volunteers. Antimicrob. Agents Ch. 8, 684–687. doi: 10.1128/AAC.8.6.684

PubMed Abstract | CrossRef Full Text | Google Scholar

Gadhamsetty, S., Marée, A., Beltman, J., and de Boer, R. (2014). A general functional response of cytotoxic T lymphocyte-mediated killing of target cells. Biophys. J. 106, 1780–1791. doi: 10.1016/j.bpj.2014.01.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghoneim, H., Thomas, P., and McCullers, J. (2013). Depletion of alveolar macrophages during influenza infection facilitates bacterial superinfections. J. Immunol. 191, 1250–1259. doi: 10.4049/jimmunol.1300014

PubMed Abstract | CrossRef Full Text | Google Scholar

Grant, E., Quiñones-Parra, S., Clemens, E., and Kedzierska, K. (2016). Human influenza viruses and CD8+ T cell responses. Curr. Opin. Virol. 16, 132–142. doi: 10.1016/j.coviro.2016.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Graw, F., and Regoes, R. (2009). Investigating CTL mediated killing with a 3D cellular automaton. PLoS Comput. Biol. 5:e1000466. doi: 10.1371/journal.pcbi.1000466

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutenkunst, R., Waterfall, J., Casey, F., Brown, K., Myers, C., and Sethna, J. (2007). Universally sloppy parameter sensitivities in systems biology models. PLoS Comput. Biol. 3:e189. doi: 10.1371/journal.pcbi.0030189

PubMed Abstract | CrossRef Full Text | Google Scholar

Halloran, S., Wexler, A., and Ristenpart, W. (2012). A comprehensive breath plume model for disease transmission via expiratory aerosols. PLoS ONE 7:e37088. doi: 10.1371/journal.pone.0037088

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayden, F., Fritz, R., Lobo, M., Alvord, W., Strober, W., and Straus, S. (1998). Local and systemic cytokine responses during experimental human influenza A virus infection. Relation to symptom formation and host defense. J. Clin. Investig. 101:643.

PubMed Abstract | Google Scholar

Holder, B., Liao, L., Simon, P., Boivin, G., and Beauchemin, C. (2011a). Design considerations in building in silico equivalents of common experimental influenza virus assays. Autoimmunity 44, 282–293. doi: 10.3109/08916934.2011.523267

PubMed Abstract | CrossRef Full Text | Google Scholar

Holder, B., Simon, P., Liao, L., Abed, Y., Bouhy, X., Beauchemin, C., et al. (2011b). Assessing the in vitro fitness of an oseltamivir-resistant seasonal A/H1N1 influenza strain using a mathematical model. PLoS ONE 6:e14767. doi: 10.1371/journal.pone.0014767

PubMed Abstract | CrossRef Full Text | Google Scholar

Holte, S., Melvin, A., Mullins, I., Tobin, N., and Frenkel, L. (2006). Density-dependent decay in HIV-1 dynamics. JAIDS 41, 266–276. doi: 10.1097/01.qai.0000199233.69457.e4

PubMed Abstract | CrossRef Full Text | Google Scholar

Jao, R., Wheelock, E., and Jackson, G. (1970). Production of interferon in volunteers infected with Asian influenza. J. Infect. Dis. 121, 419–426. doi: 10.1093/infdis/121.4.419

PubMed Abstract | CrossRef Full Text | Google Scholar

Kakizoe, Y., Nakaoka, S., Beauchemin, C., Morita, S., Mori, H., Igarashi, T., et al. (2015). A method to determine the duration of the eclipse phase for in vitro infection with a highly pathogenic SHIV strain. Sci. Rep. 5:10371. doi: 10.1038/srep10371

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, T., Sun, J., and Braciale, T. (2011). T cell responses during influenza infection: getting and keeping control. Trends Immunol. 32, 225–231. doi: 10.1016/j.it.2011.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Larson, E., Dominik, J., Rowberg, A., and Higbee, G. (1976). Influenza virus population dynamics in the respiratory tract of experimentally infected mice. Infect. Immun. 13, 438–447.

PubMed Abstract | Google Scholar

Le, D., Miller, J., and Ganusov, V. (2015). Mathematical modeling provides kinetic details of the human immune response to vaccination. Front. Cell. Infect. Microbiol. 4:177. doi: 10.3389/fcimb.2014.00177

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., and Handel, A. (2014). Modeling inoculum dose dependent patterns of acute virus infections. J. Theor. Biol. 347, 63–73. doi: 10.1016/j.jtbi.2014.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

McMichael, A., Gotch, F., Noble, G., and Beare, P. (1983). Cytotoxic T-cell immunity to influenza. New Engl. J. Med. 309, 13–17. doi: 10.1056/NEJM198307073090103

PubMed Abstract | CrossRef Full Text | Google Scholar

Medina, R., and García-Sastre, A. (2011). Influenza A viruses: new research developments. Nat. Rev. Microbiol. 9, 590–603. doi: 10.1038/nrmicro2613

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, H., Hollenbaugh, J., Zand, M., Holden-Wiltse, J., Mosmann, T., Perelson, A., et al. (2010). Quantifying the early immune response and adaptive immune response kinetics in mice infected by influenza A virus. J. Virol. 84, 6687–6698. doi: 10.1128/JVI.00266-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, H., Xia, X., Perelson, A., and Wu, H. (2011). On identifiability of nonlinear ODE models and applications in viral dynamics. SIAM Rev. 53, 3–39. doi: 10.1137/090757009

PubMed Abstract | CrossRef Full Text | Google Scholar

Neumann, A., Lam, N., Dahari, H., Gretch, D., Wiley, T., Layden, T., et al. (1998). Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-α therapy. Science 282, 103–107. doi: 10.1126/science.282.5386.103

PubMed Abstract | CrossRef Full Text | Google Scholar

Palmer, J., Dobrovolny, H., and Beauchemin, C. (2017). The in vivo efficacy of neuraminidase inhibitors cannot be determined from the decay rates of influenza viral titers observed in treated patients. Sci. Rep. 7:40210. doi: 10.1038/srep40210

PubMed Abstract | CrossRef Full Text | Google Scholar

Pawelek, K., Huynh, G., Quinlivan, M., Cullinane, A., Rong, L., and Perelson, A. (2012). Modeling within-host dynamics of influenza virus infection including immune responses. PLoS Comput. Biol. 8:e1002588. doi: 10.1371/journal.pcbi.1002588

PubMed Abstract | CrossRef Full Text | Google Scholar

Perelson, A., Neumann, A., Markowitz, M., Leonard, J., and Ho, D. (1996). HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science 271, 1582–1586. doi: 10.1126/science.271.5255.1582

PubMed Abstract | CrossRef Full Text | Google Scholar

Pilyugin, S., Mittler, J., and Antia, R. (1997). Modeling T-cell proliferation: an investigation of the consequences of the Hayflick limit. J. Theor. Biol. 186, 117–129. doi: 10.1006/jtbi.1996.0319

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinky, L., and Dobrovolny, H. (2016). Coinfections of the respiratory tract: viral competition for resources. PLoS ONE 11:e0155589. doi: 10.1371/journal.pone.0155589

PubMed Abstract | CrossRef Full Text | Google Scholar

Reed, L., and Muench, H. (1938). A simple method of estimating fifty percent endpoints. Am. J. Epidemiol. 27, 493–497. doi: 10.1093/oxfordjournals.aje.a118408

CrossRef Full Text | Google Scholar

Reuman, P., Bernstein, D., Keefer, M., Young, E., Sherwood, J., Schiff, G. M., et al. (1989). Efficacy and safety of low dosage amantadine hydrochloride as prophylaxis for influenza A. Antivir. Res. 11, 27–40. doi: 10.1016/0166-3542(89)90018-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Simonsen, L., Fukuda, K., Schonberger, L., and Cox, N. (2000). The impact of influenza epidemics on hospitalizations. J. Infect. Dis. 181, 831–837. doi: 10.1086/315320

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A. (2017). Quantifying the therapeutic requirements and potential for combination therapy to prevent bacterial coinfection during influenza. J. Pharmacokinet. Pharm. 48:81. doi: 10.1007/s10928-016-9494-9

CrossRef Full Text | Google Scholar

Smith, A., Adler, F., McAuley, J., Gutenkunst, R., Ribeiro, R., McCullers, J., et al. (2011a). Effect of 1918 PB1-F2 expression on influenza A virus infection kinetics. PLoS Comput. Biol. 7:e1001081. doi: 10.1371/journal.pcbi.1001081

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A., Adler, F., and Perelson, A. (2010). An accurate two-phase approximate solution to an acute viral infection model. J. Math. Biol. 60, 711–726. doi: 10.1007/s00285-009-0281-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A., Adler, F., Ribeiro, R., Gutenkunst, R., McAuley, J., McCullers, J., et al. (2013). Kinetics of coinfection with influenza A virus and Streptococcus pneumoniae. PLoS Pathog. 9:e1003238. doi: 10.1371/journal.ppat.1003238

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A., and McCullers, J. (2014). “Secondary bacterial infections in influenza virus infection pathogenesis,” in Influenza Pathogenesis and Control, Vol. 1. Current Topics in Microbiology and Immunology, vol. 385. (Cham: Springer), 327–356.

Google Scholar

Smith, A. M., and McCullers, J. A. (2014). Secondary Bacterial Infections in Influenza Virus Infection Pathogenesis. In: Compans R., Oldstone M. (eds) Influenza Pathogenesis and Control - Volume I. Current Topics in Microbiology and Immunology, vol 385. Springer, Cham

PubMed Abstract | Google Scholar

Smith, A., McCullers, J., and Adler, F. (2011b). Mathematical model of a three-stage innate immune response to a pneumococcal lung infection. J. Theor. Biol. 276, 106–116. doi: 10.1016/j.jtbi.2011.01.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A., and Perelson, A. (2011). Influenza A virus infection kinetics: quantitative data and models. WIREs Sys. Biol. Med. 3, 429–445. doi: 10.1002/wsbm.129

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A., and Ribeiro, R. (2010). Modeling the viral dynamics of influenza A virus infection. Crit. Rev. Immunol. 30, 291–298. doi: 10.1615/CritRevImmunol.v30.i3.60

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A., and Smith, A. (2016). A critical, nonlinear threshold dictates bacterial invasion and initial kinetics during influenza. Sci. Rep. 6:38703. doi: 10.1038/srep38703

PubMed Abstract | CrossRef Full Text | Google Scholar

Taubenberger, J., and Morens, D. (2008). The pathology of influenza virus infections. Annu. Rev. Pathmechdis. Mech. Dis. 3, 499–522. doi: 10.1146/annurev.pathmechdis.3.121806.154316

PubMed Abstract | CrossRef Full Text | Google Scholar

Thompson, W., Shay, D., Weintraub, E., Brammer, L., Bridges, Cox, C. B., et al. (2004). Influenza-associated hospitalizations in the United States. J. Am. Med. Assoc. 292, 1333–1340. doi: 10.1001/jama.292.11.1333

PubMed Abstract | CrossRef Full Text | Google Scholar

Toapanta, F., and Ross, T. (2009). Impaired immune responses in the lungs of aged mice following influenza infection. Respir. Res. 10:112. doi: 10.1186/1465-9921-10-112

PubMed Abstract | CrossRef Full Text | Google Scholar

Warnking, K., Klemm, C., Löffler, B., Niemann, S., Krüchten, A., Peters, G., et al. (2015). Super-infection with Staphylococcus aureus inhibits influenza virus-induced type I IFN signalling through impaired STAT1-STAT2 dimerization. Cell. Microbiol. 17, 303–317. doi: 10.1111/cmi.12375

PubMed Abstract | CrossRef Full Text | Google Scholar