The relative strength and timing of innate immune and CD8 T-cell responses underlie the heterogeneous outcomes of SARS-CoV-2 infection

Mathematical model of within-host SARS-CoV-2 dynamics

We considered an individual infected by SARS-CoV-2. We modeled disease progression in the individual by following the time-evolution of the population of infected cells (I), the strength of the effector CD8 T-cells (E), the cytokine-mediated innate immune response (X), and tissue damage (D) (Figure 1). We considered the essential interactions between these entities (30) and constructed the following equations to describe their time-evolution:

• Download figure
• Open in new tab

Figure 1: Schematic of the mathematical model of within-host SARS-CoV-2 infection.

I represents infected cells, X represents the innate immune response, E the CD8 T-cell response and D the tissue damage. Arrows and blunt-head arrows depict positive and negative regulation, respectively. The parameters and expressions shown next to the arrows are described in the text.

Here, the infected cells follow logistic growth (30), with a per capita growth rate k₁ and carrying capacity I_max. This growth represents the infection of target cells by virions produced by infected cells (30). I_max is the maximum number of cells that can get infected, due to target cell or other limitations. The growth rate k₁ is assumed to be reduced by the innate immune response, X, with the efficacy ε_IX, due to interferon-mediated protection of target cells and/or lowering of viral production from infected cells (10). Effector cell-mediated killing of infected cells is captured by a mass action term with the second-order rate constant k₂. The proliferation and exhaustion of CD8 T-cells are both triggered by infected cells at maximal per capita rates k₃ and k₄, respectively. k_p and k_e are the levels of infected cells at which the proliferation and exhaustion rates are half-maximal, respectively. Following previous studies, we let k₃ < k₄ and k_p < k_e, so that proliferation dominates at low antigen levels and exhaustion at high antigen levels (30–32). The innate response, X, is triggered by infected cells at the per capita rate k₅ and is depleted with the first-order rate constant k₆. To assess the severity of infection, we employ D, which represents the instantaneous tissue damage, with contributions from CD8 T-cell mediated killing of infected cells, determined by αIE, and from proinflammatory cytokines, represented by βX. Inflamed tissue may recover with the first order rate constant γ.

Solving these equations would predict the time-course of the infection. We tested the model by applying it to describe available patient data of viral load changes.

Model fits patient data

A number of studies have reported viral load measurements during the course of SARS-CoV-2 infection (33, 34). In most studies, measurements begin from the time of symptom onset because the time of contracting the disease is rarely known. Because the prodromal period may vary substantially across individuals (35), measurements from symptom onset may miss the initial phases of the immune response, which can be an important determinant of disease outcome. In asymptomatic individuals, this early response clears the infection (36). We therefore sought datasets that included accurate estimates of the time of contracting the disease. Fortunately, we found such data in a study of one of the first SARS-CoV-2 transmission chains in Germany in early 2020 (37, 38). The study traced the dates of first exposure to the virus for each patient in the transmission chain (38) (supplementary text section A; supplementary table S1). Further, daily viral load data, measured in nasopharyngeal swab and sputum samples, for all patients starting from the onset of symptoms or earlier were reported (37). Data from the nasopharyngeal swabs are thought not to be the best correlates of disease outcome and severity (39). The sensitivity of SARS-CoV-2 detection in sputum is substantially higher than in nasopharyngeal swabs (40). We therefore employed data from the sputum samples in this study. We considered data from day zero to day 15 into the infection (supplementary text section A; supplementary tables S1-S3). Beyond two weeks, the humoral response is mounted in most patients (7, 25), the role of which, as mentioned above, is poorly understood (7).

We fit our model to the above viral load data, representing the dynamics of the infection and immune responses in the respiratory tract. All patients in this dataset had mild symptoms, which waned by day 7 after the first virological test. The patients were of working age and otherwise healthy. In such patients, markers of T-cell exhaustion are not significantly higher than healthy individuals and are markedly lower than severely infected patients (15). Therefore, to facilitate more robust parameter estimation, we ignored CD8 T-cell exhaustion in the present fits (by fixing k₄ = 0). Furthermore, we assumed that the viral population, V, is in a pseudo-steady state with the infected cell population, so that V ∝ I. Since, the dynamics of tissue damage (D) is dependent on but does not affect the dynamics of infected cells (I), CD8 T-cells (E) and the cytokine mediated innate response (X), in our model, we ignored D for the present fitting. This is further justified because the patients considered for fitting are mildly/moderately infected, and are expected to suffer minimal tissue damage. Because the patients were all similar, we assumed that I_max would be similar in them and proportional to V_max, the highest viral load reported across the patients. We thus fit log₁₀(I/I_max) calculated with our model to the normalized data of log₁₀(V/V_max). Our fits were not sensitive to I_max (supplementary tables S4, S5). We allowed a delay following exposure to account for the incubation period before viral replication can begin. We used a nonlinear mixed-effects modelling approach for parameter estimation (41). Our model provided good fits to the data (figure 2, first column of subplots) and yielded estimates of the parameters at the population-level (supplementary table S6) and for the individual patients (supplementary table S7).

• Download figure
• Open in new tab

Figure 2: Fits of the mathematical model to patient data.

X-axis represents time from viral exposure in all subplots. The quantity plotted on the Y-axis for all subplots in a given column is mentioned at the top of the column. The first column shows data from patients in open circles (37). These represent normalized viral loads from sputum samples. Best-fit model predictions are shown as black curves. Patient IDs as provided in Böhmer et al. (38) are in the top-left corner of each subplot. The magenta curves in the second column, the blue area plots in the third column and the green area plots in the fourth column represent the corresponding dynamics of infected cells, cytokine mediated innate immune response and effector CD8 T-cell mediated adaptive immune response, respectively. Parameter values used are listed in supplementary table S7.

To ascertain the robustness of our model and fits, we tested several variants of our model. We fit variants without the adaptive response, without the innate response, with a logistic growth formulation of the innate immune response, with the innate response amplifying the adaptive response, or combinations thereof to the same data (supplementary text section B, supplementary table S8). The fits were all poorer than the present model (figure 1, supplementary table S8). We thus employed our present model for further analysis.

Model elucidates distinct patterns of viral clearance and associated immune responses

The best-fits above yielded important insights into the underlying dynamics of disease progression and clearance. First, the time between viral exposure and noticeable escalation of the viral load, i.e., the post-exposure delay in viral replication, varied from 0.8 d to 6.6 d in the patients analyzed, with a mean of 2.7±0.8 days, reflecting the variability in the time of the establishment of systemic infection following exposure, and consistent with the variable prodromal period observed (35). The initial, possibly stochastic (42), events during the establishment of infection might be associated with the variability in the delay in viral replication. Second, our model offered an explanation of the two distinct patterns of clearance observed in the patients. Patients 1, 2, 3, 4, and 14 had a single peak in viral load (or infected cell numbers) followed by a decline of viral load leading to clearance (figure 2, second column of subplots). Patients 7, 8 and 10, in contrast, had a second peak following the first. Our model predicted these distinct patterns as arising from the temporal variation in the dynamics of the innate and CD8 T-cell responses.

The interactions between the innate response, X, and infected cells, I, in our model have signatures of the classic predator-prey system (43) with I the prey and X the predator: In the absence of X, I grows. I also triggers X, which in turn suppresses I. X declines in the absence of I. These interactions, as with the predator-prey system (43), predict oscillatory dynamics. Thus, following infection, I grows, causing a rise of X in its wake. When X rises sufficiently, it suppresses I. When I declines substantially, the production of X is diminished and X declines. This allows I to rise again and the cycle repeats. This cycle is broken in our model by CD8 T-cells, E. Viral clearance is not possible in our model without E (supplementary figure S1). When E rises, it can suppress I independently of X, breaking the cycle and allowing X to dominate I. Together, X and E can then clear the infection. In patients 7, 8, and 10, our best-fits predicted an early innate immune response and a delayed CD8 T-cell response. The second peak was thus predicted as a result of the above predator-prey oscillations that occurred before the CD8 T-cell response was mounted. In patients 3, 4, and 14, a relatively early CD8 T-cell response was predicted, which precluded the second peak. In patients 1 and 2, both the innate and CD8 T-cell responses were delayed, leaving little time for the oscillations to arise in the 15 d period of our observations. We note that interpretations of the multiple peaks in longitudinal viral load data have not been forthcoming (44). Our predictions offer a plausible interpretation.

Third, the transient but robust innate immune response predicted (figure 2, third column of subplots) is consistent with observations in mildly infected patients (45). Fourth, our prediction of the dynamics of the CD8 T-cell response, where a gradual build-up is followed by a stationary phase (figure 2, fourth column) is also consistent with observations. In mildly infected patients, SARS-CoV-2 specific T-cells were detected as early as 2-5 days post symptom onset (12). This effector population remained stable or increased for several months after clinical recovery (16, 46).

Our model thus fit the dynamics of infection in individuals showing mild symptoms and offered explanations of disease progression patterns that had remained confounding. We examined next whether the model could also describe more severely infected patients. For this, we varied different parameters in our model and assessed the resulting dynamical features of the infection.

Interplay between innate and CD8 T-cell responses underlies heterogeneous disease outcomes

We reintroduced the CD8 T-cell exhaustion term, which we had ignored in the fits above because the patients were mildly infected, and selected associated parameter estimates from previous studies (30). We ensured that this did not affect our fits above (supplementary figure S2).

Next, to estimate the severity of the disease, we also examined the dynamics of the instantaneous tissue damage (D). Typically, D rose as the infection progressed and declined as it got resolved (supplementary figure S3A, C). We reasoned that the severity of infection would be determined by the maximum tissue damage suffered and the duration for which such damage lasted. Significant damage that is short-lived or minimal damage that is long-lived may both be tolerable and lead to mild symptoms. We therefore calculated the area under the curve (AUC) of D, starting from when D ascended above its half-maximal level to the time when it descended below that level (supplementary figure S3A), as a measure of immunopathology (P) and associated disease severity. (Note that the parameters α, β, and γ, which describe the dynamics of tissue damage, are unknown constants; our results were not sensitive to their values because changing them only minimally affected the relative extents of immunopathology across different disease severity categories (supplementary text section C; supplementary figure S5).)

With this framework, we varied the strengths of the CD8 T-cell and innate responses, by changing the values of the parameters k₃ and k₅, respectively, and examined the predicted dynamical features (figure 3A). Recall that k₃ is the proliferation rate of CD8 T-cells and k₅ is the rate of generation of the innate immune response. The other parameters were fixed (supplementary table S6) at their population estimates, for which the model elicited clearance of the infection by day 14 (figure 3A, center, subplot with an arrowhead). Increasing k₅ resulted in a decrease in the peak of infected cells (figure 3A, the row of subplots with arrowhead, right to left). With decrease in k₅, the induction of the cytokine mediated antiviral innate response was substantially delayed and that corresponded to an increased number of infected cells (supplementary figure S4A, B). Clearance was achieved in all cases without substantial variation in the infection duration because of the CD8 T-cell response (figure 3A, the row of subplots with arrowhead, right to left). Decreasing k₃ weakened the CD8 T-cell response and increased the duration of the infection (figure 3A, the column of subplots with arrowhead, bottom to top). We next explored the effects of varying both k₃ and k₅ simultaneously.

• Download figure
• Open in new tab

Figure 3: Variations in innate and CD8 T-cell responses capture disease heterogeneity.

(A) Effect of simultaneous variation of parameters determining the strengths of innate and CD8 T-cell responses. The black, annotated triangles at the top and right indicate the nature and the direction of the variation of the indicated parameters. For instance, k₅ is decreased from left to right. Individual subplots show the dynamics of infected cells normalized by carrying capacity, cytokine mediated innate immune response and effector CD8 T-cell mediated adaptive immune response. The legends are provided at the top-right. The left Y-axis shows the normalized infected cell dynamics. The right Y-axis shows the other two species, i.e., cytokine-mediated innate immune response and effector CD8 T-cell-mediated adaptive immune response. The rectangular, colored patch at the top of each subplot represents the extent of immunopathology. The range of immunopathology is given by the color scale at the bottom of the figure (below 3C). At the left of the color scale, a separate legend denotes the texture used for depicting diverged immunopathology. The arrowhead on the scale indicates the immunopathology for the population parameters. The central subplot, which also contains an arrowhead, charts the simulation using the population parameters (supplementary table S6). (B, C) Similar calculations corresponding to variations in other parameters associated with the innate (B) and CD8 T-cell (C) responses. Plots with the population parameters are marked with the arrowhead. The colored patches should be interpreted using the color scale provided at the bottom.

When k₃ was high, i.e., the response of CD8 T-cells was strong, irrespective of the innate immune response, the infection was cleared within ten days (figure 3A, subplots on bottom-left and bottom-right). Associated immunopathology was nominal. These predictions were akin to asymptomatic and mild infection scenarios. An early and robust effector T-cell response has been associated with milder infections (12, 16, 46). Here, in some cases with high k₅, a blip of the viral load was observed after an initial phase of clearance. When k₅ was decreased, marking a weaker innate response, the peak viral load rose and immunopathology moderately increased. This was also observed when we decreased ε_I, which lowered the efficacy with which the innate immune response inhibits the spread of the infection (figure 3B). The latter trends associated with high k₃ and low k₅ have parallels to infected patients with robust CD8 T-cell responses but impaired innate responses, such as those harboring mutations in the genes associated with the activation of the antiviral resistance in host cells (47). Clearance was achieved in such cases due to the robust CD8 T-cell response.

When k₃ was low and k₅ was high (figure 3A, four subplots on top-left), the infection was prolonged. However, the immunopathology was lower than when both k₃ and k₅ were high. The efficient antiviral innate response controlled the initial peak of the infection. However, the slow proliferation of the effector cells delayed clearance. This scenario had parallels to the reported cases of prolonged RT-PCR positivity of viral loads (22–24). Restrained CD8 T-cell differentiation was associated with such cases (22). Delayed clearance was also realized when the parameter k_p was increased, which increased the antigen level required for significant effector T-cell proliferation (figure 3C). These predictions were consistent with observations of defects in T-cell proliferation delaying the clearance of infection (22).

When both k₃ and k₅ were low (figure 3A, four subplots on top-right), severe immunopathology along with prolonged infection with high viral load and high cytokine levels was predicted. When they were the lowest, clearance was not achieved in our simulations. Although, clearance of the infection is the predominant outcome associated with a wide range of parameter values (figure 3), our dynamical systems analysis revealed that in certain parameter regimes clearance may not result (supplementary text section D, supplementary figures S6, S7). Instead, escape from immune protection with a high level of infected cells and cytokines together with a high degree of CD8 T-cell exhaustion may occur. Such runaway trajectories were associated with high immunopathology in our model (figure 3A, top right corner, supplementary figure S3B, S3C, top right corner) and were predicted to be terminated by fatality. These trends in the model mirrored clinical features of severe COVID-19 (45), which include consistently very high viral loads, heightened proinflammatory cytokines and interferons (39, 45, 48), attenuated proliferation (13) and increased exhaustion of T-cells (13, 14, 17). The predictions also include cases where in the late phase of the infection, although the viral load in sputum shows a decline (49), and that in nasopharyngeal swab becomes undetectable (33, 39, 49), mortality results due to intolerable immunopathology.

Note that the initial pool of CD8 T-cells, E₀, was important in determining outcomes (supplementary text section D, supplementary figure S7), with a large pool leading to rapid clearance, in agreement with observations of such clearance facilitated by cross-reactive effector T cells (12, 50). The outcomes were less sensitive to the viral inoculum size (supplementary text, section E, supplementary figure S8), i.e., I₀, consistent with studies on macaques where different inoculum sizes led to comparable disease outcomes (51).

Model can recapitulate clinical data of varying disease severity across patients

In a recent study, patients were stratified by disease outcome and measurements of longitudinal viral load from their saliva were fit using cubic splines, yielding ribbons of confidence intervals on viral loads for each category (39). We digitized these ribbons and tested our model predictions against them (figure 4A-D, grey patches). The study reported data from symptom onset. We therefore added an estimated length of the prodromal period to the timepoints in order to compare our model predictions. We set this length to 4.8 d from the German transmission chain data (37, 38) (supplementary table S2), which is also consistent with other reported estimates (33). We estimated the viral load, V, from our model predictions using the pseudo-steady state approximation, V ≈ pI/c, where p is the per capita rate of viral production from infected cells and c is the per capita rate of viral clearance. We set p and c to values estimated previously (52) (supplementary text section A).

• Download figure
• Open in new tab

Figure 4: Model recapitulates viral load data from patients with different severities of infection.

(A) The grey patch, as provided in Silva et al. (39), represents the confidence interval of a cubic spline fitted to the viral load data from non-hospitalized individuals. Curves depict simulated viral load trajectories in individuals representative of those in the patch. The scale on the right represents immunopathology. The immunopathology associated with the simulated trajectories are marked on the scale using arrowheads having the same colors as the profiles. (B-D) Simulations recapitulating the viral load trajectories from moderate, severe and deceased patients, respectively. Immunopathology of deceased individuals not shown were higher than the upper limit of the scale shown here. The blue bar on the scale in (D) represents the range of immunopathology beyond which fatality is the likely clinical outcome. The parameters used for each trajectory are in supplementary tables S9-S12. Those associated with the profiles with dashed lines were used for sensitivity analyses (supplementary figures S9-S12).

The ribbons of data are not amenable to fitting. We therefore varied parameters associated with the innate and CD8 T-cell responses in our model to achieve dynamical profiles of viral load resembling the ribbons. For the purpose of these simulations we ignored both the potential adverse and favorable effects of antibodies. The ribbon for the non-hospitalized patients was associated with low viral loads (figure 4A). The peak viral load was approximately 10⁶ copies of viral RNA. This relatively low peak viral load could be captured by our model when the strength of the innate response was increased, which we achieved by increasing k₅ and/or ε_I (supplementary table S9). The duration of the infection was dependent on the CD8 T-cell response. Strong CD8 T-cell stimulation, achieved with a low value of k_p, led to rapid clearance, whereas weaker stimulation, corresponding to a higher k_p, allowed the infection to remain for an extended period. We calculated the immunopathology associated with these simulations and found it to be low (figure 4A, scale on the right; also see below).

The ribbon for patients eliciting moderate symptoms had a relatively higher viral load at the peak, reaching approximately 10⁸ copies of viral RNA (figure 4B). Parameters sufficiently close to the population estimates above allowed us to capture the dynamics for these patients (supplementary table S10). The associated immunopathology was considerably higher than the non-hospitalized patients. For severely infected patients, the viral load peak was above 10⁸, reaching as high as 10¹⁰ copies (figure 4C). We achieved this high peak viral load by lowering the strength of the innate response (decreasing k₅ and/or ε_I; see supplementary table S11). The delayed clearance could be recapitulated by lowering CD8 T-cell stimulation (increasing k_p and/or decreasing k₃). The immunopathology was higher than those calculated to capture the viral load dynamics in moderate patients (compare the scales in figure 4B and 4C). Lastly, for the deceased individuals, the peak viral load was similar to the severe patients. However, the downward trend in the viral load after the peak seen with the severely infected patients was no longer apparent (figure 4D). The viral load remained around 10⁸ RNA copies till day 30 post-exposure. A much weaker innate response (low ε_I) and a weaker CD8 T-cell response (high k_p) could generate matching profiles (supplementary table S12). The immunopathology for the deceased patients was consistently higher than the severe patients, indicating that there might be an upper limit to the extent of immunopathological tissue damage that lay somewhere between our estimates for severe and deceased patients, and crossing which mortality would almost certainly result.

Our model thus recapitulated the trends in the viral load seen in patients with different severity of infection. Furthermore, the model indicated that there should be a narrow range of immunopathology, which acts as a threshold to determine the fatal outcomes in COVID-19 (figure 4D, scale on the right).

To assess whether variations in other model parameters could also achieve the above trends, we performed a global sensitivity analysis of our parameters using one representative parameter combination from each disease outcome category as reference (figure 4A-D, dashed curves; supplementary figures S9-S12; supplementary tables S9-S12). Specifically, we calculated how sensitive our measure of immunopathology was to the parameters. We found that immunopathology was most sensitive to ε_I and k₆ for non-hospitalized and moderately symptomatic patients (supplementary figures S9, S10). For severely infected and deceased patients, I_max, k₃ and ε_I emerged as the important parameters (supplementary figures S11, S12). These results reinforce our expectations above. In mild infections, the innate immune response is strong and any variation in its strength would have the most influence on disease severity. In more severe infections, the innate and adaptive responses are both involved and the severity is therefore sensitive to variations in the strengths of both.

Model predictions thus successfully recapitulated the heterogeneous outcomes and the associated dynamical patterns of SARS-CoV-2 infection.

Study design

We constructed a mathematical model of within-host SARS-CoV-2 infection using ordinary differential equations. Next, we utilized a nonlinear mixed-effects approach to fit the model to an available clinical dataset and estimated model parameters (supplementary text, section A, B). The model was then utilized for exploring the effects of parameter variations (figure 3), recapitulating the viral load trajectories in patients stratified by disease severity (figure 4), and for sensitivity analysis (supplementary figures S9-S12).

Model construction

The equations of our model are provided in the results section. The other models which were fit to the data are described in the supplementary text (supplementary text, section B).

Parameter estimation and model selection

Published data from Wölfel et al.(37) was digitized by a custom script in the MATLAB (version R2020a) image analysis toolbox (www.mathworks.com). This dataset was further used for fitting different models using the stochastic approximation expectation maximization (SAEM) algorithm available in Monolix 2020R1 (www.lixoft.com) (supplementary text, section A). The Akaike information index (AIC) was calculated within the Monolix environment. The model with the lowest AIC was selected for further mathematical analysis (supplementary text, section B).

Fixed points and linear stability analysis

For the steady-state analysis, estimated parameter values were utilized. MATLAB (version R2020a) was used to estimate the fixed points of the system and to determine the nature of their stability. Individual fixed points and their corresponding Jacobian matrices were estimated using the Symbolic Math Toolbox (www.mathworks.com). Calculation of the eigenvalues and eigenvectors for individual fixed points yielded the nature of their stability and facilitated determination of the phase portraits (supplementary text, section D, supplementary figures S6, S7).

Recapitulating patient viral load data stratified by disease severity

The published data from Silva et al. (39), was digitized using a custom MATLAB code, using functions from the image analysis toolbox. The parameters were manually varied and the quality of the fits determined by visual inspection of the simulation profile and the confidence interval ribbons.

Sensitivity analysis

We executed variance based global sensitivity analysis (VBGSA) on the models; the details of the algorithm have been described elsewhere (94). We simultaneously varied the parameters up to 5% above and below the population parameters in Monte Carlo simulations and calculated total effect indices for the parameters (supplementary figures S9-S12).