In late December of 2019, the first cases of COVID-19, the disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), were described in the city of Wuhan in Hubei province, China (
1,
2). The virus quickly spread within China (
3). The cordon sanitaire that was put in place in Wuhan on 23 January 2020 and mitigation efforts across China eventually brought about an end to sustained local transmission. In March and April 2020, restrictions across China were relaxed (
4). By then, however, COVID-19 was a pandemic (
5).
A concerted effort has been made to retrospectively diagnose the earliest cases of COVID-19 and thus determine when the virus first began transmitting among humans. Both epidemiological and phylogenetic approaches suggest an emergence of the pandemic in Hubei province at some point in late 2019 (
2,
6,
7). The first described cluster of COVID-19 was associated with the Huanan Seafood Wholesale Market in late December 2019, and the earliest sequenced SARS-CoV-2 genomes came from this cluster (
8,
9). However, this market cluster is unlikely to have denoted the beginning of the pandemic, as COVID-19 cases from early December lacked connections to the market (
7). The earliest such case in the scientific literature is from an individual retrospectively diagnosed on 1 December 2019 (
6). Notably, however, newspaper reports document retrospective COVID-19 diagnoses recorded by the Chinese government going back to 17 November 2019 in Hubei province (
10). These reports detail daily retrospective COVID-19 diagnoses through the end of November, suggesting that SARS-CoV-2 was actively circulating for at least a month before it was discovered.
Molecular clock phylogenetic analyses have inferred the time of the most recent common ancestor (tMRCA) of all sequenced SARS-CoV-2 genomes to be in late November or early December 2019, with uncertainty estimates typically dating to October 2019 (
7,
11,
12). Crucially, though, this tMRCA is not necessarily equivalent to the date of zoonosis or index case infection (
13,
14) because coalescent processes can prune basal viral lineages before they have the opportunity to be sampled, potentially pushing SARS-CoV-2 tMRCA estimates forward in time from the index case by days, weeks, or months. For a point of comparison, consider the zoonotic origins of the HIV-1 pandemic, whose tMRCA in the early 20th century coincided with the urbanization of Kinshasa, in what is now the Democratic Republic of the Congo (
15,
16), but whose cross-species transmission from a chimpanzee reservoir occurred in southeast Cameroon, likely predating the tMRCA of sampled HIV-1 genomes by many years (
17). Despite this important distinction, the tMRCA has been frequently conflated with the date of the index case infection in the SARS-CoV-2 literature (
7,
18,
19).
Here, we combine retrospective molecular clock analysis in a coalescent framework with a forward compartmental epidemiological model to estimate the timing of the SARS-CoV-2 index case in Hubei province. The inferred dynamics during these unobserved early days of SARS-CoV-2 highlight challenges in detecting and preventing nascent pandemics.
We first explored the evolutionary dynamics of the first wave of SARS-CoV-2 infections in China. We used Bayesian phylodynamics (
20) to reconstruct the underlying coalescent processes using a Bayesian Skyline approach for 583 SARS-CoV-2 complete genomes, sampled in China between when the virus was first discovered at the end of December 2019 and the last of the non-reintroduced circulating virus in April 2020. Applying a strict molecular clock, we inferred an evolutionary rate of 7.90 × 10
−4 substitutions per site per year [95% highest posterior density (HPD): 6.64 × 10
−4 to 9.27 × 10
−4]. The tMRCA of these circulating strains was inferred to fall within a 34-day window with a mean of 9 December 2019 (95% HPD: 17 November to 20 December) (
Fig. 1). This estimate accounts for the many disparate inferred rooting orientations [see supplementary text and (
21)]. Notably, 78.7% of the posterior density postdates the earliest published case on 1 December, and 95.1% postdates the earliest reported case on 17 November. Relaxing the molecular clock provides a similar tMRCA estimate, as does applying a Skygrid coalescent approach (fig. S1). The recency of this tMRCA estimate in relation to the earliest documented COVID-19 cases obliges us to consider the possibility that this tMRCA does not capture the index case and that SARS-CoV-2 was circulating in Hubei province before the inferred tMRCA.
If the tMRCA postdates the earliest documented cases, then the earliest diverged SARS-CoV-2 lineages must have gone extinct (
Fig. 2). As these early basal lineages disappeared, the tMRCA of the remaining lineages would move forward in time (fig. S2). Thus, we interrogated the posterior trees sampled from the phylodynamic analysis to determine whether this time of coalescence had stabilized before the sequencing of the first SARS-CoV-2 genomes on 24 December 2019 or whether this process of basal lineage loss was ongoing in late December and/or early January. Notably, these basal lineages need not be associated with specific mutations, as the phylodynamic inference reconstructs the coalescent history, not the mutational history (
20).
We find only weak evidence for basal lineage loss between 24 December 2019 and 13 January 2020 (fig. S3A). The root tMRCA is within 1 day of the tMRCA of virus sampled on or after 1 January 2020 in 78.5% of posterior samples (fig. S3B). The tMRCA of genomes sampled on or after 1 January 2020 is 3 days later than the tMRCA of all sampled genomes. By contrast, the mean tMRCA does not change when considering genomes sampled on or after 1 January 2020 versus on or after 13 January 2020. This consistency indicates a stabilization of coalescent processes at the start of 2020, when an estimated total of 1000 people had been infected with SARS-CoV-2 in Wuhan (
22). Nonetheless, to account for the weak signal of a delay in reaching a stable coalescence (i.e., the point in time at which basal lineages cease to be lost), we identified the tMRCA for all viruses sampled on or after 1 January 2020 (i.e., at the time of stable coalescence) for each tree in the posterior sample.
Phylogenetic analysis alone cannot tell us how long SARS-CoV-2 could have circulated in Hubei province before the tMRCA. To answer this question, we performed forward epidemic simulations (
23). These simulations were initiated by a single index case using a compartmental epidemiological model across scale-free contact networks (mean number of contacts: 16). This compartmental model was previously developed to describe SARS-CoV-2 transmission dynamics in Wuhan (
22). This model, termed SAPHIRE, includes compartments for susceptible (S), exposed (E), presymptomatic (P), unascertained (A), ascertained (I), hospitalized (H), and removed (R) individuals. Our simulations used parameters from the time period before COVID-19 mitigation efforts, from 1 January through 22 January 2020 (table S1), based on the work of Hao
et al. (
22). We analyzed 1000 epidemic simulations that resulted in ≥1000 total infected people. These simulated epidemics had a median doubling time of 4.1 days (95% range across simulations: 2.7 to 6.7), matching premitigation incidence trends in Wuhan (table S2).
We simulated coalescent processes across the transmission network to determine the tMRCA of the virus at the end of the simulation. This approach allowed us to determine the distribution of the expected number of days between index case infection and the stable coalescence (i.e., tMRCA) (
Fig. 2). The median number of days between index case infection and this tMRCA was 8.0 days (95% range: 0.0 to 41.5 days) (
Fig. 3A). The median time between index case infection and the first person exiting the presymptomatic phase (i.e., ascertained or unascertained infection) was 5.7 days (95% range: 0.9 to 15.7 days).
As a robustness check, we also simulated epidemics with more densely and more sparsely connected contact networks (mean: 26 and 10 contacts, respectively), rescaling the per-contact transmission rate to maintain empirical epidemic growth dynamics. We also explored the effects of faster (mean: 3.1 days; 95% range: 2.0 to 5.1 days) and slower (mean: 5.3 days; 95% range: 3.6 to 7.5 days) epidemic doubling times (table S2). Slower transmission rates led to more days between the index case and the stable coalescence, but modifying the density of the contact network had minimal effect on this time interval (fig. S4).
To estimate the date of infection for the index case in Hubei province, we combined the retrospective molecular clock analysis with the forward epidemic simulations (fig. S5). We identified the stable tMRCA in the posterior trees as an anchor to the real-world calendar dates and then extended this date back in time according to the number of days between the index case infection and the time of stable coalescence from the compartmental epidemic simulations. However, a random sample of tMRCAs and days from index case infection to coalescence will not produce epidemiologically meaningful results because many of these combinations do not precede the earliest dates of reported COVID-19 cases. Therefore, we implemented a rejection sampling–based approach to generate a posterior distribution of dates of infection for the Hubei index case, conditioning on at least one individual who had progressed past the presymptomatic stage in the simulated epidemic before the date of the first reported COVID-19 case (see materials and methods and fig. S6).
In our primary analysis, we assume that 17 November represents the first documented case of COVID-19 (ascertained or unascertained in the SAPHIRE model). Under this assumption, the median number of days between index case infection and stable coalescence after rejection sampling is 37 days (95% HPD: 12 to 55 days) (
Fig. 3B). Consequently, the index case in Hubei likely contracted SARS-CoV-2 on or around 4 November 2019 (95% upper HPD: 15 October; 99% upper HPD: 7 October) (
Fig. 3C).
This time frame for the Hubei index case is robust (fig. S7). Epidemic simulations with faster or slower transmission rates and more or less densely connected contact networks produce similar date estimates. Furthermore, using the root tMRCA from the sampled posterior trees, rather than adjusting for the shifting coalescence between 24 December 2019 and 1 January 2020, produces a median date of 3 November (95% upper HPD: 14 October; 99% upper HPD: 4 October). Incorporating a relaxed molecular clock or adjusting the coalescent prior assumption (Skyline versus Skygrid) also had minimal effect.
If we enforce that there must be at least one ascertained case in our simulations before 17 November 2019, the median date of the Hubei index case is pushed back about a week to 28 October (95% upper HPD: 11 October; 99% upper HPD: 5 October) (fig. S7). However, the distinction between ascertained and unascertained in the original SAPHIRE model was meant to reflect the probability of missed diagnoses in January 2020 of the Wuhan epidemic and does not account for the investigations that resulted in retrospective diagnoses in November and December 2019.
If we discount the reported evidence of retrospective COVID-19 diagnoses throughout the end of November and instead take 1 December as representing the first confirmed case of COVID-19, then the median time between index case infection and stable coalescence after rejection sampling is 23 days (95% range: 1 to 47 days) (
Fig. 3D). Under this scenario, the index case in Hubei would have contracted SARS-CoV-2 on or around 17 November 2019 (95% upper HPD: 24 October; 99% upper HPD: 13 October) (
Fig. 3E). Similar dates are inferred with a relaxed molecular clock and conditioning on an ascertained infection by 1 December (fig. S7).
It is reasonable to postulate that the variant of SARS-CoV-2 that first emerged was less fit than the variant that spread through China and that evolutionary adaptation was critical to its establishment in humans (
12). Therefore, we simulated two-phase epidemics in which the index case was infected with a less-fit variant (i.e., half as transmissible) that went extinct, but not before giving rise to a mutant strain matching the transmission dynamics estimated in Wuhan (figs. S8 and S9). If we condition on an ascertained or unascertained COVID-19 case (due to either the original or adapted variant) by 17 November, the original variant transmits for a median of 5 days (95% HPD: 0 to 36 days) before the adapted strain emerges (
Fig. 3F); this adapted strain then circulates for a median of 28 days (95% HPD: 0 to 52 days) before reaching a stable coalescence (
Fig. 3G). In this scenario, the index case in Hubei would have likely contracted an unobserved variant of SARS-CoV-2 on 5 November 2019 (95% upper HPD: 11 October; 99% upper HPD: 25 September) (
Fig. 3H). If we again discount the reported evidence of COVID-19 in November and take 1 December as the first confirmed case, then the virus would spend less time in humans across both phases of the early epidemic (
Fig. 3, I and J), and the index case in Hubei would have acquired SARS-CoV-2 on 18 November 2019 (95% upper HPD: 20 October; 99% upper HPD: 5 October) (
Fig. 3K). As in the primary analysis, the inferred date of the index case in the two-phase epidemic was robust to varying model assumptions (fig. S10), including the amount by which viral fitness differed between the two phases (supplementary text and fig. S18).
These one-phase and two-phase epidemic simulations both suggest that the tMRCA is not representative of the emergence of SARS-CoV-2. In the primary analysis, the index case remains infected at the tMRCA (as in
Fig. 2, upper left panel) in <1.1% of simulated epidemics (table S3). In the two-phase epidemics, the index case is still infected at the tMRCA in 2.2% of simulations, likely because of an increased variance in the time between index case and tMRCA. The initial, less-fit variant persisted until the tMRCA in 17% of simulated epidemics and until 1 January 2020 in 3.7% of simulated epidemics. However, when this less-fit variant persisted, it was represented by only a single infected individual at the tMRCA; this low frequency suggests that even if this less-fit variant did exist, it could have easily been missed in early genome sequencing efforts. In the two-phase epidemic, the first ascertained (or unascertained) case was due to the less-fit variant in around two-thirds of simulated epidemics.
By anchoring our epidemic simulations to specific tMRCA estimates, we can reconstruct a plausible range for the number of SARS-CoV-2 infections before the discovery of the virus (
Fig. 4A). The median number of individuals infected with SARS-CoV-2 in our primary analysis is less than one until 4 November. The median number of infected individuals is four (95% HPD: 1 to 13) on 17 November and reaches nine (95% HPD: 2 to 26) on 1 December. These values are generally robust to model specifications, molecular clock method, and date of first COVID-19 case (table S4 and fig. S11). The two-phase epidemics tend to exhibit similar growth patterns (
Fig. 4B, fig. S12, and table S5). Notably, we do not see any evidence for an increase in hospitalizations until mid- to late December, even when we increase the virulence of the less-fit variant in the two-phase epidemics or increase the probability of hospitalization before the stable coalescence (supplementary text and figs. S13 and S14).
Empirical observation throughout the SARS-CoV-2 pandemic has shown the outsized role of superspreading events in the propagation of SARS-CoV-2 (
24–
27), wherein the average infected person does not transmit the virus. Our results suggest that the same dynamics likely influenced the initial establishment of SARS-CoV-2 in humans, as only 29.7% of simulated epidemics from the primary analysis went on to establish self-sustaining epidemics. The remaining 70.3% of epidemics went extinct (
Fig. 4C). Simulated epidemics that went extinct typically produced only 1 infection (95% range: 1 to 9) and never more than 44 infections total or 14 infections at any given time (table S2). The median failed epidemic went extinct by day 8. As the contact network became more or less densely connected, the number of epidemics that went extinct was similar: 68.3 and 69.4%, respectively (table S2). However, the percentage of extinct epidemics increased as the transmission rate decreased (80.5%) and decreased as the transmission rate increased (53.6%). In the two-phase epidemic, this original less-fit variant went extinct by day 9 (95% HPD: 2 to 52) and produced a median of one infection (95% HPD: 1 to 13) (
Fig. 4D and table S6).
The overdispersed nature of SARS-CoV-2 transmission patterns favors its persistence, as epidemics simulated over random contact networks (with the same mean number of contacts) that are not characterized by superspreading events (
27) tended to go extinct more frequently, 83.7% of the time. Furthermore, the large and highly connected contact networks characterizing urban areas seem critical to the establishment of SARS-CoV-2. When we simulated epidemics in which the number of connections was reduced by 50 or 75% (without rescaling per-contact transmissibility) to reflect emergence in a rural community, the epidemics went extinct 94.5 or 99.6% of the time, respectively.
Our results highlight the unpredictable dynamics that characterized the earliest days of the COVID-19 pandemic. The successful establishment of SARS-CoV-2 postzoonosis was far from certain, as more than two-thirds of simulated epidemics quickly went extinct. It is highly probable that SARS-CoV-2 was circulating in Hubei province at low levels in November 2019 and possibly as early as October 2019, but not earlier. Nonetheless, the inferred prevalence of this virus was too low to permit its discovery and characterization for weeks or months. By the time that COVID-19 was first identified, the virus had firmly established itself in Wuhan. This delay highlights the difficulty in surveillance for novel zoonotic pathogens with high transmissibility and moderate mortality rates.
The high extinction rates we inferred suggest that spillover of SARS-CoV-2–like viruses may be frequent, even if pandemics are rare (
28). Furthermore, the same dynamics that characterized the establishment of SARS-CoV-2 in Hubei province may have played out all over the world, as the virus was repeatedly introduced but only occasionally took hold (
29,
30). The reports of cases in December 2019 and January 2020 in France and California that did not establish sustained transmission fit this pattern (
31–
33). However, our results suggest that polymerase chain reaction evidence of SARS-CoV-2 in wastewater outside of China before November 2019 is unlikely to be valid (
34), and the suggestion of international spread in mid-November or early December 2019 should be viewed with skepticism (
35–
37), given that our results suggest that fewer than 20 people were infected with SARS-CoV-2 at this time (table S4 and fig. S11). Our results also refute claims (
38) of large numbers of patients requiring hospitalization because of COVID-19 in Hubei province before December 2019 (figs. S13 and S14). Nevertheless, SARS-CoV-2 may be detectable in archived wastewater samples or other biomaterials from Hubei province from early to mid-November 2019, should they exist, and incorporating these types of data in our model could further refine our timing estimates. Moreover, wastewater detection may present the best chance of early detection of future pandemics during the early phase of spread for which we estimate very low numbers of infections (
39).
Even though all of the earliest documented cases of COVID-19 were found in Hubei province, we cannot discount the possibility that the index case initially acquired the virus elsewhere. Nonetheless, our dating inference is insensitive to geography. Furthermore, our results suggest that if the virus first emerged in a rural community, it would have needed to migrate to an urban setting to avoid extinction. The lack of reports of COVID-19 elsewhere in China in November and early December suggests that Hubei province is the location where human-to-human transmission chains were first established.
The circumstances surrounding the emergence of SARS-CoV-2 in Hubei province remain shrouded. Although SARS-CoV-2 is repeatedly adapting to spread among humans (
40,
41), our findings do not reveal whether the virus that first emerged was less fit than the virus that spread throughout China. Nevertheless, the inferred timing of the index case is generally similar in both of these scenarios because less-fit viruses in our simulations that went extinct tended to do so very quickly. It is yet unknown whether the virus emerged directly from its animal reservoir [presumably horseshoe bats (
42,
43)] or first circulated in and possibly adapted to an intermediate host. Our estimates for the timing of the Hubei index case further distance this individual from the outbreak at the Huanan Seafood Wholesale Market. Finding the animal reservoir or hypothetical intermediate host will help to further narrow down the date, location, and circumstances of the original SARS-CoV-2 infection in humans. However, even in the absence of that information, coalescent-based approaches permit us to look back beyond the tMRCA and toward the earliest days of the COVID-19 pandemic. Although there was a pre-tMRCA fuse to the COVID-19 pandemic, it was almost certainly very short. This brief period of time suggests that future pandemics with similar characteristics to those of the COVID-19 pandemic permit only a narrow window for preemptive intervention.
Acknowledgments
We gratefully acknowledge the authors from the originating laboratories and the submitting laboratories, who generated and shared via GISAID the viral genomic sequence data on which this research is based. A complete list acknowledging the authors who submitted the data analyzed in this study can be found in data S1. We thank J. Havens, M. Sanderson, and B. Wertheim for fruitful conversations and insight.
Funding: J.O.W. acknowledges funding from the National Institutes of Health (AI135992 and AI136056). J.P. acknowledges funding from the National Institutes of Health (T15LM011271) and the Google Cloud COVID-19 Research Credits Program. M.W. was supported by the David and Lucile Packard Foundation as well as the University of Arizona College of Science. N.M. acknowledges funding from the National Science Foundation (2028040) and the Google Cloud COVID-19 Research Credits Program.
Author contributions: Conceptualization: J.O.W.; Methodology: J.P., N.M., M.W., K.S., and J.O.W.; Software: J.P., N.M., and J.O.W.; Validation: J.P. and J.O.W; Formal analysis: J.P. and J.O.W.; Investigation: J.P. and J.O.W.; Resources: J.O.W.; Data curation: J.O.W.; Writing – original draft: J.P. and J.O.W.; Writing – review & editing: J.P., M.W., K.S., J.O.W., and N.M.; Visualization: J.P. and J.O.W.; Supervision: J.O.W.; Project administration: J.O.W.; and Funding acquisition: J.O.W. and M.W.
Competing interests: J.O.W. has received funding from Gilead Sciences, LLC (completed), and the CDC (ongoing) via grants and contracts to his institution unrelated to this research.
Data and materials availability: All data used in this analysis are free to access via GISAID. BEAST XML input, FAVITES JSON, and results files are available in Dryad (
44). This work is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. To view a copy of this license, visit
https://creativecommons.org/licenses/by/4.0. This license does not apply to figures/photos/artwork or other content included in the article that is credited to a third party; obtain authorization from the rights holder before using such material.
RE: Origin of Sars-CoV-2
Prof. Zhang used an Illumina NovaSeq machine to decode this virus on the 5th Jan 2020. The CEO of Illumina, Mr. Francis D'Souza in his April '21TIME magazine article by Alice Park, states he was contacted in November 2019 and it took 2 months to decode tis virus. Furthermore he states that the authorities were already trying to decode the virus before contacting him because of a pneumonia in Wuhan. Arithmetically, we have proof from the very person who provide the decoding machine that this was this virus was causing infections by late October/ early November.
Appropriateness of forward and retrospective coalescent analysis to infer age of the SARS-CoV-2 index case
We appreciate the attention Huang and Jiang have paid to our Report; however they misunderstand our methodological approach both in concept and in practice. To recapitulate the dynamics of the Wuhan epidemic, we used the SAPHIRE compartmental model. To incorporate individuality in these simulations, we added an underlying transmission network. A preferential attachment network is necessary to appropriately capture the transmission network of viruses like SARS-CoV-2 (1). Hence, transmission rates calibrated under a compartmentalized model (like SAPHIRE) must be adjusted to accommodate the underlying network, so as to not overestimate the transmission rate per individual. As we point out in our Report (2), our timing inference was robust to the specific number of mean contacts in this transmission network.
We agree with Huang and Jiang that the network growth rate saturates in our model, particularly beyond 10 thousand simulated infections. If we had wanted to simulate the entire Wuhan epidemic, then 10 million nodes would have been appropriate. But our timing estimates were specifically focused on simulations up to the time of stable coalescence, at which there was a maximum of 35 infected nodes—well before network saturation would become a concern.
The use of rejection sampling served as an extension of the Bayesian phylodynamic inference from BEAST, which does not permit the inference of infections prior to the time of most recent common ancestor (tMRCA) (3). Huang and Jiang expressed reservations about the step where we modeled P(X,Y|Z,Dc) as proportional to C(Y)P(X,Y|Z), with C(Y) being a consistency function. Here, P(X,Y|Z) is the degree of plausibility of a pair of hypothesized dates (index case date X and first ascertained/unascertained case date Y), given that only the tMRCA date Z is known, while P(X,Y|Z,Dc) is the plausibility of the same pair of dates when we have also learned the date of a reported case (Dc). It is inconsistent with epidemiological data for the first appearance of the virus to occur after the reported case date, so learning this date rules out any such scenarios: it would be nonsensical to use a model in which learning Dc does not reduce the probability of X,Y to zero. Beyond that effect, our assumption is minimal in that it does not allow learning Dc to further impact our information state: the relative plausibilities of hypotheses that are not ruled out by the new information remain unchanged. The claim by Huang and Jiang that "X, Y were both predicted from the model and not dependent on Dc" suggests a misunderstanding of the framework within which we are working: the equation applies to any X and Y regardless of whether predicted by a model, and describing the dependence of a hypothesized (X,Y) on Dc is exactly the purpose of the equation. Our assessment of the plausibility of a hypothesis certainly does depend on whether or not we possess information rendering the hypothesis impossible. In this context, the SAPHIRE-based simulations serve as a prior distribution for the time from index case to tMRCA; the rejection sampling permits us to infer an epidemiologically- and phylogenetically-informed posterior distribution. Importantly, failure to reject epidemiological unrealistic combinations would bias the inferred date for the index case forward in time.
Lastly, Huang and Jiang misunderstand why we included the phylodynamic analysis in the first place. We maintain that the inferred tMRCA is a robust approach for connecting the stable coalescence in SAPHIRE-based simulations to real-world calendar dates. The estimated number of SARS-CoV-2 cases in late-2019/early-2020 is imprecise, and the outsized role of superspreading events makes simulating the epidemic trajectory leading to these early case counts unreliable. Hence, we developed an approach that combines epidemic simulations with phylodynamic inference to overcome these issues and provide greater precision in timing the index case.
References
1. M. J. Keeling, K. T. D. Eames, Networks and epidemic models. J. R. Soc. Interface. 2, 295–307 (2005).
2. J. Pekar, M. Worobey, N. Moshiri, K. Scheffler, J. O. Wertheim, Timing the SARS-CoV-2 index case in Hubei province. Science. 372, 412–417 (2021).
3. M. A. Suchard, P. Lemey, G. Baele, D. L. Ayres, A. J. Drummond, A. Rambaut, Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 4, vey016 (2018).
RE: Timing the SARS-CoV-2 index case in Hubei province
RE: Timing the SARS-CoV-2 index case in Hubei province
Sisi Huang 1,2 & Chao Jiang 1*
Affiliations:
1. Life Sciences Institute, Zhejiang University, Hangzhou 310058, China.
2. School of Mathematical Sciences, Zhejiang University, Hangzhou 310027, China.
*Correspondence: [email protected] (C.J.)
The statistical approaches implemented in Pekar et al. (1) for timing the index case requires further optimizations to be useful. The authors used the parameters of the SAPHIRE model in Hao et al. (2) but took a different modeling approach (GEMF model) to combine the transmission network with phylogenetic inference (3-5). For clarity, we will refer to the transmission model in Pekar et al. (1) as SAPHIRE-GEMF, which should be mathematically equivalent to the SAPHIRE model. We first discuss the impact of authors' modifications on the SAPHIRE model. The simulated curves of SAPHIRE-GEMF are highly variable and failed to resemble reality. In addition, both SAPHIRE and SAPHIRE-GEMF models can time the index case without involving tMRCA. Finally, the rejection sampling approach described in Pekar et al. is inappropriate and we show theoretical proof and simulation results.
First, the modifications made in the SAPHIRE-GEMF model may deviate the simulation from reality. SAPHIRE-GEMF and SAPHIRE models are mathematically equivalent when simulating the same dynamic equations. However, if the interaction network is unevenly connected, such as preferential attachments (6), the average number of contacts may not reflect the actual spreading process. Therefore, the inaccuracy of the simulated contact network structure could be why the authors needed to greatly reduce the transmission rate b from 1.31 to 0.385 to fit the median doubling time. However, median doubling time alone is not a robust indicator for the model as the simulated curves do not reflect reality. The indicator of the median doubling time assumes the growth of the infected population N as a function of t is exponential. The 50 (out of 200) successfully established simulations using 100,000 nodes (also inaccurate as the actual number of nodes is nearly 10 million) showed drastically slowed growth at an early stage, with a much lower peak of daily ascertained cases for the first 100 days compared to the reality and the curves simulated with 1,000,000 nodes.
Second, both SAPHIRE-GEMF and SAPHIRE models can estimate the index case without involving tMRCA. Both models had been trained to be consistent with the epidemic growth in Wuhan during the outbreak period January 1-23, 2020 (1, 2). Assuming the SAPHIRE-GEMF and SAPHIRE models were suitable during early pandemic, as Pekar et al. assumed (1), we could fit the simulated curves to real data by minimizing the mean square error to find the index case. For the SAPHIRE-GEMF model, we chose 1,000,000 nodes to prevent premature saturation of the contact network when using only 100,000 nodes. We ran 5000 epidemic simulations and 1513 simulations successfully established epidemics. We then implemented the rejection sampling approach to remove "epidemiologically meaningless" results (1). As a result, if we assume November 17 represents the first documented case of COVID-19 (ascertained or unascertained in the SAPHIRE model), the index case in Wuhan likely contracted SARS-CoV-2 around November 9, 2019 (95% upper HPD on October 25; 99% upper HPD on October 17). Similar simulations can be implemented using the SAPHIRE model, and the median date for the index case is around November 12. Therefore, both models can time the index case without involving phylogenetic analysis.
Thirdly, the rejection sampling approach is inappropriate. The authors tried to sample the index case X and the date Y of the index case from P(X,Y|Z,Dc)=1[Y consistent with Dc]P(X,Y|Z), where Z denotes the common ancestor (tMRCA) and Dc denotes the first reported case (1). However, this equation may not hold because X, Y were both predicted from the model and not dependent on Dc. Therefore, the rejection sampling subjectively selects preferred samples of the index case prior to a given date Dc. As indicated in Pekar et al., many of the estimations did not precede the earliest dates of reported COVID-19 cases (1) and were artificially neglected by rejection sampling. In other words, sampling dates before November 17 resulted in heavily biased samples from the left tail of the distribution. When including all simulations, the median dates of infection for the index case are December 4 and December 1 for SAPHIRE-GEMF and SAPHIRE, respectively. Therefore, the date of the index case could be easily pushed around depending on the preset rules of rejection sampling. As 87.2% (1320/1513) failures of the SAPHIRE-GEMF simulation results were not epidemiologically meaningful, the biases introduced by rejection sampling should be further evaluated. In the authors' analysis, combining the phylogenetic analysis with epidemic simulations would require sampling twice (sample tMRCAs and days from index case infection to tMRCA), which further exaggerated the variation.
In conclusion, while we readily appreciate the authors' efforts in timing the index case, the modeling approaches require further investigations. We show that the SAPHIRE-GEMF and the SAPHIRE models are mathematically equivalent but parameters were modified without rigorous justifications. The quick saturation of nodes (N = 100,000) led to a deviation from reality. In addition, we can time the index case using both models without involving tMRCA. Finally, rejection sampling lacks a theoretical foundation and could lead to subjective selection bias of the sampling results. Therefore, instead of rejection sampling, the SAPHIRE-GEMF model still needs significant revisions given the extremely high failure rate (96.14% overall) of simulations. The author's explanation that this high failure rate reflects how real infectious diseases originate is only one of many possibilities.
References and Notes
1. J. Pekar, M. Worobey, N. Moshiri, K. Scheffler, J. O. J. S. Wertheim, Timing the SARS-CoV-2 index case in Hubei province. 372, 412-417 (2021).
2. X. Hao et al., Reconstruction of the full transmission dynamics of COVID-19 in Wuhan. 584, 420-424 (2020).
3. F. D. Sahneh, C. Scoglio, P. J. I. A. T. o. N. Van Mieghem, Generalized epidemic mean-field model for spreading processes over multilayer complex networks. 21, 1609-1620 (2013).
4. F. D. Sahneh, A. Vajdi, H. Shakeri, F. Fan, C. J. J. o. c. s. Scoglio, GEMFsim: A stochastic simulator for the generalized epidemic modeling framework. 22, 36-44 (2017).
5. N. Moshiri, M. Ragonnet-Cronin, J. O. Wertheim, S. J. B. Mirarab, FAVITES: simultaneous simulation of transmission networks, phylogenetic trees and sequences. 35, 1852-1861 (2019).
6. A.-L. Barabási, R. J. s. Albert, Emergence of scaling in random networks. 286, 509-512 (1999).
Acknowledgments: We thank Dr. Susan P. Holmes for the valuable discussions.
Data and materials availability: Codes are available on GitHub at https://github.com/Hc1023/time-science-test.
RE: Timing the SARS-CoV-2 index case in Hubei province"
Dear Editor
Tracking the spread of SARS-CoV-2 in real time is crucial to monitoring novel zoonotic pathogens and inform public health efforts as recently highlighted by Pekar et al., ("Timing the SARS-CoV-2 index case in Hubei province") and other recent reports (1, 2, 3)
Advances in DNA sequencing technology have ushered in a new era of pan-genomics and genomic surveillance, in which traditional molecular diagnostics and genotyping methods are being enhanced and even replaced by genomics-based methods to aid epidemiologic investigations of communicable diseases (3, 4). Considering the necessity to monitor the emergence of new viral pathogens in nearly-real-time, bioinformatics tools and the combination of genomic and epidemiological data can give essential information for understanding the past and the future of a new pandemic virus . Thus, coalescent framework based on a combination of retrospective phylodynamic inference and prospective epidemiological simulations, helped by the unprecedented number of cases and more than 800K SARS-2-CoV genomes generated, allowed to determine the cryptic circulation of SARS-CoV-2 with its probable introduction between mid-October and mid-November 2019 (1). In this line, worthy of mention is the early observation made by Benvenuto et al (2) who, by employing a molecular clock model, could date back to the fourth week of November 2019 (95%HPD: September 28, 2019; December 21, 2019) the origin of this novel pathogen. This dating, which is consistent with the estimation made by Pekar et al (1) remains quite noticeable in consideration of the very limited set of data, particularly regarding the paucity of SARS-CoV-2 complete genomes available at the beginning of the epidemic, that of course may justify its particular wide 95% HPD. Pekar et al (1) hypothesis of a frequent spill-over of beta-coronaviruses from animal reservoirs, lineages extinctions and success are indeed attractive and invites to a major effort for a much-needed improvement of international surveillance of zoonotic diseases. Together those findings supports the notion that continued genomic surveillance strategies are needed to assist in the monitoring and understanding of new emerging viral pathogens to ensure substantial improvements in the control of future epidemics with spill-over potential to humans (6,7).
1. Pekar J, Worobey M, Moshiri N, Scheffler K, Wertheim JO. Timing the SARS-CoV-2 Index Case in Hubei Province. bioRxiv [Preprint]. 2020 Nov 24:2020.11.20.392126. doi: 10.1101/2020.11.20.392126. Update in: Science. 2021 Mar 18: PMID: 33269353; PMCID: PMC7709179
2. Benvenuto D, Giovanetti M, Salemi M, Prosperi M, De Flora C, Junior Alcantara LC, Angeletti S, Ciccozzi M. The global spread of 2019-nCoV: a molecular evolutionary analysis. Pathog Glob Health. 2020 Mar;114(2):64-67. doi: 10.1080/20477724.2020.1725339. Epub 2020 Feb 12. PMID: 32048560; PMCID: PMC7099638.
3. A. Rambaut, "Phylodynamic Analysis | 176 genomes | 6 Mar 2020." Virological (2020); https://virological.org/t/phylodynamic-analysis-176-genomes-6-mar-2020/356.
4. Tegally H, Wilkinson E, Giovanetti M, Iranzadeh A, Fonseca V, Giandhari J, Doolabh D, Pillay S, San EJ, Msomi N, Mlisana K, von Gottberg A, Walaza S, Allam M, Ismail A, Mohale T, Glass AJ, Engelbrecht S, Van Zyl G, Preiser W, Petruccione F, Sigal A, Hardie D, Marais G, Hsiao M, Korsman S, Davies MA, Tyers L, Mudau I, York D, Maslo C, Goedhals D, Abrahams S, Laguda-Akingba O, Alisoltani-Dehkordi A, Godzik A, Wibmer CK, Sewell BT, Lourenço J, Alcantara LCJ, Kosakovsky Pond SL, Weaver S, Martin D, Lessells RJ, Bhiman JN, Williamson C, de Oliveira T. Emergence of a SARS-CoV-2 variant of concern with mutations in spike glycoprotein. Nature. 2021 Mar 9. doi: 10.1038/s41586-021-03402-9. Epub ahead of print. PMID: 33690265.
5. Gardy JL, Loman NJ. Towards a genomics-informed, real-time, global pathogen surveillance system. Nat Rev Genet. 2018 Jan;19(1):9-20. doi: 10.1038/nrg.2017.88. Epub 2017 Nov 13. PMID: 29129921; PMCID: PMC7097748.
6. T. Bedford, A. L. Greninger, P. Roychoudhury, L. M. Starita, M. Famulare, M.-L. Huang, A. Nalla, G. Pepper, A. Reinhardt, H. Xie, L. Shrestha, T. N. Nguyen, A. Adler, E. Brandstetter, S. Cho, D. Giroux, P. D. Han, K. Fay, C. D. Frazar, M. Ilcisin, K. Lacombe, J. Lee, A. Kiavand, M. Richardson, T. R. Sibley, M. Truong, C. R. Wolf, D. A. Nickerson, M. J. Rieder, J. A. Englund, J. Hadfield, E. B. Hodcroft, J. Huddleston, L. H. Moncla, N. F. Müller, R. A. Neher, X. Deng, W. Gu, S. Federman, C. Chiu, J. S. Duchin, R. Gautom, G. Melly, B. Hiatt, P. Dykema, S. Lindquist, K. Queen, Y. Tao, A. Uehara, S. Tong, D. MacCannell, G. L. Armstrong, G. S. Baird, H. Y. Chu, J. Shendure, K. R. Jerome; Seattle Flu Study Investigators, Cryptic transmission of SARS-CoV-2 in Washington state. Science 370, 571–575 (2020). doi:10.1126/science.abc0523pmid:32913002
7. M. A. Suchard, P. Lemey, G. Baele, D. L. Ayres, A. J. Drummond, A. Rambaut, Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 4, vey016 (2018). doi:10.1093/ve/vey016pmid:29942656
Conflict of Interest:
None declared