key: cord-0037184-1tkevj8v
authors: Edholm, Christina J.; Emerenini, Blessing O.; Murillo, Anarina L.; Saucedo, Omar; Shakiba, Nika; Wang, Xueying; Allen, Linda J. S.; Peace, Angela
title: Searching for Superspreaders: Identifying Epidemic Patterns Associated with Superspreading Events in Stochastic Models
date: 2018-10-25
journal: Understanding Complex Biological Systems with Mathematics
DOI: 10.1007/978-3-319-98083-6_1
sha: f6202ca147d1ae17734cb820f9c1888848eeb0a2
doc_id: 37184
cord_uid: 1tkevj8v

The importance of host transmissibility in disease emergence has been demonstrated in historical and recent pandemics that involve infectious individuals, known as superspreaders, who are capable of transmitting the infection to a large number of susceptible individuals. To investigate the impact of superspreaders on epidemic dynamics, we formulate deterministic and stochastic models that incorporate differences in superspreaders versus nonsuperspreaders. In particular, continuous-time Markov chain models are used to investigate epidemic features associated with the presence of superspreaders in a population. We parameterize the models for two case studies, Middle East respiratory syndrome (MERS) and Ebola. Through mathematical analysis and numerical simulations, we find that the probability of outbreaks increases and time to outbreaks decreases as the prevalence of superspreaders increases in the population. In particular, as disease outbreaks occur more rapidly and more frequently when initiated by superspreaders, our results emphasize the need for expeditious public health interventions.

Ebola outbreak to date [38] . This spread might have been increased due to infected health-care workers' close contact with susceptible individuals. Additionally, burial ceremonies may increase contact with infectious deceased bodies that contain the virus. The incubation period, defined as the time of infection to onset of symptoms, ranges from 2 to 21 days [38] . Individuals can recover from Ebola; however, mortality rates range from 25 to 90%. In 2016, the WHO announced that the first vaccine trial implemented in Guinea was 100% effective [17, 37] . The recent preventive measures announced by the Centers for Disease Control and Prevention (CDC) include: reducing contacts with infected animals or bodily fluids of infected individuals, isolating infected and deceased individuals, early detection of infected individuals, and maintaining a clean environment [8] .

MERS was first identified in 2012 from an outbreak that occurred in Saudi Arabia [40] . The source of infection was identified as dromedary camels. However, most cases are not due to camel-to-human infections. MERS outbreaks among humans arise from human-to-human interactions, where many cases occur in healthcare settings with poor health prevention and control practices. In 2015, an outbreak of MERS in South Korea was driven by three SS, initiated with one SS contracting MERS during international travel. The first SS was responsible for 29 secondary infections through various clinical visits. Two subsequently infected individuals were responsible for 106 tertiary infections [39, 40] . Individuals infected with MERS can be asymptomatic, while others may experience the following symptoms: fever, coughs, shortness of breath, diarrhea, and pneumonia. Nearly, 35% of MERS cases resulted in death. While no vaccine or treatments are available, individuals are advised to maintain good hygiene when coming into contact with animals, particularly camels, such as washing hands and avoiding contact with sick animals. Additional prevention strategies include consuming thoroughly cooked and prepared animal products [39] .

Mathematical models formulated for recent outbreaks of MERS and Ebola have applied the compartmental setting with various disease stages such as susceptible, exposed, infectious, and recovered (SEIR) or performed statistical analyses to identify important parameters in spread of the disease ( [11, 23, 24] MERS and [4, 10, 16, 22] Ebola). Additional classes for asymptomatic, hospitalized, or isolated individuals were also included in MERS models [11, 23] . Time-dependent transmission parameters accounted for superspreading events (e.g., [22] [23] [24] ). Superspreading events have also been investigated with multitype branching processes by including individual heterogeneity in offspring generating functions [18, 26] . All of these models have contributed to a better understanding of the role of superspreaders in disease outbreaks. Our models incorporate the compartmental framework and apply stochastic simulations with theory from branching processes to further elucidate the role of superspreaders in disease dynamics.

In this investigation, we develop a mathematical modeling framework that incorporates the heterogeneity of hosts through differences in transmission rates to assess the role of SS in disease spread at the population level. Specifically, we aim to study the disease dynamics in a heterogeneous population consisting of SS and NS individuals, and develop a deterministic model based on ordinary differential equations (ODEs) which is expanded to a stochastic model that is implemented as a continuous-time Markov chain (CTMC) system and approximated by a multitype branching process [1, 2] . We incorporate estimated parameter values from published data of prior MERS and Ebola epidemics into our models. Next, we compute the basic reproduction number for the ODE model, and perform sensitivity analysis using Latin hypercube sampling and partial rank correlation. By varying the initial size of SS and model parameters of the CTMC model, we derive and verify analytical estimates obtained using multitype branching process approximations with model simulations to predict the probability of an epidemic outbreak. In further numerical simulations of the CTMC model, we compute sample paths, probability of outbreak, number of deaths, time to outbreak, time to peak infection, and peak number of infectious individuals. Our analyses and numerical simulations reveal how SS influence the dynamics of epidemic outbreaks, which may provide useful insight for public health interventions.

We formulate a simple modeling framework for host heterogeneity due to differences in individuals that account for either SS or NS. In particular, SS or NS may differ in transmission, transitions between disease stages, deaths, recovery, or population size. The SS and NS mix homogeneously, such as in a hospital setting (MERS) or at a large gathering such as a funeral (Ebola). Our basic modeling framework is a system of ODEs with five disease stages for SS and NS as described by the compartmental diagram in Fig. 1 and by the differential equations in (2.1), where i = 1 is NS and i = 2 is SS. The description of the model variables are summarized in Table 1 . Such types of models have been used in metapopulation settings and are referred to as multigroup models (e.g., [25, 34] ).

For the ODE model, we assume that the disease duration is short and, therefore, we do not include birth or natural death rates. In addition, we make the simplifying assumption that NS cannot become SS and vice versa. We make this assumption due to the short duration of the epidemic period and the fact that no control measures are applied (which could change the transmission patterns). In the model, the number of susceptible NS and SS are denoted by S 1 and S 2 , respectively. Susceptible individuals transition into their respective exposed classes, E 1 and E 2 , at a rate of

where β 1 is the transmission rate of the NS asymptomatic A 1 and infective I 1 classes with N 1 as the total number of NS and similarly for SS variables. The total number of individuals is

. From the exposed class, individuals transition to the asymptomatic class, A 1 and A 2 , at a rate of α 1 or α 2 . In the asymptomatic class, there is diseaseinduced mortality with rates μ A1 or μ A2 , respectively. Asymptomatic individuals do not display symptoms but are infectious. Individuals transition into the infective class at a rate of δ 1 or δ 2 , where the disease-induced mortality rates are μ I 1 or μ I 2 . [10] (9.9, 100) × 10 −5 μ ji Disease-induced death rate 0.08 [7] (0.02, 0.14) 0.09 [9] (0.075, 0.125) γ i Recovery rate 0.075 [19] (0.05, 0.1) 0.05 [31, 35] Table 2 . For Ebola parameter values, we used the outbreak in Sierra Leon in 2014 [7, 12, 19, 31] and for the MERS 2015 outbreak in South Korea [9, 10, 31, 35] . Note that these parameter values are taken from a single outbreak of MERS and Ebola, which means that they vary from other outbreaks and may present some constraints when asserting conclusions for outbreaks of the same infectious disease [14, 20] . However, the parameter values used from the two outbreaks provide an excellent baseline for our model simulations.

We compute the reproduction number for the ODE system (2.1) using the nextgeneration matrix [34] . The basic reproduction number, R 0 , is defined as the number of secondary cases produced by the introduction of a single infected individual into a fully susceptible population. If R 0 > 1, an outbreak occurs in the ODE model. We start by defining two matrices, F and V , where the F matrix represents the newly infected rates in the system, and V represents the remaining rates in the infected compartments, Eqs. (1) and (2) respectively. The matrix F −V is the Jacobian matrix of the infected compartments evaluated at the disease-free equilibrium (DFE), wherē S i = N i (0) and E i = A i = I i = R i = 0, i = 1, 2. We find the spectral radius of the matrix F V −1 (Appendix 1), which equals the basic reproduction number,

(2.

2)

The basic reproduction number has the form typical of a multigroup/stage progression model [34] . It is the sum of two basic reproduction numbers, one for each group, NS when i = 1 and SS when i = 2. In particular,

The two terms in the preceding expression represent new infections resulting from either the asymptotic stage A i or from the infectious stage I i , i = 1, 2. In addition for group i, the term β i (N i /N ) is the number of successful transmissions from an individual in stage A i (first term) or from an individual in stage I i (second term) that result in exposed individuals. The term 1/(δ i + μ A i ) is the average length of the asymptotic stage while 1/(γ i + μ I i ) is the average length of the infectious stage, and δ i /(δ i + μ A i ) is the probability of transitioning from A i to I i . For parameter values in Table 2 and for equal proportion of SS and NS, N 1 /N = 0.5 = N 2 /N, the basic reproduction number for MERS is R 0 = 2.36 and for Ebola it is R 0 = 3.75.

If the number of hosts/pathogens is sufficiently small, an ODE model is not appropriate. To that end, we utilize a continuous-time Markov chain (CTMC) model, which is continuous in time and discrete in the state space, to study the variability at the initiation of an outbreak, in time to outbreak, and in the peak level of infection. For simplicity, we use the same notation for the state variables as in the ODE model. In particular, time t ∈ [0, ∞) and the states are discrete random variables, e.g.,

The Markov property implies that the future states of the stochastic process only depend on the current states. In particular, there is an exponential waiting time between events. To formulate a CTMC, it is necessary to define the infinitesimal transition probabilities corresponding to each change (event) in the state variables. The CTMC model consists of 12 distinct events, six events for each of the groups, NS and SS. The changes and the corresponding infinitesimal transition rates are summarized in Table 3 .

The theory of multitype (Galton-Watson) branching processes has a long history (e.g., [2, 13, 36] and references therein). It has been used to approximate the dynamics of the CTMC model near the DFE and the stochastic threshold for a disease outbreak [1] [2] [3] 36] . In fact, the stochastic threshold (i.e., probability of a disease outbreak) is directly related to the basic reproduction number as defined in Table 3 State transitions and rates for the CTMC model with Poisson rates

the corresponding deterministic model (2.1) (see [3, 36] ). More specifically, if the basic reproduction is less than unity, then disease extinction occurs with probability one. In this case, the branching process is called subcritical. However, if the basic reproduction number is greater than unity, the probability of disease extinction is less than one (probability of outbreak is greater than zero) and the process is referred to as supercritical.

In what follows, we will apply a multitype branching process approximation of the CTMC model at the DFE to estimate disease extinction probability. First, let us define the offspring probability generating function (pgf) for the exposed, asymptomatic, and infectious individuals in NS and SS. Let X = (X 1 , . . . , X 6 ) := (E 1 , A 1 , I 1 , E 2 , A 2 , I 2 ) be a vector of integer-valued random variables and δ ij denote the Kronecker delta (i.e., δ ij = 1 if i = j and zero otherwise). In general, the offspring pgf for type i given X j (0) = δ ij is a function from [0, 1] 6 to [0, 1], and it takes the form:

Here, P i (k 1 , k 2 , . . . , k 6 ) is the probability that the individual of type i gives "birth" to k j individuals of type j for j = 1, 2, . . . , 6. In particular, the pgfs f i : 6) are given by:

According to the theory of multitype branching processes [5, 13] , the fixed points of the offspring pgfs give an estimate of the disease extinction probability. Let (q 1 , q 2 , q 3 , q 4 , q 5 , q 6 ) be the minimal fixed points of pgfs; that is, f i (q 1 , . . . , q 6 ) = q i for i = 1, . . . , 6. Then, an estimate of the extinction probability given

and hence the probability of an outbreak is

However, due to the simplicity of f 1 and f 4 (no deaths during stage E i ), the pgfs can be simplified. That is, x 1 = x 2 and x 4 = x 5 . Therefore, we only solve for q 2 , q 3 , q 5 , and q 6 . The expectation matrix M = (m ij ) can be shown to be directly related to the basic reproduction number [3] with m ij = ∂f j ∂x i | X=1 . We include this calculation in Appendix 2. It is known that the spectral radius of M, denoted as ρ(M), determines whether the disease extinction probability is equal to or less than the unity [2, 3, 13] . Specifically, if ρ(M) < 1, q 1 = · · · = q 6 = 1, then the extinction probability is one; if ρ(M) > 1, then there exists a unique fixed point (q 1 , · · · , q 6 ) ∈ (0, 1) 6 , and hence the extinction probability is strictly less than one. By the Threshold Theorem of reference [3] , it follows that the spectral radius of the matrix M is strictly less than one if and only if the basic reproduction number is strictly less than one. Analogous statements hold whenever the spectral radius of M is equal to one or is strictly greater than one.

We perform a sensitivity analysis on the parameters ranges given in Table 2 for the ODE models for MERS and Ebola using a uniform distribution for the values. Latin hypercube sampling (LHS), first developed by McKay et al. [29] , with the statistical sensitivity measure partial rank correlation coefficient (PRCC), performs a sensitivity analysis that explores a defined parameter space of the model. The parameter space considered is defined by the parameter intervals depicted in Table 2 . Rather than simply exploring one parameter at a time with other parameters held fixed at baseline values, the LHS/PRCC sensitivity analysis method globally explores multidimensional parameter space. LHS is a stratified Monte Carlo sampling without replacement technique that allows an unbiased estimate of the average model output with limited samples. The PRCC sensitivity analysis technique works well for parameters that have a nonlinear and monotonic relationship with the output measure. PRCC shows how the output measure is influenced by changes in a specific parameter value when the linear effects of other parameter values are removed. The PRCC values were calculated as Spearman (rank) partial correlations using the partialcorr function in MATLAB 2016. Their significances, uncorrelated p-values, were also determined. The PRCC values vary between −1 and 1, where negative values indicate that the parameter is inversely proportional to the output measure. Following Marino et al. [27] , we performed a z-test on transformed PRCC values to rank significant model parameters in terms of relative sensitivity. According to the z-test, parameters with larger magnitude PRCC values had a stronger effect on the output measures.

We start by verifying the monotonicity of the output measures. Monotonicity was observed for all parameters except μ I 2 with total SS deaths, which exhibited two monotonic ranges [0.02, 0.0278] and [0.0278, 0.14]. For non-monotonic trends, alternative methods based on decomposition of model output variances such as eFAST (extended Fourier Amplitude Sensitivity Test) can be used instead of PRCC [27] ; however, since all other parameters were monotonic, we use PRCC and just consider the two monotonic ranges of μ I 2 separately. PRCC analysis of these two ranges produces similar results. For an analysis of the monotonicity, refer to Appendix 3. Once meeting the monotonicity requirements, we proceed to utilize LHS with PRCC for both MERS and Ebola parameters. For each disease, we calculate the PRCC for the following output measures: total NS cases, total SS cases, total NS deaths, and total SS deaths. The number of total cases refers to the total number of transmission events where susceptible individuals become exposed (latently infected) individuals. For the outputs of NS/SS cases, the PRCC results were similar in both Ebola and MERS. According to the PRCC values, the β 2 and μ I 2 are significant in the model for MERS. Meanwhile, in the Ebola model, both transmission parameters are significant in the model (see Fig. 2 ). Note that β 2 is calculated from R 0 , which we will vary later in simulations.

For the CTMC model, we numerically simulate sample paths to compute the probability of an outbreak, number of deaths, time to outbreak, time to peak infection, and peak number of infectious individuals. For sample paths and probability of outbreak, we compare our results with the deterministic model. In the remainder of this analysis, we assume that the initial total population size is N(0) = 2000. Reference to infected individuals will imply the variables I 1 and I 2 , unless stated otherwise. For example, peak number of infectious individuals refers to the maximum value of The number of total cases refers to the total number of transmission events where susceptible individuals become exposed (latently infected) individuals. P -values that are greater than 0.05 are labeled as not significant (n.s.) I 1 + I 2 and the time to peak infection refers to the time t at which this maximum occurs. However, an outbreak means that the total number in classes E i , A i , and I i for both NS and SS has reached at least 50, i.e., (E i + A i + I i ) ≥ 50. In addition, we note that for the CTMC model, an outcome measure (e.g., peak values, time to peak, and number of deaths) is defined by a corresponding probability distribution and a sample path yields one outcome from the distribution.

An example of the sample paths resulting from our CTMC model is shown in Fig. 3 , for both MERS and Ebola cases. These sample paths are generally well aligned with the population average response that is captured by our ODE model (shown with a black line). However, the sample paths of the CTMC model illustrate the potential variability in timing of the peak level of infection and the peak number of infectious individuals. Note that some sample paths are not shown because in those simulations the disease becomes extinct. Also, note that the A class is not shown for Ebola (Fig. 3b) given that the asymptomatic stage is extremely short for this disease.

Next, in order to do a comprehensive comparison of the stochastic simulation and ODE model results, we probe the relationship between two model parameters-the value of R 0 (Fig. 4) as well as the fraction of the susceptible population that are in the SS class (Fig. 4 ) and a key model output: the probability of outbreak. Probability of outbreak is defined by monitoring the number of people in the E, A, and I classes and an outbreak is declared when the cumulative size of these compartments reaches the threshold value of 50. Although the value of 50 appears relatively large, it is reasonable given that we are counting the cumulative number in all three classes for a relatively large population size of 2000. For these simulations, we vary β 2 given the significant effect of this parameter on the model outputs as confirmed by the LHS analysis. We note a negative correlation between the proportion of SS in the S class and the probability of outbreak (Fig. 4) and attribute this to the fact that the value of β 2 is varied in order to maintain a constant value of R 0 (MERS, R 0 = 2.5 and Ebola, R 0 = 2.39). In other words, as the fraction of SS susceptible individuals is increased, the value of β 2 decreases and results in a reduction in the probability of outbreak (1 − q 6 ). Results in Fig. 4 are shown only for q 3 and q 6 since these outputs are similar to q 1 and q 4 , respectively. We also note that q 1 = q 2 and q 4 = q 5 and therefore exclude those plots as well.

As expected, the probability of an outbreak is dependent on the initial fraction of the population that is infected, with an increasing chance of an outbreak (Fig. 4) . Furthermore, the probability of outbreak is significantly enhanced when the initially infected population is composed of SS rather than NS individuals. We also find a strong agreement between the probability of outbreak predicted by stochastic simulations of the CTMC model and the associated branching process approximations for all of these analyses (Fig. 4a-f ). 

Utilizing our stochastic model of MERS and Ebola dynamics within a population of individuals, we next sought to investigate whether the presence of SS individuals within the population could be reflected in key metrics that capture the severity of disease outbreak: the number of deaths, time to disease outbreak, probability of outbreak, time to peak number of infections, and the peak number of infectious individuals.

We first assess the impact of SS individuals on the number of deaths that accumulate over a 150-day time frame following disease initiation. We observe a modest increase in the frequency of deaths as the size of the susceptible SS class of individuals is increased from 5 to 50% of the total population for both MERS and Ebola disease simulations (not shown). We note a higher frequency of epidemics For all subsequent simulations, we initialize the population consisting of 1000 SS and 1000 NS susceptible individuals. Most notably, there is a ten-fold increase in the frequency of deaths expected when the initial infected individual (for both MERS and Ebola) is an SS rather than an NS (Fig. 5) . The statistical significance of the difference between NS-and SS-initiated epidemics is confirmed with a Kolmogorov-Smirnov test (p < 0.001) [28] . It is clear that the distributions are bimodal. This is due to the fact that there may be only a minor outbreak (with probability q 3 or q 6 ) with none or a few deaths or a major outbreak (with probability 1 − q 3 or 1 − q 6 ) with a significant number of deaths.

We next explore the relationship between the number of deaths and the number of initially infected individuals. For each fraction of the population initially infected, there are 1000 points, one point from each of the 1000 sample paths, representing the total number of deaths over a 150-day time period. As expected, we find that as the number of initially infected NS individuals increases, the expected number of deaths increases as well (Fig. 6) . Interestingly, we find a threshold response as the fraction of initially infected SS individuals increases. As the fraction of initially infected SS individuals increases beyond 0.005 for MERS (Fig. 6b) and 0.0075 for Ebola (Fig. 6d) , we find that the simulation always gives rise to an outbreak, resulting in a maximal number of around 1000 deaths over a 150-day simulated period. We also note that there is a decrease in the variability in the number of deaths when the outbreak is initiated by an SS rather than an NS, which contributes to this threshold response. The seemingly binary response in the number of deaths resulting from a MERS or Ebola epidemic initiated by infected SS individuals who only contribute to 0.5-0.75% of the starting population is a good indication that by tracking the number of deaths in an epidemic, the presence of an SS may be predicted. Thus, while the observation that an outbreak has occurred does not necessarily suggest the existence of SS individuals in the population, the severity of the outbreak in terms of lives lost may be more suggestive of the presence of an SS, especially when the number of known initial infections is low.

Similarly, we find that the time to outbreak-where an outbreak is defined as 50 or more people in all of the E, A, and I classes-is reduced when the initial infected individual in a simulated MERS or Ebola disease situation is an SS rather than an NS (Fig. 7a, c) . We confirmed that this reduction is, indeed, statistically significant ( Fig. 8a-b) . These results also illustrate that as the fraction of susceptible SS increases the time to outbreak increases as well, which we attribute to the fact that β 2 values decrease (detailed in Figs. 7 and 8) . In Fig. 7 , each distribution is based on 10,000 sample paths, whereas in Fig. 8 , for each fraction initially infected, the time points are based on 1000 sample paths.

We also find a clear separation between the time to outbreak of an epidemic initiated by a fraction of SS versus NS infected individuals. Mean differences were significantly distinct for each percentage in (Fig. 8a-b) , p < 0.001. In fact, if a MERS or Ebola outbreak is initiated by 1.5% or more of the initial population size and these individuals are SS, then the time to outbreak is predicted to be no more than 20 days where an outbreak is defined as 2.5% of the population becoming infected (Fig. 8c-f ). The fraction of susceptible SS is increased for comparison. All results were statistically significant when p < 0.05 from t-test (two-tailed). Comparisons that are not statistically significant were denoted n.s.

Given that the time to outbreak shows a significant difference between epidemics initiated by SS versus NS individuals, we next asked whether SS-initiated epidemics will also reach peak infection in a shorter time. To investigate this, we calculated mean (± SD) of time to peak infection (in days) for MERS and Ebola, where the percent of SS varied in the susceptible population (see Fig. 9a-b) . Mean differences between the introduction of 1 infected NS (black) compared to 1 infected SS (white) were assessed separately as the percent of SS varied (e.g., 5%, 25%, and 50%) using t-tests where statistical significance was accepted when p < 0.05. For MERS, time to peak infection was slightly significantly lower for SS when 5% and 25% of the susceptible population was SS, but not significant when 50% of the population was SS (Fig. 9a) . For Ebola, time to peak infection was significantly lower for SS regardless of changes in the percent of SS in the susceptible population (Fig. 9b) . Hence, while the differences in mean time to peak infection between SS and NS-initiated epidemics are only modestly different, we find their difference to be statistically significant (Fig. 9a-b) . Thus, this confirmed that epidemics initiated by infected SS individuals reaches its peak value more quickly. We repeated the same analysis to assess mean differences in the peak number of infections. Surprisingly, we did not find a significant difference between the peak number of infections for epidemics initiated from a single infected NS versus SS individual for Ebola (Fig. 9d) . However, significant differences were observed when 5% or 50% of the susceptible population was SS for MERS.

In this investigation, we capture the dynamics of MERS and Ebola epidemics by applying both deterministic and stochastic modeling strategies. To investigate the role of SS on the epidemic dynamics and to compare our results, we keep the R 0 constant for both MERS and Ebola while varying β 2 , the transmission rate of SS. Parameter sensitivity analysis, using Latin hypercube sampling and partial rank correlation coefficient, shows that β 2 has a significant effect on all the output measures (Fig. 2) .

From Fig. 4 , we can conclude that the stochastic model simulations agree with the branching process analytical results. As the value of R 0 increases, we observe that the probability of an outbreak increases for both diseases. This result is expected since more individuals in the population are infected. The probability of an outbreak is greater for Ebola than MERS, which is due to the transmission parameters for Ebola being larger than MERS. Furthermore, these results show that if the outbreak is initiated by an SS, then the probability of an outbreak is significantly higher. Additionally, fewer SS individuals than NS individuals are sufficient to cause an outbreak irrespective of the disease (MERS or Ebola).

As an outbreak initiated by SS has a greater probability of occurrence and peaks earlier than with NS, the accumulated number of deaths is more severe in an epidemic initiated with the same proportion of SS than NS (Figs. 5 and 6). Disease severity (number of deaths) for both MERS and Ebola occurs earlier with SS than NS. Our findings agree with prior epidemiological studies on superspreading events [15, 32, 40] . For example, the 2003 outbreak of the respiratory infection SARS in Beijing found that SS had higher mortality rates, higher attack rates, and greater number of contacts in comparison to NS [40] . From a public health perspective, as SS events will be observed more frequently, intervention/prevention methods must have rapid response to reduce disease severity. For example, Wong et al. [40] suggested that several community-based efforts could have been made to reduce the number of MERS and Ebola cases in Guinea and Sierra Leone, such as tracking contacts, earlier diagnosis, treatment strategies, and community education. Effective responses to control superspreading events and reduce disease transmission in MERS and Ebola outbreaks included: "early discovery, diagnosis, intervention, and quarantine of confirmed cases." [40] . Other epidemics that are more likely to occur in hospital settings, e.g., SARS, could be controlled through hospital administrative strategies, such as reducing contact between the infected patient and healthcare workers, visitors, or other patients whose immune system may be comprised due to other infections [32] . Thus, a rapid response is needed to reduce disease severity of SS events.

Evident in Figs. 7 and 8 , when an outbreak is initiated by an SS rather than an NS, the time to outbreak is shorter and has less variability. Therefore, if the number of disease cases rises rapidly, there may be SS in the community. In this scenario, healthcare managers should search for potential SS. Similar results apply for time to peak infection, Fig. 9 . If peak infection occurs quickly, it is more likely that there is an SS in the population.

Interestingly, varying the percentage of SS in the population has little influence on the peak number of infections (Fig. 9c, d) . This is likely due to the fact that the R 0 values are held constant.

We have formulated, analyzed, and numerically simulated deterministic and stochastic epidemic models that include heterogeneity in transmission for NS and SS. We applied our models to emerging and re-emerging infectious diseases, MERS and Ebola, where the models were parameterized with data from the literature but with a fixed initial population size of 2000. There are a number of extensions and generalizations that we will consider in the future work. We assumed homogeneous mixing and only two types of classifications of individuals (NS/SS) for the entire population. Generalizing this model to include heterogeneous mixing and spatial components are key features that can provide insight on how a superspreaders can be classified. In our model, we considered inter-host variability, which naturally leads to constructing a model with intra-host variability utilizing stochastic differential equations or other types of models. In addition, variability of the pathogen on epidemic dynamics can be explored. Additionally, we will validate our models' findings against time series data, test our models' abilities to detect the presence of SS, and interpret the results for public health implementation. Finding answers to these problems will lead to our ultimate goal of constructing novel ways to quantify, characterize, and identify an SS during the initiation of an outbreak.

Women Through Research-Focused Networks" (NSF-HRD 1500481), Society for Mathematical Biology, and Microsoft Research. We give special thanks to the WAMB organizers: Ami Radunskaya, Rebecca Segal, and Blerta Shtylla. Anarina L. Murillo acknowledges that this work has been supported in part by the grant T32DK062710 from the National Institute of Diabetes and Digestive and Kidney Diseases and grant T32HL072757 from the National Heart, Lung, and Blood Institute. Omar Saucedo acknowledges that this research has been supported in part by the MBI and the grant NSF-DMS 1440386. Nika Shakiba is the recipient of the NSERC Vanier Canada Graduate Scholarship. This work was partially supported by a grant from the Simons Foundation (#317047 to Xueying Wang). In addition, we thank Texas Tech University for hosting our second WAMB group meeting and the Paul Whitfield Horn Professorship of Linda JS Allen for providing financial support. We thank the two anonymous reviewers for their helpful suggestions on the original manuscript.

In the calculations below, we denote N i (0) as N i and N(0) = N 1 (0) + N 2 (0) as N :

and

The next-generation matrix is

For each pgf, we take partial derivatives with respect to x 1 , . . . , x 6 (Jacobian matrix of the functions f 1 , . . . , f 6 ), then evaluate at x 1 , . . . , x 6 = 1 and take the transpose of the matrix. The result is the expectation matrix, M, where Eq. (3) reduces to (4) .

In addition, we create the matrix W , Eq. (5), a diagonal matrix.

and

Finally, we check that

For the MERS parameters, the graph of the output measure with μ I 2 parameter had a concave curvature implying that another sensitivity analysis maybe implemented. However, the range of the output measure is small enough that we can ignore the monotonicity. The remaining graphs for MERS (Fig. 10 ) and all the graphs for Ebola (Fig. 11 ) are all monotonic which means that we can trust the sensitivity analysis and proceed to the PRCC analysis (Fig. 12 , Tables 4, 5 and 6). 

Stochastic Population and Epidemic Models

Extinction thresholds in deterministic and stochastic epidemic models

Relations between deterministic and stochastic thresholds for disease extinction in continuous-and discrete-time infectious disease models

Estimating the reproduction number of Ebola virus (EBOV) during the 2014 outbreak in West Africa

Branching Processes

Dynamically modeling SARS and other newly emerging respiratory illnesses: past, present, and future

Mortality Risk Factors for Middle East Respiratory Syndrome Outbreak

Ebola Virus Disease) Prevention

Ebola Outbreak in West Africa-Case Counts

Transmission dynamics and control of Ebola virus disease (EVD): a review

Synthesizing data and models for the spread of MERS-CoV, 2013: key role of index cases and hospital transmission

Preliminary epidemiologic assessment of MERS-CoV outbreak in South Korea

In the garden of branching processes

A pandemic risk assessment of middle east respiratory syndrome coronavirus in Saudi Arabia

Epidemiology: dimensions of superspreading

Assessing the international spreading risk associated with the 2014 West African Ebola outbreak

Efficacy and effectiveness of an rVSVvectored vaccine expressing Ebola surface glycoprotein: interim results from the Guinea ring vaccination cluster-randomised trial

An event-based model of superspreading in epidemics

The daily computed weighted averaging basic reproduction number for MERS-CoV in South Korea

Estimating the basic reproductive ratio for Ebola outbreak in Liberia and Sierra Leone

A deterministic model for gonorrhea in a nonhomogeneous population

Spatial and temporal dynamics of superspreading events in the 2014-2015 West Africa Ebola epidemic

A dynamic compartmental model for the Middle East respiratory syndrome outbreak in the Republic of Korea: a retrospective analysis on control interventions and superspreading events

Statistical inference in a stochastic epidemic SEIR model with control intervention: Ebola as a case study

Spatial heterogeneity in epidemic models

Superspreading and the effect of individual variation on disease emergence

A methodology for performing global uncertainty and sensitivity analysis in systems biology

The Kolmogorov-Smirnov test for goodness of fit

Comparison of three methods for selecting values of input variables in the analysis of output from a computer code

Compartmental disease models with heterogeneous populations: a survey

Modeling the impact of interventions on an epidemic of Ebola in

Super-spreaders in infectious diseases

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Time from infection to disease and infectiousness for Ebola virus disease, a systematic review

The outcome of a stochastic epidemic-a note on Bailey's paper

WHO, Final trial results confirm Ebola vaccine provides high protection against disease

Ebola Virus Disease Fact Sheet

East respiratory syndrome coronavirus (MERS-CoV) Fact Sheet

the role of super-spreaders in infectious disease

Heterogeneities in the transmission of infectious agents: implications for the design of control programs