1. Computational and Systems Biology
  2. Immunology and Inflammation
Download icon

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
Research Article
  • Cited 2
  • Views 701
  • Annotations
Cite this article as: eLife 2021;10:e68864 doi: 10.7554/eLife.68864

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.

Decision letter

  1. Joshua T Schiffer
    Reviewing Editor; Fred Hutchinson Cancer Research Center, United States
  2. Sara L Sawyer
    Senior Editor; University of Colorado Boulder, United States
  3. Joshua T Schiffer
    Reviewer; Fred Hutchinson Cancer Research Center, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Thank you for submitting your work entitled "Dynamically linking influenza virus infection with lung injury to predict disease severity" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Joshua T Schiffer.

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

The reviewers all appreciated the attention to an important question and the combination of experimental data and modeling to understand immune control of influenza in mice. However, the decision to reject the manuscript was based on two major factors:

1) The use of total CD8 rather than antigen-specific or 'active' (cytokine producing) cells.

2) The lack of statistical analysis in comparing different (and quite complex) models (especially since many of the conclusions rested upon model comparison).

As indicated, all reviewers were excited by the approach to combining novel experimental and mathematical approaches and felt that the validity of the conclusions and priority of the manuscript would increase if these issues could be addressed. However, ELife requires that authors are able to complete revisions in less than two months and it was concluded that it is unlikely that these concerns could be addressed within this timeframe.

I hope that this result will not discourage your work in this area or future submission to eLife.

Essential revisions:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Reviewer #1:

This is a tour de force paper describing a set of elegant experiments specifically designed for mathematical model testing and validation. The model takes the unprecedented step of linking kinetics of viral shedding, generation of infected cells, CD8+ T cell response and development of lung injury. Scientific conclusions are justified based on the analysis. The prediction that CD8 expansion but not levels of CD8+ T cells impact time to viral elimination is interesting and novel. The idea that infected cell density lowers Tcell-mediated clearance is also interesting and suggests a bottleneck to rapid immune clearance once a threshold of severity is surpassed.

Figures 3 and 4 are particularly instructive and interesting.

Overall, the paper could be substantially improved in terms of interpretation of the results and clarification of language regarding scientific conclusions. There are also a few slightly unclear sections.

1) In the introduction and conclusion, there is no mention that the balb c murine model may not capture the pathophysiology of influenza infection in humans. In particular, the relative contributions of humoral versus cell-mediated containment of infection may differ. The concepts of original antigenic sin and differential host susceptibility which may account for substantial heterogeneity in disease severity in human adults, are not captured in mice. The authors should directly admit this limitation and specifically acknowledge that the conclusions of their mathematical model (as elegant as they are) may not be fully generalizable to human infection. In addition, cited articles regarding flu pathogenesis should be labelled according to whether they are from humans, ferrets or mice. I would strongly consider reorganizing the discussion towards summarizing what conclusions might be relevant for human infections and what experimental data could be gathered to validate model predictions.

2) A key component of influenza pathogenesis is completely neglected, which is that a majority of infected people do not develop any lung disease. Instead, infection is limited to the upper airways. Therefore, a key role of the acquired immune response may be limiting spread from the upper to the lower airway. Under the best of circumstances, the author's model is relevant to only a subset of human cases. While the experimental system does not allow assessment of progression from upper to lower tract disease, this should be acknowledged as a major limitation of the system.

3) Figure 1: It is stated multiple times in the paper that Figure 1 demonstrates heterogeneity in viral loads. This is not the case as the standard deviation merely shows the confidence in the mean. Individual data points should be plotted to demonstrate heterogeneity in viral load and CD8+ T cell counts at each time point across mice.

4) Line 108: What percentage of infections are cleared within 4h? Are these included in the means and SDs in Figure 1? What is the proposed mechanism of these aborted infections?

5) A critical component of the model that is only mentioned in passing throughout the paper is the time delay in CD8 proliferation. The time delay parameter is as critical to the model as the density dependent CD8 killing rate (see Appendix 3 Figure A4). Yet the biology underpinning this delay is given short thrift in the discussion. It would also be useful to show that a model without this assumption fails to fit the data.

6) Table 1: are all parameter fitted? Are there references from the literature confirming that some of these values are realistic?

7) The term "model ensembles" described in line 148 is never adequately defined. The subsequent section of the paper is therefore confusing. Is this essentially describing the fact that several parameters are not identifiable because they are correlated to achieve model fit? This would not negate the model's validity but should at least be stated. Does Figure 2A only include parameter sets that resulted in optimal fit to the data (<5% error from the best fit) as described in the methods? Please clarify.

8) Line 320: the authors never show that small changes in viral load can lead to major changes in disease severity. They actually could do this with the sensitivity analysis output and it would be quite useful. However, without these simulations explicitly shown, this sentence should be removed.

9) Line 386: this is a false statement. See Figure 3 in https://www.nejm.org/doi/full/10.1056/NEJMoa1716197

Reviewer #2:

General assessment: The study presents a thorough combination of different experimental measurements with mathematical modeling to link viral dynamics and disease pathology in a mouse model of influenza infection. They find a very remarkable connection between the area of lung injured, the CD8+ T cell response and the development of disease symptoms of the animals. The discovered connections could advance the understanding of influenza virus infection in mice, but, in my view, the proposed predictive value should be further corroborated by appropriate analyses. I consider this as a very interesting study that is methodologically sound and innovative. However, I would have some comments (see below) that question the significance of this work for the field as requested by eLife.

1) A major concern affects the CD8+ T cell response used for modeling. It is stated that the total number of CD8+ T cells rather than the virus specific CD8+ T cells were used as the later often show varying dynamics (line 447-451). In which way would the use of virus-specific CD8+ T cells skew the results as mentioned? It would be an interesting question if the found density-dependent CD8+ T cell mediated clearance also holds if only virus-specific CD8+ T cells are considered. In addition, a population of memory cells was explicitly included in the model to reduce the number of CD8+ T cells responsible for clearing infected cells. This might not be necessary if virus-specific cells show a different dynamic than the total lung-resident CD8+ T cell population. In this regard, also the statement that the magnitude rather than the efficacy of the CD8 + T cells controls clearance could be questioned (line 98-99), as these aspects were not separately investigated in the experimental data nor the modeling framework.

2) Although the authors investigate the impact of the magnitude of the CD8+T cell response on the recovery time using their identified model (Figure 3), they do not use/extend these analyses to show the predicted changes for disease dynamics and pathology as based on the "workflow" shown in Figure 5. I think the study would substantially benefit if the postulated connectivity between disease dynamics and pathology, and the proposed impact of these findings on "forecasting disease progression, potential complications and therapeutic efficacy (line 23)" is shown at least theoretically for some scenarios (e.g. those used in Figure 3 that affect recovery time).

3) The analysis here benefits from the possibility to measure local viral loads and CD8+ T cell populations within the lungs of mice. I consider this as being rather difficult to be done in humans, as well as finding appropriate quantitative markers for disease pathology, such as weight loss. Therefore, I do not directly see the claimed ability of this study to enhance the ability to forecast disease progression and potential complications, at least not in humans. Even for mice it would be important to know if e.g. having only a limited number of measurements on the dynamics of the CD8+ T cell response and the viral load (e.g. until day 4 or 6 if this could be measured in vivo) is sufficient to parameterize the whole model appropriately in order to predict further dynamics.

Reviewer #3:

The authors present an interesting mix of modelling and experimental work looking at CD8 T cell control of influenza virus infection in mice. This extends previous work by looking at infected cell area, and coming to some slightly different conclusions on the mechanisms of T cell control. The strength of the conclusions is not well justified, and therefore the advance of previous work seems somewhat incremental.

1) Table 1 includes 21 parameters to fit CD8 and viral load. However, figure 1 suggests there are only 24 data points. This seems rather over-parameterized? This seems very evident when the authors justify the model, because for every potential inflection of each curve, they seem to add another parameter. But are all of these justified? Does adding additional parameters improve the fit (by AIC, for example)?

2) There appears no consideration of 'significance' in comparing the fit of different models. For example, the authors state (lines 127-135) that density dependence in clearance was a better fit than in CD8 expansion (as used previously), and just refer to the shape of curves on specific days. Surely something like an AIC for overall fit would be useful in comparing different models?

3) There seems a major confusion between 'total CD8' in lung, and virus-specific CD8. These seem used interchangeably. For example, line 166-168: 10% of (antigen specific) cells survive into memory, and here the authors comment they observe 17% (of peak total CD8). In the discussion (line 262 onwards) the authors discuss CD8E – referring to literature on antigen specific cells and comparing their results to these – when they have not measured this. This is justified because cells of different specificity have different kinetics. But since the total CD8 number is being used to predict the effects of antigen-specific cells (and this is the central conclusion of the manuscript) – surely some measure of what is functional (?cytokine secreting) or antigen-specific (tetramer or ICS positive) is necessary?

4) Model comparisons with data seem very vague. For example in Figure 4c the authors state (line 221-222) "the dynamics of the damaged cells of the lung correspond precisely to the dynamics of the maximum CD8". This sounds like it may be quantitative, but the figure just looks like they both peak around day 8 and decline. The CD8 increase on day 2 – whereas lesion does not increase until later. The same is true for the claim that AUC of infection matches active infection area (Figure 4B).

5) There seems relatively little explanation of the link between viral loads and lung lesions. For example, viral load peaks on day 2 and decreases thereafter. But lung lesion size is barely detectable at this time, and increases rapidly? The authors argue that AUC of virus = active area. How does this arise? Do infected cells produce a burst of virus and remain antigen positive for some time? What resolves antigen negative cells? How does lung lesion size (and %active inactive) directly relate to viral load? The authors make arguments around these issues, and figure 5 might lead one to believe there is some modelling to link all these, but there does not appear to be a mechanistic explanation?

6) The entire section "density dependent infected cell clearance" involves a lot of speculation on recovery time, CD8 thresholds etc. This seems entirely dependent on the model formulation, and average parameters (which are widely distributed). Is there any experimental evidence to support this modelling speculation?

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Dynamically Linking Influenza Virus Infection Kinetics, Lung Injury, Inflammation, and Disease Severity" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, including Joshua T Schiffer as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Sara Sawyer as the Senior Editor.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential Revisions:

1) Please include analyses in which δ_e does not approach its boundary value.

2) Please include AICs and model fits for models which are less supported by the data.

3) Please rewrite the abstract, introduction and discussion to highlight the scientific conclusions of the paper, rather than just the substantial technical achievements of the paper (fitting a model to several longitudinal data types in a complex and potentially representative model system of influenza).

4) Please frame the paper's conclusions and limitations more specifically in reference to their potential relevance for human infection. Specifically, highlight that the system employed in this paper is akin to a primary infection model rather than re-infection. Please also make an effort to partition past literature into mouse and human studies for better clarity.

Reviewer #1 (Recommendations for the authors):

The experimental data and modeling are highly robust. The conclusions of the paper are clearly supported by the results. The sensitivity analysis is particularly impressive and suggests a system that is highly conserved across a wide parameter space. Model validation with CD8+ depletion is a nice addition that leads to interesting and surprising conclusions. The figures are highly instructive and easy to read.

An area where the paper could be improved is conveying the actual scientific conclusions more clearly and precisely with more focused review of existing literature. The relevance of the paper's conclusions for human influenza could be discussed with more careful language.

First, the mechanistic conclusions of the work could be emphasized along with the methodology of the work. At present, these are completely lacking from the abstract which somewhat blandly just says that the paper describes a model which fits to data. From my perspective, currently underemphasized and novel / interesting conclusions are that:

1) CD8+ mediated killing becomes much more rapid on a per capita basis (40000 fold increase) when infected cells dip below several hundred cells approximately 7 days post infection.

2) There is a negative correlation between infected cell clearance by innate versus CD8+ mediated mechanisms, implying that poorer initial clearance of virus may result in more effective later killing by acquired immune mechanisms.

3) Even ~80% reduction in maximal CD8E+ levels could prolong infection by 10 days though delay in attaining these threshold CD8E+ levels due to experimental or in silico CD8+ depletion only delays viral elimination by a day.

4) Most interesting and counterintuitively, CD8+ depletion allows for considerable reductions in the size of lung lesions as well as inflammation scores and degree of weight loss during primary influenza infection. This result suggests that CD8+ T cells have the potential to create significant bystander damage in the lung.

Second, the introduction and discussion continue to not differentiate whether past experimental results are from humans or mice. It is somewhat misleading to cite mouse studies without acknowledging that these are from a model that in no way captures the totality of human infection conditions. For all animal models of human infection, the strengths of the model (ability to control experimental inputs and obtain frequent measurements) are counter-balanced by lack of realism. Humans have a complex background of immunity based on past vaccination and infection, different modes of exposure and other innumerable differences. In most human infections, the degree of lung involvement is minimal. Please stipulate in the review of existing literature which papers were done in mice versus humans. Please also frame conclusions of this paper in the discussion in terms of how it may or may not be relevant to human infection.

Third, this is a primary infection model, and this point also should be emphasized. The greatest relevance of the mouse model in the paper may be for pediatric infection in humans, rather than adults who have had multiple prior influenza exposures and possibly vaccinations. Presumably CD8+ responses can be expected to be more rapid with availability of a pre-existing population of tissue resident CD8+ T cells as would occur with re-infection. The results of CD8+ depletion prior to re-infection would potentially be very different (likely harmful) in a re-infection model and this should be discussed. This is mentioned in Line 467 but is given short attention elsewhere.

Line 60: stating that other studies have had limited success is rather insulting. Please rephrase and be more specific about why this study breaks new ground.

Line 81: "viral loads in the upper respiratory tract do not reflect the lower respiratory tract environment. " Please include a citation, remove or clarify that this is a possible confounding variable in the analysis.

Line 91: define lung histomorphometry. This is a fairly novel approach for most readers.

Line 101: This is a strong statement about viral load. Unless formal correlate studies have been done in humans (which they have not), I would day "may not be correlated" or remove altogether.

Line 201: involved with what? I am not sure what this sentence means.

Line 209: I would suggest denoting a separate section to the sensitivity analysis versus the parameter fitting as the fitted correlation between δ and δ_e appears separate mechanistically from the relationship between δ and viral clearance / total # of CD8E

Line 251: Please cite the clinical correlate oof this in the discussion. Immuncompromised humans often shed influenza (and SARS CoV-2) for months. See work from Jesse Bloom's group published in eLife on this subject.

Line 321 should this read "clear infected cells from the lung?" I am confused about what this sentence means.

Figure 5D: why are the dots yellow? Is the magenta line CD8 depleted?

Line 386: Has antiviral therapy been linked with extent of radiologic lung lesions in clinical trials. This would be a very atypical clinical trial endpoint so please be more precise with language. It is possible as previously mentioned in the paper that viral load may not predict lesion size or disease severity in humans.

Line 477: add degree of immunity from prior infections as a critical variable

Reviewer #2 (Recommendations for the authors):

This is a revised version of a previously reviewed article. The authors performed extensive additional analyses to address previous concerns and issues raised by the reviewers. However, there are a few additional points which, in my view, still would need some clarification:

1. As pointed out by one of the previous reviewers, the remarkable ability of the model to basically cover every change in the dynamics observed within the data (Figure 1) could suggest that the model is overfitting the data. In response, the authors mentioned that they performed robust fitting and sensitivity analyses, and also mentioned that they performed several different model attempts to reach at their final model. However, the current information provided, as e.g. in Figure 2 showing the parameter estimates, do not seem to fully support this claim. The authors state that "the majority of parameters are well-defined" with the exception of three parameters. However, also the death rate δ_E reaches the imposed boundary for fitting, which seems to be not addressed. As this parameter controls the CD8-mediated death rate and, thus, could be critical for the argument of a density-dependent death rate, I think this should be discussed/investigated in more detail. In addition, it could be very convincing if the AICs of some of the models tested that did not fit the data (i.e. reduced models as claimed within the response), are shown within the manuscript to support their claims. In addition, just as a suggestion, approaches that do not explicitly describe the mechanistic of the CD8+ T cell response but rather describe the measured responses by a spline function could be considered as well (see Kouyos et al. PLoS Comp Biol 2010), reducing the complexity of the model by still being able to examine the relationship between CD8+ T cell response and immunopathology.

2. Lung immunohistopathology is quantified based on tissue sections of individual mice. Did the authors analyze the whole (i.e. 3D) lung for regions of active/inactive lesions or only representative 2D tissue sections? In addition, how many mice/tissue sections were analyzed at each time point (e.g. Figure B/C)? It would be interesting to know how representative a 2D tissue section as shown in Figure 4A would be for the situation in a mouse, and also how representative it would be across mice. I apologize if I missed this information in the text but I think these details should be provided.

3. In Figure 5 B and D the difference between the fitted (black) and predicted (purple) relationships does not seem to be explained within the figure legend nor the text. I have difficulties to understand how these two things are related. The same holds true for Figure 5A, where there is no information on what the purple curve represents.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your thoroughly revised article "Dynamically Linking Influenza Virus Infection Kinetics, Lung Injury, Inflammation, and Disease Severity" for consideration by eLife. Your article has been reviewed by 1 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Sara Sawyer as the Senior Editor.

Essential Revisions:

We thank you for taking the effort to making excellent revisions and thoroughly addressing all reviewer comments. We apologize for not noticing this in the last revision but the Discussion section of the paper lacks a limitations paragraph and there are indeed several limitations that are not mentioned (lack of complete realism of the mouse model as a model of human infection; the fact that the weight loss equations that relate to inflammation are not "mechanistic" and therefore provide only some insight; other arms of the immune response are not studied in depth (humoral) experimentally and this system and may be important). A brief acknowledgment of these and other limitations, and possible next steps to address them, would really strengthen the impact of this paper.

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

Author response

Essential revisions:

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

Reviewer #1:

This is a tour de force paper describing a set of elegant experiments specifically designed for mathematical model testing and validation. The model takes the unprecedented step of linking kinetics of viral shedding, generation of infected cells, CD8+ T cell response and development of lung injury. Scientific conclusions are justified based on the analysis. The prediction that CD8 expansion but not levels of CD8+ T cells impact time to viral elimination is interesting and novel. The idea that infected cell density lowers Tcell-mediated clearance is also interesting and suggests a bottleneck to rapid immune clearance once a threshold of severity is surpassed.

Figures 3 and 4 are particularly instructive and interesting.

Overall, the paper could be substantially improved in terms of interpretation of the results and clarification of language regarding scientific conclusions. There are also a few slightly unclear sections.

We appreciate the careful reading and kind words about our work. As noted, both above and below, we have added a significant amount of new data and included discussion (see below) to address the concerns raised.

1) In the introduction and conclusion, there is no mention that the balb c murine model may not capture the pathophysiology of influenza infection in humans. In particular, the relative contributions of humoral versus cell-mediated containment of infection may differ. The concepts of original antigenic sin and differential host susceptibility which may account for substantial heterogeneity in disease severity in human adults, are not captured in mice. The authors should directly admit this limitation and specifically acknowledge that the conclusions of their mathematical model (as elegant as they are) may not be fully generalizable to human infection. In addition, cited articles regarding flu pathogenesis should be labelled according to whether they are from humans, ferrets or mice. I would strongly consider reorganizing the discussion towards summarizing what conclusions might be relevant for human infections and what experimental data could be gathered to validate model predictions.

The reviewer is correct that we had previously limited the discussion of human data, but it is an important discussion point. It is possible that some features could vary, particularly given that humans and influenza viruses are highly heterogeneous and immune responses change with age, sex, race, various co-morbidities, immunologic history, etc. This heterogeneity is difficult to control for in the laboratory or clinic, and clinical data, in particular, is insufficient to dissect the contribution of each. Further, it is quite difficult and invasive to obtain lung measurements (e.g., virus or CD8s) from living humans, where even CT imaging is a rare event. This makes the animal model highly desirable and the only system to begin establishing these connections. Fortunately, the mouse model is well established, recapitulates many aspects of human influenza disease, including a striking similarity in lung pathology (noted in Lines 67-68 and 382), and provides a highly controlled, reproducible system. We agree that there is much more to do to prove the human connection, but this is well beyond the scope of this manuscript. Nevertheless, we believe that the model and results will be similar, although may need adaptations due to pre-existing immunity (even if heterosubtypic), and resident T cells are only beginning to be understood.

2) A key component of influenza pathogenesis is completely neglected, which is that a majority of infected people do not develop any lung disease. Instead, infection is limited to the upper airways. Therefore, a key role of the acquired immune response may be limiting spread from the upper to the lower airway. Under the best of circumstances, the author's model is relevant to only a subset of human cases. While the experimental system does not allow assessment of progression from upper to lower tract disease, this should be acknowledged as a major limitation of the system.

The reviewer is correct that some influenza viruses and in some individuals only a mild, upper respiratory tract (URT) infection is established and that our data is not relevant for that case given our focus on the lung. Our workflow would generally not be appropriate for URT viral loads (e.g., nasal washes) because these do not reflect the lung environment or indicate clinical severity. There is copious data showing that many influenza virus strains (e.g., 2009 pH1N1) and certain individuals (e.g., elderly, pregnancy, immunocompromised, etc.) are more susceptible to lower respiratory tract (LRT) infections and manifestations. Other factors that contribute to developing LRT infections include dose, volume of inhaled inoculum, lung function and size, etc.. LRT infections are important to understand as they place the highest burden on public health and are a major cause of hospitalizations and influenza-related mortality. Importantly, CT scans of hospitalized patients show strikingly similar pathological features as our data. To highlight this, we have included discussion and references (Lines 67-68 and 382) to reflect this data.

With respect to the question regarding migration from URT to LRT, this question is related yet distinct from this particular data and some clinical situations as it is not necessarily a qualifying event for a lung infection. From an experimental standpoint, we do have the ability to mimic the migration using a low volume inoculum, which initially stays predominantly in the nasopharynx. High dose is not technically possible in low volume inoculation, so we have that as an additional confounder. Our collaborators have done this experiment for parainfluenza infection (e.g., Burke et al. (2015)), which is quite similar to influenza. Notably, individual animals were able to be tracked in that system using bioluminescence, so the migration was visually documented. In short, that data showed that the viral loads in the lung lag by ~2 days and resulted in a lower viral load with smaller areas of infection. The reduced infection and the reduced weight is precisely in line with our findings here. Similar reductions in infected cells and weight loss occurred in the animals that were CD8-depleted (Figure 5). Thus, the connections and results here are unlikely to change in that scenario. We recently submitted a manuscript analyzing the parainfluenza data, which the reviewer will certainly find of interest.

3) Figure 1: It is stated multiple times in the paper that Figure 1 demonstrates heterogeneity in viral loads. This is not the case as the standard deviation merely shows the confidence in the mean. Individual data points should be plotted to demonstrate heterogeneity in viral load and CD8+ T cell counts at each time point across mice.

We believe the reviewer is referencing standard error of the mean (SEM) rather than standard deviation (SD). In this type of data, SEM is not appropriate and we use SD. With SD, the heterogeneity in the data is directly captured and reflected as plotted. However, the data for each individual is included in our data files and the viral loads of individual animals were previously published in Smith et al. (2018) Front Microbiol. Here, we do not include a figure of the individual data due to the inclusion of the raw data as mandated by eLife and the limited space for additional figures.

4) Line 108: What percentage of infections are cleared within 4h? Are these included in the means and SDs in Figure 1? What is the proposed mechanism of these aborted infections?

100% of the animals have undetectable virus at 4 h. The means and SD are included. However, this is not clearance per se, but rather indicates that cells were infected but not yet producing sufficient virus to be measured by TCID50.

5) A critical component of the model that is only mentioned in passing throughout the paper is the time delay in CD8 proliferation. The time delay parameter is as critical to the model as the density dependent CD8 killing rate (see Appendix 3 Figure A4). Yet the biology underpinning this delay is given short thrift in the discussion. It would also be useful to show that a model without this assumption fails to fit the data.

The reviewer is correct that the delay is important. This delay reflects the time it takes for antigen presentation, expansion in the lymph, and trafficking to/from the lymph (Lines 671-675).

6) Table 1: are all parameter fitted? Are there references from the literature confirming that some of these values are realistic?

~14 parameters are estimated using a global search algorithm (see Table 1 and Lines 686-730). Most of the parameters are unknown biologically and some encompass several processes, and thus are not available in the literature. We did compare the ones that had to biologically known entities (e.g., percent memory cells generated (see Lines 229-233)). These were in agreement, which is a good indication that our parameters are realistic.

7) The term "model ensembles" described in line 148 is never adequately defined. The subsequent section of the paper is therefore confusing. Is this essentially describing the fact that several parameters are not identifiable because they are correlated to achieve model fit? This would not negate the model's validity but should at least be stated. Does Figure 2A only include parameter sets that resulted in optimal fit to the data (<5% error from the best fit) as described in the methods? Please clarify.

Parameter ensembles or model ensembles are collections of parameters that produce solutions within a 95% confidence interval of the best-fit. This collection of solutions is the gray shading in all model fits (Figures 2 and 4) and the parameter ‘clouds’ in Figure 3 and S2-S3. Some parameters are non-identifiable (e.g., eclipse phase, k, and memory delay, tm), which is evident from visualizing the parameter ensembles and plotting the histograms (see Figures S2-S3). That is, non-identifiable parameters will return uniform-like distributions where nearly every value works equally well. Correlated parameters, on the other hand, are not considered non-identifiable. These parameters are still bounded and have well defined distributions. We have shown in our coinfection work that correlated parameters do not inhibit our ability to accurately estimate their values (see Smith (2018) Immunol Rev for further discussion). In that work, we were able to experimentally validate a correlated parameter and its 95% CI. It is quite useful to know of their correlation, however, in order to accurately interpret the results. This is why we always plot the 2D parameter ensembles.

8) Line 320: the authors never show that small changes in viral load can lead to major changes in disease severity. They actually could do this with the sensitivity analysis output and it would be quite useful. However, without these simulations explicitly shown, this sentence should be removed.

This is directly shown in the 95% confidence intervals surrounding the model predictions of the lesion data. That is, parameters within the 95% CIs simultaneously generate the viral load curve (V), the CD8 curve (E+EM+E*), the active lesion via CAUC(I2), and the inactive lesion via E/Emax (Figures 1 and 4). Thus, one should be able to directly see that small changes in viral loads (gray band around the V curve in Figure 1) can lead to larger changes in the percent of the lung infected (gray band around CAUC(I2)).

9) Line 386: this is a false statement. See Figure 3 in https://www.nejm.org/doi/full/10.1056/NEJMoa1716197

In that article, baloxavir does slightly better than oseltamivir, but neither do a great job of reducing viral loads as compared to placebo. For example, in Figure 3A, the mean is reduced but the error bars are extremely large and overlap with the placebo group. This does not indicate a robust efficacy, at least in our opinion. Nevertheless, we updated this line to read: “This is particularly important because current antivirals alleviate symptoms but do not s effectively lower viral loads”.

Reviewer #2:

General assessment: The study presents a thorough combination of different experimental measurements with mathematical modeling to link viral dynamics and disease pathology in a mouse model of influenza infection. They find a very remarkable connection between the area of lung injured, the CD8+ T cell response and the development of disease symptoms of the animals. The discovered connections could advance the understanding of influenza virus infection in mice, but, in my view, the proposed predictive value should be further corroborated by appropriate analyses. I consider this as a very interesting study that is methodologically sound and innovative. However, I would have some comments (see below) that question the significance of this work for the field as requested by eLife.

We appreciate the reviewer’s assessment of our work and concerns. To help substantiate the claims, we have added a significant amount of new data as described both above and below.

1) A major concern affects the CD8+ T cell response used for modeling. It is stated that the total number of CD8+ T cells rather than the virus specific CD8+ T cells were used as the later often show varying dynamics (line 447-451). In which way would the use of virus-specific CD8+ T cells skew the results as mentioned? It would be an interesting question if the found density-dependent CD8+ T cell mediated clearance also holds if only virus-specific CD8+ T cells are considered. In addition, a population of memory cells was explicitly included in the model to reduce the number of CD8+ T cells responsible for clearing infected cells. This might not be necessary if virus-specific cells show a different dynamic than the total lung-resident CD8+ T cell population. In this regard, also the statement that the magnitude rather than the efficacy of the CD8 + T cells controls clearance could be questioned (line 98-99), as these aspects were not separately investigated in the experimental data nor the modeling framework.

To address any potential differences between use of the total CD8s and IAV-specific CD8s, we performed additional experiments that are described at the beginning of this document and in the manuscript (see Figure 1). In short, we detailed IAV-specific CD8s, including the effector phenotype CD25-CD43+ and ensured they were present in the lung parenchyma. Because they have similar dynamics as the total, our results were not skewed. The density-dependence was still necessary to include in the model, and only select parameters were altered due (see Table 1). Plotting the IAV-specific CD8s against the % inactive lesion also showed that the nonlinearity is even stronger with these cells (Figure S4/Supplementary file 1D). In addition, we further tested our model by depleting CD8s (Figure 3). The model fit these data and highlighted the density dependence (noted in Lines 270-271). With respect to our statement about the efficacy versus the magnitude, we do investigate this through a sensitivity analysis (Appendix 3) and now through depletion of the cells. Our findings suggest that the efficacy parameter was insensitive (see Appendix 3) and that this parameter was not changed during depletion even though the population level did change (see Table S1).

2) Although the authors investigate the impact of the magnitude of the CD8+T cell response on the recovery time using their identified model (Figure 3), they do not use/extend these analyses to show the predicted changes for disease dynamics and pathology as based on the "workflow" shown in Figure 5. I think the study would substantially benefit if the postulated connectivity between disease dynamics and pathology, and the proposed impact of these findings on "forecasting disease progression, potential complications and therapeutic efficacy (line 23)" is shown at least theoretically for some scenarios (e.g. those used in Figure 3 that affect recovery time).

The reviewer highlights an interesting point that was not included in our original submission. In response, we performed an additional experiment to specifically investigate the recovery time (described in detail at the beginning of this document, in the previous point, and in the manuscript text). This experiment is slightly different than our in silico experiment where we had altered the expansion rate. Isolating and changing rates in vivo is extraordinarily difficult, but altering the population levels is relatively simple. In short, we depleted CD8s during the infection and showed how virus, CD8s (Figure 3F), and weight loss (Figure 5A) change. Because viral loads were lower by 2 d π and because CD8 depletion antibodies are known to alter various aspects of the response, we refit the model to identify potential differences in model parameters (see Figure 3 and Table S1). Our analysis suggested 3 differences: lower T0, lower x, and higher h. Plotting the relation between the activated CD8s and virus also suggested a lower x and higher h (see Figure 3G). We chose to not perform the histormorphometry for this experiment to verify the lower T0 because it would have doubled the number of animals, which was not ethically justifiable. Instead, we used the opportunity to test our workflow by predicting the changes in lung lesions and inflammation (Figure 4). Then, we tested the correlation between these and weight loss (Figure 5). The primary Hill function parameters (n’s and K’s) were unaltered and the maximum (lmax’s) was obtained directly from the estimated lesion and inflammation values. These results show that the workflow is robust, and indirectly validate the effect of depletion on T0.

3) The analysis here benefits from the possibility to measure local viral loads and CD8+ T cell populations within the lungs of mice. I consider this as being rather difficult to be done in humans, as well as finding appropriate quantitative markers for disease pathology, such as weight loss. Therefore, I do not directly see the claimed ability of this study to enhance the ability to forecast disease progression and potential complications, at least not in humans. Even for mice it would be important to know if e.g. having only a limited number of measurements on the dynamics of the CD8+ T cell response and the viral load (e.g. until day 4 or 6 if this could be measured in vivo) is sufficient to parameterize the whole model appropriately in order to predict further dynamics.

The reviewer raises interesting questions and is correct that measuring viral loads or CD8s from the human lung is quite challenging. CD8s can be measured from the blood, but these don’t directly reflect the lung environment. Although influenza does replicate well in the upper respiratory tract (URT), which is easily accessible, URT measurements are limited in their utility. This is because similar viral loads can be observed in patients with different outcomes due to the URT being an initial contact site. In some cases, patients have x-rays or CT scans performed. CT scans look strikingly similar to our lung histomorphometry (noted in Lines 67-68 and 382), which gives us hope that our findings here would be translatable to humans. Another challenge will be finding a symptom that relates in the same way as weight loss. Humans experience a range of symptoms (e.g., fever, cough, shortness of breath, fatigue, etc.), which do not necessarily have the same timescales or act as the best predictor of disease (noted in Lines 104-106). Measurement of these can also be quite subjective. These difficulties highlight the importance of animal and quantitative studies like ours to begin establishing these dynamics and connections.

The question about whether less data could be used is interesting. Our CD8 depletion data offered us the possibility to analyze the data with fewer time points (5 time points instead of 13) and fewer measurements (no lesions or inflammation scoring). Fewer time points does reduce the number of parameters that one would be able to estimate with statistical confidence. Nevertheless, the results showed that our workflow could still be used. Further, had the depletion antibody not been known to affect the rates of the CD8 response, we likely could have inferred their dynamics without data. Nevertheless, these new data and analyses do support the utility of less data.

Reviewer #3:

The authors present an interesting mix of modelling and experimental work looking at CD8 T cell control of influenza virus infection in mice. This extends previous work by looking at infected cell area, and coming to some slightly different conclusions on the mechanisms of T cell control. The strength of the conclusions is not well justified, and therefore the advance of previous work seems somewhat incremental.

We have updated our manuscript in numerous ways to better highlight the novelty and importance of this work. As noted above, we added a significant amount of data and analyses that further increase the novelty and substantiate the model, analyses, and findings. The reviewer is correct that we did follow up our own work on a simplified density-dependent model (Smith et al. (2018) Front Microbiol). We described this and highlighted that we came to the same conclusion with respect to the density-dependent clearance (see Lines 42, 167, 399-408, and Appendix 2). In addition, we added further insight into this term and how to appropriate interpret the smaller model (see Appendix 2). As the reviewer mentioned, this conclusion is distinct from some models, which assume the density-dependence is in the effector expansion. We investigated those models and clearly show that these other models do not fit our data (Appendix 1).

While this small portion of our work could be considered as incremental, there are numerous other aspects that significantly set our work apart from others with respect to both data, modeling, and insight. These include (1) quantitative data and model fits to various phenotypes of CD8s, (2) validated CD8 model and findings, (3) validated infected cell dynamics, (4) quantitative data and modeling of CD8 depletion, (5) quantitative data and modeling of lung injury, (6) quantitative data and modeling of inflammation, and (7) connection of both pathologic measurements with weight loss along with validation. Some studies, including some of our own, have attempted to predict immunopathology or damage, often using these interchangeably, but our data highlight their inaccuracies and clearly show that these are distinct. We are aware of a two studies that did attempt to follow up on our own novel suggestion to use weight loss data (see Smith and Perelson (2011) WIRES Sys Biol). However, one had (self-proclaimed) failure (Price et al. (2015) J Theor Biol) and another did not model any other feature to substantiate their equation for symptom scores (Manchanda et al. (2014) Biosystems). These are discussed in Lines 57-60 and 426-428. If there are other papers we have inadvertently missed that were published prior to the publication of our preprint (Feb 2019), we would be happy to include those. We hope this explanation and our new data and analyses more clearly show the numerous novelties and importance of this work.

1) Table 1 includes 21 parameters to fit CD8 and viral load. However, figure 1 suggests there are only 24 data points. This seems rather over-parameterized? This seems very evident when the authors justify the model, because for every potential inflection of each curve, they seem to add another parameter. But are all of these justified? Does adding additional parameters improve the fit (by AIC, for example)?

In Figure 1B, there are a total of 260 data points across 13 time points and 2 populations (virus and CD8 T cells). There are 14 estimated parameters in Table 1. Thus, the model is not overparameterized. Of note, the general rule of 1 parameter per 1 time point does not apply when multiple variables are involved, so numerous additional parameters could have also been fit without concern of overfitting. We’re unsure what is meant by “for every potential inflection of each curve, they seem to add another parameter”. We do not choose parameters based on inflection points. We build our models with extreme care and with the minimum number of terms and parameters needed to fit and explain the data. Our general model building algorithm is to develop it term-by-term while simultaneously testing various functional forms, beginning always with linear terms then moving to nonlinear terms if necessary. At the end, we have typically gone through dozens to hundreds of different model formulations. We do not include these in a manuscript because that would be extremely cumbersome and unnecessary. In addition, we perform robust fitting and sensitivity analyses (Figures 1-3, A2-A4) so that we and our readers have the knowledge needed to appropriately interpret the model and its results.

2) There appears no consideration of 'significance' in comparing the fit of different models. For example, the authors state (lines 127-135) that density dependence in clearance was a better fit than in CD8 expansion (as used previously), and just refer to the shape of curves on specific days. Surely something like an AIC for overall fit would be useful in comparing different models?

We do consider significance when evaluating models but only do so when it is statistically appropriate. We did not initially consider them when evaluating alternate models because it is not technically sound to compare model fits that use different data. Nevertheless, in response to the reviewer, we did do this by first including all data even though the alternate model was not fit to d3-5 and then excluding d3-5 even though our model was fit to these. We further examined how the models faired against the lung-specific CD8 data. In all cases, the AICs of the alternate models were statistically worse and could not be supported. We now include this analysis in Table S2/Supplementary file 1G along with a direct visual comparison of the alternate models in Figure A1.

3) There seems a major confusion between 'total CD8' in lung, and virus-specific CD8. These seem used interchangeably. For example, line 166-168: 10% of (antigen specific) cells survive into memory, and here the authors comment they observe 17% (of peak total CD8). In the discussion (line 262 onwards) the authors discuss CD8E – referring to literature on antigen specific cells and comparing their results to these – when they have not measured this. This is justified because cells of different specificity have different kinetics. But since the total CD8 number is being used to predict the effects of antigen-specific cells (and this is the central conclusion of the manuscript) – surely some measure of what is functional (?cytokine secreting) or antigen-specific (tetramer or ICS positive) is necessary?

We apologize for the confusion. To make this more clear, we added several new pieces of data and model fits. We’ll keep our response here brief, but a more detailed description in included at the beginning of this document and in the manuscript. We added data that define the difference between the total CD8, lung-specific CD8, and specific effector and memory phenotypes (Figures 1, S1). These data support our original approach, showed that our model predicted dynamics for CD8E were reflective of CD8s in the lung parenchyma, and showed how parameters change if we assume only certain phenotypes are in action (Table 1).

4) Model comparisons with data seem very vague. For example in Figure 4c the authors state (line 221-222) "the dynamics of the damaged cells of the lung correspond precisely to the dynamics of the maximum CD8". This sounds like it may be quantitative, but the figure just looks like they both peak around day 8 and decline. The CD8 increase on day 2 – whereas lesion does not increase until later. The same is true for the claim that AUC of infection matches active infection area (Figure 4B).

Figure 4 is not a model fit but rather an overlay of the model prediction and the data. We specifically use the verbiage “correspond” and “match” within the text to indicate that it is an overlay of the model and data rather than a fit. To make this more clear, we added the words “fit” and “predicted” to most figures so that the reader can clearly distinguish these. Figure 4, specifically, shows that they correspond and that the model accurately predicts these lung pathologic features. It should be clear that the active lesion is the CAUC rather than AUC. In addition, it is not simply the CD8s, which do peak at day 8, but rather the relative CD8s. To more clearly illustrate the nonlinear relation between CD8s and the inactive lesion, we also added Figure S4F:

5) There seems relatively little explanation of the link between viral loads and lung lesions. For example, viral load peaks on day 2 and decreases thereafter. But lung lesion size is barely detectable at this time, and increases rapidly? The authors argue that AUC of virus = active area. How does this arise? Do infected cells produce a burst of virus and remain antigen positive for some time? What resolves antigen negative cells? How does lung lesion size (and %active inactive) directly relate to viral load? The authors make arguments around these issues, and figure 5 might lead one to believe there is some modelling to link all these, but there does not appear to be a mechanistic explanation?

We apologize for the confusion regarding our data. The technical description of histomorphometry is provided in the Methods (Lines 622-638). We included more detail below to further explain type of data and experimental procedures:

In our histomorphometric staining, the cells are stained with an antibody against the influenza nucleoprotein (abbreviated as NP). The labeled NP+ cells are infected cells, which should not be confused with full virions or viral loads. Infected cells express all influenza proteins, so they are able to be labeled. NP is expressed in abundance and thus widely used for staining. Virions do not express internal proteins like NP on their surface, so these are not stained in the process. Further, the influenza virus is an enveloped virus, which means that the virus requires an intact, living cell to produce virions (no bursting). Viral loads are a measure of extracellular infectious virus, and we quantify them using TCID50. Because of these differences, we never claim that “AUC of virus = active area”. Rather, it is the CAUC, which is distinct from AUC, of productively infected cells (CAUC(I2)) because histomorphometry is a measurement of the infected cell area. Obviously, free virus is what causes a cell to be infected (one could consider this the “mechanism”), so virus is indirectly related to the histomorphometry measurement. The indirect nature is noted in Figure 6 by an arrow from virus to the model (to estimate the infected cell dynamics) and from the model to the histomorphometry. “Inactive” lesions are areas that were previously infected. Once a cell becomes infected, it will either die from infection or be removed by an immune cell (one could consider this the “mechanism”). Both the active and inactive areas are defined and quantified by a board-certified pathologist who is blinded from the study.

6) The entire section "density dependent infected cell clearance" involves a lot of speculation on recovery time, CD8 thresholds etc. This seems entirely dependent on the model formulation, and average parameters (which are widely distributed). Is there any experimental evidence to support this modelling speculation?

There is prior biological evidence that the rate is density-dependent, which is discussed on Lines 49-52. We did enhance this discussion and included additional references (within Lines 409-431). In addition, our analysis in Figure S4F (noted above), our data and model fit to the CD8 depletion data (Figure 3), the close comparison of our model with the lung pathology (Figure 4), and the lack of a robust fit when this density dependence is excluded (Appendix 1) provide further support. In our opinion, this is quite a bit of supportive evidence and we do not have reason to believe that it would be dependent on model formulation. Where one would need to be cautious is only fitting the CD8 peak and expecting the model to be correct.

[Editors’ note: what follows is the authors’ response to the second round of review.]

Essential Revisions:

1) Please include analyses in which δ_e does not approach its boundary value.

Below is the ensemble plot of dE and d with larger bounds on dE (no other parameters are significantly altered). This shows that d approaches zero with increasing values of dE. Mathematically, this is due to the viral loads during the first decay phase (days 2-6) having little dynamics and thus being insufficient to distinguish between two mechanisms of infected cell clearance. Biologically, it would suggest that CD8s are the only mechanism of infected cell clearance throughout the entire infection, which we believe is incorrect. Thus, in restricting dE, we are able to maintain biological relevancy. In addition, this allows us to keep our results consistent with prior publications where d has been shown to dictate the slope in the first decay phase. That is, the viral load slope is -0.2 log10 TCID50/day (Smith et al. (2018) Front Microbiol), which equates to d=0.46 (i.e., V(t)=e-dt as derived in Smith et al. (2010) J Math Biol) in the absence of other mechanisms. Generally speaking, restricting parameters to biologically meaningful ranges is a relatively common practice (e.g., k is a prime example where numerous papers have either fixed or restricted the range of this parameter). We did not include this figure in our revised version as we feel its inclusion does not add value to the manuscript, but would be happy to if the Editors or Reviewer think it would be beneficial. We did add a brief explanation about the parameter restriction in Lines 750-752.

Author response image 1

2) Please include AICs and model fits for models which are less supported by the data.

Fits and AICs were previously included in Appendix-Figure 1 and Supplementary file 1 (previously Supplementary file 1G; tables are not supported in eLife’s appendix box within LaTeX). These are referenced in Lines 194, 496, 691, 695, and 699.

3) Please rewrite the abstract, introduction and discussion to highlight the scientific conclusions of the paper, rather than just the substantial technical achievements of the paper (fitting a model to several longitudinal data types in a complex and potentially representative model system of influenza).

We updated the abstract, the last paragraph of the introduction (Lines 130-147), and portions of the Discussion (Lines 401-406, 414-416, 423-427, 446-450, 500-502) to better highlight our results.

4) Please frame the paper's conclusions and limitations more specifically in reference to their potential relevance for human infection. Specifically, highlight that the system employed in this paper is akin to a primary infection model rather than re-infection. Please also make an effort to partition past literature into mouse and human studies for better clarity.

Thank you for the suggestions. We updated the text to better highlight the human relevance and distinguish results between different host species in Lines 43, 61, 69, 74, 77, 88, 99, 101, 111, 116, 118, 120, 400-406, 410-411, 413, 435, 437, 446-450, 452, 463, 471, 485, and 846. We use the words ‘human’, ‘clinical’, and ‘patient’ to denote humans and the words ‘animal’, ‘murine’, and ‘experimental’ to denote animal models. Further distinguishment in some areas comes from the methods noted (e.g., CT scans/imaging is used/relevant for humans but not animals (this would be a microCT)). Because many features of the infection are observed in both humans and animals (see, for example, 10.1016/j.jim.2014.03.023 and 10.3390/pathogens3040845), underscoring the strong relevance of animal models to study influenza (noted in Lines 453-454), we limited specifying this in every sentence as we feel it would reduce the readability. In addition, we feel that splitting up the references mid-sentence would also reduce readability given the journal’s reference style, and have left most at the end of a sentence. The references themselves should provide a reader to easily distinguish.

Reviewer #1 (Recommendations for the authors):

The experimental data and modeling are highly robust. The conclusions of the paper are clearly supported by the results. The sensitivity analysis is particularly impressive and suggests a system that is highly conserved across a wide parameter space. Model validation with CD8+ depletion is a nice addition that leads to interesting and surprising conclusions. The figures are highly instructive and easy to read.

An area where the paper could be improved is conveying the actual scientific conclusions more clearly and precisely with more focused review of existing literature. The relevance of the paper's conclusions for human influenza could be discussed with more careful language.

Thank you for the suggestions. As mentioned above, we updated the abstract and the text to better highlight the biological conclusions in addition to the mathematical conclusions. We also included some additional explanation on the specific points below.

First, the mechanistic conclusions of the work could be emphasized along with the methodology of the work. At present, these are completely lacking from the abstract which somewhat blandly just says that the paper describes a model which fits to data. From my perspective, currently underemphasized and novel / interesting conclusions are that:

1) CD8+ mediated killing becomes much more rapid on a per capita basis (40000 fold increase) when infected cells dip below several hundred cells approximately 7 days post infection.

2) There is a negative correlation between infected cell clearance by innate versus CD8+ mediated mechanisms, implying that poorer initial clearance of virus may result in more effective later killing by acquired immune mechanisms.

3) Even ~80% reduction in maximal CD8E+ levels could prolong infection by 10 days though delay in attaining these threshold CD8E+ levels due to experimental or in silico CD8+ depletion only delays viral elimination by a day.

In our CD8 depletion data, the entire infection is altered (see Lines 270-303). That is, our results suggest that there are fewer infected cells initially infected, which leads to lower viral loads at d2 and will automatically result in fewer CD8s. However, the depletion antibody itself is known to directly alter CD8s (see Lines 282-284), so one cannot make direct, quantified conclusions about the precise reduction percentage under this experimental condition.

4) Most interesting and counterintuitively, CD8+ depletion allows for considerable reductions in the size of lung lesions as well as inflammation scores and degree of weight loss during primary influenza infection. This result suggests that CD8+ T cells have the potential to create significant bystander damage in the lung.

While bystander damage may be present, our CD8 depletion data do not directly show bystander damage. As we mention above, one important aspect of depleting CD8s prior to infection was that it reduced the viral loads early on. We believe this is likely a result of immune activation from the initial CD8 kill-off (noted in Lines 284-285). A decrease in target/infected cells will automatically reduce the number of total cells that become infected and, thus, reduce the lesioned area of the lung. This was verified by the reduced weight loss, and no modifications were made to the predictions (Figure 4) or correlation to weight loss (Figure 5).

Second, the introduction and discussion continue to not differentiate whether past experimental results are from humans or mice. It is somewhat misleading to cite mouse studies without acknowledging that these are from a model that in no way captures the totality of human infection conditions. For all animal models of human infection, the strengths of the model (ability to control experimental inputs and obtain frequent measurements) are counter-balanced by lack of realism. Humans have a complex background of immunity based on past vaccination and infection, different modes of exposure and other innumerable differences. In most human infections, the degree of lung involvement is minimal. Please stipulate in the review of existing literature which papers were done in mice versus humans. Please also frame conclusions of this paper in the discussion in terms of how it may or may not be relevant to human infection.

As mentioned above, we updated the text to better highlight the human relevance and distinguish results between different host species in Lines 43, 61, 69, 74, 77, 88, 99, 101, 111, 116, 118, 120, 400-406, 410-411, 413, 435, 437, 446-450, 452, 463, 471, 485, and 846. We use the words ‘human’, ‘clinical’, and ‘patient’ to denote humans and the words ‘animal’, ‘murine’, and ‘experimental’ to denote animal models. Further distinguishment in some areas comes from the methods noted (e.g., CT scans/imaging is used/relevant for humans but not animals (this would be a microCT)). Because many features of the infection are observed in both humans and animals (see, for example, 10.1016/j.jim.2014.03.023 and 10.3390/pathogens3040845), underscoring the strong relevance of animal models to study influenza (noted in Lines 453-454), we limited specifying this in every sentence as we feel it would reduce the readability. In addition, we feel that splitting up the references mid-sentence would also reduce readability given the journal’s reference style, and have left most at the end of a sentence. The references themselves should provide a reader to easily distinguish.

Third, this is a primary infection model, and this point also should be emphasized. The greatest relevance of the mouse model in the paper may be for pediatric infection in humans, rather than adults who have had multiple prior influenza exposures and possibly vaccinations. Presumably CD8+ responses can be expected to be more rapid with availability of a pre-existing population of tissue resident CD8+ T cells as would occur with re-infection. The results of CD8+ depletion prior to re-infection would potentially be very different (likely harmful) in a re-infection model and this should be discussed. This is mentioned in Line 467 but is given short attention elsewhere.

We added text to highlight that we are studying a primary infection in Lines 46, 130, and 176-177. In general, primary infections may not necessarily always equate to children as any novel strain or strain novel to that individual may act as a primary infection. It is also feasible that waning immunity would appear similar to the dynamics of a primary infection.

Because CD8 depletion would also deplete out resident and other T cells (>99% efficiency as noted in Line 273), one would expect the exact same results as we showed here. Thus, it would not be the appropriate experimental design to study recall responses. As we mentioned in our prior response, we might expect some adjustments to account primed responses and added text that highlights this in Lines 424-427 and 498.

Line 60: stating that other studies have had limited success is rather insulting. Please rephrase and be more specific about why this study breaks new ground.

We did not mean to be insulting and reworded Lines 66-67 to “…but have had not yet found the appropriate mathematical relation with the available data.” Of note, even the authors of Price et al. noted their inability to capture the dynamics: “Some trajectories, notably activated macrophages and epithelial damage are not well captured by the model, suggesting that the immunophenotype we selected for active macrophages may not be accurate, and that using animals’ weight as a proxy for epithelial damage may not be appropriate”. The novelty of our study and the gaps in the field are stated throughout the introduction and discussion.

Line 81: "viral loads in the upper respiratory tract do not reflect the lower respiratory tract environment. " Please include a citation, remove or clarify that this is a possible confounding variable in the analysis.

We’ve added several citations from both human and animal studies to Lines 90-91.

Line 91: define lung histomorphometry. This is a fairly novel approach for most readers.

This was defined in our prior revision and remains in Lines 102-104.

Line 101: This is a strong statement about viral load. Unless formal correlate studies have been done in humans (which they have not), I would day "may not be correlated" or remove altogether.

We updated Line 110 to “…may not be directly correlated…”.

Line 201: involved with what? I am not sure what this sentence means.

We were referencing effector-mediated killing and memory generation as noted earlier in the sentence. We updated Lines 210-212 to clarify: “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”

Line 209: I would suggest denoting a separate section to the sensitivity analysis versus the parameter fitting as the fitted correlation between δ and δ_e appears separate mechanistically from the relationship between δ and viral clearance / total # of CD8E

We appreciate the suggestion but have chosen to leave this part of the text unaltered.

Line 251: Please cite the clinical correlate oof this in the discussion. Immuncompromised humans often shed influenza (and SARS CoV-2) for months. See work from Jesse Bloom's group published in eLife on this subject.

The suggested article has been added to Line 474.

Line 321 should this read "clear infected cells from the lung?" I am confused about what this sentence means.

Line 332 has been updated to “…of the infected areas within the lung”.

Figure 5D: why are the dots yellow? Is the magenta line CD8 depleted?

We inadvertently left off the explanation, but added clarification to the caption and text (Lines 371374 and 380-384). The yellow markers are interstitial inflammation while the white markers are alveolar inflammation. The magenta markers/lines are the CD8 depletion prediction.

Line 386: Has antiviral therapy been linked with extent of radiologic lung lesions in clinical trials. This would be a very atypical clinical trial endpoint so please be more precise with language. It is possible as previously mentioned in the paper that viral load may not predict lesion size or disease severity in humans.

To our knowledge, CT images are not typically taken in many contexts for a variety of clinical reasons (e.g., cost, exposure to the patient, etc.), but antivirals have been linked to reductions in disease severity (e.g., see 10.1056/NEJMoa1716197). In that particular line, we mention that minor reductions in viral load are paired with more significant reductions in disease/symptom, which is reported in the referenced clinical and experimental data.

Line 477: add degree of immunity from prior infections as a critical variable

This has been added to Line 478.

Reviewer #2 (Recommendations for the authors):

This is a revised version of a previously reviewed article. The authors performed extensive additional analyses to address previous concerns and issues raised by the reviewers. However, there are a few additional points which, in my view, still would need some clarification:

1. As pointed out by one of the previous reviewers, the remarkable ability of the model to basically cover every change in the dynamics observed within the data (Figure 1) could suggest that the model is overfitting the data. In response, the authors mentioned that they performed robust fitting and sensitivity analyses, and also mentioned that they performed several different model attempts to reach at their final model. However, the current information provided, as e.g. in Figure 2 showing the parameter estimates, do not seem to fully support this claim. The authors state that "the majority of parameters are well-defined" with the exception of three parameters. However, also the death rate δ_E reaches the imposed boundary for fitting, which seems to be not addressed. As this parameter controls the CD8-mediated death rate and, thus, could be critical for the argument of a density-dependent death rate, I think this should be discussed/investigated in more detail.

We appreciate the reviewer noting the great fits, which are the accumulation of several years of model development and parameter estimation. As we mentioned in our previous response to the reviewer’s concern about overfitting, we have copious data to fit the model without concern for overfitting. We believe the reviewer may be confusing over-/under-fitting with parameter identifiability and parameter correlations. These three parameter estimation features are not equivalent nor do the latter two equate to any issue in fitting or in making robust conclusions. Being fortunate enough to experimentally validate our models, we have previously shown that accurate conclusions and parameter estimates can be made in the face of parameter correlations (see our prior body of work on influenza-pneumococcal infections, which is summarized in Smith (2018) Immunol Rev). Again here, we have successfully validated our model with several pieces of data. Generally speaking, nearly every model-data pairing will result in correlated parameters and some parameters that are non-identifiable. While we could fix non-identifiable parameters to a chosen value, we would obtain the same results as if we fit them due to their non-identifiability (i.e., every values works equally well). We prefer to let these parameters fit for times when this turns out to not be the case.

With respect to dE, this parameter is only correlated to d as one might expect (i.e., both are infected cell clearance parameters). As we noted at the beginning of this document, we have previously shown that d can be directly estimated from the slope of the viral load data during the first decay phase (see Smith et al. (2010) J Math Biol for the original derivation). Allowing larger values of dE results in d approaching zero and would suggest that the CD8s are the only mechanism of infected cell clearance, which we believe to be untrue. We included a short explanation about the parameter restriction is in Lines 750-752.

In addition, it could be very convincing if the AICs of some of the models tested that did not fit the data (i.e. reduced models as claimed within the response), are shown within the manuscript to support their claims.

Fits and AICs were previously included in Appendix-Figure 1 and Supplementary file 1 (previously Supplementary file 1G; tables are not supported in eLife’s appendix box within LaTeX). These are referenced in Lines 194, 496, 691, 695, and 699.

In addition, just as a suggestion, approaches that do not explicitly describe the mechanistic of the CD8+ T cell response but rather describe the measured responses by a spline function could be considered as well (see Kouyos et al. PLoS Comp Biol 2010), reducing the complexity of the model by still being able to examine the relationship between CD8+ T cell response and immunopathology.

We appreciate the suggestion, and the reviewer is correct that splines are one possibility. However, we find these to be inferior due to their non-mechanistic structure, and have chosen to not explore this path because our model is quite simplistic while simultaneously being mechanistic and accurate.

2. Lung immunohistopathology is quantified based on tissue sections of individual mice. Did the authors analyze the whole (i.e. 3D) lung for regions of active/inactive lesions or only representative 2D tissue sections? In addition, how many mice/tissue sections were analyzed at each time point (e.g. Figure B/C)? It would be interesting to know how representative a 2D tissue section as shown in Figure 4A would be for the situation in a mouse, and also how representative it would be across mice. I apologize if I missed this information in the text but I think these details should be provided.

We analyzed 2D tissue sections (Line 660) from 5 animals (Lines 596 and 668). The choice of image per group was provided by the pathologist, who is blinded from the study, as a representative image. Other animals may have lesions located in different parts of the lung (as the exact location of cells becoming infected will happen by chance), but the percent infected areas are relatively similar. The error bars on the quantified data in Figure 4 can be directly used to infer how these change between animals. The complete experimental details are in Lines 653669.

3. In Figure 5 B and D the difference between the fitted (black) and predicted (purple) relationships does not seem to be explained within the figure legend nor the text. I have difficulties to understand how these two things are related. The same holds true for Figure 5A, where there is no information on what the purple curve represents.

Our apologies that we inadvertently left the explanation out of the caption and the text. Figures 5B, 5D are now explained in the caption and text (Lines 371-374 and 380-384).

[Editors’ note: what follows is the authors’ response to the second round of review.]

Essential Revisions (for the authors):

We thank you for taking the effort to making excellent revisions and thoroughly addressing all reviewer comments. We apologize for not noticing this in the last revision but the Discussion section of the paper lacks a limitations paragraph and there are indeed several limitations that are not mentioned (lack of complete realism of the mouse model as a model of human infection; the fact that the weight loss equations that relate to inflammation are not "mechanistic" and therefore provide only some insight; other arms of the immune response are not studied in depth (humoral) experimentally and this system and may be important). A brief acknowledgment of these and other limitations, and possible next steps to address them, would really strengthen the impact of this paper.

We added discussion points to Lines 407-408, 446-447, 450-451, and 456. Many limitations and potential future studies were previously acknowledged in Lines 426-429, 442-446, 447-450, 451-455, 484-487, 502-503, 520-522, 539-540, 547-549, and 552-555. We apologize that we did not point these out more vehemently in our prior revision.

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

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).

Senior Editor

  1. Sara L Sawyer, University of Colorado Boulder, United States

Reviewing Editor

  1. Joshua T Schiffer, Fred Hutchinson Cancer Research Center, United States

Reviewer

  1. Joshua T Schiffer, Fred Hutchinson Cancer Research Center, United States

Publication 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

  • 701
    Page views
  • 122
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

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)

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

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

Further reading

    1. Computational and Systems Biology
    Christopher P Mancuso et al.
    Research Article Updated

    Environmental disturbances have long been theorized to play a significant role in shaping the diversity and composition of ecosystems. However, an inability to specify the characteristics of a disturbance experimentally has produced an inconsistent picture of diversity-disturbance relationships (DDRs). Here, using a high-throughput programmable culture system, we subjected a soil-derived bacterial community to dilution disturbance profiles with different intensities (mean dilution rates), applied either constantly or with fluctuations of different frequencies. We observed an unexpected U-shaped relationship between community diversity and disturbance intensity in the absence of fluctuations. Adding fluctuations increased community diversity and erased the U-shape. All our results are well-captured by a Monod consumer resource model, which also explains how U-shaped DDRs emerge via a novel ‘niche flip’ mechanism. Broadly, our combined experimental and modeling framework demonstrates how distinct features of an environmental disturbance can interact in complex ways to govern ecosystem assembly and offers strategies for reshaping the composition of microbiomes.

    1. Computational and Systems Biology
    Michael S Lauer, Deepshikha Roychowdhury
    Research Article Updated

    Previous reports have described worsening inequalities of National Institutes of Health (NIH) funding. We analyzed Research Project Grant data through the end of Fiscal Year 2020, confirming worsening inequalities beginning at the time of the NIH budget doubling (1998–2003), while finding that trends in recent years have reversed for both investigators and institutions, but only to a modest degree. We also find that career-stage trends have stabilized, with equivalent proportions of early-, mid-, and late-career investigators funded from 2017 to 2020. The fraction of women among funded PIs continues to increase, but they are still not at parity. Analyses of funding inequalities show that inequalities for investigators, and to a lesser degree for institutions, have consistently been greater within groups (i.e. within groups by career stage, gender, race, and degree) than between groups.