Combining culture probability data
To estimate the association between viral load and culture probability, we used data previously described by Wölfel (
19) and Perera (
20). Four other datasets could not be included because Ct values were not converted to viral loads (
35,
46,
61,
62). The data from the study by van Kampen
et al. (
63) were not included because they differed (by viral load of ~1.0) from the data used for the current analysis (
97); this is likely due to a combination of factors including many patients who were in critical or immunocompromised condition, a high proportion of samples obtained from the lower respiratory tract (including late in the infectious course), and likely differences in cell culture trials. It is unsurprising that these data result in a shifted viral load/culture probability curve, and we excluded them because our focus was largely on first positive RT-PCR results from the upper respiratory tract, including from many subjects who were PAMS. [See (
97) for a figure comparing the plot of the van Kampen dataset to the two we used.] To calculate the expected culture probability, by age (as in
Fig. 2D) or by day from peak viral load (as in
Fig. 4C), we combined the estimated viral loads (
Figs. 2A and
4B) with the results of the regression of culture probability shown in
Fig. 2C. We used posterior predictions from the age regression model, which reflect the variation of viral load within age groups, to estimate culture probabilities by age. For instance, to obtain the culture probability for a specific age and group, we look up the estimated (expected) viral load for that group, add an error term according to the estimated error variance, and, using the association shown in
Fig. 2C, determine the expected culture probability. We used expected time courses (i.e., the model’s best guess for a time course) to estimate culture probability time courses.
Estimating viral load time course
Each RT-PCR test in our dataset has a date, but no information regarding the suspected date of subject infection or onset of symptoms (if any). Although determining the day of peak viral load for a single person based on a series of dated RT-PCR results would not in general be feasible because of individual variation, data from a large enough set of people would enable the inference of a clear and consistent model of viral load change over time with very few assumptions.
We included a single leading and/or trailing negative RT-PCR result, if dated within 7 days of the closest positive RT-PCR. To produce a model of typical viral load decline on a reasonable single-infection time scale, we excluded subjects whose full time series contains positive RT-PCRs spread over a period exceeding 30 days. Such time series may be attributable to contamination, to later swabbing that picks up residual RNA fragments in tonsillar tissue (
66), or to re-infection (
67–
69), or they may represent atypical infection courses (such as in immunocompromised or severely ill elderly patients) (
70). We excluded data from subjects with an infection delimited by both an initial and a trailing negative test when there was only a single positive RT-PCR result between them.
We estimated the slopes for a model of linear increase and then decline of log10(viral load). To compensate for the absence of information regarding time of infection, we also estimated the number of days from infection to the first positive test for each participant, so as to position the observed time series relative to the day of peak viral load. The analysis was implemented in two ways. Initially, simulated annealing was used to find an optimized fit of the parameters, minimizing a least-squares error function. Second, a Bayesian hierarchical model estimated subject-specific time courses, imputed the viral load assigned to each initial or trailing negative test, and captured effects of age, gender, clinical status, and RT-PCR system with model parameters. We tested both methods on data subsets ranging from subjects with at least three to at least nine RT-PCR results. The two methods produced results that were in generally good agreement (table S5). The finer-grained Bayesian approach appears more sensitive than the simulated annealing; its results, for subjects with at least three RT-PCR results, are those described in the main text.
Simulated annealing approach: A simulated annealing optimization algorithm (
71) was used to adjust the time series for each subject slightly earlier or later in time, by amounts drawn from a normal distribution with mean 0.0 and standard deviation 0.1 days. The error function was the sum of squares of distances of each viral load from a viral load decline line whose slope was also adjusted as part of the annealing process. In the error calculation, negative test results were assigned a viral load of 2.0, in accordance with our SARS-CoV-2 assay limit of detection and sample dilution (
19). The initial slope of the decline line was set to –2.0 and was varied using
N(0, 0.01). A second, optional, increase line initialized with a slope of 2.0, adjusted using an
N(0, 0.01) random variable, was included in the error computation if the day of a RT-PCR test was moved earlier than day zero (the modeled day of peak viral load). The height of the intercept (i.e., the estimated peak viral load) between the increase line (if any) and the decline line was also allowed to vary randomly [starting value 10.0, varied using
N(0, 0.1)]. The full time series for each subject was initialized with the first positive result positioned at day 2 +
N(0.0, 0.5) after peak viral load. The random-move step of the simulated annealing modified either of the two slopes or the intercept, each with probability 0.01, otherwise (with probability 0.97) one subject’s time series was randomly chosen to be adjusted earlier or later in time. After the simulated annealing stage, each time series was adjusted to an improved fit (when possible) based on the optimized increase and decline lines. Linear regression lines were then fitted through the results occurring before and after the peak viral load (
x = 0) and compared to the lines with slopes optimized by the simulated annealing alone. This final step helped to fine-tune the simulated annealing, in particular sometimes placing a time series much earlier or much later in time after it had stochastically moved initially in a direction that later (when the increase and decline line slopes had converged) proved to be suboptimal. The slopes of the lines fitted via linear regression after this final step were in all cases very similar (generally ±0.1) to those produced by the initial simulated annealing step. The final adjustments can be regarded as a last step in the optimization, using a steepest-descent movement operator instead of an uninformed random one. A representative optimization run for subjects with at least three RT-PCR results is shown in fig. S12.
Bayesian approach: The Bayesian analysis of viral load time course implements the same basic model, and additionally estimates associations of model parameters with covariates age, gender, B.1.1.7 status, and clinical status, estimates subject-level parameters (slope of log
10 viral load increase, peak viral load, slope of log
10 viral load decrease) as random effects, and accounts for effects of PCR system and test center types with random effects. To estimate the number of days from infection to the first test (henceforth “shift”), we constrained the possible shift values from –10 to 20 days and used a uniform prior on the support. In contrast to the other subject-level parameters, we estimated subject-level shifts independently (i.e., without a hierarchical structure). Figure S7 shows the placement in time of individual viral loads after shifting for subjects with RT-PCR results from at least 3 days. Model parameters changed gradually when subsets of subjects with an increasing minimum number of RT-PCR results, from three to nine, were examined (fig. S11 and table S5). The viral load assigned to negative test results (which may include viral loads below the level of detection) is estimated with a uniform prior on the support from –Inf to 3 (see also the caption of fig. S7). Using prior predictive simulations, we specified (weakly) informative priors for this analysis. This analysis was implemented in Stan (
72), as described in (
97).
Checking convergence of the model parameters showed that although 99.3% of all parameters converged with an R-hat value below 1.1, some subject-level parameters of 118 subjects (among 4344 subjects with at least three RT-PCR results) showed R-hat values between 1.1 and 1.74. Inspection of these parameters showed that these convergence difficulties were due to observed time courses that could arguably be placed equally well at the beginning or a later stage of the infection. Figure S16 shows a set of 81 randomly selected posterior predictions, to give an impression of time-series placement; fig. S17 shows the 49 participants with the parameters with the highest R-hat values. Although the high R-hat values could be removed by using a mixture approach to model shift for these participants, in light of their low frequency we retained the simpler model to avoid additional complexity. Alternatively, constraining the shift parameter to negative numbers would also improve R-hat values for these subjects, at the cost of the additional assumption that infections are generally not detected weeks after infection.
Sensitivity analysis: In addition to examining the viral load time series of subjects with RT-PCR results on at least 3 days, we tested both approaches on data from subjects with results from a minimum of 4 to 9 days. Given the degree of temporal viral load variation seen in other studies (
18–
20,
35,
41,
46,
63,
73,
74) and in our own data, our expectation was that a relatively high minimum number of results might be required before reliable parameter estimates with small variance would be obtained, but this proved not to be the case. The simulated annealing approach was tested with a wide range of initial slopes and intercept heights as well as seven different methods for the initial placement of time series. In general, maximum viral load and decline slopes were robust to data subset and initial time-series location, although there was variation in the length of the time to peak viral load, depending on how early in time the time series were initially positioned, the initial slopes of the increase and decrease lines and height of the maximum viral load. This is as expected, as the settings of these parameters can be used to bias the probability that a time series is initially positioned early or late in time and how difficult it is for it to subsequently move to the other side of the peak viral load at day zero. Table S5 shows parameter values for both approaches on the various data subsets.
Onset of shedding: We define the onset of shedding as the time point at which the increasing viral load crosses zero of the log
10 y axis—that is, when just one viral particle was estimated to be present. Because the estimated time of infection depends on the estimated peak viral load and the slope with which viral load increases, the data should optimally include multiple pre-peak viral load test results for each individual. If, as in the current dataset, only a subset of subjects have test results from pre-peak viral load, a hierarchical modeling approach still allows calculating subject-level estimates. Intuitively, this approach uses data from all subjects to calculate an average slope parameter for increasing viral load. In addition, it models subject-level parameters as varying around the group-level parameter. To further refine the estimation of slope parameters, the model also uses the age (see fig. S10), gender, and clinical status as covariates. Because negative test results could be false negatives, viral loads for these tests are imputed (with an upper bound of 3). Subject-level peak viral load and declining slope are modeled with the same approach. More generally, using a hierarchical model and shrinkage priors for the effects of covariates results in more accurate predictions in terms of expected squared error (
75) compared to analyzing each subject in isolation, but the overall improvement introduces a slight bias toward the group mean, resulting in an underestimation of the true variability of subject-level parameters. This is especially the case if, as in the current dataset, subject-level data are sparse.
Onset of symptoms: The 317 onset-of-symptoms dates for hospitalized patients were collected as part of the Pa-COVID-19 study, a prospective observational cohort study at Charité–Universitätsmedizin Berlin (
76,
77), approved by the local ethics committee (EA2/066/20), conducted according to the Declaration of Helsinki and Good Clinical Practice principles (ICH 1996), and registered in the German and WHO international clinical trials registry (DRKS00021688).
RE: Estimating infectiousness throughout SARS-CoV-2 infection course
Dear Editor,
The modeling of infectiousness during the course of SARS-CoV-2 infection based on 4,344 cases presented by Jones et al. [1] is problematic. The conclusions of the authors - namely that "onset of infectious virus shedding, peak viral load and the rates of viral load increase and decline" can be inferred from their model and can be used to compare "between groups of subjects and virus strains"- are not well founded.
Scientific standards with respect to design, conduct and reporting of the study have been largely ignored, e.g. no clear hypothesis has been formulated and the study design is unclear. For example, "Study Composition" only provides scant information on the 25,381 test-positives during the study period February 24, 2020 to April 2, 2021, giving rise to 4,344 subjects with at least 3 RT-PCR tests used for complex statistical modelling. Test-positives were derived from three different sources: 1) 9,519 hospital patients who "tested positive in an in-patient context at any point in their infection", 2) 6,110 test-positives (out of a total of 6,159 of whom 49 were excluded because of later hospitalisation) detected in the 24 walk-in test centers in Berlin, and 3) 9,752 "other" subjects selected from unclear sources.
Stratification of first RT-PCR test results according to source of cases by age (table 1 and figure 1) shows that the majority of people in the younger age groups (0-24 years) were tested only once at "other" settings, with 5% of 0-4, 6% of 5-9 and 4% of 10-14 year-olds being tested more than 3 times. Young and middle-aged adults (age 25-34 and 35-44) were mostly derived from walk-in test centers. Test-positive adults from age group 45-64 originated from hospitals, test centers, and other sources in equal parts. Increasing adult age was related to being tested more frequently, ranging from 7% of cases age 25-34, 9% of cases age 35-44, 12% of cases age 45-54, 18% of cases age 55-64 to 32% (N=2,492/7,851) of cases age 65 and above being tested more than three times. Hence, more than half of the 4,344 cases used for modelling the course of infectiousness in the paper by Jones et al. were hospitalised individuals mostly between age 65 and 84. The remaining 1,852 cases are not further described with respect to timing and reason for repeated testing, for example because of suspicion of a new variant.
Details on RT-PCR testing, including gene targets, as well as viral culture settings are also lacking. For clarification, we refer to the INSTAND quality assurance assessment (Ringversuch) among participating laboratories in Germany.[2] A given log10 of 6 of the E-gene viral particle concentration and a similar log10 ORF gene viral particle concentration, corresponded to measured Ct-values of 25.5 and 25.4, respectively, by the Roche cobas system. Each log10 difference of 0.5 corresponded to a Ct value difference of 1.65. Applying these INSTAND figures to the results presented by the authors in table 1, the corresponding mean Ct values of hospitalised cases are 25 and above in all age groups, indicating low infectivity at first RT-PCR test result. The proportion of test positives in the hospitalised group is generally low (hospitalisation for other reasons?), ranging from 1% to below 3% across different age groups. In contrast, test-center cases from different age groups show a higher positivity rate of 3%-7% and much lower mean Ct-values ranging from below 25 to 22.5, indicating a higher risk of being infectious. For comparison, in our recently published analysis of 162,457 first RT-PCR test results from a defined source population measured by one single laboratory [3,4] using exclusively the Roche cobas system, Ct-values were highly variable and the detection of test-positives was dependent on age, calendar time, and the related national test strategy. Between calendar weeks 10-19 and 45-49, PCR-testing of mostly symptomatic persons yielded a higher positive rate of up to 8%, while primarily testing asymptomatic persons during the summer months, produced a low positive rate ranging from 0.4% to 1.8%. The ONS coronavirus infection community survey in the UK recently also reported a wide variation of Ct-values and a strong influence of calendar time on positivity rate at first RT-PCR test, ranging from 0.01% (May to September 14, 2020) to 2.2% (December 22, 2020 to January 4, 2021) and declining again to 0.3% in March 2021.
The model by Jones et al. suggests that low viral load at time of first measurement may be compatible with a pre-symptomatic state of super-spreaders with exponential growth of viral particles within the next hours up to 5 days. Besides the described methodological shortcomings implicating strong bias and uncontrolled confounding, the underlying assumption of a uniform pattern of infectiousness of SARS-CoV-2 is unfounded.
Since positive results of a first RT-PCR test are the basis for calculating "incidence" of SARS-CoV-2 infection, on which far-reaching non-pharmaceutical interventions such as isolation, quarantine, lock-down measures, and school closures are based and legally enforced, first RT-PCR positive test results are at the core of pandemic decision making and monitoring. Real-world evidence of SARS-CoV-2 infectiousness compared to modelling "infectiousness throughout SARS-CoV-2 infection course" becomes increasingly contradictory.
Andreas Stang, MD, MPH∗
Institute of Medical Informatics, Biometry and Epidemiology,
University Hospital Essen, Germany; School of Public Health,
Department of Epidemiology, Boston University, Boston, USA
Karl-Heinz Jöckel, PhD
Institute of Medical Informatics, Biometry and Epidemiology,
University Hospital Essen, Germany
Angela Spelsberg, MD, SM
Tumorzentrum Aachen e.V., Pauwelsstraße 30, 52074 Aachen,
Germany
Paul Cullen, MD, MSc
MVZ Labor Münster Hafenweg GmbH, Hafenweg 9-11; 48155,
Münster, Germany
Ulrich Keil, MD, PhD
Institute of Epidemiology and Social Medicine, University of Münster,
Albert Schweitzer Campus 1, 48149 Münster, Germany
∗Corresponding author:
Prof. Andreas Stang, MD, MPH, Institut für Medizinische Informatik, Biometrie und Epidemiologie, Universitätsklinikum Essen, Germany. Tel.: +49 201 723 77 201; fax: +49 201 723 77 333
E-mail addresses:
[email protected] (A. Stang) (ORCID: 0000-0001-6363-9061)
[email protected] (K.-H. Jöckel)
[email protected] (A. Spelsberg)
[email protected] (P. Cullen)
[email protected] (U. Keil)
[1] T. C. Jones et al., Science 10.1126/science.abi5273 (2021).
Paul Cullen has received speaker's fees from Roche Diagnostics. None of the other authors declares any conflict of interest.[2] INSTAND Quantitative Bezugsproben zur Verbesserung der Vergleichbarkeit und Bewertung von Laborergebnissen zum Virusgenom-Nachweis von SARS-CoV-2. Information zur Testdurchführung und Anwendung der quantitativen Bezugsproben incl. Anhang: Auswertung der Ergebnisse vom 1. Versand (03.11.2020) und 2. Versand (17.11.2020), 3. Versand vom 15.01.2021. INSTAND e.V. Berlin, 20210118g Begleitheft - quantitative Bezugsproben 1 und 2 - SARS-CoV-2.docx
[3] A. Stang et al. The performance of the SARS-CoV-2 RT-PCR test as a tool for detecting SARS-CoV-2 infection in the population. J Infection 2021, in press
[4] https://www.medrxiv.org/content/10.1101/2021.05.06.21256289v1
[5] S. Walker et al. Ct threshold values, a proxy for viral load in community SARS-CoV-2 cases, demonstrate wide variation across populations and over time. medRxiv preprint DOI: https://doi.org/10.1101/2020.10.25.20219048