Dynamically linking influenza virus infection kinetics, lung injury, inflammation, and disease severity

  1. Margaret A Myers
  2. Amanda P Smith
  3. Lindey C Lane
  4. David J Moquin
  5. Rosemary Aogo
  6. Stacie Woolard
  7. Paul Thomas
  8. Peter Vogel
  9. Amber M Smith  Is a corresponding author
  1. Department of Pediatrics, University of Tennessee Health Science Center, United States
  2. Department of Anesthesiology, Washington University School of Medicine, United States
  3. Flow Cytometry Core, St. Jude Children's Research Hospital, United States
  4. Department of Immunology, St. Jude Children’s Research Hospital, United States
  5. Department of Pathology, St. Jude Children's Research Hospital, United States

Abstract

Influenza viruses cause a significant amount of morbidity and mortality. Understanding host immune control efficacy and how different factors influence lung injury and disease severity are critical. We established and validated dynamical connections between viral loads, infected cells, CD8+ T cells, lung injury, inflammation, and disease severity using an integrative mathematical model-experiment exchange. Our results showed that the dynamics of inflammation and virus-inflicted lung injury are distinct and nonlinearly related to disease severity, and that these two pathologic measurements can be independently predicted using the model-derived infected cell dynamics. Our findings further indicated that the relative CD8+ T cell dynamics paralleled the percent of the lung that had resolved with the rate of CD8+ T cell-mediated clearance rapidly accelerating by over 48,000 times in 2 days. This complimented our analyses showing a negative correlation between the efficacy of innate and adaptive immune-mediated infected cell clearance, and that infection duration was driven by CD8+ T cell magnitude rather than efficacy and could be significantly prolonged if the ratio of CD8+ T cells to infected cells was sufficiently low. These links between important pathogen kinetics and host pathology enhance our ability to forecast disease progression, potential complications, and therapeutic efficacy.

Introduction

Over 15 million respiratory infections and 200,000 hospitalizations result from influenza A viruses (IAVs) each year (Thompson et al., 2004; Simonsen et al., 2000; Taubenberger and Morens, 2008; Medina and García-Sastre, 2011). 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 (Duan and Thomas, 2016; Moskophidis and Kioussis, 1998; Mauad et al., 2010; Rygiel et al., 2009; La Gruta et al., 2007). High viral loads can also play a role in disease progression (Boon et al., 2011; de Jong et al., 2006), but these do not always correlate with the strength of the host response or with disease severity (Granados et al., 2017; Marathe et al., 2016; Toapanta and Ross, 2009; Smith et al., 2019; Gao et al., 2013). 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.

During IAV infection in both humans and animals, viral loads increase rapidly for the first 1–2 days of infection before reaching a peak (e.g., as in Smith et al., 2018; Baccam et al., 2006; Miao et al., 2010; Toapanta and Ross, 2009; Carrat et al., 2008; Srivastava et al., 2009; Gao et al., 2013). If the host is infected with a novel strain or has no prior immunity, viral loads in the lung then begin to decline, first slowly (sometimes <1 log10 TCID50/d ) then rapidly (> 4-5 log10 TCID50/d) (Smith et al., 2018; Gao et al., 2013). We previously quantified this biphasic viral decline with a mathematical model, which indicated that the rate of infected cell clearance increases as the density of infected cells decreases (Smith et al., 2018). The timing of the second, rapid viral decay phase coincides with the expansion of CD8+ T cells, which are the primary cell responsible for clearing infected cells and resolving the infection (Zhang and Bevan, 2011; Chen et al., 2018; McMichael et al., 1983; La Gruta and Turner, 2014; Kreijtz et al., 2011), and, to a lesser extent, neutralizing antibodies (Chen et al., 2018; Kreijtz et al., 2011; Fang and Sigal, 2005; Wang et al., 2015; McMichael et al., 1983) and cytotoxic CD4+ T cells (Wilkinson et al., 2012). For the CD8+ T cell response, in particular, it remains unclear whether the efficacy of these cells is dictated by their own density (Graw and Regoes, 2009; Wiedemann et al., 2006), infected cell density (Merrill, 1982; Caramalho et al., 2009; Halle et al., 2016), or both (Gadhamsetty et al., 2014). While quantifying dynamically changing CD8+ T cell efficacy is difficult in vitro and in vivo, the question is ripe for in silico investigation. Several modeling studies have described CD8+ T cell-mediated infected cell clearance for various viral infections in humans and animals, including IAV, HIV, and LCMV (e.g., as in Ganusov et al., 2011; De Boer and Perelson, 1998; Perelson and Bell, 1982; Miao et al., 2010; Cao et al., 2016; Lee et al., 2009; Price et al., 2015; Graw and Regoes, 2009; Merrill, 1982; Gadhamsetty et al., 2014; Baral et al., 2019; Bonhoeffer et al., 2000; Conway and Perelson, 2015). Some of these studies have also attempted to link this response and other immune responses to inflammation or disease severity (Price et al., 2015; Baral et al., 2019; Manchanda et al., 2014), but have not yet found the appropriate mathematical relation with the available data. In addition, for IAV infections in particular, varied efficiency of the CD8+ T cell response throughout the course of infection, their early antigen-specific and lung-resident responses (Yoon et al., 2007; McGill and Legge, 2009; Keating et al., 2018), and the balanced consequences on viral loads and host pathology have not yet been investigated in detail.

A better understanding of infected cell clearance may also yield insight into the damage induced to the lung and the ensuing immunopathology during IAV infection. In general, widespread alveolar disease is observed in patients who succumb to the infection (Bautista et al., 2010; Shieh et al., 2010; Gao et al., 2013), and CT scans show bilateral and multi-focal areas of lung damage (e.g. as in Kohr et al., 2010; Li et al., 2011; Ajlan et al., 2009; Soto-Abraham et al., 2009; Gill et al., 2010; Li et al., 2012). Further, hospitalized patients that died as a result from infection with the 2009 H1N1 influenza virus had large numbers of CD8+ T cells present in their lung tissue (Mauad et al., 2010). Large pulmonary populations of CD8+ T cells contribute to lung injury by targeting IAV-infected cells (reviewed in Duan and Thomas, 2016) and inducing ‘bystander damage’ to uninfected epithelial cells (van de Sandt et al., 2017). However, their population levels do not necessarily indicate active infected cell removal or immunopathology as their extended presence in the lung also contributes to surveillance and repair (Matheu et al., 2013; Sun et al., 2009). Macrophages and neutrophils also persist in the lung following viral clearance, and their role in inflammation and tissue damage is well documented (reviewed in Watanabe et al., 2019; La Gruta et al., 2007). A favorable outcome requires a balance between immune-mediated protection and pathology, and the viral-immunological landscape that drives disease and the markers that distinguish mild from severe influenza are complex. This is particularly true in humans (Tang et al., 2019) where obtaining high quality data from the lung is challenging, and viral loads in the upper respiratory tract do not always reflect the lower respiratory tract environment (Feikin et al., 2017; Ai et al., 2020; Wong et al., 2020; Li et al., 2012; Ruuskanen et al., 2011).

The accumulation of damage to the epithelium during IAV infection, either from virally-induced cell lysis or immune-mediated effects, is relatively understudied. We and others have modeled the lung damage and inflammation inflicted during pulmonary infections (e.g., as in Smith et al., 2011b; Reynolds et al., 2006; Dunster et al., 2014; Ramirez-Zuniga et al., 2019; Baral et al., 2019) but did so without sufficient data to constrain the model formulations. The difficulty in measuring the dynamics of infected cells and in establishing how damage correlates to specific host responses has been the primary impediment. However, recent technological advances, including the use of reporter viruses (Manicassamy et al., 2010; Tran et al., 2013; Karlsson et al., 2015; Luker and Luker, 2010) and lung histomorphometry (Marathe et al., 2016; Marathe et al., 2017; Sartorius et al., 2007; Zachariadis et al., 2006; Boyd et al., 2020) within animal models, have provided an opportunity to acquire these types of measurements. Whole lung histomorphometry, which is broadly defined as a quantitative microscopic measurement of tissue shape, has recently increased in use due to the ability to directly stain, visualize, and quantify areas of infected cells. 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 inflammation and damage, which are generally thought to be markers of disease severity. In addition, it should bolster the development of robust, predictive computational models, which have historically lacked constraint to these types of data.

In general, disease severity may not be directly correlated to viral loads or 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) (Carrat et al., 2008; Lee and Lee, 2010; Han et al., 2018). Some symptoms (e.g., fever) are cytokine mediated (Monto et al., 2000), but respiratory symptoms often last longer and can remain following viral clearance (Lee and Lee, 2010; Eccles, 2005). This is particularly true when there is scarring of the lung tissue (Li et al., 2012). 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 (Trammell and Toth, 2011; Bouvier and Lowen, 2010; Parzych et al., 2013). Animal weights generally drop slowly in the first several days during an IAV infection and more rapidly as the infection begins to resolve (Srivastava et al., 2009; Huang et al., 2017). 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 (Smith et al., 2018). Because weight loss can occur following resolution of the infection, immune-mediated pathology is thought to play a role (Lauder et al., 2011; Xu et al., 2004; Ostler et al., 2002; Parzych et al., 2013; Duan and Thomas, 2016). Host and pathogen factors, such as age, viral proteins, and inoculum size, have also been shown to influence disease progression (Toapanta and Ross, 2009; Lu et al., 2018; Smith et al., 2019; McAuley et al., 2007). While the precise contribution of different factors to IAV-associated disease and mortality remain elusive, having tools that can link these with disease severity and decipher correlation from causation would improve our ability to effectively predict, understand, and treat the disease.

To gain deeper insight into the dynamics of viral resolution during primary pulmonary influenza infection and investigate the connection between viral loads, damage, inflammation, and disease severity, we developed and validated a mathematical model that describes viral kinetics and CD8+ T cell-mediated infected cell clearance. The model accurately predicted the dynamics of effector and memory CD8+ phenotypes, agreed with our previous results that infected cells are cleared in a density-dependent manner with CD8 efficiency rapidly accelerating by over 48,000 times in 2 d, and illuminated tradeoffs between the innate and adaptive immune responses. Our model predicted that infection duration is dependent on the magnitude of CD8+ T cells rather than their efficacy, which we verified by depleting CD8+ T cells. Quantitative whole lung histomorphometry showed that infected areas increase during the first 6 d of infection before clearing rapidly, that the relative number of CD8+ T cells corresponded to the cleared area, and corroborated the model-predicted infected cell dynamics. These data also showed that the infection-induced lung injury is distinct from inflammation, and that each correlates to different immune responses and could be predicted using independent mathematical equations. In addition, the data revealed nonlinear connections between these two pathological measurements and disease severity. Together, our 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 the dynamics of viral loads, infected cells, lung injury, inflammation, and disease severity.

Results

Virus and CD8+ T cell kinetics

In animals infected with 75 TCID50 PR8, virus rapidly infects cells or is cleared within 4 hr pi and is undetectable in all animals at this time (Figure 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 log10 TCID50/d) (Smith et al., 2018). In the second phase (7–8 d pi), virus is cleared rapidly with a loss of 4–5 log10 TCID50 in 1–2 d (average of –3.8 log10 TCID50/d) (Smith et al., 2018). CD8+ T cells remain at their baseline level from 0–2 d pi, where ∼15% are located in the circulating blood, ∼75% in the lung vasculature, and ∼10% in the lung parenchyma (Figure 1—figure supplement 1) with the majority in the lung as CD43CD127+. IAV-specific CD8+ T cells that are recently primed (CD25+CD43+) or effector cells (CD25CD43+) (Keating et al., 2018) begin appearing in the lung and increase slightly from 2–3 d pi but remain at low levels. These cells are thought to proliferate within the lung at least once by 4 d pi (McGill and Legge, 2009), and their populations briefly contract (3–5 d pi) before expanding rapidly (5–8 d pi). During the expansion phase, >95% of recovered cells are in the lung (∼70% in the parenchyma and ∼30% in the vasculature) with CD25+CD43+ and CD25CD43+ as the predominant phenotypes. This is in accordance with prior studies that showed these phenotypes are IAV-specific and that their expansion dynamics in the lung were not altered by removing blood-borne CD8+ T cells from the analysis (Anderson et al., 2012; Keating et al., 2018).

Figure 1 with 2 supplements see all
Schematic and fit of the CD8+ T cell viral kinetic model.

(A) Schematic of the CD8+ T cell model in Equation (1)-(6). In the model, target cells (T) are infected at rate βV. Infected cells (I1) undergo an eclipse phase and transition to become productively-infected cells (I2) 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(I2,E) by effector CD8+ T cells (E; denoted CD8E). The dashed lines represent interactions between infected cells and CD8E. Initial CD8E influx (ξ(E)=ξ/(KE+E)) is proportional to infected cells and is limited by CD8E with half-saturation constant KE. CD8E expansion (η) occurs proportional to infected cells τE days ago. Memory CD8+ T cell (EM; denoted CD8M) generation occurs at rate ζ and proportional to CD8E τM days ago. (B) Fit of the CD8+ T cell model (Equation (1)-(6)) to virus and total CD8+ T cells from the lungs of mice infected with 75 TCID50 PR8 (10 mice per time point). The total number of CD8+ T cells is E^=E+EM+E^0. (C) Total CD8+ T cells in the lung parenchyma (gray circles) and overlay of the model predicted values (E+EM). (D) Fit of the model to virus, CD25CD43+ CD8+ T cells (cyan diamonds; E), and CD43CD127+CD103+CD69+ CD8+ T cells (blue diamonds; EM) (five mice per time point). The solid and dashed black lines are the optimal solutions and the gray shading is are the model solutions using parameter sets within the 95% CIs. Parameters are given in Table 1. Data are shown as mean ± standard deviation.

CD8+ T cell expansion corresponds to the second viral decay phase with sixty percent of mice clearing the infection by 8 d pi and the other forty percent by 9 d pi (Figure 1). Most CD43+CD8+ T cell phenotypes decline following viral clearance (8–10 d pi) but do not return to their baseline level by 12 d pi. Long-lived antigen-specific memory phenotypes down regulate CD43 (Jones et al., 1994; Harrington et al., 2000; Onami et al., 2002) and gradually increase substantially beginning at 9 d pi with most as CD25, which is qualitatively similar to other studies (Harrington et al., 2000). At 12 d pi, ∼55% of CD8+ T cells remain in the lung parenchyma and ∼20% in the circulating blood indicating exit from the lung (Figure 1—figure supplement 1).

Viral kinetic model with density-dependent CD8+ T cell-mediated clearance

We previously described the influenza viral load kinetics and biphasic decline during primary pulmonary infection using a density-dependent (DD) model (Smith et al., 2018), which assumes that the rate of infected cell clearance increases as the density of infected cells decreases (i.e., δd(I2)=δd/(Kδ+I2)). Because the rapid decay of virus is thought to be due to the clearance of infected cells by CD8+ T cells, it is unknown if early CD8+ T cell presence contributes to infected cell clearance, and no model has captured the entire CD8+ T cell time course, we developed a model that describes the dynamics of these cells and their efficiency in resolving the infection (Equation (1)-(6); Figure 1A). The model includes equations for effector (E, denoted CD8E) and memory (EM, denoted CD8M) CD8+ T cells, and two mechanisms of infected cell clearance. The first mechanism is from unspecified, innate mechanisms, which is relatively constant (δ) and primarily acts during the first viral decay phase (2–6 d pi). The second is the CD8E-mediated infected cell clearance, which occurs at a rate that increases as the density of infected cells decreases (δE(I2,E)=δEE/(KδE+I2)) and primarily acts during the rapid, second viral decay phase (7–8 d pi). Excluding this density dependence entirely resulted in a significant and premature decline in viral loads, which disagreed with the experimental data. We also tested whether the density dependence could be included in the CD8+ T cell expansion rate rather than in the infected cell clearance rate (see Appendix 1) as other models have done (e.g., as in Baral et al., 2019; Bonhoeffer et al., 2000; Conway and Perelson, 2015). This modification yielded a close fit to the CD8+ T cell data at 6–10 d pi but not at the early time points (Appendix 1—figure 1). In addition, the viral load data was underestimated at 7 d pi causing the solution to miss the rapid decline between 7–8 d pi and result in a statistically poorer fit. Thus, retaining the density-dependence in the rate of infected cell clearance most accurately captured the entire dataset. The model includes terms for the initial CD8E influx at 2–3 d pi (ξI2/(KE+E)) and for CD8E expansion (ηEI2(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 the initial CD8E influx is regulated by their own population (i.e., ξ(E)=ξ/(KE+E)). In both terms, the increase is proportional to the number of infected cells. Although memory CD8+ T cells were not the primary focus here, it was necessary to include the CD8M population because CD8+ T cells are at a significantly higher level at 10–12 d pi than at 0 d pi (Figure 1B). Fitting the model simultaneously to viral loads and total CD8+ T cells from the (non-perfused) lungs of infected animals illustrated the accuracy of the model (Figure 1B). The resulting parameter values, 95% confidence intervals (CIs), ensembles, and histograms are given in Table 1, Figure 2, and Figure 2—figure supplements 12.

Figure 2 with 2 supplements see all
Parameter ensembles and histograms.

Parameter ensembles (A) and histograms (B) resulting from fitting the CD8+ T cell viral kinetic model (Equation (1)-(6)) to viral titers and total CD8+ T cells from mice infected with 75 TCID50 PR8. (A) The rates of infected cell clearance by non-specific mechanisms (δ) and by CD8E (δE) are slightly negatively correlated. Correlations were also present between the rates of CD8E clearance (dE), CD8E expansion (η), and CD8M generation (ζ). The axes limits reflect imposed bounds. Additional ensemble plots are in Figure 2—figure supplements 12. (B) The histograms show that the majority of parameters are well-defined with the exception of the eclipse phase transition rate (k), one of the half-saturation constants (KE), and the CD8M generation delay (τM).

Table 1
CD8+ T cell model parameters.

Parameters and 95% confidence intervals obtained from fitting the CD8+ T cell model (Equation (1)-(6)) to viral titers and total CD8+ T cells (‘Total CD8’) or viral titers, CD25CD43+CD8+ T cells, and CD43CD127+CD103+CD69+CD8+ T cells (‘Specific CD8E,M Phenotypes’) from mice infected with 75 TCID50 PR8. CD8E and CD8M denote effector (E) and memory (EM) CD8+ T cells, respectively. The total number of CD8+ T cells is E^=E+EM+E0^ and is denoted by CD8.

ParameterDescriptionUnitsTotal CD8Specific CD8E,M phenotypes
Value95% CIValue95% CI
βVirus infectivityTCID50-1d-16.2 × 10–5[5.3 × 10–6, 1.0 × 10–4]3.7 × 10–5[1.1 × 10–5, 9.4 × 10–5]
kEclipse phase transitiond-14.0[4.0, 6.0]5.1[4.0, 6.0]
pVirus productionTCID50cell-1d-11.0[5.8 × 10–1 × 1.1 × 102]1.5[0.73, 13.6]
cVirus clearanced-19.4[5.6, 9.5 × 102]12.1[5.8, 17.5]
δInfected cell clearanced-12.4 × 10–1[1.0 × 10–1, 6.6 × 10–1]3.0 × 10–1[1.9 × 10–1, 5.9 × 10–1]
δEInfected cell clearance by CD8EcellsCD8E-1d-11.9[3.3. × 10–1, 2.0]5.7[1.7, 8.5]
KδEHalf-saturation constantcells4.3 × 102[1.0 × 102, 2.9 × 105]1.3 × 102[1.0 × 101, 8.6 × 102]
ξCD8E infiltrationCD8E2 cell1 d12.6 × 104[1.3 × 102, 8.7 × 104]2.9 × 103[8.0 × 102, 1.7 × 104]
KEHalf-saturation constantCD8E8.1 × 105[1.0 × 103, 1.0 × 106]2.2 × 106[1.1 × 106, 8.3 × 106]
ηCD8E expansioncell-1d-12.5 × 10–7[1.6 × 10–8, 6.7 × 10–7]3.5 × 10–7[2.3 × 10–7, 5.2 × 10–7]
τEDelay in CD8E expansiond3.6[2.1, 5.9]3.3[2.6, 3.8]
dECD8E clearanced-11.0[5.1 × 10–2, 2.0]1.1[3.3 × 10–1, 2.5]
ζCD8M generationCD8M CD8E1 d12.2 × 10–1[1.0 × 10–2, 9.4 × 10–1]3.2 × 10–2[1.0 × 10–2, 2.2 × 10–1]
τMDelay in CD8M generationd3.5[3.0, 4.0]3.3[2.1, 5.3]
E^0Baseline CD8 or CD8ECD8 or CD8E4.2 × 105[3.3. × 105, 5.3 × 105]6.0 × 102[3.1 × 102, 1.8 × 102]
E^M0Baseline CD8MCD8M--6.5 × 102[3.2 × 102, 1.5 × 103]
T(0)Initial uninfected cellscells1 × 107-1 × 107-
I1(0)Initial infected cellscells75-75-
I2(0)Initial infected cellscells0-0-
V(0)Initial virusTCID500-0-
E(0)Initial CD8ECD8E0-0-
EM(0)Initial CD8MCD8M0-0-

Plotting the model predicted dynamics for the lung-specific CD8+ T cells (CD8L=E+EM) illustrated the accuracy of the model in predicting their dynamics within the lung parenchyma without fitting to these data (Figure 1C). One benefit of using the total CD8+ T cells is that the model automatically deduces the dynamics of effector-mediated killing and memory generation without needing to specify which phenotypes might be involved in these processes as they may be dynamically changing. However, the rates of expansion, contraction, and infected cell clearance may be different if only certain phenotypes are engaged. Thus, we examined whether the model could fit the dynamics of the predominant effector (CD25CD43+) and memory (CD43CD127+CD103+CD69+) phenotypes. Re-fitting the model to these data suggested that no changes to model formulation were needed and there were only small alterations to select CD8-specific parameter values (Figure 1D; Table 1).

Plotting the model ensembles revealed a correlation between the two infected cell clearance parameters (δ and δE; Figure 2B), which represent the efficacy of the non-specific immune response and the CD8+ T cell response, respectively. Performing a sensitivity analysis showed that the viral load dynamics do not change substantially when these parameters are increased or decreased (Appendix 2—figure 1) . However, decreasing the rate of non-specific infected cell clearance (i.e., lower δ) resulted in a significant increase in the number of CD8E due to the small increase in the number of infected cells (Figure 1). Even with a larger CD8E population, recovery was delayed by only ~0.1 d. Given the correlation between δ and δE (Figure 2B), a more efficient CD8E 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 (Appendix 2—figure 1; Smith et al., 2018), 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 CD8E expansion (η) (Figure 1, Appendix 1—figure 1; 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 CD8E expansion (η) and clearance (dE) were slightly correlated, indicating a balance between these two processes (Figure 2A). This correlation and the sensitivity of η produced model dynamics that were also sensitive to changes in the CD8E clearance rate (dE) (Appendix 3—figure 2). As expected, the rates of CD8M generation (ζ) and CD8E clearance (dE) were correlated (Figure 2A). It has been estimated that approximately 5–10% of effector CD8+ T cells survive to become a long-lasting memory population (Kaech et al., 2002). Despite the inability to distinguish between CD8E and CD8M in the total CD8+ T cell data, the model predicts that 17% of CD8E transitioned to a memory class by 15 d pi. When considering only the CD25CD43+/– effector and memory phenotypes, the model estimates this value to be ~7%. 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 included in Appendix 2.

Density-dependent infected cell clearance

Given the accuracy of the model, we next sought to gain further insight into the nonlinear dynamics of CD8+ T cell-mediated infected cell clearance. Plotting the clearance rate (δE(I2,E)=δEE/(KδE+I2)) as a function of infected cells (I2) and CD8E (E) (Figure 3) confirmed that there is minimal contribution from CD8E-mediated clearance to viral load kinetics or infected cell kinetics prior to 7 d pi (Figure 3A–C, markers a–b). At the initiation of the second decay phase (7 d pi), the clearance rate is δE(I2,E)=3.5/d (Figure 3A–B, marker c). As the infected cell density declines toward the half-saturation constant (KδE=4.3×102cells), the clearance rate increases rapidly to a maximum of 4830/d (Figure 3A–C, markers d–g). The model predicts that there are 6 × 105 infected cells remaining at 7 d pi, which can be eliminated by CD8E in 6.7 h.

Density-dependent infected cell clearance by CD8+ T cells and their impact on recovery time.

(A) The rate of CD8E-mediated infected cell clearance (δE(I2,E)=δEE/(KδE+I2)) plotted as a function of infected cells (I2) and effector CD8+ T cells (E; CD8E). The colored markers (denoted a–g) indicate the infected cell clearance rate that corresponds to different time points during the infection for the best-fit solution. (B) Values of δE(I2,E) for the indicated time points associated with the markers a–g. (C) Corresponding locations of the various δE(I2,E) values (markers a–g) on the best-fit solution of the CD8+ T cell model for virus (V), infected cells (I2), and CD8E (E). (D) Solutions of the CD8+ T cell model (Equation (1)-(6)) for virus (V) and total CD8+ T cells (E^) using the best-fit parameters (black line) and when varying the CD8E 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 E^max relative to E^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 (I2=KδE; gray circles) and the infection duration (time where log10V=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 (log10V=0). (F) Fit of the CD8+ T cell model (Equation (1)-(6)) to viral loads and CD8+ T cells (magenta diamonds) following depletion at −2, 0, 3, and 7 d pi (magenta arrows). The best model (Supplementary file 2) resulted in fewer target cells (T0), a lower CD8E influx (ξ), and a higher CD8E expansion rate (η). All other parameters were fixed to the best-fit value in Table 1. The solid lines are E^=E+EM+E^0 and the dashed lines are E for the cases where CD8+ T cells were depleted (magenta) and where they were not depleted (black). (G) Comparison of the log10 ratio of virus to CD69+CD8+ T cells with and without CD8+ T cell depletion (magenta and black, respectively). All data are shown as mean ± standard deviation.

To explore how recovery time is altered by varying CD8E levels, we examined the resulting dynamics from increasing or decreasing the rate of CD8E expansion (η). When η was increased by 50%, the CD8E population increased by a substantial 670% (Appendix 3—figure 1). However, this was insufficient to significantly shorten the infection (8.4 d versus 7.8 d). The infection duration could be reduced if CD8E expansion began earlier (i.e., smaller τE; Appendix 3—figure 2). Although recovery is not significantly expedited by a larger CD8E population, our model predicted that the infection would be dramatically prolonged if these cells are sufficiently depleted (Figure 3D–E and Figure 1). This in silico analysis revealed a bifurcation in recovery time, such that the infection is either resolved within ∼15 d pi or may last up to ∼45 d if CD8E are below a critical magnitude required to resolve the infection (Figure 3D–E). The critical number of total CD8+ T cells needed for successful viral clearance was E^maxcrit=7.4×105 CD8, which was 39.2% of the maximum number of CD8+ T cells obtained from the best-fit solution (i.e., E^max=1.9×106 CD8). This corresponds to 17% of CD8E (i.e., Emaxcrit=2.3×105 CD8E; Figure 3D–E). The model analysis indicated that decreasing the total number of CD8+ T cells by as little as 0.1% from this critical level (i.e., 39.2% to 39.1%) lengthened the infection from 15 d pi to 25 d pi (Figure 3D).

Dynamical changes during CD8+ T cell depletion

To further test our model formulation and identify how viral load kinetics are altered with dynamically changing CD8+ T cells, we infected groups of mice with 75 TCID50 and depleted these cells at −2 d, 0 d, 3 d, and 7 d pi with an anti-CD8α antibody (clone 2.43; Figure 3F). CD8+ T cells were reduced with over 99% efficiency and only 1.3% remained 2 d after depletion in the absence of infection (i.e., at 0 d pi; Figure 3F). CD8+ T cells remained >1 log10 lower throughout the infection. Unexpectedly, the corresponding viral loads were significantly lower at 2 d pi (p=0.02) and consistently lower at other time points early in the infection (4 d pi (p=0.1) and 6 d pi (p=0.2)). As expected and predicted by our mathematical model, viral loads were significantly higher at 8 d pi (4.68 log10 TCID50 compared to 1.42 log10 TCID50; p=0.01). By 10 d pi, a sufficient number of CD8+ T cells were present and all animals had cleared the infection (Figure 3F). Interestingly, the number of CD8+ T cells at 10 d pi was only slightly higher than at 8 d pi (5.25 log10 versus 5.02 log10; p=0.064) further highlighting the density-dependent dynamics described above.

Given that viral loads were lower at early time points and that the anti-CD8α antibody is known to cause concentration-dependent changes in CD8E differentiation, activation, and efficacy (Cross et al., 2019) in addition to resulting in death of the cells that would initiate activation of and removal by other immune cells, we did not expect our model to match the data without modulation of parameter values. However, we did expect that no changes to the model formulation would be required. In total, we tested >30 ‘models’ where 1–4 parameters were altered and found one model that was significantly better according to the AIC (Supplementary file 1). In that model, the initial number of target cells (T0) was 2.5x lower (~4 × 106 versus 1 × 107 cells), the rate of initial CD8E influx (ξ) was 2x lower (1.3 × 104 versus 2.6×104CD8E2 cell1 d1), and the rate of CD8E expansion (η) was 4x higher (1 × 10–6 versus 2.5×107cell1 d1). The second best model had a lower cost but was penalized by an additional parameter. That model suggested similar results but replaced the effect on T0 with a combination of a lower virus production rate (p) and higher non-specific infected cell clearance rate (δ). Both models resulted in fewer infected cells and, thus, have approximately equivalent interpretations.

To further examine these findings, we plotted the ratio of virus to activated (CD69+) CD8+ T cells (Figure 3G). At 2 d pi, the ratio was similar for the depleted and mock-treated groups (p=0.41) suggesting that the depletion-induced reduction in virus was proportional to the reduction in activated CD8+ T cells. However, at 4 d pi, the ratio in the depleted groups was significantly higher (p=0.0028) for otherwise similar levels of virus. This signifies that there were disproportionately low numbers of CD8+ T cells and, thus, reduced influx (i.e., lower ξ). By 6 d pi, the ratios were again similar (p=0.25) implying a higher expansion of these cells (i.e., higher η). Investigation into the predicted target cell reduction is included below.

Modeling lung injury dynamics

To investigate the dynamics of infected cells and their clearance by CD8+ T cells, we quantified these cells and the progression and resolution of lung injury using serial whole lung histomorphometry (Figure 4A). Antigen-positive areas of the lung (‘active’ lesions) were first detectable at 2 d pi (Figure 4A–B), which coincides with the peak in viral loads. The infected area continued to increase in a nonlinear manner until 6 d pi, whereas viral loads remained high until 7 d pi (Figure 1B). At this time, resolution of the infection began and the infected area declined at a rate of ∼28.7%/d between 6–7 d pi (Figure 4—figure supplement 1). Few to no infected cells were present at 8 d pi (Figure 4A). Correspondingly, virus was undetectable in most animals by 8 d pi (Figure 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 (I2). Plotting the percent active lesion against the CAUC of I2 shows that the model accurately reflects the cumulative infected cell dynamics and, thus, the infection progression within the lung (Figure 4B). Plotting the CAUC of I2 for all parameter sets in the 95% CIs and from fitting the model to specific phenotypes further illustrates the accuracy by showing that the heterogeneity in the histomorphometry data is captured (Figure 4B), which is larger than the heterogeneity in viral loads (Figure 1B). Plotting the predicted active lesion dynamics for the model parameterized to the data where CD8+ T cells were depleted suggested that there was a ∼22% reduction in the active lesion (Figure 4B).

Figure 4 with 2 supplements see all
Histomorphometry and inflammation of the IAV-infected lung.

(A) Whole lung sections with histomorphometry showing the areas of influenza NP-positive ‘active’ lesions (red), previously infected ‘inactive’ lesions with minimal antigen-positive debris (green), or mixed active and inactive regions (orange) throughout the infection. Representative images from each group are shown. (B) Percent active lesion (red squares) plotted together with the cumulative area under the curve (CAUC) of the predicted infected cell dynamics (I2) obtained from fitting the CD8+ T cell model. The linear decline in the active lesion (−28.7%/d; see Figure 4—figure supplement 1) was used to estimate the decline after 6 d pi. (C) Percent inactive lesion (green squares) plotted together with the percent maximum CD8E (E/Emax) obtained from fitting the CD8+ T cell model. (D) The total lesion (blue squares) is the addition of the active and inactive lesions. To include all measurements on the same scale, the CAUC of I2 was multiplied by a scaling factor of 14.2% per 1 × 107 cells, and the percent maximum CD8E was multiplied by a scaling factor of 0.46%. (E) Fit of Equation (7) to the alveolar (white triangles) and interstitial (yellow triangles) inflammation scores. The solid black, dashed black, and solid magenta lines are the curves generated using the best-fit parameters obtained from fitting the model to the total CD8+ T cells, the CD25CD43+ and CD43CD127+CD103+CD69+ CD8+ T cells, and the total CD8+ T cells during CD8 depletion, respectively. The gray shading are the curves generated using the 95% CI parameters from fitting the model to the total CD8+ T cells. All data are shown as mean ± standard deviation.

Antigen-negative, previously-infected or damaged areas of the lung (‘inactive’ lesions) are evident beginning at 5 d pi (Figure 4A,C). This resolution of the infection continued from 5-8 d pi, causing a 15.1%/d increase in the area of inactive lesions (Figure 4—figure supplement 1). Following this, healing of the injured lung is apparent as the inactive lesioned area declines (−14.5%/d from 9–10 d pi; Figure 4C and Figure 4—figure supplement 1). These dynamics generally parallel the CD8+ T cell dynamics but are nonlinearly correlated (Figure 4—figure supplement 1). Fitting a line to the CD8+ T cell data from 5-8 d pi indicated that the influx rate of all phenotypes is 4.94 × 105 cells/d, of lung-specific phenotypes is 3.97 × 105 cells/d, and of the CD25CD43+ effector phenotype is 1.98 × 105 cells/d (Figure 4—figure supplement 1). Thus, the model estimates that, on average, every 100,000 total, lung-specific, or CD25CD43+ CD8+ T cells clear ~3.1%, ~3.8%, or ~7.6% of the infected areas within the lung, respectively. During the CD8+ T cell contraction phase, a similar linear regression analysis suggested that these cells decline at rates of ~ 4.13 × 105 CD8/d (total), ~ 2.35 × 104 CD8/d (lung-specific), and ~ 3.83 × 104 CD8/d (CD25CD43+) (Figure 4—figure supplement 1). Similar to the relation discussed above, the dynamics of the damaged areas of the lung corresponded to the dynamics of the percent maximum CD8E (i.e., E/Emax) in the model (Figure 4C). Our model suggested a ∼37% reduction in the inactive lesion during CD8+ T cell depletion (Figure 4C). Adding the predicted dynamics for the active and inactive lesions agrees with the dynamics of the total lesion, but is slightly underestimated when using the model fit to CD25CD43+ CD8+ T cells (Figure 4D). In addition, the predicted total lesion was reduced for the case where CD8+ T cells were depleted (Figure 4D).

Modeling lung inflammation dynamics

In addition to measuring virus-induced lung damage, lung inflammation was scored (Figure 4E). Both alveolar and interstitial inflammation begin to increase at 2 d pi with the sharpest increase between 1–3 d pi. Inflammation continues to increase until 6 d pi with a maximum score of 3.5–4.0 out of 5. Resolution of inflammation was not evident during the time course of our data, which concluded at 10 d pi. This is in contrast to the lung damage inflicted by the virus, which begins declining at 8 d pi and shows that ∼15% of the damage was repaired by 10 d pi. Inflammation was linearly correlated to inflammatory macrophages and to neutrophils (Figure 4—figure supplement 2), which were excluded from the model. However, the knowledge that the model’s predicted infected cell dynamics are accurate and that these cells can produce cytokines and chemokines that attract cells like macrophages and neutrophils suggested that our model could be used to estimate inflammation. To do this, we fit Equation (7) to the inflammation scores while keeping all other parameters fixed to their best-fit values (Table 1). The results suggested that the equation captured the inflammation dynamics, and that the contribution from the initial infected cell class (I1) was α1=4.27per 107cells/d and the contribution from the productively infected cells (I2) was α2=0.87per 107cells/d. For the case where CD8+ T cells were depleted, Equation (7) estimated the inflammation scores to be reduced to ∼1.5 (Figure 4E).

Weight loss relates nonlinearly to lung injury and inflammation

To monitor disease progression, weight loss was measured daily throughout the course of infection (Figure 5). 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; Figure 5A). Weight loss was reduced during CD8+ T cell depletion (Figure 5A). Plotting weight loss against the percent total (active and inactive) lesioned area of the lung shows the similarity in their dynamics (Figure 5A) and revealed a nonlinear relation (Figure 5B). To further quantify their relationship, we fit the saturating function L(W)=lmaxWn/(Wn+Kwn) to these data, where L is the percent total lesioned area, W is percent weight loss, lmax is the maximum rate of the interaction, Kw is the half-saturation constant, and n is the Hill coefficient. This function provided a close fit to the data (R2=0.92; black line in Figure 5B) with best-fit parameters lmax=39.7% total lesioned area, Kw=2.58% weight loss, and n=5.24. Plotting the estimated percent total lesion during CD8+ T cell depletion (Figure 4D) together with the measured weight loss (Figure 5A) also supported the nonlinear relation. Using the same best-fit parameters (Kw=2.58% weight loss and n=5.24) together with the predicted maximum percent total lesion (lmax=17.73% total lesioned area) showed the accuracy of this function (magenta line in Figure 5D). Weight loss was also nonlinearly correlated to the alveolar and interstitial inflammation scores (Figure 5C), although the relation was slightly different compared to that for the lung injury data. Fitting the same function to these data independently for alveolar and interstitial inflammation also provided a close fit (R2=0.98 and R2=0.97; black lines in Figure 5D). Best-fit parameters for alveolar inflammation were lamax=3.63 score, Kaw=1.95% weight loss, and na=3.65. Best-fit parameters for interstitial inflammation were limax=3.40 score, Kiw=1.96% weight loss, and ni=3.15. For the CD8-depleted case, plotting the estimated inflammation (Figure 4E) together with the measured weight loss (Figure 5A) again supported the nonlinear relation. In addition, using the same best-fit parameters (Kaw=1.95% weight loss, na=3.65; Kiw=1.96% weight loss, ni=3.15) together with the predicted maximum inflammation score (limax=1.41 score) showed the accuracy (magenta line in Figure 5D).

Weight loss dynamics and its relation to lung injury and inflammation.

(A) The percent total (active and inactive) lesion (blue squares) plotted together with the percent weight loss (black diamonds) to illustrate their similar dynamics. (B) Fit of a saturating function (L(W)=lmaxWn/(Kwn+Wn); black line) to the mean percent total lesioned area (L) and mean weight loss (W) for all time points (black circles). The best-fit parameters were lmax=39.7% total lesioned area, Kw=2.6% weight loss, and n=5.2. The magenta circles are the predicted percent total lesioned area (Figure 4D) with the corresponding weight loss during CD8 depletion, and the magenta line is the prediction using the best-fit parameters (Kw=2.6% weight loss and n=5.2) together with the maximum predicted percent total lesion (lmax=17.73% total lesioned area). (C) The alveolar (white triangles) and interstitial (yellow triangles) inflammation plotted together with the percent weight loss (black diamonds). (D) Fit of a saturating function (black line) to the mean alveolar (white circles) and interstitial (yellow circles) inflammation scores and mean weight loss for all time points. The best-fit parameters for alveolar inflammation (white circles; black line) were lamax=3.63 score, Kaw=1.95% weight loss, and na=3.65, and for interstitial inflammation were limax=3.40 score, Kiw=1.96% weight loss, and ni=3.15. The magenta circles are the predicted percent inflammation score (Figure 4E) with the corresponding weight loss during CD8 depletion, and the magenta line is the prediction using the best-fit parameters (Kaw=1.95% weight loss, na=3.65; Kiw=1.96% weight loss, ni=3.15) together with the maximum predicted inflammation (lmax=1.41 score).

Discussion

Influenza A virus infections pose a significant threat to human health, and it is crucial to understand and have tools that can predict how the virus spreads within the lower respiratory tract, how specific immune responses contribute to infection control, and how these relate to disease progression. Although it has been difficult to directly relate these features and obtain high quality data from the lower respiratory tract in humans, we circumvented the challenge by pairing comprehensive experimental data with robust mathematical models and analyses. Our iterative model-driven experimental approach (Smith, 2018a; Smith, 2018b) revealed important dynamic relations between virus, infected cells, immune cells, lung damage, inflammation, and disease severity (summarized in Figure 6). 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.

Summary of the connections between the kinetics of virus, infected cells, CD8+ T cells, lung injury, lung inflammation, and disease severity.

Summary of the relations between the dynamics of virus, infected cells, CD8+ T cells, lung lesions, lung inflammation, and weight loss established by our analysis. Given that viral loads and weight loss are the most easily measured variables, our mathematical model (Equation (1)-(6)) could be used to estimate the kinetics of infected cells, CD8+ T cells, and inflammation. The cumulative area under the curve (CAUC) of the productively infected cell dynamics (I2) yields an estimate of the percent lung infected (active lesion) while the predicted relative CD8E dynamics (E/Emax) yields an estimate of the percent lung resolved (inactive lesion). The total amount of lung involved (% lung infected and % lung resolved) or the inflammation scores can then be used to estimate weight loss through the functions LL(W) and LI(W). These connections could be reversed and weight loss used to estimate viral load kinetics.

Our histomorphometric data and analyses provided a robust quantification of the extent of influenza virus infection within the lung and the dynamics of lung injury and inflammation. The images look strikingly similar to CT scans and postmortem histopathology from infected patients, which also show extensive bilateral damage and cytopathic effects in epithelial cells (Kohr et al., 2010; Li et al., 2011; Ajlan et al., 2009; Soto-Abraham et al., 2009; Gill et al., 2010; Li et al., 2012; Shieh et al., 2010). In addition, our results agree with weekly CT imaging of hospitalized patients infected with the 2009 H1N1 influenza virus showing that the damage worsens during the first week of illness and that cell proliferation begins in the second week (Li et al., 2012). Further investigating the utility of CT imaging could enhance the translatability of our findings. The model’s 95% CI predictions of the lesioned area nicely captured the heterogeneity in the data (Figure 4), which is remarkable given the minimal variability in the virus and CD8+ T cell data and model fit (Figure 1). The differential heterogeneity suggests that small changes to lung viral loads can induce significant changes in disease severity. Indeed, this has been observed in both humans and animals when antivirals are administered (Toniolo Neto, 2014; Baloxavir Marboxil Investigators Group et al., 2018; Treanor et al., 2000; Moscona, 2005; Gubareva et al., 2000; Koshimichi et al., 2018) or when host responses are experimentally modulated (Channappanavar et al., 2016; Szretter et al., 2009; Iwasaki and Nozima, 1977). It was also observed in our CD8 depletion data, where our model predicted a significant reduction in the lesioned area (Figure 4) while viral loads were only slightly lower (Figure 3). This may indicate that viral loads are generally not a reliable measure of infection severity, and the knowledge should aid future studies seeking to estimate therapeutic efficacy or the effects of immune modulation. There was a slight over-estimation of the inactive lesion early in the infection (1–5 d pi, dark gray shading; Figure 4) when using the fit to the total CD8+ T cells, which was abrogated when using the CD25-CD43+CD8+ T cells. However, these cells underestimated the total lesion at 6–7 d pi, suggesting contributions from other CD8+ phenotypes. In general, modeling different CD8+ T cell phenotypes may help clarify their in vivo functions and better define the rates of CD8-mediated infected cell clearance. Although our results suggest that this is not necessarily required, they may be more critical to consider when predicting the response of preexisting and potentially cross-reactive immunity particularly given that memory resident T cells are only beginning to be understood (Zens et al., 2016; Slütter et al., 2017; Van Braeckel-Budimir et al., 2018; Wakim et al., 2015; Uddbäck et al., 2021).

The new knowledge about the extent of infection within the lung also uncovered the nonlinear connection between disease severity and various measures of lung pathology (Figure 5). This discovery is significant because it suggests there is a unique link to the extent of both lung injury and inflammation, which are often used interchangeably yet have distinct dynamics (Figure 4) and correlate to different immune responses (Figure 4—figure supplements 12). In the same vein, both measurements can be estimated using the dynamics of infected cells (Figure 4) but with independent mathematical relations (CAUC of I2 versus Equation (7)). Our CD8 depletion data (Figure 3) and other experimental studies using histomorphometry (Marathe et al., 2016) corroborate a relationship between weight loss and lung pathology during IAV infection. For example, animals treated with antivirals in various conditions (single or combination therapy and in immunocompetent or immunosuppressed hosts) demonstrated that, although viral loads are not always significantly reduced, the antiviral-induced reductions in weight loss were paired with decreased infected areas of the lung (Marathe et al., 2016). Further examining lung injury and inflammation kinetics and their connection to weight loss in various experimental infection settings (e.g., different ages or sex, under therapy, other strains or viruses, etc.) and deciphering how innate and humoral responses contribute to these measurements should improve our predictive capabilities. In addition, the nonlinearity of the connections warrants further investigation as it could indicate cooperativity of multiple mechanisms. Studies such as these may be particularly helpful in understanding the exacerbated morbidity and mortality in elderly individuals, where there are lower viral loads, increased weight loss, symptoms, and/or mortality within animal models (Toapanta and Ross, 2009; Smith et al., 2019). Because murine models are imperfect systems with important immunologic and physiologic distinctions from humans (Mestas and Hughes, 2004), uncovering how our results translate to human infection with novel or circulating influenza viruses and identifying the corresponding, tractable human symptom that could be a robust predictor of disease would be ideal. However, this could be arduous due to inherent human and viral heterogeneity and contributions from co-morbidities and complex immunologic histories.

Despite their limitations, animal models are used to make viral and immunological discoveries and to develop and test therapeutics and vaccines (reviewed in Thangavel and Bouvier, 2014; Margine and Krammer, 2014; Smith and McCullers, 2014), and this work demonstrates a significant potential for the easily obtained weight loss data to be used, analyzed, and interpreted within both modeling and experimental studies. In our data and others’ animal data, there is a spike in weight loss, symptom score, and/or inflammation at ∼6 d pi (Music et al., 2014; Kugel et al., 2009; Lu et al., 2018; Huang et al., 2017; Figure 5). Our data suggests that this may be due to CD8E activity while the infection continues to spread. The increase in animal weights starting at 7 d pi is concurrent with the decline in infected lesions (Figure 5). Interestingly, there was no corresponding decrease in inflammation, and the spike was not apparent in the CD8 depletion data where severity was reduced. The subsequent increase in weight loss following infection resolution (8–9 d pi; Figure 4C) could be attributed to CD8+ T cell-mediated pathology and/or to ongoing inflammation. The tighter correlation to inflammation might suggest the latter; however, some human and animal studies indicate that large numbers of CD8+ T cells pose a risk of acute lung tissue injury (Duan and Thomas, 2016; Moskophidis and Kioussis, 1998; Mauad et al., 2010; Rygiel et al., 2009; La Gruta et al., 2007). According to our model predictions, excessive CD8+ T cell numbers may augment disease progression yet do not improve recovery time (Appendix 3—figure 2). Instead, an earlier onset of CD8E proliferation (i.e., smaller τE) would be required to significantly shorten the infection (Appendix 3—figure 2). This aligns with evidence that hosts with adaptive immune responses primed by vaccine or prior infection recover more rapidly (Duan et al., 2015; Weinfurter et al., 2011). While higher CD8+ T cell numbers have little impact on viral kinetics, our model and data agree with clinical and experimental studies from a wide range of host species that impaired CD8+ T cell responses can prolong an IAV infection (Wang et al., 2015; Miao et al., 2010; McMichael et al., 1983; Yap and Ada, 1978; Wells et al., 1981; Bender et al., 1992). In some scenarios, such as in immunocompromised hosts, virus can persist for up to several weeks if CD8+ T cell-mediated clearance is unsuccessful (Figure 3; Wang et al., 2015; Eichelberger et al., 1991; Hou et al., 1992; Xue et al., 2017; Memoli et al., 2014). The bifurcation in recovery time revealed by the model suggests that this may occur when the number of lung-specific CD8+ T cells are less than 2.31 × 105 cells (Figure 3D–E). Our CD8 depletion data show that the precise number will be dependent on other infection variables. Minimally, we would expect this number to vary depending on parameters like the dose, rate of virus replication, degree of prior immunity, and/or the infected cell number and lifespan, which has been noted in another modeling study that detailed similar bifurcating behavior (Yates et al., 2011). Although some previously published models also suggested delayed resolution with depleted CD8+ T cell responses (Dobrovolny et al., 2013; Miao et al., 2010; Cao et al., 2016; Lee et al., 2009), this bifurcation has not been observed and their estimated delays in recovery do not amount to the long-lasting infections in immunodeficient patients (Wang et al., 2015; Eichelberger et al., 1991; Hou et al., 1992). Our model’s ability to capture the dynamics when CD8+ T cells are depleted is encouraging, and the data are an important reminder that experimental modulation of cell populations (e.g., through depletion or genetic knockouts) is complicated and that the data from such systems should be interpreted cautiously. It also brings into question studies that have used these types of data without validated mathematical models or quantification of other immunological variables.

Our data and analyses provided strong support that the infected cell density impacts the rate at which they are cleared by effector CD8+ T cells (CD8E) (Figures 34, Figure 4—figure supplement 1). Although CD8+ T cell dynamics can be somewhat replicated when assuming the density dependence lies within their expansion, that assumption could not recover the viral load dynamics (Appendix 1—figure 1). Discriminating between these mechanisms is difficult in vivo, but the ability of the model in Equation (1)-(6) to capture the entire time course of CD8E dynamics in multiple scenarios is compelling. Regardless of the mechanism, we first detected this density-dependence in a model that excluded specific immune responses (see Appendix 2) (Smith et al., 2018). Simple models like that one are useful to capture nonlinearities, but they cannot distinguish between different mechanisms.

Several factors may contribute to the density-dependent change in the rate of CD8E-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 Smith et al., 2018; Li and Handel, 2014; Le et al., 2014; Graw and Regoes, 2009; Gadhamsetty et al., 2014; Pilyugin and Antia, 2000; Smith et al., 2011b). When CD8E interact with infected cells, a complex is formed for ∼20–40 min (Mempel et al., 2006; Deguine et al., 2010; Zagury et al., 1975; Yannelli et al., 1986; Perelson and Bell, 1982). Because CD8E could not interact with other infected cells during this time, the global rate of infected cell clearance would be lowest when infected cells outnumber CD8E. In addition, contact by a single CD8E may be insufficient to remove an infected cell (Halle et al., 2016). Infected cell clearance is more frequently observed after interactions with multiple CD8E, with an average of 3.9 contact events either serially or simultaneously (Halle et al., 2016). 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 of times to induce cell death. However, as CD8E accumulate and the density of infected cells decreases (Figure 3A), the probability of simultaneous interactions will increase. These interactions may also be influenced by CD8E movement, where their velocity slows at the height of infected cell clearance (6–8 d pi) (Lambert Emo et al., 2016; Matheu et al., 2013). This should reduce the time required to remove an infected cell and, thus, result in a higher efficiency. Moreover, 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., in tuberculosis granulomas [Ulrichs et al., 2004; Russell et al., 2009]) and has been suggested in agent-based models of influenza virus infection (Levin et al., 2016).

In addition to illuminating the connections between viral kinetics and pathology, the histomorphometric data validated the model’s infected cell dynamics (Figure 4B). The dynamics of susceptible and infected cells throughout the infection and the accuracy of the target cell limited approximation used within influenza viral kinetic models have been questioned for several years (Smith and Perelson, 2011; Smith, 2018b; Beauchemin and Handel, 2011; Ahmed et al., 2017; Gallagher et al., 2018; Boianelli et al., 2015; Handel et al., 2018). The ability of the model to accurately predict the histomorphometry and CD8+ T cell depletion data corroborates the use of this approximation, which assumes a limited number of available target cells and describes their decline by infection only. Although the slowing of viral loads beginning at 2 d pi could be due to a variety of innate immune-mediated mechanisms (e.g., macrophages, neutrophils, type I interferons, etc.), adding more complex dynamics to our model was unnecessary to describe the data. Further, the model and data agree that there are few infected cells during the time when viral loads are growing most rapidly (0–2 d pi; Figures 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 (Smith et al., 2010), which has helped numerous studies interpret their results (Miao et al., 2010; Smith et al., 2011a; Smith, 2018b; Holder and Beauchemin, 2011). The data also supports the model’s hypothesis that that there is minimal clearance of infected cells prior to CD8E expansion (Figure 4). The knowledge of the model’s accuracy and of the spatiality in 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 theoretical predictions like the ones presented here is pivotal to elucidating the mechanisms of infection spread and clearance (Smith, 2018b). Examining other infections (e.g., coronaviruses), modifications to the dynamics (e.g., lethal doses), and the connection between lung measurements and those more easily acquired from upper respiratory tract will help refine the dynamical links between virus, host immune responses, and disease severity and identify their generalizability. Determining the factors that influence disease severity is vital to understanding the disproportionate mortality in at-risk populations (e.g., elderly) and to improving therapeutic design. This is particularly important because current antivirals alleviate symptoms but do not always effectively lower viral loads (Toniolo Neto, 2014; Baloxavir Marboxil Investigators Group et al., 2018; Treanor et al., 2000; Moscona, 2005; Gubareva et al., 2000; Koshimichi et al., 2018). The predictive capabilities of validated models like the one here should prove useful in forecasting infection dynamics for a variety of scenarios. These tools and analyses provide a more meaningful interpretation of infection data, new ways to interpret weight loss data in animal models, and a deeper understanding of the progression and resolution of the disease, which will undoubtedly aid our ability to effectively combat influenza.

Materials and methods

Ethics statement

Request a detailed protocol

All experimental procedures were performed under protocols O2A-020 or 17–096 approved by the Animal Care and Use Committees at St. Jude Children’s Research Hospital (SJCRH) or the University of Tennessee Health Science Center (UTHSC), respectively, under relevant institutional and American Veterinary Medical Association (AVMA) guidelines. All experimental procedures were performed in a biosafety level two facility that is accredited by the American Association for Laboratory Animal Science (AALAS).

Mice

Adult (6-week-old) female BALB/cJ mice were obtained from Jackson Laboratories (Bar Harbor, ME) or Charles River Laboratories (Wilmington, Massachusetts). Mice were housed in groups of five mice in high-temperature 31.2 cm × 23.5 cm × 15.2 cm polycarbonate cages with isolator lids (SJCRH) or in 38.2 cm × 19.4 cm × 13.0 cm solid-bottom polysulfone individually ventilated cages (UTHSC). Rooms used for housing mice were maintained on a 12:12 hr light:dark cycle at 22 ± 2°C with 50% humidity in the biosafety level two facility at SJCRH (Memphis, TN) or UTHSC (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; SJCRH) or Teklad LM-485 Mouse/Rat Sterilizable Diet (Envigo, Indianapolis, IN; UTHSC) 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 SJCRH or UTHSC.

Infectious agents

Request a detailed protocol

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

Infection experiments

Request a detailed protocol

The viral infectious dose (TCID50) was determined by interpolation using the method of Reed and Muench, 1938 using serial dilutions of virus on Madin-Darby canine kidney (MDCK) cells. Mice were intranasally inoculated with 75 TCID50 PR8 diluted in 100 μl of sterile PBS. In a subset of animals, CD8+ T cells were depleted by IP injection at days −2, 0, 3, and 7 pi with 100 μg of the rat anti-CD8α antibody (clone 2.43) that was purified from ATCC hybridoma (per manufacturer instructions) in 250 μl of PBS. Depletion efficiency was confirmed in the lung and spleen as described below. Mice were weighed at the onset of infection and each subsequent day to monitor illness and mortality. Mice were euthanized if they became moribund or lost 30% of their starting body weight. For viral load and total CD8+ T cell quantification, experiments were repeated three times and in each facility to ensure reproducibility and two complete experiments (ten animals per time point) were used for these studies. For all other experiments, the data was compared to prior results in addition to being repeated at select time points to ensure reproducibility, and one complete experiment (five animals per time point) was used. For pathology scoring and histomorphometry quantitation, five animals per time point were used. Power was calculated using G*Power.

Lung harvesting for viral and cellular dynamics

Request a detailed protocol

For total CD8+ T cell quantification, mice were euthanized by CO2 asphyxiation. To distinguish blood-borne CD8+ T cells from those in the lung parenchyma, deeply anesthetized mice (5% isoflurane) were retro-orbitally injected with 3 μg of anti-CD45 antibody (PerCP, clone 30-F11, Biolegend) 3 min prior to euthanasia (Anderson et al., 2012; Keating et al., 2018). Mice were then euthanized by 33% isoflurane inhalation and their lungs perfused with 10 ml PBS prior to removal. For all experiments, lungs were then 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, 500xg for 7 min. The supernatants were used to determine the viral titers (TCID50) by serial dilutions on MDCK monolayers. Following red blood cell lysis, cells were washed in MACS buffer (PBS, 0.1M EDTA, 0.01M HEPES, 5 mM EDTA and 5% heat-inactivated FBS). Cells were then counted with trypan blue exclusion using a Cell Countess System (Invitrogen, Grand Island, NY) and prepared for flow cytometric analysis as indicated below.

Lung titers

Request a detailed protocol

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 Figure 1B, provided in Source data 1, and were previously published and utilized for calibration of the density-dependent model (see Appendix 2; Smith et al., 2018).

Flow cytometric analysis

Request a detailed protocol

Flow cytometry (LSRII, BD Biosciences, San Jose, CA (SJCRH) or ZE5 Cell Analyzer, Bio-Rad, Hercules, CA (UTHSC)) was performed on the cell pellets after incubation with 200 μl of a 1:2 dilution of Fc block (human-γ globulin) on ice for 30 min, followed by viability (Biolegend, Zombie Violet Fixable Viability Kit) and surface marker staining with anti-mouse antibodies. For total CD8+ T cell, macrophage, and neutrophil quantification, we used antibodies CD3e (BV786, clone 145–2 C11, Biolegend), CD4 (PE-Cy5, clone RM4-5, BD Biosciences), CD8α (BV605, clone 53–6.7, BD Biosciences), Ly6C (APC, clone HK1.4, eBioscience), F4/80 (PE, clone BM8, eBioscience), CD11c (eFluor450, clone N418, eBioscience), CD11b (Alexa700, clone M1/70, BD Biosciences), MHC-II (FITC, clone M5/114.15.2, Invitrogen), CD49b (APCe780, clone DX5, eBioscience), and Ly6G (PerCp-Cy5.5, clone 1A8, Biolegend). The data were analyzed using FlowJo 10.4.2 (Tree Star, Ashland, OR) where viable cells were gated from a forward scatter/side scatter plot and singlet inclusion. Following neutrophil (Ly6Ghi) and subsequent macrophage (CD11chiF4/80hi) exclusion, CD8+ T cells were gated as CD3e+DX5CD4CD8+.

To distinguish blood borne CD8+ T cells from those in the lung parenchyma, we used antibodies CD3e (BV786, clone 145–2 C11, Biolegend), CD4 (FITC, clone RM4-5, Biolegend), CD8α (BV605, clone 53–6.7, BD Biosciences), B220 (APCe780, clone RA3-6B2, eBioscience), CD49b (APCe780, clone DX5, eBioscience), CD62L (PE-Cy7, clone MEL-14, Biolegend), CD69 (PE, clone H1.2F3, Biolegend), CD44 (PE-Dazzle594, clone IM7, Biolegend), CD25 (Alexa700, clone PC61, Biolegend), CD43 (APC, clone 1B11, Biolegend), CD127 (PE-Cy5, clone A7R34, Biolegend), and CD103 (BV711, clone 2E7, Biolegend) or CD314 (NKG2D) (BV711, clone CX5, BD Bioscience). The data were analyzed using FlowJo 10.6.2 (Tree Star, Ashland, OR) where viable cells were gated from a forward scatter/side scatter plot, singlet inclusion, and viability dye exclusion. Blood borne CD8+ T cells were gated as CD45+ and those in the lung parenchyma as CD45-. The total in each sub-population was gated as CD3e+B220DX5CD4CD8+. IAV-specific CD8+ T cells were then gated as CD25+CD43+ (recently activated), CD25CD43+ (effector) (Keating et al., 2018), and CD43CD127+CD103+CD69+ (long-lived memory). Expression of CD44, CD69, CD62L, and NKG2D were also assessed to ensure appropriate classification.

For all experiments, the absolute numbers of CD8+ T cells were calculated based on viable events analyzed by flow cytometry as related to the total number of viable cells per sample. The data are provided in Source data 1. We use the kinetics of the total number of CD8+ T cells (non-perfused lung; ‘total’), the total number of CD8+ T cells in the lung parenchyma (perfused lung; ‘lung-specific’), and virus-specific CD8+ T cell phenotypes in the lung parenchyma as defined by surface staining (described above). We chose this approach because the use of tetramers yields epitope-specific cell dynamics that vary in time and magnitude (e.g., as in Toapanta and Ross, 2009; Keating et al., 2018) and would complicate the model and potentially skew the results. In addition, tetramer staining is not available for all epitopes in BALB/cJ mice and they may underestimate the dynamics of IAV-specific cells (Keating et al., 2018). The gating strategies are shown in Figure 1—figure supplement 2.

Lung immunohistopathologic and immunohistochemical (IHC) evaluation

Request a detailed protocol

The lungs from IAV-infected mice were fixed via intratracheal infusion and then immersion in 10% buffered formalin solution. Tissues were paraffin embedded, sectioned, and stained for influenza virus using a primary goat polyclonal antibody (US Biological, Swampscott, MA) against influenza A, USSR (H1N1) at 1:1000 and a secondary biotinylated donkey anti-goat antibody (sc-2042; Santa Cruz Biotechnology, Santa Cruz, CA) at 1:200 on tissue sections subjected to antigen retrieval for 30 min at 98°C. The extent of virus spread was quantified by capturing digital images of whole-lung sections (two dimensional) stained for viral antigen using an Aperio ScanScope XT Slide Scanner (Aperio Technologies, Vista, CA) then manually outlining defined fields. Alveolar areas containing virus antigen-positive pneumocytes were highlighted in red (defined as ‘active’ infection), whereas lesioned areas containing minimal or no virus antigen-positive debris were highlighted in green (defined as ‘inactive’ infection). Lesions containing a mix of virus antigen-positive and antigen-negative pneumocytes were highlighted in orange (defined as ‘mixed’ infection). The percentage of each defined lung field was calculated using the Aperio ImageScope software. Pulmonary lesions in HE-stained histologic sections were assigned scores based on their severity and extent. Representative images from each group and quantitative analyses (five animals/group) of viral spread and lung pathology during IAV infection are shown in Figure 4, and provided in Source data 1.

CD8+ T cell model

Request a detailed protocol

To examine the contribution of CD8+ T cells to the biphasic viral load decay, we expanded the density-dependent (DD) model (see Appendix 2) to include two mechanisms of infected cell clearance (non-specific clearance (δ) and CD8+ T cell-mediated clearance (δE(I2,E))) and two CD8+ T cell populations: effector (E, denoted CD8E) and memory (EM, denoted CD8M) CD8+ T cells. The model is given by Equation (1)-(6).

(1) dTdt=-βTV
(2) dI1dt=βTV-kI1
(3) dI2dt=kI1-δI2-δE(I2,E)I2
(4) dVdt=pI2-cV
(5) dEdt=ξ(E)I2+ηEI2(t-τE)-dEE
(6) dEMdt=ζE(t-τM)

In this model, target cells become infected with virus at rate βV per day. Once infected, cells enter an eclipse phase (I1) before transitioning at rate k per day to a productively-infected state (I2). These infected cells produce virus at rate p TCID50 per infected cell per day, and virus is cleared at rate c per day. Virus-producing infected cells (I2) are cleared by non-specific mechanisms (e.g., apoptosis and/or innate immune responses) at a constant rate δ per day. Innate immune responses were excluded from the model because the viral load data is linear during the time where they act (2–7 d pi) and, thus, additional equations cannot improve the fit. Productively infected cells are cleared by CD8E at rate δE(I2,E)=δEE/(KδE+I2) per day, where the rate of infected cell clearance is δE/KδE per CD8E per day and KδE is the half-saturation constant. The CD8E-mediated clearance rate (δE(I2,E)) is dependent on the density of infected cells and is similar to the infected cell clearance term in the DD model (see Equation (A11); Smith et al., 2018). Similar density-dependent forms have also been used in models that describe the CD8+ T cell response to other virus infections (De Boer and Perelson, 1998; Li and Handel, 2014; Gadhamsetty et al., 2014). Models that exclude this density-dependence were examined, but these models resulted in a statistically poor fit to the data as defined by the Akaike Information Criteria (AIC) (Supplementary file 2). This is due in part to the increase of CD8+ T cells from 3-5 d pi (Figure 1). We did examine a model that excluded these dynamics and included CD8+ T cell expansion as a density-dependent function (e.g., ηEI2(t-τE)/(KI+I2)) while keeping a linear rate of CD8-mediated infected cell clearance (δEEI2) (see Appendix 1). This model was adapted from other published models (Conway and Perelson, 2015; Bonhoeffer et al., 2000; Baral et al., 2019) and can produce similar dynamics from 6-10 d pi, but it was not statistically supported by our data (Supplementary file 1). We further tested the model in Baral et al., 2019, but this model could not fit our data and lacked statistical support (Supplementary file 1).

The model assumes that the initial CD8E influx in the lung is proportional to infected cells at rate ξ(E)=ξ/(KE+E) CD8E per cell per day, which is down-regulated by the CD8E already present in the lung. The associated half-saturation constant is KE. Similar terms for CD8E regulation have been used in modeling HIV infections (De Boer and Perelson, 1998; Müller et al., 2001) and in models that examined CD8+ T cell proliferation mechanisms (De Boer and Perelson, 1995). We also examined whether their influx is proportional to infected cells at a delay (i.e., ξ(E)I(t-τI)). While this modification better captured the initial increase in CD8+ T cells, the additional parameter was not supported. In our model, CD8E expansion occurs at rate η per infected cell per day with time delay τE. This term accounts for local CD8E proliferation in the lung (McGill and Legge, 2009; Lawrence et al., 2005) and migration of CD8E from secondary lymphoid organs (Zhang and Bevan, 2011; Bedoui and Gebhardt, 2011; Wu et al., 2011; Kang et al., 2011). The delay may signify the time it takes CD8E to become activated by antigen presenting cells, differentiate, proliferate, and/or migrate to the infection site. The lung CD8E population declines due to cell death and/or emigration at rate dE per day. These cells transition to CD8M (EM) at rate ζ CD8M per CD8E per day after τM days. The model schematic and fit to the viral load, total, effector, and memory CD8+ T cell data are in Figure 1. All code is provided in Source code 1.

Inflammation model

Request a detailed protocol

To estimate the alveolar and interstitial inflammation without modeling other cell classes (e.g., macrophages and neutrophils), we assumed that inflammation in the lung (LI) was proportional to the infected cells (I1 and I2) according to the equation,

(7) dLIdt=α1I1+α2I2,

where α1,2 has units of score/cell/d and defines the inflammation score contribution from each infected cell class. A decay term was excluded because inflammation does not resolve on the timescale of our data.

Parameter estimation

Request a detailed protocol

Given a parameter set θ, the cost C(θ) was minimized across parameter ranges using an Adaptive Simulated Annealing (ASA) global optimization algorithm (Smith et al., 2018) to compare experimental and predicted values of virus (V; log10 TCID50/lung) and log10 total CD8+ T cells/lung (E^=E+EM+E^0, where E^0 is the initial number of CD8+ T cells at 0 d pi), or log10 TCID50/lung virus (V), log10 effector CD8+ T cells/lung (E^=E+E^0), and log10 memory CD8+ T cells/lung (E^M=EM+E^M0). The cost function is defined by,

C(θ)=i,j(V(θ,ti)vi,j)2+i,j(E^M(θ,ti)mi,j)2+sE[i,j(E^(θ,ti)ei,j)2+iγi(E^(θ,ti+1)E^(θ,ti1)ti+1ti11γijei+1,jei1,jti+1ti1)2],

where (ti,vi,j) is the viral load data, (ti,ei,j) is the total or effector CD8+ T cell data, (ti,mi,j) is the memory CD8+ T cell data, and V(θ,ti), E^(θ,ti), and EM^(θ,ti) are the corresponding model predictions. Here, sE=(vmax-vmin)/(emax-emin) is a scaling factor, and γi=Ji+1Ji-1 where Ji is the number of observations at time ti. Errors of the log10 data were assumed to be normally distributed. To explore and visualize the regions of parameters consistent with the model, we fit Equation (1)-(6) to 2000 bootstrap replicates of the data. If the fit was within χ2=0.05 of the best-fit and the CD8 derivative was not a statistical outlier as determined by the function isoutlier, then the bootstrap was considered successful (Smith et al., 2011a; Smith et al., 2013; Smith et al., 2018). For each best-fit estimate, we provide 95% confidence intervals (CI) obtained from the bootstrap replicates (Table 1). Calculations were performed either in MATLAB using a custom built ASA algorithm (Smith et al., 2018) or in Python using the simanneal package (Perry, 2018) followed by a L-BFGS-B (Byrd et al., 1995; Zhu et al., 1997) deterministic minimization through SciPy’s minimize function. MATLAB ode15s and dde23 or SciPy integrate.ode using lsoda and PyDDE (Cairns, 2008) were used as the ODE and DDE solvers.

Estimated parameters in the CD8+ T cell model included the rates of virus infection (β), virus production (p), virus clearance (c), eclipse phase transition (k), non-specific infected cell clearance (δ), CD8E-mediated infected cell clearance (δE), half-saturation constants (KδE and KE), CD8E infiltration (ξ), CD8E expansion (η), delay in CD8E expansion (τE), CD8clearance (dE), CD8M generation (ζ), delay in CD8M generation (τM), and the baseline number of CD8+ T cells (E^0, M^0). 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., 2018). The rate of infection (β) was allowed to vary between 10-6-10-1TCID50-1d-1, and the rate of virus production (p) between 10-1-103TCID50cell-1d-1. Bounds for the virus clearance rate (c) were 1d-1 (t1/2=16.7h) and 103d-1 (t1/2=1min). To ensure biological feasibility, the lower and upper bounds for the eclipse phase transition rate (k) were 4-6 d-1 as done previously (Smith et al., 2018).

The rate of non-specific infected cell clearance (δ) was given limits of 0.05–1 d-1. The CD8E-mediated infected cell clearance rate (δE) varied between 0.01-2cellsCD8E-1d-1, and the associated half-saturation constant (KδE) was bounded between 101–106 cells. The upper bound of δE was chosen to maintain the convergence of δ to nonzero values and consistency with prior results where δ approximates the slope of the first decay phase (Smith et al., 2010). Bounds for the rate of CD8E infiltration (ξ) were 102106CD8E2 cell1 d1, and bounds for the half-saturation constant (KE) were 103-107 CD8E. The CD8E expansion rate (η) varied between 10–8–10–6 cell–1 d–1, and the delay in CD8E expansion (τE) between 2–6 d. The rate of CD8E clearance (dE) had limits of 0.05–2 d–1. The rate of CD8M generation (ζ) varied between 0.011 CD8M CD8E1 d1, and the delay in CD8M generation (τM) varied between 3–4 d (total CD8+ T cell fit) or 2–6 d pi (effector and memory CD8+ T cell fit). Larger bounds were examined for this parameter; however, the parameter was non-identifiable in the total CD8+ T cell fit and a small range was required for convergence. Bounds for the baseline number of CD8+ T cells (E^0, E^M0) were set to the upper and lower values of the data at 0 d pi (3.0×105-5.3×105CD8 (total), 2.5×102-1.4×103CD8E (effector), 3.4×102-8.0×102CD8M (memory)).

The initial number of target cells (T(0)) was set to 107 cells (Smith et al., 2011a). The initial number of infected cells I1(0) was set to 75 cells to reflect an initial dose of 75 TCID50 (Smith et al., 2018). We previously found that estimating I1(0), fixing V(0)=75TCID50, or estimating V(0) did not improve the fit and could not be statistically justified (Smith et al., 2018). The initial number of productively infected cells (I2(0)), the initial free virus (V(0)), and the initial number of CD8E (E(0)) and CD8M (EM(0)) were set to 0. All other parameter estimations were done as described in the text.

Linear regression

Request a detailed protocol

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 4—figure supplement 1.

Cumulative area under the curve

Request a detailed protocol

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

Appendix 1

Alternate CD8+ T cell Models

We examined alternate formulations of the CD8+ T cell model to further investigate the density-dependence in the CD8+ T cell response. Rather than assuming that CD8E-mediated clearance of infected cells is dependent on their density, the model in Equation (A1)-(A6) assumes that the rate of CD8E expansion is dependent on the density of infected cells. Similar models have been used to study CD8+ T cell responses during HIV infection (Conway and Perelson, 2015; Bonhoeffer et al., 2000) and other viral infections (Baral et al., 2019). The differences between this alternate model and the CD8+ T cell model (Equation (1)-(6)) are in bold.

(A1) dTdt=-βTV
(A2) dI1dt=βTV-kI1
(A3) dI2dt=kI1-δI2-δ𝐄𝐚𝐄𝐈𝟐
(A4) dVdt=pI2-cV
(A5) dEdt=η𝐄𝐚𝐊𝐄𝐚+𝐈𝟐I2(t-τE)-dEE
(A6) dEMdt=ζE(t-τM)

When a linear CD8E-mediated infected cell clearance rate is included, the CD8E dynamics between 3–5 d pi cannot be replicated and returned a statistically worse fit via AIC (Supplementary file 1). However, because these cells may not have effector functions and contribute to infected cell clearance at this time (see Discussion), we excluded these data and the term ξI2/(KE+E) when fitting Equation (A1)-(A6) to the viral load and total CD8+ T cell data (Appendix 1—figure 1). The CD8 dynamics are similar to those generated by CD8+ T cell model. However, the alternate model underestimates the data at day seven and the sharp decline between 7–8 d. While it is inappropriate to directly compare these models due to the varying number of data points, we assessed their goodness of fit with and without contributions from the data at 3–5 d pi (Supplementary file 1). Regardless of the data inclusion or exclusion, the alternate model had a higher AIC in all contexts and, thus, is not statistically justifiable. Similarly, we examined the model in Baral et al., 2019 (Equation (A7)-(A8)), which also could not replicate our data (Appendix 1—figure 1).

(A7) dIdt=κbI(1-IImax)-δEbIE
(A8) dEdt=ηbIEKEb+IdEbIEKDb+I
Appendix 1—figure 1
Fit of alternate CD8+ T cell models.

Fit of the CD8+ T cell model (Equation (1)-(6); solid black line) compared to the fit of two alternate CD8+ T cell models (Equation (A1)-(A6) (solid blue line) and Equation (A7)-(A8) (dashed blue line)) to virus and CD8+ T cells (excluding 3–5 d pi; white squares) from the lungs of mice infected with 75 TCID50 PR8 (10 mice per time point). Resulting parameter values were δEa=4.02×106 CD8E1d1, ηEa=3.12×105CD8E/d, and KEa=9.53×105 infected cells, and κb=3.33d-1, Imax=1.0 ×106 infected cells, δEb=1.53×10-5CD8E-1d-1, ηb=1.49d-1, KEb=9.5 infected cells, dEb=0.95d-1, and KDb=1.90 ×106 infected cells. All other parameters are in Table 1, and the AICs are in Supplementary file 1. Data are shown as mean ± standard deviation.

Appendix 2

Comparison of the Density-Dependent and CD8+ T cell Models

We previously developed a density-dependent (DD) viral kinetic model, which describes the biphasic decline of viral loads without inclusion of specific host responses (Smith et al., 2018). This model tracks four populations: susceptible epithelial (‘target’) cells (T), two classes of infected cells (I1 and I2), and virus (V) (Smith et al., 2018).

(A9) dTdt=-βTV
(A10) dI1dt=βTV-kI1
(A11) dI2dt=kI1-δd(I2)I2
(A12) dVdt=pI2-cV

Briefly, in the DD model, virus-producing infected cells (I2) are cleared according to the function δd(I2)=δd/(Kδ+I2), where δd/Kδ is the maximum per day rate of infected cell clearance and Kδ is the half-saturation constant (Appendix 2—figure 1). All other terms are common to the CD8 T+ cell model (Equation (1)-(6)). The DD model provides a close fit to the viral load data in Figure 1B and replicates the biphasic viral load decline while excluding the dynamics of specific immune responses (Smith et al., 2018). Unsurprisingly, the CD8+ T cell model is also capable of reproducing the biphasic viral load decay (Figure 1B and Appendix 2—figure 1). In that model, infected cell clearance is split into terms for non-specific clearance (δ) and CD8E-mediated clearance (δE(I2,E)=δEE/(Kδ+I2)) (Appendix 2—figure 1).

Because the CD8+ T cell model is more mechanistic than the DD model, most of the correlations between the parameters common to both models (i.e., the rates of virus infectivity (β), virus production (p), and virus clearance (c)) were reduced (Appendix 2—figure 1A). In addition, the correlations between the infected cell clearance parameters (δd and Kδ or δE and KδE) and between the rate of virus infectivity (β) and their ratios (δd/Kδ or δE/KδE) were abolished (Figure 2—figure supplement 1–2). There was a negative correlation between the infected cell clearance parameters (δ and δE; Figure 2B), which may reflect the connection between the efficacy of early immune mechanisms and the CD8+ T cell response. This result is in line with experimental evidence that the innate immune responses modulate the activation of adaptive immunity (Iwasaki and Medzhitov, 2015; Luster, 2002; Iwasaki and Medzhitov, 2010; Le Bon and Tough, 2002; Jain and Pasare, 2017).

The differences in model structure between the two models yielded changes in parameter sensitivity and model behavior during the rapid viral clearance phase (Appendix 2—figure 1). In the DD model, the most sensitive parameter is the infected cell clearance, δd (Appendix 2—figure 1). A 50% decrease in this parameter resulted in a ∼7 d delay in viral resolution (Appendix 2—figure 1; Smith et al., 2018). In the CD8+ T cell model, however, viral resolution is delayed by <1 d if the CD8E-mediated infected cell clearance parameter (δE) is reduced by 50% (Appendix 2—figure 1, Appendix 3—figure 1). The rates of CD8E expansion (η) and decay (dE) are sensitive and, thus, significantly influence the viral resolution kinetics (Figures 12). A 50% decrease in η results in a ∼6 d delay in recovery (Appendix 2—figure 1, Appendix 3—figure 1) whereas a 48% decrease in η prolongs the infection by ∼30 d (Figure 3D–E). This bifurcation in recovery time is a unique feature of the CD8+ T cell model.

Appendix 2—figure 1
Parameter behavior of the density-dependent model and the CD8+ T cell model.

(A) Comparison of parameters that were common between the density-dependent model (gray, Equation (A9)-(A12)) and the CD8+ T cell model (purple, Equation (1)-(6)). Correlations were 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) was not well-defined in either model. (B) In the density-dependent model (gray), the viral kinetics and the infection duration were sensitive to small changes in the infected cell clearance parameter (δd). This parameter was well-defined with a narrow 95% CI. In the CD8+ T cell model (purple), changing the CD8E-mediated infected cell clearance parameter (δE) had little impact on viral kinetics or CD8+ T cell kinetics. However, these kinetics were most sensitive to changes in the rate of CD8E expansion (η), which was well-defined with a narrow 95% CI.

Appendix 3

Regulation of the CD8+ T cell Response

To further understand the regulation of the CD8+ T cell response, we examined the 2-D parameter ensembles (Figure 2A–C, Figure 2—figure supplements 12) and the results from the sensitivity analysis (Appendix 3—figure 12). Overall, few parameters were correlated. There was an expected, although small, positive correlation between the rate of CD8E infiltration (ξ) and the associated half-saturation constant (KE) (Figure 2—figure supplement 12), which represents the coordination between CD8E recruitment and the processes that prevent an overabundance of these cells. Likewise, a negative correlation was detected between the rate of initial CD8E influx (ξ) and the initial number of CD8+ T cells (E^0) (Figure 2—figure supplements 12). The influx rate (ξ) was also positively correlated with the delay in CD8E expansion (τE) (Figure 2—figure supplements 12). The rates of CD8E expansion (η) and decay (dE) are correlated (Figure 2C), indicating a balance between these two processes. This correlation was expected and reflects the coordination of mechanisms that regulate CD8+ T cell numbers, which may be necessary to limit excessive immunopathology while still resolving the infection (Moskophidis and Kioussis, 1998; Duan and Thomas, 2016; La Gruta et al., 2007). Further, because of this correlation and the sensitivity of η (Appendix 3—figure 1), the CD8+ T cell kinetics are sensitive to changes in dE (Appendix 3—figure 1). However, increasing the decay rate had less impact on the viral load kinetics, comparatively. Because dE is correlated with both η and the rate of CD8M generation (ζ) (Figure 2C), it naturally follows that η and ζ are correlated (Figure 2—figure supplements 12). Changing the rates of virus infectivity (β), production (p), or clearance (c) had little effect on viral load or CD8+ T cell kinetics (Appendix 3—figure 1).

Appendix 3—figure 1
Sensitivity of the CD8+ T cell model.

Solutions of the CD8+ T cell model (Equation (1)-(6)) with the indicated parameter (β, k, p, c, δ, δE, KδE, ξ, KE, or η) increased (red) or decreased (blue) 50% from the best-fit value (Table 1). CD8E denotes effector CD8+ T cells.

Appendix 3—figure 2
Sensitivity of the CD8+ T cell model.

Solutions of the CD8+ T cell model (Equation (1)-(6)) with the indicated parameter (τE, dE, ζ, τM, or E^0) increased (red) or decreased (blue) 50% from the best-fit value (Table 1). CD8E and CD8M denote effector and memory CD8+ T cells, respectively.

Data availability

All data generated or analyzed during this study are included in the manuscript and supporting files. Source data files have been provided.

References

  1. Software
    1. Cairns BJ
    (2008)
    Pydde, version 0.2.2
    Pydde.
    1. Hou S
    2. Doherty PC
    3. Zijlstra M
    4. Jaenisch R
    5. Katz JM
    (1992)
    Delayed clearance of Sendai virus in mice lacking class I MHC-restricted CD8 T cells
    Journal of Immunology 149:1319–1325.
    1. Iwasaki T
    2. Nozima T
    (1977)
    Defense mechanisms against primary influenza virus infection in mice: I. the roles of interferon and neutralizing antibodies and Thymus dependence of interferon and antibody production
    Journal of Immunology 118:256–263.
    1. Jones AT
    2. Federsppiel B
    3. Ellies LG
    4. Williams MJ
    5. Burgener R
    6. Duronio V
    7. Smith CA
    8. Takei F
    9. Ziltener HJ
    (1994)
    Characterization of the activation-associated isoform of CD43 on murine T lymphocytes
    Journal of Immunology 153:3426–3439.
    1. Perelson AS
    2. Bell GI
    (1982)
    Delivery of lethal hits by cytotoxic T lymphocytes in multicellular conjugates occurs sequentially but at random times
    Journal of Immunology 129:2796–2801.
  2. Software
    1. Perry M
    (2018)
    Simanneal, version 0.4.1
    Simanneal.
    1. Trammell RA
    2. Toth LA
    (2011)
    Markers for predicting death as an outcome for mice used in infectious disease research
    Comparative Medicine 61:492–498.
    1. Wells MA
    2. Albrecht P
    3. Ennis FA
    (1981)
    Recovery from a viral respiratory infection. I. Influenza pneumonia in normal and T-deficient mice
    Journal of Immunology 126:1036–1041.
    1. Yannelli JR
    2. Sullivan JA
    3. Mandell GL
    4. Engelhard VH
    (1986)
    Reorientation and fusion of cytotoxic T lymphocyte granules after interaction with target cells as determined by high resolution cinemicrography
    Journal of Immunology 136:377–382.

Article and author information

Author details

  1. Margaret A Myers

    Department of Pediatrics, University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Software, Formal analysis, Investigation, Visualization, Writing - original draft
    Contributed equally with
    Amanda P Smith
    Competing interests
    No competing interests declared
  2. Amanda P Smith

    Department of Pediatrics, University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Margaret A Myers
    Competing interests
    No competing interests declared
  3. Lindey C Lane

    Department of Pediatrics, University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  4. David J Moquin

    Department of Anesthesiology, Washington University School of Medicine, St. Louis, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  5. Rosemary Aogo

    Department of Pediatrics, University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Formal analysis, Writing - review and editing
    Competing interests
    No competing interests declared
  6. Stacie Woolard

    Flow Cytometry Core, St. Jude Children's Research Hospital, Memphis, United States
    Contribution
    Formal analysis
    Competing interests
    No competing interests declared
  7. Paul Thomas

    Department of Immunology, St. Jude Children’s Research Hospital, Memphis, United States
    Contribution
    Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  8. Peter Vogel

    Department of Pathology, St. Jude Children's Research Hospital, Memphis, United States
    Contribution
    Data curation, Formal analysis, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
  9. Amber M Smith

    Department of Pediatrics, University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    amber.smith@uthsc.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7092-6904

Funding

National Institute of Allergy and Infectious Diseases (AI139088)

  • Margaret A Myers
  • Amanda P Smith
  • Lindey C Lane
  • Rosemary Aogo

National Institute of Allergy and Infectious Diseases (AI125324)

  • Margaret A Myers
  • Amanda P Smith
  • Lindey C Lane
  • David J Moquin
  • Amber M Smith

National Institute of Allergy and Infectious Diseases (AI100946)

  • Amber M Smith

American Lebanese Syrian Associated Charities (Internal Funding)

  • Margaret A Myers
  • Amanda P Smith
  • Lindey C Lane
  • David J Moquin
  • Stacie Woolard
  • Paul Thomas
  • Peter Vogel
  • Amber M Smith

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by NIH grants AI100946, AI125324, and AI139088 and ALSAC. We thank Alan Perelson for his helpful comments, and Robert Michael for technical assistance.

Ethics

Animal experimentation: All experimental procedures were performed under protocols O2A-020 or 17-096 approved by the Animal Care and Use Committees at St. Jude Children's Research Hospital (SJCRH) or the University of Tennessee Health Science Center (UTHSC), respectively, under relevant institutional and American Veterinary Medical Association (AVMA) guidelines. All experimental procedures were performed in a biosafety level 2 facility that is accredited by the American Association for Laboratory Animal Science (AALAS).

Version history

  1. Preprint posted: February 19, 2019 (view preprint)
  2. Received: March 28, 2021
  3. Accepted: July 14, 2021
  4. Accepted Manuscript published: July 20, 2021 (version 1)
  5. Version of Record published: August 17, 2021 (version 2)

Copyright

© 2021, Myers et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,717
    views
  • 391
    downloads
  • 32
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Margaret A Myers
  2. Amanda P Smith
  3. Lindey C Lane
  4. David J Moquin
  5. Rosemary Aogo
  6. Stacie Woolard
  7. Paul Thomas
  8. Peter Vogel
  9. Amber M Smith
(2021)
Dynamically linking influenza virus infection kinetics, lung injury, inflammation, and disease severity
eLife 10:e68864.
https://doi.org/10.7554/eLife.68864

Share this article

https://doi.org/10.7554/eLife.68864

Further reading

    1. Computational and Systems Biology
    2. Genetics and Genomics
    Weichen Song, Yongyong Shi, Guan ning Lin
    Tools and Resources

    We propose a new framework for human genetic association studies: at each locus, a deep learning model (in this study, Sei) is used to calculate the functional genomic activity score for two haplotypes per individual. This score, defined as the Haplotype Function Score (HFS), replaces the original genotype in association studies. Applying the HFS framework to 14 complex traits in the UK Biobank, we identified 3619 independent HFS–trait associations with a significance of p < 5 × 10−8. Fine-mapping revealed 2699 causal associations, corresponding to a median increase of 63 causal findings per trait compared with single-nucleotide polymorphism (SNP)-based analysis. HFS-based enrichment analysis uncovered 727 pathway–trait associations and 153 tissue–trait associations with strong biological interpretability, including ‘circadian pathway-chronotype’ and ‘arachidonic acid-intelligence’. Lastly, we applied least absolute shrinkage and selection operator (LASSO) regression to integrate HFS prediction score with SNP-based polygenic risk scores, which showed an improvement of 16.1–39.8% in cross-ancestry polygenic prediction. We concluded that HFS is a promising strategy for understanding the genetic basis of human complex traits.

    1. Computational and Systems Biology
    Qianmu Yuan, Chong Tian, Yuedong Yang
    Tools and Resources

    Revealing protein binding sites with other molecules, such as nucleic acids, peptides, or small ligands, sheds light on disease mechanism elucidation and novel drug design. With the explosive growth of proteins in sequence databases, how to accurately and efficiently identify these binding sites from sequences becomes essential. However, current methods mostly rely on expensive multiple sequence alignments or experimental protein structures, limiting their genome-scale applications. Besides, these methods haven’t fully explored the geometry of the protein structures. Here, we propose GPSite, a multi-task network for simultaneously predicting binding residues of DNA, RNA, peptide, protein, ATP, HEM, and metal ions on proteins. GPSite was trained on informative sequence embeddings and predicted structures from protein language models, while comprehensively extracting residual and relational geometric contexts in an end-to-end manner. Experiments demonstrate that GPSite substantially surpasses state-of-the-art sequence-based and structure-based approaches on various benchmark datasets, even when the structures are not well-predicted. The low computational cost of GPSite enables rapid genome-scale binding residue annotations for over 568,000 sequences, providing opportunities to unveil unexplored associations of binding sites with molecular functions, biological processes, and genetic variants. The GPSite webserver and annotation database can be freely accessed at https://bio-web1.nscc-gz.cn/app/GPSite.