Abstract
Extensive virological testing is central to SARS-CoV-2 containment, but many settings face severe limitations on testing. Group testing offers a way to increase throughput by testing pools of combined samples; however, most proposed designs have not yet addressed key concerns over sensitivity loss and implementation feasibility. Here, we combine a mathematical model of epidemic spread and empirically derived viral kinetics for SARS-CoV-2 infections to identify pooling designs that are robust to changes in prevalence, and to ratify losses in sensitivity against the time course of individual infections. Using this framework, we show that prevalence can be accurately estimated across four orders of magnitude using only a few dozen pooled tests without the need for individual identification. We then exhaustively evaluate the ability of different pooling designs to maximize the number of detected infections under various resource constraints, finding that simple pooling can identify up to 20 times as many positives compared to individual testing with a given budget. We illustrate how pooling affects sensitivity and overall detection capacity during an epidemic and on each day post infection, finding that sensitivity loss is mainly attributed to individuals sampled at the end of infection. Crucially, we confirm that our theoretical results can be accurately translated into practice using pooled human nasopharyngeal specimens. Our results show that accounting for variation in sampled viral loads provides a nuanced picture of how pooling affects sensitivity to detect epidemiologically relevant infections. Using simple, practical group testing designs can vastly increase surveillance capabilities in resource-limited settings.
Introduction
The ongoing pandemic of SARS-CoV-2, a novel coronavirus, has caused over 24 million reported cases of coronavirus disease 2019 (COVID-19) and 800,000 reported deaths between December 2019 and August 2020. (1) Although wide-spread virological testing is essential to inform disease status and where outbreak mitigation measures should be targeted or lifted, sufficient testing of populations with meaningful coverage has proven difficult. (2–7) Disruptions in the global supply chains for testing reagents and supplies, as well as on-the-ground limitations in testing throughput and financial support, restrict the usefulness of testing– both for identifying infected individuals and to measure community prevalence and epidemic trajectory. While these issues have been at the fore in even the highest-income countries, the situation is even more dire in low income regions of the world. Cost barriers alone mean it is often simply not practical to prioritize community testing in any useful way, with the limited testing that exists necessarily reserved for the healthcare setting. These limitations urge new, more efficient, approaches to testing to be developed and adopted both for individual diagnostics and to enable public health epidemic control and containment efforts.
Group or pooled testing offers a way to increase efficiency by combining samples into a group or pool and testing a small number of pools rather than all samples individually. (8–10) For classifying individual samples, including for diagnostic testing, the principle is simple: if a pool tests negative, then all of the constituent samples are assumed negative. If a pool tests positive, then the constituent samples are putatively positive and must be tested again individually (Fig. 1A). Further efficiency gains are possible through combinatorial pooling, where, instead of testing every sample in every positive pool, each sample can instead be represented across multiple pools and potential positives are identified based on the pattern of pooled results (Fig. 1B). (9,10)
Simple pooling designs can also be used to assess prevalence without individual specimen identification (Fig. 1C). It has already been shown that the frequency of positive pools can allow estimation of the overall prevalence. (11) Crucially however, we show here that prevalence estimates can be greatly honed by considering quantitative viral loads measured in each positive pool, rather than simply using binary results (positive / negative). In short, the viral (RNA) load measurement from a pool is proportional to the sum of the (diluted) viral loads from each positive sample in the pool. Thus, here we show how evaluating the viral loads greatly improves potential efficiency gains in prevalence estimates by providing crucial information on the estimated number of positive samples in the pool – when the expected distribution of viral loads across specimens is known, which is easily measured empirically in a given lab. (12,13) Although this approach requires more complex statistical methods, the efficiency gains for public health surveillance can be large, and simplifying templates can be produced to improve ease of use and access to these types of analyses. The outcome is a highly efficient method for estimating population prevalence and enabling robust public health surveillance where it was previously out of reach.
Whilst the literature on theoretically optimized pooling designs for COVID-19 testing has grown rapidly, formal incorporation of biological variation (i.e., viral loads) and incorporation of general position along the epidemic curve, has received little attention. (14–17) Test sensitivity for example is not a fixed value, but depends on viral load, which can vary by many orders of magnitudes across individuals and over the course of an infection. (18–20) This large variation within a single infection affects sensitivity to detect infections at different points in the disease course, which has implications for appropriate intervention and the interpretation of a viral load measurement from a sample pool.
Here, we comprehensively evaluate designs for pooled testing of SARS-CoV-2 whilst accounting for epidemic dynamics and variation in viral loads arising from viral kinetics and extraneous features such as sampling variation. We demonstrate efficient, logistically feasible pooling designs for individual identification (i.e., diagnostics) and prevalence estimation (i.e., population surveillance). To do this, we use realistic simulated viral load data at the individual level over time, representing the entire time course of an epidemic to generate synthetic data that reflects the true distribution of viral loads in the population at any given time of the epidemic. We then used these data to derive optimal pooling strategies for different use cases and resource constraints in-silico. Finally, we demonstrate the approach using discarded de-identified human nasopharyngeal swabs initially collected for diagnostic and surveillance purposes.
Results
Modelling a synthetic population to assess pooling designs
To identify optimal pooling strategies for distinct scenarios, we required realistic estimates of viral loads across epidemic trajectories. We developed a population-level mathematical model of SARS-CoV-2 transmission that incorporates empirically measured within-host virus kinetics, and used these simulations to generate population-level viral load distributions representing real data sampled from population surveillance, either using nasopharyngeal swab or sputum samples (Fig. 2). Full details are provided in Materials and Methods and Supplementary Material 1, sections 1-4. These simulations generated a synthetic, realistic epidemic with a peak daily per incidence of 19.5 per 1000 people, and peak daily prevalence of RNA positivity (viral load greater than 100 virus RNA copies per ml) of 265 per 1000 (Fig. 2E). We used these simulation data to evaluate optimal group testing strategies at different points along the epidemic curve for diagnostic as well as public health surveillance, where the true viral loads in the population is known fully.
Improved testing efficiency for estimating prevalence
We developed a statistical method to estimate prevalence of SARS-CoV-2 based on cycle threshold (Ct) values measured from pooled samples (Materials and Methods), potentially using far fewer tests than would be required to assess prevalence based on number of positive samples identified. We used our synthetic viral load data to assess inferential accuracy under a range of sample availabilities and pooling designs. Across the spectrum of simulated pools and tests we found that simple pooling allows accurate estimates of prevalence across at least four orders of magnitude, ranging from 0.02% to 20%, with up to 400-times efficiency gains (i.e., 400 times fewer tests) than would be needed without pooling (Fig. 3). For example, in a population prevalence study that collects ∼2,000 samples, we accurately estimated infection prevalences as low as 0.05% by using only 24 total qPCR tests (i.e., 24 pools of 96 samples each; Fig. 3A; Fig. S1). Importantly, because the distribution of Ct values may differ depending on the sample type (sputum vs. swab), the instrument, and the phase of the epidemic (growth vs. decline, Fig. S2), in practice, the method should be calibrated to viral load data (i.e., Ct values) specific to the laboratory and instrument (which can differ from one laboratory to the next) and the population under investigation.
Estimation error arises in two stages: sample collection effects, and as part of the inference method (Fig. 3B). Error from sampling collection became less important with increasing numbers of positive samples, which occurred with increasing population prevalence or by increasing the total number of tested samples (Fig. 3B; Fig. S2). At very low prevalence, small sample sizes (N) risk missing positives altogether or becoming biased by false positives. We found that accuracy in prevalence estimation is greatest when population prevalence is greater than 1/N and that when this condition was met, partitioning samples into more pools always improves accuracy (Fig. S2). In summary, very accurate estimates of prevalence can be attained using only a small fraction of the tests that would be needed in the absence of pooling.
Pooled testing for individual identification
We next analyzed effectiveness of group testing for identifying individual sample results at different points along the epidemic curve with the aim of identifying simple, efficient pooling strategies that are robust to a range of prevalences (Fig. 1A&B). Using the simulated viral load data described in Materials and Methods, we evaluated a large array of pooling designs in silico (Table S1). Based on our models of viral kinetics and given a PCR limit of detection of 100 viral copies per ml, we first estimated a baseline sensitivity of conventional (non-pooled) PCR testing of 85% during the epidemic growth phase (i.e., 15% of the time we sample an infected individual with a viral load greater than 1 but below the LOD of 100 viral copies per mL, Fig. 4A), which largely agrees with reported estimates. (21,22) This reflects sampling during the latent period of the virus (after infection but prior to significant viral growth) or in the relatively long duration of low viral titers during viral clearance.
Sensitivity of pooled tests, relative to individual testing, is affected by the dilution factor of pooling and by the population prevalence – with lower prevalence resulting in generally lower sensitivity as positives are diluted into many negatives (Fig. 4A). The decrease in sensitivity is roughly linear with the log of the dilution factor employed, which largely depends on the number and size of the pools and, for combinatorial pooling, the number of pools that each sample is placed into (Fig. S3A-C).
There is a less intuitive relationship between sensitivity and prevalence as it changes over the course of the epidemic. Early in an epidemic there is an initial dip in sensitivity for both individual and pooled testing (Fig. 4A). Early during exponential growth of an outbreak, a random sample of infected individuals will be sampled closer to their peak viral load, while later on there is an increasing mixture of newly infected with individuals with lower viral loads at the tail end of their infection. We found that this means at peak prevalence, sensitivity of pooled testing increases as samples with lower viral loads, which would otherwise be missed due to dilution, are more likely to be ‘rescued’ by coexisting in the same pool with high viral load samples and thus get individually retested (at their undiluted concentration) during the validation stage. During epidemic decline, fewer new infections arise over time and therefore a randomly selected infected individual is more likely to be sampled during the recovery phase of their infection, when viral loads are lower (Fig. S4D). Overall sensitivity is therefore lower during epidemic decline, as more infected individuals have viral loads below the limit of detection; during epidemic growth (up to day 108), overall sensitivity of RT-PCR for individual testing is 85%, whereas during epidemic decline (from day 168 onward) it is 60% (Fig. S5A). Sensitivity of RT-PCR for individual testing was ∼75% across the whole epidemic. We note that in practice, sensitivity is likely higher than estimated here, because individuals are not sampled entirely at random. Together, these results describe how sensitivity is affected by the combination of epidemic dynamics, viral kinetics, and pooling design when individuals are sampled randomly from the population.
We find that on average the majority of false negatives arise from individuals sampled seven days or more after their peak viral loads, or around seven days after what is normally considered symptom onset (∼75% in swab samples during epidemic growth; ∼96% in swab samples during epidemic decline; ∼68% in sputum during epidemic growth). Importantly, only ∼3% of false negative swab samples arose from individuals tested during the first week following peak viral load during epidemic growth, and only ∼1% during epidemic decline – (peak titers usually coincide with symptom onset) – and thus most false negatives are from individuals with the least risk of onward transmission (Fig. S3D&E).
As mentioned above, the lower sensitivity of pooled testing is counterbalanced by gains in efficiency. When prevalence is low, efficiency is roughly the number of samples divided by the number of pools, since there are rarely putative positives to test individually. However, the number of validation tests required will increase as prevalence increases, and designs that are initially more efficient will lose efficiency (Fig. 4B). In general, we find that at very low population prevalence the use of fewer pools each with larger numbers of specimens offers relative efficiency gains compared to larger numbers of pools, as the majority of pools will test negative. However, as prevalence increases, testing a greater number of smaller pools pays off as more validations will be performed on fewer samples overall (Fig. 4B). For combinatorial designs with a given number of total samples and pools, splitting each sample across fewer pools results in a modest efficiency gains (dashed versus solid lines in Fig.4B).
To address realistic resource constraints, we integrated our analyses of sensitivity and efficiency with limits on daily sample collection and testing capacity to maximize the number of positive individuals identified (see Materials and Methods). We analyzed the total number of samples screened and the fold increase in the number of positive samples identified relative to individual testing for a wide array of pooling designs evaluated over a period of 50 days during epidemic spread (days 40-90 where point prevalence reaches ∼2.5%; Fig. 4C&D). Because prevalence changes over time, the number of validation tests may vary each day despite constant pooling strategies. Thus, tests saved on days requiring fewer validation tests can be stored for days where more validation tests are required.
Across all resource constraints considered, we found that effectiveness ranged from one (when testing every sample individually is optimal) to 20 (i.e., identifying 20x more positive samples on a daily basis compared with individual testing within the same budget; Fig. 4D). As expected, when capacity to collect samples exceeds capacity to test, group testing becomes increasingly effective. Simple pooling designs are most effective when samples are in slight excess of testing capacity (2-8x), whereas we find that increasingly complex combinatorial designs become the most effective when the number of samples to be tested greatly exceeds testing capacity. Additionally, when prevalence is higher (i.e., sample prevalence from 1.03% to 9.90%), the optimal pooling designs shift towards combinatorial pooling, and the overall effectiveness decreases – but still remains up to 4x more effective than individual testing (Fig. S6). Our results were qualitatively unchanged when evaluating the effectiveness of pooling sputum samples, and the optimal pooling designs under each set of sample constraints were either the same or very similar (Fig. S7). Furthermore, we evaluated the same strategies during a 50-day window of epidemic decline (days 190-250) and found that similar pooling strategies were optimally effective, despite lower overall sensitivity as described above (Fig. S5).
Pilot and validation experiments
We validated our pooling strategies using anonymized clinical nasopharyngeal swab specimens. To evaluate simple pooling across a range of inputs, we diluted 5 nasopharyngeal clinical swab samples with viral loads of 89000, 12300, 1280, 140 and 11 viral copies per ml, respectively, into 23 negative nasopharyngeal swab samples (pools of 24). Further details are provided in Materials and Methods, and Supplementary Material 1, sections 7&8. The results matched the simulated sampling results: the first three pools were all positive, the fourth was inconclusive (negative on N1, positive on N2), and the remaining pool was negative (Fig. 5A, Table S2). These results are as expected because the EUA approved assay used has a limit of detection of ∼100 virus copies per ml, such that the last two specimens fall below the limit of detection given a dilution factor of 24 (i.e. 0.46 and 5.8 virus copies per ml once pooled).
We next tested combinatorial pooling, first using only a modest pooling design. We split 48 samples, including 1 positive, into 6 pools with each sample spread across three different pools. The method correctly identified the three pools containing the positive specimen (Fig. 5B, Table S2). One negative sample was included in the same 3 pools as the positive sample; thus, 8 total tests (6 pools + 2 validations) were needed to accurately identify the status of all 48 samples, a 6x efficiency gain, which matched our expectations from the simulations.
We next performed two larger validation studies (Materials and Methods, and Supplementary Material 1, section 8). To validate combinatorial pooling, we used anonymized samples representing 930 negative and 30 distinct positive specimens (3.1% prevalence), split across 10 batches of 96 specimens each (Table S3). For each batch of 96, we split the specimens into 6 pools and each specimen was spread across 2 pools (Fig. 5C, Table S4). For this combinatorial pooling design and prevalence, our simulations suggest that we would expect to identify 26 out of 30 known positives (87%) and would see a 2.81x efficiency gain – using only 35% of the number of tests compared to no pooling. We identified 24 of the 30 known positives (80%) and, indeed, required 35% fewer tests (341 vs 960, a 2.8x efficiency gain).
To further validate our methods for prevalence estimation, we created a large study representing 2,304 samples with a (true) positive prevalence of 1%. We aimed to determine how well our methods would work to estimate the true prevalence using 1/48th the number of tests compared to testing samples individually. To do this, we randomly assigned 24 distinct positive samples into 48 pools, with each pool containing 48 samples (Table S3; to create the full set of pools, we treated some known negatives as distinct samples across separate pools). By using the measured viral loads detected in each of the pools, our methods estimated a prevalence of 0.87% (compared to the true prevalence of 1%) with a bootstrapped 95% confidence interval of 0.52% - 1.37% (Fig. 5D), and did so using 48x fewer tests than without pooling. This level of accuracy is in line with our expectations from our simulations. Notably, the inference algorithm applied to these data used viral load distributions calibrated from our simulated epidemic, which in turn had viral kinetics calibrated to samples collected and tested on another continent, demonstrating robustness of the training procedure.
Discussion
Our results show that group testing for SARS-CoV-2 can be a highly effective tool to increase surveillance coverage and capacity when resources are constrained. For prevalence testing, we find that fewer than 40 tests can be used to accurately infer prevalences across four orders of magnitude, providing large savings on tests required. For individual identification, we determined an array of designs that optimize the rate at which infected individuals are identified under constraints on sample collection and daily test capacities. These results provide pooling designs that maximize the number of positive individuals identified on a daily basis, while accounting for epidemic dynamics, viral kinetics, viral loads measured from nasopharyngeal swabs or sputum, and practical considerations of laboratory capacity.
While our experiments suggest that pooling designs may be beneficial for SARS-CoV-2 surveillance and identification of individual specimens, there are substantial logistical challenges to implementing theoretically optimized pooling designs. Large-scale testing without the use of pooling already requires managing thousands of specimens per day, largely in series. Pooling adds complexity because samples must be tracked across multiple pools and stored for potential re-testing. These complexities can be overcome with proper tracking software (including simple spreadsheets) and standard operating procedures in place before pooling begins. Such procedures can mitigate the risk of handling error or specimen mix-up. In addition, expecting laboratories to regularly adapt their workflow and optimize pool sizes based on prevalence may not be feasible in some settings. (8,16) A potential solution is to follow a simple, fixed protocol that is robust to a range of prevalences. Supplementary Material 2 provides an example spreadsheet guiding a technician receiving 96 labeled samples to create 6 pools, enter the result of each pool and be provided a list of putative positives to be retested.
Enhancing efficiency at the expense of sensitivity must be considered depending on the purpose of testing. For prevalence testing, accurate estimates can be obtained using very few tests if individual identification is not the aim. For individual testing, although identifying all positive samples that are tested is of course the objective, increasing the number of specimens tested when sacrificing sensitivity may be a crucially important tradeoff. This tradeoff is particularly pertinent because the specimens most likely to be lost due to dilution are those samples with the lowest viral loads already near the limit of detection (Fig.S3D&E). Although there is a chance that the low viral load samples missed are on the upswing of an infection – when identifying the individual would be maximally beneficial – the asymmetric course of viral titers over the full duration of positivity means that most false negatives would arise from failure to detect late-stage, low-titer individuals who are less likely to be infectious. (20) Optimal strategies and expectations of sensitivity should also be considered alongside the phase of the epidemic and how samples are collected, as this will dictate the distribution of sampled viral loads. For example, if individuals are under a regular testing regimen or are tested due to recent exposure or symptom onset, then viral loads at the time of sampling will typically be higher, leading to higher sensitivity in spite of dilution effects.
Testing throughput and staffing resources should also be considered. If a testing facility can only run a limited number of tests per day, it may be preferable to process more samples at a slight cost to sensitivity. Back-logs of individual testing can result in substantial delays in returning individual test results, which can ultimately defeat the purpose of identifying individuals for isolation - potentially further justifying some sensitivity losses. (23,24) Choosing a pooling strategy will therefore depend on target population and availability of resources. For testing in the community or in existing sentinel surveillance populations (e.g., antenatal clinics), point prevalence is likely to be low (<0.1 - 3%), which may favor strategies with fewer pools. (6,25–27) Conversely, secondary attack rates in contacts of index cases may vary from <1% to 17% depending on the setting (e.g., casual vs. household contacts), (28–30) and may be even higher in some instances, (31,32) favoring more pools. These high point prevalence sub-populations may represent less efficient use cases for pooled testing, as our results suggest that pooling for individual identification is inefficient once prevalence reaches 10%. However, group testing may still be useful if testing capacity is severely limited – for example, samples from all members of a household could be tested as a pool and quarantined if tested positive, enabling faster turnaround than testing individuals. This approach may be even more efficient if samples can be pooled at the point of collection, requiring no change to laboratory protocols.
Our modelling results have a number of limitations and may be updated as more data become available. First, our simulation results depend on the generalizability of the simulated Ct values, which were based on viral load data from symptomatic patients. Although some features of viral trajectories, such as viral waning, differ between symptomatic and asymptomatic individuals, population-wide data suggest that the range of Ct values do not differ based on symptom status. (20,33) Furthermore, we have assumed a simple hinge function to describe viral kinetics. Different shapes for the viral kinetics trajectory may become apparent as more data become available. Nonetheless, our simulated population distribution of Ct values is comparable to existing data and we propagated substantial uncertainty in viral kinetics parameters to generate a wide range of viral trajectories. For prevalence estimation, the MLE framework requires training on a distribution of Ct values. Such data can be available based on past tests from a given laboratory, but care should be taken to use a distribution appropriate for the population under consideration. For example, training the virus kinetics model on data skewed towards lower viral loads (as would be observed during the tail end of an epidemic curve) may be inappropriate when the true viral load distribution is skewed higher (as might be the case during the growth phase of an epidemic curve). Nevertheless, we used our simulated distribution of Ct, which were fit to virus kinetics in published reports from distinct labs from across the world and obtained highly accurate results throughout. Thus, despite the limitations just mentioned, this shows that the virus kinetics models are quite robust and may not, in practice, require new fitting to the individual laboratory or population. In addition, while we assume that individuals are sampled from the population at random in our analysis, in practice samples that are processed together are also typically collected together, which may bias the distribution of positive samples among pools.
We have shown that simple designs that are straightforward to implement have the potential to greatly improve testing throughput across the time course of the pandemic. These principles likely also hold for pooling of sera for antibody testing, which remains an avenue for future work. There are logistical challenges and additional costs associated with pooling that we do not consider deeply here, and it will therefore be up to laboratories and policy makers to decide where these designs are feasible. Substantial coordination will therefore be necessary to make group testing practical but investing in these efforts could enable community screening where it is currently infeasible and provide epidemiological insights that are urgently needed.
Materials and Methods
Simulation model of infection dynamics and viral load kinetics
We developed a population-level mathematical model of SARS-CoV-2 transmission that incorporates realistic within-host virus kinetics. Full details are provided in Supplementary Material 1, sections 1-4, but we provide an overview here. First, we fit a viral kinetics model to published longitudinally collected viral load data from nasopharyngeal swab and sputum samples using a Bayesian hierarchical model that captures the variation of peak viral loads, delays from infection to peak and virus decline rates across infected individuals (Fig. 5A&B; Fig. S8). (19) By incorporating estimated biological variation in virus kinetics, this model allows random draws each representing distinct within-host virus trajectories. We then simulated infection prevalence during a SARS-CoV-2 outbreak using a deterministic Susceptible-Exposed-Infected-Removed (SEIR) model with parameters reflecting the natural history of SARS-CoV-2 (Fig. 5D). For each simulated infection, we generated longitudinal virus titers over time by drawing from the distribution of fitted virus kinetic curves, using distributions derived using either nasopharyngeal swab or sputum data (Fig. S4). All estimated and assumed model parameters are shown in Table S5, with model fits shown in Fig. S8. Posterior estimates and Markov chain Monte Carlo trace plots are shown in Fig. S9 and Fig. S10. We accounted for measurement variation by: i) transforming viral loads into Ct values under a range of Ct calibration curves, ii) simulating false positives with 1% probability, and, importantly, iii) simulating sampling variation. We assumed a limit of detection (LOD) of 100 RNA copies / ml.
Estimating prevalence from pooled test results
We adapted a statistical (maximum likelihood) framework initially developed to estimate HIV prevalence with pooled antibody tests to estimate prevalence of SARS-CoV-2 using pooled samples. (12,13) The framework accounts for the distribution of viral loads (and uncertainty around them) measured in pools containing a mixture of negative and potentially positive samples. By measuring viral loads from multiple such pools, it is possible to estimate the prevalence of positive samples without individual testing. See Supplementary Material 1, section 5 for full details.
We evaluated prevalence estimation under a range of sample availabilities (N total samples; N=288 to ∼18,000) and pooling designs. We varied the pool size of combined specimens (n samples per pool; n=48, 96, 192, or 384) and the number of pools (b=6, 12, 24, or 48). For each combination, we estimated the point prevalence from pooled tests on random samples of individuals drawn during epidemic growth (days 20-120) and decline (days 155-300). Because the data is realistic but simulated, we used ground truth prevalence in the population and, separately, in the specific set of samples collected from the overall population to assess accuracy of our estimates (see for example Fig. 3B). We calculated estimates for 100 entirely distinct epidemic simulations.
Pooled tests for individual sample identification
Using the same simulated population, we evaluated a range of simple and combinatorial pooling strategies for individual positive sample identification (Supplementary Material 1, section 6). In simple pooling designs, each sample is placed in one pool, and each pool consists of some pre-specified number of samples. If a pool tests positive, all samples that were placed in that pool are retested individually (Fig. 1A). For combinatorial pooling, each sample is split into multiple, partially overlapping pools (Fig. 1B). (9,10) Every sample that was placed in any pool that tested negative is inferred to be negative, and the remaining samples are identified as potential positives. Here, we consider a very simple form of combinatorial testing, where identified potential positive samples are individually tested in a validation stage.
A given pooling design is defined by three parameters: the total number of individuals to be tested (N); the total number of pools to test (b); and the number of pools a given sample is included in (q). For instance, if we have 50 individuals (N) to test, we might split the 50 samples into four pools (b) of 25 samples each, where each sample is included in two pools. Note that, by definition, in simple pooling designs each sample is placed in one pool (q=1).
To identify optimal testing designs under different resource constraints, we systematically analyzed a large array of pooling designs under various sample and test kit availabilities. We evaluated different combinations of between 12 and ∼6,000 available samples/tests per day. The daily testing capacity shown is the daily average, though we assume that there is some flexibility to use fewer or more tests day to day (i.e., there is a budget for period of time under evaluation).
For each set of resource constraints, we evaluated designs that split N samples between 1 to 96 distinct pools, and with samples included in q=1 (simple pooling), 2, 3, or 4 (combinatorial pooling) pools (Table S1). To ensure robust estimates (especially at low prevalences of less than 1 in 10,000), we repeated each simulated pooling protocol at each time point in the epidemic up to 200,000 times.
In each scenario, we calculated: i) the sensitivity to detect positive samples when they existed in the pool; ii) the efficiency, defined as the total number of samples tested divided by the total number of tests used; iii) the total number of identified true positives (total recall); and iv) the effectiveness, defined as the total recall relative to individual testing.
Pilot experiments
For validation experiments of our simulation-based approach, we used fully de-identified, discarded human nasopharyngeal specimens obtained from the Broad Institute of MIT and Harvard. In each experiment, sample aliquots were pooled before RNA extraction and qPCR and pooled specimens were tested using the Emergency Use Authorization (EUA) approved SARS-CoV-2 assay performed by the Broad Institute CLIA laboratory. The protocol is described in full detail in Supplementary Material 1, section 7. Specifics of each pooling approach is detailed alongside the relevant results and in Supplementary Material 1, section 8.
Data Availability
Raw data generated in this study are available in Table S2 and Table S6. Simulated data can be regenerated using the accompanying code.
Author contributions
B.C., J.H., A.R. and M.M. conceived the study. B.C. and J.H. performed the simulations and modelling. B.B., M.H, M.C, J.B, and B.S performed the pooled testing. B.C., J.H, A.R, and M.M. analyzed results and wrote the manuscript with input from all authors.
Competing interests
A.R. is a co-founder and equity holder of Celsius Therapeutics, an equity holder in Immunitas, and until July 31, 2020, was an SAB member of ThermoFisher Scientific, Syros Pharmaceuticals, Asimov, and Neogene Therapeutics. From August 1, 2020, A.R. is an employee of Genentech.
Code availability
The SEIR model, viral kinetics model and MCMC were implemented in R version 3.6.2. The remainder of the work was performed in Python version 3.7. The code used for the MCMC framework is available at: https://github.com/jameshay218/lazymcmc. All other code and data used are available at: https://github.com/cleary-lab/covid19-group-tests.
Data availability
Raw data generated in this study are available in Table S2 and Table S6. Simulated data can be regenerated using the accompanying code.
Supplementary Material 1: additional materials and methods
1. Model of infection dynamics
We implemented a deterministic compartmental model to describe the increase and subsequent decline of SARS-CoV-2 infection incidence and prevalence. The model captured the following progression: individuals begin fully susceptible to infection (S); become exposed but not yet infectious (E); become infectious (I); and are finally removed (R). Our aim was to capture changes in incidence over time (increasing incidence rate before the peak and declining incidence thereafter) and the resulting population-level distribution of viral loads, and we therefore did not model additional complexity such as deaths, pre-symptomatic transmission or age-stratified outcomes. Although incorporating these mechanisms may alter the simulated epidemic dynamics for a given set of transmission parameters, we do not expect them to impact inference from the pooled testing analyses.
The model was defined by the following set of ordinary differential equations: where β is the transmission rate, scaled to give a basic reproductive number R0 of 2.5; 1/σ is the incubation period, assumed to have a mean of 6.4 days; 1/γ is the infectious period assumed to have a mean of 7 days; and N is the population size, set to 12,500,000 to generate at least 10,000,000 infected individuals. All model parameters and assumed values are described in Table S5. Using this model, we simulated per capita daily infection probability (i.e., the daily incidence) and prevalence.
2. Viral load kinetics
Following infection (tinc) and after a latent period (tg), viral load was assumed to increase exponentially, peaking a short time (tp) after the end of the latent period and before the onset of symptoms (to). Viral load then decreased monotonically to undetectable levels, defined by a time parameter (tw) giving the number of days post symptom onset before crossing the limit of detection. These assumptions are equivalent to linear growth and decay on a log scale. We chose to model viral load waning with respect to symptom onset rather than peak viral load, as almost all available time-series viral load data are presented with respect to symptom onset time. The log10 viral load over time, v(t), was given by: where tinf is the time of infection; tinc is the incubation period for symptom onset; α is the peak viral load; tg is the latent period before viral growth; tp is the time taken to reach peak viral load post; and tw is the number of days from symptom onset to becoming undetectable. Note that all individuals are assigned a symptom onset time regardless of whether they show symptoms or not, as the parameter is used here to describe viral kinetics.
3. Fitting the viral kinetics model
Time-series viral load data were obtained from a case series of 9 hospitalized COVID-19 patients from a single hospital in Munich, Germany. (34) These data provide regular measurements of log10 RNA copies per nasopharyngeal swab and per ml of sputum. Data were extracted using a web plot digitizer (https://automeris.io/WebPlotDigitizer/). We used a random-effects model to infer individual-level viral kinetics parameters alongside population-level distributions. Under this model, each individual had their own parameter values for tinc, α, tg, tp and tw. α and tw drawn from a multivariate normal distributed with estimated means , standard deviations and correlation (Table S5).
Although the viral load measurements were taken after the onset of symptoms, we are interested in modelling the entire time course from the time of infection, including uncertainty in the incubation period. We therefore estimated tinc for each individual by placing a strongly informative log-normal prior on the incubation period based on previous estimates. (35) Because no observations were made before symptom onset, the posterior distribution for this parameter matched the prior.
The model was fit separately to the swab and sputum data using a Markov chain Monte Carlo (MCMC) framework, using a likelihood that assumes normally distributed observation error with estimated standard deviation σobs (Fig. S8). The likelihood function accounted for truncation at the lower limit of detection (LOD) by integrating over the left-hand tail of the normal cumulative density function for model-predicted viral loads less than the LOD. Furthermore, Wölfel et al. present data with a limit of detection of 0 log10 copies / swab or / ml, but a limit of quantification of 2 log10 copies. We therefore considered that observed viral loads between 0 and 2 had unknown between 0 and 2. We placed uniform priors on and , half Cauchy priors on σα and with scale parameters of 1, a uniform prior on between −1 and 1 and informative uniform priors on tg, and tp.
Code to regenerate the MCMC fitting is available in the accompanying GitHub repository. (36) Briefly, an adaptive Metropolis-Hastings algorithm was run for 3,000,000 iterations with the first 1,000,000 iterations discarded as burn in to obtain posterior estimates for the model parameters (Fig. S9). Convergence was assessed visually based on trace plots of 3 independent chains (Fig. S10). All estimated parameters had an effective sample size (ESS) greater than 200 and .
4. Simulation framework
We combined the estimates from the transmission model and viral kinetics model above to simulate viral loads for a population of 12,500,000 individuals over time. We performed separate simulations for swab and sputum data. For each individual, we simulated a time of infection (or no infection) from the time-varying probability of infection, p(t). For each infected individual, an incubation period was drawn from a log normal distribution with mean (of the natural log) = 1.621 and standard deviation (of the natural log) = 0.418 as estimated previously. (35) This incubation period was used as a reference point for tw, but we did not explicitly model symptomatic vs. asymptomatic individuals. Next, we drew viral kinetics parameters for each individual from the joint posterior distribution. Finally, we simulated the time-varying viral load of each infected individual, starting at the time of infection. Because of the high variance of the drawn kinetics parameters, some simulated viral loads reached very high values, and we therefore truncated all simulated viral loads at 16 log10 virus copies per ml. These simulations provide a population of viral loads incorporating realistic variability in viral load kinetics and between-individual parameters.
5. Estimating prevalence from pooled test results
In this section, we adapt a previously described maximum likelihood estimate (MLE) framework to estimate prevalence of positive samples, p, amongst N samples using b pooled tests, with each pooling containing . (12,13) Briefly, the MLE is determined by computing the conditional probability of an observed RT-qPCR result given that there were k positive samples in the pool, and integrating over all values of k. To compute these conditional distributions, we estimate an empirical distribution of cycle threshold (Ct) values that follows from the distribution of viral loads in infected individuals derived from a subset of the simulated viral load measurements described above. This is conceptually similar to using the empirical distribution of all real measurements from a population taken to date.
Let fk,n = P(Y| k, n) represent the conditional probability of the observed Ct values in each pool, Y = [y1, …, yb], given that there were k positive samples in the pool and n-k negative samples. Computing these distributions will be the central task of generating prevalence estimates.
To begin, we model Ct values as a function of the underlying viral loads: where v is the viral load (in copies per ml). To capture variance in the LOD between labs or across batches within a lab, we let the intercept, y0, be normally distributed with unknown mean and standard deviation. We assume uniform priors over each y0 (μ ∼ Uniform(37,39), σ ∼ Uniform(0.1,1)). These ranges and the log-linear relationship are based on a calibration of Ct values to a dilution curve of positive control SARS-CoV-2 RNA samples (data not shown) and are consistent with prior measurements calibrated to SARS coronavirus. (37)
fk,n depends not only on the distribution of Ct values arising from a given viral load, but also on the distribution of viral loads in infected and uninfected individuals. Here we assume that the viral load in uninfected individuals is always zero (we do, however, allow for false positive PCR results, as discussed below). Given access to the true distribution of viral loads, fk,n can be generated by the convolution of k density functions. Since we do not have access to the true distributions, we approximate fk,n through a kernel density estimate (KDE) of empirical convolutions. Specifically, for a given value of n and each value of k from 1 to n, we generate 10,000 random combinations of k viral loads, where zj are the randomly sampled viral loads. For each i, we sample μ ∼ Uniform(37,39), σ ∼ Uniform(0.1,1) and y0 ∼ N(μ, σ) and then generate Ct values:
We then approximate fk,n = KDE(Xk,n) with a bandwidth parameter of 0.1.
Several adjustments to these density approximations are needed to account for false positives and undefined Ct values (i.e., undetected viral RNA). To allow for Ct values that might arise from false positives, we define f0,n = f1,n r, where r is an assumed false positive rate of PCR (here, we set r = 0.2% and allow this to be misspecified when simulating false positives below). When viral RNA is undetected and k > 1, we let where a Ct value of 40 is used as a typical limit of detectable RNA. When k =0, f0,n(undetected) = 1 − r.
Following the model of Zenios and Wein, (12) the MLE of prevalence, p, is defined as: where n is the number of samples per pool, b, is the number of pools, and yi is the Ct value observed in pool i. We calculate P(k | yi, p) using the conditional densities above: substituting the appropriate values of f when the Ct value is undefined. Note that each of the observed Ct values is adjusted to account for n-fold dilution. Finally, we find pestimated through an iterative expectation-maximization algorithm, where at each iteration
Using this model, we evaluated prevalence estimation across time in our simulated populations. We first partition the simulated population into training and testing sets. Observed viral loads in 20% of individuals (selected at random) are used to train the KDEs above. Furthermore, training data were divided into growth phase (days 20 to 120) or decline phase (days 155 to 300) data, allowing us to ascertain the effect of training and estimating on data from consistent or inconsistent phases of the epidemic. All remaining prevalence estimation and analyses are done on the testing population. We simulate viral loads in b pools of samples as for i from 1 to b, where Poisson sampling is used to model the sampling of a small volume from each swab. From these we simulate Ct values, Y. At each time point, we randomly select N individuals, partition them into b pools, and simulate the pooled viral loads and Ct values, Y, with y0 ∼ N(38.5,1) sampled randomly in each trial and applying the transformation described previously. To capture false positive PCR results at a rate of 1%, with 1% probability in each pool we add a viral load of , with zj selected at random from all positive samples with on a given date. We then compute the true prevalence in the population, ptrue, the prevalence in the population of N samples, psampled, and the prevalence estimated from b tests, pestimated. To quantify accuracy we calculate the mean and standard deviation of these values at each time point across 100 random trials and summarize differences between pestimated.and either ptrue or psampled with the normalized root mean squared error (NRMSE; either or ; Fig. S2).
To demonstrate the importance of calibrating the correct empirical distribution of Ct values, Fig. S2 shows the accuracy of the framework in recovering the true population prevalence during the growth and decline phases of the epidemic when the model is calibrated to Ct values measured during the growth and decline phases. Estimation accuracy is highest when the calibrated Ct values were measured during the same phase as is being measured.
6. Group testing for sample identification
We evaluated an array of pooling designs for individual identification as described in the main text. For any given design we can simulate pooled testing results on the simulated population described above. For each time point, we randomly select N individuals from the population. Each individual is then randomly assigned to q out of b pools. For each pool, we then calculate the net sampled viral load, using the viral loads zj for each individual included in poolj, where Poisson sampling is used to model a small volume sampled from each swab. If vi > LOD the pool is assigned a positive value, otherwise the pool is negative (here, with LOD=100). To allow for false positive PCR results at a rate of 1%, with probability 1% we set the result of the pool to be positive, regardless of the viral load.
Simulations were run by repeating the random pooling, pooled testing, and decoding procedures described in the main text. At each point in time, 200,000 trials are run selecting N individuals at random in each trial. We then record the number of validation tests, s, in each trial, corresponding to the number of putative positive samples. Average efficiency at a point in time, t, is then calculated as . Average sensitivity is calculated as where Npositive is the total number of infected individuals (viral load >1) that were sampled across all trials and Nputative is the number of such individuals who were identified as putative positives. From these values we calculate Total Recallt =Efficiencyt Sensitivityt.
Based on these results we evaluated a large array of group testing designs (Table S1) under a series of constraints on number of samples collected and number of reactions run per day (Fig. 4C). For a given design, D: (N, b, q) we calculate the average number of tests run, stotal = b + s, and the total number of times that D could be run on a daily basis, We then calculate the effective testing capacity on a daily basis, cD = N rD SensitivityD, and average this value over days 40 to 90 (Fig. 4C) or 80 to 108 (Fig. S6) of the simulated epidemic. Notably, when sample collection is limiting, cD =samples collected ∗ Sensitivity(sample collected, sample collected,1), corresponding to individual testing. Finally, the optimal design for a given combination of constraints is the design with the greatest average effective testing capacity. In Fig. 4C and Fig. S6 we show the relative effectiveness, corresponding to effective capacity relative to individual testing. Fig. S7 shows the same results as Fig. 4 and Fig. S6, but when samples were simulated based on viral trajectories using the sputum data. Fig. S5 shows the same results again, but for the decline phase of the epidemic.
7. Broad Institute CRSP SARS-CoV-2 Real-time Reverse Transcriptase (RT)-PCR Diagnostic Assay
RNA is isolated from nasopharyngeal and oropharygneal specimens in 50ul VTM/UTM using MagMAX™-96 Viral RNA Isolation Kits (Thermo Fisher Scientific, AMB18365), performed on an Agilent Bravo Liquid Handling Platform running VWorks (Build 11.4.0.1233), and is reverse transcribed to cDNA and subsequently amplified in the Applied Biosystems® ViiA7 Real-Time PCR Instrument with QuantStudio version 1.3 software. In the process, the probe anneals to a specific target sequence located between the forward and reverse primers. During the extension phase of the PCR cycle, the 5’ nuclease activity of Taq polymerase degrades the probe, causing the reporter dye to separate from the quencher dye, generating a fluorescent signal. With each cycle, additional reporter dye molecules are cleaved from their respective probes, increasing the fluorescence intensity. Fluorescence intensity is monitored at each PCR cycle by Applied Biosystems® ViiA7 Real-Time PCR System with QuantStudio version 1.3 software. The software allows the fluorescence intensity to be monitored at each PCR cycle to allow for the qualitative detection of nucleic acid from SARS-CoV-2.
To create the final pools for our larger validation studies (10 batches of 96 for sample identification, and 48 pools of 48 samples for prevalence estimation), we treated a total of 96 known negative specimens as distinct samples across batches (in the case of sampled identification study) or pools (for the prevalence study). Each positive sample in a batch or pool was a distinct sample, and only used in one batch or pool. For the prevalence study, only pools with one or more positive samples were tested using the assay above (the remaining pools were assigned an “undefined” Ct value). All pools in all batches were tested in the identification study, regardless of whether they contained a positive sample.
8. Selecting pool compositions for large-scale sample identification validation
To form pools, we put each of the N = 96 individuals into q = 2 out of b = 6 pools (A-F) by cycling through the following ordered list of pool pairs: AB, CD, EF, BC, DF, AE, BD, AF, CE, BE, CF, AD, BF, DE, AC. Namely, individual 1 is put in pools A and B, individual 2 in C and D, and so on; after individual 15, we return to the beginning of the list and cycle through again. Finally, since each pair of pools is assigned to multiple individuals, we reordered the individuals by what pair of pools they were put in. This final reordering simplifies the presentation to one more conveniently arranged for manual pipetting.
A strength of the procedure is that it is simple and flexible; it can easily be carried out for any number of individuals. Moreover, the resulting design here has pools with 32 individuals each, so the final pool sizes are balanced. In addition, the first 6 pool pairs are each assigned to 7 individuals and the remaining 9 pairs are each assigned to 6 individuals, so the pairs of pools are also approximately balanced. Indeed, the ordered list of pool pairs above was chosen to achieve this balance; the list uses each pair of pools once and the five consecutive sets of three pairs all use each pool once.
Acknowledgements
We thank Benedicte Gnangnon, Rounak Dey, Xihong Lin, Edgar Doriban, Rene Niehus and Heather Shakerchi for useful discussions. Work was supported by the Merkin Institute Fellowship at the Broad Institute of MIT and Harvard (B.C.), by the National Institute of General Medical Sciences (#U54GM088558; J.H. and M.M.); and by an NIH DP5 grant (M.M.), and the Dean’s Fund for Postdoctoral Research of the Wharton School and NSF BIGDATA grant IIS 1837992 (D.H.).