Clusters of patients infected with Middle East respiratory syndrome (MERS) coronavirus continue to occur in countries throughout the Middle East, where the virus is thought to be endemic in camels (Kayali and Peiris, 2015) . While rare, countries elsewhere in the world experience importations from infected individuals traveling from the endemic region (Carias et al., 2016) . Most identified importations of MERS from travelers have not resulted in documented transmissions in the destination country ; however, the recent large cluster of 186 infected patients stemming from a single introduction in the Republic of Korea (ROK) (Korea Centers for Disease Control and Prevention, 2015) demonstrated that explosive outbreaks are possible.
The ROK outbreak, combined with a non-negligible likelihood of further exportations of MERS from Middle Eastern countries (Carias et al., 2016) , is cause for continued concern for importation of MERS to other countries. For public health officials requiring quantitative assessment of the risk posed by incoming infected travelers, it is important to have a nuanced understanding of the full spectrum of possible outcomes, especially when they are highly variable (Fisman et al., 2014) ; modeling studies can play an important role in this regard.
Recent studies Kucharski and Althaus, 2015; Chowell et al., 2015) have quantified the variability implied by different data sets of MERS cluster sizes resulting from importation of cases. These analyses found that the data are potentially consistent with high transmission variability associated with the occurrence of superspreading events, similar to what was observed during severe acute respiratory syndrome (SARS) outbreaks in 2003 (Lloyd-Smith et al., 2005) . These studies quantified transmission probabilities using a negative binomial offspring distribution within a branching process outbreak model, assuming that every infected individual transmits with an average of R 0 transmissions and dispersion parameter k, where k < 1 implies high over-dispersion (Lloyd-Smith et al., 2005) .
In this paper, we extend the results of the above work to allow the reproductive number R to vary across subsequent generations of transmissions during an outbreak. The ROK outbreak consisted of a large number of transmissions from the initial traveler and from a few patients in the next transmission generation. was an extremely rapid decrease in transmissions such that the entire outbreak was extinguished after three total generations of transmission following the introduction (Korea Centers for Disease Control and Prevention, 2015) . This type of differential transmissibility before vs. after implementation of control measures has also been observed during localized outbreaks of SARS (Lloyd-Smith et al., 2005; Wallinga and Teunis, 2004) and Ebola (Toth et al., 2015; Shuaib et al., 2014) .
A simple way to model a post-control change in average transmissibility is to use one parameter for the reproduction number in early generations (R 0 ) and another for later generations (R c , or post-control reproductive number), as assumed in several previous modeling studies of observed outbreaks and public health response for different diseases (Lloyd-Smith et al., 2005; Wallinga and Teunis, 2004; Toth et al., 2015; Chowell et al., 2004) . We hypothesized that a model allowing this type of switch would produce a substantially better fit to the data from outbreak clusters caused by MERS importations. Given results from our previous work assessing Ebola importation risk (Toth et al., 2015) , we also hypothesized that this model might produce substantially different results for the risk of a very large outbreak compared to a model assuming a single reproductive number across all transmission generations.
We developed a data set of cluster sizes from MERS importations to countries entirely outside of the Arabian Peninsula (Table 1) ; we excluded data from Jordan, the Kingdom of Saudi Arabia, Kuwait, Oman, Qatar, the United Arab Emirates, and Yemen, countries where it was not always clear whether the initial or subsequent cases within clusters acquired infection from exposure to MERS cases or animals (camels). The data were extracted from World Health Organization reports (World Health Organization, 2015) as well as published accounts of individual clusters (Yavarian et al., 2015; Puzelli et al., 2013; Abroug et al., 2014 ; The Health Protection Agency U. K. Novel Coronavirus Investigation team, 2013). Our data set consists of 31 importation events, of which 23 resulted in no confirmed or suspected transmissions (clusters of size 1) and the other 8 resulted in clusters of size 2-186. Following Nishiura et al. (2015) , we also recorded the total number of generations of transmission that occurred after the introduction.
For each generation of transmission, we assumed a negative binomial offspring distribution with parameter set Â i = (R i , k i ), where i is the generation of transmission (i = 0 from the initial traveler). This assumption results in the following equations. Each row represents a unique individual infected traveler to the indicated country. a Cluster size includes the initial infected traveler and any subsequent infected persons epidemiologically linked to that traveler; a cluster of size 1 indicates no known transmission from the traveler in the destination country. b Transmission generations are the maximum number of transmission links from an infected person in the cluster back to the initial traveler.
First, the probability that x independent cases in generation i produce a total of y cases in generation i + 1 is
Next, given n independent introductions (generation 0), the joint probability of a cluster of total size j consisting of exactly g generations of transmission, under parameter set Â = (Â 0 , Â 1 , Â 2 , Â 3 ), is
We used the above equations to evaluate ten different models. In Model 0, we assumed constant parameter values across all generations of transmission, i.e.,Â 0 = Â 1 = Â 2 = Â 3 = (R, k). In Models 1a, 1b, and 1c, we assumed the initial patient transmitted with reproductive number R 0 and dispersion parameter k 0 , and all subsequent patients transmitted with a post-control reproductive number R c and dispersion parameter k c i.e., Â 0 = (R 0 , k 0 ) ; Â 1 = Â 2 = Â 3 = (Rc, k c) . Because we found that allowing k c to range freely in the optimization scheme resulted in wide uncertainty (due to few multi-generation clusters in the data), we chose to test three differ-ent assumptions for this parameter. In Model 1a, we assumed that k c = k 0 ; in Model 1b, we assumed k c = 1, a special case in which the negative binomial distribution reduces to the geometric distribution; and in Model 1c, we assumed infinite k c , another special case in which the negative binomial distribution reduces to the Poisson distribution.
In Models 2a, 2b, and 2c, we assumed that the reproductive number and dispersion parameter switched from R 0 to R c and k 0 to k c after two generations of transmission, i.e., Â 0 = Â 1 = (R 0 , k 0 ) ; Â 2 = Â 3 = (Rc, k c) , and made the same three assumptions regarding k c as described above. In summary,
For each of these three parameterizations Â, we quantified the likelihood of observing the 31 clusters of size j m extinguished after g m generations using the formula L = 31 m=1 q Â (1, j m , g m) .
We compared the maximum likelihood fits under the three models using the Akaike information criterion (AIC), which evaluates model parsimony in determining statistical support for the hypothesized difference in transmission across outbreak generations (Blumberg et al., 2014) .
We also developed a model extension to test the robustness of our results against the possibility that there were additional MERS exportations outside the Arabian Peninsula causing clusters that were not detected. If undetected clusters exist, the data set in Table 1 might be biased toward larger cluster sizes, as smaller clusters presumably would be more likely to go undetected.
To quantify the implications of undetected clusters, we made the following assumptions for this part of the analysis. Let N be the number of undetected clusters, and u be the probability that an individual infected patient outside the Arabian Peninsula goes undetected. We assumed that if any one patient in a cluster was detected with MERS, then the entire cluster was detected, due to the intensive contact tracing that would be initiated after the first detection. Under those assumptions, the probability that a cluster of size j would go undetected is u j . We also assumed that transmission among patients in an undetected cluster was governed by the R 0 , k 0 parameters only, under every model, because the presumed mechanism for shifting to R c , k c (implementation of transmission control measures) would only be relevant if detection occurred.
The new likelihood L N for a given test value of N undetected clusters is comprised of the product of the joint probabilities that each of the 31 clusters was of the given size and number of generations and was detected, times the probability that N clusters were unob-served; this latter factor includes the probabilities for undetected outbreaks of any size.
We estimated the infinite sum using a partial sum that had converged to six decimal places. The likelihood was maximized for N = 31 and N = 93, representing scenarios where 50% and 75% of importation clusters were undetected, respectively, over the parameters u, R 0 , k 0 , and R c if applicable, for Models 0, 1b, and 2b.
Model 0 produced an MLE of R = 0.87 (95% CI: 0.46-1.90) and k = 0.035 (0.016-0.069). Model 1a, assuming a change in the reproductive number after the first generation of transmission, produced an MLE of R 0 = 5.2 (1.7-29.9), R c = 0.19 (0.05-0.53), k 0 = k c = 0.068 (0.031-0.14). Model 1b produced R 0 = 5.5 (1.7-35.4), R c = 0.14 (0.05-0.29), and k 0 = 0.061 (0.025-0.13). Model 1c produced R 0 = 5.5 (1.7-35.7), R c = 0.14 (0.05-0.27), and k 0 = 0.061 (0.025-0.13).
Model 2a, assuming the change occurred after the second generation of transmission, resulted in estimates of R 0 = 2.0 (1.0-6.7), R c = 0.064 (0.007-0.27), and k 0 = k c = 0.078 (0.034-0.16). Model 2b produced R 0 = 2.2 (1.0-6.8), R c = 0.060 (0.008-0.18), and k 0 = 0.076 (0.032-0.16). Model 2c produced R 0 = 2.2 (1.0-6.8), R c = 0.060 (0.008-0.17), and k 0 = 0.076 (0.032-0.16).
Each version of Models 1 and 2 produced an MLE with substantially higher likelihood and lower AIC than Model 0, the two-parameter (R, k) model previously implemented Kucharski and Althaus, 2015; Chowell et al., 2015) . Of these, models 2b and 2c produced the lowest AIC value (Table 2) ; we chose Model 2b to represent an optimal model under this criterion.
We compared the risk assessment implications of the optimal model against those of other models. The optimal model produces a higher probability of smaller outbreaks across one or two generations of transmission, but a much lower probability of very large outbreaks or of outbreaks exceeding several transmission generations (Table 3) .
The results under the assumption of undetected clusters (Table 4) show that Model 2b is still optimal according to AIC, although the change in AIC compared to Model 0 becomes smaller as the number of assumed undetected clusters increases. Also, as the number of assumed undetected clusters increases, the optimal model's estimate of "worst-case" outbreak sizes at the 0.1% or 0.01% probability level move closer to those of the simpler Model 0 ( Fig. 1 panels A, C, E) . However, the optimal model still produces much lower estimates of the probability of outbreaks lasting several generations across all assumptions for undetected clusters (Fig. 1  panels B, D, F) .
We have considered a simple method to assess the statistical support for differential transmission in earlier versus later generations after a new introduction of MERS, based only on outbreak data for the sizes of transmission clusters and total number of transmission generations that produced them. This method demonstrated strong statistical support for assuming a higher reproductive number in earlier generations after a MERS introduction in a non-endemic area.
Projections from the optimal model have important implications for assessing the risk posed by new introductions of MERS. Compared to previous assessments Kucharski a For Model 0, the reproductive number R is the average number of transmissions from each individual regardless of the transmission generation; for Models 1a, 1b, and 1c, the initial reproductive number R0 and dispersion parameter k0, apply to the initial traveler only (generation 0), and the post-control reproductive number Rc and dispersion parameter kc apply to any infected persons in generations ≥1; for Models 2a, 2b, and 2c, R0 and k0 apply for both generations 0 and 1, and Rc and kc apply for generations ≥2.
b Parameters were optimized according to the shown maximal log likelihood. c AIC = Akaike information criterion, used to determine the optimal model (Model 2b represents an optimal model, with lowest AIC value).
Risk assessment implications of each model.
Probability of >(10, 100, 500, 1000) total transmissions Probability of >(1, 2, 3, 4, 5) generations of transmission Model 0 None (3.9%, 1.0%, 0.3%, 0.1%) (11%, 4.5%, 2.6%, 1.7%, 1.2%) Model 1a 1 (12%, 1.5%, 0.007%, 0.00002%) (26%, 3.6%, 0.6%, 0.12%, 0.02%) Model 1b 1 (11%, 1.6%, 0.011%, 0.00006%) (24%, 3.3%, 0.5%, 0.07%, 0.01%) Model 1c 1 (11%, 1.6%, 0.011%, 0.00006%) (24%, 3.3%, 0.4%, 0.06%, 0.009%) Model 2a 2 (11%, 1.6%, 0.008%, 0.00002%) (23%, 14%, 0.8%, 0.05%, 0.003%) Model 2b 2 (11%, 2.0%, 0.018%, 0.00011%) (23%, 14%, 0.9%, 0.05%, 0.003%) Model 2c 2 (11%, 2.0%, 0.018%, 0.00011%) (23%, 14%, 0.8%, 0.05%, 0.003%)
Probabilities of exceeding selected numbers of total transmissions/generations of transmission after a single importation of Middle East respiratory syndrome, under three different models. Model 3 was the optimal model given the data in Table 1 , according to criterion summarized in Table 2 .
Sensitivity analysis -results of fitting models to the cluster data given that portion of importation clusters were undetected. Chowell et al., 2015) that were similar to those from our Model 0, our optimal Model 2b suggests a higher probability of moderately sized outbreaks (e.g., on the order of 10 total transmissions) across one or two total generations of transmission, but a much lower probability of outbreaks significantly larger than the one in ROK or of outbreaks of any size lasting several generations. These conclusions are robust to assuming that only 50% of MERS importations outside the Arabian Peninsula have been detected. If the non-detection rate is much higher than 50%, then our optimal model would produce closer estimates to previous models for the probability of very large outbreaks, but the conclusion that outbreaks are less likely to last several generations than previous predictions is robust to high rates of non-detection.
The results from all the models we fit to the data suggest very high transmission variability from the index patient (and perhaps also from subsequent patients, depending on the model), as the MLE for the parameter k 0 was less than 0.1 for each model, which indicates even higher over-dispersion than what was estimated for SARS (Lloyd-Smith et al., 2005) . The MLE value of k 0 was even lower in the analyses assuming there were undetected clusters, as undetected clusters were likely small, making the ROK outbreak even more extreme compared to the average. The implications of very high initial variability are 1) a high probability of no transmissions from the index patient, even if R 0 > 1; and 2) a relatively high probability of a superspreading event, i.e., an unusually large number of transmissions, if any do occur. For example, using the MLE (R 0 , k 0 ) from our optimal Model 2b ( Table 2) there would be 77% chance of no transmissions from the initial traveler, but a 5% chance of more than 12 transmissions and a 1% chance of more than 40 transmissions from the initial traveler.
For public health officials in countries anticipating further introductions of MERS-CoV from travelers, it is important to anticipate the non-negligible possibility of an explosive outbreak in early generations of transmission driven by superspreading. There are several reasons that superspreading might occur from an infected individual, including unusually high levels of viral shedding, long length of infectious period, or high numbers of person-to-person contacts, particularly when numerous contacts coincide with peak timing of infectiousness and/or if contacts have unusual susceptibility, such as hospital patients. Investigations of the MERS superspreading events in ROK suggest that patient symptoms (frequent and vigorous coughing) during close proximity with many others in crowded hospital areas contributed to unusually high numbers of transmissions from certain individuals (Oh et al., 2015) .
While the potential for superspreading exists, our results also suggest that a prompt public health response in the early stages of a new outbreak, with efforts to prevent further transmission similar to what has been implemented previously, would most likely Table 1 . Panels C and D assuming 50% of importation clusters were undetected. Panels E and F assuming 75% of importation clusters were undetected. reduce the risk of a very large or long-lasting outbreak to negligible levels. Compared to projections from our optimal model, previously published models extrapolate higher probability of MERS outbreaks that are larger or longer-lasting than what occurred in ROK, but those models did not fully incorporate the rapid decline in transmission rate that was achieved in later generations of the ROK outbreak once it had been identified. Nonetheless, any modelbased extrapolation beyond the data is subject to potentially wide uncertainty and should be interpreted with caution.
Regardless of the true risk posed by infected travelers, the key elements of a coordinated strategy to mitigate new out-breaks of MERS, as with any emerging infection, are continued awareness, targeted surveillance strategies based on importation risk from travelers, appropriately detailed travel histories of ill patients, pre-positioned availability of laboratory diagnostics, and a strong public health response once a potential case is suspected or recognized.
of Infectious Disease Agent Study Grant 5U01GM070694-11, as well as support with resources and the use of facilities at the Department of Veterans Affairs Salt Lake City Informatics, Decision-Enhancement and Analytic Sciences Center (IDEAS 2.0) -Health Services Research and Development #CIN 13-414. The views expressed in this article are those of the authors and do not necessarily reflect the position or policy of the Department of Veterans Affairs or the United States government.