Comparison of Estimated Effectiveness of Case-Based and Population-Based Interventions on COVID-19 Containment in Taiwan

Key Points Question What are the explanations for the initial success of COVID-19 control in Taiwan, a country that has one of the lowest per capita incidence and mortality rates in the world? Findings In this comparative effectiveness research study that used detailed epidemiologic and contact tracing data, neither case-based interventions (including contact tracing and quarantine) or population-based interventions (including social distancing and facial masking) alone were estimated to have been sufficient to contain COVID-19. The combination of case-based and population-based interventions was needed. Meaning The combination of case-based interventions with population-based interventions with wide adherence may explain the success of COVID-19 control in Taiwan.


eMethods 1. Estimation of the interval parameters
The incubation period, the onset-to-isolation interval, and the serial interval (time interval between symptom onsets in an infector-infectee pair) were estimated from the dates (or date intervals if the exact date is unknown) of exposure, symptom onset, isolation of the confirmed SARS-CoV-2 cases in Taiwan. We estimated the distributions of these interval parameters using a Bayesian framework 1 that deals with the situation where the exact date is uncertain for both ends of the interval (i.e., doubly interval-censored). This method uses a hierarchical model to estimate the exact two ends of the time interval for each individual as well as the overall distribution of the interval among the population. In the case of incubation period estimation, the exact infection and onset times were unobserved and uniform priors bounded by individual exposure and onset windows were assumed. For each individual the exact time points of infection and symptom onset are modelled as follows.
is the exact time of infection, which was not observed. � , , , � represents the time window within which transmission could have occurred, and was obtained based on patients' travel/exposure history in case investigation reports. Likewise, is the unobserved exact time of symptom onset and � , , , � is the time window within which symptom onset could occurred. is the individual incubation period, and was assumed to follow a gamma distribution with the shape ( ) and the scale ( ). We applied flat exponential priors to and as follows.
(1/1,000) ∼ (1/1,000) To estimate the serial interval, we simply replaced the exposure-to-onset quantities with onset-to-onset quantities between the primary and secondary cases. To estimate the onset-to-isolation interval, we simply replaced the exposure-to-onset quantities with onset-to-isolation quantities Specifically, we added an offset of 1 days for both the serial interval and the onset-to-isolation interval (i.e. add 1 days on each interval ) because there are negative values generated by the observed windows. We included these observations for that negative serial intervals and negative onset-to-isolation intervals are plausible in practice. Particularly, negative onset-to-isolation intervals occurred due to contact tracing and active surveillance measures, and negative serial intervals occurred when the incubation period of the index case happened to be much longer than the incubation period of the secondary case.
For each estimation, we ran 4 Markov chain Monte Carlo (MCMC) chains using the Non-U-Turn Sampler in Stan. 2 Each chain contains 1,000 warm-up iterations and 500 samples, rendering to a total of 2,000 samples. The credible intervals are obtained from the 2.5 th and 97.5 th percentiles in the posterior predictive simulations of the gamma distributions.

Overview
We adopted the structure of the branching process model developed by Hellewell, et al., 3 which in essence consists of two components: 1) the branching process, and 2) the intervention model. The branching process simulates the disease transmission dynamics in the early phase of an outbreak. The intervention model considers the effects of various case-based interventions for COVID-19 control. Specifically, case detection, contact tracing, and quarantine of close contacts are considered. The parameters and their values for this dynamical transmission model are listed in Table 1.

The branching process
The branching process simulates growing transmission trees starting from some given initial cases. This process is implemented by the following steps.
i. Initiate with initial active cases. ii. For each active symptomatic case, determine the onset time and the testing time by drawing an individual incubation period and an onset-to-isolation interval from their probability distributions (estimated from case series data). The testing time for active asymptomatic cases are set to infinite. iii. For each active case (both symptomatic and asymptomatic), draw the number of newly infected cases from the distribution of secondary cases. iv. For each active case (both symptomatic and asymptomatic), apply the intervention model (see eMethods 2.3 for more details) and determine the period of quarantine and isolation. v. For each secondary case, determine the infection time by drawing an individual generation interval (infection times in an infector-infectee pair). vi. For each transmission pair, determine whether the transmission is realized, or prevented by comparing the infection time of the secondary case and the period of quarantine/ isolation of the index case. vii. Deactivate the index cases, and activate the realized secondary cases. viii. Repeat step ii~vii, until there is no active case, or when the maximal number of generations is reached.
The distributions of the incubation period and onset-to-isolation interval were estimated from case series data with the method described in eMethods 1. We assumed a negative binomial distribution for the secondary case distribution, which is governed by the reproduction number ( ) and the dispersion parameter ( ). We assume that = 0 = 2.5 in the counterfactual scenario where no interventions and behavioral changes were in effect. In the real-world setting should equal , denoting for the reproduction number under the effect of population-based interventions. and were estimated by fitting the dynamical model to the observed cluster sizes (see eMethods 3.3). The generation interval distribution was assumed to be a skewed normal distribution centered at each index case's onset time to avoid the inconsistent length of incubation periods and generation intervals. This parameterization also makes the proportion of pre-symptomatic transmission ( ) an explicit parameter in our model. The standard deviation of the generation interval ( ), and the proportion of pre-symptomatic transmission ( ) were estimated by fitting the dynamical model to the observed serial intervals (see eMethods 3.2). Regarding the possibility of asymptomatic infection, we assumed a fixed probability of being asymptomatic ( = 0.4), and fixed relative transmissibility ( = 0.5) for all infections.

The intervention model
The intervention model is composed of a set of rules which determine whether and when active cases and their contacts are quarantined or isolated and hence unable to transmit disease. Specifically, case detection, contact tracing, and quarantine of close contacts were implemented as described in the following.
i. Case detection: Each active, untraced, and symptomatic case was tested with probability , and was immediately isolated if tested positive. The secondary cases generated during the incubation period (presymptomatic period) and the onset-to-isolation interval cannot be prevented. ii. Contact tracing: Each secondary case (except for initial introductions) was ascertained as a close contact (successfully traced) of the detected index case with probability . If a secondary case was ascertained and showed symptoms on the time of contact tracing (the onset time of the detected index case), the secondary case was immediately isolated, unless the case had already been detected and isolated. Case detection plus contact tracing was able to prevent the transmission during the onset-to-isolation interval, but not during the incubation period. iii. Quarantine of close contacts: Each active and traced case (regardless of the presence of symptoms) was immediately quarantined at the time of being traced. If the case develops symptoms during the quarantine period, he/she will be immediately isolated. Only the combination of detection, tracing, and quarantine was able to prevent pre-symptomatic transmissions and transmissions from asymptomatic cases.
Note that asymptomatic cases were never detected, traced or isolated, but could be quarantined. We also assumed perfect compliance with isolation and quarantine order, and all transmissions were prevented during that period. eFigure 2 gives examples that illustrate the effects of these interventions on the prevention of disease transmission. eFigure 3 provides the one-way sensitivity analysis of all the parameters in this model.

Overview and the Sequential Monte Carlo (SMC) algorithm
The observed serial intervals and size distribution of stuttered transmission chains were utilized to calibrate the dynamical model and to estimated related parameters. We note that the reproduction number ( ) and the dispersion parameter ( ) only affect the number of infections, not the temporal relationship between successive generations, e.g. the distribution of the serial interval. In contrast, the probability of pre-symptomatic transmission ( ) and the standard deviation of the generation interval ( ) affect the number of secondary infections because they interact with case-based interventions, such that pre-symptomatic transmissions are hardly prevented by case-based interventions. Therefore, we first calibrated the model and estimated and using the observed serial intervals. Then the resulting model was used for further calibration to the observed cluster sizes, and for the estimation of and .
Since the likelihood of the dynamical model is intractable, we used a sequential Monte Carlo algorithm to obtain posterior distributions of parameters of interest. This algorithm was used to fit branching process-based dynamical models in a similar context. 4,5 . The algorithm started from a population of 1,000 parameter sets drawn from the prior distributions. Data were simulated with the branching process model parameterized by these parameter sets, and the distance between the simulated and empirical data was measured by Kolmogorov-Smirnov (KS) statistics.
In each round of iteration, the parameter set was resampled, perturbed, and passed on until the criteria of convergence were met. We ran 4 chains of SMC algorithms to generate a total of 4,000 posterior samples for inference.
The steps of an SMC algorithm are as follows: i. Initiation: generate a population of 1,000 parameter sets by Latin Hypercube sampling from prior distributions. ii. Simulation: simulate data points with each parameter set and the dynamical model. iii. Evaluation: calculate the KS statistics. iv. Evolution: resample a new population of 1,000 parameter sets from the current population weighted by 1/KS 2 . v. Mutation: perturb each new parameter set by up to 10%. vi. Repeat step ii~v until the median KS statistics of the population is less than 0.05 and is within 10% of each of the previous two rounds.

Model calibration with observed serial intervals
The probability of pre-symptomatic transmission ( ) and the standard deviation of the generation interval ( ) were estimated by fitting the model to the observed serial interval distribution. Wide uniform priors were assigned to both parameters. The posterior mean estimates are 0.55 (95%CrI 0.41-0.68) and 2.70 (95%CrI 1.88-3.76), for and , respectively. eTable 1 lists the values of priors, posterior estimates and other fixed parameters used in this calibrating stage. eFigure 4 presents the posterior distributions, the convergence plots, and the posterior predictive of the serial interval.

Model calibration with observed cluster sizes
The reproduction number under population-based interventions only ( ), and the dispersion parameter ( ) were estimated in this stage, by fitting the model to observed cluster sizes. Wide uniform priors were assigned to the parameters. The posterior mean estimates are 1.30 (95%CrI 1.03-1.58) and 29.21 (95%CrI 6.28, 50.00), for and , respectively. eTable 2 lists the values of priors, posterior estimates and other fixed parameters used in this calibrating stage. eFigure 5 presents the posterior distributions, the convergence plots, and the posterior predictive of the cluster sizes. The limited amount of observations for cluster sizes make it hard to infer the dispersion parameter, as the posterior distribution of did not converge well (eFigure 5B). However, the posterior distribution of did converge to consistent estimates (eFigure 5A) which had accounted for the uncertain in simultaneously.

Estimation of the time-varying reproduction numbers
The time-varying reproductive number ( ) for SARS-CoV-2 and influenza were estimated using Wallinga and Teunis method, 6 also known as the "case reproduction number". 7 This method attributes the transmission events and assigns the value of to the cohort of primary cases at time . Since represents the transmissibility of primary cases, it explains the future incidence, and can reflect the subsequent impacts of events after specific time points. 8 . Practically, this method estimates the transmission probabilities between every possible transmission pair, according to their observed serial intervals. The probability that case was infected by case ( ) is given as is calculated by normalizing the likelihood of case infecting case by the sum of the likelihood from all possible infector cases ≠ . � − � is the transmission likelihood quantifying how well the observed serial interval (the onset time difference between case and , − ) fits the serial interval distribution of ascertained transmission pairs. The effective reproduction number of case is by definition (the expected number of secondary infections) the sum of all the transmission probabilities where case is the infector.

= �
We then summarized the 's into the time-varying reproduction numbers ( ) by calculating the 7-day moving averages according to their onset time. The confidence intervals were calculated by the 2.5th and 97.5th percentiles in the of 100 simulated transmission trees from the 's matrix, as in Cori,et al. 5 For SARS-CoV-2, we directly used the daily incidence based on the symptom onset date to estimate . For influenzae, the weekly incidence from two different data sources were used, including the notified influenzae cases with severe complications in the National Notifiable Disease Surveillance System, and the influenzaee-like illness (ILI) consulting rate in the out-patient and emergency departments. The ILI consultation rate was further multiplied by the positive rate of influenzae according to the laboratory surveillance data. Cubic spline smoothing was used to disaggregate the weekly-basis data into daily-basis incidence. 9 To avoid the problem of right truncation, we used the influenzae data through the end of March to estimate the Rt until Feb 21st (one month after the first reported case of Covid-19 in Taiwan). Another key input to estimation is the distribution of the serial interval. For SARS-CoV-2, the serial interval was estimated using the ascertained transmission pairs in our case series data. For seasonal influenzae, we assumed the mean and standard deviation of the serial interval to be 3.6 and 1.6 according to a previous study 10 . In addition to the Wallinga and Teunis method, we also estimated the instantaneous (real-time) reproduction number Rt by Cori et al. for the influenzae analysis, and the results were shown in eFigure 6.

eFigure 1. Conceptual framework for the transmission model.
The transmission model is a stochastic branching process model that simulates the epidemic based on given input parameters (in a forward fashion). On the other hand, the transmission model could be used to estimate the parameter values through fitting/calibrating the model output to the observed data (in a backward fashion). A two-step approach was used in this study. First, the stochastic model was calibrated to the observed serial interval distribution to estimate the proportion of pre-symptomatic transmission and the standard deviation of generation interval (eMethods 3.2). Second, the calibrated model from the first step was used to estimate and project the potential effectiveness of case-based interventions in the absence of population-based interventions (the arrow from R 0 to R c ). In this case, R 0 is the input parameter and R c is the output from the model. The model was also used to estimate potential effectiveness of population-based interventions and joint interventions (the arrow from R p to R pc ). In this case, R p is the (unknown) input parameter and R pc is the output from the model that depended on R p . By fitting the cluster size distribution from the model to the cluster size distribution observed in contact tracing, R p and R pc were jointly estimated. R 0 , basic reproduction number (assumed); R c , effective reproduction number with case-based interventions (estimated); R p : effective reproduction number with population-based interventions (estimated); R pc , effective reproduction number with both case-based and population-based interventions.

eFigure 2. Examples of the effects of case detection, contact tracing, and quarantine.
Case A and A* demonstrate the effect of mere detection, which can only prevent the transmission once the active cases are tested and isolated. That is, the active cases can transmit the disease during their incubation and delay of case detection. Case B and B* demonstrate the effect of detection plus tracing (without quarantine), where case B was successfully traced and onset within a buffer period. Therefore, case B was immediately isolated when the source was detected. Detection plus tracing can prevent transmission during the delay of case detection. Case C and C* demonstrate the combined effect of detection, contact tracing, and quarantine. Only in this scenario that transmission during the incubation period can be prevented. We assume that there is no delay between testing and isolation, the buffer of contact tracing to be one day, and the same effect of quarantine and isolation. Besides, asymptomatic cases are never detected or traced but could be quarantined and have lower transmissibility.

eFigure 3. Tornado diagrams from the one-way sensitivity analysis on the effects of casebased interventions.
The effective reproduction number under case-based interventions ( ). The blue bars represent the change in the measured outcome when the corresponding parameter value decreased; the red bars represent the change when the parameter value increased. The tuning ranges of the parameters are shown next to ends of the bars. eFigure 6. The incidence and instantaneous reproduction number (Rt) of influenzae in Taiwan, 2018-2020.
(A) Estimates from the notified number of severe influenzae with complications. (B) Estimates from the overall influenzae cases derived from influenzae-like illness consultation rate and the positive rate from laboratory testing for influenzae (eMethods). The gray bars represent the number of weekly incidence cases, and the red lines represent the Rt with 95% confidence intervals in the shaded area. Thirty days before and after (yellow and blue background) January 21, the date of the first SARS-CoV-2 infection confirmed in Taiwan