Disease Emergence in Multi-Patch Stochastic Epidemic Models with Demographic and Seasonal Variability.

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.

Introduction

Seasonal fluctuations in disease incidence have been known since the time of Hippocrates, about 400 B.C. (Fisman 2007). Seasonal changes in temperature, humidity, and rainfall affect pathogen survival, host behavior, and host immune function, which in turn affect successful transmission (Grassly and Fraser 2006; Altizer et al. 2006). 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., (Allen et al. 2007; 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 multi-patch 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 with stability of the disease-free state has been verified for (Bacaër 2007; Zhang and Zhao 2007). Interestingly, numerical examples have demonstrated that 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 from the underlying nonautonomous ODE model serves as a threshold for disease extinction in the corresponding stochastic time-nonhomogeneous model. If , the probability of disease extinction is equal to 1, but when , the probability of disease extinction is periodic with probability 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 . 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:and , where and is the number of infected individuals in patch j for at the initial time with 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.

ODE Multi-Patch Model

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 , infected individuals as , and recovered individuals as , . We assume there are births and natural deaths as well as disease-related deaths. The total population size in each patch equals and the total population size equals . For each patch j, parameter represents the birth and natural death rate, parameter is the disease-related death rate, and parameter is the recovery rate. The transmission and dispersal rates are assumed to be time-periodic functions. The disease transmission rate is frequency-dependent, , and the dispersal rate from patch j to patch k is , and . That is,where represents the dispersal rate for either susceptible, infected or recovered individuals. These functions are continuous for with period .

The multi-patch model with a total of n patches takes the following form:where , , , and with all parameters nonnegative such that and for . In general,If there are no disease-related deaths, , 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). Assumefor . In the absence of infection, , the equilibrium population sizes within each of the patches are constant, i.e., . The unique DFE is , . 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 results in the linear systemThe basic reproduction number, , 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):andwhere . 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 as the patch reproduction number which in this case is simply the average of the transmission rate, denoted as , multiplied by the average duration of the infection,

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, can be found from the following linear periodic system:The spectral radius of the fundamental matrix solution of system, evaluated at , equals one, then (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 , and , where until a value of is reached such that . Then, . 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 and is low-risk if (Allen et al. 2007).

Time-Nonhomogeneous Stochastic Process

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, ,for . The changes in the random variables are during a small interval of time . The rates in the ODE model are used to define the infinitesimal transition probabilities. They are summarized in Table 1.

Table 1

Infinitesimal transition probabilities for the stochastic n-patch model

Description	Transition	Probabilities
Infection in patch j	S_j→S_j-1,I_j→I_j+1	β_j(t)i_jn_js_jΔt+o(Δt)
Healthy birth in patch j	S_j→S_j+1	μ_jn_jΔt+o(Δt)
Healthy death in patch j	S_j→S_j-1	μ_js_jΔt+o(Δt)
Healthy dispersal from patch j to k	S_j→S_j-1,S_k→S_k+1	djkS(t)s_jΔt+o(Δt)
Healthy dispersal from patch k to j	S_j→S_j+1,S_k→S_k-1	dkjS(t)s_kΔt+o(Δt)
Infected death in patch j	I_j→I_j-1	(μ_j+α_j)i_jΔt+o(Δt)
Recovery in patch j	I_j→I_j-1,R_j→R_j+1	γ_ji_jΔt+o(Δt)
Infected dispersal from patch j to k	I_j→I_j-1,I_k→I_k+1	djkI(t)i_jΔt+o(Δt)
Infected dispersal from patch k to j	I_j→I_j+1,I_k→I_k-1	dkjI(t)i_kΔt+o(Δt)
Recovered death in patch j	R_j→R_j-1	μ_jr_jΔt+o(Δt)
Recovered dispersal from patch j to k	R_j→R_j-1,R_k→R_k+1	djkR(t)r_jΔt+o(Δt)
Recovered dispersal from patch k to j	R_j→R_j+1,R_k→R_k-1	dkjR(t)r_kΔt+o(Δt)
Sum of Changes		Σ(t)Δt+o(Δt)

For example, a new infection in patch j at time t results in a change in the random variables, and with probabilityThe lower case letters , , and denote the values of the random variables , , and , respectively.

A simple example of two patches illustrates the dynamics of the stochastic time-nonhomogeneous process when either a disease outbreak or disease extinction occurs. In Fig. 1, two sample paths are graphed for . 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.

Fig. 1

(Color Figure Online) Two sample paths are graphed with one infected individual introduced into patch 1, . Patch 1 is high-risk with and patch 2 is low-risk with . 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 , and

Branching Process Approximation

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

Table 2

Infinitesimal transition probabilities for the infected states in the branching process approximation of the stochastic n-patch model

Event	Description	Transition	Probabilities
1	Infection in patch j	I_j→I_j+1	β_j(t)i_jΔt+o(Δt)
2	Infected death in patch j	I_j→I_j-1	(μ_j+α_j)i_jΔt+o(Δt)
3	Recovery in patch j	I_j→I_j-1	γ_ji_jΔt+o(Δt)
4	Infected dispersal from patch j to k	I_j→I_j-1,I_k→I_k+1	djkI(t)i_jΔt+o(Δt)
5	Infected dispersal from patch k to j	I_j→I_j+1,I_k→I_k-1	dkjI(t)i_kΔt+o(Δt)
Total	Sum of Changes		Σ_I(t)Δt+o(Δt)

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 from the ODE model serves as a threshold for the branching process approximation. In particular, they showed if , 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.e., and . If , then the asymptotic probability of disease extinction is 1, but if , thenIn 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 aswhere the notation , , is the number infected at time t and , , 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 is the standard unit vector in .

Next, we define the offspring pgfs, one for each of the random variables . The offspring pgfs for this time-continuous process are distinct from those in a discrete-time Markov chain (DTMC). The number of infections generated by in the time-continuous process is assumed to occur in a small period of time rather than during the entire infectious period as in a DTMC. Given and for , one of four events defined in Table 2, if any, occurs in a small period of time . 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 that includes these four events is defined bywhere , Each pgf has the property that and is nonlinear in u and positive on the set . 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. Letwhere . In the branching process, the random variables are assumed to be independent, resulting inThis is a reasonable assumption if the initial number of infected individuals is small. Let and the remaining infected states , Then, for , 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 . 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 functionfor (Bacaër and Ait Dads 2014). In particular, the probability of an outbreak given , equals

Numerical Methods

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 in the differential equation (13) and for simplicity, we fix , , so thatThe preceding differential equation is solved (backward in time) for sufficiently large so that there is convergence to a periodic solution,Another change of variable yields .

For each numerical example, we estimate 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 , chosen sufficiently small, such that during each time step, only one, if any, of the events in Table 1 occurs (). The accuracy of the Monte Carlo approximation is tested by using a sequence of progressively smaller time steps. A total of sample paths are simulated for a given set of initial conditions at a specific time . Each sample path continues until a time when either or . If the total infected population reaches 50, then it is assumed that there is an outbreak. The proportion of sample paths q out of that reach zero is an estimate for the probability of extinction, and the remaining proportion that reach 50 is an estimate of the probability of an outbreak.

Two Patches

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

Fig. 2

(Color Figure Online) Compartmental diagram for the two-patch model

Sinusoidal functions are assumed for the transmission and dispersal rates with period : . The amplitude of these rates are and , respectively. The dispersal rates in each patch are the same, except for the value of which may be either positive or negative while the transmission rates differ in the magnitude of . If the values and 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.

Fig. 3

(Color Figure Online) An illustration when patches 1 and 2 are either 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 . 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 or and but keep , constant, i.e., . Each patch has the same recovery and death rates but differs in the average transmission rates, , . Patch 1 is high-risk with patch reproduction number and patch 2 is low-risk with (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.

Table 3

Baseline Parameter values for the two- and three-patch models where the units are per season

Description	Parameters	Two patches, i=1,2	Three Patches, i=1,2,3
Total population size	N_i	2000	2000
Natural death rate	μ_i	0, 0.01	0, 0.01
Disease-related death rate	α_i	0, 0.5	0, 0.5
Recovery rate	γ_i	6-μ_i-α_i	2-μ_i-α_i
Average transmission rates	β¯_i	18, 1.2	10, 1, 0.5
Average dispersal	d¯ijℓ,ℓ=S,I,R	2	2
Amplitude of transmission rate	δ_i	0, 0.8	0, 0.8
Amplitude of dispersal rate	δijℓ,ℓ=S,I,R	0, ±0.8	0, ±0.8
Length of period	T	4	4

Fig. 4

(Color Figure Online) Two patches with one high-risk patch, , and one low-risk patch,

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 and (0, 1), respectively. The branching process estimate is verified by checking thirteen different times with a Monte Carlo simulation of sample paths of the nonhomogeneous stochastic process. The results from the Monte Carlo simulations show good agreement with the branching process approximation.

Fig. 5

(Color Figure Online) Two-patch model with the periodic transmission (black solid curve) and dispersal (red dashed curve) rates, graphed in each of the top two (a)–(d), and probability of an outbreak, graphed in each of the bottom two (a)–(d) for patches 1 and 2 when or (0, 1), respectively. In a constant transmission with periodic dispersal, b constant dispersal with periodic transmission, c transmission synchronized with dispersal and d transmission desynchronized with dispersal. Thirteen different introductions of one infected individual into either patch 1 or 2 at times are checked using sample paths of the Monte Carlo simulation of the full nonhomogeneous process (Num Sim, circles). The parameter values are given in Table 3

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 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 ), the patch j, the location of the infectious individual and the value of the transmission rate 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 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 is not a good predictor of disease severity in these cases as the largest value of has the lowest average probability of an outbreak.

Table 4

Average probability of an outbreak for patch 1 and for patch 2 and the threshold value for the two-patch system when and

Two-patch	β_j	djkI	R_0	Patch 1	Patch 2
(A)	Constant	Periodic	2.45	0.59	0.15
(B)	Periodic	Constant	2.46	0.44	0.13
(C)	Periodic	Synchronized	2.41	0.45	0.16
(D)	Periodic	Desynchronized	2.58	0.44	0.09
Base	Constant	Constant	2.41	0.59	0.17

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) equalsThe probability of an outbreak is close to 1 when and but drops rapidly to for . 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. 6

(Color Figure Online) For synchronous dispersal as in 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, , respectively ). The blue circles are the results of the Monte Carlo simulation of the nonhomogeneous stochastic process with sample paths at The parameter values are given in Table 3

Next, we examine how an increase in the dispersal rate affects the probability of an outbreak (Fig. 7). As the average dispersal rate 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).

Fig. 7

(Color Figure Online) The probabilities of an outbreak for initial conditions (patch 1, solid curves) or (0, 1) (patch 2, dashed curves) are graphed when the average dispersal rate is increased , , 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 as increases

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, , The dispersal parameters for the two-patch model arewith remaining parametersfor 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 time-nonhomogeneous process, the initial conditions for the population sizes are equal to the disease-free population sizes. That is, takes the value of in Fig. 8 for When infection is introduced, or (0, 1), then , and , .

Fig. 8

(Color Figure Online) Disease-free solution of the two-patch model is graphed for the ODE two-patch model (black solid and red dashed curves) and one sample path of the stochastic time-nonhomogeneous process (green curves). Initial conditions are and . Parameter values are summarized in the list (16)

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; at the initiation of an epidemic.

Fig. 9

(Color Figure Online) Periodic probability of a disease outbreak for two patches with oscillating population sizes. Parameter values are summarized in the lists (16) and (17)

Three Patches

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., (Fig. 10). Patch 1 is high-risk with , and patches 2 and 3 are low-risk with and . The second type of connection is called a circular connection (CC), i.e., . 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., . The arrangement is changed in the third case, the second patch has and . The periodic functions for dispersal and transmission are given in Eq. (15) with parameter values in Table 3.

Fig. 10

(Color Figure Online) Three patches with either full connection (FC), circular connection (CC) or bidirectional connection (BC) containing one high-risk patch and two low-risk patches, and (FC and CC) or and (BC)

The periodic probabilities of an outbreak with FC are illustrated in Fig. 11 when initial conditions , 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, . 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).

Fig. 11

(Color Figure Online) The periodic probabilities of an outbreak for FC with the four cases of dispersal and transmission (a)–(d) as in Fig. 5. The initial conditions are and (0, 0, 1) for patches 1, 2, and 3, respectively. Patch 1 is high-risk, , and patches 2 and 3 are low-risk, and . Parameter values are given in Table 3

Fig. 13

(Color Figure Online) The periodic probabilities of an outbreak for the four cases of dispersal and transmission rates, (a)–(d) as in Fig. 5, with bidirectional connection (BC). The initial conditions are , and (0, 0, 1) for patches 1, 2, and 3, respectively. Patch 1 is high-risk, , and patches 2 and 3 are low-risk, and

Fig. 14

(Color Figure Online) The periodic probabilities of an outbreak for the four cases of dispersal and transmission rates, (a)–(d) as in Fig. 5, with circular connection (CC). The initial conditions are and (0, 0, 1) for patches 1, 2, and 3, respectively. Patch 1 is high-risk, , and patches 2 and 3 are low-risk, and

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 , but the average probabilities are generally lower than for constant transmission and dispersal.

Table 5

Average probabilities of an outbreak in three patches with FC for the four cases in Fig. 11

FC	β_j	djkI	R_0	Patch 1	Patch 2	Patch 3
				R_01=5	R_02=0.5	R_03=0.25
(A)	Constant	Periodic	2.79	0.64	0.34	0.32
(B)	Periodic	Constant	2.75	0.52	0.31	0.30
(C)	Periodic	Synchronized	2.60	0.55	0.33	0.32
(D)	Periodic	Desynchronized	3.21	0.53	0.25	0.23
Base	Constant	Constant	2.60	0.63	0.36	0.35

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 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).

Table 6

Average probabilities of an outbreak in three patches with BC and CC for the four cases (A)–(D) corresponding to Figs. 13 and 14

Case	BC				CC
	R_0	Patch 1	Patch 2	Patch 3	R_0	Patch 1	Patch 2	Patch 3
		R_01=5	R_02=0.25	R_03=0.5		R_01=5	R_02=0.5	R_03=0.25
(A)	3.29	0.70	0.29	0.18	3.09	0.68	0.21	0.33
(B)	3.28	0.58	0.27	0.17	3.08	0.56	0.19	0.31
(C)	3.17	0.61	0.29	0.20	2.92	0.59	0.22	0.33
(D)	3.61	0.57	0.21	0.12	3.50	0.55	0.14	0.24
Base	3.17	0.69	0.31	0.20	2.92	0.67	0.23	0.36

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) (top), (ii) , (ii) , (iv) (bottom), . Here, we assume natural and disease-related death rates are so that . As the recovery rate increases from (i) to (iv), the average duration of the infection shortens, each patch reproduction number 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 when transmission is decreasing.

Fig. 12

(Color Figure Online) Probability of an outbreak with synchronized dispersal in three patches with FC, BC, and CC as in Fig. 10. Initial conditions in patch 1 are (solid black curve), in patch 2, (red dashed curve) and in patch 3, (blue dotted curve). Values of the recovery rate in rows (i) , (ii) , (iii) and (iv) , with natural and disease-related death rates set to zero . All other parameters are in Table 3

Table 7 is a summary of the average probabilities of an outbreak and the values of 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, . 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 equalswhere . 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.

Table 7

Average probability of an outbreak in three patches for FC, BC, and CC with synchronized dispersal and transmission, corresponding to Fig. 12

Patch	FC				BC				CC
	1	2	3	R_0	1	2	3	R_0	1	2	3	R_0
(i)	0.86	0.74	0.72	8.46	0.88	0.69	0.64	9.51	0.88	0.66	0.72	8.86
(ii)	0.75	0.55	0.54	4.59	0.78	0.50	0.41	5.44	0.77	0.44	0.54	4.97
(iii)	0.55	0.33	0.32	2.60	0.61	0.29	0.20	3.17	0.59	0.22	0.33	2.92
(iv)	0.28	0.12	0.12	1.53	0.36	0.12	0.06	1.85	0.33	0.07	0.14	1.74

Summary and Conclusion

Seasonal changes that affect the pathogen life cycle or the within-host and between-host 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 time-nonhomogeneous 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 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 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.