Densely sampled viral trajectories suggest longer duration of acute infection with B.1.1.7 variant relative to non-B.1.1.7 SARS-CoV-2

Ethics

Residual de-identified viral transport media from anterior nares and oropharyngeal swabs collected from players, staff, vendors, and associated household members from a professional sports league were obtained from BioReference Laboratories. In accordance with the guidelines of the Yale Human Investigations Committee, this work with de-identified samples was approved for research not involving human subjects by the Yale Internal Review Board (HIC protocol # 2000028599). This project was designated exempt by the Harvard IRB (IRB20-1407).

Study population

The data reported here represent a convenience sample including team staff, players, arena staff, and other vendors (e.g., transportation, facilities maintenance, and food preparation) affiliated with a professional sports league. Clinical samples were obtained by combined swabs of the anterior nares and oropharynx administered by a trained provider. Viral concentration was measured using the cycle threshold (Ct) according to the Roche cobas target 1 assay. For an initial pool of 298 participants who first tested positive for SARS-CoV-2 infection during the study period (between November 28^th, 2020 and January 20^th, 2021), a diagnosis of “novel” or “persistent” infection was recorded. “Novel” denoted a likely new infection while “persistent” indicated the presence of virus in a clinically recovered individual. A total of 65 individuals (90% male) had novel infections that met our inclusion criteria: at least five positive PCR tests (Ct < 40) and at least one test with Ct < 35. Seven of these individuals were infected with the B.1.1.7 variant as confirmed by genomic sequencing.

Genome sequencing and lineage assignments

RNA was extracted from remnant nasopharyngeal diagnostic specimens and used as input for SARS-CoV-2 genomic sequencing as previously described.⁶ Samples were sequenced on the Oxford Nanopore MinION. Consensus sequences were generated using the ARTIC Network analysis pipeline⁷ and samples with >80% genome coverage were included in analysis. Individual SARS-CoV-2 genomes were assigned to PANGO lineages using Pangolin v.2.1.8.⁸ All viral genomes assigned to the B.1.1.7 lineage were manually examined for representative mutations.⁹

Converting Ct values to viral genome equivalents

To convert Ct values to viral genome equivalents, we first converted the Roche cobas target 1 Ct values to equivalent Ct values on a multiplexed version of the RT-qPCR assay from the US Centers for Disease Control and Prevention.¹⁰ We did this following our previously described methods.³ Briefly, we adjusted the Ct values using the best-fit linear regression between previously collected Roche cobas target 1 Ct values and CDC multiplex Ct values using the following regression equation:

Here, y_i denotes the i^th Ct value from the CDC multiplex assay, x_i denotes the i^th Ct value from the Roche cobas target 1 test, and ε_i is an error term with mean 0 and constant variance across all samples. The coefficient values are β₀ = –6.25 and β₁ = 1.34.

Ct values were fitted to a standard curve in order to convert Ct value data to RNA copies. Synthetic T7 RNA transcripts corresponding to a 1,363 b.p. segment of the SARS-CoV-2 nucleocapsid gene were serially diluted from 10⁶-10⁰ RNA copies/μl in duplicate to generate a standard curve¹¹ (Supplementary Table 1). The average Ct value for each dilution was used to calculate the slope (−3.60971) and intercept (40.93733) of the linear regression of Ct on log-10 transformed standard RNA concentration, and Ct values from subsequent RT-qPCR runs were converted to RNA copies using the following equation:

Here, [RNA] represents the RNA copies /ml. The log₁₀(250) term accounts for the extraction (300 μl) and elution (75 μl) volumes associated with processing the clinical samples as well as the 1,000 μl/ml unit conversion.

Model fitting

For the statistical analysis, we removed any sequences of 3 or more consecutive negative tests to avoid overfitting to these trivial values. Following our previously described methods,³ we assumed that the viral concentration trajectories consisted of a proliferation phase, with exponential growth in viral RNA concentration, followed by a clearance phase characterized by exponential decay in viral RNA concentration.¹² Since Ct values are roughly proportional to the negative logarithm of viral concentration¹³, this corresponds to a linear decrease in Ct followed by a linear increase. We therefore constructed a piecewise-linear regression model to estimate the peak Ct value, the time from infection onset to peak (i.e. the duration of the proliferation stage), and the time from peak to infection resolution (i.e. the duration of the clearance stage). The trajectory may be represented by the equation

Here, E[Ct(t)] represents the expected value of the Ct at time t, “l.o.d” represents the RT-qPCR limit of detection, δ is the absolute difference in Ct between the limit of detection and the peak (lowest) Ct, and t_o, t_p, and t_r are the onset, peak, and recovery times, respectively.

Before fitting, we re-parametrized the model using the following definitions:

• ΔCt(t) = l.o.d. – Ct(t) is the difference between the limit of detection and the observed Ct value at time t.

• ω_p = t_p - t_o is the duration of the proliferation stage.

• ω_c = t_r - t_p is the duration of the clearance stage.

We constrained 0.25 ≤ ω_p ≤ 14 days and 2 ≤ ω_p ≤ 30 days to prevent inferring unrealistically small or large values for these parameters for trajectories that were missing data prior to the peak and after the peak, respectively. We also constrained 0 ≤ δ ≤ 40 as Ct values can only take values between 0 and the limit of detection (40).

We next assumed that the observed ΔCt(t) could be described the following mixture model: where E[ΔCt(t)] = l.o.d. - E[Ct(t)] and λ is the sensitivity of the q-PCR test, which we fixed at 0.99. The bracket term on the right-hand side of the equation denotes that the distribution was truncated to ensure Ct values between 0 and the limit of detection. This model captures the scenario where most observed Ct values are normally distributed around the expected trajectory with standard deviation σ(t), yet there is a small (1%) probability of an exponentially distributed false negative near the limit of detection. The log(10) rate of the exponential distribution was chosen so that 90% of the mass of the distribution sat below 1 Ct unit and 99% of the distribution sat below 2 Ct units, ensuring that the distribution captures values distributed at or near the limit of detection. We did not estimate values for λ or the exponential rate because they were not of interest in this study; we simply needed to include them to account for some small probability mass that persisted near the limit of detection to allow for the possibility of false negatives.

We used a hierarchical structure to describe the distributions of ω_p, ω_r, and δ for each individual based on their respective population means μ_ωp, μ_ωr, and μ_δ and population standard deviations σ_ωp, σ_ωr, and σ_δ such that

We inferred separate population means (μ_•) for B.1.1.7- and non-B.1.1.7-infected individuals. We used a Hamiltonian Monte Carlo fitting procedure implemented in Stan (version 2.24)¹⁴ and R (version 3.6.2)¹⁵ to estimate the individual-level parameters ω_p, ω_r, δ, and t_p as well as the population-level parameters σ*, μ_ωp, μ_ωr, μ_δ, σ_ωp, σ_ωr, and σ_δ. We used the following priors:

Hyperparameters:

Individual-level parameters:

The values in square brackets denote truncation bounds for the distributions. We chose a vague half-Cauchy prior with scale 5 for the observation variance σ*. The priors for the population mean values (μ_•) are normally distributed priors spanning the range of allowable values for that parameter; this prior is vague but expresses a mild preference for values near the center of the allowable range. The priors for the population standard deviations (σ_•) are half Cauchy-distributed with scale chosen so that 90% of the distribution sits below the maximum value for that parameter; this prior is vague but expresses a mild preference for standard deviations close to 0.

We ran four MCMC chains for 1,000 iterations each with a target average proposal acceptance probability of 0.8. The first half of each chain was discarded as the warm-up. The Gelman R-hat statistic was less than 1.1 for all parameters. This indicates good overall mixing of the chains. There were no divergent iterations, indicating good exploration of the parameter space. The posterior distributions for μ_δ, μ_ωp, and μ_ωr, were estimated separately for individuals infected with B.1.1.7 and non-B.1.1.7. These are depicted in Figure 1 (main text). Draws from the individual posterior viral trajectory distributions are depicted in Supplementary Figure 1. The mean posterior viral trajectories for each individual are depicted in Supplementary Figure 2.

Checking for influential outliers

To examine whether the posterior distributions for the B.1.1.7-infected individuals reflected the influence of a single outlier, we re-fit the model seven times, omitting one of the B.1.1.7 trajectories each time. The inferred parameter values were fairly consistent, though omitting either of two of the B.1.1.7 cases (cases 5 and 6 in Supplementary Table 2). yields an infection duration with a 90% credible interval that overlaps with that of the non-B.1.1.7 90% credible interval for infection duration.

• Download figure
• Open in new tab

Supplementary Figure 1. Ct values for 65 individuals with estimated viral trajectories.

Each pane depicts the recorded Ct values (points) and derived log-10 genome equivalents per ml (log(ge/ml)) for a single person during the study period. Points along the horizontal axis represent negative tests. Time is indexed in days since the minimum recorded Ct value (maximum viral concentration). Individuals with confirmed B.1.1.7 infections are depicted in red. Non-B.1.1.7 infections are depicted in blue. Lines depict 100 draws from the posterior distribution for each person’s viral trajectory.

• Download figure
• Open in new tab

Supplementary Figure 2. Mean posterior viral trajectories for each person in the study.

Lines depict the posterior mean viral trajectory specified by the posterior mean proliferation time, mean clearance time, and mean peak Ct. Trajectories are aligned temporally to have the same peak time. B.1.1.7 trajectories are depicted in red, non-B.1.1.7 in blue.