Modelling upper respiratory viral load dynamics of SARS-CoV-2

Fitting a regression model to the viral load data

We fitted two types of models to the viral load data from the first 15 days from onset of symptoms (see Methods), using only subjects with at least three samples during this time period (models fitted to 870 samples from 155 patients from 16 studies). The first was a linear regression model, fitted to log-transformed viral loads. This model included patient- and study-specific random effects, which captured variation from the average behaviour observed across 16 datasets. There was variation in the peak and slope of the viral loads across different studies (Figure 2, Supplementary Table 3). The variation in the (log-transformed) peak viral load varied over several orders of magnitude. Some of this variation was explained by study-specific differences: the standard deviation of the between-study variation in the peak viral load was 0.84 i.e. nearly one order of magnitude (Supplementary Table 3). This could be due to a number of factors, such as the method of sample collection, quantification method or characteristics of included patients. Furthermore, viral loads in several studies were estimated using an averaged standard curve, which introduces some uncertainty into the magnitude of the viral load. However, the inclusion of study-specific random effects allows such data to be appraised alongside data from other studies.

• Download figure
• Open in new tab

Figure 2: Viral load data and mixed-effects regression model.

Data from all 16 studies used in the regression modelling (numbered as in Table 1), showing samples taken within the first 15 days of symptom onset. We fitted a regression model to the data, with study-specific random-effects for the peak viral load and rate of decline (slope). The solid lines show the posterior mean behaviour for each study, with the shaded areas showing the 95% credible intervals. The dashed line, which is the same in each panel, is the average trajectory across the 16 studies. The 95% credible interval for the averaged trajectory is shown by the grey shaded region. Population-level parameters for this model are shown in Supplementary Table 3.

Within this regression framework, we incorporated information on age, disease severity and sex (Methods) to see if the goodness of fit could be improved. We added fixed effects for these three variables, both separately and in combination. We separately examined age as a continuous variable or a dichotomous variable (stipulating whether patients were 60 years of age or over). The goodness of fit to the data did not vary appreciably: here we report results for age as a dichotomous variable. As we had relatively few (29) samples from patients with severe disease in the subset of the data considered here, we pooled patients with moderate disease and severe disease together. The inclusion of the fixed effects for age, sex, or severity did not improve the model fit (Table 2). In other words, the wide variation between individuals observed in the viral load dynamics could not be explained by the inclusion of these variables (Figure 3).

• Download figure
• Open in new tab

Figure 3: Inclusion of fixed effects in the regression models.

All regression models included study- and patient-specific random effects for the peak and slope (i.e. rate of decline) of the viral load. We then added fixed effects, both separately and in combination, to see if the model fit could be improved (Table 2). These fixed effects were: age, sex, and severity of disease. Here we show results for the three models containing one fixed effect (left: severity; middle: age; right: sex). The inclusion of severity, age, or sex did not improve the goodness of fit (Table 2). In the modelling here, age was included as a binary variable (under or over 60 years of age).

We used simulation-based methods (Methods) to estimate the power of our analysis to detect different effect sizes for severity, age and sex on peak viral load (Supplementary Figure 4). For severity and sex, our analysis has 80% power to detect a difference of around 1.1 in the log10 viral load (about a 12-fold difference in viral load), whereas power was lower for detecting the effect of age. Importantly the detectable differences are considerably smaller than the inter-individual variation in viral load at any given time point (Supplementary Table 3), indicating that these are not the major determinants of viral load.

Fitting a mechanistic model to the viral load data

In addition to modelling viral load decline using regression models, we also developed a mechanistic model which we fitted to the dataset. We elected to keep the model relatively simple due to lack of identifiability between more complex model structures. We represented the multi-faceted immune response to the infection via two components. First, the exponential growth of viral load in the initial phase of infection is brought under control by an early immune response. Subsequent to this, the infection is gradually cleared by a late immune response. The early immune response is stimulated by a high load of infected cells and starts to block viral replication and the invasion of susceptible cells. The late response requires a maturation phase before it becomes effective and is therefore, more representative of the adaptive immune response. However, we do not attempt to fully distinguish between innate and adaptive responses in this model, as the interplay between them is complex.

To guide the model fitting (see Methods), we made the assumption that both the peak viral load and the activation of the early immune response should roughly be concurrent with the onset of symptoms. The vast majority of subjects in the dataset were only under observation after the onset of symptoms, which means we are unable to infer the rate of exponential growth during the initial phase of the infection. For each subject, we fitted two parameters in the mechanistic model, holding all other parameters fixed (see Methods). One of these free parameters governs the density of infected cells required to activate the early response, while the other determines the rate at which the late response clears infected cells and, therefore, the rate at which the viral load declines. As we did for the regression modelling, we used a nested random effect structure, with study- and patient-specific random offsets for both of these parameters (Methods). In addition, we incorporated data points that fell below the limit of detection in each study by accounting for censoring in the likelihood (Methods).

This approach allowed us to characterise the average time-course of an infection at the population-level i.e. after removing study- and patient-level offsets (Figure 4, Supplementary Table 4), as well as showing the dynamics for the study specific fits to the data. There was no significant relationship between the subject-specific random effects for either of the model parameters and maximum disease severity (Supplementary Figure 5), consistent with the outcome of the preceding regression modelling (Table 2). Similarly, there was no significant relationship between the subject-specific random effects and subject age or sex (not shown). Supplementary Figure 6 shows the 95% credible intervals for the study-specific random effects, showing the between-study variation in both the peak viral load (panel a) and the rate at which the infection is cleared in the URT (panel b). Supplementary Figure 7 shows modelled viral loads for each of the 155 subjects used to fit the mechanistic model.

• Download figure
• Open in new tab

Figure 4: Averaged trajectories obtained from the fitted mechanistic model.

In this plot we show the average trajectory predicted for each study (coloured lines), generated using the median value used for the initial viral load at t=0 in each case. We also show the average trajectory across all the studies, indicated by the black line. The dark grey shaded area indicates the 95% credible interval for this average trajectory. The light grey area accounts for the variation observed around the average trajectory (generated using samples from s, as defined in in Equation 8, and calculating the 95% prediction interval for the population-level dynamics). The fit to data from Study 3 is not shown, as this study only contained one patient, which means one cannot distinguish between study- and patient-specific random effects. The opaque black circles are the data points from the 16 studies used to fit the model. For illustrative purposes, viral samples that were negative for virus are set to 1 viral copy per ml (i.e. 0 on the log-scale). The results from the mechanistic model presented here were obtained using 1500 samples from the posterior distribution, with the median trajectories plotted.

Relating viral load dynamics to SARS-CoV-2-specific immune responses

To demonstrate the potential to link viral load dynamics to SARS-CoV-2-specific adaptive immune responses in individual subjects, we made use of detailed data on neutralising antibody and T-cell responses from the study by Tan et al. (7). We compared this immunological data with the patient-specific late immune response parameter fitted in our mechanistic model in these 12 patients, to attempt to gain insight into within-host mechanisms that might control viral load (Figure 5). Specifically, we assessed the total interferon-γ (IFN-γ) T-cell response to SARS-CoV-2 peptides measured in an ELISPOT assay, and results from a surrogate virus antibody neutralization assay. Simple logistic curves were fitted to the data from each subject’s measured immune responses, and the area under each curve (AUC) was used as a measure of the magnitude of each subject’s immune responses (Figure 5b and 5c). Overall, no significant correlation was observed between the modelled late immune response and the total interferon-γ (IFN-γ) T-cell response (r=0.39, p value = 0.208, see Figure 5d) or the neutralizing antibody response (r=0.07, p value = 0.831, see Figure 5e). However, we noted two subjects (subjects 1 and 12) displayed quite distinct immune responses to the others, characterized by absent neutralizing antibody, and initially high but subsequently declining total interferon-γ (IFN-γ) T-cell responses, possibly indicating a qualitatively different response to SARS-CoV-2. When these two subjects were removed from the analysis, a much stronger correlation was observed between the modelled late immune response and surrogate virus antibody neutralization (r=0.79, p value = 0.006).

• Download figure
• Open in new tab

Figure 5: Paired viral load and immune response dynamics.

These panels show data from 12 patients, reported by Tan et al. (7). Panel (a): viral load measurements (points) and modelled viral load trajectories from the mechanistic model (black lines show the posterior means, shaded areas are the 95% credible intervals). The coloured symbols indicate the severity score recorded for each subject (on the WHO severity scale). Panel (b): measured total T cell response (blue symbols), rescaled by the largest observed measurement. A logistic curve was fitted through the points for each patient, to facilitate the area under the curve (AUC, blue shaded area) calculation. To calculate both AUCs (antibody & T cell) we used only the first 15 days after symptom onset, as this was the time period used to fit our models. Panel (c): neutralising antibody response (purple dots). A logistic curve was fitted through the points for each patient, to facilitate the area under the curve (AUC, purple shaded area) calculation. Panel (d): relationship between the calculated AUCs of the T cell responses and modelled patient-specific immune responses (p value = 0.208). Panel (e): relationship between the calculated AUCs of the antibody responses and modelled patient-specific immune responses. Two patients failed to mount an antibody response which neutralised virus. The correlation between the patient-specific response and the AUC is much stronger when these patients are not included (p value = 0.006, compared to p value = 0.831 when all 12 subjects are considered). We note that for subjects 4 and 8 (open circles in panels d and e), fewer than 3 viral load measurements were available, meaning their fitted parameters may shrink to the study mean.

Data

To collect data on viral load dynamics, we searched PubMed and medRxiv for studies that recorded longitudinal viral load data from individuals with symptomatic COVID-19 infections. Searches were carried out between May 20^th 2020 and February 11^th 2021. In particular, we searched for studies which reported the timing of symptom onset for each patient, which we used to temporally align samples from different patients. We identified 53 studies, data from 5 of which could be extracted from the publication or preprint. We contacted the authors of 48 studies by email to request access to patient-level data, and we received data from 14 of these studies. 2 of the 19 studies were dropped from this analysis, as they contained little or no longitudinal data for URT samples. Studies for which data was obtained are summarised in Table 1: all remaining studies identified by our literature search are summarised in Supplementary Table 1. We extracted as much information as possible on the severity of illness experienced by the patients and extracted demographic information on the patients where available. Two authors independently studied the severity information and matched the descriptions to the WHO clinical progression scale (34). Some articles contained detailed information on the course of disease for each patient, but this was not available in all the studies. We considered scores of 1-3 to represent mild disease, 4-6 to represent moderate disease and scores above 6 to represent severe disease. This is slightly different to the patient state descriptors of the WHO scale, because the majority of studies did not provide sufficient detail about the method of oxygen delivery to allow us to distinguish between scores of 5 and 6. In some studies, the severity of symptoms was recorded at multiple timepoints, over the course of the infection. Here, we use the term ‘severity’ to describe the maximum severity of disease experienced by each patient.

Although some studies with a long follow-up period demonstrate that some patients can remain PCR-positive for virus for well over a month (39,40), several studies (16,17) have demonstrated that very few people are infectious, based on URT swabs, beyond 10 or 14 days after symptom onset. Data collected after this period is likely to reflect RNA debris that remains in the URT. Therefore, our models of viral load are fitted to data points within the first 15 days of symptoms. This means that we have fewer data points to fit to (Table 1), but we believe that the meaningful dynamics, in terms of how the host controls the infection, are found in this time period.

Viral load Quantification

In all studies real-time polymerase chain reaction (PCR) assays were used to quantify URT virus either as Ct values, or as viral copies /mL (V). As discussed in Ref. (41), calculation of viral copies per ml from Ct values requires the use of a ‘standard curve’, which is calibrated to the experimental set-up in a particular laboratory using reference samples. In general, these curves have the form:

Here a and b are positive numbers that fully specify the viral load for a given Ct value. This relationship indicates that there is a linear relationship between log-transformed viral loads and the raw Ct values, with higher Ct values representing lower viral loads. The units of V vary between studies (e.g. viral copies per ml, per swab, or per 1000 cells), all studies collected here that have quantified viral load have used viral copies per unit of volume. We collected as many different standard curves as possible, from the studies included in this analysis (42,43) and from other papers in the COVID-19 literature (33,44–46), to understand the variation observed. From these, we determined an ‘averaged’ standard curve (Supplementary Figure 1), using the mean observed values for parameters a and b, which we used to estimate the viral loads for studies for which only Ct values were available. This enabled us to pool data from all 17 studies.

Models

We sought to explain variation in viral load (either its peak value or its rate of decline over time) among patients, due to (e.g.) age, sex, and severity of disease. We did this using two types of models, a linear mixed effects regression approach and a mechanistic model, which took the form of a system of first-order differential equations. Both models were fitted to viral load data from within 15 days of symptom onset, using only subjects with least three samples taken during this time period (models fitted to 870 samples from 155 patients from 16 studies). Therefore, as indicated in Table 1, no data from study 13 was used to fit the models.

Linear regression models

Bayesian regression models were fitted using RStan (47), with some of the analysis carried out using the rethinking package (48). Linear models were fitted to log-transformed viral loads, with random effects for each patient and study applied to the parameters determining the peak viral load, which we assumed coincides with the onset of symptoms, and rate of its decline over time. Samples for which no virus was detected were treated as being below limit of detection (LOD), rather than truly virus negative. We show the general form of the regression model here, where L is the likelihood of each data point, to illustrate the random-effect structure and the censoring of the data:

This model is relatively simple: the log-transformed viral loads are captured by a linear model, with the intercept describing the viral load at the time of symptom onset, and the slope capturing the rate at which the viral load subsequently declines. The log-transformed viral load for a given patient is normally distributed around a modelled trajectory, which is described by μ, with a standard deviation given by σ. For data points where no virus was detected, the likelihood takes a different form. We calculate the likelihood as the probability of the viral load being below the LOD, writing N_CDF as the cumulative distribution function of the normal distribution. We allow for the fact that the LOD is study-specific (indicated in Figure 2). All remaining terms in Equation 2 define the prior distributions for the parameters. Both the slope and the intercept have a parameter that captures the average behaviour across all patients in all studies (a₀ and b₀ respectively). Random effects capture patient- and study-specific deviations from the average behaviour. Both the patient- and study-specific random effects are normally distributed with zero mean: here we write N₂ to indicate the bivariate normal distribution. We allow correlation between the random effects on the intercept and slope at each level, as indicated by the presence of parameters ρ_study and ρ_patient. For these parameters, we use the LKJ distribution (49) as a prior.

Within this framework, fixed effects could then be added to the model, to assess whether (e.g.) severity of disease, or the sex of patient affected the typical viral load trajectory observed. To do this, we simply add terms to μ in Equation 1. Here we illustrate this for sex, making female the reference category:

Here variable ‘Male’ is equal to 1 if the patient is male (0 otherwise), and we have also added zero-mean priors for the new parameters. As mentioned above, severity was expressed on the WHO scale, which takes values from 1-10. As 1 is asymptomatic, severity in our dataset is limited from 2 (mild symptoms, not hospitalised) to 10 (dead). As the dataset contains relatively few samples from patients with severe disease (score 7 to 10 on the WHO scale), we chose to make two groups from these categories: Mild (2 or 3) and Moderate or Severe (4-10). When severity is included in the model, μ becomes:

Here variable ‘ModSev’ is equal to 1 for patients with moderate or severe disease (0 otherwise). In most studies, patient age is given in years, as an integer value. Some studies expressed age of patients in decades e.g. 10-20 years of age. For these studies, we used the midpoint of the age range given for each patient. We grouped age into two groups: _ <60 years of age and ≥60 years of age. In the regression modelling, we added a fixed effect for age, making the youngest age group the default. Here, μ takes the form:

Here variable ‘Age’ is equal to 1 for patients aged 60 or over (equal to 0 otherwise). We also looked at including subject age as a continuous variable, but the goodness of fit did not change appreciably. We added the 3 terms (age, severity, sex) to the regression models both separately and in combination. We assessed the goodness of fit using the Watanabe Akaike Information Criterion (WAIC) (49), with the best-fitting model having the lowest WAIC value. As the WAIC for each candidate model is estimated from a finite sample, its standard error was calculated using the rethinking package (48) to appreciate the uncertainty in its value. These standard errors are useful when appraising the differences in WAIC between candidate models. The relative goodness of fit of a given model can also be appraised by calculating its Akaike weight among the set of all considered models. This can be interpreted as the probability that this model, out of the set of models considered, would provide the best fit to new (i.e. out of sample) data (49).

We use simulation-based estimation of statistical power to assess our capacity to detect a difference in viral load dynamics due to one of the three factors (sex, age, severity of disease) assessed here. We generated synthetic datasets of the same size as the one considered here, with study- and patient-specific random effects of the same magnitude (Supplementary Table 3). For simplicity, we generated datasets with an equal number of samples per subject, and the sampling times were randomly generated from a uniform distribution. When generating the data, we assumed that the peak viral load was influenced by one of the three aforementioned factors. We then ran the regression analysis, to see if the modelled effect could be detected. To reduce computation time, we here used frequentist regression via the lme4 package in R (50), with a p value < 0.05 for the included fixed effect indicating a significant finding. Generating 1000 synthetic datasets, the statistical power can be estimated as the percentage of datasets for which the regression analysis found a significant effect. Supplementary Figure 4 shows how the statistical power varies with the magnitude of the modelled effect, which is here expressed as a fold difference in the peak viral load. The results of the power analysis suggest that we were not well powered to measure relatively small differences in peak viral load (under 10-fold difference). The power to measure sex-specific differences in peak viral load was slightly higher than the power to measure severity- or age-specific differences, as we have a better balance between samples from male and female samples than we do between samples in under 60s and over 60s, or samples from subjects from mild disease versus samples from subjects with moderate or severe disease.

Mechanistic model

We also developed a mechanistic model to describe the viral dynamics. We started modelling the infection 5 days before the onset of symptoms. At this time, an initial dose of virus V₀ is introduced. These virions start invading susceptible cells at rate β. Infected cells then produce more virus at rate p, which can then invade subsequent susceptible cells. Meanwhile, free virions are cleared at rate γ. In this way, the viral population grows exponentially. In the model, this growth is brought under control by an early phase of the immune response, and then the infection is gradually cleared by a later phase of the immune response. We do not seek to specify the immunological mechanisms contributing to these phases of the immune response. In each infection, the exponential growth slows as the population of infected cells approaches a certain value (I_max). This reflects a combination of two within-host mechanisms: the depletion of susceptible cells in the URT and the effect of an early phase immune response which is triggered at a high level of infection. Since these mechanisms may be linked, i.e. the immune response may modify the susceptibility of target cells, we do not attempt to distinguish between them, or to explicitly model the population of susceptible cells.

After the exponential growth is brought under control, the infection is cleared by the late immune response. In this model, this is triggered by a certain density of infected cells but requires a maturation stage before it becomes effective at clearing infected cells. We write the system of equations as:

The late immune response, represented here by variables (A₁, A₂, A₃), is stimulated by the presence of infected cells. A Hill function is used to make the adaptive response dimensionless, and to rescale its value (daily input into compartment A₁ scales between 0 and 1). Parameter I₅₀ (here fixed to a value of 1000 cells) determines the magnitude of the density of infected cells required to stimulate this response. It requires time to mature (governed by rate parameter k_c, which we fix at 0.33day^-1), meaning that only stage A₃ is able to clear infected cells. Once mature, the late response does not wane in this model, as we are only interested in its ability to clear an infection, not how long it persists after the infection has been resolved. In this way the late response recapitulates features of adaptive cell-mediated and humoral immunity, without needing to specify their relative contributions. Parameter k_a represents the maximum rate at which the late response can clear the infection, and A₅₀ governs the magnitude of the adaptive response required effectively clear infected cells. We note that the magnitude of the late immune response has been set to be of order 1 for convenience, and this determines the magnitude of A₅₀.

For the vast majority of patients for whom we have data, only samples taken after the onset of symptoms are recorded. This means that we are unable to estimate the duration of the incubation period, the initial dose of virus that causes the infection, or the rate at which the virus reproduces before being acted upon by the immune response. We elect to model a five-day incubation period and fix the rate of exponential growth to be the same for all patients, setting β = 0.8 ml day^-1virion^-1, γ = 13day^-1, p = 80day^-1, with these values informed by parameter values chosen in published models (25,26,28). This is a simplification: it is unclear whether, in reality, the time from infection to symptom onset varies by age, sex or disease severity. For each patient the infection starts at Day -5 (in the dataset, Day 0 is the day of symptom onset, not the day of infection onset) with initial conditions given by (I = 0, V = V₀, A₁ = 0, A₂ = 0, A₃ = 0) i.e. the infection starts with an initial viral dose, no infected cells, and no prior exposure to the virus (late immune response completely inactive).

The initial viral dose, , for each patient is chosen with the idea in mind that the viral load should peak at or close to the time of symptom onset. Ideally, , should be fitted simultaneously with the free parameters for each patient, but this has proved challenging to date, especially given that peak viral load should coincide approximately with symptom onset but may also not be observed. Briefly, , is chosen by considering the linear subsystem of equations that govern the growth of the virus prior to the activation of the immune response to the infection. As this system of equations is linear, it can be solved exactly (we use a matrix exponential). We require that be chosen so that, just prior to the onset of symptoms, an uncontrolled infection would reach the maximum viral load observed for patient i in the dataset. This ensures that the peak viral load modelled for each patient is controlled by the early-stage immune response, as it is too soon for the adaptive response to be effective at clearing infected cells.

When fitting to the data, we allowed the values of k_a and I_max to vary and we fixed all the other parameters. As we did for the regression modelling, we included random effects in the model, to account for different behaviour in different patients and across different studies. In the case where we have P patients taken from S studies, we write:

This means that and represents the population-level average value of each parameter. The system of equations was solved numerically in R, using the dopri function from the dde package (51). We write the modelled viral load trajectory at time t for patient i in study j as V_ij(t). Samples from the posterior distribution of the fitted parameter were obtained using Markov Chain Monte Carlo methods, via the R package drjacoby (52). As with the regression modelling, the form of the likelihood for each data point depends on whether the sample is above the study-specific LOD or is recorded as being virus-negative. We write D_ijk to indicate the k^th viral load sample, taken t_k after the onset of symptoms, from patient i in study j. The likelihood for each data point has the form:

The global likelihood was then obtained by multiplying together the likelihoods for each data point. We provided prior distributions for the parameters to be fitted. These were: