Time-varying transmission heterogeneity of SARS and COVID-19 in Hong Kong

Transmission heterogeneity is a notable feature of the severe acute respiratory syndrome (SARS) and coronavirus disease 2019 (COVID-19) epidemics, though previous efforts to estimate how heterogeneity changes over time are limited. Using contact tracing data, we compared the epidemiology of SARS and COVID-19 infection in Hong Kong in 2003 and 2020-21 and estimated time-varying transmission heterogeneity (k t ) by tting negative binomial models to offspring distributions generated across variable observation windows. k t uctuated over time for both COVID-19 and SARS on a continuous scale though SARS exhibited signicantly greater (p < 0.001) heterogeneity compared to COVID-19 overall and in-time. For COVID-19, k t declined over time and was signicantly associated with increasingly stringent non-pharmaceutical interventions though similar evidence for SARS was inconclusive. Underdetection of sporadic COVID-19 cases led to a moderate overestimation of k t , indicating COVID-19 heterogeneity of could be greater than observed. Time-varying or real-time estimates of transmission heterogeneity could become a critical indicator for epidemic intelligence in the future.


Abstract
Transmission heterogeneity is a notable feature of the severe acute respiratory syndrome (SARS) and coronavirus disease 2019 (COVID-19) epidemics, though previous efforts to estimate how heterogeneity changes over time are limited. Using contact tracing data, we compared the epidemiology of SARS and COVID-19 infection in Hong Kong in 2003 and 2020-21 and estimated time-varying transmission heterogeneity (k t ) by tting negative binomial models to offspring distributions generated across variable observation windows. k t uctuated over time for both COVID-19 and SARS on a continuous scale though SARS exhibited signi cantly greater (p < 0.001) heterogeneity compared to COVID-19 overall and in-time.
For COVID-19, k t declined over time and was signi cantly associated with increasingly stringent nonpharmaceutical interventions though similar evidence for SARS was inconclusive. Underdetection of sporadic COVID-19 cases led to a moderate overestimation of k t , indicating COVID-19 heterogeneity of could be greater than observed. Time-varying or real-time estimates of transmission heterogeneity could become a critical indicator for epidemic intelligence in the future.

Main
Superspreading appears a distinct feature of coronavirus transmission during past and recent human outbreaks of novel coronaviruses including the ongoing coronavirus disease 2019 (COVID- 19) pandemic (1)(2)(3)(4)(5)(6)(7)(8). The effective reproductive number (Re), and time-varying effective reproduction number (Re t ), are essential epidemiological measures used to quantify rates of pathogen transmission and predict the progression of epidemics, however, these measures often omit heterogeneity in individual level transmission and superspreading effects. Typically, disease transmission is modelled according to a Poisson distribution where individual variance in transmission is equal to Re. However, this distribution is inappropriate when modelling datasets that feature superspreading because the variance greatly exceeds the expected mean. The negative binomial distribution is a convenient model to measure transmission heterogeneity when parameterized by the mean (µ), equivalent to Re and dispersion parameter, k, where the variance is equal to µ + µ 2 /k (9). For low values of k (0 < k < 1), the corresponding distribution is long tailed but concentrated around zero. Such overdispersion indicates the likelihood of superspreading events (9, 10) which can vastly alter epidemic dynamics and optimal intervention strategies (11). For example, highly overdispersed transmissions are often associated with large clusters of infection, thus interventions that aim to limit superspreading can have a greater marginal impact on Re (12)(13)(14). Unlike Re (and Re t ) however, k is traditionally interpreted as a xed characteristic immutable to nonpharmaceutical interventions (NPIs), and few studies have investigated potential temporal changes in k (15,16) and the relative effect of NPIs on transmission heterogeneity over time. Using detailed contact tracing data, we presented a novel approach to estimating transmission heterogeneity over time on a continuous scale, i.e., time-varying overdispersion (k t ) and other related statistics. We applied our approach to investigate and compare the epidemiology and temporal transmission heterogeneity with two coronavirus datasets from Hong Kong notable for superspreading: SARS in 2003, and the ongoing COVID-19 pandemic (1, 17).

Comparative transmission dynamics of SARS and COVID-19
in Hong Kong Figure 1 shows the epidemic curve of SARS in 2003 and COVID-19 in Hong Kong between Jan 23, 2020, and Apr 5, 2021. During the three epidemic waves of COVID-19 transmission, the majority of the local cases (n = 6,112/9,216, 66.34%) were not traceable to superspreading events (de ned as greater than six secondary cases per infector case (1)). This was in contrast to SARS where the majority of the cases matched to contact tracing data (n = 848/1,295, 65.48%) were linked to superspreading events either as cases exposed directly to source of infection (primary SSE exposure cases) or subsequent generations traced to past events (other SSE associated cases). The mean overall estimate of Re t for SARS was higher at 1.5 (range = 0.45, 5.8) compared to 1.1 (range = 0.37, 2.3) for COVID-19, though the range was notably wider for SARS ( Fig. 2A). Consistent with the known kinetics of both pathogens, the mean serial interval for all resolved transmission pairs was 5.4 days for COVID-19 and 12.1 days for SARS. Temporally, the SARS epidemic was marked by a substantial rise in Re t during the early phase of the outbreak in late February 2003 that preceded the peak in late March and coincided with a consistent decline and eventual plateau of Re t throughout the nal phase of the outbreak ( Fig. 2A). This compared to COVID-19 where Re t was marginally distributed around 1 but never peaked as high or declined as low as SARS during periods of local transmission ( Fig. 2A).
The overall dispersion parameter for SARS in Hong Kong was k = 0.04 (Table 1). For COVID-19, we estimated a higher overall k of 0.20. From this we calculated the proportion of cases responsible for 80% of onwards transmissions (Prop80) as 4.2% for SARS in 2003, where 87.4% of cases did not transmit to anyone. In comparison, 14.2% of COVID-19 cases were responsible for 80% of onward infections in Hong Kong while fewer cases, 69.9%, did not transmit to anyone (Table 1, Supplementary Table 1). Table 1 Overall estimates* of dispersion parameter (k) and proportion of cases responsible for infecting 80% (Prop80) and 0% (Prop0, i.e., no one) of cases for SARS and COVID-19 in Hong Kong.   Table S2). Similar declines were inconclusive during the rst COVID-19 wave seemingly due to the wider con dence intervals; a result of smaller sample sizes (Fig. 1). For SARS, we found that k t and Prop80 t appeared to signi cantly increase over time ( Fig. 2B Figure   S1). For each epidemic (excluding COVID-19 wave one), given Re t and k t , we showed that the probability of epidemic extinction increased over time for both SARS and COVID-19 as expected given that each epidemic was locally eliminated (Fig. 2D, Supplementary Table S2,).
Temporal variations in k t and Prop80 t for COVID-19 appeared to partially correlate temporally with large superspreading events and clusters, particularly when case numbers were low and / or towards the declining phase of each epidemic (Fig. 1, Fig. 2B-C, Supplementary Figure S2). For example, during the third COVID-19 wave, a sharp decline in k t and Prop80 t coincided with a large gym cluster (n = 155) detected in mid-March 2021. Conversely, Prop80 t and k t were mostly insensitive to change following detection of the largest COVID-19 case cluster in Hong Kong during the initial epidemic rise in wave three.
This cluster however was primarily associated with attendance across at least 21 dance and ballroom venues (n = 732), such that transmissions were unlikely to be caused by a single source case, while 48.9% (n = 355/732) of cases were secondary or tertiary cases etc. without direct exposure to the venues such as family members. Other apparent spikes in k t and Prop80 t could be seen throughout each COVID-19 wave that partially coincided with large superspreading events ( Fig. 1, Fig. 2B-C, Supplementary Figure   S1). For SARS, observed k t and Prop80 t were in contrast relatively monotonic following the initial introduction of a case from mainland China which triggered the rst clusters linked to the Metropole Hotel  Table S4). This however was not the case for SARS where the relationship was inconclusive (Fig. 3, Supplementary Table  S4). For COVID-19, increasing NPI stringency appeared to be associated with signi cant declines in k t and Prop80 t , which was not seen for SARS as with Re t (Fig. 3B, Supplementary Table S5). We did however nd a signi cant association between increasing NPI stringency and decreasing Re t and increased probability of extinction as expected for both COVID-19 and SARS (Supplementary Figure S3, Supplementary Table S5).

Imperfect observation bias and underdetection
To investigate the potential effect of imperfect observation of COVID-19 cases on our estimates of k t, we considered hypothetical worst-case scenarios involving under and over observation of chain-terminating cases in the data. To assess under-observation, we incorporated additional unobserved hypothetical cases as chain terminating singletons into the observed data. Here, the k t estimated from the base case tended towards overestimation as expected ( Fig. 4 left panel) i.e., the original k t was higher than the k t estimated from the sensitivity analysis on the hypothetical scenarios (Supplementary Figure S4). When testing a worst-case scenario where a constant 50% of singletons cases remained undetected, the k t appeared to be overestimated with a median of 0.1 (90% CI: 0.05, 0.31) while the estimate was reduced to 0.03 (0.01, 0.08) under a more plausible scenario where the underdetection was 10% (Supplementary   Table S6).
Next, to assess the upper boundary of potential over-observation, we excluded all observed singletons cases from the data while again testing a constant 10% -50% underdetection rate as before. Here, k t shifted upwards, as expected, compared to the estimate from the base case demonstrating the upper range for potential overestimation of k t (Fig. 4 right panel). Excluding the hypothetical worst-case (50% underdetection), the marginal distribution of k t for COVID-19 remained signi cantly higher than SARS (Supplementary Figure S5, Supplementary Table S7). Reassuringly, the absolute magnitude of bias was small, even in these worst-case-scenario analyses.

Discussion
This study presented for the rst time a measurement of time-varying transmission heterogeneity (here k t and Prop80 t ) and compared two signi cant beta-coronavirus epidemics in Hong Kong: SARS in 2003 and the ongoing COVID-19 pandemic ( Fig. 1 and Fig. 2). Few studies have demonstrated changing transmission heterogeneity between temporally distinct COVID-19 epidemic periods (15,16), however we presented the variation in transmission heterogeneity on a continuous scale (daily) during periods of sustained local transmission with minimal to no international introductions which would otherwise confound local inference. We found that measures of transmission heterogeneity uctuated over time and were partially correlated with large superspreading events, though measures for SARS exhibited signi cantly greater (p < 0.001) heterogeneity and were less variable over time compared to COVID-19 ( Fig. 2, Supplementary Figure S2). Furthermore, k t and Prop80 t declined over the course of each epidemic wave. Counter to expectations, declines in k t and Prop80 t were signi cantly associated with increasingly stringent NPIs that largely targeted settings for potential superspreading events such as venue closures ( Fig. 3D-E). A similar effect was observed in a recent study which saw overall estimates of k declined after strict shelter-in-place measures were introduced (16). We found that k t and Prop80 t were also signi cantly associated with Re t (Fig. 3A-B). Previous studies suggested this might be because population-wide NPIs, such as physical distancing and universal masking (that also reduced Re t ), effectively limited the number of random contacts in the community and thus increased the overall proportion of chain terminating cases, i.e., cases without onward infection, which is a hallmark feature of transmission heterogeneity (13,18). Interpretation of k t therefore is not simply proportional to the occurrence of superspreading events, though we did partially observe this, but individual transmissions overall. Thus, increasing measures of heterogeneity over time, including in response to NPIs, can indicate improving epidemic control as transmissions, including superspreading events, are prevented before they occur (13,19). This is also the rst study to our knowledge to estimate overall overdispersion for SARS during the 2003 epidemic in Hong Kong, which was notable for the occurrence of large superspreading events. Our most conservative estimate of transmission heterogeneity (k = 0.05, 95% CI: 0.04, 0.07) however is lower than the only previous estimate for SARS (k = 0.16, 90% CI: 0.11, 0.64) using data from Singapore (9), demonstrating markedly greater transmission heterogeneity for SARS in Hong Kong (Table 1,   Supplementary Table 1). Such difference between our results and others is perhaps indicative of a unique pathogen-population dynamic at the time. One major strength of our study therefore is the shared population setting, meaning comparisons between transmission heterogeneity for SARS and COVID-19 better accounts for potential confounders if comparing outbreaks between different populations not considering population changes between the epidemics. We have also re ned previous estimates of COVID-19 transmission heterogeneity in Hong Kong by greatly increasing the number of transmission pairs available for analysis (n = 4,697 pairs minimum compared to n = 169 pairs for earlier estimates), producing additional evidence that favors greater levels of transmission heterogeneity for COVID-19 than previously found (Table 1, 12-17% vs 15-24% previously) (1). Our results are closer, but still higher than global estimates that showed < 10% of COVID-19 cases were responsible for 80% of transmissions (2,3,7,21,22). This consistent difference between overall estimates of heterogeneity between Hong Kong and elsewhere again could indicate a distinct pathogen-population dynamic in Hong Kong, though this might also be due to potential underdetection which we showed could lead to a moderate overestimation of k t (Fig. 4).
Our study has some limitations. For COVID-19, no variants of concern (VOCs) sustained local transmission in Hong Kong during the study period, therefore our conclusions were based on originally dominant lineages and therefore might not be generalizable to VOCs, including those with altered transmissibility e.g., Delta (B.1.617.2) (23, 24), Omicron (B.1.1.529) (25, 26) and other novel VOCs that may emerge in the future. This study also relied on contact tracing data largely generated through both active and passive surveillance in order to construct transmission pairs from which offspring distributions and measures of transmission heterogeneity were generated, and thus assumed all cases were observed. However as shown in our study, and consistent with past studies, imperfect observation and underdetection could lead to moderate overestimation of k (including k t ), thus underestimating the true degree of heterogeneity in transmission especially if singleton cases were over or underrepresented in the data. (Sporadic, chain-terminating cases are 0's in the offspring distribution; holding all other things equal, a higher number of observed 0's is a hallmark feature of greater overdispersion.) (18). If contact tracing fails to identify all valid epidemiological links, singleton cases are overrepresented. Alternatively, they could be underrepresented if singleton cases are harder to detect than large clusters which are more likely to observed (18). Observation biases are particularly problematic when the probability of observing epidemiological links between cases is low. However, in Hong Kong, COVID-19 observation probabilities are assumed to be relatively high given that each epidemic wave was effectively controlled throughout 2020 and 2021 via strict testing and trace, isolate and quarantine (20,27), which makes us con dent that relatively few infections are missing from our data, and the vast majority of observed COVID-19 cases were traced (Fig. 1). Regardless, under various worst case hypothetical scenarios, we showed that if 50% of COVID-19 cases in Hong Kong were undetected, the marginal density of k t for COVID-19 was no longer signi cantly higher compared to SARS. However, when the same worst-case rate of underdetection was similarly applied to the overall estimates of k, heterogeneity remained signi cantly greater for SARS than for COVID-19 (k = 0.04, 95% CI: 0.03, 0.06 vs. k = 0.1, 95% CI: 0.08, 0.12 respectively).
Finally for SARS, it is likely that most cases were observed during the epidemic due to the relative severity, and thus the estimates of k t for SARS were less likely to be overestimated due to underdetection compared with COVID-19. Instead, because 26.2% con rmed SARS cases could not be matched to available contact tracing data and were therefore excluded from all offspring-related analyses, which could result in bias in either direction.

Characterization of clusters and transmission pairs
For COVID-19, we characterized case clusters and transmissions pairs con rmed between Jan 23, 2020, and Apr 5, 2021, using data provided by the Centre for Health Protection in Hong Kong, which was available for all 9,195 local cases detected during the study period. This period covered the epidemic waves of local transmission of COVID-19 in Hong Kong. Clusters linked to known imported index cases were excluded to limit inference of NPI effects to local measures only i.e., excluding travel related measures. Characterization, including the generation of transmission pairs, was performed as per our previous study with complete details of the construction process also available in the supplementary methods (1). Overall, we constructed two datasets each for COVID-19 and SARS sensitivity analyses. First for COVID-19, we identi ed 4,697 resolved transmission pairs (comprising 6,531 cases) where the speci c infector and infectee could be identi ed (the primary dataset). Second, we additionally included cases whose infector was unknown but who could be traced to a common exposure source, i.e., as infectees linked directly to an exposure setting (hereafter denoted as unresolved transmission pairs & the sensitivity dataset), increasing the number of pairs to 6,196 pairs (comprising 7,640 cases). For SARS in 2003, contact tracing was performed and data provided by the Department of Health, and could be matched to demographic data for 1,293 (73.6%) of the 1,755 cases con rmed within Hong Kong. The 462 SARS cases with that could not be matched were excluded from all subsequent analyses. Generation of transmission pairs was the same as for COVID-19, though we additionally reviewed publicly available reports to identify the con rmed index cases of large clusters where contact tracing data was ambiguous. As before, we constructed two datasets for sensitivity resulting in 472 transmission pairs (comprising 552 cases) for SARS where the speci c infector and infectee was known, and 1,061 pairs (comprising 1,160 cases) when including infectees linked to common exposure settings.

Estimating transmission heterogeneity of SARS and COVID-19
Heterogeneity in transmission can be modelled according to a negative binomial distribution and corresponding dispersion parameter k t to the individual numbers of secondary infections (offspring distribution), where the mean parameter, μ, is equivalent to the effective reproductive number (Re). Following the approach described by Lloyd Smith et al., (2) we generated offspring distributions for each of the primary and sensitivity datasets for SARS and COVID-19 by counting the number of secondary cases per primary case or exposure setting for all resolved and unresolved pairs. Cases without traceable onward transmission, including cases without a known source (unlinked singleton cases), were included in each distribution as chain terminating cases i.e., zero secondary cases, excluding those SARS cases without contact tracing data as before. Negative binomial distributions were tted by maximum likelihood (ML) estimation to the primary and sensitivity offspring distributions independently using tdistrplus v1.1-3 xing Re = 1 for control interpretation of k between groups and datasets. Uncertainty around each point estimate of k was generated by parametric bootstrap resampling (n=1,000). Final overall estimates were a composite of the primary (upper bound) and sensitivity analyses (lower bound) calculated as the median and 95% quantiles of the combined bootstrapped distributions (See Supplementary Table S1 for independent results). To provide an intuitive comparison of overall transmission heterogeneity between SARS and COVID-19 on a linear scale, we calculated the proportion of individuals responsible for 80% of onward transmission (Prop80, by analogy to the 20/80 rule) as per the formula below given Re = 1 (3): Estimating time-varying transmission heterogeneity To estimate time-varying k, i.e., k t , negative binomial distributions were t by ML estimation to offspring distributions subset within dynamic sliding windows based on a minimum sample size and xed window period forecast from each of the known infector's onset date, or con rmation date in the case of asymptomatic infection. When the infector was unknown, the infectees onset date or con rmation date when asymptomatic was used. For the dynamic windows we tested periods with a minimum of n = (10, 20, …, 50) paired cases. However, to avoid extending the dynamic windows beyond reasonable time periods when samples were low, windows were limited to one or two weeks regardless of the minimum sample size. Similarly, for the xed windows, we subset offspring distributions into both one and twoweek sliding windows. To simplify parameter estimation in time and control against biases in contact tracing data when case numbers were high, we estimated the median Re t using EpiNow2 (4) for each of the dynamics and xed sliding window periods which was then xed as the input as parameter μ when estimating k t . Uncertainty around each estimate was again generated by parametric bootstrap resampling. A nal composite estimate of k t was calculated by day as the median and 90% quantiles of the combined bootstrapped distribution of each of the variable window periods from both the primary and sensitivity datasets. We generated 90% CIs here instead of 95% for k t because more extreme values of k were most di cult to estimate accurately with small sample sizes (5) as a result of sub-setting the observed data. This also allowed easier visualization of k through time and matched with the results as presented for SARS in Lloyd-Smith et al (2). With both Re t and k t known in time, we calculated the instantaneous probability of epidemic extinction given the probability generating function of a negative binomial offspring distribution (2) denoted by : Where the probability that an epidemic has gone extinct by the n th generation is Complete details of the estimation procedure are available in the supplementary methods.

Temporal associations of changing transmission heterogeneity and NPIs
We conducted simple linear regression analyses to identify associations between changing transmission heterogeneity over time and temporally with Re t and NPI stringency given the nal estimates of k t and Prop80 t . We performed 1000 independent regression analyses each randomly sampling within the combined bootstrapped distribution of k t and Prop80 t n estimates equal to the number of days in each epidemic period with valid data. Signi cance was determined as the mean p-value of the 1000 estimates and adjusted by applying Bonferroni's correction for multiple comparisons. Because k t exhibited a skewed distribution (0 < k → ∞), values were transformed using the natural logarithm before linear regressions were conducted. The relative stringency of NPIs over time was calculated for COVID-19 and SARS based on the Oxford Government response tracker indices (level 1: <40; level 2 : 40-50; level 3: 50-60; level 4: 60-70; level 5: >70) (6). We compared the distributions of k t between COVID-19 and SARS by Welch Two Sample t-test, performing 1000 independent comparisons sampled from the marginal bootstrapped distribution. Because k t exhibits a skewed distribution (0 < k → ∞), estimates were transformed using the natural logarithm before t-tests were conducted Measuring bias due to incomplete observation and underdetection of COVID-19.
In order to investigate the magnitude and boundaries of possible bias due to the presence of an overburdened contact tracing system and asymptomatic infection, we considered worst-case scenarios that could cause us to underestimate or overestimate k t . First, we could underestimate k t if chainterminating cases were overrepresented in our data (e.g., if epidemiological links between cases were frequently missed in contact tracing). To assess this worst-case underestimation bound, we repeated the main analysis excluding unlinked singleton cases observed in the community ("without singletons" analysis). Second, we could overestimate k t if chain-terminating cases were systematically unobserved (e.g., if contact tracing made it easier to identify linked cases than singletons). We tested multiple worstcase hypothetical scenarios to investigate the overestimation boundary of k t given extremes of potential observation bias by additionally incorporating unobserved singletons as chain terminating cases given assumed rates of underdetection (between 10 and 50% unobserved cases) as a proportion of the observed cases within each sliding window period prior to ML estimation of k. Together this conceptualizes the reasonable boundaries of k t shown in Fig. 4. As before, we compared each the resulting hypothetical distributions of log-transformed k t to the observed distribution for SARS by Welch Two Sample t-test, performing 1000 independent comparisons sampled from the marginal bootstrapped distribution.

CODE AVAILABILITY
The code used for analysis is publicly available at https://github.com/dcadam/.... Figure 1 Epidemic curve of 9,237 locally acquired and con rmed COVID-19 cases, and 1,755 con rmed SARS cases in Hong Kong colored by source relative to known superspreading events. Arrows indicate the date of detection of the largest linked superspreading events and their primary setting including total number of linked cases overall including traceable onward infections.