# The role of permanently resident populations in the two-patches SIR model with commuters

Alain Rapaport, Ismail Mimouni

► To cite this version:

Alain Rapaport, Ismail Mimouni. The role of permanently resident populations in the two-patches SIR model with commuters. Bulletin of Mathematical Biology, 2023, 85 (3), 10.1007/s11538-022-01111-6. hal-03872579

HAL Id: hal-03872579

https://hal.inrae.fr/hal-03872579

Submitted on 25 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The role of permanently resident populations in the two-patches SIR model with commuters

Alain Rapaport¹* and Ismail Mimouni¹

¹*MISTEA, Université Montpellier, INRAE, Institut Agro, place Viala, Montpellier, 34060, France.

*Corresponding author(s). E-mail(s): alain.rapaport@inrae.fr; Contributing authors: ismail.mimouni@yahoo.com;

## Abstract

We consider a two-patches SIR model where communication occurs thru commuters, distinguishing explicitly permanently resident populations from commuters populations. We give an explicit formula of the reproduction number, and show how the proportions of permanently resident populations impact it. We exhibit non-intuitive situations for which allowing commuting from a safe territory to another one where the transmission rate is higher can reduce the overall epidemic threshold and avoid an outbreak.

**Keywords:** SIR model, reproduction number, patches models, commuters

**MSC Classification:** 92D30, 34D20, 90C31

## 1 Introduction

Since the pioneer work of Kermack and McKendrick [26], the SIR model has been very popular in epidemiology, as the basic model for infectious diseases with direct transmission (e.g. [1]). It retakes great importance nowadays due to the recent coronavirus pandemic. While early models were not spatialized, the importance of accounting for spatial heterogeneity has been often reported in the literature (see, e.g. [2, 22, 24, 25, 28, 36]). However, different mechanisms come into play to explain the spatial spreading of a disease. Although diffusion appears to be a natural process to describe the local propagation of
2 *The role of permanently resident populations*

an infectious agent among a population, which leads to models with partial
differential equations [33], it appears to be not well suited for describing long
distance spreading. In particular, transportation between cities comes into the
picture as a major source of rapid spreading among non-homogeneous pop-
ulations [3, 5, 10, 29, 30, 32, 37, 38, 43]. Meta-populations or multi-patches
models are then more appropriate to describe the spatial characteristics of the
propagation [4, 6, 17, 41, 42], as already well considered in ecology [20, 31].
These models require a precise description of the movements between patches,
which are most of the time assumed to be linear and thus encoded into a con-
nection matrix [4, 6]. Typically one obtains a system of ordinary differential
equations on a graph, which couples the communication dynamics with the
epidemiological one.

For diseases spreading among human populations living in different cities,
*commuters* (individuals housing in a city, traveling regularly for short periods
in a neighboring city, and coming back to their home city) play a crucial role
in the disease propagation among territories [22–24, 32, 43]. Such coupling
between patches have been already considered in the literature, distinguishing
among populations $N_i$ attached to a city *i* the sub-population $N_{ii}$ present
in its permanent housing from other sub-populations $N_{ij}$ temporary present
in another city *j* ≠ *i* (it can be also seen as multi-groups models as in [11,
18, 21]). However, such models explicitly assume that the whole population
housing in a given city can potentially commute to another one. We believe
that this is not always fully realistic and that a sub-population that never
(or very rarely) moves to another city should be distinguished from the sub-
population that visits at a regular basis another city. Therefore, we consider
an extension of such models, which explicitly takes into consideration two
kinds of movement: an Eulerian one which describes the flow between patches
that mixes populations, and an Lagrangian one which assigns home locations
of individuals, as described in the more general framework [9]. The study of
this extension, which has not yet been analyzed analytically in the literature,
to our knowledge, and how it impacts the disease spreading, is the primary
objective of the present work. For this purpose, we establish an analytical
expression of the reproduction number (as the epidemic threshold formerly
introduced and analyzed in [14–16, 39]) for the two patches case (that is also
valid for the particular case when the whole populations travel, for which the
exact expression of the reproduction number has not been yet provided in the
literature).

We also had in mind to consider heterogeneity among territories when
disease transmission differs from one city to another one. Typically, non-
pharmaceutical interventions (such as reducing physical distance in the
population) could be applied with different strength in each city, providing dis-
tinct transmission rates. When one territory being isolated presents a higher
reproduction number than the other territory, it can be considered as a core
*group* in the epidemiological terminology [8, 19], and commuters contribute
then to spread the epidemics in both territories. We aim at analyzing more precisely how the proportions of commuters in each city can increase or decrease the overall reproduction number. Intuitively, one may believe that the best way to reduce the spreading is to encourage commuters from the city with the lowest transmission rate not to travel to the other city, and on the opposite to encourage as much as possible commuters from the other city to spend time in the safer city. Indeed, we shall see that this is not always true... The second objective of the present work is thus to study the minimization of the epidemic threshold of the two-patches model with respect to these proportions, depending on the commuting rates. This analysis can potentially serve for decisions making to prevent epidemic outbreak (as in [27] for instance).

The paper is organized as follows. In the next section, we present the complete model in dimension 18 and give some preliminaries. Section 3 is devoted to the analysis of the asymptotic behavior of the solutions of the model. We give and demonstrate an explicit expression of the reproduction number, introducing four relevant quantities $q_{ij}$ ($i, j = 1, 2$). In a corollary, we also give an alternative way of computation, which is useful in the following. In Section 4, we study the minimization of the reproduction number with respect to the proportions of commuters in each patch. Finally, Section 5 gives a numerical illustration of the results, considering two territories with intrinsic basic reproduction numbers lower and higher than one. We depict the relative sizes of the permanently resident populations that can avoid the outbreak of the epidemic depending on the commuting rates, and discuss the various cases. We end with a conclusion.

## 2 The model

We follow the modeling of commuters proposed in [23] between two patches (such as cities or territories), but here we consider in addition that a part of the population in each patch do not commute (the *permanently resident* sub-population). We consider populations of size $N_i$ whose home belongs to a patch $i \in \{1, 2\}$, structured in three groups:

1. permanently resident, being all the time in patch $i$, whose population size is denoted $N_{ir}$,

2. commuters to patch $j$, but located in patch $i$ at time $t$, of population size denoted $N_{ii}$,

3. commuters to patch $j$ and located in patch $j$ at time $t$, of population size denoted $N_{ij}$.

We shall denote $N_{ic} = N_{ii} + N_{ij}$ the size of the total population of commuters with their home in patch $i$. The individuals commutes to patch $j$ at a rate $\lambda_i$ with a return rate $\mu_i$. For each group $g \in \{ir, ii, ij\}$ we denote by $S_g, I_g, R_g$ the sizes of susceptible, infected and recovered sub-populations. This modeling implicitly assumes that at any time there is no individual out the territories, that is traveling time is negligible. This assumption is therefore only valid
4 *The role of permanently resident populations*

for adjoining territories with short transportation times (by train, road...). It
would not be valid between distant territories connected for example by boat
or plane with non-negligible crossing times. In this case, it would be necessary
to consider additional nodes of in-transit populations, as it has been considered
for example in [12, 34] or in [35] where distance between nodes are explicitly
taken into consideration. This would of course complicates the model and its
study.

We consider the SIR model assuming that the recovery parameter $\gamma$ is
identical everywhere while the transmission rate $\beta_i$ depends on the patch $i$ but
is identical among each group. Typically lifestyle and hygienic measures may
differ between two cities, implying different values of $\beta$. Moreover, if two cities
are on both sides of the border between two countries, the strength of non-
pharmaceutical interventions are likely to be different, as is was for instance
the case between European countries during the SARS-2 outbreak. The model
is written as follows (with $i \neq j$ in {1,2}).

$$
\dot{S}_{ir} = -\beta_i S_{ir} \frac{I_{ir} + I_{ii} + I_{ji}}{N_{ir} + N_{ii} + N_{ji}},
$$

$$
\dot{I}_{ir} = \beta_i S_{ir} \frac{I_{ir} + I_{ii} + I_{ji}}{N_{ir} + N_{ii} + N_{ji}} - \gamma I_{ir},
$$

$$
\dot{R}_{ir} = \gamma I_{ir},
$$

$$
\dot{S}_{ii} = -\beta_i S_{ii} \frac{I_{ir} + I_{ii} + I_{ji}}{N_{ir} + N_{ii} + N_{ji}} - \lambda_i S_{ii} + \mu_i S_{ij},
$$

$$
\dot{I}_{ii} = \beta_i S_{ii} \frac{I_{ir} + I_{ii} + I_{ji}}{N_{ir} + N_{ii} + N_{ji}} - \gamma I_{ii} - \lambda_i I_{ii} + \mu_i I_{ij},
$$

$$
\dot{R}_{ii} = \gamma I_{ii} - \lambda_i R_{ii} + \mu_i R_{ij},
$$

$$
\dot{S}_{ij} = -\beta_j S_{ij} \frac{I_{jr} + I_{jj} + I_{ij}}{N_{jr} + N_{jj} + N_{ij}} + \lambda_i S_{ii} - \mu_i S_{ij},
$$

$$
\dot{I}_{ij} = \beta_j S_{ij} \frac{I_{jr} + I_{jj} + I_{ij}}{N_{jr} + N_{jj} + N_{ij}} - \gamma I_{ij} + \lambda_i I_{ii} - \mu_i I_{ij},
$$

$$
\dot{R}_{ij} = \gamma I_{ij} + \lambda_i R_{ii} - \mu_i R_{ij}
$$

Parameters $\lambda_i$, $\mu_i$ represent switching rates of populations $i$, leaving home and returning. This modeling implicitly assumes that movements between territories are not synchronized, as often considered in multi-city models (see e.g. [5, 6, 10, 22–24, 30, 36, 37, 41, 42]). Note that we also consider, in all generality, that commuting is asymmetrical (i.e. $\lambda_1$ and $\lambda_2$ may be different, as well as $\mu_1, \mu_2$). Typically, each territory may offer different activities that attract commuters from the other territory, and thus different mean sojourn times. One can check that the population sizes $N_{ir}$, and $N_{ic}$ are constant. Moreover $N_{ii}$, $N_{ij}$ fulfill the system of equations
$$ \begin{cases} \dot{N}_{ii} = -\lambda_i N_{ii} + \mu_i N_{ij}, \\ \dot{N}_{ij} = \lambda_i N_{ii} - \mu_i N_{ij} \end{cases} $$

whose solutions verify

$$ \lim_{t \to +\infty} N_{ii}(t) = \bar{N}_{ii} := \frac{\mu_i}{\lambda_i + \mu_i} N_{ic}, \quad \lim_{t \to +\infty} N_{ij}(t) = \bar{N}_{ij} := \frac{\lambda_i}{\lambda_i + \mu_i} N_{ic} \quad (1) $$

We shall assume that populations are already balanced at initial time i.e. that one has $N_{ii} = \bar{N}_{ii}$, $N_{ij} = \bar{N}_{ij}$ (constant). For simplicity, we shall drop the notation $^-$ in the following, and denote

$$ N_{ip} := N_{ir} + N_{ii} + N_{ji} $$

which represents the (constant) size of the total population present in patch i.

# 3 The epidemic threshold

We denote the vectors

$$ I = (I_{1r}, I_{11}, I_{12}, I_{2r}, I_{22}, I_{21})^{\top}, \quad S = (S_{1r}, S_{11}, S_{12}, S_{2r}, S_{22}, S_{21})^{\top} $$

and consider the state vector

$$ X = \begin{bmatrix} I \\ S \end{bmatrix} $$

which belongs to the invariant domain

$$ \mathcal{D} := \{ X \in \mathbb{R}_+^{12} ; \mathbb{M} X \leq \mathbb{N} \} $$

where $\mathbb{N}$ is the vector

$$ \mathbb{N} = (N_{1r}, N_{11}, N_{12}, N_{2r}, N_{22}, N_{21})^{\top} $$

and $\mathbb{M}$ the 6 × 12 matrix which consists in the concatenation of the identity matrix $\mathbb{I}_6$ of dimension 6 × 6

$$ \mathbb{M} = [\mathbb{I}_6, \mathbb{I}_6] $$

The disease free equilibrium is defined as

$$ X^* = \begin{bmatrix} 0 \\ \mathbb{N} \end{bmatrix} $$

Let $\mathcal{R}_i$ be the intrinsic reproduction number in the patch i (i.e. when there is no connection between patches), that is
6 *The role of permanently resident populations*

$$
\mathcal{R}_i := \frac{\beta_i}{\gamma}.
$$

We give now an explicit expression of the epidemic threshold when the two
patches communicates via commuters.

**Proposition 1** Let

$$
\mathcal{R}_{1,2} := \frac{q_{11} + q_{22} + \sqrt{(q_{22} - q_{11})^2 + 4q_{12}q_{21}}}{2} \quad (2)
$$

where

$$
\left\{
\begin{aligned}
q_{11} &= \mathcal{R}_1 \left( \frac{N_{1r}}{N_{1p}} + \frac{N_{11}}{N_{1p}} \frac{\gamma+\mu_1}{\gamma+\lambda_1+\mu_1} + \frac{N_{21}}{N_{1p}} \frac{\gamma+\lambda_2}{\gamma+\lambda_2+\mu_2} \right) \\
q_{22} &= \mathcal{R}_2 \left( \frac{N_{2r}}{N_{2p}} + \frac{N_{22}}{N_{2p}} \frac{\gamma+\mu_2}{\gamma+\lambda_2+\mu_2} + \frac{N_{12}}{N_{2p}} \frac{\gamma+\lambda_1}{\gamma+\lambda_1+\mu_1} \right) \\
q_{21} &= \mathcal{R}_1 \left( \frac{N_{11}}{N_{1p}} \frac{\lambda_1}{\gamma+\lambda_1+\mu_1} + \frac{N_{21}}{N_{1p}} \frac{\mu_2}{\gamma+\lambda_2+\mu_2} \right) \\
q_{12} &= \mathcal{R}_2 \left( \frac{N_{12}}{N_{2p}} \frac{\mu_1}{\gamma+\lambda_1+\mu_1} + \frac{N_{22}}{N_{2p}} \frac{\lambda_2}{\gamma+\lambda_2+\mu_2} \right)
\end{aligned}
\right.
\quad (3)
$$

Then, one has the following properties.

1. If $\mathcal{R}_{1,2} > 1$, then $X^*$ is unstable.

2. If $\mathcal{R}_{1,2} < 1$, then $X^*$ is exponentially stable with respect to the variable<fnref n="1" /> I.

3. If $\mathcal{R}_1 = \mathcal{R}_2 := \mathcal{R}$, then $\mathcal{R}_{1,2} = \mathcal{R}$.

*Proof* Write the dynamics of *X* as *Ẋ* = *f*(*X*). The Jacobian matrix *J* of *f* at *Ẋ* is
of the form

$$
J = \begin{bmatrix} A & 0 \\ \star & B \end{bmatrix} \text{ with } A = F - V
$$

where

$$
F = \begin{bmatrix}
\beta_1 \frac{N_{1r}}{N_{1p}} & \beta_1 \frac{N_{1r}}{N_{1p}} & 0 & 0 & 0 & \beta_1 \frac{N_{1r}}{N_{1p}} \\
\beta_1 \frac{N_{11}}{N_{1p}} & \beta_1 \frac{N_{11}}{N_{1p}} & 0 & 0 & 0 & \beta_1 \frac{N_{11}}{N_{1p}} \\
0 & 0 & \beta_2 \frac{N_{12}}{N_{2p}} & \beta_2 \frac{N_{12}}{N_{2p}} & \beta_2 \frac{N_{12}}{N_{2p}} & 0 \\
0 & 0 & \beta_2 \frac{N_{2r}}{N_{2p}} & \beta_2 \frac{N_{2r}}{N_{2p}} & \beta_2 \frac{N_{2r}}{N_{2p}} & 0 \\
0 & 0 & \beta_2 \frac{N_{22}}{N_{2p}} & \beta_2 \frac{N_{22}}{N_{2p}} & \beta_2 \frac{N_{22}}{N_{2p}} & 0 \\
\beta_1 \frac{N_{21}}{N_{1p}} & \beta_1 \frac{N_{21}}{N_{1p}} & 0 & 0 & 0 & \beta_1 \frac{N_{21}}{N_{1p}}
\end{bmatrix},
$$

<footnote n="1">¹We refer to [40] for the definition of partial stability.</footnote>
$$V = \begin{bmatrix} \gamma & 0 & 0 & 0 & 0 & 0 \\ 0 & \gamma + \lambda_1 & -\mu_1 & 0 & 0 & 0 \\ 0 & -\lambda_1 & \gamma + \mu_1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \gamma & 0 & 0 \\ 0 & 0 & 0 & 0 & \gamma + \lambda_2 & -\mu_2 \\ 0 & 0 & 0 & 0 & -\lambda_2 & \gamma + \mu_2 \end{bmatrix}$$

and

$$B = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\lambda_1 & \mu_1 & 0 & 0 & 0 \\ 0 & \lambda_1 & -\mu_1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -\lambda_2 & \mu_2 \\ 0 & 0 & 0 & 0 & \lambda_2 & -\mu_2 \end{bmatrix}$$

Note that $F$ is a non-negative matrix and $V$ is a non-singular M-matrix. We recall (see for instance from [39]) that one has the property

$$\max \mathrm{Re}(\mathrm{Spec}(A)) \le 0 \iff \rho(FV^{-1}) \le 1$$

The computation of the matrix $M := FV^{-1}$ gives the following expression

$$M = \begin{bmatrix} \mathcal{R}_1 \frac{N_{1r}}{N_{1p}} & \mathcal{R}_1 \frac{N_{1r}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_1 \frac{N_{1r}\mu_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & 0 & \mathcal{R}_1 \frac{N_{1r}\lambda_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{1r}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\ \mathcal{R}_1 \frac{N_{11}}{N_{1p}} & \mathcal{R}_1 \frac{N_{11}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_1 \frac{N_{11}\mu_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & 0 & \mathcal{R}_1 \frac{N_{11}\lambda_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{11}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\ 0 & \mathcal{R}_2 \frac{N_{12}\lambda_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}}{N_{2p}} & \mathcal{R}_2 \frac{N_{12}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_2 \frac{N_{12}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\ 0 & \mathcal{R}_2 \frac{N_{2r}\lambda_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{2r}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{2r}}{N_{2p}} & \mathcal{R}_2 \frac{N_{2r}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_2 \frac{N_{2r}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\ 0 & \mathcal{R}_2 \frac{N_{22}\lambda_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{22}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{22}}{N_{2p}} & \mathcal{R}_2 \frac{N_{22}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_2 \frac{N_{22}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\ \mathcal{R}_1 \frac{N_{21}}{N_{1p}} & \mathcal{R}_1 \frac{N_{21}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_1 \frac{N_{21}\mu_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & 0 & \mathcal{R}_1 \frac{N_{21}\lambda_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{21}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} \end{bmatrix}$$

Let us consider the diagonal matrix

$$D := \begin{bmatrix} \mathcal{R}_1 \frac{N_{1r}}{N_{1p}} & & & & & \\ & \mathcal{R}_1 \frac{N_{11}}{N_{1p}} & & & & \\ & & \mathcal{R}_2 \frac{N_{12}}{N_{2p}} & & & \\ & & & \mathcal{R}_2 \frac{N_{2r}}{N_{2p}} & & \\ & & & & \mathcal{R}_2 \frac{N_{22}}{N_{2p}} & \\ & & & & & \mathcal{R}_1 \frac{N_{21}}{N_{1p}} \end{bmatrix}$$
and the matrix $Q = D^{-1}MD$, whose computation gives the expression

$$Q = \begin{bmatrix}
\mathcal{R}_1 \frac{N_{1r}}{N_{1p}} & \mathcal{R}_1 \frac{N_{11}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}\mu_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & 0 & \mathcal{R}_2 \frac{N_{22}\lambda_2}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{21}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\
\mathcal{R}_1 \frac{N_{1r}}{N_{1p}} & \mathcal{R}_1 \frac{N_{11}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}\mu_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & 0 & \mathcal{R}_2 \frac{N_{22}\lambda_2}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{21}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\
0 & \mathcal{R}_1 \frac{N_{11}\lambda_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{2r}}{N_{2p}} & \mathcal{R}_2 \frac{N_{22}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{21}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\
0 & \mathcal{R}_1 \frac{N_{11}\lambda_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{2r}}{N_{2p}} & \mathcal{R}_2 \frac{N_{22}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{21}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\
0 & \mathcal{R}_1 \frac{N_{11}\lambda_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{2r}}{N_{2p}} & \mathcal{R}_2 \frac{N_{22}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{21}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\
\mathcal{R}_1 \frac{N_{1r}}{N_{1p}} & \mathcal{R}_1 \frac{N_{11}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \mathcal{R}_2 \frac{N_{12}\mu_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & 0 & \mathcal{R}_2 \frac{N_{22}\lambda_2}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \mathcal{R}_1 \frac{N_{21}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)}
\end{bmatrix}$$

The matrix $Q$ is non-negative and irreducible. By Perron-Frobenius Theorem (see for instance [7]), this matrix admits a unique positive eigenvector (up to a scalar multiplication) that corresponds to the simple (positive) eigenvalue $\ell = \rho(Q) = \rho(M)$.

Note that the rank of $Q$ is two. We posit

$$Y = (1, 1, 0, 0, 0, 1)^{\top}, \quad Z = (0, 0, 1, 1, 1, 0)^{\top}$$

and define $Q_Y$, $Q_Z$ the first and third lines, respectively, of the matrix $Q$. Then, for any vector $X \in \mathbb{R}^6$, on has one has $QX = (Q_Y X)Y + (Q_Z X)Z$. We look for an positive eigenvector $X$ of the form $X = \alpha Y + (1-\alpha)Z$ with $\alpha \in (0, 1)$. One has then

$$\begin{aligned}
QX &= \alpha QY + (1 - \alpha)QZ \\
&= \alpha((Q_Y Y)Y + (Q_Z Y)Z) + (1 - \alpha)((Q_Y Z)Y + (Q_Z Z)Z) \\
&= (\alpha(Q_Y Y) + (1 - \alpha)(Q_Y Z))Y + (\alpha(Q_Z Y) + (1 - \alpha)(Q_Z Z))Z \quad (4)
\end{aligned}$$

On the other hand, as $X$ is an eigenvector, one has

$$QX = \ell X = \alpha \ell Y + (1 - \alpha) \ell Z \quad (5)$$

The vectors $Y$ and $Z$ being orthogonal, one obtains from (4)-(5) the conditions

$$\begin{cases}
\alpha Q_Y Y + (1 - \alpha) Q_Y Z = \alpha \ell \\
\alpha Q_Z Y + (1 - \alpha) Q_Z Z = (1 - \alpha) \ell
\end{cases}
\quad (6)$$

Let $r = \frac{1-\alpha}{\alpha}$. Eliminating $\ell$ in the two previous equations, $r$ is the positive solution of the polynomial

$$r^2 Q_Y Z + r(Q_Y Y - Q_Z Z) - Q_Z Y = 0$$

and $\ell = Q_Y Y + rQ_Y Z$. One obtains the expression of the eigenvalue

$$\ell = \frac{Q_Y Y + Q_Z Z + \sqrt{(Q_Y Y - Q_Z Z)^2 + 4(Q_Y Z)(Q_Z Y)}}{2}$$

Finally, from the expression of $Q$, one gets

$$q_{11} = Q_Y Y = \mathcal{R}_1 \left( \frac{N_{1r}}{N_{1p}} + \frac{N_{11}}{N_{1p}} \frac{\gamma+\mu_1}{\gamma+\lambda_1+\mu_1} + \frac{N_{21}}{N_{2p}} \frac{\gamma+\lambda_2}{\gamma+\lambda_2+\mu_2} \right)$$

$$q_{22} = Q_Z Z = \mathcal{R}_2 \left( \frac{N_{2r}}{N_{2p}} + \frac{N_{22}}{N_{2p}} \frac{\gamma+\mu_2}{\gamma+\lambda_2+\mu_2} + \frac{N_{12}}{N_{1p}} \frac{\gamma+\lambda_1}{\gamma+\lambda_1+\mu_1} \right)$$

$$q_{21} = Q_Z Y = \mathcal{R}_1 \left( \frac{N_{11}}{N_{1p}} \frac{\lambda_1}{\gamma+\lambda_1+\mu_1} + \frac{N_{21}}{N_{2p}} \frac{\mu_2}{\gamma+\lambda_2+\mu_2} \right)$$

$$q_{12} = Q_Y Z = \mathcal{R}_2 \left( \frac{N_{12}}{N_{2p}} \frac{\mu_1}{\gamma+\lambda_1+\mu_1} + \frac{N_{22}}{N_{2p}} \frac{\lambda_2}{\gamma+\lambda_2+\mu_2} \right)$$

and thus $\ell = \mathcal{R}_{1,2}$, which is exactly $\rho(M)$.
i. When $\mathcal{R}_{1,2} > 1$, the matrix $A$ has at least one eigenvalue with positive real part and the matrix $J$ as well. The equilibrium $X^*$ is thus unstable on $D$.

ii. When $\mathcal{R}_{1,2} < 1$, the matrix $A$ is Hurwitz, but $X^*$ is not an hyperbolic equilibrium. However, on can write the dynamics of the vector $I$ as an non-autonomous system

$$ \dot{I} = g(t, I) := \begin{pmatrix} \beta_1 S_{1r}(t) \frac{I_{1r} + I_{11} + I_{21}}{N_{1p}} - \gamma I_{1r} \\ \beta_1 S_{11}(t) \frac{I_{1r} + I_{11} + I_{21}}{N_{1p}} - (\gamma + \lambda_1) I_{11} + \mu_1 I_{12} \\ \beta_2 S_{12}(t) \frac{I_{2r} + I_{22} + I_{12}}{N_{2p}} + \lambda_1 I_{11} - (\gamma + \mu_1) I_{12} \\ \beta_2 S_{2r}(t) \frac{I_{2r} + I_{22} + I_{12}}{N_{2p}} - \gamma I_{2r} \\ \beta_2 S_{22}(t) \frac{I_{2r} + I_{22} + I_{12}}{N_{2p}} - (\gamma + \lambda_2) I_{22} + \mu_2 I_{21} \\ \beta_1 S_{21}(t) \frac{I_{1r} + I_{11} + I_{21}}{N_{1p}} + \lambda_2 I_{22} - (\gamma + \mu_2) I_{21} \end{pmatrix} $$

Note that this dynamics is cooperative and as for any $t \ge 0$ one has $S_{ij}(t) \le N_{ij}$ for $ij \in \{1r, 11, 12, 2r, 22, 21\}$, one gets

$$ g(t, I) \le \bar{g}(I) := AI, \quad I \ge 0 $$

Therefore, any solution $I(\cdot)$ of $\dot{I} = g(t, I)$ with $I(0) = I_0 \ge 0$ verifies $0 \le I(t) \le \bar{I}(t)$ for any $t \ge 0$, where $\bar{I}(\cdot)$ is solution of the linear dynamics $\dot{I} = \bar{g}(\bar{I})$ with $\bar{I}(0) = I_0$. As $A$ is Hurwitz, we conclude that $X^*$ is exponentially stable with respect to $I$, which proves point ii.

iii. For the particular case $\mathcal{R}_1 = \mathcal{R}_2 := \mathcal{R}$, the transpose of the matrix $M$ writes

$$ M^\top = \mathcal{R} \begin{bmatrix} \frac{N_{1r}}{N_{1p}} & \frac{N_{11}}{N_{1p}} & 0 & 0 & 0 & \frac{N_{21}}{N_{1p}} \\ \frac{N_{1r}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{11}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{12}\lambda_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{2r}\lambda_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{22}\lambda_1}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{21}(\gamma+\mu_1)}{N_{1p}(\gamma+\lambda_1+\mu_1)} \\ \frac{N_{1r}\mu_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{11}\mu_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{12}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{2r}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{22}(\gamma+\lambda_1)}{N_{2p}(\gamma+\lambda_1+\mu_1)} & \frac{N_{21}\mu_1}{N_{1p}(\gamma+\lambda_1+\mu_1)} \\ 0 & 0 & \frac{N_{12}}{N_{2p}} & \frac{N_{2r}}{N_{2p}} & \frac{N_{22}}{N_{2p}} & 0 \\ \frac{N_{1r}\lambda_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{11}\lambda_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{12}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{2r}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{22}(\gamma+\mu_2)}{N_{2p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{21}\lambda_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} \\ \frac{N_{1r}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{11}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{12}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{2r}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{22}\mu_2}{N_{1p}(\gamma+\lambda_2+\mu_2)} & \frac{N_{21}(\gamma+\lambda_2)}{N_{1p}(\gamma+\lambda_2+\mu_2)} \end{bmatrix} $$

One can check that one has $M^\top U = \mathcal{R}U$ where $U = (1, 1, 1, 1, 1)^\top$. As $U$ is a positive vector, we deduce from the Perron-Frobenius Theorem that one has $\rho(M) = \rho(M^\top) = \mathcal{R}$, which ends the proof. $\square$

**Remark 1** More generally, the next-generation matrix $M = FV^{-1}$ can be shown to have a rank equal to the number $n$ of patches, and that its Perron vector can be expressed as a convex combination of a family of orthogonal vectors in the image of $M$. This implies that the positive eigenvalue of $M$ (i.e. the reproduction number) is also the positive eigenvalue of the $n$-dimensional positive matrix given by the decomposition of the image of this vectors by the matrix $M$.

Alternatively, one may consider the epidemic spread in a virgin population as a Markov process, to determine the expected numbers of secondary cases in each patch, and obtain this $n \times n$ matrix, as described in [13]. This method consists in a first-step
10 The role of permanently resident populations

analysis by determining the mean residence times of an infected individual of each group in each of the patches. Then, for a given patch the expected numbers of new infected present in each path are given by the products of the mean residence times by the transmission rate, averaged by the constant distribution given in (1).

This explains why the formula (2) takes the expression of a root of the characteristic polynomial of a 2 by 2 matrix.

Remark 2 The explicit expression (2) of the epidemic threshold given in Proposition 1 is also relevant in absence of permanently resident populations, which has not been yet provided explicitly in the literature (up to our knowledge).

**Corollary 2** One has

$$ \min(\mathcal{R}_1, \mathcal{R}_2) \leq \mathcal{R}_{1,2} \leq \max(\mathcal{R}_1, \mathcal{R}_2). $$

*Proof* Denote by $M(\mathcal{R}_1, \mathcal{R}_2)$ the matrix $FV^{-1}$ for the parameters $\mathcal{R}_1, \mathcal{R}_2$, and let $\mathcal{R}_- := \min(\mathcal{R}_1, \mathcal{R}_2)$, $\mathcal{R}_+ := \max(\mathcal{R}_1, \mathcal{R}_2)$. From the expression of the non-negative matrices $M$, one gets

$$ M(\mathcal{R}_-, \mathcal{R}_-) \leq M(\mathcal{R}_1, \mathcal{R}_2) \leq M(\mathcal{R}_+, \mathcal{R}_+) $$

which implies (see for instance [7]) the inequalities

$$ \rho(M(\mathcal{R}_-, \mathcal{R}_-)) \leq \rho(M(\mathcal{R}_1, \mathcal{R}_2)) \leq \rho(M(\mathcal{R}_+, \mathcal{R}_+)) $$

and thus

$$ \mathcal{R}_- \leq \mathcal{R}_{1,2} \leq \mathcal{R}_+. $$

□

Alternatively, the number $\mathcal{R}_{1,2}$ can be determined as follows.

**Corollary 3** Assume $\mathcal{R}_2 > \mathcal{R}_1$. Then, one has

$$ \mathcal{R}_{1,2} = \alpha \mathcal{R}_1 + (1 - \alpha) \mathcal{R}_2 \qquad (7) $$

where $\alpha \in [0, 1)$ is the smallest root of the polynomial

$$ P(\alpha) = \alpha^2(\mathcal{R}_2 - \mathcal{R}_1) - \alpha(\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21}) + q_{12} $$

*Proof* One can check, from expressions (3), that one has $q_{11} + q_{21} = \mathcal{R}_1$ and $q_{22} + q_{12} = \mathcal{R}_2$. Then, from (6), one get

$$ \mathcal{R}_{1,2} = l = \alpha \mathcal{R}_1 + (1 - \alpha) \mathcal{R}_2 \qquad (8) $$

where $\alpha$ is a root of the polynomial $P$ obtained from (6) by eliminating $l$, that is

$$ P(\alpha) = \alpha^2(\mathcal{R}_2 - \mathcal{R}_1) - \alpha(\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21}) + q_{12} $$

From Corollary 2, we know that $\alpha$ belongs to $[0, 1]$. Note that one has $P(0) = q_{12} \geq 0$ and $P(1) = -q_{21} \leq 0$. Therefore, when $\mathcal{R}_2 - \mathcal{R}_1 > 0$, $P$ admits exactly one root
in [0, 1) and another one in [1, →). However, if α = 1 one should have q$_{21}$ = 0 and thus λ$_{1}$ = 0, µ$_{2}$ = 0, which implies N$_{11}$ = N$_{1c}$, N$_{12}$ = 0, N$_{22}$ = 0, N$_{21}$ = N$_{2c}$ . Then, one obtains q$_{11}$ = R $_{1}$, q$_{22}$ = R 2 and from the expression (2) on gets R 1 , 2 = max( R $_{1}$, R $_{2}$) = R 2 which contradicts α = 1. We conclude that α belongs to [0 , 1) and is thus the smallest root of P . □

Remark 3 When there is no communication between patches (that is N$_{1}$$_{r}$ = N$_{1}$$_{p}$ = N$_{1}$ , N$_{2}$$_{r}$ = N$_{2}$$_{p}$ = N$_{2}$ ), one has q$_{21}$ = 0 and q$_{12}$ = 0. If R 2 > R $_{1}$, resp. R 1 > R $_{2}$, one has α = 0, resp. α = 1, which gives

$$
\mathcal{R}_{1,2} = \max(\mathcal{R}_1, \mathcal{R}_2).
$$

We look now for a characterization of the minimum value of the threshold R 1 , $_{2}$.

4 Minimization of the epidemic threshold

In this section, we assume that the mixing is fast compared to the recovery rate (as its is often considered in the literature), which amounts to have numbers λ$_{i}$, µ$_{i}$ large compared to γ . Our objective is to study how the proportions of commuters in the populations impact the value of R 1 , $_{2}$.

Given R $_{1}$, R $_{2}$, we consider the approximation ˜ R 1 , 2 of the threshold R 1 , 2 which consists in keeping γ = 0 in the expressions (3). For convenience, we posit the numbers

$$
\eta_i := \frac{\lambda_i}{\lambda_i + \mu_i} \in (0, 1) \quad (i = 1, 2)
$$

One has a first result about the variations of ˜ R 1 , 2 with respect to N$_{1}$$_{c}$, N$_{2}$$_{c}$.

**Proposition 4** Fix parameters N$_{i}$, β$_{i}$, γ, λ$_{i}$, µ$_{i}$ ( i = 1 , 2) such that R 2 > R $_{1}$.

1. For any N$_{1}$$_{c}$ ∈ (0 , N$_{1}$ ), the map N$_{2}$$_{c}$ ↦→ ˜ R 1 , $_{2}$( N$_{1}$$_{c}$, N$_{2}$$_{c}$ ) is decreasing.

2. The map N$_{1}$$_{c}$ ↦→ ˜ R 1 , $_{2}$( N$_{1}$$_{c}$, N$_{2}$$_{c}$ ) is increasing at ( N$_{1}$$_{c}$, N$_{2}$$_{c}$ ) when

$$
\eta_2(1 - \eta_2)N_{2c} > (1 - \eta_1)(N_2 - \eta_2 N_{2c}) \quad (9)
$$

3. The map N$_{1}$$_{c}$ ↦→ ˜ R 1 , $_{2}$( N$_{1}$$_{c}$, N$_{2}$$_{c}$ ) is increasing, resp. decreasing, at ( N$_{1}$$_{c}$, N$_{2}$$_{c}$ ) if the numbers A and B are negative, resp. positive, where

$$
\begin{align*}
A := \mathcal{R}_2 \frac{\frac{N_2}{2} - \eta_1(\frac{1}{2} - \eta_1)N_{1c} - (\frac{3}{2} - \eta_2)\eta_2 N_{2c}}{N_2 - \eta_2 N_{2c} + \eta_1 N_{1c}} \\
\qquad - \mathcal{R}_1 \frac{\frac{N_1}{2} - (\frac{3}{2} - \eta_1)\eta_1 N_{1c} - \eta_2(\frac{1}{2} - \eta_2)N_{2c}}{N_1 - \eta_1 N_{1c} + \eta_2 N_{2c}},
\end{align*}
$$
12 The role of permanently resident populations

$$
\begin{align*}
B := \mathcal{R}_2 \frac{(1 - \eta_1)(N_2 - \eta_2 N_{2c}) - \eta_2(1 - \eta_2)N_{2c}}{(N_2 - \eta_2 N_{2c} + \eta_1 N_{1c})^2} & \\
& - \mathcal{R}_1 \frac{(1 - \eta_1)(N_1 + \eta_2 N_{2c}) + \eta_2(1 - \eta_2)N_{2c}}{(N_1 - \eta_1 N_{1c} + \eta_2 N_{2c})^2}
\end{align*}
$$

Proof Following Corollary 3, one has

$$
\tilde{\mathcal{R}}_{1,2} = \tilde{\alpha}\mathcal{R}_1 + (1 - \tilde{\alpha})\mathcal{R}_2 \quad (10)
$$

where $\tilde{\alpha}$ is the smallest root of the polynomial

$$
\tilde{P}(\alpha) = \alpha^2(\mathcal{R}_2 - \mathcal{R}_1) - \alpha(\mathcal{R}_2 - \mathcal{R}_1 + \tilde{q}_{12} + \tilde{q}_{21}) + \tilde{q}_{12}
$$

where $\tilde{q}_{12}, \tilde{q}_{21}$ are the approximations of $q_{12}, q_{21}$ defined in (3). Let us note that one
can write $N_{ii} = (1 - \eta_i)N_{ic}$, $N_{ij} = \eta_i N_{ic}$ (for $j \neq i$) and also $N_{ip} = N_i - \eta_i N_{ic} +
\eta_j N_{jc}$, which leads to the following expressions of $\tilde{q}_{12}, \tilde{q}_{21}$

$$
\tilde{q}_{21} = \mathcal{R}_1 \frac{(1 - \eta_1)\eta_1 N_{1c} + \eta_2(1 - \eta_2)N_{2c}}{N_1 - \eta_1 N_{1c} + \eta_2 N_{2c}}, \quad \tilde{q}_{12} = \mathcal{R}_2 \frac{\eta_1(1 - \eta_1)N_{1c} + (1 - \eta_2)\eta_2 N_{2c}}{N_2 - \eta_2 N_{2c} + \eta_1 N_{1c}} \quad (11)
$$

For simplicity, we shall drop the notation $^{−}$ in the rest of the proof. Note that $\alpha$
being the smallest root of P, it verifies

$$
\alpha < \frac{\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21}}{2(\mathcal{R}_2 - \mathcal{R}_1)} \quad (12)
$$

Let us differentiate the equality $P(\alpha) = 0$ with respect to $q_{12}$ and $q_{21}$:

$$
2\alpha \frac{\partial \alpha}{\partial q_{12}} (\mathcal{R}_2 - \mathcal{R}_1) - \frac{\partial \alpha}{\partial q_{12}} (\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21}) - \alpha + 1 = 0
$$

$$
2\alpha \frac{\partial \alpha}{\partial q_{21}} (\mathcal{R}_2 - \mathcal{R}_1) - \frac{\partial \alpha}{\partial q_{21}} (\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21}) - \alpha = 0
$$

which gives

$$
\frac{\partial \alpha}{\partial q_{12}} = \frac{1 - \alpha}{\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21} - 2\alpha(\mathcal{R}_2 - \mathcal{R}_1)}
$$

$$
\frac{\partial \alpha}{\partial q_{21}} = \frac{-\alpha}{\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21} - 2\alpha(\mathcal{R}_2 - \mathcal{R}_1)}
$$

Then, one can write

$$
\frac{\partial \alpha}{\partial N_{ic}} = \frac{\partial \alpha}{\partial q_{12}} \frac{\partial q_{12}}{\partial N_{ic}} + \frac{\partial \alpha}{\partial q_{21}} \frac{\partial q_{21}}{\partial N_{ic}} = \frac{(1-\alpha) \frac{\partial q_{12}}{\partial N_{ic}} - \alpha \frac{\partial q_{21}}{\partial N_{ic}}}{\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21} - 2\alpha(\mathcal{R}_2 - \mathcal{R}_1)} \quad (i=1,2)
$$

and from inequality (12), we obtain that the signs of the derivatives $\frac{\partial \alpha}{\partial N_{ic}}$ are given
by the sign of the numbers

$$
\sigma_i := (1 - \alpha) \frac{\partial q_{12}}{\partial N_{ic}} - \alpha \frac{\partial q_{21}}{\partial N_{ic}} \quad (i = 1, 2) \qquad (13)
$$
We begin by the dependency with respect to $N_{2c}$. One has first

$$ \frac{\partial q_{12}}{\partial N_{2c}} = \mathcal{R}_2 \eta_2 \frac{(1 - \eta_2)(N_2 + \eta_1 N_{1c}) + \eta_1(1 - \eta_1)N_{1c}}{(N_2 + \eta_1 N_{1c} - \eta_2 N_{2c})^2} > 0 $$

Note that one has

$$ q_{21} = \frac{\mathcal{R}_1(N_2 + \eta_1 N_{1c} - \eta_2 N_{2c})}{\mathcal{R}_2(N_1 - \eta_1 N_{1c} + \eta_2 N_{2c})} q_{12} \quad (14) $$

and thus

$$ \frac{\partial q_{21}}{\partial N_{2c}} = \frac{\mathcal{R}_1(N_2 + \eta_1 N_{1c} - \eta_2 N_{2c}) \partial q_{12}}{\mathcal{R}_2(N_1 - \eta_1 N_{1c} + \eta_2 N_{2c}) \partial N_{2c}} - \frac{\mathcal{R}_1 \eta_2 (N_1 + N_2)}{\mathcal{R}_2 (N_1 - \eta_1 N_{1c} + \eta_2 N_{2c})^2} q_{12} $$

Then, one gets the inequality

$$ \sigma_2 > \left(1 - \alpha - \alpha \frac{\mathcal{R}_1(N_2 + \eta_1 N_{1c} - \eta_2 N_{2c})}{\mathcal{R}_2(N_1 - \eta_1 N_{1c} + \eta_2 N_{2c})} \right) \frac{\partial q_{12}}{\partial N_{2c}} $$

On another hand, one gets from $P(\alpha) = 0$ the inequality

$$ (1 - \alpha)q_{12} - \alpha q_{21} = \alpha(1 - \alpha)(\mathcal{R}_2 - \mathcal{R}_1) > 0 $$

and with (14)

$$ (1 - \alpha)q_{12} - \alpha q_{21} = \left(1 - \alpha - \alpha \frac{\mathcal{R}_1(N_2 + \eta_1 N_{1c} - \eta_2 N_{2c})}{\mathcal{R}_2(N_1 - \eta_1 N_{1c} + \eta_2 N_{2c})}\right) q_{12} > 0 $$

We then conclude that $\sigma_2$ is positive, and from (10) we deduce that the map $N_{2c} \mapsto \mathcal{R}_{1,2}$ is decreasing. This proves the point i.

We study now the dependency with respect to $N_{1c}$. A calculation of the partial derivative gives

$$ \frac{\partial q_{12}}{\partial N_{1c}} = \mathcal{R}_2 \eta_1 \frac{(1 - \eta_1)(N_2 - \eta_2 N_{2c}) - \eta_2(1 - \eta_2)N_{2c}}{(N_2 - \eta_2 N_{2c} + \eta_1 N_{1c})^2} \quad (15) $$

and

$$ \frac{\partial q_{21}}{\partial N_{1c}} = \mathcal{R}_1 \eta_1 \frac{(1 - \eta_1)(N_1 + \eta_2 N_{2c}) + \eta_2(1 - \eta_2)N_{2c}}{(N_1 - \eta_1 N_{1c} + \eta_2 N_{2c})^2} > 0 \quad (16) $$

When $\frac{\partial q_{12}}{\partial N_{1c}} < 0$, we can conclude that $\sigma_1$ is negative and $\mathcal{R}_{1,2}$ is thus increasing with respect to $N_{1c}$. This condition is equivalent to (9). This proves the point ii. When this last condition is not satisfied, having $\frac{\partial q_{12}}{\partial N_{1c}} < \frac{\partial q_{21}}{\partial N_{1c}}$ with $\alpha > \frac{1}{2}$ is another sufficient condition to obtain $\sigma_1 < 0$ from expression (13). However, having $\alpha > \frac{1}{2}$ amounts to have $P(\frac{1}{2}) > 0$, that is

$$ \frac{\mathcal{R}_2 - \mathcal{R}_1}{4} - \frac{\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21}}{2} + q_{12} > 0 $$

or equivalently

$$ \frac{\mathcal{R}_2}{2} - q_{12} < \frac{\mathcal{R}_1}{2} - q_{21} $$

One can check that this last condition is equivalent to $A < 0$ and that the condition $\frac{\partial q_{12}}{\partial N_{1c}} < \frac{\partial q_{21}}{\partial N_{1c}}$ is equivalent to $B < 0$. In the same manner, having $A > 0$ and $B > 0$ implies $\alpha < \frac{1}{2}$ and $\frac{\partial q_{12}}{\partial N_{1c}} > \frac{\partial q_{21}}{\partial N_{1c}}$, which is a sufficient condition to have $\sigma_1 > 0$, and thus $\mathcal{R}_{1,2}$ increasing with respect to $N_{1c}$. This proves the point iii. $\square$
14 The role of permanently resident populations

This result suggests that the map $N_{1c} \mapsto \tilde{\mathcal{R}}_{1,2}(N_{1c}, N_{2c})$ is not necessarily monotonic, differently to the map $N_{2c} \mapsto \tilde{\mathcal{R}}_{1,2}(N_{1c}, N_{2c})$. We show now that the possibilities of its variations are limited.

**Proposition 5** Under hypotheses of Proposition 4, for each $N_{2c} \in (0, N_2)$ the map $N_{1c} \mapsto \tilde{\mathcal{R}}_{1,2}(N_{1c}, N_{2c})$ possesses one of the three properties

1. it is decreasing on $(0, N_1)$,

2. it is increasing on $(0, N_1)$,

3. there exists $N_{1c}^* \in (0, N_1)$ such that it is decreasing on $(0, N_{1c}^*)$ and increasing on $(N_{1c}^*, N_1)$.

*Proof* Fix $N_{2c} \in (0, N_2)$. If the map $N_{1c} \mapsto \tilde{\mathcal{R}}_{1,2}(N_{1c}, N_{2c})$ is not monotonic, there exists $\hat{N}_{1c} \in (0, N_1)$ such that $\frac{\partial \tilde{\mathcal{R}}_{1,2}}{\partial N_{1c}}(\hat{N}_{1c}, N_{2c}) = 0$. For simplicity, we shall drop the notation $\sim$ in the rest of the proof. Following the proof of Proposition 4, one has $\mathcal{R}_{1,2} = \alpha \mathcal{R}_1 + (1 - \alpha) \mathcal{R}_2$ with

$$ \frac{\partial \alpha}{\partial N_{1c}} = \frac{(1 - \alpha) \frac{\partial q_{12}}{\partial N_{1c}} - \alpha \frac{\partial q_{21}}{\partial N_{1c}}}{\mathcal{R}_2 - \mathcal{R}_1 + q_{12} + q_{21} - 2\alpha(\mathcal{R}_2 - \mathcal{R}_1)} := \frac{\sigma_1}{\nu} $$

where $\nu > 0$. Therefore, one has $\frac{\partial \alpha}{\partial N_{1c}} = 0$ and $\sigma_1 = 0$ at $N_{1c} = \hat{N}_{1c}$, and thus

$$ \left. \frac{\partial^2 \alpha}{\partial N_{1c}^2} \right|_{N_{1c}=\hat{N}_{1c}} = \left. \frac{\frac{\partial \sigma_1}{\partial N_{1c}}}{\nu} \right|_{N_{1c}=\hat{N}_{1c}} = \left. \frac{(1-\alpha) \frac{\partial^2 q_{12}}{\partial N_{1c}^2} - \alpha \frac{\partial^2 q_{21}}{\partial N_{1c}^2}}{\nu} \right|_{N_{1c}=\hat{N}_{1c}} $$

From expressions (15) and (16), a calculation of the partial derivatives gives

$$ \frac{\partial^2 q_{12}}{\partial N_{1c}^2} = \frac{-2\eta_1 \frac{\partial q_{12}}{\partial N_{1c}}}{N_1 - \eta_1 N_{1c} + \eta_2 N_{2c}}, \quad \frac{\partial^2 q_{21}}{\partial N_{1c}^2} = \frac{2\eta_1 \frac{\partial q_{21}}{\partial N_{1c}}}{N_2 - \eta_2 N_{2c} + \eta_1 N_{1c}} $$

where $\frac{\partial q_{21}}{\partial N_{1c}} > 0$ and from $\sigma_1 = 0$ one gets $\frac{\partial q_{12}}{\partial N_{1c}} > 0$ for $N_{1c} = \hat{N}_{1c}$. Finally, one obtains

$$ \frac{\partial^2 \mathcal{R}_{1,2}}{\partial N_{1c}^2}(\hat{N}_{1c}, N_{2c}) = -(\mathcal{R}_2 - \mathcal{R}_1) \frac{\partial^2 \alpha}{\partial N_{1c}^2}(\hat{N}_{1c}, N_{2c}) < 0 $$

Consequently, any extremum of the map $N_{1c} \mapsto \mathcal{R}_{1,2}(N_{1c}, N_{2c})$ is a local minimizer, which implies that this map has at most one local minimizer. $\square$

Finally, we give conditions for which the minimization of the threshold $\mathcal{R}_{1,2}$ presents a trichotomy.

**Proposition 6** Let parameters $\beta_i, \gamma$ be such that $\mathcal{R}_2 > \mathcal{R}_1$ and assume that $N_1, N_2$ satisfy $N_1 \mathcal{R}_2 > N_2 \mathcal{R}_1$. Then, provided that $\gamma$ is small enough compared to $\lambda_i$ and $\mu_i$, the function $(N_{1c}, N_{2c}) \mapsto \mathcal{R}_{1,2}(N_{1c}, N_{2c})$ admits an unique minimum at $(N_{1c}^*, N_{2c}^*)$ with $N_{2c}^* = N_2$. Moreover, one has the following properties.

1. $N_{1c}^* = 0$ if $\eta_2 > 1 - \eta_1$,

2. $N_{1c}^* = N_1$ if $\eta_1$ and $\eta_2$ are sufficiently small,
3. there exists $\eta_1, \eta_2$ for which $N_{1c}^* \in (0, N_1)$.

*Proof* We first show that the announced properties are satisfied for the approximate function $\tilde{\mathcal{R}}_{1,2}$.

From Propositions 4 and 5, we know that $\tilde{\mathcal{R}}_{1,2}$ admits a unique minimum at $(\tilde{N}_{1c}, \tilde{N}_{2c})$ with $\tilde{N}_{2c} = N_2$. For $N_{2c} = N_2$, the condition (9) simply writes $\eta_2 > 1 - \eta_1$ which implies from point ii. of Proposition 4 that one has $\tilde{N}_{1c} = 0$ when this condition is fulfilled. This shows that point 1 is verified for the function $\tilde{\mathcal{R}}_{1,2}$.

One obtains the limits

$$ \lim_{\eta_1, \eta_2 \to 0} A = \frac{\mathcal{R}_2}{2} - \frac{\mathcal{R}_1}{2} > 0, \quad \lim_{\eta_1, \eta_2 \to 0} B = \frac{\mathcal{R}_2}{N_2} - \frac{\mathcal{R}_1}{N_1} > 0 $$

which show that numbers A and B are positive when $\eta_1, \eta_2$ are small, and thus one has $\tilde{N}_{1c} = N_1$ from point iii of Proposition 4. This shows that point 2 is verified for the function $\tilde{\mathcal{R}}_{1,2}$.

Take now any $N_{1c} \in (0, N_1)$. When $\eta_2 > 1 - \eta_1$, one has $\frac{\partial \tilde{\mathcal{R}}_{1,2}}{\partial N_{1c}}(N_{1c}, N_2) > 0$, and for $\eta_1, \eta_2$ small, $\frac{\partial \tilde{\mathcal{R}}_{1,2}}{\partial N_{1c}}(N_{1c}, N_2) < 0$ is verified. Then, by continuity of the function $\tilde{\mathcal{R}}_{1,2}$ with respect to parameters $\eta_1, \eta_2$, one deduce that the existence of values $\hat{\eta}_1, \hat{\eta}_2$ for which $\frac{\partial \tilde{\mathcal{R}}_{1,2}}{\partial N_1}(N_{1c}, N_2) = 0$. As the function $\tilde{\mathcal{R}}_{1,2}$ cannot have more than a local extremum (see Proposition 5), we deduce that $N_{1c}$ realizes the minimum of the function $N_{1c} \mapsto \tilde{\mathcal{R}}_{1,2}(N_{1c}, N_2)$ when $\eta_1 = \hat{\eta}_1$ and $\eta_2 = \hat{\eta}_2$. This shows that point 3 is verified for the function $\tilde{\mathcal{R}}_{1,2}$.

Finally, note that the exact threshold $\mathcal{R}_{1,2}$ amounts to replace in the expression of $\tilde{q}_{12}, \tilde{q}_{21}$ the numbers $\eta_i$ by $\frac{\lambda_i + \gamma}{\lambda_i + \mu_i + \gamma}$, which is continuous with respect to $\gamma$ and equal to $\eta_i$ for $\gamma = 0$. By continuity of $\tilde{\mathcal{R}}_{1,2}$ with respect to $\tilde{q}_{12}, \tilde{q}_{21}$, we deduce that uniqueness of the minimizer of $\mathcal{R}_{1,2}$ and properties 1. to 3. are also fulfilled by the function $(N_{1c}, N_{2c}) \mapsto \mathcal{R}_{1,2}$, provided that $\gamma$ is small enough. $\square$

# 5 Numerical illustration

We consider two territories of same population size $N = N_1 = N_2$ with different transmission rates such that one has $\mathcal{R}_1 < 1 < \mathcal{R}_2$ (values are given in Table 1). Typically, some precautionary measures (such as social distance) are taken in the first territory so that the disease cannot spread in this territory if it is closed, while the epidemic can spread in the second territory in absence of communication with territory 1. We aim at studying how the epidemic can die out when commuting occur between territories, depending on the proportions of resident in each population, denoted

$$ p_i := \frac{N_{ir}}{N} = 1 - \frac{N_{ic}}{N}, \quad (i = 1, 2) $$

(in other words, how to obtain $\mathcal{R}_{1,2} < 1$ playing with $p_1, p_2$). Note that when $N_1 = N_2$, the threshold $\mathcal{R}_{1,2}$ depends on the proportions $p_1, p_2$ independently of $N$.
16 *The role of permanently resident populations*

<table><caption>**Table 1** Characteristics numbers of the epidemic</caption><thead><tr><th>γ</th><th>β<sub>1</sub></th><th>β<sub>2</sub></th><th>R<sub>1</sub></th><th>R<sub>2</sub></th></tr></thead><tbody><tr><td>0.3</td><td>0.27</td><td>0.33</td><td>0.9</td><td>1.1</td></tr></tbody></table>

Conditions of Proposition 6 are satisfied provided that commuting parameters $λ_i$, $μ_i$ are large enough. We have considered three sets of these parameters, given in Table 2, that correspond to the three possible situations depicted in Proposition 6.

<table><caption>**Table 2** Three sets of commuting parameters</caption><thead><tr><th>case</th><th>λ<sub>1</sub></th><th>μ<sub>1</sub></th><th>λ<sub>2</sub></th><th>μ<sub>2</sub></th><th>η<sub>1</sub></th><th>η<sub>2</sub></th></tr></thead><tbody><tr><td>A</td><td>10</td><td>10</td><td>10</td><td>1</td><td>0.5</td><td>0.9090909</td></tr><tr><td>B</td><td>10</td><td>100</td><td>10</td><td>100</td><td>0.009901</td><td>0.009901</td></tr><tr><td>C</td><td>10</td><td>10</td><td>10</td><td>70</td><td>0.5</td><td>0.125</td></tr></tbody></table>

The approximate expression $\tilde{R}_{1,2}$ turns out to be a very good approximation of the exact value $R_{1,2}$, even in case A for which $γ$ is not so small compared to $μ_2$ (see Table 3).

<table><caption>**Table 3** Quality of the approximation $\tilde{R}_{1,2}$</caption><thead><tr><th>case</th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td>max<sub>p<sub>1</sub>,p<sub>2</sub></sub> |R&#x0303;<sub>1,2</sub> - R<sub>1,2</sub>|</td><td>1.9 10<sup>-3</sup></td><td>1.4 10<sup>-4</sup></td><td>6 10<sup>-4</sup></td></tr></tbody></table>

Figures 1, 2, 3 show families of curves $p_1 \mapsto R_{1,2}$ for different values of $p_2 \in [0, 1]$. One can observe that these curves possess the properties given by Propositions 4 and 5:

- they are either decreasing, increasing or decreasing down to a minimum and then increasing,

- they are ordered and the lower one is obtained for $p_2 = 0$ (i.e. $N_{2c} = N_2$).

This last feature is intuitive: the more there are commuters from territory 2 (that spend time in territory 1 where the conditions of transmission disease is lower), the less the epidemic spreads. A way to reduce the value of $R_{1,2}$ is thus to encourage commuting towards territory 1 (whatever are the commuting rates). However, the role of the resident population in territory 1 is far less intuitive because it does depend on the commuting rates.

1. In case A, commuters from territory 2 return more rarely to home than commuters from territory 1 do. The condition of point 1 of Proposition 6 is fulfilled. Then, the threshold $R_{1,2}$ can be made small (and below 1) when the proportion of resident in territory 1 is high i.e. when the inhabitants of territory 1 are encouraged not to commute.
2. In case B, both commuters return rapidly to their home. This means that the numbers of commuters from one territory present in the other one at a given time is low. Then the condition of point 2 of Proposition 6 is fulfilled. Here, it is better to encourage inhabitants of territory 1 to commute to the other territory where the disease spreads yet more easily... which is counter-intuitive at first sight. Indeed, commuters do not spend much time in the other territory, and therefore heuristically have less time to encounter and transmit the disease...

3. In case C, commuters from territory 2 return more rapidly to home than commuters from territory 1 do, on the opposite of case A. Conditions of points 1 and 2 of Proposition 6 are not fulfilled here and we are in an intermediate situation for which point 3 of Proposition 6 occurs. It is theoretically possible to have $\mathcal{R}_{1,2} < 1$ on the condition that the proportion of commuters of territory 1 is well balanced.

Finally, this example shows that changing only the return rates $\mu_1, \mu_2$ allows to obtain the three possible scenarios, but other changes could also exhibit them.

<figure><img src="image_2.png" /><figcaption>Fig. 1 $\mathcal{R}_{1,2}$ as a function of $p_1$ in case A (each curve corresponds to a value of $p_2 \in [0, 1]$)</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** $\mathcal{R}_{1,2}$ as a function of $p_1$ in case B (each curve corresponds to a value of $p_2 \in [0, 1]$)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** $\mathcal{R}_{1,2}$ as a function of $p_1$ in case C (each curve corresponds to a value of $p_2 \in [0, 1]$)</figcaption></figure>
# 6 Conclusion

In this work, we have been able to provide an explicit expression of the reproduction number, although the model is in dimension 18. This expression has allowed us to study its minimization with respect to the proportions of permanently resident populations in each patch. We discovered a trichotomy of cases, with some counter intuitive situations. In each case, it is always beneficial to have commuters traveling to a safer city where the transmission rate is lower. However, for the safer city, three situations occurs:

- either it is better to avoid commuting to the other city,

- or on the opposite encouraging commuting to the more risky city reduces the reproduction number,

- and in a third case there exists an optimal intermediate proportion of commuters of the safer city which minimizes the epidemic threshold.

In some sense, the permanently resident populations, which have been ignored in former modeling, can play an hidden role in an epidemic outbreak. This is illustrated on an example for which only right proportions of commuters (or permanently resident) avoid the outbreak. This suggests that counter-intuitive situations may also occur when considering networks with more than two nodes. The present study focuses on the reproduction number and how it can be reduced. The impacts of resident proportions on other epidemiological characteristics, such as the peak level or the finite size, may be the matter a future work. The extension of the present results to more general networks is also a future perspective.

# Acknowledgments

The authors are grateful for the support of the French platform MODCOV19, and the Algerian Government for the PhD grant of Ismail Mimouni.

The authors thank the anonymous referee to let us know the alternative approach to obtain the expression of the reproduction number, mentioned in Remark 1.

# References

[1] ANDERSON, R.M. AND MAY, R.M. *Infectious Diseases of Humans*, Oxford University, 1991.

[2] ANGULO, J.J, TAKIGUTI, C.K, PEDERNEIRAS, C.A, CARVALHO-DE-SOUZA, A.M, OLIVEIRA-DE-SOUZA, M.C. AND MEGALE, P. *Identification of pattern and process in the spread of a contagious disease*. Social Science & Medicine, 13D(3):183-9, 1979.

[3] ARINO, J., JORDAN, R. AND VAN DEN DRIESSCHE, P. *Quarantine in a multi-species epidemic model with spatial dynamics*. Mathematical
20 *The role of permanently resident populations*

Biosciences, 206(1), 46–60, 2007.

[4] ARINO, J. *Diseases in Metapopulations*. In Modeling and Dynamics of Infectious Diseases, World Scientific, 64–122, 2009.

[5] ARINO, J AND VAN DEN DRIESSCHE, P. *A multi-city epidemic model*. Mathematical Population Studies, 10(3), 175–193, 2003.

[6] ARINO, J AND VAN DEN DRIESSCHE, P. *Disease spread in metapopulations*. Fields Institute Communications, 48, 2006.

[7] BERMAN, A. AND PLEMMONS, R.J. *Nonnegative Matrices in the Mathematical Sciences*, SIAM, 1994.

[8] BRUNHAM, R. *Core group theory: a central concept in STD epidemiology*. Venereology, 10(1), 34–5, 8–9, 1997.

[9] CITRON, D.T., GUERRA, C.A., DOLGERT, A.J., WU, S.L., HENRY, J.M., SÁNCHEZ C, H. M. AND SMITH, D.L. *Comparing metapopulation dynamics of infectious diseases under different models of human movement*. Proceedings of the National Academy of Sciences, 118(18), e2007488118, 2021.

[10] CHEN, Y., YAN, M., AND XIANG, Z., *Transmission Dynamics of a Two-City SIR Epidemic Model with Transport-Related Infections*. Journal of Applied Mathematics, ID 764278, 2014.

[11] CLANCY, D. *Carrier-borne epidemic models incorporating population mobility*. Mathematical Biosciences, 132(2), 185–204, 1996.

[12] COLIZZA, V., BARRAT, A., BARTHELEMY, M. AND VESPIGNANI, A. *The role of the airline transportation network in the prediction and predictability of global epidemics*. Proc. Natl. Acad. Sci. U.S.A. 103, 2015–2020, 2006.

[13] DIEKMANN, O., HEESTERBEEK, H. AND BRITTON T. *Mathematical Tools for Understanding Infectious Disease Dynamics*. Princeton University Press, 2013.

[14] DIEKMANN, O., HEESTERBEEK, J.A.P. AND METZ, J.A.J. *On the definition and computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations*. Journal of Mathematical Biology, 28, 365–382 1990.

[15] DIEKMANN, O., HEESTERBEEK, J.A.P. AND ROBERT M.G., *The construction of next-generation matrices for compartmental epidemic models*. Journal of the Royal Society Interface, 7, 873–885, 2007.
[16] DHIRASAKDANON, T., THIEME, H.R. AND VAN DEN DRIESSCHE, P. *A sharp threshold for disease persistence in host metapopulations*. Journal of Biological Dynamics, 1(4), 363–378, 2007.

[17] GAO, D. *How Does Dispersal Affect the Infection Size?* SIAM Journal on Applied Mathematics, 80, 2144–2169, 2007.

[18] GUO, H., LI, M. AND SHUAI, Z. *Global stability of the endemic equilibrium of multigroup SIR epidemic models*. Canadian Applied Mathematics Quarterly, 14(3), 259–284, 2006.

[19] HADELER, K.P. AND CASTILLO-CHAVEZ, C. *A core group model for disease transmission*, Mathematical Biosciences, 128 (1-2), 41–55, 1995.

[20] HANSKI, I. *Metapopulation Ecology*. Oxford University Press, 1999.

[21] IGGIDR, A., SALLET, G. AND TSANOU, B. *Global Stability Analysis of a Metapopulation SIS Epidemic Model*. Mathematical Population Studies, 19(3), 115–129, 2012.

[22] KEELING, M. J., BJØRNSTAD, O. N. AND GRENFELL, B. T. *Metapopulation dynamics of infectious diseases*. Chapter 17 in I. Hanski, and O. Gaggiotti (eds.). Ecology, Genetics, and Evolution of Metapopulations, Elsevier, 415–445, 2004.

[23] KEELING, M.J. AND ROHANI, P. *Estimating Spatial Coupling in Epidemiological Systems: a Mechanistic Approach*. Ecology Letters, 5, 20–29, 2002.

[24] KEELING, M.J. AND ROHANI, P. *Modeling Infectious Diseases in Humans and Animals*. Princeton University Press, 2007.

[25] KELLY, M.R., TIEN, J.H., EISENBERG, M.C AND LENHART, S. *The impact of spatial arrangements on epidemic disease dynamics and intervention strategies*. Journal of Biological Dynamics, 10, 222–49, 2016.

[26] KERMACK, W. AND MCKENDRICK, A. *A contribution to the mathematical theory of epidemics*. Proceedings of the Royal Society, A 115, 700–721, 1927.

[27] KNIPL, D. *A new approach for designing disease intervention strategies in metapopulation models*. Journal of Biological Dynamics, 10(1), 71–94, 2016.

[28] LI, J., XIANG, T. AND HE, L. *Modeling epidemic spread in transportation networks: A review*. Journal of Traffic and Transportation Engineering, 8(2), 139–152, 2021.
[29] LIPSHTAT, A., ALIMI, R. AND BEN-HORIN, Y. *Commuting in metapopulation epidemic modeling*. Scientific Reports, 11, 15198, 2021.

[30] LIU, X. AND STECHLINSKI, P. *Transmission dynamics of a switched multi-city model with transport-related infections*. Nonlinear Analysis: Real World Applications, 14(1), 264–279, 2013.

[31] MACARTHUR, R.H. AND AND WILSON, E.O. *The Theory of Island Biogeography*, Princeton University Press, 2001.

[32] MPOLYA, E.A., YASHIMA, K., OHTSUKI, H AND SASAKI, A. *Epidemic dynamics of a vector-borne disease on a villages-and-city star network with commuters*. Journal of Theoretical Biology, 343, 120–26, 2014.

[33] MURRAY, J.D. *Mathematical Biology II: Spatial Models and Biomedical Applications*. Springer, 2003.

[34] PATIL, R., DAVE, R., PATEL, H., SHAH, V.M., CHAKRABARTI, D. AND BHATIA, U. *Assessing the interplay between travel patterns and SARS-CoV-2 outbreak in realistic urban setting*. Applied Network Science 6, 4, 2021.

[35] RUAN, Z., WANG, C., MING HUI, P. AND LIU, Z. *Integrated travel network model for studying epidemics: Interplay between journeys and epidemic*. Scientific Reports 5, 11401, 2015.

[36] SATTENSPIEL, L. AND DIETZ, K. *A structured epidemic model incorporating geographic mobility among regions*. Mathematical Biosciences, 128 (1-2), 71–91, 1995.

[37] TAKEUCHI, Y., LIU, X. AND CUI, J. *Global dynamics of SIS models with transport-related infection*. Journal of Mathematical Analysis and Applications, 329(2), 1460-71, 2007.

[38] TOCTO-ERAZO, M.R., OLMOS-LICEAGA, D. AND MONTOYA-LAOS, J.A. *Effect of daily human movement on some characteristics of dengue dynamics*. Mathematical Biosciences, 332, 108531, 2021.

[39] VAN DEN DRIESSCHE, P. AND WATMOUGH, J. *Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission*. Mathematical Biosciences, 180 (1-2), 29–48, 2002.

[40] VOROTNIKOV, V.I *Partial Stability and Control*, Birkhäuser, 1998.

[41] WANG, W. AND MULONE, G. *Threshold of disease transmission in a patch environment*. Journal of Mathematical Analysis and Applications, 285(1), 321–335, 2003.
[42] WANG, W. AND ZHAO, X.Q. *An epidemic model in a patchy environment*. Mathematical Biosciences,190(1), 97-112, 2004.

[43] YIN, Q., WANG, Z., XIA, C., DEHMER, M., EMMERT-STREIB, F. AND JIN, Z. *A novel epidemic model considering demographics and intercity commuting on complex dynamical networks*. Applied Mathematics and Computation, 386, N. 125517, 2020.