Multi-patch epidemic models with partial mobility, residency, and demography

The emergence and re-emergence of infectious diseases has been a global cause of concern in the past few decades. Previous research in the field has revealed that human connectivity and mobility behaviour play a major role in the spreading of an infectious disease. In this work, we propose multi-patch models that take into account the effects

of human mobility on the evolution of disease dynamics in a multi-

Introduction
In the current century, knowledge about the risks of the hosts coupled with easier and faster communication of disease risk information by health officials is very important in enabling individuals to prevent the spread of infectious diseases. Another important component of disease risk assessment by individuals and hosts is the ability of public health officials to measure and properly communicate, in a timely manner, real or perceived information on disease risk [3]. One way to conduct such analysis is through the use of meta-population studies [1,4,5]. This can help assess the effect of movement of individuals between different regions or cities (in general referred as patches) as far as the spread of an infectious disease is concerned, taking into account the disease risk in each of the cities/patches. A meta-population is network of patches among which individuals can travel. This approach has been used in the study of the spread and risks of various infectious diseases.
Most infectious diseases, for instance, SARS, SARS-COV-2, H1N1 etc can be spread among individuals through their movements, known as mobility. At the forefront of global concern are diseases that can be spread through various forms of mobility including travel and trade [3,6]. Consequently, it has increasingly attracted the interest of theoreticians and epidemiologists to quantify the impact of mobility on infectious disease dynamics [7]. Indeed, epidemiologists are tasked to understand the cause of a disease, to predict its course, and to develop approaches capable of preventing and controlling outbreaks/infections. The potential impact of mobility and dispersal of population in propagating and fanning disease spread between neighboring environments cannot be underestimated. In fact, mobility of individuals in face of an infectious disease can lead to an epidemic or even a pandemic.
The potential threat of mobility has been the focus of many research works which attempt to quantify and measure impact of mobility on disease spread within different environments. For instance, [8] proposed a patchy SIS model with standard incidence to asses how population dispersal affects the distribution and overall size of infection in the environment. Through simulation, the author revealed that fast diffusion decreases the basic reproduction number but has a potential of increasing the total infection size. Another fundamental finding is that patches with higher infection risks have higher disease prevalence in the face of human dispersal for the two-patch case, than for more than two patches. Clearly, only a fraction of people are frequent travelers while most people travel occasionally or never at all. Using this observation, [8] developed an SIS multi-patch epidemic model which captured the traveling frequency and showed, using both analytical and numerical procedures, that the model tends to underestimate the infection risk when the difference in human travel frequency is not captured and not clearly distinguished. The model led to the conclusion that travelers residing in regions with high transmission rates are at the highest risk of infection. However, the work by [8] did not incorporate the time that the proportion of the traveling individuals spend in the patches they have traveled to.
[9] considered a two-patch SIR epidemic model with state dependent residence-time matrix to asses the impact of virtual dispersal on disease transmission dynamics and optimal resource allocation in the two patches. They established that patch-specific optimal control strategies reduce the number of infections per patch. The use of state dependent residence time matrix was also one of the goals of [1] in their SIS multipatch model, although their focus was not on optimal resource allocation but a comparison of local and global disease dynamics from a constant and state-dependent residence time matrix.
Multipatch epidemic models have also been used to study the effects of media coverage on the the propagation of infectious disease in a multi-patch environment. For instance, [10] considered a patchy SIS model with varying transmission rates to investigate the effects of media coverage and human mobility on the transmission dynamics of infectious diseases among patches. They computed the basic reproduction number R 0 of their model and established that the Disease Free Equilibrium (DFE) is Globally Asymptotically Stable (GAS) when R 0 ≤ 1 and that when the travel rates of susceptible and infectious individuals are equal, then the disease is uniformly persistent and there exists a unique and GAS endemic equilibrium when R 0 > 1. [11] constructed two SIS epidemic models in a patchy environment; one with and the other without linear recruitment. They established that both models have the same basic reproduction number R 0 which they used to asses the global stability of the models. Moreover, they established that even though R 0 is an important threshold value, other factors can play an important role in the prediction and control of infectious diseases. In particular, their results showed that the variation of total population number can cause the disease to be more threatening and difficult to control even when mobility of infected and/or susceptible population tend to zero. [12] presented some analytical results for an SIS epidemic model that describes the spread of an infectious disease in a population of individuals characterized by travel between n-cities. They computed and obtained explicit bounds on R 0 in the presence of homogeneous interaction between individuals in each city. Their numerical simulations revealed that the disease dies out when R 0 ≤ 1 and remains endemic when R 0 > 1. The same observation was made by [13] who considered a SEIRS epidemic model which captures cross infection between several species and track both the current patch and the patch in which an individual usually resides, to investigate disease spread in an environment divided into p patches. [14] proposed a multi-patch SIR-SI model to study the effects of human travel on the host vector dynamics of dengue disease in different patches. They showed that the disease dies out if R 0 < 1. Otherwise, the disease remains endemic if R 0 > 1. Their results showed that controlling mobility from disease dominant patches to patches where the disease is less prevalent may result in controlling the disease in the low-disease prevalent patches while the disease remains more prevalent in the high-disease prevalent patches.
Multi-patch epidemic models have also been used to fit COVID-19 data. For instance, [15] constructed a generalized n-patch SEIR epidemic model with quarantine and hospitalization. They concluded that appropriate quarantine strategy and control of the migration rate are important to reduce the spread of covid-19 after numerically fitting their model to COVID-19 data from three patches: Hubei, Chongqing and Hunan provinces in china. They also established that their model has a unique DFE which is GAS when R 0 < 1 and is unstable when R 0 > 1. [16] studied the global stability of a multi-patch waterborne disease epidemics model in a multi-patch network. Their model encapsulates the dynamics of the susceptiple, asymptomatic and symptomatic individuals as well as the bacteria dynamics at interacting patches. The patches are connected by human mobility. They constructed Lyapunov functions to study the global stability properties of the disease free and endemic equilibrium of their model. Transmission of Ebola virus in a multi-patch environment was studied by [17]. The set-up was such that humans carrying the pathogen can move between patches, transmitting the disease in any patch. Their model captures the dynamics of people who are susceptible, infectious, isolated, deceased but not yet buried and still infectious as well as the dynamics of the pathogen at interacting patches. The author constructed Lyapunov functions to establish global stability properties of both the disease free and endemic equilibrium of their model. [18] constructed a SEIRS epidemic model for a spatially heterogeneous environment with patches in which human mobility rates between patches depends on disease statuses. They computed certain bounds on R 0 and showed that the DFE of their model is GAS when R 0 < 1. Using two patch example, they showed that propagation dynamics of a disease can greatly be influenced by travel of infectious individuals in a patchy environment. [19] modelled the transmission dynamics of tuberculosis in Africa using a patcy model. they computed R 0 and demonstrated that the DFE is GAS when R 0 < 1 and unstable when R 0 > 1, a condition which leads to the existence of one unique endemic equilibrium.
Due to the ever increasing of human population together with advances and efficiency in transportations, the connectivity that we are constantly gaining between previously separated parts of the world is one of the main factors that promotes an infectious agent spreading. Therefore, understanding the patterns of interactions between humans, between humans and vectors and more importantly, the pattern of the movement of individuals is very important to develop health and public policies [3]. [20], by modeling mobility of individuals between a large urban center and smaller satellite centers using a SIR epidemic model, discovered that most of the behaviour of the whole system is governed by the large urban center. In their analysis, they also concluded that the population size of a satellite city and its connectivity to the main urban center significantly contributes to the disease evolution of the whole system. [1] examined the effects of mobility on disease levels in an environment divided into n patches using SIS and SIR models. They introduced mobility residence matrix which quantifies the effective contact time that individuals have with the residents of the patches they move into. [2] proposed a multi-patch and multi-group epidemic model which captures the heterogeneity of the epidemiological groups in n patches. They used this model to asses the effects the mobility of susceptible, infected, exposed, and recovered individuals on the diseases levels in the various patches they move into/from. Inherent in the models is the assumption that all individuals travel and have a similar visiting times to other patches. This assumption is clearly not realistic. In all of these works, non has considered mobility of individuals in real time. All of them, except [8] models the effects of virtual human dispersal on disease dynamics in patchy environments. In other words, the existing literature captures the proportion of the population that travels and the time time they spend in the patch they have moved to using one parameter; the residence time matrix, which other literature calls the travel rates.
Our contributions are as follows. Using SEIRS model, We have generalized the existing multi-patch SIS and SIR epidemic models in [1] in several important directions. First, we separate the mobility and residence time for each patch, thus leading to a more practical and general model. Second, we have performed a complete rigorous analysis of our proposed model accompanied by extensive numerical results. In particular, we have shown that the model is locally Lipschitz and is biologically well posed. Under the proposed formulation, we compute the global and patch specific basic reproduction numbers R 0 as a function of the mobility parameters and the residence time matrix matrix. We perform global analysis where we have ascertained that in the presence of mobility, the disease dynamics purely depends on the global basic reproduction number R 0 while the patch specific reproduction numbers R i 0 controls Patch i disease dynamics in the absence of mobility. We have shown that the DFE is GAS if R 0 ≤ 1 and that there exists a unique EE which is asymptotically stable when R 0 > 1. Third, a novel definition of the patch reproduction number R i,i 0 has also been introduced where we have shown, theoretically and through numerical simulation, that there is a close relationship between R 0 and R i,i 0 as regards to the persistence and extinction of the disease. Our newly defined patch reproduction numbers obeys the expected theoretical behaviour of reproduction numbers in this regard as is discussed in the text. Fourth, we have also introduced and analyzed a more generalized epidemiological-statusmobility-dependent SEIRS epidemic model which encapsulates the idea that Patch i disease levels can be determined when a proportion of the population in the respective Patch Susceptible, Exposed, Infected and Recovered epidemiological groups travel and spend a proportion of their time in the other patch.
The rest of this paper is organized as follows: In Section 2, we derive the generalized n-Patch SEIRS model with mobility, residency and demography after which, we explore the biological well-posedness of the model. Section 3 introduces, computes and analyzes relations of the global and different local basic reproduction numbers. In Section 4, we consider the global stability of the DFE, the existence and local asymptotic stability of a unique endemic equilibrium (EE). Section 5 explores, through small scale and large scale simulation, the effects of heterogeneity in mobility and residence times on the patch and global environment disease dynamics. In Section 6, we introduce and analyze a more generalized epidemiological-statuses-mobility-dependent SEIRS model. Finally, Section 7 concludes the paper.
2 An n-Patch SEIRS model with mobility, residency and demography In this section, we first devise an n-Patch SEIRS model taking into account mobility, residency and demography. We then rigorously study the wellposedness of the model, the disease free equilibrium (DFE), existence of the unique endemic equilibrium (EE), and global stability of the equilibria.

Epidemic Model
Consider a population with n patches, with {N i } n i=1 individuals. The disease evolution of each patch is assumed to follow the SEIRS model with mobility and demography. We assume that a proportion α i of patch i residents move to and spend p ij ∈ [0, 1] fraction of their time in patch j. In [1], p ij was used as a "proxy" to measure the contact time residents in patch i move and stay in patch j. In their analysis p ij has the inherent assumption that all patch i individuals move to patch j and have p ij contact time with residents of patch j. This assumption, though greatly simplifying the analysis, is not realistic. Clearly, there is always a percentage of (in fact most) individuals remaining at home as the others travel for whatever reasons. Introducing the mobility parameter α i is our contribution to address this shortcoming. The combination of α i and p ij enables us to model not only any fraction traveling out of a patch but also different sub-fractions spending time on other patches. Therefore, if we consider tha individuals move independently, there are α i N i p ij patch i residents in patch j at time t, and there are (1 − α i )N i patch i residents who do not move out from their patch.
Since p ij is the fraction of unit of time (days, for instance) that an individual i, who travels out of his patch, spends it in patch j, then n j=1 p ij = 1 for i = 1, . . . , n.
The effective population in patch j is given as Total population of people who stays in patch j + n k=1 α k p kj N k Total population of people who move to patch j Similarly, the number of infected individuals in patch j is The proportion of infected individuals in patch j is therefore given as Notice that the number of susceptible individuals from patch i who are currently in patch j are α i p ij S i . Multiplying this with the effective proportion of infected individuals in patch j and patch j risk/rate of infection β j , we obtain the number of S i who is infected per unit time in patch j. It follows that the total number of susceptible individuals of patch i who are infected due to traveling reads #Si infected per unit Due to individuals in i, who travels out of its patch, can also spend a fraction of their time in i, the sum in (1) does not exclude j = i.
Since total number of susceptible remaining in patch i is (1 − α i )S i , the number of patch i susceptible infected per unit time while remaining in patch i is #Si infected per unit time while remaining in patch Therefore, the total number of infected susceptible individuals in patch i per unit time is given as #Si infectedd per unit time while remaining in patch i + #Si infected per unit time while traveling We assume that the dynamics of each patches follow the SEIRS model with mobility and demography. That is, we incorporate vital dynamics such as birth (recruitment), natural and disease induced death rates, and the concept of waning immunity. Waning immunity simply means that the recovered individuals loose immunity and become susceptible once again. We thus obtain a patchy SEIRS model with mobility, residency and demography as where we have definedp kj = α k p kj and i ∈ [n] := {1, . . . , n}. The parameters used in the model are described in Table 1. Per capita rate at which the exposed individuals in patch i becomes infectious

Well-posedness of the model
A dynamical system or a system of differential equations is said to be wellposed if ∀t, the system has a unique solution that depends continuously on the data. An additional condition for a dynamical system describing population dynamics is that, given any non-negative initial conditions, the solutions of the system should remain positive at all times [21]. Our proof for the wellposedness of the proposed model (3) relies on the Picard-Lipchitz theorem, which we now state.
Lemma 2.1 (Picard-Lipschitz theorem, [22]). Let t 0 ∈ R and x 0 ∈ R n . Consider the following initial valued problem where f : R × R n → R n is Lipschitz continuous with respect to x on Then there exists δ > 0 such that the system (4) has a unique solution for any t ∈ [t 0 − δ, t 0 + δ] and the solution depends continuously on [t 0 , x 0 ].
Corollary 2.1. Given an initial condition (S 0 , I 0 , I 0 , R 0 ) at t 0 . There exists δ > 0 such that (5) has a unique solution continuous with respect to Proof The result is a direct consequence of Lemma 2.1. We only need to show that the right hand side of the system (5) is locally Lipschitz continuous, but this is obvious as the right hand side is the composition of Lipschitz functions. □ We next show that system (3) is biologically well-posed. That is, given any non-negative initial conditions, the solutions of the system remains positive at all times.
Theorem 2.1. Suppose that initial conditions at t 0 = 0 for system (3) are non-negative. The states E i (t), I i (t) and R i (t) ∀i ∈ [n] := {1, 2, · · · , n} remains non-negative, while S i (t) and N i (t) remains positive at all times. Furthermore, the patch total population N i (t), i ∈ [n], is bounded above by max N i (0), Λ i /µ i .
Proof The proof is inspired by [18]. Suppose S i , for some i ∈ [n], vanishes at t ≥ 0 before other variables. We conclude from the first equation of (3) that S i (t) ≥ 0 at all times sinceṠ i (t) = Λ i + τ i R i ≥ 0. Similarly, if E i , for some i ∈ [n], vanishes at t ≥ 0 before other variables, then from the second equation of (3), we havė and thus E i (t) can never be negative. In the same vein, it is easy to show that I i (t) ≥ 0 and R i (t) ≥ 0 at all times. On the other hand, from the first equation of system (3), we have, for any i ∈ [n], which in turn implies S i (t) > 0 as long as S i (0) > 0. Therefore, N i (t) > 0, ∀t ≥ 0 as long as S i (0) > 0. Now, summing all equations of system (3) From the last result, the environment total population Now we define the patch population vector as N = N 1 N 2 · · · N n T and adding all equations in (5) altogether we obtain dynamics of the vector patch population vector aṡ It is easy to see that the set is a global compact attracting positively invariant set for the system in Equation (5). Solutions of the system (5) will thus approach to and stay in D SEIRS . Recall that the disease free equilibrium (DFE) of system (5) is the point where the population is free of disease. The below characterization of the DFE is straightforward.
Proposition 2.1. The DFE of system (5) is given as

Global and Patches' Reproduction numbers
Reproduction numbers are important threshold parameters in mathematical epidemiology as their values help us determine if an infectious disease outbreaks or dies out. This section derives the expressions for the global and local basic reproduction numbers R 0 and propose some reproduction numbers considering the sub populations, and establish some important relation between them and R 0 .

Basic reproduction number R 0
The global behaviour and disease transmission dynamics of system such as (5) can typically be characterized by the basic reproduction number R 0 . Following [23,24] we use the next generation to compute this crucial threshold quantity. Clearly, the infected compartments in our model are E and I, corresponding to the second and third equations of the system (5). We decompose the right hand side of these equations as F (E, where F (E, I) is the rate of appearance of new infections and V(E, I) is the outgoing rate of transferring already infected individuals. Let F and V be the Jacobians (with respect to (E, I) of F (E, I) and V(E, I), respectively) at the diseases free equilibrium (DFE) P 0 .
The inverse V −1 is given as and the next generation matrix (NGM) reads where we have defined The global basic reproduction number is therefore given by the spectral radius of the block matrix in the first row first column of F V −1 . That is, It can be shown (see, e.g., [24] and the references therein) that when R 0 > 1, the infection can spread to a large fraction of the population, and afterwards go into extinction, or remain endemic in the global population (or any of its respective patches). In the next sections we present three reproduction numbers that describe the number of cases that individuals in a patch can produce and describe their relations with the global basic reproduction number.

Only-local basic reproduction number
In the absence of mobility the patches act independently and then, for i = 1, 2, · · · , n, the Patch i has basic reproduction numbers given by and the disease remains endemic in the population if R i 0 > 1 for at least one i. The global basic reproduction number in the presence of mobility R 0 (8) can be bounded by the local basic reproduction numbers R i 0 [see e.g. 18,25] in the special case when the demographic and epidemiological parameters of the patches differ only in the transmission rates β i , i = 1, 2, · · · , n.
n} with the epidemiological parameters of the patches only differing in their contact rates β i , i = 1, 2, · · · , n, then Proof Without loss of generality, we consider Based on our assumption that the epidemiological parameters of the patches only differ in their contact rates, we have that G = P * diag(β)diag(P * T 1) −1 P * T . Now we examine the sum of the columns of the next generation matrix F V −1 in (6). Let C 1 be the sum of the first column, then with C 1 being the first element of C. Define r = κ (κ+µ)(γ+ψ+µ) , then we have where the inequality follows from (10). We therefore have that In a similar fashion, we have . . , n. For j = n + 1, . . . , 2n, it is also straightforward to verify that Since the spectral radius of a non-negative matrix is bounded below and above by the minimum and maximum column sums respectively [26], we have that We can write the NGM as where r ij can be interpreted as the number of expected new cases in patch j that an infective individual in i produces. One may attempt to define the patch i reproduction numbersR i 0 as the number of new cases that an infective individual from patch i can produce (in early stages of the outbreak). Hencẽ R i 0 = n j=1 r ij , i = 1, . . . , n, and R 0 = ρ(M) as usual. The problem with this definition of individual patch reproduction numbers is that while it is easily computed, we cannot always determine if R 0 is greater or smaller than the critical value 1. Indeed, we have the following observations.
Since the NGM is a nonnegative matrix, that in several applications can be considered irreducible, using the Perron-Frobenius theorem we have The Perron-Frobenious theorem provides simple bounds for the spectral radius of a matrix based on its row (column) sums, however it cannot specify if R 0 is above or below threshold if only some of theR i 0 are greater than 1. In these cases we can appeal to sharper bound for the spectral radius. Although various proposals have been obtained, one attraction they can have is that they are given in simple operation of the matrix elements. On this direction, a theorem that allows knowing the location of all the eigenvalues is the Gershgorin circle theorem that in contrast to the Perron-Frobenious, is valid for reducible or irreducible matrices. This theorem states that given Sharper lower bounds for the spectral radius of a real matrix are presented in [27] that use the traces of the matrix and its square. [28] also provide bound for the eigenvalues and when the matrix is nonegative, the lower bound coincides with the optimal lower bound given in terms of n, the traces of the matrix and its square. We present this result.
Theorem 3.2 (Lemma 7 and and Theorem 8 in [28]). Let A ⩾ 0 with a := tr A and b := tr A 2 being real numbers satisfying a 2 ≤ nb then is optimal (i.e., the best possible lower bound for r, using only n, a, and b) if and only if l ⩾ cb/n. For this matrix we have that the row sums are 0.3961 and 1.1819 and we cannot conclude that R 0 = ρ(M) > 1. Based on the Gershgorin cirle theorem we have that R 0 ∈ (−0.0447, 0.4150) ∪ (0.7411, 1.16309) that contains values above 1. However, using Theorem 3.2 we know that R 0 > 1.01082.

R i,i 0 reproduction number
Some recent works, [see, e.g., 1, 2] resorted to the i th diagonal element of the Next Generation Matrix (NGM) as a measure of the reproduction in patch i, but did not provide a clear rational and intuition to this selection. Moreover, in some cases, at least in the numerical simulation study, contradicted the intuitive behaviour of the basic reproduction number in regard to the persistence (or the extinction) of the disease. Indeed, [1] stated that R 0,i > 1 is necessary but not a sufficient condition for the persistence of the disease in Patch i.
In this section we utilize the standard definition of the NGM, represented as in (12), to define a new local reproduction number R i,i 0 for patch i. This reproduction number corresponds to the number of expected new cases in Patch i originated by a typical case in i either by direct contact or after the infection has transited to individuals in other patches. We also establish the relations between this individual patch reproduction number and the global reproduction number.
This reproductive number has the motivation of regarding Patch i as of the primary interest and the rest of individuals in the patches as vectors [29]. For example, in a two patches population with irreducible NGM, the number of new cases generated by a single infective v in the patch 1 (P1) correspond to those individuals in P1 (⃝) who are directly infected by v plus all those individuals from P1 directly infected by secondary cases in P2 (□), for which the chain of infections has not yet included individuals in P1. See Figure 1. For this example, the number of new cases that we expect a typical infective originates in P1, that we denote as R 1,1 0 , is R 1,1 0 = r 11 + r 12 r 21 + r 12 r 22 r 21 + r 12 r 22 r 22 r 21 + r 12 r 22 r 22 r 22 r 21 + · · · (14) = r 11 + r 12 r 21 where (15) follows when r 22 < 1. We see that even if all the elements in M are smaller than 1, R 1,1 0 could be greater than 1. This scenario corresponds to the outbreak due to the initial infectious cases in P1. If r 22 > 1 and r 12 , r 21 > 0, it is easy to see that the increasing number of infectious cases in P2 will keep contributing to the infected cases in P1. Clearly, if either r 12 or r 21 are null, then all the new cases in P1 are originated only by direct infections, and in that case R 1,1 0 = r 11 . When the population consists of n patches, the result can be generalized to where (17) holds if all the eigenvalues of the i th order principal submatrix M[−i] (the NGM without its row and column i) are between -1 and 1. We have defined r i: (17) can be written in terms of its eigenvalues and eigenvectors as where U is the matrix whose columns are the eigenvectors of M[−i] and D is the diagonal matrix with elements [D] jj = 1/(1 − λ j ), j = 1, . . . , n.
We next discuss the relations between R 0 and {R i,i 0 }. We begin by recalling some useful definitions and results, which will be useful in establishing our main results in this section. 26,30]). Let A ∈ R n×n be a non-negative matrix and B ∈ R (n−1)×(n−1) be the principal sub-matrix of A. Then, ρ(B) ≤ ρ(A). If A ∈ R n×n is irreducible and B ̸ = A, then ρ(B) < ρ(A).

Theorem 3.3. (Lemma 3 and Theorem 4 in [33]) If
A is a non-negative irreducible matrix, the Perron root of the generalized Perron complement, We now present our main results together with their proofs.
Proof Immediate from the Theorem 3.3. R i,i 0 can be thought of as a generalized Perron complement of M T [α] in M T with t = 1, β = {i} and α = ⟨n⟩ \ β. Now, for we have that the Neumann series in (16) diverges and since r i:−i , r −i:i > 0, R i,i 0 is unbounded from above. □ From this proposition the following corollary follows immediately.
To include the reducible case we formulate the next proposition.
Proof Let M be a nonnegative reducible matrix. Then there exists a permutation matrix P such that we obtain the normal form of a reducible matrix where M 1 , M 2, . . . , Ms are irreducible matrices.
Proof If for some i we have R i,i 0 > 1 and 1. R 0 > 1, then we have the conclusion by Proposition 3.1.

that also contradicts our assumption.
Hence the first case is the only one possible. □

GAS of the DFE, existence of a unique EE and it's local asymptotic stability
In this section we study the stability of the DFE of System (5) under the assumption that there there is no disease induced deaths (i.e., ψ = 0). In this case, the dynamics of the total population of System (5) becomeṡ We therefore have that lim t→∞ N = Λ⊙ 1 µ = N . Note that our system (5) satisfies the asymptotically autonomous condition in [35,36] and its limit equation is given by Similar to the proofs for the original system (5) (see section 2.2), we can show that limit system (22) is well-posed and its DFE is P 0 = (N , 0, 0, 0).

GAS of the DFE
We now study the stability of the DFE P 0 of the original system (5) via the stability of the DFE P 0 of the limit system (22). We use the same next generation approach (see section 3.1) to obtain the basic reproduction number of model (22). In this case, the F matrix remains the same while V now reads The basic reproduction number of the limit system (22) is then given as where we recall from (7): Theorem 4.1. Suppose P * P * T is irreducible, then the following hold: 1. if R 0 ≤ 1, the disease free equilibrium (DFE) is globally asymptotically stable 2. if R 0 > 1, the DFE is unstable Proof The second assertion is clear from [24]. For the first assertion, we adapt the proof in [2] for our system (22). Since −F V −1 is nonnegative and irreducible, R 0 = ρ − F V −1 is the largest positive eigenvalue of −F V −1 and its corresponding left eigenvector (z E , z I ) is positive, owing to the Perron-Frobenius theorem, [see 26]. In particular Due to the fact that −V −1 is nonnegative if we define then L(E, I) ≥ 0 and the derivative along any trajectory of (22) is given aṡ Since diag(S) ≤ diag(N ), we havė Thus, by the LaSalle's invariance principle, any solution of (22) approaches the largest invariant set A in the inverse imageL −1 (0). If R 0 < 1, thenL(E, I) = 0 if and only if I = E = 0, which is disease free, and thus A = {P 0 }. In other words, P 0 is globally asymptotically stable.
If R 0 = 1, a simple algebra manipulation giveṡ Thus,L(E, I) = 0 if and only if I = 0 or S = N . Either of the cases leads to the same conclusion that A = {P 0 }. That is, P 0 is also globally asymptotically stable when R 0 = 1. □ The following result is a direct consequence of Theorem 4.1 and results in [35,36].

Existence and asymptotic stability of the unique EE
We first show that for R 0 > 1 the original system (5) has a unique endemic equilibrium (EE) using its limit equation (22). To achieve this, we follow the approach in [2] to construct an auxiliary system which facilitates the uniqueness of the EE of System (22). In the following we denote by ≫ the componentwise strict inequality for two vectors.
Furthermore, E * is a fixed point of g(r) if and only if E * is an equilibrium of the systemṙ = F (r) := diag(κ + µ) g(r) − r .
Proof The second assertion is obvious and we sketch the (simple but lengthy) proof of the first assertion. Let (S * , E * , I * , R * ) be an endemic equilibrium point (I * ≫ 0) of system (22), then The last two equations of (26) provide the expressions for I * and R * in terms of E * , while the first two equations yields S * in terms of E * . Now, substituting I * , S * , and R * into the first two equations of System (26), and with some simple algebra manipulations we arrive at that is, E * is a fixed point of g. □ We now state and prove a theorem on the existence and uniqueness of EE of the limit system (22). Theorem 4.2. Assuming that P * P * T is irreducible, then the limit system (22) has a unique EE whenever R 0 > 1.
Proof Proposition (4.1) shows that the existence of the unique EE of system (22), when R 0 > 1, is equivalent to the existence of the unique non-trivial equilibrium of system (25). It is therefore sufficient to show the latter. We have and thus which is a Metzler matrix, since diag(N − Lr)P * diag(β)diag(P * T N ) −1 P * T J is nonnegative, and is a diagonal matrix. Additionally, F ′ (r) is irreducible owing to the irreducibility of P * . Therefore, F (r) is strongly monotone (see, e.g., [37] and the references therein). We observe that if a ≫ b, then F ′ (a) < F ′ (b). Hence, the map F ′ : R n → R n × R n is strictly antimonotone. Next notice that 0 is an equilibrium for (25) and Since has at least one positive eigenvalue and hence the origin 0 is unstable (see [38, page 310]). Therefore, according to [39, Theorem 6.1], there is a unique equilibrium different from 0 of the system (25), and hence a unique EE for the limit system (22)). □ We next show that the unique EE of the limit system (22) is asymptotically stable. We follow the procedure presented in [40] and adapt the proof in [41] to come up with the proof for our model. Proof Since we are dealing with the limit system, it suffices to consider the last three equations of system (22). That is, after replacing S = N − K where K = E + I + R, our system now becomes We obtain the positive EE points by solving for E * , I * and R * from the system from which we obtain that can be written as where is an irreducible matrix owing to the assumption of the irreducibility of P * P * T .
Notice that E * is the Perron-Frobenius vector of P . For any z ∈ C, we denote ℜz as the real part of z. Let Y 0 ∈ C 3n be the eigenvector corresponding to the eigenvalue {z ∈ C : ℜz ≥ 0} of the Jacobian of system (31) evaluated at the endemic equilibrium P * = (E * , I * , R * ). Using the Krasnosel'kiis trick in [42], we show the asymptotic stability of our system by showing that the linearized system has no solution of the form is the eigenvector corresponding to the eigenvalue z such that ℜz ≥ 0 , then we have We thus have that We eliminate V and W from the system by substituting for their values which we obtain from Equations (35b) and (35c), into Equation (35a), of which, after some rearrangements, we obtain Making the transformation
We assume that ℜz ≥ 0 and let ϑ(z) be the minimum of the real parts of the components of the vector ψ(z) and the magnitude |Ũ | = (|Ũ 1 |, · · · , |Ũn|). The magnitude of 37 gives 1 + ϑ(z) |Ũ | ≤ B|Ũ | Let ω be the minimum number such that By Theorem (4.2), E * is a strictly positive vector and hence ω < ∞. Now, sequentially using Equations (40) and (34) in Equation (39) we have We have from condition (38) that if ℜz ≥ 0 then ϑ(z) > 0. Equation (41) then leads to a contradiction of the assumption that ω is the minimum number such that 40 holds. We therefore conclude that ℜz < 0. This proves the (local) asymptotic stability of the EE. □ The following result is a direct consequence of Theorem 4.2, Theorem 4.3 and results in [35,36].
Corollary 4.2. If R 0 > 1, then there exists a unique (locally) asymptotically stable EE for the original system (5).

Numerical results
This section is dedicated to small scale and large scale numerical simulation study of the proposed model. The first subsection deals with numerical comparison of our model with that from [1]. The second sub-section deals with a detailed and in-depth numerical study on the role of mobility and patch residence times between highly risky and less risky patches. The numerical studies in the first and second subsections are based on a simple two-patch case with the last sub-section delving into large scale simulation of how mobility of individuals between the geostatistical units calledÁrea Geoestadística Básicas (AGEBs), influences disease dynamics in Hermosillo, the state of Sonora, Mexico. For the two-patch case, we have used the parameters β 1 = 2.0, β 2 = 0.5, γ 1 = γ 2 = 1/14, κ 1 = 1/15, κ 2 = 2/25, Λ 1 = Λ 2 = 9, µ 1 = 1/9, µ 2 = 1/8, τ 1 = 1/10, τ 2 = 1/20 and ψ 1 = ψ 2 = 0.003 and N 1 = Λ1 µ1 = 81 and N 2 = Λ2 µ2 = 72. The initial conditions used in the simulation were E 1 (0) = 10, E 2 (0) = 15, I 1 (0) = 5, I 2 (0) = 7, R 1 (0) = 0, R 2 (0) = 0, S 1 (0) = 66 and S 2 (0) = 50. There is no specific reason why these values were chosen. They are hypothetical and are meant to be used for illustrative purposes. The contact parameter values β 1 = 2.0 and β 2 = 0.5 were however chosen because we are interested in the assessment of the effects of mobility of individuals on the evolution of disease in high and low risk patches and on the global environment. These parameter values will be used throughout the simulation study in this subsection unless stated otherwise.

Comparison of actual mobility and virtual mobility
Using two patches, we compare the performance of our model with actual mobility against that proposed by [1] with virtual dispersal. As mentioned in the introduction, the model by [1] was constructed under the inherent assumption that the whole population can travel from their patch and spend homogeneous proportion of their time in any other patch. In order to account for the practical situations in which only a fraction of population traveling, we extend the model in [1] by introducing partial mobility through the parameter α i , i = 1, 2, · · · , n whose significance has been discussed in the introduction. Our model is thus more general in the sense that if we substitute α i = 1.0, i = 1, 2, · · · , n in our model and maintain the residence times P = (p ij ) j=1,··· ,n i=1,··· ,n , our model reduces to that of [1]. The infection curves from the two models are presented in Figure ( Figure (4) shows the global infection dynamics obtained by summing the infection dynamics in Figures (2) and (3). Observe that in both graphs, different combinations of mobility describe heterogeneity in the infection curves that would otherwise have been impossible to discern using the model in [1] (green curve with α i = 1.0). More importantly, our model generally reveals that a larger population spending a larger proportion of their time in a highly risky patch results into high levels of infections in that patch. This can be observed across Figures (2) and (3). In other words, it would be impossible to decipher the contribution of the mobility parameters α i if it was left coupled with the residence times as in [1].

Small scale numerical study: The two-patch case
For the numerical simulation of the proposed patchy SEIRS models, we have considered a simple two-patch case for which, our intention is to perform an in-depth analysis of the effects of the mobility parameters α i and the residence time matrix P on the evolution of the disease in each patch and global environment based on the basic reproduction numbers of the models. The idea is to asses the effects of mobility between low-risk and high-risk patches on disease levels in the respective patches and the overall environment.  To begin with, it suffices to analyse the disease dynamics in the two patches in the absence of mobility. The patches are thus in isolation and act independently. In this case, Patch 1 basic reproduction number is R 1 0 = 4.0423 and that of Patch 2 is R 2 0 = 0.9908. In order to understand the infection dynamics at the beginning of the infection period and in the long-time behavior (t → ∞), we have plotted the infection curves for t = 20 and t = 1000 in the various plots. Figure 5 shows the evolution of the disease in Patch 1 and Patch 2 in the absence of mobility. We can see that Patch 1 infections in Patch 1 starts increasing at about t = 2 up until t = 17.5 after which the disease remains endemic in the patch. This is an expected result since R 1 0 = 4.0423 > 1. Patch 2 disease dynamics portrays an opposite scenario from those of Patch 1 since its infections curve shows a decline in infections levels from the onset until the disease goes into extinction. This too is expected since R 2 0 = 0.9908 ≤ 1. Figure 14 shows the basic reproduction number computed for different combinations of the infection parameters β 1 and β 2 and in the absence of mobility. The figure reveals that in the absence of mobility, the disease will die out in both patches if approximately β i < 0.6, i = 1, 2. Otherwise, the disease will remain endemic in both patches.
In order to asses the effects of mobility between the two patches on disease dynamics, we elected to combine larger values of α 1 (first fixed at α 1 = 1 and then varying in decreasing order) with varying values of α 2 (increasing in order of magnitude). We combined these with the time that the moving proportion of individuals take in the patch that they have moved to before they go back to their own patch. We consider two instances: a proportion of individuals move and take half (p 12 = p 21 = 0.5) or all (p 12 = p 21 = 1.0) of their time in the patch they have moved to. The disease evolution in the two patches with different mobility patterns are shown in Figures 6, 8, 10 and 12 while the global environment infection dynamics (sum of Patch 1 and Patch 2 infection counts) are given in Figures 7, 9, 11 and 13. The patch specific and global basic reproduction numbers corresponding to the above figures are given in Tables 2 and 3. We make the following observations.  Figure 10 where there is a slight deviation so that there is a low infection level when α 1 = 1.0, α 2 = 1.0 and p 12 = p 21 = 1.0. The high infection levels can be confirmed from Table 2 and Table 3 where we see the global basic reproduction numbers R 0 increasing with the combination of fixed and then decreasing α 1 with increasing α 2 . The figures also reveal that there is an increased level of infection when the proportion of individuals travel and take all their time in the other patch compared to when they take only 50% of their time in the patch they have moved to. This can be visualized from the infection curves in Figure 10 and Figure 12 which are somewhat widely spread out as compared to the infection curves in Figure 6 and   2. As can be seen in Figure 10, infections in Patch 1 generally increases at the beginning the remains endemic with an exception of Patch 1 infection curves when α 1 = 1.0 with α 2 increasing and p 12 = p 21 = 1.0. Whereas Patch 1 infections levels generally increases before stabilizing at endemic levels, infections in Patch 2 generally decreases first then remains endemic. This could be attributed to the fact that Patch 1 is high risk and Patch 2 is low risk so that there is a high contact rate in Patch 1 leading to increased rate of infection at the beginning of the infection period before mobility of individuals (coupled with the high contact rate) maintains the infections at endemic levels. The same argument holds for low risk patch 2 where the low contact rate results into a decrease in the infection levels at the beginning of the infection period: then mobility within the Patch maintains the disease at endemic levels. 3. The disease remains endemic in all the patches. This is expected and is confirmed by the global basic reproduction numbers R 0 given in Table 2 and Table 3 whose values are all greater that 1. Notice that the global basic reproduction numbers given in Table 2 and Table 3 increases both when fixed α 1 = 1 is combined with increasing α 2 values and when decreasing α 1 values are combined with increasing α 2 values. However, the increase of R 0 is larger in the case when individuals move and spend 100%, as opposed to when they spend 50% of their time in the other patch. This is observed from about 0.6 differences of R 0 values down the R 0 columns of Table 2 and Table 3 Figure 7 and Figure 11. This trend is repeated when we combine reducing α 1 with increasing α 2 values whereby we observe that the lowest α 1 = 0.3 combined with the highest α 2 = 0.8 gives the largest R 0 both when p 12 = p 21 = 0.5 and p 12 = p 21 = 1.0 and subsequently, we observe the highest global infection curves as portrayed in Figure 9 and Figure 13. The large increment in R 0 values when p 12 = p 21 = 1.0 have been captured by the corresponding increment in the global infection counts displayed by Figure 11 and most conspicuously, by the large spacing in the global infection curves in Figure  13. 4. Generally, since Patch 1 is higher risk, its infection levels are higher when individuals take 50% of their time in the other Patch as compared to when individuals take all (100%) of their time in the other patch. 5. The local patch specific basic reproduction numbers R i,i 0 in the presence of mobility given in Table 2 and Table 3 reveals that in all instances, R i,i 0 > 1 whenever R 0 > 1 so that the disease remains persistent in Patch i, i = 1, 2. This is in line with Proposition (3.1) and consequently, with the expected theoretical behaviour of reproduction numbers. It is also noteworthy that R i,i 0 increases with increase in R 0 > 1. This may not be clearly visible in the cases when R i,i 0 → +∞ but we noticed that in such instances, the rate of divergence of R i,i 0 increased. This implies that in the instances where R i,i 0 → +∞, the values of r ii > 1 in the NGM were increasing with increase in R 0 . Moreover, it can be noted in the tables that both R i,i 0 and R 0 > 1 increased with increase in mobility from Patch 2 to Patch 1, signifying   increasing levels of contact and infection in Patch 1 (which is a highly risky patch) and consequently in the global environment. The distinct scenario which appears in red in Table 2 is the case where the NGM was reducible. Even in this case, it is noteworthy that the contrapositivity of Proposition (3.2) is fulfilled with the global R 0 being equal to the largest R i,i 0 , which is R 2,2 0 in this case. The relationship between R i,i 0 and R 0 > 1 suggests that in order to have a complete understanding of the patch and global disease dynamics and evolution in the presence of mobility, it is important not to use R i,i 0 in isolation but to also consider the global basic reproduction number R 0 as well. Table 2 Patch specific and global R 0 for different combinations of α 1 = 1 and varying α 2  Table 3 Patch specific and global R 0 for different combinations of varying α 1 and α 2 p 12 = p 21 = 0.5 Notice that the above analysis on the effects of mobility on patch and global environment disease dynamics is not very exhaustive as we have only considered R 0 for specific combinations of the values of α 1 and α 2 for both p 12 = p 21 = 0.5 and p 12 = p 21 = 0.5. In order to have a complete and exhaustive view about the values of R 0 and consequently, the evolution of the disease in the specific patches and the global environment, we elected to compute R 0 from the combination of different values of α 1 and α 2 , both when p 12 = p 21 = 0.5 and p 12 = p 21 = 0.1. These values are plotted as both heatmaps and surface plots in Figure 17 and Figure 18. The surface plots are just to help us view the R 0 values at different angles. Our intention is to help the reader to have a complete and exhaustive panoramic visualization of the R 0 values obtained from the various combination of mobility parameters.
Both graphs in Figure 17 and Figure 18 confirm our assertion above that the disease will remain endemic in both patches regardless of the proportion of individuals moving from Patch 1 to Patch 2 and from Patch 2 to Patch 1 and the time they take in the patches they have moved to. In addition, Figure  17 Figures  18 (a) and (b) reveal that when individuals travel and take all of their time in the patch they have traveled to, a combination of α 1 < 0.6 and α 2 > 0.5 values results into low values of R 0 and hence low infection levels. It is therefore evident that there are high infection levels when individuals travel and take a longer period (100% in this case) of their time in the other patch. This confirms our previous assertion. In fact, we see from Figure 18(a) and (b) that a combination of values of α 1 ≥ 0 with α 2 < 0.4 and α 1 > 0.6 with α 2 > 0.4 results into larger R 0 values and consequently signifying high infection levels when p 12 = p 21 = 1.0. We thus have a larger combination of mobility patterns that results into high infection levels when individuals travel and take all of their time in the patch they have moved to than when they take half of their time in the patches they are visiting. Simply put: there is a larger red area in both graphs in Figure 18 than in both graphs of Figure 17.    the patches they have to. We can observe from both figures that when α 1 is decreasing and α 2 is increasing, we obtain lower R 0 values with a combination of increasing β 1 and decreasing β 2 (up-to some constant) values. The graphs reveals that the line separating the larger and smaller R 0 values is non-linear in the case when individuals travel and take half of their time in the other patch and is linear when the travelers take all of their time in the patches they are visiting. Generally, in both cases, a combination of small β 1 (< 0.7) and β 2 (< 0.5) values should give R 0 < 1 values thereby leading to low infection rates in the specific patches and hence the global environment.  Figure 15 shows surface plots of R 0 values when varying α 1 and α 2 values while fixing p 12 = 0.1, p 21 = 0.9 (Figure 15 (a)) and fixing p 12 = 0.9, p 21 = 0.1 (Figure 15 (b)). We can see from Figure  15 (a) that R 0 values are high when a small proportion of Patch 1 individuals move and spend a small proportion (10%) of their time in Patch 2 and a small proportion of Patch 2 individuals move and spend a large proportion (90%) of their time in Patch 1. This combination of mobility scenarios thus results into a high infection levels in the high-risk Patch 1 and subsequently in the global environment. The values of R 0 decreases, but not as fast, with an increase in the proportion of Patch 1 individuals moving to (and spending 10% of their time in ) Patch 2 and a small proportion of (0%) of individuals moving from Patch 2 to (and spending 90% of their time in) high-risk Patch 1. Figure 15 (a) further reveals that R 0 decreases with a combination of increasing values of α 1 and increasing values of α 2 , up-to α 2 = 0.6, where we achieve the lowest value of R 0 which starts to increase for values of α 2 > 0.6. This reveals that an increasing proportion of Patch 1 individuals moving and spending 10% of their time in Patch 2 and an increasing proportion (up-to 60%) of Patch 2 population moving and spending 90% of their time in Patch 1 results into decreasing infection levels in the specific patches (especially Patch 2) and subsequently in the global environment. It is noteworthy from Figure 15 (a) that R 0 > 1 regardless of any mobility combination patterns. The disease therefore remains endemic in both Patch 1 and Patch 2 and subsequently in the global environment. A combination of large proportion of Patch 1 population moving and spending 10% of their time in Patch 2 and large proportion of Patch 2 population moving and spending 90% of their time in Patch 1 results into high global R 0 values and thus infection levels. However, these R 0 values (and thus infection levels) are not as high as when a low proportion of Patch 1 population moves and to (and spends 10% of their time in) Patch 2 and a low proportion of Patch 2 population moves to (and spends 90% of their time in) Patch 1. Figure 15 (b) shows that a combination of a low proportion of Patch 1 population who move to (and spend 90% of their time in) Patch 2 with whichever proportion of Patch 2 individuals who move to (and spend 10% of their time in) Patch 1 results into the largest values of R 0 and consequently, the highest levels of infections in Patch 1 and the global environment. Put differently, we see that the global R 0 decreases for the combination of increasing values of α 1 and any values of α 2 . We recall that Patch 1 is a high risk patch and hence having more people in Patch 1 will result into high infections at the beginning of the infection period before the mobility and the amount of time taken within the Patch maintains the infections levels and endemic levels. Notice that as in Figure 15 (a), Figure 15 (b) shows that R 0 > 1 for any combinations of the values of α 1 and α 2 . The disease thus remains endemic in both Patches and in the global environment. However, contrary to Figure 15 (a) which shows that the lowest R 0 = 3.2, the mobility scenarios captured by Figure 15    We only need to realize that whereas in Figure 15 we are varying α 1 and α 2 and fixing p 12 = 0.1, p 21 = 0.9 (Figure 15 (a)) and p 12 = 0.9, p 21 = 0.1 (Figure 15 (b)), in Figure 16, we are varying p 12 and p 21 and fixing α 1 = 0.1, α 2 = 0.9 (Figure 16 (a)) and α 1 = 0.9, α 2 = 0.1 (Figure 16 (b)). This implies that we can either fix p 12 = p 1 and p 21 = p 2 and vary α 1 and α 2 or fix α 1 = p 1 and α 2 = p 2 and vary p 12 and p 21 and obtain similar surface plots with similar interpretations.

Large scale simulation studies
Our interest is to compare if there are any structural differences in the forward simulations obtained by using randomly generated residence matrix and mobility parameters with estimated ones from real data. This comparison is important as it validates the performance and the use of the proposed model in a real rather than theoretical world. For this purpose, we consider forward map simulations for patches of sizes n ∈ {103, 203, 303, 403, 503}. The   residence-time matrix P and the mobility parameter α for each of the n sized patches were randomly simulated.
As an application to the real world situation, the residence matrix P and the mobility parameter α for the AGEBs in Hermosillo, Sonora Mexico were estimated using the Brownian bridge technique. We refer the reader to the work by [43] in order to have an idea of how this estimation was carried out. The total number of AGEBs (which we consider as patches in this analysis) for which  these parameters were estimated were n = 503. This explains why we have considered forward simulations for patch sizes of n ∈ {103, 203, 303, 403, 503}.
Our intention is to see if any structural differences exists between the forward simulations obtained using randomly generated residence-time matrix and mobility parameters for a patch of size n = 503 and using the estimated residence-time and mobility parameter values from the AGEBs in Hermosillo.
In order to achieve our comparison objective, we only give a plot of the epidemiological curves for each patch for n = 503. Figure 21 (a) shows the results from randomly generated residence-time matrix while Figure 21   shows the same curves obtained by using estimated P and α from the AGEBs in Hermosillo, Sonora, Mexico. A comparison of the two graphs reveals that there are structural differences between the respective curves. That is, the respective  curves do not have structural uniformity with those generated by estimated P and α from the AGEBs being somewhat non-homogeneous. Of importance to note from Table 4 is that the global basic reproduction numbers for the random P and α for n = 503 was R 0 = 7.2618 while that for the estimated P and α from the AGEBs was R 0 = 9.7775. As such, we have the disease remaining endemic in all the patches and the global environment (since R 0 > 1) for both cases. This can respectively be verified from the infected curves in Figure 21 (a) and Figure 21 (b). Notice the large difference in the global basic reproduction numbers for the two cases. This confirms the heterogeneity and variability of the epidemiological curves as can be seen in Figure 21 (b) compared with the somewhat homogeneous and nob-variable epidemiological curves in Figure 21 (a). Figure 22 shows a comparison of the global curves (sum of the susceptible, exposed, infected and recovered for each patch) for varying patch sizes, the AGEBs and a plot of the single patch SEIRS model. Figure 21 are plots of the global curves for t = 50 and t = 150 respectively. The graphs reveal that the infections in the global environment peaks (save for the single SEIRS model) at t = 25, slightly reduces and then remains endemic. The endemicity of the infections in the global environment is due to the fact that R 0 > 1 in the considered cases as can be seen in Table 4. The slight structural difference in the curves obtained by randomly simulated and estimated P and α has been confirmed in Figure 21. It can be seen that the global curves corresponding to these scenarios are somewhat comparably different with only the susceptible global curve obtained from the estimated parameters being slightly above the global susceptible curve obtained from the randomly generated ones. The rest of the global exposed, infected and recovered curves obtained from the estimated P and α are slightly below the same global curves obtained from the randomly generated P and α. There are higher counts as n increases. This is expected since the global curves are obtained from the addition of the individual patch points, which should be high when n is large. Theorem 3.1 holds in this large scale simulation case as it was confirmed that all the global basic reproduction numbers for all n ∈ {103, 203, 303, 403, 503} as in Table 4 are bounded below and above respectively by the minimum and maximum Patch i, i = 1, 2, · · · , n basic reproduction numbers R i 0 when the patches are acting as isolated units, that is, when the patches are not connected through mobility.  Fig. 22 Global curves from the AGEBS in Hermosillo, the varying n-sized patches obtained from the randomly generated residence matrix/mobility parameters and the single patch model.

The generalized Multi-Patch SEIRS model
The SEIRS model that we have analyzed above is formulated on the premise that Patch i i = 1, 2, · · · n population is moving as a whole group without any discrimination and separation of the moving population into their specific epidemiological or clinical groups. In this section, we derive a generalization of the already analyzed SEIRS model where we assume that a proportion of each Patch i epidemiological population groups S i , E i , I i and R i move and spend a proportion of their time in Patch j, j = 1, 2, · · · n. This is to say that we now discriminate and capture the effects of the mobility of each of the population belonging to different clinical statuses on Patch j, j = 1, 2, · · · n disease dynamics. We therefore introduce and capture heterogeneity in the mobility of Patch i different epidemiological groups.

Derivation of the model
Let us denote S i , E i , I i , R i as the number of susceptible, exposed, infected and recovered people in patch i. Let α S i , α E i , α I i and α R i respectively be the proportion of susceptible, exposed, infected and recovered population moving out of Patch i into Patch j. Further, let s ij , e ij , i ij and r ij be the proportion of time spent respectively by Patch i susceptible, exposed, infected and recovered population in Patch j. The defined proportions of time should be such that As such, we assign different factors (proportion of time) for different subgroup thereby accounting for heterogeneity. These factors consider that only a proportion of the population move out from their own community. Further, by assigning different factors to each subgroup, heterogeneity in their movement is captured. Based on the above definitions, we observe that there are patch i residents in patch j on average at time t and that there are patch i residents that do not move to patch j. The number of infected people in patch j is where the first term represents the total infected people moving into patch j, the second term represents the total infected population that stays at patch j. The second term can be further simplified using the identity n j=1 i ij = 1. Therefore, the effective infected people in patch j is given as: The effective population in patch j is therefore given as: where the second equality is obtained after using the identities in Equation (42) on the first equality. The first part on of the RHS of the second equality represents the remaining population in patch j, whereas the second part represents people coming into j. Therefore, proportion of infected individuals in patch j is: Number of susceptible from patch i who are currently in patch j is α S i S i s ij and hence we have the S i population infected per unit time in patch j as The effective population in Patch i (the sum of those who remain in Patch i and those who visit from Patch j) is: Hence the effective proportion of infected individuals in patch i is: But notice that the total S i population that remain in patch i is: Therefore, S i population infected per unit time in patch i is given as: Therefore, total S i infected per unit time based on Equation (43) and Equation (44) is: Now incorporating the dynamics of SEIRS model, we get the following ODE system: where Z i and Z j are as defined previously. The definitions of the parameters used in the above model are in Table 5.  Per capita natural death rate in Patch i γ i Per capita recovery rate of individuals in Patch i τ i Per capita loss of immunity rate ψ i Per capita disease induced death rate of Patch i κ i Per capita rate at which the exposed individuals in patch i becomes infectious The above system can also be written in matrix form as follows:                   Ṡ = Λ − diag(S)s ⋆ diag(β)diag −1 (s ⋆T S + e ⋆T E + i ⋆T I + r ⋆T R)i * t I −diag(µ)S + diag(τ )Ṙ E = diag(S)s ⋆ diag(β)diag −1 (s ⋆T S + e ⋆T E + i ⋆T I + r ⋆T R)i * t I −diag(κ + µ)Ė I = diag(κ)E − diag(γ + ψ + µ)İ R = diag(γ)I − diag(τ + µ)R where,

The basic reproduction number and mathematical analysis of the model
Using the next generation approach, we obtain the next generation matrix of system (46) as Hdiag(κ)diag −1 ((κ + µ) ⊙ (γ + ψ + µ)) Hdiag −1 (γ + ψ + µ) 0 0 where H = diag(N )s ⋆ diag(β)diag −1 (s ⋆T N )i ⋆T From which we obtain the basic reproduction number as R Ge 0 = ρ(Hdiag(κ)diag −1 ((κ + µ) ⊙ (γ + ψ + µ))) The local Lipschitz continuity and positive boundedness of system (45) can be shown using a similar procedure for the system (3). Theoretical and simulation studies on the effects of heterogeneity in mobility and patchiness on the basic reproduction number R Ge 0 is still an ongoing work. However, analysis of the global stability of equilibria of system (46) in its current form is quite challenging. In fact, such models exhibits intricate non-linearity that may lead to the existence of multiple endemic equilibria (see, e.g., [44]). It may therefore be meaningful to explore an analysis of global stability of equilibria for some particular cases of the generalized system (46). In this sense, we consider a particular case where we ignore the epidemiological-statuses-dependent mobility. That is, we assume that s ⋆ = e ⋆ = i ⋆ = r ⋆ . Further, as in the previous global stability analysis, we assume that there is no disease induced mortality, in which case, the total population dynamics is given aṡ N = Λ − µ ⊙ N and hence we have that lim t→∞ N = Λ ⊙ 1 µ = N . Consequently, an application of a mix of the above two assumptions and the theory of asymptotically autonomous systems for triangular systems as discussed by [36] and [45], the system (46) is asymptotically equivalent to               Ṡ = Λ − diag(S)s ⋆ diag(β)diag −1 (s ⋆T N )i * t I − diag(µ)S +diag(τ )Ṙ E = diag(S)s ⋆ diag(β)diag −1 (s ⋆T N )i * t I − diag(κ + µ)Ė I = diag(κ)E − diag(γ + ψ + µ)İ R = diag(γ)I − diag(τ + µ)R We denote the basic reproduction number of system (47) as R GA 0 and obtain its expression, using the next generation approach, as R GA 0 = ρ(diag(N )s ⋆ diag(β)diag −1 (s ⋆T N )i * t diag(κ)diag −1 ((κ + µ) ⊙ (γ + µ))).
Theorem (6.1) gives the global stability of the DFE and the existence of a unique EE of system (47). The proof of the first and second parts of Theorem (6.1) is respectively similar to that of Theorem (4.1), Theorem (4.2) and Theorem (4.3).

Conclusion
We have generalized the existing multi-patch SEIRS epidemic model [1] in several important directions. First, we introduce the proportion α i of Patch i individuals who move and the time p ij that the moving Patch i individuals take in Patch j. Second, we have performed a quite complete analysis of the proposed model. In particular, we have shown that the model is locally Lipschitz and is biologically well posed. Under the proposed formulation, we compute the global and patch specific basic reproduction numbers R 0 as a function of the mobility parameters α i and the residence time matrix elements p ij . We perform global analysis where we have ascertained that in the presence of mobility, the disease dynamics purely depends on the global basic reproduction number R 0 while the patch specific reproduction numbers R i 0 controls Patch i disease dynamics in the absence of mobility. We have shown that the DFE is GAS if R 0 ≤ 1 and that there exists a unique EE which is asymptotically stable when R 0 > 1. Third, a novel definition of the patch reproduction number R i,i 0 has also been introduced where we have shown, theoretically and through numerical simulation, that there is a close relationship between R 0 and R i,i 0 as regards to the persistence and extinction of the disease. Our newly defined patch reproduction numbers obeys the expected theoretical behaviour of reproduction numbers in this regard as is discussed in the text. Fourth, we have also introduced and analyzed a more generalized epidemiological-statusmobility-dependent SEIRS epidemic model which encapsulates the idea that Patch i disease levels can be determined when a proportion of the population in the respective Patch Susceptible, Exposed, Infected and Recovered epidemiological groups travel and spend a proportion of their time in the other patch.
We have carefully formulated and run a simple two-patch forward map simulation in order to asses the role of the mobility parameter α, the residence matrix P, and the environmental risk parameters β i playing in the evolution of the disease in the two patches. In our simulation study, we have shown that the synergy between these parameters can lead to the disease remaining endemic or going into extinction in the individual patches. For example, we have ascertained that a larger proportion α i of Patch i individuals moving and spending more of their time p ij in a sufficiently high risk Patch j (β j ≫ β i ) may lead to the disease remaining endemic in Patch j. The simulation results reveal that there are high infection levels in a high risk patch regardless of the proportion of individuals moving and the proportion of time those moving spend in low risk patch. Our formulation and analysis has revealed that heterogeneity of mixing distributions in terms of mobility and residence times are vital on the spread of infectious diseases.
The performance of the proposed model in the face of a large scale forward simulation scenario has also been considered in this research. We have combined this large scale simulation analysis with an application of the model in a real world scenario where we have used an estimated residence-time matrix P and mobility parameter α from AGEBs in the city of Hermosillo, the state of Sonora, Mexico. A comparison of the forward map simulations obtained from use of randomly generated and estimated P and α reveals some structural differences and variability. Consequently, real data should be incorporated to accurately simulate and forecast of an ongoing infectious disease.