The model of primary viral infection is a coupled system of ordinary and delay differential equations, consisting of three major components (see Figure 1 for a schematic diagram). Equations (1)–(3) describe the process of infection of target cells by influenza virus and are a major component in almost all models of viral dynamics in the literature. Equations (4) and (5) model IFN-mediated innate immunity (45, 46). Third, adaptive immunity including CD8⁺ T cells and B cell-produced antibodies for killing infected cells and neutralizing influenza virus, respectively, are described by equations (6)–(11).

In further detail, equation (1) indicates that the change in viral load (V) is controlled by four factors: the production term (p_VI) in which virions are produced by infected cells (I) at a rate of p_V (37, 45, 47); the viral natural decay/clearance (δ_VV) with a decay rate of δ_V; the viral neutralization terms (κ_SVA_S and κ_LVA_L) by antibodies (both a short-lived antibody response A_S driven by, e.g., IgM, and a longer-lived antibody response A_L driven by, e.g., IgG and IgA (12, 38)), and a consumption term (βVT) due to binding to and infection of target cells (T). In equation (2), the term g_T(T + R)(1 − (T + R + I)/T₀) models logistic regrowth of the target cell pool (46). Both target cells (T) and resistant cells (R, those protected due to IFN-induced antiviral effect) can produce new target cells, with a net growth rate proportional to the severity of infection, 1 − (T + R + I)/T₀ (i.e., the fraction of dead cells). T₀ is the initial number of target cells and the maximum value for the target cell pool (34). Target cells (T) are consumed by virus (V) due to binding (β′VT), the same process as βVT. Note that β and β′ have different measurement units due to different units for viral load (V) and infected cells (I). As mentioned earlier, the innate response may trigger target cells (T) to become resistant (R) to virus at a rate ϕFT. Resistant cells lose protection at a rate ρ (45). This process also governs the evolution of virus-resistant cells (R) in equation (5).

Equation (3) describes the change of infected cells (I). They increase due to the infection of target cells by virus (β′VT) and die at a (basal) rate δ_I. Two components of the immune response increase the rate of killing of infected cells. IFN-activated NK cells kill infected cells at a rate κ_NIF (6, 45, 46, 48). Effector CD8⁺ T cells (E)—produced through differentiation from naive CD8⁺ T cells C_n in equation (6)—kill at a rate κ_EIE. Of note, our previous work has suggested that models of the innate response containing only IFN-induced resistance for target cells (state R; equation (5)), while able to maintain a population of healthy uninfected cells, still control viral kinetics through target cell depletion. While it remains possible that these models may be able to reproduce features of the viral reexposure data (43, 46), that work also demonstrated that inclusion of IFN-activated NK cells (term κ_NIF) provides a natural explanation for the observed viral reexposure data.

Equation (4) models the innate response, as mediated by IFN (F). IFN is produced by infected cells at a rate p_F and decays at a rate δ_F (46).

Equation (6) models stimulation of naive CD8⁺ T cells (C_n) into the proliferation/differentiation process by virus at a rate β_CnV/(V + h_C), where β_Cn is the maximum stimulation rate and h_C indicates the viral load (titV) at which half of the stimulation rate is achieved. Note that this formulation does not capture the process of antigen presentation and CD8⁺ T cell activation, but rather is a simple way to establish the essential coupling between the viral load and the rate of CD8⁺ T cell activation in the model (49). In equation (7), the production of effector CD8⁺ T cells (E) is assumed to be an “advection flux” induced by a delayed virus stimulation of naive CD8⁺ T cells [the first term on the right-hand side of equation (7)]. The delayed variables, V (t − τ_C) and C_n(t − τ_C), equal zero when t < τ_C. The introduction of the delay τ_C is to phenomenologically model the delay induced by both naive CD8⁺ T cell proliferation/differentiation and effector CD8⁺ T cell migration and localization to the site of infection for antiviral action (42, 50, 51). The delay also captures the experimental finding that naive CD8⁺ T cells continue to differentiate into effector T cells in the absence of ongoing antigenic stimulation (49, 52). The multiplication factor epCτC indicates the number of effector CD8⁺ T cells produced from one naive CD8⁺ T cell, where p_C is the average of effector CD8⁺ T cell production rate over the delay period τ_C. The exponential form of the multiplication factor is derived based on the assumption that cell differentiation and proliferation follow a first-order advection–reaction equation. Effector CD8⁺ T cells decay at a rate δ_E.

Similar to CD8⁺ T cells, equations (8) and (9) model the proliferation/differentiation of naive B cells, stimulated by virus presentation at rate β_BnV/(V + h_B). Stimulation subsequently leads to the production of plasma B cells (P) after a delay τ_B. The multiplication factor epBτB indicates the number of plasma B cells produced from one naive B cell, where p_B is the production rate. Plasma B cells secrete antibodies, which exhibit two types of profiles in terms of experimental observation: a short-lived profile (e.g., IgM lasting from about day 5 to day 20 postinfection) and a longer lived profile (e.g., IgG and IgA lasting weeks to months) (12, 38). These two antibody responses are modeled by equations (10) and (11), wherein different rates of production (P_S and P_L) and consumption (δ_S and δ_L) are assumed.

The model contains 11 equations and 30 parameters (see Table 1). This represents a serious challenge in terms of parameter estimation and clearly prevents a straightforward application of standard statistical techniques. To reduce uncertainty, a number of parameters were taken directly from the literature, as per the citations in Table 1. The rest were estimated (as indicated in Table 1) by calibrating the model against the published data from Miao et al. (38) who measured viral titer, CD8⁺ T cell counts, and IgM and IgG antibodies in laboratory mice (exhibiting a full immune response) over time during primary H3N2 A/Hong Kong/X31 influenza A virus infection (see Ref. (38) for a detailed description of the experiment). The approach to estimating the parameters based on Miao et al.’s data is provided in Supplementary Material, and the estimated parameter values are given in Table 1. Note that the data were presented in scatter plots in the original paper (38), while we presented the data here in mean ± SD at each data collection time point for a direct comparison with our mean-field mathematical model.

For model simulation, the initial condition is set to be (V, T, I, F, R, C_n, E, B_n, P, A_S, A_L) = (V ₀, T₀, 0, 0, 0, 100, 0, 100, 0, 0, 0) unless otherwise specified. The initial target cell number (T₀) was estimated by Petrie et al. (53). We estimate that of order 100 cells (resident in the spleen) are able to respond to viral infection (C_n) (personal communication, N. LaGruta, Monash University, Australia). Note that 100 naive CD8⁺ T cells might underestimate the actual number of naive precursors that could respond to all the epitopes contained within the virus but does not qualitatively alter the model dynamics and predictions (see Section 3 where the naive CD8⁺ T cell number is varied between 0 and 200, where the upper bound is sufficient to show the model’s full range of behaviors). In the absence of further data, we also use this value for the initial naive B cell number (B_n), but again this choice does not qualitatively alter the model predictions. The numerical method and code (implemented in MATLAB, version R2014b, the MathWorks, Natick, MA, USA) for solving the model are provided in Supplementary Material.

Clinical influenza A (H7N9) patient data were used to test our model predictions on the relationship between CD8⁺ T cell number and recovery time. The data were collected from 12 surviving patients infected with H7N9 virus during the first wave of infection in China in 2013. (Raw data are provided in Data Sheet S1 in Supplementary Material; see the paper of Wang et al. (21) for details of data collection; this study was reviewed and approved by the SHAPHC Ethics Committee.) Note that the clinical data were scarce for some patients. For those patients, we have assumed that the available data are representative of the unobserved values in the neighboring time period. For each patient, we took the average IFNγ⁺ CD8⁺ T cell number in 10⁶ peripheral blood mononuclear cells (PBMC) for the period from day 8 to day 22 (or the recovery day if it comes earlier) postadmission as a measure of the effector CD8⁺ T cell level. This period was chosen a priori as it roughly matches the duration of the CD8⁺ T cell profile, and clinical samples were frequently collected in this period. The average CD8⁺ T cell count was given by the ratio of the total area under the data points (using trapezoidal integration) to the number of days from day 8 to day 22 (or the recovery day if it comes earlier). For those patients for whom samples at days 8 and/or 22 were missing, we specified the average CD8⁺ T cell level at the missing time point to be equal to the value from the nearest sampled time available.

We first analyze the model behavior to establish a clear understanding of the model dynamics. Figure 2 shows solutions (time series) for the model compartments (viral load, CD8⁺ T cells, and IgM and IgG antibody) calibrated against the murine data from the study by Miao et al. (38). Solutions for the remaining model compartments are shown in Figure 3. The model (with both innate and adaptive components active) prevents the depletion of target cells (see Figure 3 wherein over 50% of target cells remain during infection) and results in a minor loss of just 10–20% of healthy epithelial cells (i.e., the sum of target cells (T) and virus-resistant cells (R); see Figure S1 in Supplementary Material). We note the important difference between prevention of target cell depletion on the one hand and maintenance of healthy cells on the other hand. To be compatible with heterologous reexposure studies (43, 46), a model must not only maintain the population of healthy cells (as many of the aforementioned models do) but must also prevent depletion of target cells to enable infection on rapid reexposure. Otherwise, if T is driven low and R high, while the healthy cell population will be maintained, infection on reexposure may still be blocked. In our model, the primary driver for the maintenance of the target cell pool during acute viral infection is a timely activation of the innate immune response (Figure S2 in Supplementary Material), indicating that our model improves upon previous models where viral clearance was only achieved through depletion of target cells (a typical solution shown in Figure S2B in Supplementary Material).

The modeled viral dynamics exhibits three phases, each dominated by the involvement of different elements of the immune responses (Figure 4). Immediately following infection (0–2 days postinfection) and prior to the activation of the innate (and adaptive) immune responses, virus undergoes a rapid exponential growth (Figure 4A). In the second phase (2–5 days postinfection), the innate immune response successfully limits viral growth (Figure 4A). In the third phase (4–6 days postinfection), adaptive immunity (antibodies and CD8⁺ T cells) is activated and viral load decreases rapidly, achieving clearance. Figures 4B,C demonstrate the dominance of the different immune mechanisms at different phases. In Figure 4B, models with and without immunity are indistinguishable until day 2 (shaded region), before diverging dramatically when the innate and then adaptive immune responses influence the dynamics. In Figure 4C, models with and without an adaptive response only diverge at around day 4 as the adaptive response becomes active. We have further shown that this three-phase property is a robust feature of the model, emergent from its mathematical structure and not a property of fine tuning of parameters (see Figure S3 in Supplementary Material). Importantly, it clearly dissects the periods and effect of innate immunity, extending on previous studies of viral infection phases where the innate immune response was either ambiguous or ignored (12, 38, 55).

As reviewed by Dobrovolny et al. (39), a number of in vivo studies have been performed to dissect the contributions of CD8⁺ T cells and antibodies (12, 17, 18, 44, 56). We use the findings of these studies to validate our model, by testing how well it is able to reproduce the experimental findings (without any further adjustment to parameters). Although the determination of the role of CD8⁺ T cells is often hindered by co-inhibition of both CD8⁺ T cells and the long-lived antibody response (e.g., using nude mice), it is consistently observed that antibodies play a dominant role in final viral clearance, while CD8⁺ T cells are primarily responsible for the timely killing of infected cells and so indirectly contribute to an increased rate of removal of free virus toward the end of infection (13, 17, 18). Furthermore, experimental data demonstrate that a long-lived antibody response is crucial for achieving complete viral clearance, while short-lived antibodies are only capable of driving a transient decrease in viral load (12, 44). We find that our model (with parameters calibrated against Miao et al.’s data (38)) is able to reproduce these observations:

Virus can rebound in the absence of long-lived antibody response (see Figure 5; Figure S4 in Supplementary Material).

In addition, the model also predicts a rapid depletion of naive CD8⁺ T cells after primary infection (see Figure 3), which represents a full recruitment of naive CD8⁺ T cell precursors. This result may be associated with the experimental evidence suggesting a strong correlation between the naive CD8⁺ T cell precursor frequencies and effector CD8⁺ T cell magnitudes for different pMHC-specific T cell populations (57). Note that in Figure 5, no adjustments to the model (e.g., to the vertical scale) were made; its behavior is completely determined by the calibration to the aforementioned murine data (38), and so these findings represent a (successful) prediction of the model.

In summary, we have demonstrated that our model—with parameters calibrated against murine data (38)—exhibits three important phases characterized by the involvement of various immune responses. Advancing on previous models, our model does not rely on target cell depletion and successfully reproduces a multitude of behavior from knockout experiments where particular components of the adaptive immune response were removed. This provides us with a well-tested platform, in which all major components of the immune response have been included and tested, with which to now make predictions on the effect of the cellular adaptive response on viral clearance.

Having established that our model is (from a structural point of view) biologically plausible and that our parameterization is capable of reproducing varied experimental data under different immune conditions (i.e., knockout experiments), we now study how the cellular adaptive response influences viral kinetics in detail. We focus on the key clinical outcome of recovery time, defined in the model as the time when viral titer first falls below 1 EID₅₀/ml, the minimum value detected in relevant experiments (e.g., Figure 2).

Time series of the viral load show that the recovery time decreases as the initial naive CD8⁺ T cell number (C_n) increases (Figure S4 in Supplementary Material). With that in mind, we now examine how recovery time is associated with the clinically relevant measure of effector CD8⁺ T cell level during viral infection. With an increasing initial level of naive CD8⁺ T cells, the average level of effector CD8⁺ T cells over days 6–20 increases linearly (Figure 6A), while the recovery time decreases in an approximately exponential manner (Figure 6B). Combining these two effects gives rise to an approximately exponential relation between the level of effector CD8⁺ T cells and recovery time (Figure 6C). Note that the exponential/linear fits shown in the figures are simply to aid in interpretation of the results. They are not generated by the viral dynamics model.

If varying the delay for naive CD8⁺ T cell activation and differentiation, τ_C, while keeping the naive CD8⁺ T cell number fixed (at the default value of 100), we find that the average level of effector CD8⁺ T cells is exponentially related to the delay, while the recovery time is dependent on the delay in a piecewise linear manner (see Figure S5 in Supplementary Material). Nevertheless, the combination still leads to an approximately exponential relationship between the level of effector CD8⁺ T cells and recovery time (Figure S5C in Supplementary Material), which is almost identical to that of varying naive CD8⁺ T cells (Figure 6D). We also examine the sensitivity of the exponential relationship to other model parameters generally accepted to be important in influencing the major components of the system, such as the viral production rate p_V, IFN production rate p_F, and naive B cell number. We find that the exponential relationship is robust to significant variation in all of these parameters (see Figures S6 and S7 in Supplementary Material; the result for varying naive B cell number is shown in the last section of Section 3). These results suggest that a higher level of effector CD8⁺ T cells is critical for early recovery, consistent with experimental findings (58).

Finally, and perhaps surprisingly given our model has been calibrated purely on data from the mouse, a strikingly similar relationship as shown in Figure 6C is found in clinical data from influenza A (H7N9) virus-infected patients (Figure 7). Excluding one patient (No. a79 in Data Sheet S1 in Supplementary Material; the exclusion is considered further in Section 4), average IFNγ⁺ CD8⁺ T cells and recovery time are negatively correlated (Spearman’s ρ = − 0.8368, p = 0.0013) and well captured by an exponential fit with an estimated offset (see Figure 7 caption for details). The exponential relationship (observed in both model and data) has features of a rapid decay for relatively low/intermediate levels of effector CD8⁺ T cells and a strong saturation for relatively high CD8⁺ T levels, implying that even with a very high level of naive CD8⁺ T cells, recovery time cannot be reduced below a certain value (in this case, estimated to be approximately 17 days). Of note, the exponential relationship (i.e., the scale of CD8⁺ T cell level or recovery time) is only a qualitative one, as we have no way to determine the scaling between different x-axis measurement units, nor adjust for particular host and/or viral factors that differ between the two experiments (i.e., H3N2-infection in the mouse (38) versus H7N9 infection in humans (21)).

In addition to naive CD8⁺ T cells, memory CD8⁺ T cells (established through previous viral infection) may also significantly affect recovery time due to both their rapid activation on antigen stimulus and faster replication rate (54, 59–61). To study the role of memory CD8⁺ T cells, we must first extend our model. As we are only concerned with how the presence of memory CD8⁺T cells influences the dynamics, as opposed to the development of the memory response itself, the model is modified in a straightforward manner through addition of two additional equations that describe memory CD8⁺ T cell (C_m) proliferation/differentiation: (12)dCmdt=−βCmVV+hCmCm, (13)dEmdt=βCmV(t−τCm)V(t−τCm)+hCmCm(t−τCm)e(pCmτCm)−δEEm.

Accordingly, the term κ_EIE in equation (3) is modified to κ_EI(E + E_m). The full model and details on the choice of the additional parameters are provided in Supplementary Material. Note that the model component, C_m, may include different populations of memory CD8⁺ T cells, including those directly specific to the virus and those stimulated by a different virus but which provide cross-protection (62, 63).

Figure 8A shows how the preexisting memory CD8⁺ T cell number (C_m) changes the exponential relationship between naive CD8⁺ T cells and recovery time. Importantly, as the number of memory CD8⁺ T cells increases, the recovery time decreases for any level of naive CD8⁺ T cells and the exponential relationship remains. For patients with a relatively low level of naive CD8⁺ T cells (i.e., on the left of Figure 8A), the extent of reduction in the recovery time is greater than that for patients with a relatively high level of naive CD8⁺ T cells (i.e., on the right of Figure 8A). This suggests that people with a lower level of naive CD8⁺ T cells may benefit more through induction of memory CD8⁺ T cells, emphasizing the potential importance for taking prior population immunity into consideration when designing CD8⁺ T cell-based vaccines (64).

The above result is based on the assumption that the initial memory CD8⁺ T cell number upon reinfection is independent of the number of naive CD8⁺ T cells available during the previous infection. However, it has also been found that the stationary level of memory CD8⁺ T cells is usually maintained at about 5–10% of the maximum antigen-specific CD8⁺ T cell number during primary viral infection (8, 65). This indicates that people with a low naive CD8⁺ T cell number may also develop a low level of memory CD8⁺ T cells following infection. In consequence, such individuals may be relatively more susceptible to viral reinfection (66). This alternative and arguably more realistic relationship between the numbers of naive and memory CD8⁺ T cells is simulated in Figure 8B where memory CD8⁺ T cell levels are set to 5% of the maximum of the effector CD8⁺ T cell level. Results suggest that, on viral reinfection, preexisting memory CD8⁺ T cells are able to significantly improve recovery time except for hosts with a very low level of naive CD8⁺ T cells (Figure 8B). This is in accordance with the assumption that a smaller naive pool leads to a smaller memory pool and in turn a weaker shortening in recovery time. Although the model suggests that the failure of memory CD8⁺ T cells to protect the host is unlikely to be observed (because of the approximately 30-fold increase in the size of the memory pool relative to the naive pool), the failure range may be increased if the memory pool size is much smaller (modulated by, say, changing 5 to 1% in the model). Therefore, for people with a low naive CD8⁺ T cell number, the level of memory CD8⁺ T cells may be insufficient and prior infection may provide very limited benefit, further emphasizing the opportunity for novel vaccines that are able to induce a strong memory CD8⁺ T cell response to improve clinical outcomes.

Antibodies appear at a similar time as effector CD8⁺ T cells during influenza viral infection and may enhance the reduction in the recovery time in addition to CD8⁺ T cells. By varying the naive B cell number B_n (as a convenient, but by no means unique, way to influence antibody level), we find that increasing the antibody level shortens the recovery time regardless of the initial naive CD8⁺ T cell number, leaving the exponential relation largely intact (Figure 9). A slight saturation occurs for the case in which levels of both naive B cells and CD8⁺ T cells are low. Moreover, variation in naive B cell number also results in a wider variation in recovery time for a lower naive CD8⁺ T cell level, suggesting that people with a lower level of naive CD8⁺ T cells may, once again, receive a more significant benefit (in terms of recovery time) through effective induction of an antibody response via vaccination.