Mathematical Modelling of Multi-Region Spread of Epidemic

SAIR model of growth of an epidemic is extended to a system of many interacting regions. Interactions are described by exchange of populations between various regions. Diﬀerences caused by the exchange of susceptible population, in addition to infected individuals are noted. It is shown that initial phase of the epidemic is governed by a linear system. Analysis of linear system shows that a fundamental mode gets established that governs the spatial proﬁle of the spread of the disease.


INTRODUCTION
Covid-19 pandemic has given rise to renewed interest in the Mathematical Modelling of spread of a contagious disease, Yang and Wang [2020], Lobato et.al.
[2021], Agrawal et.al. [2021], Peter et.al. [2021]. This is an old topic, dating back to the work of Kermack and Meckendrick [1927] who proposed the SIR model. In this model the population Π of a region under study is divided in three groups, namely Susceptible (S), Infected (I) and Removed (R). It is postulated that the disease spreads through a contact between a susceptible person and an infected individual with the probability β ′ dt Π in an infinitesimal time interval dt. Thus the number of people getting infected in a time dt is given by β ′ SIdt Π , where S and I denote the number of susceptible and infected individuals at time t. The infected individuals are transferred to the Removed (R) group at the rate γIdt, either through recovery or death. One then sets up non-linear, first order, ordinary differential equations (ODE) describing how people go over from one group to another. These equations merely state the conservation of the total number of people in three groups. Their non-linearity makes it difficult to analyse them theoretically and hence are generally studied numerically. In what follows we will set β = β ′ /Π which is equivalent to the normalization Π = 1.
This basic model has been improved by various authors by including additional groups to account for the specific characteristics of a particular disease. Thus for example in the SEIR model, Li and Muldowney [1995], one introduces another group of Exposed (E) people On contacting the disease a susceptible individual moves to Exposed (E) group for a certain incubation period before moving to the Infected (I) group. The SAIR model was introduced by Robinson and Stilianakis [2013] to account for Asymptomatic (A) individuals. Susceptible individuals can contact disease by a contact with either the Asymptomatic person or an infected one and join the Asymptomatic group. Asymptomatic individuals either recover from the infection and directly move to Removed (R) group or move to Infected (I) group. The Infected persons move to the removed (R) group either through recovery or death. Once again one can set up ODE describing conservation of people in these four groups.
Most of the studies consider only one region, which can be a town, city, state or a country. Few authors have considered the spatial and temporal development of the epidemic which is the subject of this paper. Thus e.g. Noble [1974], added diffusion terms to the balance equations for both Susceptible and Infected persons of basic SIR model. This converts these two equations to time-dependent partial differential equations (PDE) Recently Besse and Faye [2021] used diffusion equation to account for migration of (only) infected individuals on connected graphs, a system of cities connected by a transportation network. This converts the infected (I) population equation for each node (city) to a non-linear PDE, coupling the neighbouring parts of the transportation network The equations for Susceptible (S) and Removed (R) groups remain the non-linear ODE. A much simpler approach was followed by Zakary, Richik and Elmouki [2017] who considered a discrete time evolution of Multi-region SIR model, allowing for the migration of only infected individuals. We follow this approach but over continuous time, using classical mathematical physics methods. Further we will examine the consequences of allowing movement of all, the Susceptible, Asymptomatic carriers and Infected people to different regions, i.e. with SAIR model.
In the next section we first briefly recall the basic SAIR model and then set out its generalization to multi-region case. We write down the equations allowing for migration of (i) Susceptibles Asymptomatic and Infected persons, (ii) of Asymptomatic and Infected individuals only and note the differences in the resulting equations. This leads to a set of 2N coupled. non-linear equations for the susceptible and asymptomatic populations of N regions. Remaining 2N equations governing the populations of infected and removed groups of N regions are linear equations where the regions are also decoupled. Then in section 3 we observe that in the initial stages of the epidemic the infected populations, asymptomatic or those showing symptoms are much smaller than the susceptible populations. This allows us to analyse a linearised model and gain some insight.
Our aim in this paper is to proceed analytically as far as possible. We obtain the modes of this linearized model and show that a dominant fundamental mode exists with all non-negative elements. This fundamental mode grows much more rapidly than other modes. We then set up the development of the solution of non-linear equations by a perturbation series expansion. In section 4, we illustrate our theoretical conclusions by numerical computations. We consider a system of just two regions when all migration from one region is to the only other region. This results in a particularly simple form of equations that can be easily solved. It is seen that model predicts all the features of the phenomenon, in quantitative terms, which one feels intuitively. Lastly in section 5 we state our conclusions.

Basic Equations
As mentioned in the Introduction, the SAIR model divides the population of a region i in four groups, namely, Susceptible (S), Asymptomatic (A), Infected (I) and Removed (R). Let these symbols also denote the number of individuals in the region. One can easily write down the equations governing the time evolution of these four groups as Here β A , β I denote the rate at which an Asymptomatic or Infected individual causes infection in a Susceptible person. γ A denotes the rate at which an asymptomatic person recovers from the disease and moves directly to the removed group while α measures the rate at which asymptomatic persons move to infected group after showing symptoms of the disease. γ I measure the rate of removal of infected persons by recovery (fraction µ) and death with fraction (1 − µ). For simplicity we will consider the simplified SAIR model when β A = β I = β and γ A = γ I = γ though this assumption can be relaxed without any additional difficulties. Eq. (1) then reduce to We observe that first two equations are non-linear and coupled while the remaining three are linear and are essentially driven by a "source" αA. We now describe the extension of Eq. (2) to multi-region case.

Model (a) Movement of Susceptible, Asymptomatic and Infected groups
Let us now consider a collection of N regions. Let us assume that a fraction ξ i,j of the population of the region j is visiting the region i at any given time.
We will assume that this fraction is time-independent. For simplicity we will also assume that same fraction is applicable to susceptible, asymptomatic and infected groups, though this assumption can be easily relaxed. We also denote ξ i = N j=1,j̸ =i ξ j,i . Let S i , A i , I i and R i denote the susceptible, asymptomatic, infected and removed populations in the i-th region, i = 1, 2, .......N . Let us also for the moment assume that the people who show symptoms of the disease i.e. the infected individuals, isolate themselves and avoid further contacts. Then the time evolution of the S i is governed by the equation Eq.
(3) states that the rate of change of susceptible population in the region i consists of two parts. A fraction (1 − ξ i )S i catches infection locally from the asymptomatic people of the same region whose number is (1 − ξ i )A i and rate β i,i . Thus this rate is (1 − ξ i ) 2 β i,i A i S i . Another set of susceptible is affected by the fraction ξ i,j of asymptomatic people of region j who were visiting the region i, and that totals to β i,j ξ i,j A j (1−ξ i )S i . The second term accounts for the fraction ξ j,i of susceptible of region i contacting infection in some other region j, from the asymptomatic carriers of that region (rate β i,j ξ j,i A j (1 − ξ j )S i ) as well as from the asymptomatic individuals visiting region j, including those from region i. We Separate the contribution from region i, set ζ 2 i = N j=1,j̸ =i ξ 2 j,i and define ξ k,i η i,k = N j=1,j̸ =i,j̸ =k ξ j,i ξ j,k . With this notation we can write Eq. (3) as This is an extension of a slightly modified form of SAIR model. A straight forward extension of SAIR model of Eq. (2) wherein asymptomatic carriers and infected individuals both come in contact of susceptible population can also be worked out and is given below. As we will see there is no material difference between these two versions of generalised SAIR models.
The growth of Asymptomatic people A i of region i is given by It is assumed that on contacting the infection a person moves initially to the asymptomatic group. A fraction recovers at the rate γ i A i while another fraction moves to the Infected category I i at the rate α i A i . Hence other three equations of the set, Eq. (2), are generalised to These are three linear equations for the region i, de-coupled from other regions and driven by a source term α i A i . It is seen that Eqs. (4) and (5) are a set of 2N coupled, non-linear equations and one cannot treat region i in isolation. It is not very hard to numerically solve Eqs. (4), (5) and (6) if the number of regions N is not too large. Fractions ξ i,j depend upon the geographical boundaries of the regions i, j and the number of people of region j that regularly visit region i for work, study and other purposes. This data can be inferred from normal monitoring of daily commute of people.
A straight forward generalization of Eq.(2) is as follows. Define M i = A i + I i , i = 1, 1, ...N . Using same arguments as before we arrive at the following equations in place of Eqs. (4), (5) and (6) and Remarks made for Eqs. (4), (5) and (6) are also applicable to Eqs. (7), (8) and (9). Eqs. (7) and (8) are a set of 2N coupled set of non-linear, first order equations while Eq. (9), a set of three linear first order ordinary differential equations determining I i , R i and D i for the region i, isolated from other regions.
Since the coupling of different regions is through Eqs. (7) and (8), we will concentrate only on these two equations. We also note that the treatment of Eqs. (4) and (5) is practically the same. It merely needs a redefinition of the parameters β i , γ i .

Model (b) Movements of only Asymptomatic and Infected groups
Many authors study the effect of movement of infected groups (Asymptomatic and Symptomatic) only. The susceptible population remains confined to their home region. With the notation introduced earlier this model leads to the following differential equations in place of Eqs. (7) and (8).
Eq. (9) remains unchanged and is written here for completeness For completeness let us also write down the solutions of Eq. (12). Thus we have and Eqs. (12), (13), (14) and (15) are applicable to both models (a) and (b). We thus have to study Eqs. (7) and (8) (or equivalently Eqs. (10) and (11)) and obtain M i (t), i = 1, 2, ...N and we proceed to do that in the next section. We may add here that although we consider Eqs. (7) and (8), most of the analysis is also applicable to Eqs. (10) and (11). Any difference will be pointed out as we go along.

Linearization and Solution of Eqs. (7) and (8) by Perturbation Expansion
In the initial stages of the epidemic, the number of infected individuals M i (t) is small in all regions. The Susceptible populations S i is close to their initial fixed values. We can therefore assume S i = S i,0 + λs i (t) where the change s i (t) is small while S i,0 are constants independent of time. We have introduced a perturbation parameter λ. We expand both M i and s i as power series in λ and write Eqs. (7) and (8) are modified to Substituting the expansions, Eq. (16), in Eqs. (17) and (18) and equating the coefficients of λ n on both sides of equation we have for n = 0, 1, 2, .... and (20) In the zeroth order approximation is a set of N coupled, first order, linear differential equations with constant coefficients. It can be solved easily by reducing it to a set of decoupled equations which we will outline below. Eq. (23) does not have the functions s let us now consider Eq. (22) and cast it in a matrix form. Let us define a Matrix H with matrix elements Time independent N × N , real square matrix H has N eigenvalues ν 1 , ν 2 , ....ν N , real or complex, repeated or distinct. Complex eigenvalues occur in conjugate pairs and lead to oscillatory solutions. We will assume that it has N eigenvectors. Let P be an N × N matrix whose columns are the eigenvectors of H. The case of H having fewer eigenvectors can also be treated by including generalized eigenvectors. That will slightly modify the treatment given below. Pre-multiplying Eq. (25) by the time independent matrix P −1 we get where Λ is a diagonal matrix with entries ν 1 , ν 2 , ....ν N along the main diagonal. The differential equations in Eq. (26) are easily solved and we have where e Λt is a diagonal matrix with entries e ν1t , e ν2t , .....e ν N t along the main diagonal. Thus we have x (0) (0) is the column vector of known initial values of M i , i = 1, 2, ...N at t = 0. Eq. (28) describes the time evolution of infections as a superposition of eigenvectors of the matrix H. Each eigenvector evolves at the rate e νit , with its eigenvalue ν i . If the real part ℜ(ν i ) of the eigenvalue is negative the eigenvector decays with time or grows if ℜ(ν i ) > 0. Consequences of these are examined in next subsection. We observe from Eq. (23) that to compute s dt. These can be found easily when we notice where I is the identity matrix. We now consider the evaluation of higher terms in the perturbation series, i.e evaluation of n − th order terms M i (t) of previous orders k = 0, 1, 2, ...n − 1. An expression for b (n) (t) can be written down as follows. Let S (k) (t) be a sequence of N × N diagonal matrices with diagonal elements (s Thus the vector b (n) (t) can be computed if the vectors x (n−1−k) (t) and matrices S (k) (t) of earlier terms in the series expansion are known. For the moment we will assume that these are known, though we have to devise methods for computing diagonal matrices S (k) (t), k = 0, 1, 2, ..n − 1 Pre-multiplying Eq. (30) by the matrix P −1 we get and the solution of Eq. (32) is given by (33) We now turn our attention to the computation of the matrices S (k) (t), k = 0, 1, 2, ..n−1. Let s (k) (t), denote N dimensional vector with components (s  It is clear that this vector immediately yields the diagonal matrix S (k) (t). We already have an expression for the vector s (0) from Eqs. (23) and (29). Thus we have . Hence we have Substituting for x (n) from Eq. (33) we have We can now successively compute M (n) i.e. all terms in perturbation expansion. We form the matrix H, given by Eq. (24), and find its eigenvalues (ν 1 , ν 2 , ...., ν N ) and the corresponding eigenvectors and obtain the matrices Λ, P and P −1 . as defined in Eq. (26). We can then obtain the vector x (0) (t) from Eq. (28) and known initial values of M i , i = 1, 2, ...N at t = 0. We also obtain the vector s (0) (t) from Eq. (34) and then compute the diagonal matrix S (0) . We can then compute the vectors b (1) , x (1) and s (1) from Eqs. (31), (33) and (37) for the case of n = 1.. Knowing these vectors for n = 1 we can compute the corresponding vectors for n = 2 and so on.

Model (b)
The analysis given above is also applicable to Model (b) wherein Susceptible population is confined to their home region i.e. the starting equations are Eqs. (10) and (11). For perturbation expansions, Eqs. (17) and (18) are replaced by Using the expansion, Eq. (16), in these equations we obtain the following equations for s for n = 0, 1, 2, .... and Matrix formulation presented earlier (for zeroth and higher order approximations) is applicable to Eqs. (40) and (41) if we replace the matrix H with the matrix H b whose matrix elements are given by Here δ i,j is Kronecker's δ function, being unity when i = j and zero otherwise. Thus for example Eqs. (25) and (26) are modified to where Λ b and P b are the matrices of eigenvalues and eigenvectors of the matrix H b respectively. Rest of the analysis, Eqs. (27) to (41) is same as before after replacing the matrices H, Λ, P and P −1 with H b , Λ b , P b and P −1 b .

The Fundamental Mode
We now examine the eigenvalue problem of matrices H and H b more closely as it is clear from above analysis that the eigenvalues of these matrices play a crucial role in the growth of the epidemic.
Our first observation is that if the fractions ξ i,j are symmetric i.e. ξ i,j = ξ j,i then the matrices H and H b are real symmetric matrices. All their eigenvalues are real and they have N distinct eigenvectors. Thus P and P b exist and are orthogonal matrices and P −1 = P T , the transpose of P. Similarly Since all the eigenvalues are real, one of them is the largest. Corresponding eigenvector constitutes the fundamental mode if all its elements are non-negative. This mode grows much faster than other eigenvectors and dominates after a few initial time steps. We will assume that the eigenvalues are ordered i.e. ν 1 > ν 2 ≥ ν3, .....ν N In case some of the eigenvalues are complex then their ordering is according to their real part.
We also observe that all fractions ξ i,j ≥ 0. Thus the off-diagonal elements of H and H b are non-negative. If in addition the diagonal elements Let us assume that H and H b are Irreducible. Then we have from Perron-Frobenius theorem that these matrices have a largest positive eigenvalue, with one and only one corresponding eigenvector whose elements are non-negative. This eigenvector is the dominant fundamental mode which grows faster than all other modes. We note that Irreducible, non-negativity, Eq. (45), is a sufficient condition for the existence of fundamental mode. It is not necessary. Thus e.g. if all the parameters γ i = γ, i = 1, 2, ...N then also fundamental mode exists. This follows from the fact that in this situation Thus the matrices G and G b have a largest positive eigenvalue and a corresponding non-negative eigenvector. The eigenvectors of G, G b are also the eigenvectors of H, H b . Thus H, H b have the same eigenvector which is their fundamental mode. It is easy to find this fundamental mode in some simple cases. Let us assume that all coefficients β i,j = β are equal, γ i = γ and all the initial Susceptible populations S i,0 = S 0 are also equal. Further we assume that the fractions ξ i,j = ξ j,i are symmetric. Let z = (z 1 , z 2 , ...z N ) denote an eigenvector of H corresponding to an eigenvalue ν which is explicitly written We also note that under the conditions stated above, before Eq. (47), the matrix H b also has same dominant eigenvector and same eigenvalue as H. Remaining eigenvalues and eigenvectors, however, are different.
Other eigenvalues of H and H b are also of interest. They determine the rate at which the fundamental mode gets established. As noted earlier, if ℜ(ν i ) < 0, i = 1, 2, ..N then all these modes decay while if ℜ(ν i ) > 0 for some i, then that mode also grows. In that case the eigenvalue separation ν 1 − ℜ(ν i ), between the fundamental eigenvalue and the particular eigenvalue is an important parameter. This parameter determines if the fundamental mode gets established before the non-linear effects become significant.
Before leaving this section we observe that main advantage of the above analysis is theoretical. By analysing initial phase of the epidemic we can identify some characteristics of the solutions that can help us to decide which model fits the, observed data better. Thus e.g. we can find the influence of allowing movement of susceptible population in addition to that of infected people. We also expect that the relative proportion of infections in different regions will be similar to the fundamental mode, even when non-linear effects become important, whatever be the initial distribution. In fact powerful computer codes are available that yield accurate numerical solutions of the coupled differential equations. The role of perturbation methods for obtaining numerical results is limited.

Two Region Problem
We now apply above considerations to a simple two region problem. Since there are only two regions, people going out of region 1 are only going to region 2. Thus ξ 2,1 = ξ 1 , ξ 1,2 = ξ 2 . It is seen that in this case η 1,2 = η 2,1 = 0 and ζ 2 1 = ξ 2 1 , ζ 2 2 = ξ 2 2 . To simplify the problem still further we will assume β 1,1 = β 1,2 = β 2,2 = β and γ 1 = γ 2 = γ. We will also assume that the initial susceptible populations in two regions are equal i.e. S 1,0 = S 2,0 = S 0 . Eqs. (7) and (8) then take a simple form Matrix H is a real symmetric matrix. It has two real eigenvalues ν 1 , ν 2 given by the equation where the coefficients a, b, c are identified as the elements of the matrix H + Γ in Eq. (51). H has two real eigenvectors and the orthogonal matrix P (of eigenvectors) and its inverse are given by Here p 0 , p 1 are the components of the normalised fundamental eigenvector. The eigenvalue ν 1 > ν 2 and the corresponding eigenvector is the fundamental mode. In general the eigenvalues ν 1 , ν 2 and corresponding eigenvectors depend on the parameters ξ 1 , ξ 2 . These also conform to the transformation mentioned above. However in case ξ 1 = ξ 2 = ξ we have ν 1 = βS 0 − γ, independent of ξ 1 , ξ 2 while ν 2 = βS 0 (1 − 4ξ + 4ξ 2 ) − γ. In this case we have c = a and r = 2b. The matrices P and P −1 then reduce to Another interesting case is when ξ 2 = (1 − ξ 1 ). In this case the eigenvalue ν 1 = 2βS 0 [ξ 2 1 + ξ 2 2 ] − γ > βS 0 − γ and ν 2 = −γ. The matrix P of eigenvectors and its inverse P −1 are again given by Eq. (54). Thus the fundamental mode is same as that for the case ξ 1 = ξ 2 but is attained earlier in time as the separation of eigenvalues is larger.
If we consider the model (b) then we have the equations in place of Eqs. (49) and (50). Linearised treatment deals with the matrix H b which in this case reduces to The matrix H b is not symmetric. It has two eigenvalues ν 3 = βS 0 − γ and ν 4 = βS 0 [1 − (ξ 1 + ξ 2 )] − γ. Its matrix of eigenvectors P b is not orthogonal. P b and P b −1 are given by Clearly the first eigenvalue ν 3 > ν 4 and the vector (ξ 2 , ξ 1 ) T is the fundamental mode. It is interesting to note that the eigenvalue ν 3 is independent of ξ 1 , ξ 2 , but it is the second eigenvector that does not depend on these parameters.

Numerical Solutions
We solved non-linear equations, Eq. (49) and (50) We varied the parameters ξ 1 , ξ 2 over a fairly wide range with values of 0.1, 0.3, 0.5, 0.7 and 0.9, considered three values of β = 0.15, 0.2, 0.25 but kept γ = 0.1. We observe that the roles of two regions 1 and 2 are interchangeable. We note that the solution of initial value problems governed by differential equations is quite sensitive to the number of digits retained in computations. We observed that retaining six digits in computations gives large errors in the number of infected persons. If eight digits are retained in all computations, the errors are still significant. Results reported in tables 4 -8 were obtained with retaining 16 digits in all computations.
Similarly ν 3 > 0 and ν 4 is generally negative. Hence only the fundamental mode grows with time while the second eigenvector decays. Further, ν 1 ≥ ν 3 This implies that the movement of Susceptible population leads to a more rapid growth of disease, compared to the case of movement of infected people (asymptomatic or otherwise) only. We also note that even a slight increase in ν 1 (ν 3 does not vary with ξ 1 , ξ 2 ) leads to a substantial increase in the number of infections because of the exponential behaviour.
Eigenvalue separation (ν 1 −ν 2 ) or (ν 3 −ν 4 ) increases with βS 0 , being directly proportional to this parameter. This suggests that fundamental mode is established earlier for higher values of β. Variation of these two separations with ξ 1 , ξ 2 is mixed. Sometimes we have (ν 1 − ν 2 ) > (ν 3 − ν 4 ) while at other times reverse is true. Thus in some cases fundamental mode is attained earlier for model (a) while in other cases model (b) reaches this equilibrium first. Numerical computations confirm these trends though non-linear effects, due to reduction of S 1 (t), S 2 (t), S 3 (t) and S 4 (t) from their initial value S 0 , also affect the results. Typically the linear model is a fair approximation when S 1 (t), S 2 (t), S 3 (t) and S 4 (t) decrease only by a few percentage points from their initial value S 0 .
We present our results for some typical values of ξ 1 , ξ 2 in tables 4, 5, 6 and 7 for the smallest value of β = 0.15 considered by us. It is in this case the fundamental mode is attained more slowly and is accompanied by a larger reduction of S i (t), i = 1, 2, 3, 4. Further we have chosen those cases when initial infections are all confined to region 1, while the fundamental mode predicts a much higher infections in region 2. Thus is a situation when non-linear effects are maximum during approach to equilibrium of linear models.
We begin by discussing table 4 where we tabulate S 1 (t), S 2 (t), S 3 (t) and S 4 (t) as well as M 1 (t), M 2 (t), M 3 (t) and M 4 (t) for ξ 1 = 0.3, ξ 2 = 0.1 and t ∈ (0, 150) days for the initial conditions M 1 (0) = M 3 (0) = 0.002 and M 2 (0) = M 4 (0) = 0.0. In this case the fundamental eigenvalue ν 1 = 0.0590833 of model (a) is somewhat larger than ν 3 = 0.05 but the second eigenvalue ν 2 = −0.0490833 is much lower than the eigenvalue ν 4 = −0.01. Thus it takes longer for model (b) to attain equillibrium than model (a). Numerical results show that it takes about 60 days for the fundamental mode to be established for model (a), i.e. infections are distributed in two regions in proportion to the fundamental eigenvector. In this time the susceptible population S 1 (t), S 2 (t), decreases by about 6%, 8% respectively. We see that at t = 60 the infected population M 1 (t) = 0.0197329 and M 2 (t) = 0.0272415. Their sum grows by a factor of 23.5 from its initial value of 0.002. However the zeroth order solution, Eq. (59), predicts M 1 (t) = 0.023181 and M 2 (t) = 0.032088. Clearly the zeroth order approximation is not a good measure of growth of the disease because of non-linear effects of reduction in S 1 (t), S 2 (t). In this interval about 4.2% of the population move to the "Removed" group by recovery or death. The ratio M 1 (t)/M 2 (t) = 1.3805 at t = 60 though is closer to 1.413 predicted by fundamental eigenvector. For model (b) the infected population M 3 (t) = 0.00967471, M 4 (t) = 0.0247318 at t = 60 does not compare well with their Eq. (60) estimate of 0.010866, 0.029305 respectively while the growth factor is 17.2. The susceptible populations S 3 (t), S 4 (t) decrease by 3.24%, 7.1% respectively. This is in conformity with the eigenvalue ν 3 being somewhat less than ν 1 . The percentage of population that moves to Removed group in this time interval is 3.5%. The effect of larger second eigenvalue ν 4 = −0.01 also shows up in the actual calculations. The ratio M 4 (t)/M 3 (t) = 2.556 at t = 60 as against the value 3.0 predicted by the eigenvector. This is partly due to slower decay of the second eigenvector and partly due to highly skewed fundamental eigenvector, three times in region 2 than in region 1. Reduction in S 4 (t) in this time interval also contributes. Overall we observe that the movements of susceptible increases the infections significantly, though not by an order of magnitude. It also leads to more even spatial distribution of fundamental mode. Further this relative proportion of two regions is maintained much longer, beyond the applicability of linearized model. This suggests that one need not solve for different initial distributions of infected people. It is sufficient to consider initial distribution as given by fundamental mode.
Above conclusions are corroborated by comparing them with results for the initial conditions M 1 (0) = 0.0002, M 2 (0) = 0.0 in table 5. Since the initial infections are an order of magnitude less, the non-linear effects are much less than in table 4. Thus at t = 60 the susceptible populations S 1 (t), S 2 (t), S 3 (t) and S 4 (t) all reduce by less than 1%. It is now seen that M 1 (t) = 0.00227942, M 2 (t) = 0.00320050 are much closer to their Eq. (59) estimates of 0.0023181, 0.0032088 and the ratio M 2 (t)/M 1 (t) = 1.404 is much closer to theoretical value 1.413. However, more glaring difference is in the results of model (b). The ratio M 4 (t)/M 3 (t) = 2.681 is higher than previos value of 2.556 at t = 60 and is still increasing with time. The value of this ratio 2.850 at t = 80, is still less than theoretical value 3.0. This is because of the slower decay of second eigenvector as ν 4 = −0.01. At t = 60, the values M 3 (t) = 0.00107362 and M 4 (t) = 0.0028791 are close to their zeroth order estimates. Comparing the results of tables 4 and 5 we conclude that if the initial infections are high then non-linear effects come into play earlier and interfere with the settling of fundamental mode. In that case the solution will depend on the initial locations of infections. However if the epidemic starts from small initial infections, then its future growth is largely determined by the fundamental eigenvector and the initial location of outbreak has less significance.
We next discuss the case ξ 1 = ξ 2 = 0.3 in table 6. In this case the fundamental eigenvalues ν 1 = ν 3 = 0.05 and hence the epidemic grows at the same rate for both models (a) and (b). Further the two eigenvectors are identical for both models, the fundamental mode is symmetric in regions 1 and (2)  in S 1 (t), S 2 (t), S 3 (t) and S 4 (t). If we reduce the initial infections to 0.0002, the reduction in S 1 (t), S 2 (t), S 3 (t) and S 4 (t) is less than 1% and the observed values of M 1 (t), M 2 (t), M 3 (t), M 4 (t) are much closer to their linearized estimates..
In previous cases the fundamental eigenvector for both models (a) and (b) are of similar nature. In table 4 it is greater in region 2 while in table 6 it is equal in both regions. Moreover the second eigenvector decayed faster in model (a) than in model (b). We now present a case where this eigenvector is higher in region 2 for model (a) while it is greater in region 1 for model (b). In table 7  There is an apparent mismatch between actual M 4 (t) value being higher than zeroth order estimate. This is observed for all t < 40 and disappears by t = 50 and is accompanied by significantly lower computed M 3 (t) compared to its zeroth estimate.
Lastly we observe that if the interaction between the regions is small, each region will evolve almost independently and thhe fundamental mode will hardly ever set in. In table 8 we present the results for the case β = 0.15, ξ 1 = ξ 2 = 0.1 with reduced initial conditions M 1 (0) = 0.0002, M 2 (0) = 0.0 so as to minimize the non-linear effects. Fundamental mode is symmetric in two regions for both the models (a) and (b) and so are the eigenvalues ν 1 = ν 3 = 0.05. The second eigenvector is antisymmetric with ν 2 = −0.004, small but negative, while ν 4 = 0.02 is positive. Table 8 shows that M 1 (t) ≈ M 2 (t) only by t = 100 when they differ by less than 1%. By this time S 1 (t), S 2 (t) reduce by nearly 4%. Model (b) results show that there are more than 1% difference between M 3 (t) and M 4 (t) even for t = 160, suggesting that fundamental mode does not get established.

CONCLUSIONS
In this paper we have extended the SAIR model of spread of an epidemic to a system comprising of many regions. Interaction between various regions is described by the fraction of population of any region that is present in other regions. We considered two models, one which allows for migration of Susceptible and Aymptomatic groups and the other which confines Susceptible population to their home region. We then developed ordinary differential equations that describe the time evolution of Susceptible, Asymptomatic carriers, Infected and Removed groups of people of all the regions. It is noted that only half the equations are non-linear that are also coupled for different regions. Further in the initial stages of the evolution of disease, the number of infected and asymptomatic carriers is much smaller than the susceptible population. This allows us to treat these equations by expanding in a perturbation series involving only linear equations in all orders of approximation. The zeroth order approximation yields a simple set of coupled, linear, first order ordinary differential equations involving a constant matrix of the interactions between various regions. We cast these equations in matrix form and note that same interaction matrix governs the evolution of all higher order terms, only the source terms are different. We then observed the crucial role played by the spectrum of this matrix and noted that a fundamental mode exists which grows much faster than other eigenvectors. This fundamental mode gets established if real part of other eigenvalues is sufficiently negative. Future spatial and temporal evolution of the disease then is governed by the shape of this mode only and the location of initial outbreak does not have much significance. If some other modes also grow, though slower than the fundamental one, then non-linear effects interfere with setting in of the dominant mode. It is noted that this is the case when the regions interact weakly with each other. We then applied above analysis to a simple two region problem for which lot of work can be carried out analytically. We verified our theoretical conclusions and confirmed that restricting the movement even of susceptible population only does lead to a slower growth of the epidemic. It was also found that once the fundamental mode sets in, subsequent time evolution roughly maintains the spatial shape much longer, well beyond the applicability of linearized model.
An apparent limitation of this analysis is that the spread of disease during commute from one region to another is neglected. It is well known that every major city has a well developed public transportation system where people come in close contact with each other. Likewise different cities and countries are linked through a travel network. This network makes a major contribution to spread of disease. It should be noted that this transport network can be regarded as another region with a population that equals the number of people normally using the network. The exchange of population from transport region to other regions and vice versa is high so that sum of people that go to (and from) other regions almost equals its entire population.