Multi-Patch and Multi-Group Epidemic Models: A New Framework

TL;DR

This paper introduces a comprehensive multi-patch, multi-group epidemic model incorporating host mobility based on epidemiological status, deriving the basic reproduction number and analyzing how heterogeneity and mobility influence disease spread.

Contribution

It presents a novel framework for modeling complex host interactions and mobility patterns, including explicit derivation of the basic reproduction number in this context.

Findings

Reproduction number increases with the number of patches and mobility complexity.

Explicit formula for $ ext{R}_0$ when mobility matrices are rank one.

Mobility patterns significantly impact disease prevalence and $ ext{R}_0$.

Abstract

We develop a multi-patch and multi-group model that captures the dynamics of an infectious disease when the host is structured into an arbitrary number of groups and interacts into an arbitrary number of patches where the infection takes place. In this framework, we model host mobility that depends on its epidemiological status, by a Lagrangian approach. This framework is applied to a general SEIRS model and the basic reproduction number $\mathcal{R_0}$ is derived. The effects of heterogeneity in groups, patches and mobility patterns on $\mathcal{R_0}$ and disease prevalence are explored. Our results show that for a fixed number of groups, the basic reproduction number increases with respect to the number of patches and the host mobility patterns. Moreover, when the mobility matrix of susceptible individuals is of rank one, the basic reproduction number is explicitly determined and was found to be independent of the latter. The cases where mobility matrices are of rank one capture important modeling scenarios. Additionally, we study the global analysis of equilibria for some special cases. Numerical simulations are carried out to showcase the ramifications of mobility pattern matrices on disease prevalence and basic reproduction number.

Tables

Table 1. Table 1: Description of the parameters used in System ( 1 ).

Parameters	Description
${\Lambda }_i$	Recruitment of the susceptible individuals in Group $i$
${\beta }_j$	Instantaneous risk of infection in Patch $j$
${\mu }_i$	Per capita natural death rate of Group $i$
${\nu }_i$	Per capita rate at which latent in Group $i$ become infectious
${\gamma }_i$	Per capita recovery rate of Group $i$
$m_{ij}$	Proportion of time susceptible individuals of Group $i$ spend in Patch $j$
$n_{ij}$	Proportion of time latent individuals of Group $i$ spend in Patch $j$
$p_{ij}$	Proportion of time infectious individuals of Group $i$ spend in Patch $j$
$q_{ij}$	Proportion of time recovered individuals of Group $i$ spend in Patch $j$
${\eta }_i$	Per capita loss of immunity rate
${\delta }_i$	Per capita disease induced death rate of Group $i$.

Equations

$$\frac{{\sum }_{k=1}^u{p}_{kj}{I}_k}{{\sum }_{k=1}^u(m_{kj}{S}_k+n_{kj}{E}_k+p_{kj}{I}_k+q_{kj}{R}_k)}$$

$${\dot{S}}_i={\Lambda }_i-\sum _{j=1}^v{\beta }_j{m}_{ij}{S}_i\frac{{\sum }_{k=1}^u{p}_{kj}{I}_k}{{\sum }_{k=1}^u(m_{kj}{S}_k+n_{kj}{E}_k+p_{kj}{I}_k+q_{kj}{R}_k)}-{\mu }_i{S}_i+{\eta }_i{R}_i$$

$${\dot{E}}_i=\sum _{j=1}^v{\beta }_j{m}_{ij}{S}_i\frac{{\sum }_{k=1}^u{p}_{kj}{I}_k}{{\sum }_{k=1}^u(m_{kj}{S}_k+n_{kj}{E}_k+p_{kj}{I}_k+q_{kj}{R}_k)}-({\nu }_i+{\mu }_i)E_i$$

$${\dot{I}}_i={\nu }_i{E}_i-({\gamma }_i+{\mu }_i)I_i$$

$${\dot{R}}_i={\gamma }_i{I}_i-({\eta }_i+{\mu }_i)R_i$$

$$\left\{\begin{array}[]{llll}\dot{S}_{i}=\Lambda_{i}-\sum_{j=1}^{v}\beta_{j}m_{ij}S_{i}\dfrac{\sum_{k=1}^{u}p_{kj}I_{k}}{\sum_{k=1}^{u}(m_{kj}S_{k}+n_{kj}E_{k}+p_{kj}I_{k}+q_{kj}R_{k})}-\mu_{i}S_{i}+\eta_{i}R_{i},\\[8.53581pt] \dot{E}_{i}=\sum_{j=1}^{v}\beta_{j}m_{ij}S_{i}\dfrac{\sum_{k=1}^{u}p_{kj}I_{k}}{\sum_{k=1}^{u}(m_{kj}S_{k}+n_{kj}E_{k}+p_{kj}I_{k}+q_{kj}R_{k})}-(\nu_{i}+\mu_{i})E_{i}\\ \dot{I}_{i}=\nu_{i}E_{i}-(\gamma_{i}+\mu_{i}+\delta_{i})I_{i}\\ \dot{R}_{i}=\gamma_{i}I_{i}-(\eta_{i}+\mu_{i})R_{i}\end{array}\right.$$

$$\left\{\begin{array}[]{llll}\dot{\mathbf{S}}={\mathbf{\Lambda}}-\textrm{diag}(\mathbf{S})\mathbb{M}\textrm{diag}(\mathcal{B})\textrm{diag}^{-1}(\mathbb{M}^{T}\mathbf{S}+\mathbb{N}^{T}\mathbf{E}+\mathbb{P}^{T}\mathbf{I}+\mathbb{Q}^{T}\mathbf{R})\mathbb{P}^{T}\mathbf{I}-\textrm{diag}(\mu)\mathbf{S}+\textrm{diag}(\eta)\mathbf{R}\\ \dot{\mathbf{E}}=\textrm{diag}(\mathbf{S})\mathbb{M}\textrm{diag}(\mathcal{B})\textrm{diag}^{-1}(\mathbb{M}^{T}\mathbf{S}+\mathbb{N}^{T}\mathbf{E}+\mathbb{P}^{T}\mathbf{I}+\mathbb{Q}^{T}\mathbf{R})\mathbb{P}^{T}\mathbf{I}-\textrm{diag}(\nu+\mu)\mathbf{E}\\ \dot{\mathbf{I}}=\textrm{diag}(\nu)\mathbf{E}-\textrm{diag}(\gamma+\mu+\delta)\mathbf{I}\\ \dot{\mathbf{R}}=\textrm{diag}(\gamma)\mathbf{I}-\textrm{diag}(\eta+\mu)\mathbf{R}\end{array}\right.$$

$${\dot{S}}_i={\Lambda }_i-\sum _{j=1}^v{\beta }_j{m}_{ij}{g}_i{S}_i\frac{{\sum }_{k=1}^u{p}_{kj}{I}_k}{{\sum }_{k=1}^u(m_{kj}{S}_k+n_{kj}{E}_k+p_{kj}{I}_k+q_{kj}{R}_k)}-{\mu }_i{S}_i+{\eta }_i{R}_i,$$

$${\dot{E}}_i=\sum _{j=1}^v{\beta }_j{m}_{ij}{g}_i{S}_i\frac{{\sum }_{k=1}^u{p}_{kj}{I}_k}{{\sum }_{k=1}^u(m_{kj}{S}_k+n_{kj}{E}_k+p_{kj}{I}_k+q_{kj}{R}_k)}-({\nu }_i+{\mu }_i)E_i.$$

$$\dot{\mathbf{N}}=\mathbf{\Lambda }-\mu \circ \mathbf{N}-\delta \circ \mathbf{I}\le \mathbf{\Lambda }-\mu \circ \mathbf{N}$$

$$\Omega=\left\{(\mathbf{S},\mathbf{E},\mathbf{I},\mathbf{R})\in{{\rm I\mkern-3.5muR}}^{4u}_{+}\;\mbox{ \LARGE$\mid$ }\;\mathbf{S}+\mathbf{E}+\mathbf{I}+\mathbf{R}\leq\mathbf{\Lambda}\circ\dfrac{1}{\mu}\right\}$$

$$\displaystyle\left(\begin{array}[]{c}\dot{\mathbf{E}}\\ \dot{\mathbf{I}}\end{array}\right)$$

$$=$$

$$F=D\mathcal{F}(\mathbf{E},\mathbf{I})\Bigg{|}_{\textrm{DFE}}=\left(\begin{array}[]{cccc}\textbf{0}_{u,u}&\textrm{diag}(\mathbf{S^{\ast}})\mathbb{M}\textrm{diag}(\mathcal{B})\textrm{diag}^{-1}(\mathbb{M}^{T}\mathbf{S^{\ast}})\mathbb{P}^{T}\\ \textbf{0}_{u,u}&\textbf{0}_{u,u}\end{array}\right)$$

$$V=D\mathcal{V}(\mathbf{E},\mathbf{I})\Bigg{|}_{\textrm{DFE}}=\left(\begin{array}[]{cccc}-\textrm{diag}(\mu+\nu)&\textbf{0}_{u,u}\\ \textrm{diag}(\nu)&-\textrm{diag}(\mu+\gamma+\delta)\end{array}\right)$$

$$-V^{-1}=\left(\begin{array}[]{cccc}\textrm{diag}^{-1}(\mu+\nu)&\textbf{0}_{u,u}\\ \textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma+\delta))&\textrm{diag}^{-1}(\mu+\gamma+\delta)\end{array}\right)$$

$$-FV^{-1}=\left(\begin{array}[]{cccc}Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma+\delta))&Z\textrm{diag}^{-1}(\mu+\gamma+\delta)\\ \textbf{0}_{u,u}&\textbf{0}_{u,u}\end{array}\right)$$

$$Z=\text{diag}({\mathbf{S}}^{\ast })\mathbb{M}\text{diag}(\mathcal{B}){\text{diag}}^{-1}({\mathbb{M}}^T{\mathbf{S}}^{\ast }){\mathbb{P}}^T$$

$${\mathcal{R}}_0(u,v)=\rho (Z\text{diag}(\nu ){\text{diag}}^{-1}((\mu +\nu )\circ (\mu +\gamma +\delta )))$$

$${\mathcal{R}}_0(u,v,\mathbb{P})$$

$${\tilde{\mathcal{R}}}_0^i(u,v)$$

$${\mathcal{R}}_0^i=\frac{{\nu }_i}{({\mu }_i+{\nu }_i)({\mu }_i+{\gamma }_i+{\delta }_i)}\sum _{k=1}^v{\beta }_k{p}_{ik}$$

$$e_i^T\text{diag}({\mathbf{S}}^{\ast })\mathbb{M}=(m_{i1}{S}_i^{\ast },m_{i2}{S}_i^{\ast },\dots ,m_{iv}{S}_i^{\ast })$$

$$e_i^T\text{diag}({\mathbf{S}}^{\ast })\mathbb{M}\text{diag}(\mathcal{B})=({\beta }_1{m}_{i1}{S}_i^{\ast },{\beta }_2{m}_{i2}{S}_i^{\ast },\dots ,{\beta }_v{m}_{iv}{S}_i^{\ast })$$

$$\mathbb{M}^{T}\mathbf{S}^{\ast}=\left(\begin{array}[]{c}\sum_{k=1}^{u}m_{k1}S_{k}^{\ast}\\ \sum_{k=1}^{u}m_{k2}S_{k}^{\ast}\\ \vdots\\ \sum_{k=1}^{u}m_{kv}S_{k}^{\ast}\end{array}\right)$$

$$\textrm{diag}^{-1}(\mathbb{M}^{T}\mathbf{S}^{\ast})\mathbb{P}^{T}e_{i}=\left(\begin{array}[]{c}\frac{p_{i1}}{\sum_{k=1}^{u}m_{k1}S_{k}^{\ast}}\\ \frac{p_{i2}}{\sum_{k=1}^{u}m_{k2}S_{k}^{\ast}}\\ \vdots\\ \frac{p_{iv}}{\sum_{k=1}^{u}m_{kv}S_{k}^{\ast}}\end{array}\right)$$

$$e_i^T\text{diag}({\mathbf{S}}^{\ast })\mathbb{M}\text{diag}(\mathcal{B}){\text{diag}}^{-1}({\mathbb{M}}^T{\mathbf{S}}^{\ast }){\mathbb{P}}^T{e}_i$$

$${\mathcal{R}}_0(u,v)\ge \underset{i=1,\dots ,u}{\max }{\tilde{\mathcal{R}}}_0^i(u,v)$$

$$Z=\text{diag}({\mathbf{S}}^{\ast })\mathbb{M}\text{diag}(\mathcal{B}){\text{diag}}^{-1}({\mathbb{M}}^T{\mathbf{S}}^{\ast }){\mathbb{P}}^T$$

$$z_{ij}=\frac{{\nu }_j}{({\mu }_j+{\nu }_j)({\mu }_j+{\gamma }_j+{\delta }_j)}\sum _{k=1}^v\frac{{\beta }_k{m}_{ik}{p}_{jk}{S}_i^{\ast }}{{\sum }_{l=1}^u{m}_{lk}{S}_l^{\ast }}\forall 1\le i,j\le u.$$

$$\underset{j}{\min }{z}_j\le {\mathcal{R}}_0(u,v)\le \underset{j}{\max }{z}_j$$

Full text

Mathematics Subject Classification: 92D25, 92D30

Keywords:

Multi-Patch, Multi-Group, Mobility, Heterogeneity, Residence Times, Global Stability.

1 Introduction

The role of heterogeneity in populations and their mobility have long been recognized as driving forces in the spread of infectious diseases [1, 17, 37, 43]. Indeed, populations are composed of individuals with different immunological features and hence differ in how they can transmit or acquire an infection at a given time. These differences could result from demographic, host genetic or socio-economic factors [1]. Populations also move across different geographical landscapes, importing their disease history with them either by infecting or getting infected in the host/visiting location.

While the concept of modeling epidemiological heterogeneity within a population goes back to Kermack and McKendrick in modeling the age of infection [31], the approach gained prominence with Yorke and Lajmonivich’s seminal paper [32] on the spread of gonorrhea, a sexually transmitted disease. An abundant and varied literature have followed on understanding the effects of “superspreaders ” which are core groups on the disease dynamics [11, 13, 29, 28, 47] or related multi-group models [10, 20, 22, 25, 35, 40, 43] (and the references therein). Similarly, spatial heterogeneity in epidemiology has been extensively explored in different settings. Continuum models of dispersal have been investigated through diffusion equations [34] whereas islands models have been dealt through metapopulation approach [2, 3, 4, 26, 27, 41, 42], defined here as continuous models with discrete dispersal.

Although the importance and the complete or partial analysis of these two types of heterogeneities have been studied separately in the aforementioned papers, little attention has been given to the simultaneous consideration of groups and spacial heterogeneities. Moreover, previous studies on multi-group rely on differential susceptibility in each group through the WAIFW (Who Acquires Infection From Whom [1]) matrices which, we argue, are difficult to quantify. Similarly, in metapopulation (Eulerian) settings, the movement of individuals between patches is captured in terms of flux of population, making it nearly impossible to track the life-history of individuals after the interpatch mixing.

In this paper, we introduce a general modeling framework that structures populations into an arbitrary number of groups (e.g. demographic, ethnic or socio-economic grouping). These populations, with different health statuses, spend certain amounts of time in an arbitrary number of locations, or patches, where they could get infected or infect others. Each patch is defined by a particular risk of infection tied to environmental conditions of each patch. This approach allows us to track individuals of each group over time and to avoid the use of differential susceptibility of individuals or groups, which is theoretically nice but practically difficult to assess. The likelihood of infection depends both on the time one spends (in a particular patch) and the risk associated with that patch. Moreover, we incorporate individuals’ behavioral decisions through differential residence times. Indeed, individuals of the same group spend different amounts of time in different areas depending on their epidemiological conditions. We also considered two cases of the general framework, that are particularly important from modeling standpoint: when the susceptible and/or infected individuals of different groups have proportional residence times in different patches. That is, when the mobility matrix of susceptible (or infected) individuals, $\mathbb{M}$ (or $\mathbb{P}$) is of rank one. In these cases, we obtain explicit expressions of the basic reproduction number in terms of mobility patterns. It turns out that if $\mathbb{M}$ is of rank one and stochastic, the basic reproduction number is independent of the mobility patterns of susceptible host.

In short, we address how group heterogeneity, or groupness, patch heterogeneity, or patchiness, mobility patterns and behavior each alter or mitigate disease dynamics. In this sense, our paper is a direct extension of [7, 8, 9, 12] but also other studies that capture dispersal through Lagrangian approaches – in which it is possible to track host movement after the interpatch mixing – [15, 26, 38, 39] and a recent paper [19] that investigates the effects of daily movements in the context of Dengue.

The paper is organized as follows. Section 2 explains the model derivation, states the basic properties and the computation of the basic reproduction number $\mathcal{R}_{0}(u,v)$ for $u$ groups and $v$ patches. Section 3 investigates the role of patch and group heterogeneity on the basic reproduction number, and how dispersal patterns alter $\mathcal{R}_{0}(u,v)$ and the disease prevalence. Section 5 is devoted to the existence, uniqueness and stability of equilibria for the considered system under certain conditions. Finally, Section 6 is dedicated to concluding remarks and discussions.

2 Derivation of the model

We consider a population that is structured in an arbitrarliy many $u$ groups interacting in $v$ patches. We consider a typical disease captured by an SEIRS structure. Naturally, $S_{i}$, $E_{i}$, $I_{i}$ and $R_{i}$ are the susceptible, latent, infectious and recovered individuals of Group $i$ respectively. The population of each group is denoted by $N_{i}=S_{i}+E_{i}+I_{i}+R_{i}$, for $i=1,\dots,u$. Individuals of Group $i$ spend on average some time in Patch $j$, $j=1,\dots,v$. The susceptible, latent, infected and recovered populations of group $i$ spend $m_{ij}$, $n_{ij}$, $p_{ij}$ and $q_{ij}$ proportion of times respectively in Patch $j$, for $j=1,\dots,v$. At time $t$, the effective population of Patch $j$ is $N_{j}^{\textrm{eff}}=\sum_{k=1}^{u}(m_{kj}S_{k}+n_{kj}E_{k}+p_{kj}I_{k}+q_{kj}R_{k})$. This effective population of Patch $j$ describes the temporal dynamics of the population in Patch $j$ weighted by the mobility patterns of each group and each epidemiological status. Of this patch population, $\sum_{k=1}^{u}p_{kj}I_{k}$ are infectious. The proportion of infectious individuals in Patch $j$ is therefore,

[TABLE]

Susceptible individuals of Group $i$ could be infected in any Patch $j$, $j=1,\dots,v$ while visiting there. Hence, the dynamics of susceptible of Group $i$ is given by:

[TABLE]

where $\Lambda_{i}$ denotes a constant recruitment of susceptible individuals of Group $i$, $\mu_{i}$ the natural death rate, $\beta_{j}$ the risk of infection and $\eta_{i}$ the immunity loss rate. The patch specific risk vector $\mathcal{B}=(\beta_{j})_{1\leq j\leq v}$ is treated as constant. However, in Subsection 5.2, we also considered the case when this risk depends on the effective population size.

The latent individuals of Group $i$ are generated through infection of susceptible and decreased by natural death and by becoming infectious at the rate $\nu_{i}$. Hence the dynamics of latent of Group $i$, for $i=1,\dots,u$, is given by:

[TABLE]

The dynamics of infectious individuals of Group $i$ is given by

[TABLE]

where $\gamma_{i}$ is the recovery rate of infectious individuals. Finally, the dynamics of recovered individuals of Group $i$ is:

[TABLE]

The complete dynamics of $u$-groups and $v$-patches SEIRS epidemic model is given by the following system:

[TABLE]

The description of parameters in Model (1) is given in Table 1. These parameters are composed of three set of parameters: ecological/environmental (number of patches $v$ and their risk $\mathcal{B}$), epidemiological (Recruitment, death rates, recovery rate, etc) and behavioral (mobility matrices) parameters. A schematic description of the flow is given in Fig 1.

Model (1) could be written in the compact form,

[TABLE]

where $\mathbf{S}=[S_{1},S_{2},\dots,S_{u}]^{T}$, $\mathbf{E}=[E_{1},E_{2},\dots,E_{u}]^{T}$, $\mathbf{I}=[I_{1},I_{2},\dots,I_{u}]^{T}$ and $\mathbf{R}=[R_{1},R_{2},\dots,R_{u}]^{T}$. The matrices $\mathbb{M}=(m_{ij})_{\begin{subarray}{c}1\leq i\leq u,\\ 1\leq j\leq v\end{subarray}}$, $\mathbb{N}=(n_{ij})_{\begin{subarray}{c}1\leq i\leq u,\\ 1\leq j\leq v\end{subarray}}$, $\mathbb{P}=(p_{ij})_{\begin{subarray}{c}1\leq i\leq u,\\ 1\leq j\leq v\end{subarray}}$ and $\mathbb{Q}=(q_{ij})_{\begin{subarray}{c}1\leq i\leq u,\\ 1\leq j\leq v\end{subarray}}$ represent the residence time matrices of susceptible, latent, infectious and recovered individuals respectively. Moreover, ${\mathbf{\Lambda}}=[\Lambda_{1},\Lambda_{2},\dots,\Lambda_{u}]^{T}$, $\mathcal{B}=[\beta_{1},\beta_{2},\dots,\beta_{v}]^{T}$, $\mu=[\mu_{1},\mu_{2},\dots,\mu_{u}]^{T}$, $\nu=[\nu_{1},\nu_{2},\dots,\nu_{u}]^{T}$, $\gamma=[\gamma_{1},\gamma_{2},\dots,\gamma_{u}]^{T}$, $\delta=[\delta_{1},\delta_{2},\dots,\delta_{u}]^{T}$ and $\eta=[\eta_{1},\eta_{2},\dots,\eta_{u}]^{T}$.

Model (2) brings added value to the existing literature in the following ways:

The structure of the host population is different and independent from the patches where the infection takes place. Indeed, in the previous epidemic models describing human dispersal or mixing (Eulerian or Lagrangian), hosts’ structure unit and the geographical landscape unit, be it group or patch, is the same and homogeneous in term of transmission rate. Our model captures added heterogeneity in the sense that we decouple the structure of the host to that of patches. For instance, our framework fits well for nosocomial diseases (hospital-acquired infections), where the hospitals could be treated as patches and host’s groups as gender or age (see [18, 30] for the effects of gender and age on nosocomial infections). 2. 2.

In our formulation, there is no need to measure contacts rates, a difficult task for nearly all diseases that are not either sexually transmitted or vector-borne. Each patch is defined by its specific risk of infection that could be tied to environmental or hygienic conditions. Hence, susceptibility is not individual-based nor group-based as in classical formulation of multi-group models (the contact matrices in these type of models are known as WAIFW, i.e., Who Acquires Infection From Whom [1]), but a patch specific risk. In fact, our framework is capable of capturing a wide-range of modeling scenarios, including group-susceptibility. Indeed, if $g_{i}$ is the risk of infection of Group $i$, $i=1,2,\dots,u$, it suffices to replace $S_{i}$ by $g_{i}S_{i}$ in only the infection terms in (1). That is, the dynamics of susceptible and latent hosts, for $i=1,2,\dots,u$ will be:

[TABLE]

and

[TABLE]

For the sake of simplicity, we considered the case where all host groups have the same risk of infection, though all the results obtained in this paper hold without this simplification.

The risk in each patch may be fixed, as in Model (2), or variable and dependent of the effective patch population (See Subsection 5.2). The prospect of infection is tied to the environmental risk and time spent in that environment. This fits, for example, pandemic influenza in schools and, again, the nosocomial infections (length of stay in hospitals and their corresponding risks). These residences times and patch related risks are easier to quantify than contact rates. This paper extend earlier results in [7, 9]. 3. 3.

The model allows individuals of different groups to move across patches without losing their identities. This approach allows a more targeted control strategy for public health benefit. Therefore, the model follows a Lagrangian approach and generalize [7, 9, 15, 26, 38, 39]. 4. 4.

There are different mobility patterns depending on the epidemiological class of individuals. This allows us to highlight and assess the effects of hosts’ behavior through social distancing and their predilection for specific patches on the disease dynamics. Although the differential mobility have been considered in an Eulerian setting [41, 46], its incorporation in a Lagrangian setting is new and is an extension of [7, 9, 15, 19, 26, 38, 39] (for which mobility is independent of hosts’ epidemiological class). 5. 5.

In this framework, we consider only patches where the infection takes place (hospitals, schools, malls, etc) whereas previous models suppose that the patches are distributed over the whole space. In short, the mobility matrices are not assumed to be stochastic.In this case, a natural condition on the mobility matrices arises: $\mathbb{X}\mathbf{1}\leq\mathbf{1}$, for $\mathbb{X}\in\{\mathbb{M},\mathbb{N},\mathbb{P},\mathbb{Q}\}$, where $\mathbf{1}$ is the vector whose components are all equal to unity. These conditions stem from the fact that the added proportion of time spend in all patches cannot be more that $100\%$. However, as pointed out by a reviewer, the stochasticity of the mobility matrices is not really restrictive. Indeed, as we are considering an arbitrary number of patches, we can, without loss of generality, add an additional patch within which individuals spent “the rest of their time” and where no infection takes place in it. That is, $\beta_{v+1}=0$.

We denote by $\mathbf{N}$ the vector of populations of each group. The dynamics of the population in each group is given by the following:

[TABLE]

where $\circ$ denotes the Hadamard product. Thus, the set defined by

[TABLE]

is a compact attracting positively invariant for System (2).

The disease-free equilibrium (DFE) of System (2) is given by $(\mathbf{S}^{\ast},\mathbf{0},\mathbf{0},\mathbf{0})$ where $\mathbf{S}^{\ast}=\mathbf{\Lambda}\circ\dfrac{1}{\mu}$.

Remark 2.1**.**

If the susceptible or infected individuals do not go to the patches where the infection takes place, either due to intervention strategy or social distancing, that is when the residence time matrices $\mathbb{M}$ or $\mathbb{P}$ are the null matrix (the susceptible individuals do not spend any time in the considered patches), the disease does not spread and eventually dies out.

We compute the basic reproduction number following [16, 44]. By decomposing the infected compartments of (2) as a sum of new infection terms and transition terms,

[TABLE]

The Jacobian matrix at the DFE of $\mathcal{F}(\mathbf{E},\mathbf{I})$ and $\mathcal{V}(\mathbf{E},\mathbf{I})$ are given by:

[TABLE]

and,

[TABLE]

Hence, we obtain

[TABLE]

The basic reproduction number is the spectral radius of the next generation matrix

[TABLE]

where

[TABLE]

Finally, the basic reproduction number for $u$ groups and $v$ patches is given by

[TABLE]

The disease-free equilibrium is asymptotically stable whenever $\mathcal{R}_{0}(u,v)<1$ and unstable if $\mathcal{R}_{0}(u,v)>1$ [16, 44].

3 Effects of heterogeneity on the basic reproduction number

In this section, we investigate the effects of patchiness, groupness and mobility on the basic reproduction number. More particularly, how the basic reproduction number changes its monotonicity with respect to the number of patches, groups and mobility patterns of individuals.

The following theorem gives the monotonicity of the basic reproduction with respect the residence times patterns of the infected individuals.

Theorem 3.1**.**

The basic reproduction number $\mathcal{R}_{0}(u,v)$ is a nondecreasing function with respect to $\mathbb{P}$, that is, the infected individuals movement patterns.

Proof.

Recall that $\mathcal{R}_{0}(u,v)=\rho(Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma+\delta)))$ where,

$Z=\textrm{diag}(\mathbf{S^{\ast}})\mathbb{M}\textrm{diag}(\mathcal{B})\textrm{diag}^{-1}(\mathbb{M}^{T}\mathbf{S^{\ast}})\mathbb{P}^{T}$. The matrix $Z$ is linear in $\mathbb{P}$ and has all non-negative entries. We consider the order relation for the matrices as follows: $A\leq B$ if $a_{ij}\leq b_{ij}$, for all $i$ and all $j$, where $a_{ij}$ and $b_{ij}$ are entries of $A$ and $B$ respectively. Also, $A<B$ if $A\leq B$ and $A\neq B$. Hence, since the Perron-Frobenius theorem [5] (Corollary 1.5, page 27) guarantees that for any positives matrices $A$ and $B$ such that $A\geq B\geq 0$, then $\rho(A)\geq\rho(B)$, we deduce that, for any matrix $\mathbb{P}^{\prime}\geq\mathbb{P}$,

[TABLE]

∎

The variation in monotonicity of $\mathcal{R}_{0}(u,v)$ with respect to the residence times patterns of susceptible individuals, that is $\mathbb{M}$, is more complicated and difficult to assess in general and even in some more restrictive particular cases (see Remark 3.2).

Hereafter, we define two bounding quantities tied to the global basic reproduction number:

[TABLE]

and,

[TABLE]

It is worthwhile noting that $\mathcal{R}_{0}^{i}=\mathcal{R}_{0}(1,v)$. That is, $\mathcal{R}_{0}^{i}$ is also the basic reproduction number of the global system in presence of one group only, namely the $i^{th}$, spread over $v$ patches. $\mathcal{R}_{0}^{i}$ could be seen as a group specific “reproduction number”.

The quantity $\tilde{\mathcal{R}}_{0}^{i}(u,v)$ could be heuristically seen as the sum of the average number of cases produced by an infected of group $i$ over all patches, in presence of other groups.

In the following theorem, we explore how the general basic reproduction number $\mathcal{R}_{0}(u,v)$ is tied to these specific reproduction numbers and whether it increases or decreases when the number of patches and/or groups changes. An underlying assumption in the following theorem is that when adding patches, the proportion of time spent in the existing patches remain exactly the same.

Theorem 3.2**.**

We have the following inequalities:

$\displaystyle\max\left\{\max_{i=1,\dots,u}\tilde{\mathcal{R}}_{0}^{i}(u,v),\min_{i=1,\dots,u}\mathcal{R}_{0}^{i}\right\}\leq\mathcal{R}_{0}(u,v)\leq\max_{i=1,\dots,u}\mathcal{R}_{0}^{i}$** 2. 2.

$\mathcal{R}_{0}(u,v)\geq\mathcal{R}_{0}(1,v)\geq\mathcal{R}_{0}(1,1).$** 3. 3.

For a fixed number of groups $u$, $\mathcal{R}_{0}(u,v)\geq\mathcal{R}_{0}(u,v^{\prime})$ where $v$ and $v^{\prime}$ are integers such that $v\geq v^{\prime}$.

Proof.

1. We prove first that $\mathcal{R}_{0}(u,v)\geq\displaystyle\max_{i=1,\dots,u}\tilde{\mathcal{R}}_{0}^{i}(u,v)$ and then $\displaystyle\min_{i=1,\dots,n}\mathcal{R}_{0}^{i}\leq\mathcal{R}_{0}(u,v)\leq\max_{i=1,\dots,n}\mathcal{R}_{0}^{i}$.

Let $e_{i}$ the $i-$th vector of the canonical basis of $\mathbb{R}^{4u}$. We have

[TABLE]

It follows that,

[TABLE]

We also have

[TABLE]

Since $\mathbb{P}^{T}e_{i}$ is the $i-$th column of $\mathbb{P}^{T}$, we obtain:

[TABLE]

Hence, the diagonal elements of $\mathbb{M}\textrm{diag}(\mathcal{B})\textrm{diag}(\mathbb{M}^{T}\mathbf{S^{\ast}})^{-1}\mathbb{P}^{T}$ is given by

[TABLE]

This implies that, for all $i=1,\cdots,v$, $\tilde{\mathcal{R}}_{0}^{i}(u,v)$ is a diagonal element of the next generation matrix. Since the spectral radius of a matrix is the greater or equal to its diagonal elements, we can conclude that $\mathcal{R}_{0}(u,v)\geq\tilde{\mathcal{R}}_{0}^{i}$ for all $i=1,\cdots,u$. This implies that

[TABLE]

It remains to prove that $\displaystyle\min_{i=1,\dots,u}\mathcal{R}_{0}^{i}\leq\mathcal{R}_{0}(u,v)\leq\max_{i=1,\dots,u}\mathcal{R}_{0}^{i}$. The basic reproduction number is given by $\mathcal{R}_{0}(u,v)=\rho(Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma+\delta)))$ where

[TABLE]

It can be shown that the elements of this matrix are the following:

[TABLE]

If $\mathbb{M}\mathbb{P}^{T}$ is irreducible, the matrix $Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma+\delta)))$ is irreducible, and therefore its spectral radius satisfy the Frobenius’ inequality ([24], Theorem 8.1.22, page 492):

[TABLE]

where $z_{j}=\sum_{i=1}^{u}z_{ij}$ and $z_{ij}$ are given by (13). We have:

[TABLE]

Hence,

[TABLE]

The relations (12) and (14) imply the desired inequality.

2. By using the inequality proved in the first part, we have:

[TABLE]

Finally, we have:

[TABLE]

3. Let $u$ a fixed number of groups. We would like to prove that $\mathcal{R}_{0}(u,v)\geq\mathcal{R}_{0}(u,v^{\prime})$ for any $v\geq v^{\prime}$. Since, $\mathcal{R}_{0}(u,v)=\rho(Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma+\delta)))$ and the number of groups is fixed, the epidemiological parameters remain the same for any number of patches. Hence, it remains to compare $Z_{v}$ and $Z_{v^{\prime}}$ where $Z$ is the part of the next generation matrix that depends on the number of patches.

For $v$ patches, we have

[TABLE]

For $v^{\prime}$ patches,

[TABLE]

Hence, for $v\geq v^{\prime}$, we have clearly $Z_{v}^{ij}\geq Z_{v^{\prime}}^{ij}$. Hence, thanks to Perron-Frobebenius’ theorem, we conclude that $\mathcal{R}_{0}(u,v)\geq\mathcal{R}_{0}(u,v^{\prime})$. ∎

Remark 3.1**.**

• •

The inequality in Item 3 of Theorem 3.2 is independent of the risk of infection in the additional patches.

• •

If the residence times network configuration changes due the newly added patches, the increasing property of the basic reproduction number with respect to the number of patches (Item 3 of Theorem 3.2) may not hold. This is an interesting avenue to exploring the monotonicity of $\mathcal{R}_{0}$ and/or the dynamics of the disease.

We investigate relevant modeling scenarios where the expression of the general basic reproduction number for $u$ patches and $v$ patches, $\mathcal{R}_{0}(u,v)$, could be explicitly obtained. In the rest of the paper, we use \langle x\,\mbox{\Large\mid}\,y\rangle to denote the canonical scalar product.

Theorem 3.3**.**

If the susceptible residence times matrix $\mathbb{M}$ is of rank one, an explicit expression of $\mathcal{R}_{0}$ is given by

[TABLE]

where $\xi\in{{\rm I\mkern-3.5muR}}^{u}$ is such that $\mathbb{M}=\xi^{T}m$, with $m\in{{\rm I\mkern-3.5muR}}^{v}$. Moreover, if the matrix $\mathbb{M}$ is stochastic, we have:

[TABLE]

Proof.

If the susceptible residence times matrix $\mathbb{M}$ is of rank one, it exist a vector $\xi\in{{\rm I\mkern-3.5muR}}^{u}$ and a vector $m\in{{\rm I\mkern-3.5muR}}^{v}$ such that $\mathbb{M}=\xi m^{T}$. We have the following:

[TABLE]

Hence,

[TABLE]

and

[TABLE]

We deduce that the non-zero diagonal block of the next generation matrix could be written as:

[TABLE]

This matrix is clearly of rank 1, since it could be written as $wz^{T}$ where $w\in{{\rm I\mkern-3.5muR}}^{u}$ and $w\in{{\rm I\mkern-3.5muR}}^{v}$. Hence, its unique non zero eigenvalue is

[TABLE]

or, equivalently,

[TABLE]

Now, if $\mathbb{M}$ is of rank one and stochastic, that is , $\sum_{j=1}^{v}m_{ij}=1$, for all $i=1,\dots,u$, it is not difficult to show that $\xi=\mathbf{1}$, where $\mathbf{1}$ is the vector whose components are all equal to unity. This leads to

[TABLE]

∎

Remark 3.2**.**

If the residence times matrix of susceptible individuals, that is $\mathbb{M}$, is of rank one and stochastic, the basic reproduction number is independent of $\mathbb{M}$.

It is worthwhile noting that there is a special case for which the result of Remark 3.2 holds even if the matrix $\mathbb{M}$ is not stochastic but only of rank one and sub-stochastic. Indeed, by adding a new patch $v+1$ with $\beta_{v+1}=0$ where the hosts of different groups spend “the rest of their times”, the new mobility matrices will become the matrices $\tilde{\mathbb{M}}=(\mathbb{M},\mathbb{M}^{\prime}),\tilde{\mathbb{N}}=(\mathbb{N},\mathbb{N}^{\prime}),\tilde{\mathbb{P}}=(\mathbb{P},\mathbb{P}^{\prime})$ and $\tilde{\mathbb{Q}}=(\mathbb{Q},\mathbb{Q}^{\prime})$, where $\mathbb{M}^{\prime},\mathbb{N}^{\prime},\mathbb{P}^{\prime}$ and $\mathbb{Q}^{\prime}$ are column vectors. The new mobility matrices are stochastic and $\mathcal{R}_{0}(u,v,\mathbb{M},\mathbb{P})=\mathcal{R}_{0}(u,v+1,\tilde{\mathbb{M}},\tilde{\mathbb{P}})$ since $\beta_{v+1}=0$. Hence, if $\mathbb{M}$ and $\tilde{\mathbb{M}}$ are of rank one, the basic reproduction number is still independent of $\mathbb{M}$. In this case, the matrix $\mathbb{M}$ could be expressed as $\mathbf{1}m^{T}$ with $\sum_{j=1}^{v}m_{j}<1$. Thus, there is a special case when $\mathbb{M}$ is rank 1, yet sub-stochastic, and the reproduction number does not depend on $\mathbb{M}$.

From a modeling standpoint, the rank one condition of $\mathbb{M}$ (i.e., $\mathbb{M}=\xi m^{T}$ with $\xi\in{{\rm I\mkern-3.5muR}}^{u}$ and $m\in{{\rm I\mkern-3.5muR}}^{v}$) can be interpreted as follows:

$\bullet$ The ratio of the proportions of time spent in any given patch by susceptible individuals belonging to two different groups, is identical. Indeed, for any given group $i$, the ratio of the proportion of time spent in any given patch by susceptible individual is:

[TABLE]

which is independent of $i$. Moreover, if $\mathbb{M}$ is stochastic, we deduce that the susceptible of each group spend the exact proportion of time in any given patch, since $\sum_{k=1}^{v}m_{k}=1$.

$\bullet$ A straightforward case that stems from the previous point is whenever there is one patch and multiple groups; or when there are multiple patches and one group.

Similar remarks hold when the matrix $\mathbb{P}$ is of rank one, which is dealt in the next theorem.

Theorem 3.4**.**

If the infected residence times matrix $\mathbb{P}$ is of rank one, an explicit expression of $\mathcal{R}_{0}$ is given by

[TABLE]

where $\alpha\in{{\rm I\mkern-3.5muR}}^{u}$ and $p\in{{\rm I\mkern-3.5muR}}^{v}$ are such that $\mathbb{P}=\alpha p^{T}$. Moreover, if $\mathbb{P}$ is stochastic,

[TABLE]

Proof.

If the susceptible residence times matrix $\mathbb{P}$ is of rank one, there exists a vector $p\in{{\rm I\mkern-3.5muR}}^{v}$ and $\alpha\in{{\rm I\mkern-3.5muR}}^{u}$ such that $\mathbb{P}=\alpha p^{T}$. The next generation matrix is:

[TABLE]

which is of rank one since it could be written as $xy^{T}$ where $x=\textrm{diag}(\mathbf{S^{\ast}})\mathbb{M}\textrm{diag}(\mathcal{B})\textrm{diag}^{-1}(\mathbb{M}^{T}\mathbf{S^{\ast}})p$ and $y=\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma+\delta))\alpha$. Hence, its unique non zero eigenvalue is,

[TABLE]

If $\mathbb{P}$ is stochastic, we can show that $\alpha=\mathbf{1}$ and hence,

[TABLE]

which is the desired result. ∎

The condition of rank one of the matrices $\mathbb{M}$ and $\mathbb{P}$, when both matrices are stochastic, means that the susceptible and infected individuals of different groups spend the same proportion of time in each and every patch. When the matrices are not stochastic, the rank one condition means that the proportion of times spent by susceptible or infected individuals of different groups in each patch are proportional. That is, there exists $\alpha_{j}$ such that $m_{ij}=\alpha_{j}m_{i}$ for all $1\leq i,j\leq u$.

4 Simulations

In this section, we run some numerical simulations for 2 groups and 3 patches in order to highlight the effects of differential residence times and to illustrate the previously obtained theoretical results. To that end, unless otherwise stated, the baseline parameters of the model are chosen as follows:

[TABLE]

[TABLE]

[TABLE]

Although the values of $\beta_{j}$ are chosen throughout this section, for convenience, to be between 0 and 1, they need only to be nonnegative. We begin by simulating the dynamics of Model 2 when the basic reproduction number is below or above unity. Figure 2 shows the dynamics of infected individuals of Group 1 (Fig. 2(a)) and Group 2 (Fig. 2(b)). The disease persists in both groups when $\mathcal{R}_{0}>1$ ( Fig. 2(a) and Fig. 2(b), dotted red and dashed green curves) while it dies out when $\mathcal{R}_{0}<1$ ( Fig. 2(a) and Fig 2(b), solid blue and dash-dotted black curves).

Figure (3) displays how the equilibrium value of infected individuals of Group 1 and Group 2 change with respect to residence times of infected and susceptible of group 1 in Patch 1, that is $m_{11}$ and $p_{11}$. For instance, in Fig. 3(a), the disease burden in Group 1 ($\bar{I}_{1}$) is moderately low for all values of $m_{11}$ as long as $p_{11}$, the residence times of Group 1’s infected into Patch 1, is below 0.3, even if Patch 1 is the riskiest patch with $\beta_{1}=0.35$. However, this prevalence level is more marked when $m_{11}\geq 0.4$ and $p_{11}\geq 0.5$. The heatmap of $\bar{I_{1}}$ with respect to $m_{12}$ and $p_{21}$ shows similar patterns. We decided not to display this figure. Fig 3(b) shows the changes in the values of infected in Group 2 ($\bar{I}_{2}$) due to movement patterns of susceptible and infected of Group 1 ($m_{11}$ and $p_{11}$) when their own movement patterns are fixed ($m_{21}=0.6$ and $p_{21}=0.4$).

Fig 4 gives an overview of the dynamics of the basic reproduction number with respect of mobility patterns of susceptible and infected individuals of Group 1 in Patch 1 and Patch 2. Fig 4(a) shows that $m_{11}$ and $p_{11}$ could bring $\mathcal{R}_{0}$ from bellow unity to above unity. Particularly, if $m_{11}\geq 0.4$, then $\mathcal{R}_{0}>1$, which lead to the persistence of the disease. Also, $\mathcal{R}_{0}$ is much higher when $m_{11}\geq 0.7$ and $p_{11}\geq 0.2$. Fig 4(b) shows how $\mathcal{R}_{0}$ varies when the movement of infected and susceptible of Group 1 in Patch 2 change.

In Figure (5), we revisit the variability of the basic reproduction number with respect of mobility patterns of susceptible and infected individuals of Group 1 (Fig 4). However, we obtain a clear picture on how it changes. Indeed, Fig. 5(a)) suggests that $\mathcal{R}_{0}$ increases with respect to $m_{11}$ and $m_{12}$; and $p_{11}$ and $p_{12}$ (Fig. 5(b)). However, $\mathcal{R}_{0}$ increases much faster with respect to $p_{11}$ than to $m_{11}$. Moreover, Fig 5(b) confirms also the result of Theorem 3.1, which states that the basic reproduction number increases with respect of $p_{ij}$, that is the movement patterns of infected individuals.

Fig (6) showcases that, for a fixed number of groups (3 in this case), the basic reproduction number increases as the number of patches increases, and that independently of the values of the risk of infection of the added patches. This figure, also confirms our the theoretical result in Item 3 of Theorem 12. It also shows a linear monotonicity of $\mathcal{R}_{0}(u,v)$ with respect to $\mathbb{P}$. Other values of $\beta$s than those of Fig. 6 exhibit similar patterns.

5 Global stability of equilibria

The global stability of equilibria for the general Model (2) happens to be very challenging. In fact, for models with such intricated nonlinearities, it is shown in [25] that multiple endemic equilibria may exist. In this section, we explore the global stability of equilibria for some particular cases of the general model.

5.1 Identical Mobility and no disease induced mortality

In this subsection, we suppose that the host mobility to the patches is independent of the epidemiological status and that we neglect the disease induced mortality. In this case, the dynamics of the total population is given by

[TABLE]

Hence, $\displaystyle\lim_{t\to\infty}\mathbf{N}=\frac{\mathbf{\Lambda}}{\mu}:=\bar{\mathbf{N}}$. By using the theory of asymptotic systems [14, 45], System (2) is asymptotically equivalent to:

[TABLE]

Model (16) generalizes models considered in [9]. Let us denote $\mathcal{R}_{0}^{\textrm{Eq}}(u,v)$ the corresponding basic reproduction number of Model (16). Its expression is

[TABLE]

The following theorem gives the global stability of the disease free equilibrium.

Theorem 5.1**.**

Whenever the host-patch mobility configuration $\mathbb{M}\mathbb{M}^{T}$ is irreducible, the following statements hold:

If $\mathcal{R}_{0}^{\textrm{Eq}}(u,v)\leq 1$, the DFE is globally asymptotically stable (GAS). 2. 2.

If $\mathcal{R}_{0}^{\textrm{Eq}}(u,v)>1$, the DFE is unstable.

Proof.

Let $(w_{E},w_{I})$ a left eigenvector of $Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma))$ corresponding to $\rho(Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma)))$ where

[TABLE]

Hence,

[TABLE]

Since $\mathbb{M}\mathbb{M}^{T}$ is irreducible, the matrix $Z\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma))$ is irreducible. This implies that $(w_{E},w_{I})\gg 0$.

We consider the Lyapunov function

[TABLE]

The derivative of $V(\mathbf{E},\mathbf{I})$ along trajectories of (16) is

[TABLE]

where $\tilde{w}_{E}=w_{E}\textrm{diag}^{-1}(\mu+\nu)+w_{I}\textrm{diag}(\nu)\textrm{diag}^{-1}((\mu+\nu)\circ(\mu+\gamma))$ and $\tilde{w}_{I}=w_{I}\textrm{diag}^{-1}(\mu+\gamma)$, or equivalently $(\tilde{w}_{E},\tilde{w}_{I})=(w_{E},w_{I})(-V^{-1})$.

Since $\textrm{diag}(\mathbf{S})\leq\textrm{diag}(\mathbf{S^{\ast}})$ and $\mathbf{S^{\ast}}=\bar{\mathbf{N}}$, we obtain (denoting $\mathbb{I}$ the identity matrix),

[TABLE]

We consider first the case when $\mathcal{R}_{0}^{\textrm{Eq}}(u,v)<1$. Let $\mathcal{E}$ be an invariant set contained in $\Omega$, where $\dot{V}(\mathbf{E},\mathbf{I})=0$. This set is reduced to the origin of ${{\rm I\mkern-3.5muR}}^{2u}$ (i.e., $(\mathbf{E},\mathbf{I})=(0,0)$). This, combined to the invariance of $\mathcal{E}$, leads to $\mathbf{R}=0$ and $\mathbf{S}=\mathbf{S}^{\ast}$. Hence, the only invariant set contained in $\Omega$, such that $\dot{V}(\mathbf{E},\mathbf{I})=0$, is reduced to the DFE. Hence, by LaSalle’s invariance principle [6, 33], the DFE is globally asymptotically stable on $\Omega$. Since $\Omega$ is an attracting set, we conclude that the DFE is GAS on the positive orthant ${{\rm I\mkern-3.5muR}}^{4u}_{+}$.

When $\mathcal{R}_{0}^{\textrm{Eq}}(u,v)=1$, we can show that

[TABLE]

Therefore, as above, LaSalle’s invariance principle allows to conclude.

The instability of the DFE when $\mathcal{R}_{0}^{\textrm{Eq}}(u,v)>1$ follows from [16, 44]. ∎

The following theorem provides the uniqueness of the endemic equilibrium.

Theorem 5.2**.**

If $\mathcal{R}_{0}^{\textrm{Eq}}(u,v)>1$, Model (16) has a unique endemic equilibrium.

The proof of this theorem is similar to that of Theorem 5.3 in the next subsection.

5.2 Effective population size dependent risk

So far, the risk associated with each patch is represented by the constant vector $\mathcal{B}$. However, in some cases, it is more appropriate to assume that the risk of catching a disease depends on the size of the population or crowd, that is the effective population size in each patch. In this subsection, we suppose that the risk of infection in each patch $j$ is linearly proportional to the effective population size, that is $N_{j}^{\textrm{eff}}=\sum_{k=1}^{u}(m_{ij}S_{i}+n_{ij}E_{i}+p_{ij}I_{i}+q_{ij}R_{i})$. Hence,

[TABLE]

Hence, the rate at which susceptible individuals are infected in Patch $j$ is, therefore

[TABLE]

Therefore, in this settings, the dynamics of the population in different epidemiological classes take the form:

[TABLE]

System (36) could be written in a compact form as follows:

[TABLE]

Clearly, System (37) is a particular case of System (2) when the transmission term takes a modified density-dependent form. Positivity and boundedness properties of solutions of System (2) hold for those of System (37). The basic reproduction number of Model (37), denoted by $\mathcal{R}_{0}^{\textrm{DD}}(u,v)$ is:

[TABLE]

We explore the properties of steady state solutions. The following result gives the global stability of the DFE. Its proof is similar to that of Theorem 5.1.

Corollary 5.1**.**

Whenever the host-patch mobility configuration $\mathbb{M}\mathbb{P}^{T}$ is irreducible, the following statements hold:

If $\mathcal{R}_{0}^{\textrm{DD}}(u,v)\leq 1$, the DFE is globally asymptotically stable. 2. 2.

If $\mathcal{R}_{0}^{\textrm{DD}}(u,v)>1$, the DFE is unstable.

The proof of the existence and uniqueness of the endemic equilibrium (EE) for Model (37) is done in two steps, by carefully crafting a new auxiliary system whose EE uniqueness is tied to that of Model (37).

Let

[TABLE]

We have the following lemma,

Lemma 5.1**.**

Model (37) has a unique endemic equilibrium if the function

[TABLE]

has a unique fixed point.

Proof.

Let ($\mathbf{\bar{S}},\mathbf{\bar{E}},\mathbf{\bar{I}},\mathbf{\bar{R}}$) an equilibrium point of System (37) with $\mathbf{\bar{I}}\gg 0$. This equilibrium satisfies the following system:

[TABLE]

We can easily see that $\mathbf{\bar{R}}=A\mathbf{\bar{E}}$ and $\mathbf{\bar{I}}=L\mathbf{\bar{E}}$, where $A$, $L$ and $K$ are as defined in (38). Thus, $\mathbf{\bar{I}}\gg 0$ implies that $\mathbf{\bar{E}}\gg 0$ and $\mathbf{\bar{R}}\gg 0$.

Hence, System (39) could be written only in terms of $\mathbf{\bar{S}}$ and $\mathbf{\bar{E}}$, that is:

[TABLE]

Let $x=\textrm{diag}^{-1}(\mu)\Lambda-\mathbf{\bar{S}}$ and $y=\mathbf{\bar{E}}$. Since $\mathbf{\bar{S}}\in\Omega$, it is clear that $x\geq\mathbf{0}$ and $y\geq\mathbf{0}$. Expressing the system (40) into new variables, we obtain:

[TABLE]

[TABLE]

where

[TABLE]

It follows from (41b) that $f(x,y)=\textrm{diag}(\nu+\mu)\,y$, and hence (41a) implies that $x=Ky$ where

[TABLE]

After some algebraic manipulations, it could be shown that $K>0$. Combining the fact that $x=Ky$ and (41b), it follows that (41), and subsequently (39), could be written in the single vectorial equation:

[TABLE]

where

[TABLE]

Thus, Model (37) has a unique endemic equilibrium $\mathbf{\bar{I}}\gg 0$ if and only if $g(y)$ has a unique fixed point $\bar{y}\gg 0$. The desired result is achieved. ∎

Next, we present another lemma whose proof is straightforward:

Lemma 5.2**.**

The function $g(y)$ has a fixed point $\bar{y}$ if and only if $\bar{y}$ is an equilibrium of $\dot{y}=F(y)$ where

[TABLE]

The proof of this lemma is straightforward.

Theorem 5.3**.**

Under the assumption that the host-patch mobility configuration $\mathbb{M}\mathbb{P}^{T}$ is irreducible, Model (37) has a unique endemic equilibrium whenever $\mathcal{R}_{0}^{\textrm{DD}}(u,v)>1$.

Proof.

By using Lemma 5.1 and Lemma 5.2, the uniqueness of EE for Model (37) is equivalent to the uniqueness of an EE for this system

[TABLE]

when $\mathcal{R}_{0}^{\textrm{DD}}(u,v)>1$. Therefore, we will prove that the auxiliary system (43) has an unique EE. In fact, we will prove that this equilibrium is globally attractive if $\mathcal{R}_{0}^{\textrm{DD}}(u,v)>1$. The proof of the latter is based on Hirsch’s theorem [23], by using elements of monotone systems. The Jacobian of the vector field $F(y)$ is:

[TABLE]

where $\mathbb{I}$ is the identity matrix. The matrix -\textrm{diag}(\nu+\mu)\,\mathbb{I}-\textrm{diag}\Big{(}\mathbb{M}\textrm{diag}(\mathcal{B})\,\mathbb{P}^{T}\,Ly\Big{)}\,K is a diagonal matrix and $\textrm{diag}(S^{*}-Ky)\mathbb{M}\textrm{diag}(\mathcal{B})\,\mathbb{P}^{T}\,L$ is a nonnegative matrix (since $S^{*}-Ky=\mathbf{\bar{S}}$). It follows that $F^{\prime}(y)$ is Metzler and is irreducible since $\mathbb{M}\,\mathbb{P}^{T}$ is. Therefore, System (43) is strongly monotone. Moreover, it is clear that the map $F^{\prime}:{{\rm I\mkern-3.5muR}}^{u}\longrightarrow{{\rm I\mkern-3.5muR}}^{u}\times{{\rm I\mkern-3.5muR}}^{u}$ is antimonotone. Also, $F(0_{{{\rm I\mkern-3.5muR}}^{u}})=0_{{{\rm I\mkern-3.5muR}}^{u}}$ and $F^{\prime}(0_{{{\rm I\mkern-3.5muR}}^{u}})=g^{\prime}(0_{{{\rm I\mkern-3.5muR}}^{u}})-\mathbb{I}=\textrm{diag}(S^{*})\mathbb{M}\textrm{diag}(\mathcal{B})\,\mathbb{P}^{T}\,L-\mathbb{I}$. Since $\rho(g^{\prime}(0_{{{\rm I\mkern-3.5muR}}^{u}}))=\rho(\textrm{diag}(S^{*})\mathbb{M}\textrm{diag}(\mathcal{B})\,\mathbb{P}^{T}\,L)=\mathcal{R}_{0}^{\textrm{DD}}(u,v)>1$, we deduce that $F^{\prime}(0_{{{\rm I\mkern-3.5muR}}^{u}})$ has at least a positive eigenvalue and therefore $0_{{{\rm I\mkern-3.5muR}}^{u}}$ is unstable. Therefore, System (43) has unique equilibrium $\bar{y}\gg 0_{{{\rm I\mkern-3.5muR}}^{u}}$, which is globally attractive, due to Hirsch’s theorem [23] (Theorem 6.1). We conclude that Model (37) has a unique endemic equilibrium whenever $\mathcal{R}_{0}^{\textrm{DD}}(u,v)>1$. ∎

Note that with the choice of $\mathbb{P}=\mathbb{M}\textrm{diag}^{-1}(\mathbb{M}^{T}\bar{\mathbf{N}})$ and $\delta=0$, System (37) is exactly System (16). Therefore, their solutions have the same asymptotic behavior.

6 Conclusion and discussions

Heterogeneity in space and social groups are often studied separately and sometimes interchangeably in the context of disease dynamics. Moreover, in these settings, susceptibility of the infection is based on group or individual. In this paper, we propose a new framework that incorporates heterogeneity in space and in group for which the structure of the latter is independent from that of the former. We define patch as a location where the infection takes place, which has a particular risk of infection. This risk is tied to environmental or hygienic or economic conditions that favors the infection. The likelihood of infection in each patch depends on both the risk of the patch and the proportion of time each host spend in that environment. We argue that this patch-specific risk is easier to assess compared to the classical differential susceptibility or WAIFW matrices. Human host is structured in groups, where a group is defined as a collection of individuals with similar demographic, genetic or social characteristics. In this framework, the population of each patch at time $t$ is captured by the temporal mobility patterns of all host groups visiting the patches, which in turn depends on the host’s epidemiological status.

Under this framework, we propose a general SEIRS multi-patch and multi-group model with differential state-host mobility patterns. We compute the basic reproduction number of the system with $u$ groups and $v$ patches, $\mathcal{R}_{0}(u,v)$, which depends on the mobility matrices of susceptible, $\mathbb{M}$, and infected, $\mathbb{P}$. The disease persists when $\mathcal{R}_{0}>1$ and dies out from all patches when $\mathcal{R}_{0}(u,v)<1$ (Fig. 2), when $\mathbb{M}\mathbb{P}^{T}$ is irreducible. When this matrix is not irreducible, the disease will persist or die out in all patches of the subsystem for which the configuration group-patch is irreducible and will be decoupled from the remaining system.

We systematically investigate the effects of heterogeneity in mobility patterns, groups and patches on the basic reproduction and on disease prevalence. Indeed, we have shown that, if the epidemiological parameters are fixed, the basic reproduction number is an increasing function of the entries of infected hosts’ movement matrix (e.g. Theorem 3.1). Also, if the number of groups is fixed, an increase in the number of patches increases the basic reproduction number (e.g. see Theorem 3.2). Explicit expressions of the basic reproduction numbers are obtained when the mobility matrices $\mathbb{M}$ and $\mathbb{P}$ are of rank one. That is, when, for all groups, susceptible (and infected) individuals’ residence times in all patches are proportional (Theorems 3.3 and 3.4). It turns out that if the susceptible residence time matrix is of rank one and stochastic, the basic reproduction number is independent of $\mathbb{M}$. Moreover, we also show that if $\mathbb{M}$ is of rank one, its stochasticity is sufficient but not necessary for the basic reproduction number to be independent of $\mathbb{M}$. However, if the infected residence time matrix $\mathbb{P}$ is of rank one, stochastic or otherwise, the basic reproduction number still depends on the infected movement patterns.

The patch-specific risk vector $\mathcal{B}$ could also depend on the effective population size. We explored the case when this dependence is linear, that is when, for each patch $j$, $\beta_{j}(N_{j}^{\textrm{eff}})=\beta_{j}\sum_{k=1}^{u}(m_{kj}S_{k}+n_{kj}E_{k}+p_{kj}I_{k}+q_{kj}R_{k})$. In this case, the transmission term of our model is captured by a density dependent incidence. Moreover, we show that this case is isomorphic to the general model, where the mobility patterns of host does not dependent on the epidemiological class, that is when $\mathbb{M}=\mathbb{N}=\mathbb{P}=\mathbb{Q}$. We prove that, in this case the disease free equilibrium is globally asymptotically stable whenever $\mathcal{R}_{0}\leq 1$ while an unique endemic equilibrium exists if $\mathcal{R}_{0}>1$.

We suspect that the disease free equilibrium is globally asymptotically stable whenever $\mathcal{R}_{0}\leq 1$ for Model (2), where the patch-specific risk is constant. A similar remark holds for the global stability of the endemic equilibrium of Model (37) and Model (2) when $\mathcal{R}_{0}>1$. This is still under investigation. Further areas of extensions of this study include more general forms of the patch-specific risks and when mobility patterns reflect the choices that individuals make at each point in time. These choices are based on maximizing the discounted value of an economic criterion à la [21, 36].

Acknowledgements

We are grateful to two anonymous referees and the handling editor Dr. Gabriela Gomes for valuable comments and suggestions that led to an improvement of this paper. We also thank Bridget K. Druken for the careful reading and constructive comments. A. Iggidr acknowledges the partial support of Inria in the framework of the program STIC AmSud (project MOSTICAW).