The Nonlinear Relations that Predict Influenza Viral Dynamics, CD8

Influenza A viruses cause a significant amount of morbidity and mortality. Understanding how the infection is controlled by host immune responses and how different factors influence severity are critical to combat the infection. During infection, viral loads increase exponentially, peak, then decline until resolution. The viral decline is often biphasic, which we previously determined is a consequence of density-dependent infected cell clearance. The second, rapid clearance phase corresponds with the infiltration of CD8 + T cells, but how the rate changes with infected cell density and T cell density is unclear. Further, the kinetics of virus, infected cells, and CD8 + T cells all contribute to disease severity but do not seem to be directly correlated. Thus, we investigated the relations between viral loads, infected cells, CD8 + T cells, lung pathology, and disease severity/symptoms by infecting mice with influenza A/PR8, simultaneously measuring virus and CD8 + T cells, and developing and calibrating a kinetic model. The model predicted that infection resolution is sensitive to CD8 + T cell expansion, that there is a critical T cell magnitude below which the infection is significantly prolonged, and that the efficiency of T cell-mediated clearance is dependent on infected cell density. To further examine the latter finding and validate the model’s predicted dynamics, we quantified infected cell kinetics using lung histomorphometry. These data showed that the area of lung infected matches the predicted cumulative infected cell dynamics, and that the area of resolved infection parallels the relative CD8 + T cell magnitude. Our analysis further revealed a nonlinear relationship between disease severity (i.e., weight loss) and the percent of the lung damaged. Establishing the predictive capabilities of the model and the critical connections that map the kinetics of virus, infected cells, CD8 + T cells, lung pathology, and disease severity during influenza virus infection aids our ability to forecast the course of infection, disease progression, and potential complications, thereby providing insight for clinical decisions.


Author Summary
Influenza A viruses infect millions of people each year. An understanding of how virus growth and host responses impact disease progression is critical to identify disease-specific markers that help predict hospitalization, complications, and therapeutic efficacy. To establish these relations, we developed and validated a mathematical model that accurately forecasts the kinetics of virus, infected cells, and CD8 + T cells. We discovered that the rate of infected cell removal by CD8 + T cells increases as infected cells decline, that there is a critical number of CD8 + T cells below which recovery is prolonged, and that recovery time depends on the number of CD8 + T cells rather than their efficiency. Further, examining lung pathology showed that the area of the lung infected and the area of the lung resolved parallel the model's predicted cumulative infected cell and CD8

Introduction
Over 15 million respiratory infections and 200,000 hospitalizations result from influenza A viruses (IAVs) each year [1][2][3][4]. The incidence and severity of IAV infections increases when new strains emerge and/or when there is a lack of prior immunity. A robust immune response is crucial for resolving viral infections, but immune-mediated pathology can exacerbate disease [5][6][7][8][9]. High viral loads also play a role in disease progression [10], but these do not always correlate with the strength of the host response or with disease severity [11][12][13][14]. An understanding of how viral loads, host immune responses, and disease progression are related is critical to identify disease-specific markers that may help predict hospitalization or other complications.
A better understanding of infected cell clearance may also yield insight into the damage induced to the lung during IAV infection. In general, widespread alveolar disease is observed in patients who succumb to the infection [36]. Further, in hospitalized patients who died as a result from infection with the 2009 H1N1 influenza virus had large numbers of CD8 + T cells present in their lung tissue [7]. Large pulmonary populations of CD8 + T cells may result in 'bystander damage' to uninfected epithelial cells [37] in addition to targeting IAV-infected cells. The accumulation of damage to the epithelium during IAV infection, either from virally-induced cell lysis or immune-mediated effects, is relatively understudied. This is in part due to the difficulty in measuring the dynamics of infected cells and in establishing how damage correlates to different 3/35 host responses. Recent technological advances, including the use of reporter viruses [38][39][40][41] and lung histomorphometry [12,[42][43][44], have provided an opportunity to acquire these types of measurements. However, even with these techniques, quantitative data over the course of infection is not currently available. Having data such as these should help reconcile potential nonlinearities in infected cell clearance and provide insight into the accumulated lung damage, which we hypothesize is a marker of disease severity.
In general, measures of disease severity do not seem to be directly correlated to viral loads or to specific immunological components. In humans infected with IAV, viral loads typically increase prior to the onset of systemic symptoms, which peak around 2-3 d post-infection (pi) [18,45,46].
Respiratory symptoms often last longer and can remain following viral clearance [45]. In the murine model, weight loss is used as an indicator of disease progression and severity, where greater weight loss corresponds to more severe disease [47][48][49]. Animal weights typically drop slowly in the first several days during an IAV infection and more rapidly as the infection begins to resolve [19,50]. This is unlike viral load dynamics in these animals, which increase rapidly during the first 0-3 d pi then remain relatively constant prior to resolution [15]. Because weight loss often occurs following resolution of the infection, immune-mediated pathology is thought to play a role [5,49,[51][52][53]. Host and pathogen factors, such as age, viral proteins, and inoculum size, can also influence disease progression [13,14,54,55]. While the causes of IAV-associated disease and mortality remain elusive, this gap in knowledge impairs our ability to effectively predict, understand, and treat the disease.
To gain deeper insight into the dynamics of viral resolution and investigate the connection between viral loads and disease severity, we simultaneously measured viral loads and CD8 + T cells daily from the lungs of BALB/cJ mice infected with influenza A/Puerto Rico/8/34 (H1N1) (PR8). We then developed a model that describes their kinetics to explore the mechanisms and dynamics of CD8 + T cell influx and their efficiency in removing virus-infected cells. The model verified our previous results that infected cells are cleared in a density-dependent manner. We further determined that infection duration is dependent on the magnitude of CD8 + T cells rather than their efficacy. Exploring these findings through quantitative whole lung histomorphometry corroborated the predicted infected cell dynamics and the direct relation between infected cell clearance and CD8 + T cell expansion. In addition, these data revealed a nonlinear connection between disease severity, as measured by weight loss, and the amount of the lung impacted by the infection. The data, model, and analyses provide a robust quantification of the density-dependent nature of CD8 + T cell-mediated clearance, and the critical connections between these cells and 4/35 the dynamics of viral loads, infected cells, lung pathology, and disease severity.

Results
Virus and CD8 + T cell kinetics In animals infected with 75 TCID 50 PR8, virus rapidly infects cells or is cleared within 4 h pi (Fig 1). Virus then increases exponentially and peaks after ∼2 d pi. Following the peak, virus enters a biphasic decline. In the first phase (2-6 d pi), virus decays slowly and at a relatively constant rate (0.2 log 10 TCID 50 /d) [15]. In the second phase (7-8 d pi), virus is cleared rapidly with a loss of 4 − 5 log 10 TCID 50 in 1-2 d (average of −3.8 log 10 TCID 50 /d) [15]. CD8 + T cells remain at their baseline level from 0-2 d pi before they infiltrate the lung tissue and increase slightly at 2-3 d pi. The population briefly contracts (3-5 d pi) before expanding rapidly ( Viral kinetic model with CD8 + T cell-mediated clearance We previously described the viral load kinetics and biphasic decline using the density-dependent (DD) model in Eq (1)-(4), which assumes that the rate of infected cell clearance increases as the density of infected cells decreases (i.e., δ d (I 2 ) = δ d /(K δ + I 2 )) [15]. Because the second viral decay phase is thought to be due to the clearance of infected cells by CD8 + T cells, we developed a model to describe the dynamics of these cells and their efficiency in resolving the infection (Eq (5)-(10); Fig 1A). The model includes equations for effector (E, denoted CD8 E ) and memory (E M , denoted CD8 M ) CD8 + T cells, and two mechanisms of infected cell clearance.
During the first viral decay phase (2-6 d pi), the rate of infected cell clearance by non-specific mechanisms is relatively constant (δ). During the rapid, second viral decay phase (7-8 d pi), CD8 E -mediated infected cell clearance occurs at a rate that increases as the density of infected cells decreases (δ E (I 2 , E) = δ E E/(K δ E + I 2 )). Excluding this density dependence resulted in a significant and premature decline in viral loads, which disagreed with the experimental data.
Although memory CD8 + T cells are not the focus here, it was necessary to include the CD8 M population because CD8 + T cells are at a significantly higher level at 10-12 d pi than at 0 d 5/35 Parameters and 95% confidence intervals obtained from fitting the CD8 + T cell model (Eq (5)-(10)) to viral titers and CD8 + T cells from mice infected with 75 TCID 50 PR8. CD8 E and CD8 M denote effector (E) and memory (E M ) CD8 + T cells, respectively. The total number of CD8 + T cells isÊ = E + E M +Ê 0 and is denoted by CD8. pi ( Fig 1B). The model includes terms for CD8 E infiltration (ξI 2 /(K E + E)), which accounts for the initial increase at 2-3 d pi, and for CD8 E expansion (ηEI 2 (t − τ E )), which accounts for the larger increase between 5-8 d pi. To capture the contraction of CD8 + T cells between these times (3-5 d pi), it was necessary to assume that CD8 E infiltration is reduced by their own population (i.e., ξ(E) = ξ/(K E + E)). In both terms, the increase is proportional to the number of infected cells. Fitting the model simultaneously to viral loads and CD8 + T cells from the lungs of infected animals illustrated the accuracy of the model (Fig 1B). The resulting parameter values, ensembles, histograms, and 95% confidence intervals (CIs) are given in Table 1 and Figs 2 and S1-S2.
Plotting the model ensembles indicated that the behavior of the virus-specific parameters (i.e., the rates of virus infectivity (β), virus production (p), virus clearance (c), and eclipse phase transition (k)) in the CD8 + T cell model were similar to the results obtained fitting the DD model to the viral load data [15] (Fig 2A). In addition, they revealed a correlation between the 6/35 two infected cell clearance parameters (δ and δ E ; Fig 2B), which represent the efficacy of the non-specific and CD8 + T cell immune responses, respectively. Performing a sensitivity analysis showed that the viral load dynamics do not change substantially when these parameters are increased or decreased. However, decreasing the rate of non-specific infected cell clearance (i.e., lower δ) resulted in a significant increase in the number of CD8 E due to the small increase in the number of infected cells (Fig S3). Even with a larger CD8 E population, recovery was delayed by only ∼ 0.1 d. Given the correlation between δ and δ E (Fig 2B), a more efficient CD8 E response (i.e., higher δ E ) may be able to overcome this short delay in resolution. The lack of sensitivity to changes in the infected cell clearance parameters is in contrast to the DD model, where the viral dynamics were most sensitive to perturbations in δ d (Fig S5) [15], which encompasses multiple processes. With CD8 + T cells explicitly included in the model, the infection duration was most sensitive to changes in the rate of CD8 E expansion (η) (Figs S3 and S5; discussed in more detail below).
Examining the parameter ensembles and sensitivity analysis also yielded insight into how other model parameters affect the CD8 + T cell response. The rates of CD8 E expansion (η) and clearance (d E ) were slightly correlated, indicating a balance between these two processes ( Fig 2C).
This correlation and the sensitivity of η produce model dynamics that were also sensitive to changes in the CD8 E clearance rate (d E ) ( Fig S4). As expected, the rates of CD8 M generation (ζ) and CD8 E clearance (d E ) were correlated ( Fig 2C). It has been estimated that approximately 5-10% of effector CD8 + T cells survive to become a long-lasting memory population [56]. Despite the inability to distinguish between CD8 E and CD8 M in the data, the model predicts that 17% of CD8 E transitioned to the memory class by 15 d pi. Additional insight into the regulation of the CD8 + T cell response, results from the model fitting, and a comparison of the DD model and the CD8 + T cell model are in the Supporting Information.

Density-dependent infected cell clearance
To gain further insight into the nonlinear dynamics of CD8 ) as a function of infected cells (I 2 ) and CD8 E (E) (Fig 3). This confirmed that there is minimal contribution from CD8 E -mediated clearance to viral load kinetics or infected cell kinetics prior to 7 d pi ( At the initiation of the second decay phase (7 d pi), the clearance rate is δ E (I 2 , E) = 3.5/d ( Fig 3A-B, marker c). As the infected cell density declines towards the half-saturation constant 7/35 (K δ E = 4.3 × 10 2 cells), the clearance rate increases rapidly to a maximum of 4830/d (Fig 3A- To explore how recovery time is altered by varying CD8 E levels, we examined the resulting dynamics from increasing or decreasing the rate of CD8 E expansion (η). When η was increased by 50%, the CD8 E population increased by a substantial 670% ( Fig S3)

Lung histomorphometry validates infected cell dynamics
To investigate the dynamics of infected cells, we quantified these cells and the infection progression and resolution within the lung using serial whole lung histomorphometry ( Fig 4A). Antigenpositive areas of the lung ("active" lesions) were first detectable at 2 d pi ( Fig 4A-B), which coincides with the peak in viral loads. The infected area continued to increase until 6 d pi, whereas viral loads remained high until 7 d pi (Fig 1B). At this time, resolution of the infection began and ∼28.7% of the infected area was cleared between 6-7 d pi ( Fig S6A). Few to no infected cells were present by 8 d pi ( Fig 4A). Correspondingly, virus was undetectable in most animals by 8 d pi (Fig 1B). Because the percent active lesion is a reflection of the influenza-positive cells, we examined whether the CD8 + T cell model accurately predicted these dynamics. In the model, the accumulated infection is defined by the cumulative area under the curve (CAUC) of the productively infected cells (I 2 ). Plotting the percent active lesion against the CAUC of I 2 showed that the model accurately predicts the cumulative infected cell dynamics and, thus, the 8/35 infection progression within the lung (Fig 4B).
Antigen-negative, previously-infected or damaged areas of the lung ("inactive" lesions) are evident beginning at 5 d pi (Fig 4A,C). This resolution of the infection continued from 5-8 d pi, causing a 14.6%/d increase in the area of inactive lesions (Fig S6B). Following this, healing of the lung is apparent as the inactive lesioned areas decline (-14.5%/d from 9-10 d pi; Figs 4C and S6B). These dynamics parallel the CD8 + T cell dynamics (Fig 1B). Fitting a line to the CD8 + T cell data from 5-8 d pi indicated that the influx rate is 4.7 × 10 5 cells/d ( Fig S6C).
Thus, every 100,000 CD8 + T cells clear ∼ 3.1% of the lung. During CD8 + T cell contraction phase (8-10 d pi; Fig 1B), a similar linear regression analysis suggested that these cells decline at a rate of ∼ 3.3 × 10 5 CD8/d (Fig S6C). Similar to the relation discussed above, the dynamics of the damaged areas of the lung corresponded precisely to the dynamics of the percent maximum CD8 E (i.e., E/E max ) in the model (Fig 4C).

Weight loss predicts area of lung affected
To monitor disease progression, weight loss was measured daily throughout the course of infection ( Fig 4D). During the first 5 d pi, animals gradually lost ∼4% of their initial weight. This was followed by a sharper drop (8%) at 6 d pi. Animal weights increased slightly at 7 d pi (∼6%) before reaching peak weight loss (10-14%) at 8 d pi. Following virus resolution, the animals' weights began to restore as the inactive lesions resolved (9-10 d pi; Fig 4D). We hypothesized that these weight loss dynamics may reflect the area of the lung impacted by the infection. Indeed, plotting weight loss together with the percent total (active and inactive) lesioned area of the lung revealed the similarity in their dynamics (Fig 4D-E). To further quantify their relationship, we plotted the percent weight loss against the percent total lesioned area and observed a nonlinear relation. Thus, we fit the saturating function L = l max W n /(W n + K n w ) to these data, where L is the percent total lesioned area, W is percent weight loss, l max is the maximum rate of the interaction, K w is the half-saturation constant, and n is the Hill coefficient. This function provided a close fit to the data (R 2 = 0.9; Fig 4E) with best-fit parameters l max =39.7% total lesioned area, K w =2.6% weight loss, and n = 5.2.

Discussion
Influenza A virus infections pose a significant threat to human health, and it is crucial to understand how the virus spreads within the respiratory tract, how specific immune responses 9/35 contribute to infection control, and how these relate to disease progression. Although it has been difficult to directly relate these features, we circumvented the challenge by pairing comprehensive, experimental data with robust mathematical models and analyses. Our iterative model-driven experiment approach [57,58] revealed important dynamic relations between virus, infected cells, CD8 + T cells, lung damage, and weight loss (Fig 5). Identifying these nonlinear connections allows for more accurate interpretations of viral infection data and significant improvement in our ability to predict disease severity, the likelihood of complications, and therapeutic efficacy.
Our histomorphometric data provided the first quantification of the spread of influenza virus infection within the lung. These data allowed us to validate our model's predicted infected cell dynamics and, thus, confirm that their density impacts the rate at which they are cleared by effector CD8 + T cells (CD8 E ) (Fig 4). We first detected this density-dependence in a model that excluded specific immune responses (i.e., Eq (1)- (4)) [15]. That model was a mathematically elegant way to capture the nonlinearity in the viral load decline, but it could not distinguish between different mechanisms of viral clearance. Here, we verified that the second viral clearance phase reflects CD8 + T cell-mediated clearance, determined that their efficiency increases as more infected cells are removed, and identified the critical level needed for a timely recovery (Fig 3).
Several factors may contribute to the density-dependent change in the rate of CD8 E -mediated clearance. One possibility is that the slowed rate of clearance at high infected cell densities is due to a "handling time" effect, which represents the time required for an immune cell to remove an infected cell (e.g., as in [15,24,29,[59][60][61][62]). When CD8 E interact with infected cells, a complex is formed for ∼20-40 min [32,[63][64][65][66]. Because CD8 E cannot interact with other infected cells during this time, the global rate of infected cell clearance would be lowest when infected cells outnumber CD8 E . In addition, contact by a single CD8 E may be insufficient to remove an infected cell [28]. Infected cell clearance is more frequently observed after interactions with multiple CD8 E , with an average of 3.9 contact events either serially or simultaneously [28]. Thus, the high density of infected cells early in the infection reduces the probability that a single virus-infected cell would be targeted a sufficient number times to induce cell death. However, as CD8 E accumulate and the density of infected cells decreases (Fig 3A), the probability of simultaneous interactions will increase. This should reduce the handling time required to remove of an infected cell and, thus, result in a higher efficiency. While this is a likely explanation, it is possible that spatial limitations also contribute, such that a high infected cell density may hinder CD8 + T cells from reaching infected cells on the interior of the infected tissue. Crowding of immune cells at the periphery of infected lesions has been observed in other infections (e.g.,
The new knowledge about how the infection spreads throughout the lung also uncovered a novel connection between the percent weight loss of an infected animal (i.e., disease severity) and the percent total area of the lung impacted by the virus infection (Fig 4D-E). This discovery is significant because it suggests that disease severity is linked to the amount of lung involvement, which we showed is directly connected to infected cell and CD8 E kinetics (Fig 4B-C) and indirectly connected to viral load kinetics (Fig 5). Other studies utilizing histomorphometry data, although not quantitatively, provide support for the relationship between weight loss and lung involvement during IAV infection [12]. For example, animals treated with antivirals in various conditions (single or combination therapy and in immunocompetent or immunosuppressed hosts) demonstrated reduced weight loss that corresponded to reduced infected areas of lung [12,42].
Examining the histomorphometric kinetics and the connection between weight loss and lung pathology in various infection settings should improve our predictive capabilities. This may be particularly helpful in understanding the exacerbated morbidity and mortality in elderly individuals. In the animal model, elderly mice typically have lower viral loads compared to adult mice yet experience more weight loss, symptoms, and/or mortality [13,14].
The nonlinear relations derived from the histomorphometric analysis further provide a way to interpret and utilize weight loss data, which is the most readily available type of data and commonly used as a measure of disease severity in infection models (Figs 4-5). In our data and others', there is a spike in weight loss, symptom score, and/or inflammation at ∼6 d pi [50,54,70,71] (Fig 4D). Our data suggests that this is due to CD8 E infiltration and activity while the infection continues to spread (Fig 4). Utilizing this observation and the function describing the correlation between the percent lung lesion and weight loss (L(W ); Fig 4E), the kinetics of the active and inactive lesioned areas can be directly estimated from weight loss data ( Fig 5). Taking this one step further, the dynamics of the CAUC of the infected cells and the relative CD8 E can be obtained. Finally, the CD8 + T cell model (Eq (5)-(10)) can be employed to estimate the viral load and CD8 + T cell dynamics (Fig 5). These connections can be reversed to estimate disease severity (i.e., weight loss) from viral load data.
Weight loss increases following infection resolution (8-9 d pi), which is likely attributed to CD8 + T cell-mediated pathology (Fig 4C). This is supported by other studies that suggest a large number of CD8 + T cells poses a risk of acute lung tissue injury [5][6][7][8][9]. According to our model predictions, excessive CD8 + T cell numbers may augment disease progression yet do not 11/35 improve recovery time (Fig S3). Instead, an earlier onset of CD8 E proliferation (i.e., smaller τ E ) would be required to significantly shorten the infection (Fig S4). This aligns with evidence that hosts with adaptive immune responses primed by vaccine or prior infection recover more rapidly [72,73]. While higher CD8 + T cell numbers have little impact on viral kinetics, the model agrees with clinical and experimental studies from a wide range of host species that impaired CD8 + T cell responses can prolong an IAV infection [17,22,[74][75][76][77]. That is, virus can persist for up to several weeks if CD8 + T cell-mediated clearance is unsuccessful (Fig 3; [74, 78, 79]). The bifurcation in recovery time revealed by the model suggests that this may occur when the total number of CD8 + T cells are <39.2% of their maximum (Fig 3D-E)). This number is expected to vary depending on parameters like the rate of virus replication and the infected cell lifespan, which has been noted in another modeling study that detailed similar bifurcation behavior [80].
Although some previously published models also predict delayed resolution with depleted CD8 + T cell responses [17,33,34,81], this bifurcation has not been observed and their estimated delays in recovery do not amount to the long-lasting infections in immunodeficient patients [74,78,79].
In addition to illuminating the connections between virus spread, virus clearance, the associated pathology, and the severity of disease, the histomorphometric data validated the model's infected cell dynamics (Fig 4B). The dynamics of susceptible and infected cells throughout the infection and, thus, the accuracy of the target cell limited approximation used within influenza viral kinetic models have been questioned for several years [58,[82][83][84][85][86][87]. The data here corroborate the use of this approximation, which does not define what limits the availability of target cells.
The data also agree with the model that there are few infected cells during the time when viral loads are growing most rapidly (0-2 d pi; Figs 1B and 4B). We previously used this information to derive approximations for the model and gain deeper insight into how each parameter influences the kinetics [88], which has been utilized by numerous studies [15,17,89,90]. Further, the data support the model's hypothesis that that there is minimal clearance of infected cells prior to CD8 E infiltration (Fig 3A-B). The knowledge of the model's accuracy and of the spatial growth throughout the lung should aid investigation into the mechanisms that limit virus growth during the early stages of the infection.
Employing targeted model-driven experimental designs to examine and validate the predictions of models like the one presented here is pivotal to elucidating the mechanisms of infection spread and clearance [57,58]. Determining the factors that influence disease severity/weight loss is the first step to understand the disproportionate mortality in at-risk populations (e.g., elderly) and to improve therapeutic design. This is particularly important because current antivirals alleviate 12/35 symptoms but do not effectively lower viral loads [91][92][93][94][95][96]. Validated models like the one here have enormous predictive capabilities and should prove useful in forecasting infection dynamics for a variety of scenarios. These tools and analyses provide a more meaningful interpretation of infection data and a deeper understanding of the progression and resolution of the disease, which will undoubtedly aid our ability to make robust clinical decisions.

Use of experimental animals
All experimental procedures were approved by the Animal Care and Use Committee at St. Jude

Children's Research Hospital under relevant institutional and American Veterinary Medical
Association guidelines and were performed in a biosafety level 2 facility that is accredited by AALAAS.

Mice
Adult (

Infection experiments
The viral infectious dose (TCID 50 ) was determined by interpolation using the method of Reed and Muench [97]

Lung titers
For each animal, viral titers were obtained using serial dilutions on MDCK monolayers and normalized to the total volume of lung homogenate supernatant. The resulting viral loads are shown in Fig 1B and were previously published and utilized for calibration of the densitydependent model (Eq (1)-(4); see below) [15].

Lung histomorphometry and immunohistochemistry
The lungs from IAV infected mice were fixed via intratracheal infusion and then immersion in

Density-dependent model
We previously developed a density-dependent (DD) viral kinetic model that describes the biphasic decline of viral loads [15]. This model tracks 4 populations: susceptible epithelial ("target") cells (T ), two classes of infected cells (I 1 and I 2 ), and virus (V ) [15].
In this model, target cells become infected with virus at rate βV per day. Once infected, target cells enter an eclipse phase (I 1 ) before transitioning at rate k per day to a productively-infected state (I 2 ). These infected cells produce virus at rate p TCID 50 per infected cell per day. Virus is cleared at rate c per day and virus-producing infected cells (I 2 ) are cleared according to the function δ d (I 2 ) = δ d /(K δ + I 2 ), where δ d /K δ is the maximum per day rate of infected cell clearance and K δ is the half-saturation constant. This model provides a close fit to the viral load data in Fig 1B and replicates the biphasic viral load decline [15].

CD8 + T cell model
To examine the contribution of CD8 + T cells to the biphasic viral load decay, we expanded the DD model (Eq (1)-(4)) to include two mechanisms of infected cell clearance (non-specific clearance (δ) and CD8 + T cell-mediated clearance (δ E (I 2 , E))) and two CD8 + T cell populations: effector (E, denoted CD8 E ) and memory (E M , denoted CD8 M ) CD8 + T cells. The model is given by Eq (5)- (10).
In this model, virus-producing infected cells (I 2 ) are cleared by non-specific mechanisms 16/35 (e.g. apoptosis and innate immune responses) at a constant rate δ per day and by CD8 E at rate δ E (I 2 , E) = δ E E/(K δ E + I 2 ) per day, where the rate of infected cell clearance is δ E /K δ E per CD8 E per day and K δ E is the half-saturation constant. The CD8 E -mediated clearance rate (δ E (I 2 , E)) is dependent on the density of infected cells and is similar to the infected cell clearance term in the DD model (see Eq (3)) [15]. Similar density-dependent forms have also been used in models that describe the CD8 + T cell response to other virus infections [31,59].
Models that exclude this density-dependence were examined, but these models resulted in a poor fit to the data. The model assumes that CD8 E infiltrate the lung proportional to infected cells at rate ξ(E) = ξ/(K E + E) CD8 E per cell per day, which is down-regulated by the CD8 E already present in the lung. The associated half-saturation constant is K E . Similar terms for CD8 E regulation have been used in modeling HIV infections [31,98] and in models that examine CD8 + T cell proliferation mechanisms [99].

Parameter estimation
Given a parameter set θ, the cost C(θ) was minimized across parameter ranges using an Adaptive Simulated Annealing (ASA) global optimization algorithm [15] to compare experimental and predicted values of log 10 TCID 50 /lung virus (V ) and of log 10 total CD8 + T cells/lung (Ê = E + E M +Ê 0 , whereÊ 0 is the initial number of CD8 + T cells at 0 d pi). The cost function is defined by is the viral load data, (t i , e i,j ) is the CD8 + T cell data, and V (θ, t i ) and E(θ, t i ) are the corresponding model predictions. Here, s E = (v max − v min )/(e max − e min ) is a scaling factor, and γ i = J i+1 J i−1 where J i is the number of observations at time t i . Errors of the log 10 data were assumed to be normally distributed. To explore and visualize the regions of parameters consistent with the model, we fit Eq (5)- (10) to 2000 bootstrap replicates of the data. If the fit was within χ 2 = 0.05 of the best-fit, then the bootstrap was considered successful [15,89,105]. For each best-fit estimate, we provide 95% confidence intervals (CI) obtained from the bootstrap replicates (Table 1). All calculations were performed in Python using the simanneal package [106] followed by a L-BFGS-B [107,108] deterministic minimization through SciPy's minimize function. SciPy integrate.ode using lsoda and PyDDE [109] were used as the ODE and DDE solvers.
The rate of non-specific infected cell clearance (δ) was given limits of 0 The initial number of target cells (T (0)) was set to 10 7 cells [15,89,105]. The initial number of infected cells I 1 (0) was set to 75 cells to reflect an initial dose of 75 TCID 50 [15]. We previously found that estimating I 1 (0), fixing V (0) = 75 TCID 50 , or estimating V (0) did not improve the fit and could not be statistically justified [15]. The initial number of productively infected cells (I 2 (0)), the initial free virus (V (0)), and the initial number of CD8 E (E(0)) and CD8 M (E M (0)) were set to 0.

Linear regression
The function polyfit in MATLAB was used to perform linear regression on the percent active lesioned area, the percent inactive lesioned area, and the CD8 + T cells during the expansion phase (5-8 d pi) and the contraction phase (9-10 d pi). Linear fits are shown in Figure S6.

Area under the curve
The function cumtrapz in MATLAB was used to estimate the cumulative area under the curve (CAUC) for the infected cells (I 2 ) for the best-fit model solution.

Author contributions
Conceptualization: Amber M. Smith  (T ) are infected at rate βV . Infected cells (I 1 ) undergo an eclipse phase and transition to become productively-infected cells (I 2 ) at rate k. Virus (V ) is produced by infected cells at rate p and is cleared at rate c. Infected cells are cleared at rate δ by non-specific mechanisms and at rate δ E (I 2 , E) by effector CD8 + T cells (E; denoted CD8 E ). The dashed lines represent interactions between infected cells and CD8 E . CD8 E infiltration (ξ(E) = ξ/(K E + E)) is proportional to infected cells and is limited by CD8 E with half-saturation constant K E . CD8 E expansion (η) occurs proportional to infected cells τ E days ago. Memory CD8 + T cell (E M ; denoted CD8 M ) generation occurs at rate ζ and proportional to CD8 E τ M days ago. (B) Fit of the CD8 + T cell model (Eq (5)-(10)) to virus and CD8 + T cells from the lungs of mice infected with 75 TCID 50 PR8 (10 mice per time point). The total number of CD8 + T cells iŝ E = E + E M +Ê 0 . The solid black line is the optimal solution and the gray shading is the minimum and maximum of the model solution using parameter sets within the 95% CIs. Parameters are given in Table 1. Data are shown as mean ± standard deviation.  [15] or the CD8 + T cell viral kinetic model (Eq (5)-(10)) to viral titers and CD8 + T cell counts from mice infected with 75 TCID 50 PR8. (A) Comparison of parameters that were consistent between the DD model (red) and the CD8 + T cell model (black). Correlations are evident between parameters relating to the rates of virus infectivity (β), virus production (p), and virus clearance (c). However, the strength of the correlation was significantly reduced in the CD8 + T cell model. The eclipse phase parameter (k) is not well-defined in either model. (B) The rates of infected cell clearance by non-specific mechanisms (δ) and by CD8 E (δ E ) are slightly negatively correlated. (C) Additional correlations were present between the rates of CD8 E clearance (d E ), expansion (η), and memory generation (ζ). The axes limits reflect imposed bounds. Additional ensemble plots are in Figs S1-S2. , and CD8 E (E). (D) Solutions of the CD8 + T cell model (Eq (5)-(10)) for virus (V ) and total CD8 + T cells (Ê) using the best-fit parameters (black line) and when varying the CD8 E expansion rate (η; magenta lines) to illustrate how different total CD8 + T cell magnitudes alter infection duration. The magenta lines are solutions from when the percentÊ max relative toÊ max from the best-fit solution was 42% (solid line), 39.2% (dash-dotted line), 39.1% (dashed line), or 37% (dotted line). (E) The time at which infected cells reach the half-saturation constant (I 2 = K δ E ; gray circles) and the infection duration (time where log 10 V = 0; black diamonds) are shown for the various CD8 + T cell magnitudes. The gray line between these points is the time required to eliminate K δ E infected cells and achieve complete resolution of the infection (log 10 V = 0). and disease severity. Viral loads and weight loss are the two most easily measured pieces of data. Our analysis relates these through mathematical models. That is, given viral loads, the CD8 + T cell model (Eq (5)-(10)) can be used to predict the kinetics of infected cells and CD8 + T cells. The cumulative area under the curve (CAUC) of the predicted infected cell dynamics (I 2 ) yields an estimate of the percent lung infected (active lesion) while the predicted relative CD8 E dynamics (E) yield an estimate of the percent lung resolved (inactive lesion). The total amount of lung involved (% lung infected and % lung resolved) can then be used to estimate weight loss through the function L(W ). These connections can be reversed and weight loss used to predict viral load kinetics.