key: cord-0839386-6i5u3mgq
authors: Nipa, Kaniz Fatema; Allen, Linda J. S.
title: Disease Emergence in Multi-Patch Stochastic Epidemic Models with Demographic and Seasonal Variability
date: 2020-11-24
journal: Bull Math Biol
DOI: 10.1007/s11538-020-00831-x
sha: e237257a20bb47f08f402af5530142f4effa8874
doc_id: 839386
cord_uid: 6i5u3mgq

Factors such as seasonality and spatial connectivity affect the spread of an infectious disease. Accounting for these factors in infectious disease models provides useful information on the times and locations of greatest risk for disease outbreaks. In this investigation, stochastic multi-patch epidemic models are formulated with seasonal and demographic variability. The stochastic models are used to investigate the probability of a disease outbreak when infected individuals are introduced into one or more of the patches. Seasonal variation is included through periodic transmission and dispersal rates. Multi-type branching process approximation and application of the backward Kolmogorov differential equation lead to an estimate for the probability of a disease outbreak. This estimate is also periodic and depends on the time, the location, and the number of initial infected individuals introduced into the patch system as well as the magnitude of the transmission and dispersal rates and the connectivity between patches. Examples are given for seasonal transmission and dispersal in two and three patches.

Seasonal changes also affect the movement of susceptible and infected hosts and impact the spatial spread of disease in humans and in domestic and wild animal populations (Altizer et al. 2006) . Some of the many infectious diseases that are impacted by seasonal variations include seasonal influenza, avian influenza (AI), malaria, cholera, dengue, plague, foot-and-mouth disease, Ebola, Middle East respiratory syndrome (MERS), and severe acute respiratory syndrome (SARS), e.g., (Breban et al. 2009; Brown et al. 2014; Camitz and Liljeros 2006; Endo and Nishiura 2018; Fisman 2007; Gao et al. 2019; Parham and Michael 2010; Keeling 2005; McLennan-Smith and Mercer 2014; Schmidt et al. 2017; Vaidya and Wahl 2015) . The emergence of the new coronavirus disease named COVID-19 and the worldwide pandemic has raised concerns about seasonal or sporadic outbreaks after the pandemic period (Kissler et al. 2020) . It is the goal of the present investigation to study infectious disease outbreaks in a multi-patch setting when transmission within and movement between the patches are seasonal. Our modeling framework and methods with knowledge about travel patterns and seasonal trends will provide insight about effective methods for control or prevention of disease outbreaks.

Many modeling studies have focused on disease spread among discrete patches with no seasonal variations. These studies have been based on ordinary differential equations (ODEs), e.g., Arino 2009; Arino and van den Driessche 2003b, a, 2006; Baguette et al. 2014; Cosner et al. 2009; Gao and Ruan 2012; Gao et al. 2019; Kelly Jr. et al. 2016 ) and on stochastic models, e.g., (Ball 1991; Ball and Clancy 1993; Clancy 1994; Breban et al. 2009; Lahodny Jr and Allen 2013; McCormack and Allen 2007; Milliken 2017; Neal 2012) . The ODE multi-patch epidemic models have established upper and lower bounds on the basic reproduction number, existence and stability of an endemic equilibrium and conditions for disease persistence or extinction under various dispersal/movement assumptions. In a cholera ODE multi-patch model, Kelly Jr. et al. (2016) showed that targeting vaccination to specific patches is effective in controlling the disease. In stochastic time-homogeneous models, the effect of dispersal on epidemic final size and on probability of an outbreak has also been studied in multipatch settings.

The effects of seasonal variation on disease dynamics have been studied in ODE models, e.g., (Bacaër and Guernaoui 2006; Bacaër 2007; Gao et al. 2014; Jin and Lin 2018; McLennan-Smith and Mercer 2014; Mitchell and Kribs 2017; Posny and Wang 2014; Schwartz and Smith 1983; Suparit et al. 2018; Vaidya and Wahl 2015; Wang et al. 2012; Wang and Zhao 2017a, b; Zhang and Zhao 2007; Wolf et al. 2006 ) and in stochastic models, e.g., (Bacaër and Ait Dads 2014; Billings and Forgoston 2018; Breban et al. 2009; Jin and Lin 2018; Keeling et al. 2001; Lin et al. 2015; Nipa and Allen 2020) . The ODE models are nonautonomous and include seasonal variation through periodic coefficients. In these models, the existence of a basic reproduction number R 0 with stability of the disease-free state has been verified for R 0 < 1 (Bacaër 2007; Zhang and Zhao 2007) . Interestingly, numerical examples have demonstrated that R 0 for the nonautonomous model differs from the autonomous model when coefficients are fixed at their average values. Recently, Bacaër and Ait Dads (2014) showed that the basic reproduction number R 0 from the underlying nonautonomous ODE model serves as a threshold for disease extinction in the corresponding stochastic time-nonhomogeneous model. If R 0 < 1, the probability of disease extinction is equal to 1, but when R 0 > 1, the probability of disease extinction is periodic with probability < 1. Nipa and Allen (2020) applied these results to estimate the probability of disease extinction in a stochastic vector-host epidemic model. This investigation extends the previous research by combining multi-patch epidemic models with seasonal variation in a stochastic setting. In particular, we formulate a stochastic multi-patch epidemic model with seasonal variability in transmission and dispersal rates. Through multi-type branching process approximation and application of the backward Kolmogorov differential equation, an estimate for the probability of disease extinction is derived. It is shown that the probability of disease extinction and disease outbreak is periodic and less than 1 when R 0 > 1. In particular, the probability of disease extinction and disease outbreak is dependent on the time, the location, and the number of initial infected individuals i introduced at time τ into the patch system:

. . , i n ) and i j is the number of infected individuals in patch j for j = 1, . . . , n at the initial time τ ∈ [0, T] with T the total length of the seasons (generally one year) and n the total number of patches.

In the next section, a nonautonomous ODE multi-patch epidemic model is introduced which serves as a framework for the stochastic time-nonhomogeneous model in Sect. 3. Approximation of the time-nonhomogeneous process by a multi-type branching process and the backward Kolmogorov differential equation leads to estimates of the periodic probability of disease extinction. Several numerical examples with two patches in Sect. 4 illustrate the dependence of the probability of an outbreak on time, the location, and the initial number of infected individuals. Examples with three patches in Sect. 5 show the importance of patch arrangement and connectivity on the probability of a disease outbreak. In the concluding section, we summarize our results and suggest areas for further research.

We consider a susceptible-infected-recovered (SIR) multi-patch epidemic model with the dispersal of individuals between patches. The term dispersal is used to represent seasonal movement or migration rather than human movement associated with daily commuter travel, e.g., (Arino and van den Driessche 2003b; Cosner et al. 2009 ). Before defining the stochastic process, we define the underlying nonautonomous ODE model.

Let n be the total number of patches. For patch j, susceptible individuals are denoted as S j , infected individuals as I j , and recovered individuals as R j , j = 1, 2, . . . , n. We assume there are births and natural deaths as well as disease-related deaths. The total population size in each patch equals N j = S j + I j + R j , j = 1, 2, . . . , n and the total population size equals N = n j=1 N j . For each patch j, parameter μ j represents the birth and natural death rate, parameter α j is the disease-related death rate, and parameter γ j is the recovery rate. The transmission and dispersal rates are assumed to be time-periodic functions. The disease transmission rate is frequency-dependent, β j (t), and the dispersal rate from patch j to patch k is d jk (t), j, k = 1, 2, . . . , n and j = k. That is,

where = S, I , R represents the dispersal rate for either susceptible, infected or recovered individuals. These functions are continuous for t ∈ (−∞, ∞) with period T > 0.

The multi-patch model with a total of n patches takes the following form:

If there are no disease-related deaths, α j = 0, then the total population size N is constant. However, with dispersal and disease, the population size within each patch may change over time.

We make some simplifying assumptions to ensure there exists a unique disease-free equilibrium (DFE). Assume

for = S, I , R. In the absence of infection, I j ≡ 0, the equilibrium population sizes within each of the patches are constant, i.e., N 1 = N 2 = · · · = N n . The unique DFE isS j = N j = N /n, j = 1, . . . , n. The simplifying assumption in Eq.

(3) implies movement to and from a patch is balanced or symmetric. This assumption does not hold, for example, if there is directed movement such as asymmetric connectivity as in Chen et al. (2020) . We relax this assumption for an example in Sect. 4. Linearizing the system (2) about the DFE for the n infected states I = (I 1 , . . . , I n ) T results in the linear systemİ

The basic reproduction number, R 0 , depends on the two matrices, F(t) and V (t). Matrix F(t) represents new infections, a diagonal matrix, and matrix V (t) represents other transitions (dispersal, recovery and deaths):

and

where

. For the nonautonomous system, explicit solutions for the basic reproduction number can be calculated only in special cases, e.g., Mitchell and Kribs (2017) ; Wang and Zhao (2008) . For example, if there is no dispersal between patches, the basic reproduction number for patch j can be computed as follows:

We refer to R 0 j as the patch reproduction number which in this case is simply the average of the transmission rate, denoted asβ j , multiplied by the average duration of the infection, 1/(μ j + α j + γ j ).

In the more general case, Assumptions (A1)-(A7) (Appendix A) ensure that the basic reproduction number exists for the system (2) (Wang and Zhao 2008) . Theoretically, R 0 can be found from the following linear periodic system:

The spectral radius of the fundamental matrix solution of system, evaluated at t = T, equals one, then λ = R 0 (Theorem 2.1 p. 704, (Wang and Zhao 2008) ), i.e.,

If an analytical solution cannot be obtained, numerical methods as well as computer algebra systems can be used to compute the spectral radius of this matrix (Klausmeier 2008; Posny and Wang 2014) . In the following examples, we compute the results numerically by choosing a sequence of values λ i , i = 1, 2, . . . and 0 < 1, where 

Healthy birth in patch j S j → S j + 1 μ j n j Δt + o(Δt)

It should be noted that the basic reproduction number does not depend on the dispersal of S or R, only the dispersal rate of I . We use the terminology that patch j is high-risk if R 0 j > 1 and is low-risk if R 0 j < 1 ).

The ODE model serves as a framework for formulation of the time-nonhomogeneous stochastic process. For simplicity, the same notation is used for the variables as in the ODE model. The random variables are discrete and time is continuous, t ∈ [0, ∞),

The changes in the random variables are ±1 during a small interval of time Δt. The rates in the ODE model are used to define the infinitesimal transition probabilities. They are summarized in Table 1 . For example, a new infection in patch j at time t results in a change in the random variables, ΔS j (t) = S j (t +Δt)− S j (t) = −1 and ΔI j (t) = I j (t +Δt)− I j (t) = +1 with probability , I 1 (0), R 1 (0)) = (1999, 1, 0) (S 2 (0), I 1 (0), R 2 (0)) = (2000, 0, 0). Patch 1 is high-risk with R 01 = 3 and patch 2 is low-risk with R 0 = 0.2. One sample path illustrates a major outbreak, while the second sample path illustrates rapid disease extinction. The top two panels show the dynamics of susceptible, infected, and recovered individuals in the two patches. The bottom two panels show a close-up view of the two sample paths of the infected individuals. Parameter values are β 1 (t) = 18(1 + 0.8 sin(π t/2)), β 2 (t) = 1.2(1 + 0.8 sin(π t/2)), d jk = 2(1 + 0.8 sin(π t/2)), = S, I , R, μ 1 = μ 2 = 0.01, γ 1 = γ 2 = 5.99, and

The lower case letters s j , i j , and n j denote the values of the random variables S j (t), I j (t), and N j (t), respectively. A simple example of two patches illustrates the dynamics of the stochastic timenonhomogeneous process when either a disease outbreak or disease extinction occurs. In Fig. 1 , two sample paths are graphed for (I 1 (0), I 2 (0)) = (1, 0). The top two panels show the dynamics of the SIR process in two patches when an outbreak occurs. In this example, susceptible, infected, and recovered individuals have the same periodic dispersal rates in both patches, but the periodic transmission rates differ in patches 1 and 2. Disease extinction is not visible in the top two panels, but a close-up view in the bottom two panels shows the infected population in each patch. In one sample path, there is no outbreak, the number of infected individuals reaches zero by time 0.25 (approximately 3 weeks), whereas in the other sample path there is a disease 

outbreak. It is clear that a major outbreak occurs in patch 1, where the number of infected individuals reaches a level of over 500 individuals, and a smaller outbreak is visible in patch 2. There is only a single outbreak as the susceptible populations are depleted in both patches, and the infected population sizes reach zero in both patches before a second outbreak occurs. That is, the populations in both patches become disease-free.

A multi-type branching process approximation is applied to the n infected states I j in the time-nonhomogeneous process near the DFE, where S j = N j . The infinitesimal transition rates in the branching process approximation during a small period of time Δt are summarized in Table 2 .

The changes in the infected states are related to the transitions that occur in the linearized ODE system, Eq. (4). Unlike the nonhomogeneous stochastic process defined in Table 1 , the random variables for the branching process are not bounded above,

Bacaër and Ait Dads (2014) applied the forward Kolmogorov differential equation to show that the basic reproduction number R 0 from the ODE model serves as a threshold for the branching process approximation. In particular, they showed if R 0 ≤ 1, then the multi-type branching process with the same periodic rates has an asymptotic periodic probability of disease extinction that depends on the time of introduction τ and the number of infected individuals i = (i 1 , . . . , i n ), i.e., I j (τ ) = i j ≥ 0 and I (τ ) = i. If R 0 ≤ 1, then the asymptotic probability of disease extinction is 1, but if R 0 > 1, then In addition, the asymptotic probability of a disease outbreak can be approximated as follows:

We describe another method to obtain the estimate for the probability of disease extinction via probability generating functions (pgfs) and the backward Kolmogorov differential equation (Nipa and Allen 2020) . Denote the transition probability for the multi-type branching process as

where the notation m = (m 1 , . . . , m n ), I j (t) = m j ≥ 0, j = 1, . . . , n is the number infected at time t and i = (i 1 , . . . , i n ), I j (τ ) = i j , j = 1, . . . , n is the number infected at time τ . The backward Kolmogorov differential equation can be derived from the transitions defined in Table 2 , e.g., (Athreya and Ney 2004; Harris 1963; Jagers and Nerman 1985) :

where e j is the standard unit vector in R n . Next, we define the offspring pgfs, one for each of the random variables I j . The offspring pgfs for this time-continuous process are distinct from those in a discretetime Markov chain (DTMC). The number of infections generated by I j in the timecontinuous process is assumed to occur in a small period of time Δt rather than during the entire infectious period as in a DTMC. Given I j (τ ) = 1 and I k (τ ) = 0 for k = j, one of four events defined in Table 2 , if any, occurs in a small period of time Δt. An infected individual in patch j (1) transmits the infection to another individual in patch j or (2) dies (naturally or from the infection) or (3) recovers from the infection or (4) moves to another patch k. The pgf for I j that includes these four events is defined by

where u = (u 1 , . . . , u n ) ∈ [0, 1] n , j = 1, . . . , n. Each pgf has the property that f j (u, τ )| u=1 = 1 and f j is nonlinear in u and positive on the set [0, 1] n . The difference between this pgf and the pgf for the model with constant parameters (equation (12) in Lahodny Jr and Allen (2013)) is that the pgf (10) is dependent on τ .

In addition to the offspring pgfs for each of the infected states, we define the generating function for the multi-type branching process consisting of all n infected states. Let (11) where u = (u 1 , . . . , u n ). In the branching process, the random variables I j are assumed to be independent, resulting in

This is a reasonable assumption if the initial number of infected individuals is small. Let I j (τ ) = i j and the remaining infected states I k (τ ) = 0, k = j. Then, for i = (0, . . . , i j , . . . , 0), differentiating Eq. (12) with respect to τ , using the definition (11), applying the identity (9), and substituting the pgf (10), leads to a differential equation for G e j . That is,

where

(A more detailed derivation of (13) in the case of two patches is given in Appendix B.) An estimate of the probability of extinction follows from (11) and (12):

The asymptotic probability of disease extinction is a periodic function

for τ ∈ [0, T] (Bacaër and Ait Dads 2014). In particular, the probability of an outbreak given I j (τ ) = i j , j = 1, . . . , n equals

In the following examples for two and three patches, we estimate the probability of an outbreak given by Eq. (14) from the differential equation (13). In particular, we make a change of variable s = t − τ in the differential equation (13) and for simplicity, we fix t = 0, P j (s) = p e j ,0 (τ, 0), so that

The preceding differential equation is solved (backward in time) for sufficiently large s > 0 so that there is convergence to a periodic solution,

For each numerical example, we estimate P ext (e j , τ ) and apply formula (14). This estimate is checked against the full nonlinear time-nonhomogeneous stochastic process based on the transition probabilities in Table 1 . A Monte Carlo approach is used with a small time step Δt, chosen sufficiently small, such that during each time step, only one, if any, of the events in Table 1 occurs (Σ(t)Δt < 1). The accuracy of the Monte Carlo approximation is tested by using a sequence of progressively smaller time steps. A total of 10 4 sample paths are simulated for a given set of initial conditions at a specific time τ . Each sample path continues until a time t s > τ when either j I j (t s ) = 0 or j I j (t s ) = 50. If the total infected population reaches 50, then it is assumed that there is an outbreak. The proportion of sample paths q out of 10 4 that reach zero is an estimate for the probability of extinction, and the remaining proportion 1 − q that reach 50 is an estimate of the probability of an outbreak.

For two patches, the SIR ODE model has a simple form, illustrated in the compartmental diagram in Fig. 2 .

Sinusoidal functions are assumed for the transmission and dispersal rates with period T:

j, k = 1, 2, j = k. The amplitude of these rates are |δ j | and |δ jk |, respectively. The dispersal rates in each patch are the same, except for the value of δ jk which may be either positive or negative while the transmission rates differ in the magnitude ofβ j . If the values δ j and δ jk have the same sign, we say the parameters for transmission and dispersal are "synchonized," but if they have opposite signs, one positive and one negative, we say they are "desynchronized." Synchronization of transmission and dispersal rates means that the times of the extrema for the two rates are the same, whereas desynchronization means the times of the extrema are reversed, i.e., high transmission rate corresponds to low dispersal rate and vice versa. Figure 3 shows an example when the two patches are synchronized or desynchronized with respect to dispersal and transmission. Parameter values for the two-patch model are given in Table 3 . The year is divided into four seasons, and the transmission and dispersal rates in Eq. (15) have a period T = 4. The time unit 1 equals three months. We assume either all susceptible, infected, and recovered individuals disperse at the same rate or only infected individuals disperse. In addition, we let μ i = 0 = α i or μ i = 0.01 and α i = 0.5 but keep R 0i , i = 1, 2 constant, i.e., R 0i =β i /(γ i + μ i + α i ). Each patch has the same recovery and death rates but differs in the average transmission rates,β i , i = 1, 2. Patch 1 is high-risk with patch reproduction number R 01 = 3 and patch 2 is low-risk with R 02 = 0.2 (Fig. 4) .

The following examples show that the probabilities of an outbreak depend on the time and location of the first infected individual and the relation between dispersal and transmission rates.

In Fig. 5 , four different cases are considered, (A)-(D), with constant transmission or dispersal and synchronized or desynchronized transmission and dispersal. In the top two panels of each of (A)-(D), the transmission and dispersal rates are graphed for patches 1 and 2 and in the bottom two panels, the periodic probabilities of an outbreak for initial conditions (I 1 (τ ), I 2 (τ )) = (1, 0) and (0, 1), respectively. The branching process estimate is verified by checking thirteen different times τ = 0, 1/3, 2/3, . . . , 4 Table 3 with a Monte Carlo simulation of 10 4 sample paths of the nonhomogeneous stochastic process. The results from the Monte Carlo simulations show good agreement with the branching process approximation. Dispersal between the two patches allows an outbreak to occur even when an infected individual is introduced into a low-risk patch, but there is a much smaller probability of an outbreak if the infection is initiated into a low-risk patch than into a high-risk patch (Fig. 5) . With desynchronized dispersal and transmission (Fig. 5D) , the probability of an outbreak when the infection starts in the low-risk patch is smaller than in cases (A), (B), or (C). Also, fluctuation in the transmission rate has a greater effect on the amplitude of the periodic outbreaks than dispersal (compare Fig. 5A, B) . In all cases, the time of the highest (lowest) transmission rate is closely related to the highest (lowest) probability of an outbreak. But the extrema of the probability of an outbreak are shifted slightly left of the extrema of the transmission rates. This phenomenon Table 4 Average probability of an outbreakP outbreak (1, 0) for patch 1 andP outbreak (0, 1) for patch 2 and the threshold value R 0 for the two-patch system when R 01 = 3 and R 02 = 0.2 was also noted by Nipa and Allen (2020) in a vector-host epidemic model. When infected individuals are first introduced, they may not transmit the infection to others immediately. Whether there is a successful transmission depends on the duration of the infectious period (average duration 1/(γ j + μ j + α j )), the patch j, the location of the infectious individual and the value of the transmission rate β j (t) at the time of transmission. The largest probability of an outbreak occurs just before the peak of transmission when the transmission rate is increasing. Table 4 records the average of the periodic probabilities of an outbreak that are displayed in Fig. 5A -D. The average is computed as follows:

For each case (A)-(D), the threshold value R 0 from the underlying nonautonomous ODE model is also calculated. For comparison purposes, the baseline case (Base) with constant transmission and dispersal rates, set at their average values, is also included. It is notable that the four cases of periodic transmission and dispersal have an overall lower average probability of an outbreak than at the baseline case. The lowest average probabilities of an outbreak occur with desynchronized transmission and dispersal. The magnitude of R 0 is not a good predictor of disease severity in these cases as the largest value of R 0 has the lowest average probability of an outbreak. We show that the branching process estimate is also a good approximation when there are several infected individuals introduced into one or both patches (Fig. 6) . For two patches, the probability of an outbreak in Eq. (14) equals

The probability of an outbreak is close to 1 when I 2 (τ ) = 3 and τ ∈ [0, 5/3] but drops rapidly to < 0.3 for τ ∈ [7/3, 3]. The example in Fig. 6 illustrates synchronous dispersal and transmission. Additional Monte Carlo simulations showed good agreement with the branching process estimate for desynchronous dispersal and transmission and for unequal population sizes in the two patches ranging from 200 to 5000 (not shown). The branching process approximation may fail if population sizes are too small. Some other limitations of the branching process approximation are discussed in the conclusion. Fig. 5C , the probability of an outbreak for one infected individual in patch 2 and either 0, 1, 2, or 3 infected individuals in patch 1 (four graphs ordered from bottom to top, have initial conditions for the high-risk patch 1, I 1 (τ ) = 0, 1, 2, 3, respectively ). The blue circles are the results of the Monte Carlo simulation of the nonhomogeneous stochastic process with 5 × 10 3 sample paths at τ = 0, 1/3, 2/3, 1, . . . , 4. The parameter values are given in Table 3 Fig. 7 (Color Figure Online) The probabilities of an outbreak for initial conditions (I 1 (τ ), I 2 (τ )) = (1, 0) (patch 1, solid curves) or (0, 1) (patch 2, dashed curves) are graphed when the average dispersal rate is increased d I jk = 1/4, 1/2, 1, 2, 4, 8, j, k = 1, 2, j = k, when transmission and dispersal rates in patches 1 and 2 are synchronized (Fig. 5C ). The arrows indicate the direction of change in the probability of an outbreak asd I jk increases

Next, we examine how an increase in the dispersal rate affects the probability of an outbreak (Fig. 7) . As the average dispersal rated I jk increases, the periodic curves for the probability of an outbreak in the two patches come closer together. That is, with a large amount of mixing of infected individuals in the two patches, whether infection is initiated into the high-risk or the low-risk patch does not matter, as the periodic probabilities of outbreaks for patches 1 and 2 are very close to each other. This phenomenon also occurs when constant dispersal rates are increased (Lahodny Jr and Allen 2013) .

In the final two-patch example, we relax the restriction imposed on dispersal, condition (3). We assume individuals disperse at different rates in the two patches, d 12 = d 21 , 

with remaining parameters

for i = 1, 2. In patch 1, transmission and dispersal are synchronized but in patch 2, they are desynchronized. The values of the patch reproduction numbers are the same as in Fig. 4 . Relaxation of the condition (3) results in oscillation of the populations in the absence of infection (Fig. 8) . In the Monte Carlo simulations of the timenonhomogeneous process, the initial conditions for the population sizes are equal to the disease-free population sizes. That is, N j (τ ) takes the value of S j (τ ) in Fig. 8 for τ ∈ [0, T]. When infection is introduced, (I 1 (τ ), I 2 (τ )) = (1, 0) or (0, 1), then S j (τ ) = N j (τ ) − I j (τ ), and R j (τ ) = 0, j = 1, 2. The branching process estimate in Fig. 9 agrees well with the Monte Carlo simulations of the time-nonhomogeneous process. The reason for this agreement is due to the assumption of frequency-dependent transmission; S i (t) ≈ N i (t) at the initiation of an epidemic.

For three patches, we consider how the connections between low-and high-risk patches affect the probability of an outbreak. For each example, we keep the same relation between transmission and dispersal rates, as in Fig. 5A -D. Three types of connection between the patches are referred to as full connection (FC), bidirectional connection (BC), and circular connection (CC).

With full connection (FC), each patch is connected to its neighbors by a one-step transition, i.e., 1 ↔ 2 ↔ 3 ↔ 1 (Fig. 10) . Patch 1 is high-risk with R 01 = 5, and (16) and (17) patches 2 and 3 are low-risk with R 02 = 0.5 and R 02 = 0.25. The second type of connection is called a circular connection (CC), i.e., 1 → 2 → 3 → 1. The patch reproduction numbers are the same as in FC. In the third type of connection, referred to as a bidirectional connection (BC), the patches are linearly ordered and connected only to nearest neighbor, i.e., 1 ↔ 2 ↔ 3. The arrangement is changed in the third case, the second patch has R 02 = 0.25 and R 03 = 0.5. The periodic functions for dispersal and transmission are given in Eq. (15) with parameter values in Table 3 .

The periodic probabilities of an outbreak with FC are illustrated in Fig. 11 when initial conditions (I 1 (τ ), I 2 (τ ), I 3 (τ )) = (1, 0, 0), (0, 1, 0), and (0, 0, 1) for patches 1, 2, and 3, respectively. (The graphs for BC and CC are in Appendix C, Figs. 13 and 14). Similar trends are observed for three patches with FC as in two patches. The largest probability of an outbreak occurs when infected individuals are introduced into the high-risk patch 1 and the lowest when introduced into low-risk patch 3 with the smallest patch reproduction number, R 03 = 0.25. Desynchronized dispersal and transmission rates in case (D) show that the peaks of an outbreak are smaller in both of the low-risk patches but larger in the high-risk patch than in cases (A), (B), or (C). The periodic transmission rate has a greater impact on the magnitude of the periodic probability of an outbreak than periodic dispersal, case (A) versus case (B).

There are similar trends in the average probabilities for cases (A)-(D) for FC as in two patches (Tables 4, 5). Case (D), desynchronized dispersal and transmission, has the largest threshold value of R 0 , but the average probabilities are generally lower than for constant transmission and dispersal.

The most important difference between FC and the two other connections is whether a low-risk patch is directly connected to the high-risk patch. The low-risk patch directly and (0, 0, 1) for patches 1, 2, and 3, respectively. Patch 1 is high-risk, R 01 = 5, and patches 2 and 3 are low-risk, R 02 = 0.5 and R 02 = 0.25. Parameter values are given in Table 3 connected to the high-risk patch (patch 3 in CC and patch 2 in BC) has a larger probability of an outbreak than the low-risk patch not directly connected. In the case of BC, infection initiated in the patch with the smallest patch reproduction number does not correspond to the lowest probability of an outbreak. More important is the direct connection to the high-risk patch (see Table 6 ). In a final example for three patches, we compare the probabilities of an outbreak in each of the three patches with synchronized dispersal and transmission, (case (C) in Figs. 11, 13, 14) for FC, BC, and CC. In Fig. 12 , the changes in the periodic probabilities of an outbreak are shown for increasing recovery rates (i) γ i = 0.5 (top), (ii) γ i = 1, (ii) γ i = 2, (iv) γ i = 4 (bottom), i = 1, 2, 3. Here, we assume natural and disease-related death rates are μ i = 0 = α i so that R 0i =β i /γ i . As the recovery rate γ i increases from (i) to (iv), the average duration of the infection shortens, each patch reproduction number R 0i decreases and the probability of an outbreak decreases. The low-risk patches 2 and 3 have the lowest probability of outbreak in all three cases with a large recovery rate. It is interesting to note that the high-risk patch in BC and CC has a larger probability of an outbreak than FC, but the low-risk patches in BC and CC have lower probabilities of outbreaks than FC. This is clearly visible in case (D) with the largest recovery rate. The lowest risk of an outbreak in Figure 12 occurs when infections are introduced during τ ∈ (2, 3) when transmission is decreasing. Fig. 10 . Initial conditions in patch 1 are (I 1 (τ ), I 2 (τ ), I 3 (τ )) = (1, 0, 0) (solid black curve), in patch 2, (I 1 (τ ), I 2 (τ ), I 3 (τ )) = (0, 1, 0) (red dashed curve) and in patch 3, (I 1 (τ ), I 2 (τ ), I 3 (τ )) = (0, 1, 0) (blue dotted curve). Values of the recovery rate in rows (i) γ i = 0.5, (ii) γ i = 1, (iii) γ i = 2 and (iv) γ i = 4, i = 1, 2, 3 with natural and disease-related death rates set to zero μ i = 0 = α i . All other parameters are in Table 3 Table 7 is a summary of the average probabilities of an outbreak and the values of R 0 for FC, BC, and CC corresponding to Fig. 12 . For comparison purposes, consider the case where the patches are isolated, not connected via dispersal, and there is no seasonal variation in transmission, β j (t) =β j . In this case, the well-known result of Whittle (1955) can be applied: the probability of an outbreak in the low-risk patches equals zero, and the probability of an outbreak in the high-risk patch equals where R 01 =β 1 /γ 1 . With isolation and no seasonality, the probabilities of an outbreak in the high-risk patch equal 0.95, 0.9, 0.8, and 0.6 in cases (i), (ii), (iii), and (iv), respectively. With dispersal connecting the high-and low-risk patches and seasonality, it can be seen that the average probabilities in Table 7 decrease in the high-risk patch, from 0.6 to 0.28 in case FC (iv), but the average probabilities of an outbreak initiated from the low-risk patches increase significantly from 0 to as high as 0.74 in case FC (i) in patch 2.

Seasonal changes that affect the pathogen life cycle or the within-host and betweenhost interactions result in seasonal outbreaks of many infectious diseases, such as seasonal influenza, AI, malaria, cholera, dengue, Ebola, and coronavirus diseases (Altizer et al. 2006; Endo and Nishiura 2018; Fisman 2007; Grassly and Fraser 2006; Kissler et al. 2020; Martinez 2018) . In this investigation, we applied a stochastic time-nonhomogeneous, multi-patch epidemic model with seasonal and demographic variability to investigate how time and location affect the probability of an outbreak when the infection is introduced into the patch system. We applied a multi-type branching process approximation to the infected classes to obtain an analytical estimate for the probability of a disease outbreak. Examples with two and three patches were used to test the branching process approximation against the full nonlinear timenonhomogeneous process. In particular, we investigated how the periodic probability of an outbreak changes with different assumptions regarding patch arrangement, connectivity, and relations between the periodic dispersal and transmission rates (cases (A)-(D) and FC, BC, and CC). We found that seasonality in dispersal and transmission affects the time and location of the greatest risk for an outbreak. If the infection is introduced into a high-risk patch (R 0 > 1) at a time when the transmission rate is large and increasing, then the probability of an outbreak is high. But if the infection is introduced at a time when the transmission rate is small and decreasing, or the dispersal rate to a low-risk patch is large, then the probability of an outbreak is low. Interestingly, seasonality in transmission and dispersal rates can, in some cases, result in a lower average risk of an outbreak than if the transmission and dispersal rates are constant (Tables 4, 5, 6). The shapes of the outbreak probabilities do not directly coincide with those of either the dispersal or transmission rates. In an ODE patch model without seasonality, developed for cholera , Kelly Jr. et al. (2016) investigated the optimal vaccination strategy. In a linear arrangement of patches, similar to BC, they found that after an outbreak occurs in a patch, the optimal vaccination strategy is to target neighboring or downstream patches. In our stochastic model, with dispersal and transmission rates changing during the season, the magnitude of these rates and whether they are increasing or decreasing at the time of the outbreak is important and makes the stochastic optimal control problem with seasonality an interesting and challenging problem.

As shown here and elsewhere, multi-type branching processes are often a good approximation method for the probability of an outbreak near the DFE, e.g., (Allen and van den Driessche 2013; Allen and Lahodny Jr 2012; Lahodny Jr and Allen 2013; Lahodny Jr et al. 2015; Milliken 2017; Nipa and Allen 2020) . The branching process approximation may fail if R 0 is too close to one, if demographic dynamics dominate the infection dynamics, or if population sizes are too small. The reasons for these failures are that the underlying assumptions of the multi-type branching process are no longer valid. The total population sizes must be sufficiently large and not deviate substantially from the disease-free state, and the initial infected population size must be small but grow to a sufficiently large size. Therefore, it is important to test the branching process approximation against the time-nonhomogeneous process.

There are many generalizations of these stochastic multi-patch epidemic models worthy of further investigation, such as additional stages or classes of infection, a mass action incidence rate, patch-dependent population densities, patch arrangements, and alternate movement and dispersal strategies. Realistic models that incorporate both demographic and seasonal variability in conjunction with case data and environmental data on specific diseases provide useful public health information on preventive and control measures such as isolation, quarantine, and travel restrictions. The impact of such measures has been evaluated in general models and in specific models applied to outbreaks of SARS, MERS, Ebola, and COVID-19, e.g., (Arino et al. 2007; Camitz and Liljeros 2006; Chinazzi et al. 2020; Dénes and Gumel 2019; Gao and Ruan 2012; Hsieh et al. 2007; Kwok et al. 2019; Parmet and Sinha 2020; Peak et al. 2020 ) and references therein. Application of our modeling framework and methods, coupled with information about travel patterns and seasonal trends on specific diseases (such as seasonal influenza, AI, dengue, malaria, cholera, Ebola, and coronavirus diseases), provide additional insight into the times and locations for quarantine and travel restrictions.

The following assumptions are taken from Wang and Zhao (2008) and applied to the ODE system (2). The system can be written as follows:

where x = (x 1 , . . . , x 3n ) = (I 1 , . . . , I n , S 1 , . . . , S n , R 1 , . . . , R n ) and

In particular, linearization of system (18) about the DFE has the forṁ

Matrices F(t) and V (t) are defined in Eqs. (5) and (6) and O is the 2n ×2n zero matrix. Letż = M(t)z andẏ = −V (t)y. The fundamental matrix solutions for these linear periodic systems are known as monodromy matrices. Let these monodromy matrices be denoted as φ M (t) and φ −V (t), respectively. The following assumptions guarantee there exists a basic reproduction number R 0 > 0 (Theorem 2.1 p. 704, Wang and Zhao (2008) ).

x) are nonnegative and continuous on R × R 3n + and continuously differential with respect to x.

x) = 0 for i = n + 1, . . . , 3n. (A5) If I = 0, then F i (t, x) = V + i (t, x) = 0 for i = 1, . . . , n. (A6) The DFE is stable if I = 0, i.e., the spectral radius ρ(φ M ( p)) < 1. (A7) The spectral radius ρ(φ −V ( p)) < 1.

Assumptions (A1)-(A7) hold for the ODE system (2). The following assumptions were applied by Bacaër and Ait Dads (2014) for the multitype branching process approximation based on the linear system (4).

where matrices F(t) and V (t) are defined in (5) and (6) (2014)).

(H1) F(t) is nonnegative with at least one entry strictly positive for t ≥ 0. From these assumptions, the following results hold for the multi-type branching process. If R 0 ≤ 1, then the multi-type branching process with the same periodic rates (defined in Table 2 ) has an asymptotic probability of extinction P ext (i 1 , . . . , i n , τ ) = 1 and if R 0 > 1, then the asymptotic probability of disease extinction is a p-periodic function satisfying 0 < P ext (i 1 , . . . , i n , τ ) < 1.

Applying a multi-type branching process approximation for two patches, the detailed steps of the derivation of the probability of extinction are derived here using generating functions and the backward Kolmogorov differential equation. For I j (τ ) = 1 and I k (τ ) = 0, j, k = 1, 2 and k = j the offspring pgf for patch j is f j (u 1 , u 2 , τ ) = β j (τ )u 2 j +μ j +α j +γ j +d I jk (τ )u k β j (τ )+μ j +α j +γ j +d I jk (τ ) , j = 1, 2.

The backward Kolmogorov differential equation in (9) for the branching process approximation of I 1 and I 2 is − ∂ p (i 1 ,i 2 ),( j,k) (τ, t) ∂τ = β 1 (τ )i 1 p (i 1 +1,i 2 ),( j,k) (τ, t) + (μ 1 + α 1 + γ 1 )i 1 p (i 1 −1,i 2 ),( j,k) (τ, t) + d I 21 (τ )i 2 p (i 1 +1,i 2 −1),( j,k) (τ, t) + β 2 (τ )i 2 p (i 1 ,i 2 +1),( j,k) (τ, t) + (μ 2 + α 2 + γ 2 )i 2 p (i 1 ,i 2 −1),( j,k) (τ, t) + d I 12 (τ )i 1 p (i 1 −1,i 2 +1),( j,k) (τ, t) − [β 1 (τ )i 1 + (μ 1 + α 1 + γ 1 )i 1 + d I 21 (τ )i 2 + β 2 (τ )i 2 + (μ 2 + α 2 + γ 2 )i 2 + d I 12 (τ )i 1 ] p (i 1 ,i 2 ),( j,k) (τ, t).

Also, the generating function G (i 1 ,i 2 ) (u 1 , u 2 , τ, t) for the two-patch model can be expressed in terms of the transition probabilities in (11):

G (i 1 ,i 2 ) (u 1 , u 2 , τ, t) = E u I 1 (t) 1 u I 2 (t) 2 |(I 1 (τ ), I 2 (τ )) = (i 1 , i 2 ) = j,k≥0 p (i 1 ,i 2 ),( j,k) (τ, t)u j 1 u k 2 , τ < t.

For j = k = 0, G (i 1 ,i 2 ) (0, 0, τ, t) = p (i 1 ,i 2 ),(0,0) (τ, t) is the probability of disease extinction. From the independent assumptions of the random variables I 1 and I 2 in the branching process approximation, it follows that G (i 1 ,i 2 ) (u 1 , u 2 , τ, t) = [G e 1 (u 1 , u 2 , τ, t)] i 1 [G e 2 (u 1 , u 2 , τ, t)] i 2 ,

where e 1 = (1, 0) and e 2 = (0, 1). For initial data (I 1 (τ ), I 2 (τ )) = (i 1 , i 2 ), the probability of extinction at time t > τ is p (i 1 ,i 2 ),(0,0) (τ, t) = [p e 1 ,(0,0) (τ, t)] i 1 [ p e 2 ,(0,0) (τ, t)] i 2 .

(23)

The periodic probabilities of an outbreak are computed from the branching process approximation and verified with the time-nonhomogeneous process (circles) for the three-patch models with bidirectional (BC) or circular connections (CC) (Figs. 13, 14) . Parameter values are given in Table 3 . . The initial conditions are (I 1 (τ ), I 2 (τ ), I 3 (τ )) = (1, 0, 0), (0, 1, 0) and (0, 0, 1) for patches 1, 2, and 3, respectively. Patch 1 is high-risk, R 01 = 5, and patches 2 and 3 are low-risk, R 02 = 0.5 and R 03 = 0.25

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

Asymptotic profiles of the steady states for an SIS epidemic patch model

Seasonality and the dynamics of infectious diseases

Modeling and dynamics of infectious diseases

The basic reproduction number in a multi-city compartmental epidemic model

A multi-city epidemic model

Disease spread in metapopulations

Quarantine in a multi-species epidemic model with spatial dynamics

Mineola Bacaër N (2007) Approximation of the basic reproduction number R0 for vector-borne diseases with a periodic vector population

On the probability of extinction in a periodic environment

The epidemic threshold of vector-borne diseases with seasonality: the case of cutaneous leishmaniasis in Chichaoua, Morocco

The pros and cons of applying the movement ecology paradigm for studying animal dispersal

Dynamic population epidemic models

The final size and severity of a generalised stochastic multitype epidemic model

Seasonal forcing in stochastic epidemiology models

The role of environmental transmission in recurrent avian influenza epidemics

Neutrality, cross-immunity and subtype dominance in avian influenza viruses

The effect of travel restrictions on the spread of a moderately contagious disease

Asymptotic profiles of the steady states for an SIS epidemic patch model with asymmetric connectivity matrix

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak

Some comparison results for multitype epidemic models

The effects of human movement on the persistence of vector-borne diseases

Modeling the impact of quarantine during an outbreak of Ebola virus disease

The role of migration in maintaining the transmission of avian influenza in waterfowl: A multisite multispecies transmission model along East Asian-Australian flyway

Seasonality of infectious diseases

A multipatch malaria model with logistic growth populations

A periodic Ross-Macdonald model in a patchy environment

Habitat fragmentation promotes malaria persistence

Seasonal infectious disease epidemiology

Impact of travel between patches for spatial spread of disease

Branching processes in periodically varying environment

Periodic solution of a stochastic sirs epidemic model with seasonal variation

Models of foot-and-mouth disease

Seasonally forced disease dynamics explored as switching between attractors

The impact of spatial arrangements on epidemic disease dynamics and intervention strategies

Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period

Floquet theory: a useful tool for understanding nonequilibrium dynamics

Epidemic models of contact tracing: systematic review of transmission studies of severe acute respiratory syndrome and Middle East respiratory syndrome

Probability of a disease outbreak in stochastic multipatch epidemic models

Estimating the probability of an extinction or major outbreak for an environmentally transmitted infectious disease

Nontrivial periodic solution of a stochastic epidemic model with seasonal variation

The calendar of epidemics: seasonal cycles of infectious diseases

Multi-patch deterministic and stochastic models for wildlife diseases

Complex behaviour in a dengue model with a seasonally varying vector population

The probability of extinction of infectious salmon anemia virus in one and two patches

A comparison of methods for calculating the basic reproductive number for periodic epidemic systems

The basic reproduction number and the probability of extinction for a dynamic epidemic model

Effects of demographic, environmental and seasonal variability on disease outbreaks in stochastic vector-host, multi-patch and dengue epidemic models

Modeling the effects of weather and climate change on malaria transmission

Covid-19-the law and limits of quarantine

Individual quarantine versus active monitoring of contacts for the mitigation of COVID-19: a modelling study

Computing the basic reproductive numbers for epidemiological models in nonhomogeneous environments

Spatiotemporal fluctuations and triggers of Ebola virus spillover

Infinite subharmonic bifurcation in an SEIR epidemic model

A mathematical model for Zika virus transmission dynamics with a time-dependent mosquito biting rate

Avian influenza dynamics under periodic environmental conditions

Threshold dynamics for compartmental epidemic models in periodic environments

Dynamics of a time-delayed Lyme disease model with seasonality

A malaria transmission model with temperature-dependent incubation period

A simple stochastic model with environmental transmission explains multi-year periodicity in outbreaks of avian flu

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

A multi-patch epidemic model with periodic demography, direct and indirect transmission and variable maturation rate

A periodic epidemic model in a patchy environment

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Acknowledgements This research is part of the PhD research of Nipa (2020) and is partially supported by

For simplicity, let G (i 1 ,i 2 ) ≡ G (i 1 ,i 2 ) (u 1 , u 2 , τ, t) and p (i 1 ,i 2 ),( j,k) ≡ p (i 1 ,i 2 ), ( j,k) (τ, t) . By differentiating G (i 1 ,i 2 ) with respect to τ in Eq. (22) and applying the identitiesSubstituting the backward Kolmogorov differential equation (20) into the right side, applying the identity (22), using offspring pgf f j , and simplifying leads toandSince the right side of the differential equation is independent of u 1 and u 2 , the differential equation also applies to G e j (0, 0, τ, t) = p e j ,(0,0) (τ, t) = p e j ,(0,0) . That is, dp e 1 ,(0,0) dτ = (β 1 (τ ) + μ 1 + α 1 + γ 1 + d I 12 (τ ))[ f 1 ( p e 1 ,(0,0) , p e 2 ,(0,0) , τ ) − p e 1 ,(0,0) ], dp e 2 ,(0,0) dτ = (β 2 (τ ) + μ 2 + α 2 + γ 2 + d I 21 (τ ))[ f 2 ( p e 1 ,(0,0) , p e 2 ,(0,0) , τ ) − p e 2 ,(0,0) ].