key: cord-1014808-v8181t5n
authors: Shakiba, Nika; Edholm, Christina J.; Emerenini, Blessing O.; Murillo, Anarina L.; Peace, Angela; Saucedo, Omar; Wang, Xueying; Allen, Linda J.S.
title: Effects of environmental variability on superspreading transmission events in stochastic epidemic models
date: 2021-03-18
journal: Infect Dis Model
DOI: 10.1016/j.idm.2021.03.001
sha: 950c22ac6c672582426d8e23a2d8e6a55ec9e8ea
doc_id: 1014808
cord_uid: v8181t5n

Superspreaders (individuals with a high propensity for disease spread) have played a pivotal role in recent emerging and re-emerging diseases. In disease outbreak studies, host heterogeneity based on demographic (e.g. age, sex, vaccination status) and environmental (e.g. climate, urban/rural residence, clinics) factors are critical for the spread of infectious diseases, such as Ebola and Middle East Respiratory Syndrome (MERS). Transmission rates can vary as demographic and environmental factors are altered naturally or due to modified behaviors in response to the implementation of public health strategies. In this work, we develop stochastic models to explore the effects of demographic and environmental variability on human-to-human disease transmission rates among superspreaders in the case of Ebola and MERS. We show that the addition of environmental variability results in reduced probability of outbreak occurrence, however the severity of outbreaks that do occur increases. These observations have implications for public health strategies that aim to control environmental variables.

Infectious diseases are one of the biggest threats to human life and a major concern to public health (Tatem et al., 2006) . The spreading of emerging and re-emerging diseases, such as Middle East Respiratory Syndrome (MERS) and Ebola virus disease (Ebola), have been attributed to population mobility and the interconnectedness of the world's heterogeneous population (Tatem et al., 2006) . Superspreader (SS) is a 5 term used to informally describe high-risk individuals who disproportionately spread infection, because they have increased contacts with susceptible individuals, they are carriers of a more contagious strain, and/or both specific and general relationships between environmental change, infectious disease and their spread (Eisenberg et al., 2007) . Hence, it is established that the environment can be a major factor in superspreading 45 events. For example, crowding, poor ventilation, improper isolation procedures, unnecessary movement of the infectious, and misdiagnosis have all been identified as environmental factors that contributed to MERS (Chan et al., 2015) and Ebola (Glynn et al., 2017) outbreaks.

Building on our prior research, we extend our stochastic models for MERS and Ebola to investigate the effect of environmental variability in transmission rates on the probability of initiating an outbreak as 50 well as the severity of outbreaks (Edholm et al., 2018) . Accounting for this variability results in temporally correlated transmission parameters that vary randomly and continuously with the environment.

In the following sections, we review briefly some background on MERS and Ebola, before we formulate two new stochastic models that include demographic and environmental variability. We discuss advantages and disadvantages of the two stochastic modeling approaches. In addition, we summarize the results of 55 our extensive numerical simulations comparing model outputs on designated measures: probability of an outbreak, time to the outbreak, number of deaths, peak number of infections reached in the outbreak, and time to the peak infection.

MERS coronavirus belongs to a family of enveloped, single-stranded RNA viruses (genus Betacoronavirus), that infect a number of different species, including humans, by predominantly causing mild selflimiting upper respiratory tract infections and other related illnesses. Other well-known viruses from the same genus include severe acute respiratory syndrome coronavirus (SARS-CoV) and severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which cause the two diseases SARS and COVID-19, respectively (Peeri et al., 2020) . MERS was first identified in 2012 from an outbreak in Saudi Arabia (Chan et al., 2015) .

MERS transmission arises from human-to-human interactions, where infected individuals can be asymptomatic or symptomatic (WHO, 2017 ). An outbreak in South Korea in 2015 with 166 cases was identified to be driven by three SS. The initial infected individual was a SS responsible for 29 secondary infections, of which two individuals were also SS who were responsible for 106 subsequent infections (WHO, 2017) .

Ebola virus is a member of the Filoviridae virus family, known as one of the emerging and re-emerging 70 zoonotic pathogens, causing acute hemorrhagic fever with high case fatality rate in humans (Beeching et al., 2014; Nii-Trebi, 2017) . Ebola was first reported in 1976 in Africa during the Ebola outbreak in the Democratic Republic of the Congo. The transmission of Ebola virus results from direct contact with bodily fluids containing virus particles from infected patients or diseased individuals, as well as from contact with contaminated objects. Epidemiological studies have revealed that family members of infected persons and 75 associated health care workers are potential spreaders of the virus since they are likely to form a link between the infected persons and susceptible individuals. Though people can recover from Ebola infection the mortality rate is still very high (WHO, 2020).

We extend our previous stochastic SEAIR model, which included demographic variability, to incorporate 80 environmental variability in the disease transmission rate (Edholm et al., 2018) . We did this by assuming an average transmission rate and applying a mean-reverting stochastic differential equation for both MERS (Figure 1(a) ) and Ebola (Figure E.1(a) in Appendix E). This led to the development of two novel models: a time-nonhomogenous stochastic process (NHP) model with discrete random variables, and a stochastic differential equation (SDE) model with continuous random variables. As with our previously published 85 CTMC model (Edholm et al., 2018) , these new models enabled us to track several key measures of disease dynamics, including: time to the outbreak, number of deaths, peak number of infections reached in the outbreak, time to the peak infection, as well as probability of an outbreak (when total infections in the population exceeds an outbreak threshold of 100) (Figure 1(b) ). Our models do not account for the emergence of new pathogens, which would have allowed for NS individuals to become SS depending on the pathogenic 90 variance. Rather, they focus on physiological and behavioral causes of superspreading.

The variables S i , E i , A i , I i and R i denote the number of susceptible, exposed, asymptomatic, infected and recovered individuals of class i, respectively, where i = 1 is NS (nonsuperspreader) and i = 2 is SS (superspreader). The underlying ODE framework for MERS can be seen in the compartmental flow diagram 95 in Figure 1 (a). The differential equations for the ODE model are given in Appendix A. Parameter N is the total population size and β i (t) is the time-dependent transmission rate for i = 1, 2. If β i (t) = β i is constant, then the basic reproduction number for the ODE and CTMC models is

where N 1 + N 2 = N . A description of the parameters and their values are summarized in Table 1 . We assume the basic reproductive number is R 0 = 2.5 and use this to parameterize the mean transmission rate 100 for superspreaders β 2 . We note that the value of R 0 = 2.5 is in the absence of control measures, in contrast to commonly reported effective reproductive numbers of 0.3-0.5 for MERS outbreaks (Majumder et al., 2014) and 1.3 to 1.9 for Ebola outbreaks, e.g., (Chowell et al., 2004; Khan et al., 2015) . In addition, we assume a fixed proportion of SS and NS (20% SS and 80% NS out of a total population size of 2000). A full exploration of the effects of initial conditions and proportion of SS in the population in our designated model outputs 105 was presented in our previous paper focusing on demographic variability (Edholm et al., 2018) . Our previous stochastic modeling with MERS and Ebola found that the SS proportion did not significantly affect time to peak infection or peak number of infections when initiated either by one infected SS or one infected NS (Edholm et al., 2018) . 

% " (')) " *# + *" %"(')," *# + *" %#('),# *# + *" %"('))" *# + *" %"(')," *# + *" %#('),# *# + *" ** ** ** Figure 1 : A time-nonhomogeneous process that incorporates environmental variability in transmission rates.

(a) A schematic representation of our compartmental model for MERS disease dynamics in which individuals pass through the susceptible (S i ), exposed (E i ), asymptomatic (A i ), infected (I i ), and recovered (R i ) classes, where i = 1 corresponds to non-superspreaders and i = 2 corresponds to superspreaders. The model incorporates demographic variability, as before Edholm et al. (2018) , and environmental variability is incorporated through the use of a mean-reverting stochastic differential equation in which the disease transmission rate (β i (t)) changes continuously over time. We assume that α 1 = α 2 , δ 1 = δ 2 , and µ A1 = µ A2 = µ I1 = µ I2 . Estimated parameters are marked by * * for assumed R 0 = 2.5 (Edholm et al., 2018) . †For a total population size of N = 2000, N 1 = 0.8N , and N 2 = 0.2N are assumed.

First, we describe how demographic variability is incorporated in the MERS NHP and SDE models in 110 Sections 3.1.1 and 3.1.2, respectively. Second, we describe how environmental variability is incorporated through the transmission rates in Section 3.1.3.

The MERS NHP model is based on the 12 discrete changes in Table 2 

The MERS SDE model is based on a diffusion approximation, where the first two moments of the process are approximated from the transition rates in Table 2 (Allen, 2007; Allen et al., 2008; Kurtz, 1970 Kurtz, , 1971 Kurtz, , 1972 . The result is a system of Itô SDEs. In particular, the SDE model has the general form:

where X is the vector of random variables. The drift vector F agrees with the underlying ODE model, the 125 diffusion matrix G accounts for demographic variability and W is a vector of independent Wiener processes.

In the derivation, to order ∆t, F (X(t), t)∆t ≈ E[∆X(t)] and G(X(t), t)G(

with ∆X(t) = X(t + ∆t) − X(t) (Allen, 2007; Allen et al., 2008) . As ∆t → 0, the SDE takes the form of model (1).

J o u r n a l P r e -p r o o f 

Rate, a(t)

In the MERS SDE model, there are 10 state variables and 12 distinct events. We use the same notation for the variables as we did for the ODE and NHP models, specifically

Also, W (t) = (W 11 (t), . . . , W 16 (t), W 21 (t), . . . , W 26 (t)) T is a vector of 12 independent one-dimensional Wiener processes. The first index i of W ij (t) represents NS (i = 1) or SS (i = 2) and the second index j is one of the six different events j = 1, . . . , 6 for NS or SS. The drift vector F (X(t), t) has dimension 10 × 1 and is given by the rates in the ODE model. The diffusion matrix G can be expressed as a 10 × 12 matrix (10=number of random variables and 12=number of events) (Allen, 2007; Allen et al., 2008) :

where O is a 5 × 6 zero matrix and C N S is a 5 × 6 matrix with the following form:

where all variables are evaluated at t and

Matrix C SS has a similar form but with the role of NS and SS interchanged (subscript 1 is replaced by 2).

The order of the six different events for NS is the same as that described in Table 2 .

The explicit Itô SDE model for demographic stochasticity in a MERS epidemic consists of a system of 10 equations:

for i = 1, 2. For simplicity, the dependence on t has been omitted from the state variables, e.g., S i ≡ S i (t).

To ensure that solutions are nonnegative and that the dynamics follow that of the CTMC when solutions 135 are near zero, we modify this model slightly. We make the assumptions that all transition rates given in Table 2 are nonnegative and that specific terms in the diffusion matrix are zero when the random variables are close to zero (Allen et al., 2020) . See details in Appendix B.

The description of the NHP and SDE models for Ebola are similar to those for the MERS models, with the exception that there is no asymptotic stage for Ebola. The Ebola models are described in Appendix C.

The similarity between the overall dynamics of the two diseases is consistent with the epidemiological compartmental models which have been previously studied, focusing on different parameters to differentiate the diseases (Wong et al., 2015; Chowell et al., 2016) . In terms of compartmental structure the difference between the two diseases is the presence of an asymptomatic class for MERS, which from previous models is not present for Ebola (Edholm et for the stochastic MERS and Ebola models is described next.

Environmental variability affects model parameters in a variety of ways (Allen et al., 2006; Cai et al., 2018; Marion et al., 2000; Truscott and Gilligan, 2003; Varughese and Fatti, 2008) . For example, environmental 150 fluctuations have been included in the carrying capacity of a zoonotic reservoir (Allen et al., 2006) , the death rate in an epidemic model (Cai et al., 2018) , and the parasite egg development in a helminth infection (Marion et al., 2000) . The parameters in these models, carrying capacity, death rate, and egg development, were 8 J o u r n a l P r e -p r o o f modeled by a mean-reverting Ornstein-Uhlenbeck (OU) process. A disadvantage of the OU process is that the parameter governed by an OU process has an asymptotic normal distribution, and therefore, nonnegativity 155 of the parameter is not preserved (Allen, 2016; Iacus, 2009) . Several alternative mean-reverting processes have been proposed that are more realistic than the OU process and preserve nonnegativity (Allen, 2016 ).

The Cox-Ingersoll-Ross (CIR) process is used to model the continuous and random environmental effects on the transmission parameters (Iacus, 2009 ). The CIR process takes the form of the following SDE:

The mean and variance of β i (t) are

The mean value ofβ i was the constant transmission parameter used for NS or SS, respectively, in the CTMC 160 models for MERS and Ebola by Edholm et al. (2018) . (See Table 1 .)

The advantages of the CIR process over the OU process are that the CIR process has solutions that are nonnegative and has a well-known and commonly-used asymptotic distribution, a gamma distribution with constant meanβ i and varianceβ i σ 2 i /(2r i ). The coefficient of variation CV (standard deviation divided by the mean) for the asymptotic gamma probability density function (pdf) equals

The asymptotic pdf depends on the parameters,β i , r i , and CV . As the value of CV increases (CV > 1) the tail of the distribution becomes fatter, meaning that there is a greater probability for the transmission rate to have large values. See Figure 2 (a). Formally, the parameter r i is the rate of return to equilibriumβ i (Iacus, 2009 ). More specifically, the parameter r i is related to the fluctuations in β i (t) over time. For a small value of r i there is greater temporal correlation between the transmission values at time t and t + ∆t so that the 170 transmission rate changes slowly even when the CV is large. Conversely, for a large value of r i the temporal correlation between two consecutive times is much smaller, and the two values can have widely differing values when the CV is large. See Figure 2 (b) which illustrates several sample paths and the asymptotic pdf.

Understanding the amount of temporal correlation associated with environmental variability is important as it affects the entire population. 175 We make the assumption that environmental variability in the transmission rates has the same impact on NS and SS, relative to their mean values. That is, the coefficient of variation, CV , and the rate of return to equilibrium, r i = r, are the same for NS and SS. The two SDEs in (3) differ in their mean valuesβ i , i = 1, 2. In summary, the two SDEs in equation (3) coupled with the NHP and SDE models discussed in the previous two sections account for demographic and environmental variability in MERS and Ebola epidemics.

Numerical methods for simulation of the NHP and SDE models are discussed in Appendix D.

185 4. Results

We tested the functionality of the SDE and NHP models and found that both models generated pre- simulated. The SDE loses prediction accuracy for outbreaks when initial infection numbers are close to zero (when differences between discrete and continuous random variables are most evident), a common feature of diffusion approximations in general, e.g., (Ovaskainen and Meerson, 2010) . Better agreement between these two models is reached as the total number of infected individuals increases (Figure 3(b) ). An example of the 

We next utilized our model to explore the impact of environmental variability on MERS and Ebola outbreak probability and severity, given outbreaks that begin with a single infected SS, NS, or one of each infected individual. Given the low initial conditions used, we utilized our NHP model based on the observations that the SDE model was less accurate at low initial conditions, as shown in Figure 3 We were able to tune the level of environmental variability using r and CV , where r impacts the speed at which the disease transmission rate reverts to the mean value and CV impacts the degree to which fluctuations in the transmission rate deviate from the mean ( Figure 2) . Therefore, over a fixed period of time, the environmental variability increases with r but the magnitude of this effect is controlled by CV . (Figure 6(a) ).

These demonstrate that spikes in the transmission rate above the mean level often trigger corresponding increases in the number of infections over time (Figure 6(b) ). We also explored the effect of initial NS and SS transmission rates on model predictions, either using the mean valueβ i , i = 1, 2, or a random value from the asymptotic gamma distribution at the first time step. We found that the trends observed were conserved across model predictions using these two rules for transmission rate initialization, validating our strategies were used for small values of r, though this was not of notable issue within the range of r values we explored.

We next sought to explore the effect of sweeping r and CV values in disease scenarios with high initial 255 numbers of infected individuals, all introduced to the population at the same time, in which the SDE model is known to give rise to reliable predictions ( Figure E.2(a) in Appendix E). In fact, we validated that the SDE model predictions agreed with those of the NHP with environmental variability ( Figure E .5 in Appendix E), with only minor differences that do not impact the trends observed. As with simulations conducted with low initial conditions, we observed that the probability of an outbreak increased with increasing r values, 260 especially with large CV values (Figure 7(a) ). Furthermore, we observed that while environmental variability (by increasing CV ) reduces the chance of both MERS and Ebola outbreaks, those that do occur experience higher severity in terms of number of deaths (Figure 7(b) ), time to outbreak (Figure 7(c) ), number of infections (Figure 7 (d)), and peak number of infections (Figure 7 (e)). Interestingly, however, we observed a non-monotonic relationship between r and measures of outbreak severity. This is in contrast to our previous 265 observation that decreasing the r value increases the outbreak severity (Figure 4 , 5). In fact, we found that a medium r value resulted in the slowest time to outbreak (Figure 7 (c)) and peak infection (Figure 7 (f)) for both MERS and Ebola. We also found that these trends were observed for conditions with equal numbers of NS and SS susceptible individuals.

Environmental changes over time, such as social behavior, temperature fluctuations, or medical quality, influence contact between individuals and successful disease transmission (Brauer and Castillo-Chavez, 2012) .

Our simulations reveal that deviations in disease transmission rates, both above and below the mean, drive significant changes in outbreak dynamics. This aligns with our previously reported sensitivity analysis predictions, which uncovered the transmission rate as particularly impactful on disease outcomes (Edholm 275 et al., 2018) . Importantly, environmental variability reduces the probability of outbreak occurrence since temporal decreases in the transmission rate reduces the chance that the outbreak threshold will be reached.

On the other hand, once an outbreak is triggered, variability in the disease transmission rate increases the severity of the disease dynamics by increasing the number of deaths and total infections, increasing the peak number of infections reached, and reducing the time to the outbreak and peak number of infections. This is a 280 result of the fact that any spikes in the transmission rate has compounding effects as the number of infectious individuals increases. These methods can be used to investigate the effects of environmental variability on transmission or other model parameter values in epidemic or ecological settings.

While this study focuses on modeling transmission dynamics rather than evaluating specific public health strategies, our findings highlight the importance of altering environmental factors to mitigate disease spread.

Masking and social distancing are two such factors that have been shown to affect the transmission rate (Eikenberry et al., 2020; . rescinded but then applied again, such as the lockdowns with respect to COVID-19 (Acuña-Zegarra et al., 2020; Della Rossa et al., 2020; Dropkin, 2020) , could be viewed as contributing to a high amplitude and low frequency of variability (small r and large CV ). In this case, our models predict that such interventions may 290 lead to more severe epidemics. On the other hand, with high frequency in variability (large r) compared to the time scale of the disease dynamics, the outcomes may be similar to those with no environmental variability.

Consequently, high variability in public health interventions may result in their reduced effectiveness. Hence, findings of this study indirectly support the potential for modifiable environmental factors, influenced by public health strategies, to control disease transmission. Further work is needed to fully evaluate intervention 295 programs and their impact on disease prevalence.

We also made the interesting observation that the severity of outbreaks does not vary monotonically as a function of the degree of environmental variability. Indeed, we found that for outbreaks that are started by a number of SS and NS individuals introduced at the same time, rather than a single SS or NS individual, there is an optimal level of environmental variability (determined by r and CV values) that results in a delayed 300 time to outbreak and time to peak infection. This suggests that modeling could play a role in predicting an optimized distribution for transmission rates that can be attained by tuning social measures in such a way as to reduce outbreak probability and severity while minimizing socioeconomic impacts of social distancing.

Our study analyzes the dynamics of two diseases, MERS and Ebola, with the same R 0 value. Qualitatively, the trends in outcome measures shown in Figures 3 to 5 appear similar as the initial conditions, r, 305 and CV parameters are varied. In Figure 4 to Figure 5 , the number of deaths, the number infected, and peak number infected are similar among the three cases of initial conditions but are higher when r is small and CV is large for both diseases. The probability of an outbreak was higher for both diseases when the initial conditions included both a superspreader and a non-superspreader. Further, the time to outbreak and time to peak infection were more sensitive to increases in r when CV was higher. While the results between 310 the two diseases appear similar despite differences in compartments, the quantitative or estimated outcomes differ between MERS and Ebola. In Figure 7 In using Ebola and MERS as a case study for the impact of environmental variability on disease transmission dynamics, our theoretical investigations lay the foundation for understanding the impact of randomly 320 varying transmission rates on outbreak probability and severity. However, the amount of variability observed in a particular epidemic (r and CV ) depends on the variability and heterogeneity in social contacts, as well as factors affecting successful transmission during the epidemic. Quantifying and identifying the factors that account for this variability present a challenging problem. Nevertheless, the modeling strategy and trends we observed, which are closely in agreement between Ebola and MERS despite differences in asymptomatic 325 disease states, are expected to apply more broadly to other communicable diseases, such as COVID-19.

Overall, our results suggest that for a set mean transmission rate, environmental variability can reduce the probability of an outbreak yet increase the severity of those outbreaks that occur. Thus, tighter control over environmental variables such as adherence to hospital disease transmission-reduction guidelines, social distancing in the general population, and other control measures may be warranted.

It is worth mentioning that, given the large number of random variables, time-nonhomogeneity, and nonlinearity, exact analytical results for our models are not feasible. Certainly, in the case of time-homogeneous stochastic processes such as continuous-time Markov chain and SDE models, methods to approximate the mean and higher-order moments of the random variables, the final size, the duration until disease extinction, and the variability about the endemic state have been applied, e.g., Britton (2010) (2010); Van Kampen (1992) . Unfortunately, many of these methods do not extend to the time-nonhomogeneous models with mean-reverting processes. However, in simple SIR and SIS models, various outcomes from the timenonhomogeneous processes can be more easily compared to the analytical approximations from the corresponding time-homogeneous processes. An interesting topic for future research is to apply our methods and 340 the preceding methods in simpler settings to investigate and compare the effects of demographic and/or environmental variability and the amount of temporal correlation with mean-reverting processes on disease outcomes. Another future extension of our models could seek to include the contributions of novel pathogenic variants as well as the behavioral aspects of SS individuals, which may make them more susceptible to infection by NS and SS individuals. 

The underlying ODE model for MERS is

J o u r n a l P r e -p r o o f

The MERS and Ebola SDE models with vector random variables X = (X 1 , . . . , X n ) T , drift vector F (X, t) = F = (F 1 , . . . , F n ) T and diffusion matrix G(X, t) = G = [G ij ], i = 1, . . . , n and j = 1, . . . , m (n =number of random variables, m =number of events) are modified slightly to form new SDE models to ensure (1) the rates in the SDE are nonnegative and (2) the SDE sample paths are nonnegative (Allen et al., 505 2020; Cresson and Sonner, 2018) .

To ensure that the SDE dynamics follow closely those of the NHP, the transition rates defined in Tables 2 and C.1 (Appendix C) are assumed to be nonnegative (Allen et al., 2020) . In particular, the transition rates a(t n ) at time t n are replaced byâ

With this change the SDE models have the form dX =F (X, t) dt +G(X, t) dW .

Next to ensure that the SDE solutions are nonnegative some of the diffusion terms are modified slightly when a random variable is close to zero. For the SDE to have nonnegative sample paths it must be the case that each termG ij in the diffusion matrix satisfyG ij Xi=0 = 0 (Cresson and Sonner, 2018) . If this condition is not satisfied for a fixed i and j, then the SDE termG ij is modified as follows:

for 0 < ≤ 1 (Allen et al., 2020) . A value = 0.1 was chosen after testing the numerical code to ensure that reducing the value of does not change the results. Applying the previous modifications, the SDE models applied to MERS and Ebola have nonnegative transition rates and nonnegative sample paths (Allen et al., 2020) :

dX =F (X, t) dt +Ĝ(X, t) dW, X(0) ≥ 0.

Unlike MERS, the asymptotic stage for Ebola is short (Chowell and Nishiura, 2014) and thus it is omitted.

The underlying Ebola model is a system of 8 ODEs, four equations for SS and four for NS:

J o u r n a l P r e -p r o o f for i = 1, 2. If β i (t) ≡ β i is constant, then the basic reproduction number for the Ebola ODE model is

The corresponding Itô SDE SEIR model for Ebola has 8 Wiener processes, corresponding to the 8 events in Table C. 1. 

Rate, a

The Ebola SDE model with demographic variability can be derived in a similar manner as for the MERS SDE model:

(C.1) for i = 1, 2. For simplicity, the dependence on t has been omitted from the state variables, e.g.,

The Ebola SDE model is modfied slightly, similar to the MERS SDE model, to ensure that the rates in 

For the MERS NHP model, the sum of all of the rates in the right-hand column of Table 2 is

J o u r n a l P r e -p r o o f

In general, for well-defined and integrable time-varying parameters, given an event at time t n , the time until the next event τ can be found by inverting the following cumulative distribution function:

where the transmission rates β i (t) are functions of t but each of the random variables in Γ(t) are fixed at the time of the last event t n , e.g., S i (t n ), E i (t n ), etc. (Ross, 2014) . In particular, the time until the next event is

where U is a uniform random variable on [0, 1] and X = (S 1 , E 1 , A 1 , I 1 , R 1 , S 2 , E 2 , A 2 , I 2 , R 2 ) T is the discrete random vector for the stochastic process. For constant β i , this method simplifies to the Gillespie algorithm (Gillespie, 1977) . Calculation via inversion of equation (D.1) is time-consuming and complex. Also, as β i (t) is a random variable, this method is not possible for our models. Instead we apply a Monte Carlo method by incrementing time by a sufficiently small positive time step ∆t to ensure that each 525 of the 12 probabilities and Γ(t)∆t are nonnegative and less than one and that only one event, if any, occurs in a given time interval ∆t. A similar method is applied to the Ebola NHP model.

It is straightforward to include the two SDE equations for β i (t) in equation (3) in the NHP and SDE models for MERS and Ebola. The two SDES in equation (3) are approximated at a discrete set of points t = 0, ∆t, 2∆t, . . .,

where ∆t is sufficiently small and η is a random value from a standard normal distribution. Each of the 10 4 sample paths are approximated at the discrete set of points.

For the NHP models, the time step is chosen sufficiently small to ensure the probabilities and their sum 530 in Tables 2 and C.1 (Appendix C) are nonnegative in the Monte Carlo simulation. For the SDE models, the Euler-Maruyama method is applied. The random vector for 10 4 sample paths is updated each time step.

The numerical methods are checked for accuracy by reducing the time steps sequentially by 1/2 to ensure there are no significant differences in the output measures with smaller time steps. The time step chosen for the numerical simulation of the NHP models was ∆t = 5 × 10 −4 and for the SDE models was ∆t = 5 × 10 −2 . 

Global transport networks and infectious disease spread

Typhoid Mary: Captive to the Public's Health

Super-spreaders in infectious diseases

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

Ebola superspreading, The Lancet Infectious Diseases

Searching for superspreaders: Identifying epidemic patterns associated with superspreading events in 20

Airborne infection, Theoretical limits of protection achievable by building ventilation

Combining environmental assessment and contact investigations to make tuberculosis screening decisions

Dynamics and control of infections transmitted from person to person through the environment

Environmental factors influencing the spread of communicable diseases

Environmental determinants of infectious disease: a framework for tracking causal links and guiding public health research

Middle East respiratory syndrome coronavirus: another zoonotic betacoronavirus causing SARS-like disease

Variability in intrahousehold transmission of Ebola virus, and estimation of the household secondary attack rate

COVID-19) epidemics, the newest and biggest global health threats: what lessons have we learned?

Ebola virus disease

Response of a deterministic epidemiological system to a stochastically varying environment

Incorporating environmental stochasticity within a biological population model

Environmental variability and mean-reverting processes

Simulation and Inference for Stochastic Differential Equations: with R Examples

Stochastic models of population extinction

To mask or not to mask: Modeling the potential for face mask use by the general public to curtail the 465 COVID-19 pandemic

Effectiveness of facemasks to reduce exposure hazards for airborne infections among general populations

Modeling behavioral change and COVID-19 containment in Mexico: A trade-off between lockdown and compliance

A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic

COVID-19 UK lockdown forecasts and R0

Stochastic epidemic models: a survey

Comparison of Markov chain and stochastic differential equation population models under higher-order moment closure approximations

Estimating variability in models for recurrent epidemics: assessing the use of moment closure techniques

Stochastic Processes in Physics and Chemistry

A note on a derivation method for SDE models: Applications in biology and viability 485 criteria

Introduction to Probability Models

The authors acknowledge the support of an American Institute of Mathematics SQuaREs grant. A.P. was supported by NSF grant DMS-1815750. N.S. was supported by a Natural Sciences and Engineering Research

Council of Canada (NSERC) Postdoctoral Fellowship. In addition, we thank Texas Tech University for hosting our third WAMB group meeting and the P.W. Horn Professorship of L.J.S.A. for providing financial support. We thank the two referees for their helpful suggestions.