# A framework for the modelling and the analysis of epidemiological spread in commuting populations

Pierre-Alexandre Bliman *¹, Boureima Sangaré†², and Assane Savadogo‡¹,²

¹Sorbonne Université, Inria, CNRS, Université Paris Cité, Laboratoire Jacques-Louis Lions UMR 7598, Equipe MUSCLEES, Paris, France

²Laboratoire de Mathématiques, Informatique et Applications, Université Nazi Boni, Bobo-Dioulasso, Burkina Faso

August 28, 2024

## Contents

<table><tr><td>1 Introduction</td><td>2</td></tr><tr><td>2 Representation of commuting populations</td><td>5</td></tr><tr><td>2.1 Partitioned-population models</td><td>6</td></tr><tr><td>2.2 Admissible switching signals</td><td>8</td></tr><tr><td>2.3 Models of commuting populations</td><td>9</td></tr><tr><td>3 Epidemiological dynamics</td><td>9</td></tr><tr><td>3.1 A class of compartmental disease transmission models</td><td>10</td></tr><tr><td>3.2 Treatment model</td><td>12</td></tr><tr><td>3.3 Multigroup model</td><td>13</td></tr><tr><td>3.4 Staged progression model</td><td>14</td></tr><tr><td>3.5 Multistrain model</td><td>14</td></tr><tr><td>3.6 Vector-host model</td><td>15</td></tr><tr><td>4 Models of epidemics in commuting populations</td><td>15</td></tr><tr><td>4.1 A class of epidemiological models in commuting populations</td><td>15</td></tr><tr><td>4.2 Disease-free states and disease-free equilibrium points</td><td>16</td></tr><tr><td>4.3 Assumptions on the models</td><td>16</td></tr><tr><td>4.4 Basic reproduction number for a single-population model</td><td>18</td></tr><tr><td>5 Basic reproduction number of epidemiological models in commuting populations</td><td>22</td></tr><tr><td>5.1 Basic reproduction number of stationary partitioned-population models</td><td>22</td></tr><tr><td>5.1.1 Ad hoc representation of the state variable for partitioned-population model</td><td>22</td></tr><tr><td>5.1.2 Computation of the Jacobian matrix</td><td>23</td></tr><tr><td>5.1.3 DFE stability for a partitioned-population model</td><td>25</td></tr><tr><td>5.2 Basic reproduction number of general epidemiological models in commuting populations</td><td>27</td></tr><tr><td>5.2.1 Equilibrium states of a periodic commuting population model</td><td>27</td></tr><tr><td>5.2.2 State variable and Jacobian matrix representation for commuting population model</td><td>27</td></tr><tr><td>5.2.3 DFE stability for a commuting population model</td><td>28</td></tr></table>

<footnote>*Corresponding author, email: pierre-alexandre.bliman@inria.fr</footnote>

<footnote>†Email: mazou1979@yahoo.fr</footnote>

<footnote>‡Email: savadogoass302@yahoo.fr</footnote>
<table>
  <tbody>
    <tr>
      <td>5.2.4 Computing the exponentials of the blocks B<sub>q</sub> in the Jacobian matrix J<sub>II,k</sub>(X<sup>*</sup>)</td>
      <td>31</td>
    </tr>
    <tr>
      <td>5.2.5 An illustrative example</td>
      <td>32</td>
    </tr>
    <tr>
      <td><b>6 Conclusion</b></td>
      <td><b>35</b></td>
    </tr>
    <tr>
      <td><b>References</b></td>
      <td><b>37</b></td>
    </tr>
  </tbody>
</table>

Abstract

In the present paper, our goal is to establish a framework for the mathematical modelling and the analysis of the spread of an epidemic in a large population commuting regularly, typically along a time-periodic pattern, as is roughly speaking the case in populous urban center. We consider a large number of distinct homogeneous groups of individuals of various sizes, called subpopulations, and focus on the modelling of the changing conditions of their mixing along time and of the induced disease transmission. We propose a general class of models in which the ‘force of infection’ plays a central role, which attempts to ‘reconcile’ the classical modelling approaches in mathematical epidemiology, based on compartmental models, with some widely used analysis results (including those by P. van den Driessche and J. Watmough in 2002), established for apparently less structured systems of nonlinear ordinary-differential equations. We take special care in explaining the modelling approach in details, and provide analysis results that allow to compute or estimate the value of the basic reproduction number for such general periodic epidemic systems.

# 1 Introduction

Generally speaking, spatial or geographic heterogeneity plays an important role in the transmission process of many infectious diseases, as well as the temporal variations of the spreading conditions. The diversity of the corresponding (time and space) scales makes it a challenge to model, simulate and analyze the evolution of such phenomena. In the present paper, our goal is to establish a framework for the mathematical modelling and the analysis of the spread of an epidemic in a large population commuting regularly, typically along a time-periodic pattern, as is roughly speaking the case in populous urban center. The latter account for host displacements of varied nature according to the context, and typically distinguish the current location of the individuals, in order to define the possible infection paths.

A powerful framework to model the spread of infectious diseases among populations that are naturally partitioned into spatial sub-units, is based on *metapopulation models*, see [22, 30] and [5, 3] for a quite complete overview. Practically, ‘a metapopulation model involves explicit movements of the individuals between distinct locations’ [3], called *patches*. The fluxes between the different patches, seen as the vertices of an associated graph, are usually described through a discrete Laplacian matrix and take place continuously. The transfers between patches are generally instantaneous, but some contributions modelled specifically the infection occurring in transportation systems, see [31, 37] and references therein. At the price of an increased number of variables, some models keep track of the ‘origin’ of each individual (e.g. through a *residency patch*), but most of them just account for the number of individuals that are present at each location at the current time, making them otherwise indistinguishable. This omission is inadequate to describe precisely the complex, but quite regular, displacement schemes that shape urban commutations. In addition, the ‘residency patch’ is not sufficient to characterize univocally the history of the encounters, which are as many opportunities of inter-infection.

On the other hand, much work has been made to integrate various heterogeneity traits relevant to describe the spread of an infection within a population. Such approaches consider structured popula- tions subdivided into *subpopulations* that differ from each other [25]. As written by H.W. Hethcote and J.W. van Ark [25], the latter ‘can be determined not only on the basis of disease-related factors such as mode of transmission, latent period, infectious period, and genetic susceptibility or resistance, but also on the basis of social, cultural, economic, demographic, and geographic factors. In a spatially heteroge- neous population the subpopulations can be different schools, neighborhoods, cities, states, countries, or continents.’ This also includes the heterogeneity of the contact rates and of the behaviours that impact the disease spread. Attempts have been made to distinguish between ‘social’ encounters, occurring be- tween individuals of the same subpopulation (on a somehow ‘regular’ or ‘predictable’ basis), and random encounters between individuals of possibly different subpopulations [39]. Also, some structured models accommodate geographical aspects, in the way of metapopulation [38].
<figure><img src="image_1.png" /><figcaption>Figure 1: A population with 22 individuals, divided into 5 subpopulations $p_i$, $i = 1, \dots, 5$. Their cardinals are respectively $|p_1| = 6$, $|p_2| = 1$, $|p_3| = 5$, $|p_4| = 3$, $|p_5| = 7$. The set of subpopulations is $\mathcal{P} = \{p_1, p_2, p_3, p_4, p_5\}$, of cardinal $|\mathcal{P}| = 5$.</figcaption></figure>

An important characteristic of the 'subpopulations' considered in this type of models is their *homo-geneity*: the members of each subpopulation are perfectly mixed, submitted to the same 'encounters' with other subpopulations, and they cannot usually change from one subpopulation to the other (so that in particular, age classes are usually not treated as subpopulations). In other words, within each sub-population, the individuals are indistinguishable from one another from the point of view of the epidemic dynamics.

Generally speaking, in these settings, the movements are rather thought of as permanent and sta-tionary, rather than time-varying, phenomena. On the other hand, the evolution of disease transmission is usually subject to time fluctuations, and in particular periodic fluctuations [45]. The latter may result from seasonal variations induced by the climatic conditions, but also from human activities at various scales (daily commuting, opening and closing of schools, vaccination programs, etc.). Efforts have been made to analyze the effects of such forcing conditions on infection spread [6, 9, 7]. Due to the dif-ficulty to analyze the behaviour of periodic nonlinear systems in general, several simplifications have been made, amounting to analyze *stationary* systems obtained from a (usually informal) averaging of the periodic effects. However, examples in [45] have established that such simplifications may produce overly optimistic or overly pessimistic estimation of the occurrence of an epidemic (overestimation or underestimation of the basic reproduction number). Modelling and analysis of time-periodic epidemic models therefore remain of great interest, compared to the use of 'residence time' stationary models. This point is already true for the determination of the basic reproduction number, deduced from the model behaviour in the vicinity of a disease-free equilibrium; but is even more necessary for important issues outside the linearity range of the infection, such as e.g. the determination of the epidemic final size [36, 8].

Commuting involves complex interactions between population groups through complicated, but mostly repetitive, spatial and temporal patterns of encounters. In the present paper, in order to account for such situations, we shift the paradigm from a *spatial* 'Eulerian' description of the population and its displacements, to a 'Lagrangian' description in terms of *meeting and mixing* of homogeneous groups followed along their evolution. To carry out this program, we introduce (a large number of) distinct *subpopulations*, which are homogeneous groups of individuals of various sizes that constitute a partition of the total population. To fix ideas, here a subpopulation represents typically the group of all individuals living in a same location A that go at the same time to another location B (which represents e.g. school, work, shops, but also means of public transportation...), and then to a third location C at the same moment, etc. A key point is that all individuals of a given subpopulation have the same 'history'; that is, all are in contact at every moment with exactly the same other subpopulations. This is similar in spirit with the homogeneity assumption in the subpopulation models mentioned above. The example of a subpopulation with five subpopulations is provided in Figure 1, and two different partitions are shown in Figure 2.

In a complementary way, we assume that in each location, the different subpopulations are *perfectly mixed*. This assumption will allow to describe the epidemiological evolution at each fixed location where a fixed mixture of subpopulations through a *unique* aggregate compartmental model. Of course, this
<figure><img src="image_2.png" /><figcaption>Figure 2: Two different partitions $\mathcal{R}_1, \mathcal{R}_2$ of the set of subpopulations $\mathcal{P}$ shown in Figure 1. On the left, the set of classes $\mathcal{Q}_1 := \mathcal{P} \setminus \mathcal{R}_1 = \{(p_1, p_2, p_3), (p_4, p_5)\}$ contains $|\mathcal{Q}_1| = 2$ elements, namely the class $q_1 := \{p_1, p_2, p_3\}$ with $|q_1| = 3$ elements and the class $q_2 := \{p_4, p_5\}$ with $|q_2| = 2$ elements. On the right, $\mathcal{Q}_2 := \mathcal{P} \setminus \mathcal{R}_2 = \{(p_1), (p_2, p_5), (p_3, p_4)\}$ contains $|\mathcal{Q}_2| = 3$ elements, namely the classes $q_1 := \{p_1\}$ with $|q_1| = 1$ elements, $q_2 := \{p_2, p_5\}$ with $|q_2| = 2$ elements and $q_3 := \{p_3, p_4\}$ with $|q_3| = 2$ elements. In both configurations all individuals from a same subpopulation are together in the same class, and one may check that the balance formula (46) given later holds in both cases: $|\mathcal{P}| = |q_1| + |q_2| + |q_3| + |q_4| + |q_5|$. This is a consequence of the fact that in a partition, every subpopulation pertains exactly to one class.</figcaption></figure>

evolution is a priori different in different locations, due to different population densities and different behaviors modifying the contagion rate (e.g. at home, in public transportation...), and also due to the meeting of different subpopulations (e.g. at work, at school, during seasonal holidays...). This framework is quite distinct from most models that use subpopulations to represent heterogeneity traits: here the main difference between the subpopulations lies in the history of their contacts with other subpopulations, and in the local conditions of infection transmission, with otherwise identical characteristics. On the other hand, the modelling of heterogeneous traits may be achieved as usual through pertinent choice of the state variable.

The implementation of this program requires two steps: first to model the alternating mixture of different subpopulations; and second to describe how the epidemic spreads in a homogeneous melange of subpopulations, in a given location during a given time interval. On the one hand, the commuting and the induced changes in the contacts between subpopulations are modelled as piecewise-constant, periodic, mixing of the latter. We abandon any reference to 'physical locations', in favour of a vision in terms of *partitions* of the set of the different subpopulations. In this sense, trains or metros are also 'locations', and 'physical locations' that are actually empty (schools at night...) do not have to appear. On the other hand, to describe how the disease spreads locally, we propose an original class of compartmental models in which the *force of infection* $\lambda$ is explicitly distinguished and infection occurs by application of the latter to some (usually susceptible, but not only) compartments. Generally speaking, the force of infection is defined as 'the per capita rate at which susceptible individuals contract the infection' [30]. While the latter is usually taken as proportional either to the number of infected or to the proportion of infected, more complex models have been used, e.g. in order to take into account protective measures and intervention policies, or to reflect the necessity of multiple contacts to transmit the disease, see [15, 33, 43, 47, 18, 18, 32, 40, 34, 17]. We capture explicitly the force of infection in the class of models proposed below. This induces a special shape of the right-hand sides of the ODE models, which are jointly *linear with respect to the state variable x and affine with respect to the force of infection* $\lambda$. Notice that we consider in general several infective forces, first of all according to every location, but also to every susceptible compartment in the case of 'multi-group population', like done e.g. in [41, Section 23.1]. This leads, to a *vector-valued* force of infection.

Of course, the value of the force of infection is as usual a (nonlinear) function $\lambda(x)$ of the state variable $x$. The framework adopted here offers the advantage that, when two subpopulations described respectively by state vectors $x_1, x_2$ come into contact (i.e. pertain to the same class in the current partition), the resulting force of infection is simply $\lambda(x_1 + x_2)$ and applies identically to each of them, with in fact no other interaction between the two subgroups.
This feature has an important consequence. Many models of mathematical epidemiology describe the evolution of the epidemiological status of a given population in terms of proportions (of susceptible individuals, of infected individuals, etc.), allowing in particular to assume constant recruitment terms. Here, due to our modelling choice, it is necessary to have the ability to sum up the values at each compartment of the different populations. As a result, the models introduced here manipulate exclusively extensive quantities, and it is usually necessary to substitute constant recruitment terms by expressions proportional to the size of each subpopulation, in order to allocate realistically the demographic effects.

Up to our knowledge, the use of periodic epidemiological model structured in subpopulations with the purpose of accounting for commuting is original. In contrast with the models mentioned above, here all subpopulations will have essentially the same behaviour with respect to the infection spread: they differ only through the history of their displacements and the local definition and intensity of the force of infection, and subpopulations that pertain to the same equivalence class of the partition at a certain time participate jointly and identically to the spreading.

We display in the present paper the proposed modelling framework and provide some key analysis results extending the approach in [42]. This has to be seen as a preliminary step towards more complex results that will fully exploit the modelling capabilities. The modelling approach based on compartmental models has produced numerous variants in mathematical epidemiology [14, 19, 30, 13]. Attempts have been made to obtain general classes regrouping related models, liable to common general analysis results. See for example multi-group models [41, 21], metapopulation models [44, 4, 3], models of differential susceptibility and infectivity [27, 12].

In a sense, the present work constitutes an attempt to unify them, and to ‘reconcile’ the classical modelling approaches with some analysis results [42, 45] established for apparently less structured and less specific systems of equations. It is an effort to contribute to rationalizing and systematizing the modelling and analysis of epidemics, in line with the work [24].

The paper is organised as follows. We display first in Section 2 the representation of the changing melange at the scale of the global population. We then propose in Section 3 a class of epidemiological models that describe locally the infectious processes, in which the force of infection is explicitly depicted. We show by ample examples that it includes various compartmental models previously shown in the literature. These two scales of description are built in a consistent way, and are put together in Section 4, which summarizes the key modelling contribution of the paper, consisting of an original class of models of epidemics in commuting populations. Adequate assumptions (numbered (B.1), (B.2), (B.3)) are presented, and it is shown how they allow to recover essentially the framework developed by van den Driessche and Watmough [42] and to compute the basic reproduction number $R_0$ of an epidemic in a given, unique, subpopulation under stationary conditions (Theorem 2). The general analysis results are then given in Section 5, indicating how to compute the basic reproduction number for the proposed general class under periodic commutations and providing some applications (Theorem 4, Corollary 5 and Lemma 1). Concluding remarks are provided in Section 6.

**Notations** A list of the notations specific to the present paper is provided in Table 1. We recall here some usual mathematical notions.

• By definition let the vector $1_n \in \mathbb{R}^n$, $n \in \mathbb{N} \setminus \{0\}$ be such that $1_n := (1 \dots 1)$. • The componentwise order relation on any Cartesian power of the set $\mathbb{R}$ is denoted as usual $\ge$. The relation $>$ means ‘larger but not equal’, and the notation $>$ is used to denote the strict ordering of every component ($0 < x \in \mathbb{R}^n$ means $0 < x_i$, $i = 1, \dots, n$). • The spectral radius of a real square-matrix $M$ (that is the maximum of the moduli of its eigenvalues) is denoted by $\rho(M)$; its stability modulus (that is the maximum real part of its eigenvalues) is denoted by $s(M)$.

## 2 Representation of commuting populations

In order to represent commuting populations, we first introduce in Section 2.1 a notion of partitioned-population model, with several subpopulations evolving in parallel, in contact or out of contact. We then introduce in Section 2.2 a way to represent the commutations. This allows to define models of commuting populations in Section 2.3, which are basically partitioned-population models where the contact pattern between subpopulations is varying with time.
<table><caption>Table 1: Nomenclature. The cardinals of the sets $P$, $Q$ are written $|P|$, $|Q|$, and the cardinal of a class $q \in Q$ is written $|q|$. Without creating confusion, the same letter $i$ is used to index the subpopulations and the compartments in the state vector of each subpopulation. The notations and explanations related to the model itself are to be found in the text.</caption><tr><td>x</td><td>State variable of the epidemiological model</td></tr><tr><td>x<sup>I</sup></td><td>Infected compartments of the state variable</td></tr><tr><td>x<sup>U</sup></td><td>Non-infected compartments of the state variable</td></tr><tr><td>λ</td><td>Vector of the forces of infection</td></tr><tr><td>n</td><td>Dimension of the state variable x</td></tr><tr><td>n<sub>I</sub></td><td>Dimension of the infected compartment state x<sup>I</sup></td></tr><tr><td>n<sub>U</sub></td><td>Dimension of the non-infected compartment state x<sup>U</sup></td></tr><tr><td>n<sub>λ</sub></td><td>Dimension of the vector of the forces of infection λ</td></tr><tr><td>X<sub>s</sub></td><td>Set of all disease-free states (X<sub>s</sub> ⊂ R<sub>+</sub><sup>n</sup>)</td></tr><tr><td>X<sub>eq</sub></td><td>Set of all disease-free equilibrium points (X<sub>eq</sub> ⊂ X<sub>s</sub>)</td></tr><tr><td>P</td><td>Set of the subpopulations</td></tr><tr><td>R</td><td>Partition of P (possibly time-dependent)</td></tr><tr><td>Q</td><td>Quotient set P \ R of the equivalence classes ('loci') of subpopulations by R</td></tr><tr><td>X</td><td>Global state variable (of dimension |P|n)</td></tr><tr><td>p</td><td>Subpopulation name</td></tr><tr><td>[p], [p]<sub>R</sub></td><td>Class of the subpopulation p</td></tr><tr><td>q</td><td>Class name</td></tr><tr><td>i</td><td>Subpopulation index (i ∈ {1, ..., |P|}) or compartment index (i ∈ {1, ..., n})</td></tr><tr><td>j</td><td>Class index</td></tr><tr><td>m</td><td>Interswitching interval number (in periodic evolution)</td></tr><tr><td>k</td><td>Interswitching interval index</td></tr></table>

## 2.1 Partitioned-population models

A basic compartmental model describing the epidemiological evolution within an isolated, homogeneous and perfectly mixed population, is given [14, 19, 30, 13] as a dynamical system on $R_+^n$ defined by a system of ordinary differential equations

$$ \dot{x} = f(x(t), \lambda(t)), \quad \lambda(t) := \lambda(x(t)), \qquad (1) $$

for $f: R_+^n \times R_+^{n\lambda} \to R^n$, $\lambda: R_+^n \to R_+^{n\lambda}$. The vector $x$ bears the *numbers of individuals* in each of the model compartments, and the vector $\lambda$ gathers *forces of infection* exerted on the components of $x$. *Frequency-dependent* transmission mechanisms are typically rendered by a function $\lambda$ positively homogeneous of degree 0, and *density-dependent* transmission mechanisms by function positively homogeneous of positive degree<fnref n="1" />, but as mentioned in Section 1, other more complex expressions may be found in the literature. The dimension $n_\lambda$ of the latter may be larger than 1, see examples in Section 3.

**Definition 1 (Solution of a partitioned-population model).** Let $\mathcal{P}$ be a finite set, called the set of subpopulations, and $\mathcal{R} \subset \mathcal{P} \times \mathcal{P}$ be an equivalence relation on $\mathcal{P}$. For any subpopulation $p \in \mathcal{P}$, the corresponding equivalence class of $p$ in the quotient set $Q := \mathcal{P}\backslash\mathcal{R}$ is denoted $[p]_{\mathcal{R}}$, simply abbreviated $[p]$ when no confusion is possible. We call $(\mathcal{P}, \mathcal{R})$-solution of model (1) any continuous function $\mathcal{X}: R_+ \to R_+^{n|\mathcal{P}|}$, $\mathcal{X} := (x_1, \dots, x_{|\mathcal{P}|})$, continuously differentiable on $R_+$ such that

$$ \dot{x}_p = f(x_p(t), \lambda_{[p]}(t)), \quad \lambda_{[p]}(t) := \lambda_{[p]}(x_{[p]}(t)), \quad x_{[p]} := \sum_{p' \in [p]} x_{p'}, \quad p \in \mathcal{P}. \qquad (2) $$

$\mathcal{P}$ represents the index of different subpopulations, while the equivalence relation $\mathcal{R}$ defines which of them may infect each other. Therefore, one may consider each equivalence class as a *location*, namely

<footnote n="1">¹Recall that (in an adequate setting) a function $g$ is called positively homogeneous of degree $d$ if $g(\alpha x) = \alpha^d g(x)$ for any positive scalar $\alpha$ and any $x$ in its domain.</footnote>
the place where are present every subpopulation pertaining to this class. However, what is done here is not exactly to define the *place* where is located a given population at a given time, but rather which are the other populations with which it interacts. In particular, this perspective suppresses the necessity of any absolute notion of (geographical) location: after a switch in the repartition of the populations, the new partition is a priori not related to the previous one, and the equivalence classes change and become incomparable without supplementary information.

In (2), the conditions of transmission of the infection at time $t$ in a certain location $[p]$ depend upon the total population present therein at that time, which is exactly what is written $x_{[p]}(t)$. These conditions may be different from one place to the other: this is rendered by the dependence of $\lambda_{[p]}$ upon $[p]$. On the other hand, the *mechanisms of the infection* themselves are the same in every location, and for this reason the function $f$ is the same at every location. Recall that, as mentioned in Section 1, the definition of $x_{[p]}$ in (2) makes it essential to place in the components of $x$ *extensive quantities* such as numbers of individuals, and not relative proportions or densities.

Another point worth noting is that the coupling of two different subpopulations $p, p' \in \mathcal{P}$ that are present together in the same place (i.e. $[p] = [p']$) is reflected exclusively through the shared force of infection.

Summing up the contributions of every population within a given class yields:

$$ \dot{x}_{[p]} = \sum_{p' \in [p]} f(x_{p'}(t), \lambda_{[p]}(t)), \quad \lambda_{[p]}(t) = \lambda_{[p]}(x_{[p]}(t)), \quad p \in \mathcal{P}. \quad (3) $$

If $f$ is linear with respect to its first argument (an assumption that will be explained and adopted in Section 3), the effective number of equations in (3) is the number of components of the partition $\mathcal{R}$, that is the cardinal $|\mathcal{Q}|$ (rather than $|\mathcal{P}|$), times the state dimension $n$ of the model (1). As a matter of fact, one then obtains by adding all equations within a given equivalence class

$$ \dot{x}_{[p]} = f(x_{[p]}(t), \lambda_{[p]}(t)), \quad \lambda_{[p]}(t) = \lambda_{[p]}(x_{[p]}(t)), \quad p \in \mathcal{P}, $$

or more simply:

$$ \dot{x}_q = f(x_q(t), \lambda_q(t)), \quad \lambda_q(t) = \lambda_q(x_q(t)), \quad q \in \mathcal{Q}. \quad (4) $$

This describes the evolution in any 'location' $q$.

Let $q \in \mathcal{Q}$ be given and assume $x_p(0) \in \mathbb{R}_+^n$ given, for any $p \in \mathcal{P}$ such that $[p] = q$. From the solution of the nonlinear problem (4), one may compute the solution of the different interacting subpopulations $p$ such that $[p] = q$, as follows.

1. Compute

$$ x_q(0) = \sum \{x_p(0) : p \in \mathcal{P}, [p] = q\}. \quad (5) $$

2. Solve (4) with the initial condition (5). This provides $x_q(t)$, and as a by product, this also provides $\lambda_q(t) = \lambda_q(x_q(t))$.

3. Then, for any $p \in \mathcal{P}$ such that $[p] = q$, one has

$$ \dot{x}_p = f(x_p(t), \lambda_q(t)), \quad x_p(0) \text{ known.} $$

Solving the Cauchy problem for this *time-varying* ordinary differential equation, one obtains the solution $x_p$ for any $p \in \mathcal{Q}$. More specifically, as $f$ is linear with respect to its first argument, the solution yields the following algebraic relation

$$ x_p(t) = \Phi_q(t)x_p(0), \quad t \ge 0, \quad (6) $$

where $\Phi_q(\cdot)$ is the fundamental matrix of the *linear*, *time-varying*, equation

$$ \dot{y} = f(y, \lambda_q(t)), \quad \lambda_q(t) = \lambda_q(x_q(t)), $$

that is by definition

$$ \dot{\Phi}_q = f(\Phi_q(t), \lambda_q(t)), \quad t \ge 0, \quad \Phi_q(0) = I_n. \quad (7) $$

By convention, the result of applying $f$ to the matrix $\Phi_q(t)$ in the right-hand side of (7), is the matrix obtained as the concatenation of the vectors that result from applying $f$ to each of the column vectors of $\Phi_q(t)$. In other words, the $i$-th column of $\Phi_q(t)$ contains the value at time $t$ of the trajectory initiated on the $i$-th vector of the canonical basis at time 0.
The previous consideration is quite important from the point of view of numerical simulation. As a matter of fact, solving (2) as a system of $|q|$ coupled equations engages $|q|n$ scalar variables for any subpopulation class $q$; while computing jointly the solution of (4) with initial condition (5) amounts to solve a system of ODEs with $n$ scalar variables, and solving (7) (in sequence or in parallel with the previous one) necessitates $n^2$ scalar variables, independently of the size of the class $q$. This number is usually much smaller than the previous one. The same observation extends directly to the piecewise autonomous models presented below in Section 2.3.

Before going further, we provide now explicit expressions of $f$ and $\lambda$ for an elementary example.

**Example 1.** We consider the simple case where (1) is an SIR model with no induced mortality, that is

$$ \dot{S} = \mu N - \beta S \frac{I}{N} - \mu S, \quad \dot{I} = \beta S \frac{I}{N} - (\gamma + \mu) I, \quad \dot{R} = \gamma I - \mu R, \quad (8) $$

where $\beta, \mu, \gamma$ are positive parameters, and $N := S + I + R$ is the total population value, invariant along every trajectory due to the fact that $\dot{N} \equiv 0$. In order to comply with the analysis framework constructed by van den Driessche and Watmough [42], we put first in the state vector $x$ the (unique) infected compartment $I$, and take

$$ x := \begin{pmatrix} I \\ S \\ R \end{pmatrix}, \quad f(x, \lambda) := \begin{pmatrix} \lambda S - (\gamma + \mu) I \\ -\lambda S + \mu(I + R) \\ \gamma I - \mu R \end{pmatrix} = \begin{pmatrix} -(\gamma + \mu) & \lambda & 0 \\ \mu & -\lambda & \mu \\ \gamma & 0 & -\mu \end{pmatrix} x, \quad \lambda := \beta \frac{I}{S + I + R}. \quad (9) $$

The previous notation showcases the linear dependence of the function $f$ with respect to its first argument $x$ and its affine dependence with respect to $\lambda$. These properties are one of the key elements of the framework introduced later in Section 3.

The situation corresponding to system (2) (or (4)) is retrieved by introducing the function $f$ in (9) together with the definition of $\lambda_q$, in such a way that, for any $p \in \mathcal{P}$, the evolution of disease analogous to system (2) is given by:

$$ x_p := \begin{pmatrix} I_p \\ S_p \\ R_p \end{pmatrix}, \quad \dot{x}_p = \begin{pmatrix} -(\gamma + \mu) & \lambda_{[p]} & 0 \\ \mu & -\lambda_{[p]} & \mu \\ \gamma & 0 & -\mu \end{pmatrix} x_p, \quad \lambda_{[p]} = \beta_{[p]} \frac{I_{[p]}}{S_{[p]} + I_{[p]} + R_{[p]}}, \quad t \ge 0, \quad p \in \mathcal{P}, \quad (10) $$

for given positive parameters $\beta_q, q \in \mathcal{Q}$. The latter may be identical (homogeneous transmission conditions), or not (heterogeneous transmission conditions). In (10), in agreement with equation (2), one has $I_{[p]} := \sum_{p' \in [p]} I_{p'}$, and similarly for $S_{[p]}$ and $R_{[p]}$.

## 2.2 Admissible switching signals

In order to introduce the notion of commuting population models in the next Section 2.3, we first define the key notion of admissible switching signals.

**Definition 2 (Admissible switching signals).** Let $\mathcal{P}$ be a finite set of subpopulations. We call admissible switching signal any function $\mathcal{R}(\cdot)$ which associates to any $t \ge 0$ an equivalence relation (or equivalently a partition) on $\mathcal{P}$, denoted $\mathcal{R}(t)$, and such that $\mathcal{R}(\cdot)$ is constant on any bounded set of $\mathbb{R}^+$, except on a finite number of points (called switches).

For any admissible switching signal $\mathcal{R}$, we denote $\mathcal{I}$ the (open) set of time instants on which $\mathcal{R}$ is continuous. Otherwise said, $\mathbb{R}_+ \setminus \mathcal{I}$ is the set of switching times.

For any admissible switching signal, there exists a finite or denumerable number of switches, and the set $\mathcal{I}$ is a (finite or denumerable) union of successive open intervals sharing common endpoints.

In order to represent commutations, our interest lies especially in switching signals showing repeating pattern. For this, we will consider in the sequel admissible switching signals $\mathcal{R}(t)$ which are $T$-periodic for some $T > 0$. Denoting $m \in \mathbb{N}$ the number of commutations of such a signal $\mathcal{R}(\cdot)$ within each period, we will represent the corresponding switching pattern by an element of some set $\Lambda_m, m > 0$, defined as

$$ \Lambda_m = \{(\tau_0, \tau_1, \dots, \tau_m) : 0 = \tau_0 < \tau_1 < \dots < \tau_m = T\}. \quad (11) $$
For convenience, we define also the $m$ unions $\mathcal{I}_k$ of open intervals, $k = 1, \dots, m$, by

$$ \forall t \ge 0, \quad t \in \mathcal{I}_k \Leftrightarrow t \in (\tau_{k-1}, \tau_k) \quad \text{mod } T. $$

One then has $\mathcal{I} = \bigcup_{k=1}^{m} \mathcal{I}_k$, for $\mathcal{I}$ defined in Definition 2.

As no commutation occurs during the intervals contained in $\mathcal{I}_k$, we may denote unambiguously

* $\mathcal{R}(\mathcal{I}_k)$ the partition in force at any $t \in \mathcal{I}_k$, instead of $\mathcal{R}(t)$;

* $[p]_{\mathcal{R}(\mathcal{I}_k)}$ the corresponding class of a subpopulation $p \in \mathcal{P}$, instead of $[p]_{\mathcal{R}(t)}$;

* and $\mathcal{Q}(\mathcal{I}_k) = \mathcal{P} \setminus \mathcal{R}(\mathcal{I}_k)$ the set of subpopulation classes on the set $\mathcal{I}_k$, instead of $\mathcal{Q}(t)$, $t \in \mathcal{I}_k$.

## 2.3 Models of commuting populations

We are finally in position to introduce the models of commuting populations. Consider the following system of equations

$$ \dot{x}_p = f(x_p(t), \lambda_{[p]_{\mathcal{R}(t)}}(t)), \quad \lambda_{[p]_{\mathcal{R}(t)}}(t) = \lambda_{[p]_{\mathcal{R}(t)}} \left( x_{[p]_{\mathcal{R}(t)}}(t) \right), \quad x_{[p]_{\mathcal{R}(t)}} := \sum_{p' \in [p]_{\mathcal{R}(t)}} x_{p'}, \quad t \in \mathcal{I}, \quad p \in \mathcal{P}, \tag{12a} $$

$$ x_p(t^+) = x_p(t^-), \quad t \in \mathbb{R}_+ \setminus \mathcal{I}, \quad p \in \mathcal{P}. \tag{12b} $$

We call epidemiological model in commuting population the system (12).

**Definition 3 (Solution of a model of commuting population).** Let $\mathcal{P}$ be a finite set of subpopulations, and $\mathcal{R}(\cdot)$ an admissible switching signal. We call $(\mathcal{P}, \mathcal{R}(\cdot))$-solution of model (1) any continuous function $\mathcal{X} : \mathbb{R}_+ \to \mathbb{R}_+^{n|\mathcal{P}|}$, $t \mapsto (x_1(t), \dots, x_{|\mathcal{P}|}(t))$, continuously differentiable on $\mathcal{I}$ that fulfils (12).

In (12), $[p]_{\mathcal{R}(t)}$ denotes for any subpopulation $p \in \mathcal{P}$ the equivalence class to which $p$ belongs at time $t$, that is the ‘place’ where it is located, together with the other members of the same equivalence class $[p]_{\mathcal{R}(t)}$. The difference between (2) and (12) is only in the fact that the partition $\mathcal{R}$ in (12) is now a function of time.

Due to the non-accumulation of the number of switches (which comes from the admissibility of $\mathcal{R}$), it is easy to show that there exists a unique $(\mathcal{P}, \mathcal{R}(\cdot))$-solution of model (1) passing through a given initial condition $\mathcal{X}(0) \in \mathbb{R}_+^{n|\mathcal{P}|}$.

If $f$ is linear with respect to its first argument, writing together the evolution of the populations within each class yields the following equation, which is the analogue of (4) for system (12):

$$ \dot{x}_q = f(x_q(t), \lambda_q(t)), \quad \lambda_q(t) = \lambda_q(x_q(t)), \quad t \in \mathcal{I}, \quad q \in \mathcal{Q}(t). $$

Here by definition we put

$$ \mathcal{Q}(t) := \mathcal{P} \setminus \mathcal{R}(t), \quad t \in \mathcal{I}. $$

The considerations on simulation issues made at the end of Section 2.1 apply directly to system (12), due to the continuity of the state variable at the switching instants.

**Example 2 (continuation of Example 1).** The example of partitioned system shown in Example 1 extends straightforwardly to a model of commuting population, using the definition of $f$ in (9), and defining, for a given admissible switching signal, the force of infection at every (‘location’) $q \in \mathcal{P} \setminus \mathcal{R}(t)$ at any given time $t \in \mathcal{I}$ by

$$ \lambda_q(x) := \beta_q \frac{I}{S + I + R}, \quad t \in \mathcal{I}, \quad q \in \mathcal{Q}(t). \tag{13} $$

# 3 Epidemiological dynamics

So far, we have introduced a description of the commutations and ‘mixing’ of different subpopulations. In the present section, we will first focus in Section 3.1 on the representation of the different steps of the epidemiological processes (infection, exposition, recovery, immunity...), following the framework introduced by van den Driessche and Watmough [42], further extended by Wang and Zhao [45] to periodic problems. As an illustration and in order to show the interest of the obtained class of models, we show in Sections 3.2 to 3.6 that it contains in particular the five relevant examples studied in [42].
## 3.1 A class of compartmental disease transmission models

To adapt the framework of [42], a first formality is to introduce a distinction between the ‘infected’ and ‘non-infected’ compartments (respectively with exponents *I* and *U*) and to separate accordingly the terms in the equations which correspond to infections. More precisely we write the local model (1) under the more precise form:

$$x := \begin{pmatrix} x^I \\ x^U \end{pmatrix}, \quad f(x, \lambda) := \begin{pmatrix} \mathcal{F}^I(x, \lambda) + \mathcal{V}^I(x, \lambda) \\ \mathcal{V}^U(x, \lambda) \end{pmatrix}, \qquad (14a)$$

$$n_I := \dim x^I, \quad n_U := \dim x^U, \quad n_I + n_U = n. \qquad (14b)$$

The components $\mathcal{F}^I(x, \lambda)$, $\mathcal{V}^I(x, \lambda)$, $\mathcal{V}^U(x, \lambda)$ are chosen according to the prescriptions made by van den Driessche and Watmough<fnref n="2" />: $\mathcal{F}^I(x, \lambda)$ gathers solely the rates of appearance of new infections in infectious compartments, and $\mathcal{V}^I(x, \lambda)$ the rates of transfer of individuals in and out of the infected compartments by any other means; while $\mathcal{V}^U(x, \lambda)$ expresses the rates of transfer of individuals in and out of the non-infected compartments.

In fact, more structure is present, as will be demonstrated in Sections 3.2 to 3.6 by the study of examples. A key feature of most classical models is that the force of infection, which measures the per capita probability of acquiring the infection [1], enters linearly in the infection terms; while the other terms present in the equations are linear with respect to the numbers of individuals in some compartments. In line with this observation, we will assume that the function $f(x, \lambda)$ is linear with respect to *x* and affine with respect to $\lambda$, and therefore that $\mathcal{F}^I(x, \lambda)$, $\mathcal{V}^I(x, \lambda)$ and $\mathcal{V}^U(x, \lambda)$ share the same properties. In particular, the nonlinear character of the system arises from the force of infection $\lambda(x)$, and from cross products between components of $\lambda(x)$ and *x*.

Due to the fact that infections appear only in presence of infected, we assume that the function $\mathcal{F}^I(x, \lambda)$ which bears the corresponding terms is linear upon $\lambda$ (and not simply affine), and consequently bilinear in $(x, \lambda)$. The terms contained in $\mathcal{F}^I$ corresponds to the appearance of new infected individuals, therefore they always take on nonnegative values. They are compensated by the disparition of other (non-infected, but also infected in case of super-infections) individuals in corresponding compartments. We also assume that the term $\mathcal{V}^I(x, \lambda)$ is linear in $x^I$ only, meaning that it does not depend upon $x^U$. This formalizes the fact that every transfer to an infected compartments coming from a non-infected one, occurs through an infection process and therefore is counted within $\mathcal{F}^I(x, \lambda)$.

Consequently we denote more specifically

$$\mathcal{F}^I(x, \lambda) := F_+^I(\lambda)x^I + F_+(\lambda)x^U, \quad \mathcal{V}^I(x, \lambda) := V^I x^I - F_-^I(\lambda)x^I, \quad \mathcal{V}^U(x, \lambda) := V^U x - F_-(\lambda)x^U \quad (14c)$$

where $V^I \in \mathbb{R}^{n_I \times n_I}$, $V^U \in \mathbb{R}^{n_U \times n_U}$, and where the functions $F_{\pm}^I(\cdot) : \mathbb{R}^{n_\lambda} \to \mathbb{R}^{n_I \times n_I}$, $F_+(\cdot) : \mathbb{R}^{n_\lambda} \to \mathbb{R}^{n_I \times n_U}$, $F_-(\cdot) : \mathbb{R}^{n_\lambda} \to \mathbb{R}^{n_U \times n_U}$ are linear. In (14c), the new infections are decomposed into infection of susceptible individuals and infection of already infected ones, respectively by the nonnegative terms $F_+(\lambda)x^U$ and $F_-(\lambda)x^I$. They result in the remotion of corresponding infected individuals. This appears through the nonpositive term $-F_-^I(\lambda)x^I$ in $\mathcal{V}^I(x, \lambda)$ in case of infection of an already infected, and through the nonpositive term $-F_-(\lambda)x^U$ in $\mathcal{V}^U(x, \lambda)$ in case of infection of a non-infected. Notice also that these two matrices have to be diagonal, while in absence of self-infections, the diagonal of $F_+^I(\lambda)$ has to be zero.

Denoting

$$F_0(\cdot) := F_+^I(\cdot) - F_-^I(\cdot), \qquad (15)$$

one is finally led to adopt the general form

$$f(x, \lambda) := \begin{pmatrix} (F_0(\lambda) + V^I)x^I + F_+(\lambda)x^U \\ V^U x - F_-(\lambda)x^U \end{pmatrix}. \qquad (16)$$

This decomposition is indeed the most general form that integrates linear dependence with respect to *x* and affine dependence with respect to $\lambda$ within the framework based on the distinction of the new infection terms introduced in [42]. To summarize the comments made below, it is reasonable to assume

<footnote n="2">²The right-hand side decomposition is denoted $f(x) = \mathcal{F}(x) - \mathcal{V}(x)$ in [42]. For simplicity, we choose here the opposite sign for the term $\mathcal{V}$. This has no major impact, including on the subsequent definition of the basic reproduction number in (36a), as the spectral radius of a matrix is an even function.</footnote>
that $F_+(\cdot), F_-(\cdot)$ take on nonnegative values with $F_-(\cdot)$ is diagonal, and that $F_0(\cdot)$ is a Metzler matrix: this will be the content of Assumption (B.1) below.

Compartmental models describe transfers of individuals between different compartments. Except for the processes of birth and death which may add or subtract individuals, the corresponding transitions are essentially lossless. As a consequence, mass balance equations are usually fulfilled in absence of specific, disease-induced, mortality, namely:

$$ \mathbf{1}_{n_I}^{\mathrm{T}} F_0(\cdot) = 0_{n_I}^{\mathrm{T}}, \qquad \mathbf{1}_{n_I}^{\mathrm{T}} F_+(\cdot) = \mathbf{1}_{n_U}^{\mathrm{T}} F_-(\cdot), \qquad (17a) $$

$$ \mathbf{1}_{n_I}^{\mathrm{T}} V^I + \mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} = 0_{1 \times n_I}, \qquad (17b) $$

$$ \mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} = 0_{1 \times n_U}. \qquad (17c) $$

These formulas express the conservation of individuals during the transfers from one compartment to others during infections, for any value of the force of infection $\lambda$ (formula (17a)); and during the non-infective transfers: formula (17b) concerns the transfers of infected individuals and (17c) the transfers of non-infected individuals. As a complementary explanation, consider the total population for system (1)-(14), which is equal to $\mathbf{1}_{n_I}^{\mathrm{T}} x = \mathbf{1}_{n_I}^{\mathrm{T}} x^I + \mathbf{1}_{n_U}^{\mathrm{T}} x^U$. The derivative of this quantity with respect to time is

$$ \begin{aligned} \mathbf{1}_{n_I}^{\mathrm{T}} \dot{x}^I + \mathbf{1}_{n_U}^{\mathrm{T}} \dot{x}^U &= \mathbf{1}_{n_I}^{\mathrm{T}} (F_0(\lambda) + V^I) x^I + (\mathbf{1}_{n_I}^{\mathrm{T}} F_+(\lambda) - \mathbf{1}_{n_U}^{\mathrm{T}} F_-(\lambda)) x^U + \mathbf{1}_{n_U}^{\mathrm{T}} V^U x \\ &= \left[ \mathbf{1}_{n_I}^{\mathrm{T}} (F_0(\lambda) + V^I) + \mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} \right] x^I \\ &\quad + \left[ \mathbf{1}_{n_I}^{\mathrm{T}} F_+(\lambda) - \mathbf{1}_{n_U}^{\mathrm{T}} F_-(\lambda) + \mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} \right] x^U, \end{aligned} $$

so that $\mathbf{1}_{n_I}^{\mathrm{T}} \dot{x}^I + \mathbf{1}_{n_U}^{\mathrm{T}} \dot{x}^U = 0_n$ when all identities in (17a), (17b), (17c) are fulfilled, and the total population is then constant along the evolution.

In presence of specific, disease-induced, mortality, these formulas have to be modified to account for this effect, which is limited to the infected compartments. Equations (17a), (17c) remain valid, while (17b) has to be replaced by

$$ \mathbf{1}_{n_I}^{\mathrm{T}} V^I + \mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} < 0_{1 \times n_I}. \qquad (17d) $$

Recall that in the previous formula, $<$ means ‘less or equal, but not equal’. As a matter of fact, the previous computation shows that one has in this case

$$ \mathbf{1}_{n_I}^{\mathrm{T}} \dot{x}^I + \mathbf{1}_{n_U}^{\mathrm{T}} \dot{x}^U = \left( \mathbf{1}_{n_I}^{\mathrm{T}} V^I + \mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} \right) x^I \le 0. $$

In any case, (17c) alone ensures that any disease-free evolution occurs under constant total population.

The inequalities presented above are included below in Assumption (B.2). Other properties are needed, which are quite natural in this setting. First, the square matrices $V^I$ and $V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}$ should be Metzler matrices, as they represent intakes from other compartments for the off diagonal terms, and outtakes to other compartments for the diagonal ones. The latter matrix accounts for the internal exchanges between non-infected compartments, governed by the equation $\dot{x}^U = V^U \begin{pmatrix} 0_{n_I} \\ x^U \end{pmatrix}$ in absence of infected. As the balance of these exchanges should be zero in order to sustain constant population levels, it is thus natural to assume that the matrix $V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}$ admits 0 as an eigenvalue, associated to the left eigenvector $\mathbf{1}_{n_U}^{\mathrm{T}}$.

Also, every non-infected compartment is assumed to receive intakes from some infected compartments, either ‘after infection’ or ‘before’ in the shape of births proportional to this compartment: otherwise it is not needed for the description of the disease spread. For this it is added in Assumption (B.2) that $\mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} > 0_{1 \times n_I}$. Recall that the symbol $>$ has been defined in the paragraph Notations at
the end of Section 1: $0 \prec x \in \mathbb{R}^n$ means $0 < x_i, i = 1, \dots, n$.

In order to exemplify the ideas behind the previous modelling choices, we first come back to the SIR model already commented. Afterwards, we study in details in Sections 3.2 to 3.6 the five epidemiological models shown as examples in [42, Section 4], with some slight adaptations. We then come back in Section 4 to the formal presentation of the proposed framework.

**Remark 1.** We find in [12, equation (1)] a general class of models with differential susceptibility and infectivity, and in [41, Section 23.1] a general class of models with multiple groups or populations (including e.g. multi-stage infection). Both show evolution governed by ODEs that are at the same time linear with respect to state variables, and affine with respect to forces of infection (which are themselves nonlinear functions of the state). With some slight modifications (in particular, making the recruitment terms proportional to the population), both classes can be embedded in the framework (14) by adequate choice of the x variable.

**Example 3** (end of Examples 1 and 2). The model (9) is embedded in the framework (14) with $x^I := I$, $x^U := (S \ R)^T$, so that $n = 3$, $n_I = 1$, $n_U = 2$, and one sees that $n_\lambda = 1$. The function $f$ defined in (9) is decomposed as in (14a) in its two, infected and non-infected, components, yielding

$$ \mathcal{F}^I(x, \lambda) := \lambda S, \quad \mathcal{V}^I(x, \lambda) := -(\gamma + \mu)I, \quad \mathcal{V}^U(x, \lambda) := \begin{pmatrix} \mu(I + R) - \lambda S \\ \gamma I - \mu R \end{pmatrix}, $$

and the decomposition (14c) with the following elements:

$$ F_0(\lambda) := 0, \quad F_+(\lambda) := \lambda \begin{pmatrix} 1 & 0 \end{pmatrix}, \quad V^I := -(\gamma + \mu), \quad (18a) $$

$$ V^U := \begin{pmatrix} \mu & 0 & \mu \\ 0 & 0 & -\mu \end{pmatrix}, \quad F_-(\lambda) := \lambda \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}. \quad (18b) $$

One sees easily that the identities in (17a) and (17c) are verified, and as well (17b), because

$$ \mathbf{1}_{n_I}^\top V^I + \mathbf{1}_{n_U}^\top V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} = -(\gamma + \mu) + \begin{pmatrix} 1 & 1 \end{pmatrix} \begin{pmatrix} \mu \\ \gamma \end{pmatrix} = 0, $$

## 3.2 Treatment model

The first example considered in [42] is a tuberculosis model with treatment proposed in [11, 16]. It possesses four compartments, namely, individuals susceptible to tuberculosis (S), exposed individuals (E), infectious individuals (I) and treated individuals (T). The forces of infection applied to susceptible and treated individuals have different rates, respectively $\beta_1 \frac{I}{N}$ and $\beta_2 \frac{I}{N}$, where $N := E + I + S + T$. The exposed individuals enter the infectious compartment at a rate $\nu$. All newborns are susceptible, and the mortality is uniform and equal to $d > 0$. The treatment rates are denoted $r_1$, resp. $r_2$, for the exposed individuals, resp. the infectious individuals. A fraction $q$ of the treatments of infectious individuals are successful, and the unsuccessfully treated infectious individuals reenter the exposed compartment. The epidemiological model is given by the following system.

$$ \left\{ \begin{aligned} \dot{E} &= \beta_1 S \frac{I}{N} + \beta_2 T \frac{I}{N} - (d + \nu + r_1)E + (1 - q)r_2 I, \\ \dot{I} &= \nu E - (d + r_2)I, \\ \dot{S} &= b(N) - dS - \beta_1 S \frac{I}{N}, \\ \dot{T} &= -dT + r_1 E + qr_2 I - \beta_2 T \frac{I}{N}. \end{aligned} \right. \quad (19) $$

In the case where the natality term is linear $b(N) = bN$, this decomposes as in (16), with state
dimension $n = 4$ and the following dynamics.

$$x^I := \begin{pmatrix} E \\ I \end{pmatrix}, \quad x^U := \begin{pmatrix} S \\ T \end{pmatrix}, \quad n_I = n_U := 2, \quad \lambda := \begin{pmatrix} \lambda_1 \\ \lambda_2 \end{pmatrix}, \quad n_\lambda := 2, \quad (20a)$$

$$F_0(\cdot) \equiv 0_{2 \times 2}, \quad F_+(\lambda) = F_-(\lambda) := \lambda_1 \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} + \lambda_2 \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}, \quad (20b)$$

$$V^I := \begin{pmatrix} -(d + \nu + r_1) & (1-q)r_2 \\ \nu & -(d + r_2) \end{pmatrix}, \quad V^U := \begin{pmatrix} b & b & b-d & b \\ r_1 & qr_2 & 0 & -d \end{pmatrix}, \quad (20c)$$

$$\lambda_1(x) = \beta_1 \frac{I}{E + I + S + T}, \quad \lambda_2(x) = \beta_2 \frac{I}{E + I + S + T}. \quad (20d)$$

One checks that $V^I$ is a nonsingular Metzler matrix, that the identity in (17a) is fulfilled and that (17c) holds when the natality and mortality rates $b$ and $d$ are equal. Whenever $b = d$, (17b) (resp. (17d)) is fulfilled when the specific mortality $d$ is zero (resp. positive).

## 3.3 Multigroup model

The second example in [42] is an $m$-group SIRS-vaccination model of Hethcote [25, 23]. The model is slightly adapted here. Using numbers of individuals instead of proportions, the model equations are as follows :

$$\left\{
\begin{aligned}
\dot{I}_i &= \sum_{j=1}^{m} \beta_{ij}(x) S_i \frac{I_j}{I_j + S_j + R_j} - (d_i + \gamma_i + \epsilon_i) I_i, \\
\dot{S}_i &= (1 - p_i) b_i (I_i + S_i + R_i) - (d_i + \theta_i) S_i + \sigma_i R_i - \sum_{j=1}^{m} \beta_{ij}(x) S_i \frac{I_j}{I_j + S_j + R_j}, \\
\dot{R}_i &= p_i b_i (I_i + S_i + R_i) + \gamma_i I_i + \theta_i S_i - (d_i + \sigma_i) R_i,
\end{aligned}
\right.
\quad (21)$$

for $i = 1, \dots, m$ where $x = (I_1 \dots I_m S_1 \dots S_m R_1 \dots R_m)^T$. The natural death rate is $d_i$, and $\epsilon_i$ the specific death rate. Recovering of the infected individuals occurs at rate $\gamma_i$, and immunity wanes at a rate $\sigma_i$. All newborns are susceptible, with group $i$ having a birth rate equal to $b_i$. Here, the dimension of the state variable is $n = 3m$, where $m$ is the number of distinct groups. Vaccination of the newborns, resp. of the susceptibles, is achieved with a rate $p_i$, resp. $\theta_i$. The incidence rates $\beta_{ij}(x)$ depend on individual behaviour and describes the amount of mixing between groups [28].

The decomposition (16) is obtained here with

$$x^I := (I_1 \dots I_m)^T, \quad x^U := (S_1 \dots S_m R_1 \dots R_m)^T, \quad n_I = n_\lambda := m, n_U := 2m, \quad (22a)$$

$$F_0(\cdot) \equiv 0_{m \times m}, \quad F_+(\lambda) := (\Lambda \quad 0_{m \times m}), \quad F_-(\lambda) := \text{diag}\{\Lambda; 0_{m \times m}\}, \quad (22b)$$

$$V^I := \text{diag}\{-(d_i + \gamma_i + \epsilon_i)\}, \quad (22c)$$

$$V^U := \begin{pmatrix} \text{diag}\{(1-p_i)b_i\} & \text{diag}\{(1-p_i)b_i - d_i - \theta_i\} & \text{diag}\{(1-p_i)b_i + \sigma_i\} \\ \text{diag}\{p_i b_i + \gamma_i\} & \text{diag}\{p_i b_i + \theta_i\} & \text{diag}\{p_i b_i - d_i - \sigma_i\} \end{pmatrix}, \quad (22d)$$

$$\Lambda(x) := \text{diag}\{\lambda_i(x)\} \in \mathbb{R}^{m \times m}, \quad \lambda_i(x) := \sum_{j=1}^{m} \beta_{ij}(x) \frac{I_j}{I_j + S_j + R_j}, \quad i = 1, \dots, m. \quad (22e)$$

Again, $V^I$ is a nonsingular Metzler matrix. The two identities in (17a) are clearly fulfilled. On the other hand, here

$$\begin{aligned}
\mathbf{1}_{n_I}^{\mathrm{T}} V^I + \mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} &= (b_1 - d_1 - \epsilon_1 \quad \dots \quad b_m - d_m - \epsilon_m), \\
\mathbf{1}_{n_U}^{\mathrm{T}} V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} &= (b_1 - d_1 \quad \dots \quad b_m - d_m \quad b_1 - d_1 \quad \dots \quad b_m - d_m),
\end{aligned}$$

so that (17c) holds when $b_i = d_i$, $i = 1, \dots, m$; and (17b), resp. (17d) also holds when $\epsilon_i = 0$, $i = 1, \dots, m$, resp. $(\epsilon_1 \dots \epsilon_m) > (0 \dots 0)$.
## 3.4 Staged progression model

The third example shown in [42] has a single susceptible compartment, but the infected individuals pass through *m* different stages of disease [26]. The model has dimension *m* + 1, with the following equations.

$$
\left\{
\begin{aligned}
\dot{I}_1 &= S \sum_{k=1}^{m-1} \frac{\beta_k I_k}{N} - (\nu_1 + d_1)I_1, \\
\dot{I}_i &= \nu_{i-1} I_{i-1} - (\nu_i + d_i)I_i, \quad i = 2, \dots, m-1, \\
\dot{I}_m &= \nu_{m-1} I_{m-1} - d_m I_m, \\
\dot{S} &= bN - bS - S \sum_{k=1}^{m-1} \frac{\beta_k I_k}{N},
\end{aligned}
\right.
$$

where the total population is $N := S + I_1 + \dots + I_m$.

Denoting $e_1$ the first vector of the canonical basis of $\mathbb{R}^m$, this model writes as (16) with

$$x^I := (I_1 \quad \dots \quad I_m)^\top, \quad x^U = S, \quad n_I = n_\lambda := m, \ n_U := 1, \qquad (23a)$$

$$F_0(\cdot) \equiv 0, \quad F_+(\lambda) := \lambda e_1, \quad F_-(\lambda) := \lambda, \qquad (23b)$$

$$V^I := \begin{pmatrix} -(\nu_1 + d_1) & 0 & \cdots & \cdots & 0 \\ \nu_1 & -(\nu_2 + d_2) & 0 & \cdots & 0 \\ 0 & \nu_2 & \cdots & \cdots & \vdots \\ \vdots & \vdots & \ddots & -(\nu_{m-1} + d_{m-1}) & 0 \\ 0 & 0 & \cdots & \nu_{m-1} & -d_m \end{pmatrix}, \quad V^U := (b\mathbf{1}_m^\top \quad 0), \qquad (23c)$$

$$\lambda(x) := \sum_{k=1}^{m-1} \frac{\beta_k I_k}{N}. \qquad (23d)$$

The matrix $V^I$ is a nonsingular Metzler matrix. Identities (17a) and (17c) are fulfilled, and (17b) as well when $d_1 = \dots = d_m = b$ (i.e. when all mortality rates are identical and equal to the natality rate). Alternatively, (17d) is fulfilled when $d_i \ge b, i = 1, \dots, m$ (every mortality rate at most equal to the natality rate).

## 3.5 Multistrain model

The fourth example from [42] is a multistrain model adapted from [16, 20]. The model has a unique susceptible compartment, but two infectious compartments corresponding to the two infectious agents, with possible super-infections by strain one of individuals already infected by strain two, giving rise to a new infection in compartment $I_1$. The model is the following 3-dimensional model, where the total population is $N := S + I_1 + I_2$.

$$
\left\{
\begin{aligned}
\dot{I}_1 &= \beta_1 \frac{I_1}{N} S - (b + \gamma_1) I_1 + \nu \frac{I_1}{N} I_2, \\
\dot{I}_2 &= \beta_2 \frac{I_2}{N} S - (b + \gamma_2) I_2 - \nu \frac{I_1}{N} I_2, \\
\dot{S} &= bN - bS + \gamma_1 I_1 + \gamma_2 I_2 - \frac{\beta_1 I_1 + \beta_2 I_2}{N} S.
\end{aligned}
\right.
$$

The setting in (16) is recovered by putting

$$x^I := \begin{pmatrix} I_1 \\ I_2 \end{pmatrix}, \quad x^U = S, \quad n_I = n_\lambda := 2, \ n_U := 1, \quad \lambda(x) = \begin{pmatrix} \lambda_1(x) \\ \lambda_2(x) \end{pmatrix}, \qquad (24a)$$

$$F_0(\lambda) := \frac{\nu}{\beta_1} \lambda_1 \begin{pmatrix} 0 & 1 \\ 0 & -1 \end{pmatrix}, \quad F_+(\lambda) := \lambda_1 \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \lambda_2 \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \quad F_-(\lambda) := \lambda_1 + \lambda_2, \qquad (24b)$$

$$V^I := \begin{pmatrix} -(b + \gamma_1) & 0 \\ 0 & -(b + \gamma_2) \end{pmatrix}, \quad V^U := (b + \gamma_1 \quad b + \gamma_2 \quad 0), \qquad (24c)$$

$$\lambda_1(x) = \beta_1 \frac{I_1}{S + I_1 + I_2}, \quad \lambda_2(x) = \beta_2 \frac{I_2}{S + I_1 + I_2}. \qquad (24d)$$
As before, $V^I$ is a nonsingular Metzler matrix and (17a), (17b) and (17c) are fulfilled. Notice that here the functions $F_+^I$ and $F_-^I$ that appear in the decomposition (14c) are not zero, so that $F_0$ is not identically zero. This is due to the fact that infections may occur for already infected individuals.

## 3.6 Vector-host model

The last example from [42] is a simplified version of a vector-host model borrowed from [20]. The latter couples a simple SIS model for the hosts with an SI model for the vectors. It has four compartments: infected hosts ($I$) and vectors ($V$), and susceptible hosts ($S$) and vectors ($M$). Inter-species transmission occurs in both directions through infective contacts between vectors and hosts. The frequency of the (infectious) contacts is proportional to the number of vectors, with different rates $\beta_s$ and $\beta_m$, depending on the direction. The equations are as follows.

$$ 
\left\{ 
\begin{aligned}
\dot{I} &= \beta_s \frac{S}{S+I} V - (b+\gamma)I, \\
\dot{V} &= \beta_m \frac{I}{S+I} M - cV, \\
\dot{S} &= b(S+I) - bS + \gamma I - \beta_s \frac{S}{S+I} V, \\
\dot{M} &= c(M+V) - cM - \beta_m \frac{I}{S+I} M,
\end{aligned}
\right. 
\qquad (25) $$

Model (25) is 4-dimensional and rewrites

$$ x^I := \begin{pmatrix} I \\ V \end{pmatrix}, \quad x^U := \begin{pmatrix} S \\ M \end{pmatrix}, \quad n_I = n_U = n_\lambda := 2, \quad \lambda := \begin{pmatrix} \lambda_1 \\ \lambda_2 \end{pmatrix}, \qquad (26a) $$

$$ F_0(\cdot) \equiv 0_{2 \times 2}, \quad F_+(\lambda) = F_-(\lambda) := \lambda_1 \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} + \lambda_2 \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}, \qquad (26b) $$

$$ V^I := \begin{pmatrix} -(b+\gamma) & 0 \\ 0 & -c \end{pmatrix}, \quad V^U := \begin{pmatrix} b+\gamma & 0 & 0 & 0 \\ 0 & c & 0 & 0 \end{pmatrix}, \qquad (26c) $$

$$ \lambda_1(x) = \beta_s \frac{V}{S+I}, \quad \lambda_2(x) = \beta_m \frac{I}{S+I}. \qquad (26d) $$

For this example too, one shows easily that $V^I$ is a nonsingular Metzler matrix and (17a), (17b) and (17c) are fulfilled.

# 4 Models of epidemics in commuting populations

For clarity, we first recall and summarize here the proposed class of models (Section 4.1). The disease-free states and in particular the disease-free equilibria are studied in Section 4.2. In Section 4.3, we state a list of adequate assumptions, which basically allows to retrieve the framework put in place by van den Driessche and Watmough [42]. It is then shown in Section 4.4 how these assumptions allow to find the basic reproduction number in the case of a single population (i.e. for $|P| = 1$), which is the case treated in the latter paper.

## 4.1 A class of epidemiological models in commuting populations

For sake of clarity, we repeat and summarize here in a compact way the ingredients of the proposed epidemiological models in commuting population.

Assume be given matrices

$$ V^I \in \mathbb{R}^{n_I \times n_I}, \quad V^U \in \mathbb{R}^{n_U \times n_U}, \qquad (27a) $$

linear maps

$$ F_0(\lambda) : \mathbb{R}^{n_\lambda} \to \mathbb{R}^{n_I \times n_I}, \quad F_+(\lambda) : \mathbb{R}^{n_\lambda} \to \mathbb{R}^{n_I \times n_U}, \quad F_-(\lambda) : \mathbb{R}^{n_\lambda} \to \mathbb{R}^{n_U \times n_U}, \qquad (27b) $$
and a function $\lambda : \mathbb{R}_+^n \to \mathbb{R}_+^{n_\lambda}$, define

$$f(x, \lambda) := \begin{pmatrix} (F_0(\lambda) + V^I)x^I + F_+(\lambda)x^U \\ V^U x - F_-(\lambda)x^U \end{pmatrix}. \qquad (27c)$$

When needed, we will use the following decomposition, deduced from (27) and quite similar to the one in [42]:

$$\mathcal{F}^I(x, \lambda) := F_+^I(\lambda)x^I + F_+(\lambda)x^U, \quad \mathcal{V}^I(x, \lambda) := V^I x^I - F_-^I(\lambda)x^I, \quad \mathcal{V}^U(x, \lambda) := V^U x - F_-(\lambda)x^U, \quad (28a)$$

where by definition

$$F_-^I(\lambda) \text{ is the main diagonal of } -F_0(\lambda), \quad F_+^I(\lambda) := F_0(\lambda) + F_-^I(\lambda). \qquad (28b)$$

The diagonal of the matrix $F_+^I(\lambda)$ is thus zero by construction.

## 4.2 Disease-free states and disease-free equilibrium points

First, as in [42], we define $\mathbf{X}_s$ to be the set of all disease-free states:

$$\mathbf{X}_s := \left\{ \begin{pmatrix} 0_{n_I} \\ x^U \end{pmatrix} : x^U \in \mathbb{R}_+^{n_U} \right\}. \qquad (29)$$

Another important subset of the state space is the set of disease-free equilibria (DFE), a subset of the previous one. Here, contrary to [42], we put in this category any equilibrium point, whatever its stability properties. In absence of infection, one has $\lambda = 0_{n_\lambda}$, and the set $\mathbf{X}_{eq} \subset \mathbf{X}_s$ of the DFEs of system (12) is given by

$$\mathbf{X}_{eq} := \{x \in \mathbf{X}_s : f(x, 0_{n_\lambda}) = 0_n\}. \qquad (30)$$

Due to the assumed property of linearity of $f$ with respect to $x$, the set $\{x \in \mathbb{R}^n : f(x, 0_{n_\lambda}) = 0_n\}$ is a vectorial space, and the set $\mathbf{X}_{eq}$ is therefore a cone. Assume $\mathbf{X}_{eq}$ is not empty, so that there exists (at least) one DFE $y^* > 0$ of (1). Then, for any given total population number of $N > 0$ (non-infected) individuals, there exists a DFE $x^* > 0$ such that $\mathbf{1}_n^\top x^* = N$, namely

$$x^* := \frac{N}{\mathbf{1}_n^\top y^*} y^*,$$

or equivalently, considering only the non-infected compartments,

$$x^{U*} := \frac{N}{\mathbf{1}_{n_U}^\top y^{U*}} y^{U*}.$$

As a matter of fact, for $x^*$ defined as before, $f(x^*, 0_{n_\lambda}) = \frac{N}{\mathbf{1}_n^\top y^*} f(y^*, 0_{n_\lambda}) = 0_n$ by linearity, so that $x^*$ is a DFE; and $\mathbf{1}_n^\top x^* := \frac{N}{\mathbf{1}_n^\top y^*} \mathbf{1}_n^\top y^* = N$.

Let us now characterize in a more precise manner the set $\mathbf{X}_{eq}$ and the condition $\mathbf{X}_{eq} \neq \emptyset$. For any $x \in \mathbf{X}_s$, one has $x^I = 0_{n_I}$ and $\lambda(x) = 0_{n_\lambda}$, so that

$$f(x, 0_{n_\lambda}) := \begin{pmatrix} (F_0(0_{n_\lambda}) + V^I)x^I + F_+(0_{n_\lambda})x^U \\ V^U x - F_-(0_{n_\lambda})x^U \end{pmatrix} = \begin{pmatrix} 0_{n_I} \\ V^U x \end{pmatrix}.$$

Therefore, the cone $\mathbf{X}_{eq}$ is characterized by the following property:

$$\mathbf{X}_{eq} = \mathbf{X}_s \cap \ker V^U. \qquad (31)$$

## 4.3 Assumptions on the models

We make the following assumptions on the objects defined in Section 4.1.
**Assumption (B.1)** (On the operators $F_0, F_±$). The operators $F_0(·), F_+(·), F_-(·)$ defined in (27b) are linear. Moreover, for any $λ ∈ R_+^{n_λ}$, the matrix $F_0(λ)$ is a Metzler matrix and $F_+(λ), F_-(λ)$ are nonnegative, with $F_-(λ)$ diagonal. Last,

$$ \mathbf{1}_{n_I}^\top F_0(\cdot) = 0_{n_I}^\top, \quad \mathbf{1}_{n_I}^\top F_+(\cdot) = \mathbf{1}_{n_U}^\top F_-(\cdot). \qquad (32) $$

**Assumption (B.2)** (On the matrices $V^I, V^U$). For the matrices $V^I, V^U$ defined in (27a),

* $V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}$ is a Metzler matrix of which $0$ is a simple eigenvalue associated to the left-eigenvector $\mathbf{1}_{n_U}^\top$;

* $V^I$ is a Metzler matrix, $V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} > 0_{n_U \times n_I}$ and
$$ -\mathbf{1}_{n_I}^\top V^I \ge \mathbf{1}_{n_U}^\top V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} > 0_{1 \times n_I}. \qquad (33) $$

**Assumption (B.3)** (On the map $λ$). The map $λ$ takes on nonnegative values, and is zero on $X_s$.

Explanations of the meaning of Assumptions (B.1) and (B.2) have been given in Section 3.1, and Assumption (B.3) requires no particular comment.

One may check without major difficulties that the examples displayed in Section 3 fit Assumptions (B.1) to (B.3). (See however an important nuance in Remark 3 on the simplicity of the eigenvalue $0$ (Assumption (B.2)), in relation with the dimension of the set of disease-free equilibrium points.)

Notice that (32) repeats (17a) without modification. Also, the link between (33) and (17b)/(17d) is evident, while (17c) is contained in the first point of Assumption (B.2). Before going further, we enlighten the relation of the previous assumptions with the framework invented by van den Driessche and Watmough [42].

**Proposition 1.** Assume $f$ defined by (27) fulfils Assumptions (B.1), (B.2), (B.3). Then the same is true for Assumptions (A.1), (A.2), (A.3), (A.4) in [42].

As will be seen below, the present setting does not contain an hypothesis on asymptotic stability as strong as [42, Assumption (A.5)]. Indeed, Assumption (B.3) allows for a unique zero eigenvalue in the spectrum of the Jacobian matrix of the system in disease-free evolution. This comes from the fact that, due to the second formula in (33), the population is constant along any disease-free evolution, see the discussion in the beginning of Section 3. (Recall that extensive quantities are considered here; see the related discussion in Section 1.) This insensitivity to the total population level generates a zero eigenvalue in the Jacobian spectrum at any equilibrium point.

*Proof of Proposition 1.*

• We first retrieve from (27) the setting of [42]. For this, let $\mathcal{F}^I(x, λ)$, $\mathcal{V}^I(x, λ)$, $\mathcal{V}^U(x, λ)$ be deduced from model (27) thanks to (28). From Assumption (B.2), one deduces that for any $λ ∈ R_+^{n_λ}$, $F_0(λ)$ is a Metzler matrix and its off-diagonal elements are nonnegative. By construction, the latter are also the off-diagonal elements of $F_+^I(λ)$, which possesses a zero diagonal. Thus $F_+^I(λ) \ge 0_{n_I \times n_I}$. One the other hand, $-F_-^I(λ)$ is a diagonal matrix that bears the diagonal elements of $F_0(λ)$. Due to Assumption (B.1), one has

$$ 0_{n_I}^\top = \mathbf{1}_{n_I}^\top F_0(\lambda) = \mathbf{1}_{n_I}^\top (F_+^I(\lambda) - F_-^I(\lambda)), $$

so that $\mathbf{1}_{n_I}^\top F_-^I(\lambda) = \mathbf{1}_{n_I}^\top F_+^I(\lambda) \ge 0$, and the diagonal matrix $F_-^I(\lambda)$ is nonnegative. Now, the decomposition of $f$ according to

$$ f(x, \lambda) = \begin{pmatrix} \mathcal{F}^I(x, \lambda) + \mathcal{V}^I(x, \lambda) \\ \mathcal{V}^U(x, \lambda) \end{pmatrix} $$

is identical to the one given in [42]. We will now establish that it fulfils the Assumptions (A.1) to (A.4) of the present paper.

• Due to Assumption (B.1), the function $\mathcal{F}^I(x, λ(x))$ in the decomposition (28), which gathers the infection terms, takes on nonnegative values. This is essentially the content of [42, Assumption (A.1)].
Assumption (A.3) in [42] states that no infection term feeds the non-infected compartments. This property is already fulfilled by the choice of the structure of model (27).

Translated in the present setting, Assumption (A.2) in [42] states that whenever $x_i = 0$ for some $i = 1, \dots, n$, then the corresponding component $V_i^I(x, \lambda(x)) \ge 0$. First let $i \in \{1, \dots, n\}$ be the index of an infected compartment, that is $i \in \{1, \dots, n_I\}$. On the one hand, due to the fact that $V^I$ is a Metzler matrix by Assumption (B.2), and that $x \ge 0$, one has $(V^I x^I)_i \ge 0$ whenever $x_i = 0$. On the other hand, the operator $F^I(\lambda)$ being diagonal, the $i$-th component of $F^I(\lambda)x^I$ is zero when $x_i = 0$. Therefore, as $V^I(x, \lambda) := V^I x^I - F^I(\lambda)x^I$, one obtains that $V_i^I(x, \lambda(x)) \ge 0$ whenever $x_i = 0$.

The same argument allows to treat the case of a non-infected compartment. Let $i \in \{1, \dots, n\}$ be the index of such a compartment, that is $i \in \{n_I + 1, \dots, n\}$, and $x \in \mathbb{R}_+^n$ such that $x_i = 0$. One has

$$V^U x = V^U \left( \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} x^I + \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} x^U \right) = V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} x^I + V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} x^U. \quad (34)$$

The $i$-th component of the first term is nonnegative due to the fact that $x^I \ge 0$ and $V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} > 0$ (Assumption (B.2)); while the $i$-th component of the second term is also nonnegative, due to the fact that $x_i = 0$ and $x^U \ge 0$ (Assumption (B.2)), $V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}$ is a Metzler matrix.

Overall, $(V^U x)_i \ge 0$. On the other hand, the operator $F_-(\lambda)$ is diagonal (Assumption (B.1)), so that the $i$-th component of $F_-(\lambda)x^U$ is zero when $x_i = 0$. Therefore, $V_i^U(x, \lambda) \ge 0$ whenever $x_i = 0$. This achieves the proof of Assumption (A.2) in [42].

Last, Assumption (A.4) in [42] states that at any disease-free state $x \in \mathbf{X}_s$, the infectious term $F^I(x, \lambda(x)) = F_+^I(\lambda)x^I + F_-(\lambda)x^U$ is equal to $0$ and that $V^I(x, \lambda(x)) = V^I x^I - F_-^I(\lambda)x^I \le 0$ (Assumption (B.1)). Let $x \in \mathbf{X}_s$, then $x^I = 0$ and all terms $F_+^I(\lambda)x^I$, $F_-^I(\lambda)x^I$ and $V^I x^I$ are thus $0$ (Assumption (B.3)), so by linearity the matrix $F_+(\lambda)$ is null, and thus $F_+(\lambda)x^U = 0$ (Assumption (B.3)), so that $V^I(x, \lambda(x)) = 0$ as well. This achieves the proof of Proposition 1. □

## 4.4 Basic reproduction number for a single-population model

The following result shows essentially that the previous setting contains the stationary, single-population, case treated in [42]. It gives the value of the basic reproduction number of model (1) with a dynamic defined by $f$ in (27). The proof consists mainly in adapting the framework of [42, Theorem 2].

For any sufficiently regular function $H(x) : \mathbb{R}^n \to \mathbb{R}^r$, we decompose the Jacobian matrix of $H$ in $\mathbb{R}^{r \times n}$ as

$$\frac{\partial H}{\partial x}(x) = \left( \frac{\partial H}{\partial x^I}(x) \quad \frac{\partial H}{\partial x^U}(x) \right), \quad \frac{\partial H}{\partial x^I}(x) \in \mathbb{R}^{r \times n_I}, \quad \frac{\partial H}{\partial x^U}(x) \in \mathbb{R}^{r \times n_U}. \quad (35a)$$

For simplicity, this will be written in the sequel with the more compact notation

$$\frac{\partial H}{\partial x}(x) = \partial_x H(x) = \left( \partial_{x^I} H(x) \quad \partial_{x^U} H(x) \right). \quad (35b)$$

**Theorem 2** (Basic reproduction number and stability for single-subpopulation epidemic models). Assume $f$ and $\lambda$ fulfil Assumptions (B.1), (B.2), (B.3). Then the following properties hold.

* System (1)-(27) admits nonzero disease-free equilibrium points.

* For any nonzero disease-free equilibrium point $x^* \in \mathbf{X}_{\text{eq}}$ of system (1)-(27), let

$$\mathcal{R}_0 := \rho(FV^{-1}), \quad (36a)$$

where $F, V \in \mathbb{R}^{n_I \times n_I}$ are given by:

$$F := \partial_{x^I} (F_+(\lambda(x))x^{U*})|_{x=x^*}, \quad V := V^I. \quad (36b)$$

The following assertions are true:

- if $\mathcal{R}_0 > 1$, then $x^*$ is unstable;
- if $\mathbb{R}_0 < 1$, then $x^*$ is marginally stable. Moreover, if the first identity in (33) holds with an equality, then $x^*$ is locally asymptotically stable.

In agreement with the notation defined in (35), the matrix $F$ of Theorem 2 is such that

$$F = \partial_{x^I} (F_+(\lambda(x))x^{U*})|_{x=x^*} = \frac{\partial (F_+(\lambda(x))x^{U*})}{\partial x^I}\Big|_{x=x^*} = \sum_{i=1}^{n_\lambda} \frac{\partial (F_+(\lambda)x^{U*})}{\partial \lambda_i}\Big|_{\lambda=0_{n_\lambda}} \frac{\partial \lambda_i(x)}{\partial x^I}\Big|_{x=x^*}.$$

### Proof of Theorem 2.

* From Assumption (B.2), the matrix $M := V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} \in \mathbb{R}^{n_U \times n_U}$ is a Metzler matrix, of which the positive vector $\mathbf{1}_{n_U}^T$ is a left-eigenvector associated to the eigenvalue 0. For large enough value of $\alpha$, $M + \alpha I_{n_U}$ is non-negative, and $\mathbf{1}_{n_U}^T$ is a left-eigenvector associated to the eigenvalue $\alpha$. The latter is thus the spectral radius of $M + \alpha I_{n_U}$ [10, Corollary 1.1.12], so that 0 is the stability modulus of $M$ and no other eigenvalue is located on the imaginary axis.

The matrix $M$ also admits a nonnegative right-eigenvector $w \in \mathbb{R}_+^{n_U}$, see [10, Theorem 1.1.1]. Now, the nonzero, nonnegative, vector $\begin{pmatrix} 0_{n_I} \\ w \end{pmatrix} \in \mathbb{R}_+^n$ is such that

$$V^U \begin{pmatrix} 0_{n_I} \\ w \end{pmatrix} = V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} w = Mw = 0_{n_U}.$$

One concludes, using (31), that

$$\begin{pmatrix} 0_{n_I} \\ w \end{pmatrix} \in \mathbf{X}_s \cap \ker V^U = \mathbf{X}_{\text{eq}},$$

which demonstrates the existence of a nonzero disease-free equilibrium of system (1)-(27).

* Generally speaking, the Jacobian matrix of the function $f(x, \lambda(x))$ for $f$ defined in (27) at any point $x \in \mathbb{R}_+^n$ is decomposed as

$$J(x) := \begin{pmatrix} J_{II}(x) & J_{IU}(x) \\ J_{UI}(x) & J_{UU}(x) \end{pmatrix}, \qquad (37)$$

$J_{II}(x) \in \mathbb{R}^{n_I \times n_I}$, $J_{IU}(x) \in \mathbb{R}^{n_I \times n_U}$, $J_{UI}(x) \in \mathbb{R}^{n_U \times n_I}$, $J_{UU}(x) \in \mathbb{R}^{n_U \times n_U}$, where

$$J_{II}(x) := \partial_{x^I} ((F_0(\lambda(x)) + V^I)x^I + F_+(\lambda(x))x^U),$$

$$J_{IU}(x) := \partial_{x^U} ((F_0(\lambda(x)) + V^I)x^I + F_+(\lambda(x))x^U),$$

$$J_{UI}(x) := \partial_{x^I} (V^U x - F_-(\lambda(x))x^U), \quad J_{UU}(x) := \partial_{x^U} (V^U x - F_-(\lambda(x))x^U).$$

Let $x^* \in \mathbf{X}_{\text{eq}}$ be a (nonzero) DFE. We will first show that at such a point, the block $J_{IU}(x^*)$ is null. As a matter of fact, one has $x^* = 0_{n_I}$, $\lambda(x^*) = 0_{n_\lambda}$, and thus the matrices $F_0(\lambda(x^*))$, $F_\pm(\lambda(x^*))$ are zero, due to the linearity of these operators, see Assumption (B.1). Therefore,

$$J_{IU}(x^*) = \partial_{x^U} (F_+(\lambda(x))x^U)|_{x=x^*} = F_+(\lambda(x^*)) + \partial_{x^U} (F_+(\lambda(x))x^{U*})|_{x=x^*} = \frac{\partial F_+}{\partial \lambda}(0_{n_\lambda}) \frac{\partial \lambda}{\partial x^U}(x^*) x^{U*}.$$

One has $\lambda(x) = 0_{n_\lambda}$ at any point $x \in \mathbf{X}_s$, so that

$$\frac{\partial \lambda}{\partial x^U}(x^*) = 0_{n_\lambda \times n_U},$$

and one concludes that at any DFE $x^*$,

$$J_{IU}(x^*) = 0_{n_I \times n_U}. \qquad (38)$$

Due to the previous property, the stability of any DFE $x^*$ is determined by the spectra of the two blocks $J_{II}(x^*)$ and $J_{UU}(x^*)$ of $J(x^*)$. Let us compute these quantities.

On the one hand, $\frac{\partial}{\partial x^U} (F_-(\lambda(x)))|_{x=x^*}$ is shown to be zero, just as $\frac{\partial}{\partial x^U} (F_+(\lambda(x)))|_{x=x^*}$ previously. Thus,

$$J_{UU}(x^*) = \partial_{x^U} (V^U x - F_-(\lambda(x))x^U)|_{x=x^*} = V^U \frac{\partial x}{\partial x^U} - F_-(\lambda(x^*)) = V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}, \quad (39)$$
which by Assumption (B.2) has zero as simple eigenvalue and possibly other eigenvalues with negative real parts.

On the other hand, using the nullity of the terms $F_0(\lambda(x))$ and $F_+(\lambda(x))$, one has

$$J_{II}(x^*) = \partial_{x^I} \left( (F_0(\lambda(x)) + V^I)x^I + F_+(\lambda(x))x^U \right) \Big|_{x=x^*} = V^I + \partial_{x^I} \left[ F_+(\lambda(x))x^{U*} \right] \Big|_{x=x^*}.$$

One checks directly that

$$\partial_{x^I} (F_+(\lambda(x))x^{U*}) \Big|_{x=x^*} = \partial_\lambda (F_+(\lambda)x^{U*}) \Big|_{\lambda=0_{n_\lambda}} \partial_{x^I} \lambda \Big|_{x=x^*},$$

where, due to the linearity of $F_+(\cdot)$, the quantity $\partial_\lambda (F_+(\lambda)x^{U*})$ does not depend upon $\lambda$. Comparing with the notations introduced in (36b), one thus obtains

$$J_{II}(x^*) = F + V. \qquad (40)$$

Notice first that the matrix $V = V^I$ is Hurwitz. As a matter of fact, by hypotheses in Assumption (B.2), one has

$$\mathbf{1}_{n_I}^\top V^I \le -\mathbf{1}_{n_U}^\top V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} < 0.$$

This, together with the fact that this matrix is Metzler, implies that $V^I$ is Hurwitz, see [35, Theorem 2.1], and in particular *invertible*.

On the other hand, by Assumption (B.1), the operator $F_+(\lambda)$ is positive. This implies that the terms $\partial_\lambda (F_+(\lambda)x^{U*})|_{\lambda=0_{n_\lambda}}$ are (constant) nonnegative matrices. Also, the terms $\partial_{x^I} \lambda|_{x=x^*}$ are nonnegative, as $\lambda(x)$ is zero at any point $x \in \mathbf{X}_s$ and takes on nonnegative values. Thus, $F$ is a *nonnegative matrix*. Arguing as in [42, Proof of Theorem 2], one then obtains that

$$s(J_{II}(x^*)) < 0 \Leftrightarrow \mathcal{R}_0 < 1 \quad \text{and} \quad s(J_{II}(x^*)) > 0 \Leftrightarrow \mathcal{R}_0 > 1, \qquad (41)$$

for the value $\mathcal{R}_0$ defined in (36), and where $s$ denotes the stability modulus, see the Notations in the end of the Introduction.

• Consider first the case where $\mathcal{R}_0 > 1$, so that $s(J_{II}(x^*)) > 0$. One has

$$s(J(x^*)) = \max\{s(J_{II}(x^*)), s(J_{UU}(x^*))\} = \max\{s(J_{II}(x^*)), 0\} = s(J_{II}(x^*)) > 0.$$

The DFE $x^*$ is then unstable.

• Consider now the case where $\mathcal{R}_0 < 1$. One has here $s(J_{II}(x^*)) < 0$, but $s(J_{UU}(x^*)) = 0$ (contrary to what happens in [42]), and $J_{UU}(x^*)$ has only 0 as an imaginary eigenvalue, which is simple. Let us first establish the (simple) stability of the DFE under the general condition (33). Let $\delta x := x - x^*$, $\delta x^I := x^I - x^{*I}$, $\delta x^U := x^U - x^{*U}$. By definition, the tangent linear equation writes

$$\delta \dot{x}^I = J_{II}(x^*) \delta x^I, \quad \delta \dot{x}^U = J_{UI}(x^*) \delta x^I + J_{UU}(x^*) \delta x^U.$$

Then exploiting the asymptotic stability of the block $J_{II}(x^*)$, for sufficiently small $\varepsilon > 0$, there exist $\zeta > 0$, $c \ge 1$, $k > 0$ such that, for any trajectory of (1)-(27) such that $\|\delta x(0)\| < \zeta$,

$$\|\delta x^I(t)\|, \|\delta \dot{x}^I(t)\| \le c e^{-kt} \|\delta x(0)\| \quad \text{for any } t \ge 0, \text{ as long as } \|\delta x(t)\| < \varepsilon. \qquad (42)$$

On the other hand, as seen in Section 3.1, one has for any trajectory

$$\mathbf{1}_{n_I}^\top \dot{x}^I + \mathbf{1}_{n_U}^\top \dot{x}^U = \left[ \mathbf{1}_{n_I}^\top (F_0(\lambda) + V^I) + \mathbf{1}_{n_U}^\top V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} \right] x^I \\ + \left[ \mathbf{1}_{n_I}^\top F_+(\lambda) - \mathbf{1}_{n_U}^\top F_-(\lambda) + \mathbf{1}_{n_U}^\top V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} \right] x^U,$$

so that in the vicinity of the DFE $x^*$, one has $\lambda(x^*) = 0_{n_\lambda}$ and

$$\mathbf{1}_{n_I}^\top \delta \dot{x}^I + \mathbf{1}_{n_U}^\top \delta \dot{x}^U = \left( \mathbf{1}_{n_I}^\top V^I + \mathbf{1}_{n_U}^\top V^U \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} \right) \delta x^I + \mathbf{1}_{n_U}^\top V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} \delta x^U.$$
Due to the first point in Assumption (B.2), the second term is zero, yielding

$$
\begin{align*}
\mathbf{1}_{n_U}^{\top} \delta \dot{x}^{U} &= -\mathbf{1}_{n_I}^{\top} \delta \dot{x}^{I} + \left( \mathbf{1}_{n_I}^{\top} V^{I} + \mathbf{1}_{n_U}^{\top} V^{U} \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} \right) \delta x^{I} \\
&= \left( -\mathbf{1}_{n_I}^{\top} J_{II} (x^{*}) + \mathbf{1}_{n_I}^{\top} V^{I} + \mathbf{1}_{n_U}^{\top} V^{U} \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} \right) \delta x^{I}.
\end{align*}
$$

As $\delta x^I$ decreases exponentially, see (42), one thus gets that,

$$
|\mathbf{1}_{n_U}^{\top} \delta x^U(t)| \leq c' \|\delta x(0)\| \quad \text{as long as } \|\delta x(t)\| < \varepsilon,
$$

for some $c' \ge 1$. Denoting $\tilde{x}^U := (x_{n_I+1} \dots x_{n-1})^\top$ and $\delta\tilde{x}^U := \tilde{x}^U - \tilde{x}^{*U}$, the evolution of $\delta\tilde{x}^U$ is then governed by the remaining eigenvalues of the matrix $V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}$ mentioned in Assumption (B.2). The latter have negative real parts, as noticed in the beginning of the present demonstration. Therefore, for some $c'' \ge 1$,

$$
\| \delta \tilde{x}^U(t) \| \leq c'' e^{-kt} \| \delta x(0) \| \quad \text{as long as} \quad \| \delta x(t) \| < \varepsilon.
$$

Take $\|\delta x(0)\| < \min \{\zeta; \frac{\varepsilon}{5c}; \frac{\varepsilon}{5c'}; \frac{\varepsilon}{5c''}\}$, then as long as $\|\delta x(t)\| < \varepsilon$, one deduces

$$
\| \delta x^I(t) \| < \frac{\varepsilon}{5}, \quad |\mathbf{1}_{n_U}^\top \delta x^U(t)| < \frac{\varepsilon}{5}, \quad \| \delta \tilde{x}^U(t) \| < \frac{\varepsilon}{5}.
$$

Using for example for $\| \cdot \|$ the $l^1$-norm, one has

$$
|\delta x_n(t)| = |\mathbf{1}_{n_U}^{\top} \delta x^U(t) - \mathbf{1}_{n_U-1}^{\top} \delta \tilde{x}^U(t)| \le |\mathbf{1}_{n_U}^{\top} \delta x^U(t)| + |\mathbf{1}_{n_U-1}^{\top} \delta \tilde{x}^U(t)| \le |\mathbf{1}_{n_U}^{\top} \delta x^U(t)| + \|\delta \tilde{x}^U(t)\| \le \frac{2\varepsilon}{5},
$$

so that

$$
\|\delta x(t)\| = \|\delta x^I(t)\| + \|\delta x^U(t)\| \le \|\delta x^I(t)\| + \|\delta \tilde{x}^U(t)\| + \|\delta x_n(t)\| < \frac{4\varepsilon}{5} < \varepsilon.
$$

Therefore $\|\delta x(t)\|$ never reaches the value $\varepsilon$, and the simple stability of the DFE $x^*$ is proved.

Assume now more specifically that the first formula in (33) holds with an equality, i.e.

$$
\mathbf{1}_{n_I}^{\top} V^{I} + \mathbf{1}_{n_U}^{\top} V^{U} \begin{pmatrix} I_{n_I} \\ 0_{n_U \times n_I} \end{pmatrix} = 0_{1 \times n_I}.
$$

Then along any trajectory of system (1)-(27), one has (see Section 3) $\mathbf{1}_n^\top \dot{x} = 0$, so that $\mathbf{1}_n^\top x(t) = \mathbf{1}_n^\top x(0)$. Replacing e.g. the coordinate $x_n(\cdot)$ by its value $\sum_{j=1}^n x_j(0) - \sum_{j=1}^{n-1} x_j(\cdot)$ yields an $(n-1)$-dimensional dynamical system, of which the vector $\tilde{x}^* := (x_1^* \dots x_{n-1}^*)$ is an equilibrium. The spectrum of the Jacobian matrix of this reduced dynamical system at this point consists of the $n-1$ eigenvalues of $J(x^*)$ with negative real parts. The point $\tilde{x}^*$ is thus a locally asymptotically stable equilibrium of the reduced system, and $x^*$ is therefore a locally asymptotically stable equilibrium of the initial system. This achieves the demonstration of Theorem 2. $\square$

**Remark 2.** The function whose Jacobian is considered in Proposition 1 is, as in [42, Assumption (A.5)], equal to the right-and side $f(x, \lambda(x))$ 'when the function $\mathcal{F}^I(x, \lambda(x))$ is set to zero'.

**Remark 3.** Consider the following system with two susceptible compartments, extending the SIR model (8) by introduction of inherited heterogeneity of susceptibility and infectivity:

$$
\dot{S}_i = \mu N_i - S_i \left( \beta_{ii} \frac{I_i}{S_i} + \beta_{ii'} \frac{I_{i'}}{N_{i'}} \right) - \mu S_i, \quad \dot{I}_i = S_i \left( \beta_{ii} \frac{I_i}{S_i} + \beta_{ii'} \frac{I_{i'}}{N_{i'}} \right) - (\gamma + \mu) I_i, \quad \dot{R}_i = \gamma I_i - \mu R_i,
$$

with $N_i = S_i + I_i + R_i$, for any $i, i'$ such that $\{i, i'\} = \{1, 2\}$. This model may be written in the preceding setting with $n_I = 2, n_U = 4$. One verifies that, with the state vector $(I_1 \ I_2 \ S_1 \ S_2 \ R_1 \ R_2)^\top$, one has (compare with (18))

$$
V^U = \begin{pmatrix} \mu & 0 & 0 & 0 & \mu & 0 \\ 0 & \mu & 0 & 0 & 0 & \mu \\ \gamma & 0 & 0 & 0 & -\mu & 0 \\ 0 & \gamma & 0 & 0 & 0 & -\mu \end{pmatrix}.
$$
Corresponding to the conservation of the two quantities $S_i+I_i+R_i$, respectively within the compartments $S_i/I_i/R_i$ for $i=1,2$, one has here

$$ (1 \ 0 \ 1 \ 0) V^U \begin{pmatrix} 0_{2 \times 4} \\ I_4 \end{pmatrix} = (0 \ 1 \ 0 \ 1) V^U \begin{pmatrix} 0_{2 \times 4} \\ I_4 \end{pmatrix} = 0_4^T. $$

These identities imply, but are stronger than the eigenvector property in the first point of Assumption (B.2), namely

$$ 1_4^T V^U \begin{pmatrix} 0_{2 \times 4} \\ I_4 \end{pmatrix} = 0_4^T. $$

Here the eigenvalue 0 of matrix $V^U \begin{pmatrix} 0_{2 \times 4} \\ I_4 \end{pmatrix}$ has thus multiplicity (at least) 2, so that Assumption (B.2) is not fulfilled and Theorem 2 does not apply directly. However, using the two conservation formulas of the quantities $N_i$, one adapts easily the demonstration, showing that two components of $\delta x^U$ are bounded instead of one, while the remaining components converge to zero. One then recovers the marginal stability property, and indeed the conclusions of Theorem 2 are still valid. This provides a substantial extension of this result.

Notice that the same phenomenon appears for the multi-group model presented in Section 3.3 and for the vector-host model in Section 3.6. For the first one, the disease-free equilibrium points may or may not present individuals in every non-infected compartments $S_i, R_i$; and in the second one, they may or may not contain both uninfected hosts $S$ and uninfected vectors $M$. Generally speaking this subtlety requires special attention when the dimension of the set $X_{eq}$ of all DFEs is at least equal to 2.

# 5 Basic reproduction number of epidemiological models in commuting populations

We gather in this section the main results of the paper. In Section 5.1 is given the basic reproduction number for a partitioned-population model (that is, without commutations). This serves as an introduction to Section 5.2, where the general case of commuting populations is treated.

## 5.1 Basic reproduction number of stationary partitioned-population models

Consider first the model of partitioned population given in equation (2), with a function $f$ as in (27). The corresponding model is stationary, and composed of $|P \setminus R|$ uncoupled ‘locations’ (classes).

Assuming the Assumptions (B.1) and (B.2) fulfilled, as well as (B.3) for every $\lambda_q, q \in Q$, it is easy to see that there exist disease-free equilibrium points to this system. As a matter of fact, equation (4) provides a description at the level of each location $q \in Q = P \setminus R$, and due to the previous assumptions, the latter admits disease-free equilibrium points $x_q^*, q \in Q$ (see Theorem 2). It is then easy to split the global population at this point and allocate proportions $\alpha_p x_q^*, 0 \le \alpha_p \le 1$, $\sum_{p \in q} \alpha_p = 1$ of the latter to each subpopulation $x_p, p \in q$, in each class $q$ in order to obtain a DFE $X^*$.

We will first introduce in Section 5.1.1 a representation of the state variable for a partitioned-population model. A difficulty is that, as we aim at obtaining tractable criteria, we can no longer use the intrinsic formulas in (2) expressed in terms of all subpopulations indistinctly: instead, one must consider the subpopulation classes and define a global ordering of the subpopulations that keep track of their contacts, in order to ease the representation of the global state variable. Once this done, one may compute the Jacobian matrix of the system (Section 5.1.2), and finally obtain in Section 5.1.3 the analysis results related to the stability of the DFE $X^*$.

### 5.1.1 Ad hoc representation of the state variable for partitioned-population model

We will now proceed to the computation of the Jacobian matrix of equation (2) at a disease-free equilibrium point. For this, we will choose a specific representation of the global state variable $X$, compatible with the previous setting. In order to picture the analysis results in terms of matrix properties, we will have to drop the intrinsic representation (2), expressed in terms of subpopulations, and argue in terms
of classes, ordered in an order having some specified properties. This complicates notably the writing of the conditions. However, this complication is more of a notational nature, than of a conceptual one.

Before displaying the new representation, notice that, in any case, the class $q_j \in \mathcal{Q}$ contains $|q_j|$ distinct subpopulations, so that the total number of subpopulations fulfils the identity

$$|\mathcal{P}| = \sum_{\{|q| : q \in \mathcal{Q}\}}. \qquad (43)$$

Also, due to the fact that $x_q = \sum_{p \in q} x_p$ for any class $q \in \mathcal{Q}$, one has for any $p, p' \in \mathcal{P}$;

$$\frac{\partial x_{[p']}}{\partial x_p} = I_n \quad \text{if} \quad [p] = [p'], \quad \frac{\partial x_{[p']}}{\partial x_p} = 0_{n \times n} \quad \text{if} \quad [p] \neq [p'], \qquad (44)$$

independently of the ordering of the components of the vector $\mathcal{X}$.

For simplicity, we will assume that the $n$-dimensional vectors $x_p$ are ordered within the $n|\mathcal{P}|$-dimensional vector $\mathcal{X}$ according to the following principles: all vectors $x_p^I$ are grouped first, then all vectors $x_p^U$; and within the vectors $x_p^I$, resp. $x_p^U$, the information of the subpopulation of a first class $q_1$ are put first, then of a second class $q_2$, and so on until $q_{|\mathcal{Q}|}$. Schematically, this corresponds to writing the state variable $\mathcal{X} \in \mathbb{R}_+^{n|\mathcal{P}|}$ as

$$\mathcal{X} = \begin{pmatrix} \mathcal{X}^I \\ \mathcal{X}^U \end{pmatrix}, \quad \mathcal{X}^I \in \mathbb{R}_+^{n_I |\mathcal{P}|}, \quad \mathcal{X}^U \in \mathbb{R}_+^{n_U |\mathcal{P}|}, \qquad (45a)$$

with

$$\mathcal{X}^{I\tau} = \left( x_{q_1, p_1}^{I\tau} \quad \dots \quad x_{q_1, p_{|q_1|}}^{I\tau} \quad x_{q_2, p_1}^{I\tau} \quad \dots \quad x_{q_2, p_{|q_2|}}^{I\tau} \quad \dots \quad x_{q_{|\mathcal{Q}|}, p_1}^{I\tau} \quad \dots \quad x_{q_{|\mathcal{Q}|}, p_{|q_{|\mathcal{Q}|}|}}^{I\tau} \right), \qquad (45b)$$

where the vector $x_{q_j, p_i}^I$, $j = 1, \dots, |\mathcal{Q}|$, $i = 1, \dots, |q_j|$, represents the $n_I$ infected compartments of the subpopulation $p_i$ in the class $q_j \in \mathcal{Q}$. The vector $\mathcal{X}^U$ is defined similarly. Notice that such a decomposition is not unique, as the ordering of the classes on the one hand, and the ordering of the subpopulations in every class on the other, are not unique. Once the ordering of the classes and the ordering of each subpopulation in its class have been chosen, then the ordering of the components in vector $\mathcal{X}$ is the corresponding lexicographic order: for any $q, q' \in \mathcal{Q}$, $p \in q, p' \in q'$, the components $x_{q,p}$ are located before the components $x_{q',p'}$ if and only if $q$ precedes $q'$, or if $q = q'$ and $p$ precedes $p'$.

**Remark 4.** With the convention exposed above, one may rewrite identity (43) as

$$|\mathcal{P}| = \sum_{j=1}^{|\mathcal{Q}|} |q_j|. \qquad (46)$$

### 5.1.2 Computation of the Jacobian matrix

With the convention of representation of the state variable $\mathcal{X}$ previously exposed, the Jacobian matrix $\mathcal{J}(\mathcal{X}^*)$ at a DFE $\mathcal{X}^*$ is written by blocks as

$$\mathcal{J}(\mathcal{X}^*) := \begin{pmatrix} \mathcal{J}_{II}(\mathcal{X}^*) & \mathcal{J}_{IU}(\mathcal{X}^*) \\ \mathcal{J}_{UI}(\mathcal{X}^*) & \mathcal{J}_{UU}(\mathcal{X}^*) \end{pmatrix}, \qquad (47)$$

for some $\mathcal{J}_{II}(\mathcal{X}^*) \in \mathbb{R}^{n_I |\mathcal{P}| \times n_I |\mathcal{P}|}$, $\mathcal{J}_{IU}(\mathcal{X}^*) \in \mathbb{R}^{n_I |\mathcal{P}| \times n_U |\mathcal{P}|}$, $\mathcal{J}_{UI}(\mathcal{X}^*) \in \mathbb{R}^{n_U |\mathcal{P}| \times n_I |\mathcal{P}|}$, $\mathcal{J}_{UU}(\mathcal{X}^*) \in \mathbb{R}^{n_U |\mathcal{P}| \times n_U |\mathcal{P}|}$. Let us assess the value of these blocks.

* Let $p \in \mathcal{P}$, then taking inspiration from the computations in the proof of Theorem 2, one has that, at any DFE $\mathcal{X}^*$, for any $p \in \mathcal{P}$,

$$\begin{aligned} & \partial_{x_p^I} \left[ (F_0(\lambda_{[p]}(x_{[p]})) + V^I) x_p^I + F_+(\lambda_{[p]}(x_{[p]})) x_p^U \right] \Big|_{\mathcal{X}=\mathcal{X}^*} \\ &= V^I + \partial_{x_p^I} \left[ F_+(\lambda_{[p]}(x)) x_p^{U*} \right] \Big|_{x=x_{[p]}^*} = V^I + \partial_{\lambda} \left[ F_+(\lambda) x_p^{U*} \right] \Big|_{\lambda=0_{n_{\lambda}}} \partial_{x_p^I} \lambda_{[p]}(x) \Big|_{x=x_{[p]}^*} . \end{aligned}$$

In coherence with the previous notations, $x_p^{U*}$ represents the uninfected components of the state vector of the subpopulation $p$ at equilibrium; and $x_{[p]}^* = \sum_{p' \in [p]} x_{p'}^*.$
**Remark 5.** As noticed before, the function $\lambda \mapsto \partial_{\lambda} [F_{+}(\lambda)x_{p}^{U*}]$ is constant, as it is the derivative of an affine map. Accordingly, here and in the sequel we will simply write

$$ \partial_{\lambda} [F_{+}(\lambda)x_{p}^{U*}] \quad \text{instead of} \quad \partial_{\lambda} [F_{+}(\lambda)x_{p}^{U*}] (0_{n_{\lambda}}). $$

Also, $x_{[p]}^* = \sum_{p' \in [p]} x_{p'}^*$, so that the preceding derivative does not depend upon the subpopulation of the class. One may thus write

$$ \partial_{x_p^I} \lambda_{[p]}(x_{[p]}^*) = \partial_{x^I} \lambda_{[p]}(x_{[p]}^*), $$

where the variable $x_p^I$ has been replaced by the dummy variable $x^I$.

On the other hand, for $p, p' \in \mathcal{P}$ such that $[p] = [p']$ but $p \neq p'$, one has

$$ \begin{aligned} & \partial_{x_{p'}^I} [(F_0(\lambda_{[p]}(x_{[p]})) + V^I)x_p^I + F_+(\lambda_{[p]}(x_{[p]}))x_p^U] \Big|_{\mathcal{X}=\mathcal{X}^*} \\ & = \partial_{x_{p'}^I} [F_+(\lambda_{[p]}(x))x_p^{U*}] \Big|_{x=x_{[p]}^*} = \partial_{\lambda} [F_+(\lambda)x_p^{U*}] \partial_{x^I} \lambda_{[p]}(x) \Big|_{x=x_{[p]}^*} \end{aligned} $$

(using again Remark 5 to introduce a dummy variable).

Last, the partial derivative $\partial_{x_{p'}^I} [(F_0(\lambda_{[p]}(x_{[p]})) + V^I)x_p^I + F_+(\lambda_{[p]}(x_{[p]}))x_p^U] \Big|_{\mathcal{X}=\mathcal{X}^*}$ is equal to $0_{n_I \times n_I}$ whenever $[p] \neq [p']$.

One shows with the same techniques that, for any $p, p' \in \mathcal{P}$,

$$ \partial_{x_{p'}^U} [(F_0(\lambda_{[p]}(x_{[p]})) + V^I)x_p^I + F_+(\lambda_{[p]}(x_{[p]}))x_p^U] \Big|_{\mathcal{X}=\mathcal{X}^*} = 0_{n_I \times n_U}. $$

Last, one has for any $p \in \mathcal{P}$,

$$ \partial_{x_p^U} (V^U x_p - F_-(\lambda_{[p]}(x_{[p]}))x_p^U) \Big|_{\mathcal{X}=\mathcal{X}^*} = V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}, $$

while for any $p, p' \in \mathcal{P}$ such that $p \neq p'$,

$$ \partial_{x_{p'}^U} (V^U x_p - F_-(\lambda_{[p]}(x_{[p]}))x_p^U) \Big|_{\mathcal{X}=\mathcal{X}^*} = 0_{n_U \times n_U}. $$

• With the help of the previous computations and adopting now the component ordering defined in (45), one may proceed to compute the three blocks $\mathcal{J}_{II}(\mathcal{X}^*)$, $\mathcal{J}_{IU}(\mathcal{X}^*)$ and $\mathcal{J}_{UU}(\mathcal{X}^*)$ in the decomposition (47). The matrix $\mathcal{J}_{II}(\mathcal{X}^*)$ is defined as

$$ \mathcal{J}_{II}(\mathcal{X}^*) = \mathrm{diag} \{ \mathcal{B}_{q_1}(\mathcal{X}^*); \dots; \mathcal{B}_{q_{|Q|}}(\mathcal{X}^*) \}, \quad (48a) $$

where, for any $q \in \mathcal{Q}$, the block $\mathcal{B}_q \in \mathbb{R}^{n_I|q| \times n_I|q|}$, is given by

$$ \mathcal{B}_q := (I_{|q|} \otimes V^I) + \begin{pmatrix} \partial_{x^I} [F_+(\lambda(x))x_{q,1}^{U*}] & (x_q^*) \\ \vdots & \\ \partial_{x^I} [F_+(\lambda(x))x_{q,|q|}^{U*}] & (x_q^*) \end{pmatrix} (\mathbf{1}_{|q|}^{\top} \otimes I_{n_I}). \quad (48b) $$

In coherence with the previous notations, in (48b) $x_{q,i}^{U*}$, $i = 1, \dots, |q|$, is the equilibrium value of the subpopulation numbered $i$ in the $q$ class/location; and $x_q^* = \sum \{x_p^* : p \in \mathcal{P}, [p] = q\} = \sum_{i=1}^{|q|} x_{q,i}^*$ is the total population present in that class. One also has:

$$ \mathcal{J}_{IU}(\mathcal{X}^*) = 0_{n_I | \mathcal{P} | n_U | \mathcal{P}|}, \quad (48c) $$

$$ \mathcal{J}_{UU}(\mathcal{X}^*) = I_{|\mathcal{P}|} \otimes V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix}, \quad (48d) $$

where $\otimes$ represents the Kronecker product. It is important to notice at this point that the 4th block $\mathcal{J}_{UI}(\mathcal{X}^*)$ has a block diagonal structure, with $|\mathcal{Q}|$ blocks, respectively of size $|q_1|, \dots, |q_{|Q|}|$.
Let us comment on that structure of the Jacobian matrix. First, formula (48a) shows a block-diagonal structure that corresponds to the class structure: there may be no influence whatsoever between two subpopulations whose classes (or locations) are different. Then (48b) unveils the interinfluence between the infected populations at the same location. The first term therein, namely $(I_{|q|} \otimes V^I)$, originates from the passage from a compartment to another within every given subpopulation $x_p$; while the second one — the product of two matrices of respective dimensions $|q|n_I \times n_I$ and $n_I \times |q|n_I$ — comes from the interinfluence between the different subpopulations in the same class $q$, which results from the fact that the force of infection $\lambda$ depends upon the total population $x_q = \sum_{p \in q} x_p$ at this location. The $n_I$ lines of this second term that correspond to a given subpopulation $p_{q,i}$, $i \in \{1, \dots, |q|\}$, are

$$ \partial_{x^I} [F_+(\lambda(x))x_{q,i}^{U*}] (x_q^*)(\mathbf{1}_{|q|}^{\top} \otimes I_{n_I}). $$

They depend only upon two quantities: on the one hand upon the equilibrium value $x_{q,i}^{U*}$ of this subpopulation; and on the other hand upon the global sum $x_q^{U*} - x_{q,i}^{U*}$ of the other subpopulations in the class $q$ (and not individually of every other subpopulation in the class $q$).

Formula (38) is the analogue of (48c) in the simple case of a unique subpopulation. It is a consequence of the fact that the equilibrium considered is disease-free. Last, formula (48d) testifies of the same type of diagonal structure than the first term in (48b), which gives in $\mathcal{J}_{II}(\mathcal{X}^*)$ a term $(I_{|\mathcal{P}|} \otimes V^I)$ similar to (48d).

### 5.1.3 DFE stability for a partitioned-population model

The following result treats the case of several subpopulations distributed in different locations in a permanent manner. This extension of Theorem 2 is essentially a rewriting, due to the fact that the partition of the subpopulations is unchanged along the time. It is put here for didactic reasons, as an intermediate towards the fully general case where the partition changes with respect to time, treated below in Theorem 4.

**Theorem 3** (Basic reproduction number and stability for epidemic models with fixed population partition). Let $\mathcal{P}$ be a finite set, $\mathcal{R}$ a partition of $\mathcal{P}$ and $\mathcal{Q} := \mathcal{P} \setminus \mathcal{R}$ the corresponding set of equivalence classes. Let $f$ and $\lambda_q, q \in \mathcal{Q}$, fulfilling Assumptions (B.1), (B.2), (B.3). Then the following properties hold.

* System (2)-(27) admits nonzero disease-free equilibrium points. The set of equilibrium points of this system is exactly the Cartesian product $\mathbf{X}_{\text{eq}}^{\mathcal{P}}$, where $\mathbf{X}_{\text{eq}}$, defined in (30), is the set of equilibrium points for an isolated subpopulation.

* For any nonzero disease-free equilibrium point $\mathcal{X}^* \in \mathbf{X}_{\text{eq}}^{\mathcal{P}}$ of system (2)-(27), let

$$ \mathcal{R}_{0,q} := \rho(F_q V^{-1}), \quad q \in \mathcal{Q}, \qquad (49a) $$

where $F_q, V \in \mathbb{R}^{n_I \times n_I}$ are given by:

$$ F_q := \partial_{x^I} (F_+(\lambda_q(x))x_q^{U*})|_{x=x_q^*}, \quad V := V^I. \qquad (49b) $$

Let

$$ \mathcal{R}_0 := \max_{q \in \mathcal{Q}} \mathcal{R}_{0,q}. \qquad (50) $$

The following assertions are true:

- if $\mathcal{R}_0 > 1$, then $\mathcal{X}^*$ is unstable;

- if $\mathcal{R}_0 < 1$, then $\mathcal{X}^*$ is marginally stable. Moreover, if the first identity in (33) holds with an equality, then $\mathcal{X}^*$ is locally asymptotically stable.

Recall that the notation $\partial_{x^I}$ is defined in (35).

While apparently complicated, the result enunciated in Theorem 3 is in fact quite natural. On the one hand, formula (50) states that the basic reproduction number of the system is the largest of the
basic reproduction numbers of the (sub-)systems describing the infection spread in each independent class (location) $q \in Q$. On the other hand, due to perfect mixing within every class, formula (49) states that this value may be computed for any class $q \in Q$ as the basic reproduction number corresponding to a unique population at the equilibrium point $x_q^* = \sum_{p \in q} x_p^*$. In this respect, it is important to notice that $F_q$ in (49b) is an $n_I \times n_I$ matrix, while the block $B_q$ corresponding to the class $q$ of the Jacobian matrix defined in (48b) is a $n_I|q| \times n_I|q|$ matrix.

Notice that when $Q$ contains a unique class gathering all subpopulations, then Theorem 3 reduces to the case, treated in Theorem 2, of a unique subpopulation. This is a consequence of the fact that what matters to portray the asymptotic evolution at the scale of a class $q \in Q$, is exclusively the total population $x_q$ present therein.

### Proof of Theorem 3.

* Due to (2), at equilibrium, every subpopulation $p \in P$ evolves uncoupled, according to

$$ \dot{x}_p = f(x_p, 0_{n_\lambda}). $$

It is therefore clear that the set of equilibrium points of (2)-(27) is exactly the Cartesian product $X_{eq}^{\mathcal{P}}$.

* Let us now consider the Jacobian matrix. As $\mathcal{J}_{IU}(\mathcal{X}^*)$ is null, one has

$$ s(\mathcal{J}(\mathcal{X}^*)) = \max\{s(\mathcal{J}_{II}(\mathcal{X}^*)), s(\mathcal{J}_{UU}(\mathcal{X}^*))\}, $$

and one is led to consider the two diagonal blocks $\mathcal{J}_{II}(\mathcal{X}^*)$ and $\mathcal{J}_{UU}(\mathcal{X}^*)$. We will exploit the block diagonal structure of these two matrices.

* Let us compute first the value of the dominant eigenvalues of the diagonal blocks $B_q$ of $\mathcal{J}_{II}(\mathcal{X}^*)$, $q \in Q$. For this, let us consider a subpopulation class $q$. We will assess the dominant eigenvalue of the corresponding matrix $B_q$ and show that it is equal to $\mathcal{R}_{0,q}$ defined in (49a). Let $\nu \in \mathbb{R}$ and the vector

$$ \begin{pmatrix} v_1 \\ \vdots \\ v_{|q|} \end{pmatrix} \in \mathbb{R}_+^{|q|n_I} $$

be respectively the dominant eigenvalue and a positive right-eigenvector of the matrix $B_q$ defined in (48b). In other words,

$$ \left( (I_{|q|} \otimes V^I) + \begin{pmatrix} \partial_{x^I} [F_+(\lambda_q(x))x_{q,1}^{U*}] & (x_q^*) \\ \vdots & \\ \partial_{x^I} [F_+(\lambda_q(x))x_{q,|q|}^{U*}] & (x_q^*) \end{pmatrix} (\mathbf{1}_{|q|}^\top \otimes I_{n_I}) \right) \begin{pmatrix} v_1 \\ \vdots \\ v_{|q|} \end{pmatrix} = \nu \begin{pmatrix} v_1 \\ \vdots \\ v_{|q|} \end{pmatrix}. \quad (51) $$

The map $x^U \mapsto F_+(\lambda)x^U$ is linear. Therefore, due to the fact that $x_q^* = x_{q,1}^* + \cdots + x_{q,|q|}^*$, so that $x_q^{U*} = x_{q,1}^{U*} + \cdots + x_{q,|q|}^{U*}$, one has

$$ (\mathbf{1}_{|q|}^{\top} \otimes I_{n_I}) \begin{pmatrix} \partial_{x^I} [F_+(\lambda_q(x))x_{q,1}^{U*}] & (x_q^*) \\ \vdots & \\ \partial_{x^I} [F_+(\lambda_q(x))x_{q,|q|}^{U*}] & (x_q^*) \end{pmatrix} = \sum_{i=1}^{|q|} \partial_{x^I} [F_+(\lambda_q(x))x_{q,i}^{U*}] (x_q^*) = \partial_{x^I} [F_+(\lambda_q(x))x_q^{U*}] (x_q^*). $$

The emergence of an epidemic in the class $q$ depends only upon the total population present therein and upon the local conditions of the infection, but not upon the division of this class in several subpopulations: the stability of the DFE is indeed a property of the class itself, independent of its constituting subpopulations. In coherence with this intuition, we now exhibit a dominant eigenvector at the level of the whole population. Using the previous computations, left-multiplying formula (51) by the matrix $(\mathbf{1}_{|q|}^{\top} \otimes I_{n_I}) \in \mathbb{R}^{n_I \times |q|n_I}$ yields the identity:

$$ \left( (\mathbf{1}_{|q|}^{\top} \otimes V^I) + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U*}] (x_q^*) (\mathbf{1}_{|q|}^{\top} \otimes I_{n_I}) \right) \begin{pmatrix} v_1 \\ \vdots \\ v_{|q|} \end{pmatrix} = \nu (\mathbf{1}_{|q|}^{\top} \otimes I_{n_I}) \begin{pmatrix} v_1 \\ \vdots \\ v_{|q|} \end{pmatrix}, $$
that is, for any $q \in \mathcal{Q}$,

$$ (V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U*}] (x_q^*) ) v_q = \nu v_q, \quad \text{with } v_q := v_1 + \dots + v_{|q|}. $$

The vector $v_q$ is a positive eigenvector, it is therefore associated to the dominant eigenvalue $\nu$. The matrix involved is the same than the one appearing in Theorem 2, and the same handling yields the formulas in (49), similar to (36). The basic reproduction number of the system then appears as the largest of the dominant eigenvalues of the diagonal blocks $B_q$ of $\mathcal{J}_{II}(\mathcal{X}^*)$, $q \in \mathcal{Q}$.

• Let us now consider the matrix $\mathcal{J}_{UU}(\mathcal{X}^*)$, expressed in (48d). Each of the $|\mathcal{P}|$ blocks possesses a unique eigenvector associated to the eigenvalue 0, therefore 0 is an eigenvalue of $\mathcal{J}_{UU}(\mathcal{X}^*)$ whose algebraic and geometrical multiplicities are both equal to $|\mathcal{P}|$. The evolution in each class of $\mathcal{Q}$ is independent of what happens in the other ones, and conducting the same analysis than in the proof of Theorem 2 gives the same stability result, depending on whether $\mathcal{R}_0$ given in (50) is smaller or larger than 1. This achieves the proof of Theorem 3. $\square$

## 5.2 Basic reproduction number of general epidemiological models in commuting populations

We study now the behaviour of the solutions of an epidemiological model in commuting population given by (12)-(27), for a T-periodic admissible switching signal $\mathcal{R}(\cdot)$ with interswitch durations in a given set $\Lambda_m$ defined in (11), assuming Assumptions (B.1), (B.2), (B.3) fulfilled.

We first show in Section 5.2.1 that the disease-free periodic solutions are indeed equilibrium points. We then extend in Section 5.2.2 the framework introduced in Section 5.1.1 for representing the state variable, and state the main result in Section 5.2.3, together with an important corollary. Section 5.2.4 provides a method to compute more easily the quantities involved. Finally a simple example is displayed and analyzed in Section 5.2.5, in order to make visible the different steps.

### 5.2.1 Equilibrium states of a periodic commuting population model

We are interested in modelling and analysing the effects of periodic ‘commutations’. The latter are specific variations in time of the considered system. Having a more acute look at the formulas (12) with T-periodic switching signal, one sees that the time-variation only affects the definition of the (piecewise constant) force of infection, namely the subpopulation mixing on each time interval.

Consider a disease-free evolution of system (12)-(27), that is one for which $x^I(t) = 0_{n_I}$ and $\lambda_{[p]_{\mathcal{R}(t)}}(t) = 0_{n_\lambda}$ for any $t \in \mathcal{I}$. Seen from the subpopulation $p \in \mathcal{P}$, this implies (see formula (34))

$$ \dot{x}_p = f(x_p(t), 0_{n_\lambda}) = \begin{pmatrix} 0_{n_I} \\ V^U x_p \end{pmatrix}, \quad t \in \mathcal{I}. $$

The only changes considered here are the commutations. The latter modify the mixing of the subpopulations, and possibly the value of the components of the force of infection. However, the structure of the model is not touched, and in particular the non-infective passages from a compartment to another are not modified along time. In consequence, contrary e.g. to what happens with the settings developed to analyse seasonal phenomena (see for example the papers [6, 9, 45]), the disease-free attractors are all the equilibrium points in the space $\mathbb{R}_+^{|\mathcal{P}|n}$. See more details in the first part of Theorem 4.

### 5.2.2 State variable and Jacobian matrix representation for commuting population model

We extend here the notational setting introduced in Section 5.1.1, in order to accommodate the changes of partition from a maximal subinterval of $\mathcal{I}$ to the next one that model the commutations. Denote $\mathcal{R}(\mathcal{I}_k)$ the partition in force during any maximal subinterval of $\mathcal{I}_k$. The class of subpopulation $p$ during this period of time will be denoted accordingly $[p]_{\mathcal{R}(\mathcal{I}_k)}$.

For simplicity we assume the existence of a ‘reference ordering of the subpopulations’, pre-existing the infection; and, for any union of intervals $\mathcal{I}_k$, $k = 1, \dots, m$, of another lexical ordering such as the ones presented in Section 5.1.1, having the property of attributing contiguous ranks to the subpopulations located in the same class. Denote $\mathcal{M}_k \in \mathbb{R}^{|\mathcal{P}| \times |\mathcal{P}|}$ the permutation matrix allowing to move from the reference ordering to the ordering attached to $\mathcal{I}_k$. In other words, the subpopulation with rank $i$ in the
initial ordering, is ranked $i'$ in the new one, such that $M_k e_i = e_{i'}$ during $\mathcal{I}_k$, where $e_i, e_{i'}$ are respectively the $i$-th and $i'$-th vectors of the canonical basis of $\mathbb{R}^{|P|}$.

Consider a given DFE $\mathcal{X}^* \in \mathbb{R}_+^{n|P|}$. The Jacobian matrix of system (12)-(27) is constant on any interval in $\mathcal{I}_k$. Referring to formula (47), its value when expressed according to the subpopulation ordering adapted to the partition $\mathcal{R}(\mathcal{I}_k)$, will be denoted $\mathcal{J}_k(\mathcal{X}^*)$, and the corresponding nonzero blocks $\mathcal{J}_{II,k}(\mathcal{X}^*)$, $\mathcal{J}_{UI,k}(\mathcal{X}^*)$, $\mathcal{J}_{UU,k}(\mathcal{X}^*)$. In the reference ordering, the Jacobian matrix $\hat{\mathcal{J}}_k(\mathcal{X}^*)$ is then equal to

$$ \hat{\mathcal{J}}_k(\mathcal{X}^*) = \begin{pmatrix} \hat{\mathcal{J}}_{II,k}(\mathcal{X}^*) & \hat{\mathcal{J}}_{IU,k}(\mathcal{X}^*) \\ \hat{\mathcal{J}}_{UI,k}(\mathcal{X}^*) & \hat{\mathcal{J}}_{UU,k}(\mathcal{X}^*) \end{pmatrix}, \qquad (52a) $$

$$ \hat{\mathcal{J}}_{\eta\eta',k}(\mathcal{X}^*) := (\mathcal{M}_k^\top \otimes I_{n_\eta}) \mathcal{J}_{\eta\eta',k}(\mathcal{X}^*)(\mathcal{M}_k \otimes I_{n_{\eta'}}), \quad \eta, \eta' \in \{I, U\}. \qquad (52b) $$

In particular, one has $\hat{\mathcal{J}}_{IU,k}(\mathcal{X}^*) = \mathcal{J}_{IU,k}(\mathcal{X}^*) = 0_{n_I|P| \times n_U|P|}$ at any equilibrium point $\mathcal{X}^*$.

### 5.2.3 DFE stability for a commuting population model

With the previous notations, it is now possible to state the main result of the paper.

**Theorem 4** (Basic reproduction number and stability for epidemic models in commuting populations). Let $\mathcal{P}$ be a finite set, $\mathcal{R}(\cdot)$ a T-periodic admissible switching signal with set of switching times $\mathbb{R}_+ \setminus \mathcal{I}$ according to pattern (11), and $\mathcal{Q}(t) := \mathcal{P} \setminus \mathcal{R}(t)$ the corresponding set of equivalence classes at time $t \in \mathcal{I}$. Let $f$ be given as in (27) and $\lambda_q, q \in \mathcal{Q}(\mathcal{I}_k)$, $k = 1, \dots, m$, fulfilling Assumptions (B.1), (B.2), (B.3). Then the following properties hold.

* The disease-free T-periodic solutions of system (12)-(27) are exactly the equilibrium points of the Cartesian product $\mathbf{X}_{\text{eq}}^{|P|}$, where $\mathbf{X}_{\text{eq}}$, defined in (30), is the set of equilibrium points for an isolated subpopulation.

* For any nonzero disease-free equilibrium point $\mathcal{X}^* \in \mathbf{X}_{\text{eq}}^{|P|}$ of system (12)-(27), let

$$ \mathcal{R}_0 := \rho \left( e^{\hat{\mathcal{J}}_{II,m}(\mathcal{X}^*)(\tau_m - \tau_{m-1})} \dots e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(\tau_2 - \tau_1)} e^{\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)(\tau_1 - \tau_0)} \right), \qquad (53) $$

for the $|P|n_I \times |P|n_I$ Jacobian matrices $\hat{\mathcal{J}}_{II,k}(\mathcal{X}^*)$, $k = 1, \dots, m$, defined in (52). The following assertions are true:

- if $\mathcal{R}_0 > 1$, then $\mathcal{X}^*$ is unstable;

- if $\mathcal{R}_0 < 1$, then $\mathcal{X}^*$ is marginally stable. Moreover, if the first identity in (33) holds with an equality, then $\mathcal{X}^*$ is locally asymptotically stable.

Recall that the notion of admissible switching signal and the notation $\mathcal{I}$ have been defined in Definition 2. Computing exactly $\mathcal{R}_0$ involves computing the spectral radius of a matrix of size $|P|n_I \times |P|n_I$.

**Remark 6.** The computation of $\mathcal{R}_0$ by (53) is made easier in the case where, in spite of the commutations, the infection evolves independently in completely disjoint subgroups of the population. This is indeed the case if for some subpopulations $p, p' \in \mathcal{P}$, there is no positive integer $l$ and finite sequence $p_0, p_1, \dots, p_l$ of elements of $\mathcal{P}$ such that $p_0 := p, p_l := p'$ and for any $i = 0, \dots, l-1$, there exists $k \in \{1, \dots, m\}$ such that

$$ [p_i]_{\mathcal{R}(\mathcal{I}_k)} = [p_{i+1}]_{\mathcal{R}(\mathcal{I}_k)}. $$

(Alternatively, this means that the transitive closure of the union of the relations $\mathcal{R}(\mathcal{I}_k)$, $k = 1, \dots, m$, has at least two distinct classes of equivalence.) When such disconnection occurs, the infection initially present in the subpopulation $p$ cannot move to the subpopulation $p'$ and vice versa, and the matrix product in (53) involves several diagonal blocks. As is the case for Theorem 3, the value of $\mathcal{R}_0$ then appears as the maximal value of the corresponding quantities computed in each connected component of $\mathcal{P}$ (that is, on each class of equivalence of the transitive closure of the union of the relations $\mathcal{R}(\mathcal{I}_k)$, $k = 1, \dots, m$).
### Proof of Theorem 4.

• Let us first consider a disease-free trajectory for (12)-(27). As observed previously, the coupling between equations describing the evolution of different subpopulations is manifested only through the infection process: due to (12), every subpopulation $p \in \mathcal{P}$ evolves according to

$$ \dot{x}_p = f(x_p, 0_{n_\lambda}). $$

It is therefore clear that the set of equilibrium points of (12)-(27) is exactly the Cartesian product $\mathbf{X}_{\text{eq}}^{\mathcal{P}}$.

More precisely, for any disease-free periodic trajectory, any subpopulation $p \in \mathcal{P}$ undergoes a stationary evolution, governed (see (27c)) by the evolution

$$ \dot{x}_p^U = V^U \begin{pmatrix} 0_{n_I \times n_U} \\ I_{n_U} \end{pmatrix} x_p^U. $$

By Assumption (B.2), the matrix present in this formula has 0 as simple eigenvalue, and all other eigenvalues have negative real parts. Therefore, every trajectory converges towards an equilibrium point, and no other periodic solution exists.

• As established in Section 5.2.2, the Jacobian matrix of the system is constant on every maximal subinterval of $\mathcal{I}$, and equal to $\tilde{J}_k(\mathcal{X}^*)$ on any $\mathcal{I}_k$, $k = 1, \dots, m$. Using the block-diagonal structure and arguments identical to the demonstration of Theorem 3, one shows the marginal stability of the uninfected components; while the stability of the infected components depends upon the position with respect to 1 of the spectral radius of the monodromy matrix. The latter is exactly $\mathcal{R}_0$. This achieves the proof of Theorem 4. □

We now present an important corollary of Theorem 4. The following result shows that, when the conditions of transmission are identical in the whole population (in a sense made precise in the statement), then the disease evolution is identical to that of a perfectly mixed population gathering all the subpopulations.

**Corollary 5 (Homogeneous density-independent transmission conditions).** Under the conditions of Theorem 4, let $\mathcal{X}^* \in \mathbf{X}_{\text{eq}}^{\mathcal{P}}$ be a nonzero disease-free equilibrium point of system (12)-(27) such that

$$ \forall k \in \{1, \dots, m\}, \quad \forall q, q' \in \mathcal{Q}(\mathcal{I}_k), \quad \partial_{x^I} [F_+(\lambda_q(x))x_q^{U*}] (x_q^*) = \partial_{x^I} [F_+(\lambda_{q'}(x))x_{q'}^{U*}] (x_{q'}^*) =: G_k. \quad (54) $$

Then, the basic reproduction number given in (53) is equal to

$$ \mathcal{R}_0 = \rho \left( e^{(V^I+G_m)(\tau_m-\tau_{m-1})} \dots e^{(V^I+G_1)(\tau_1-\tau_0)} \right). \quad (55) $$

A major interest of Corollary 5, when it applies, is of course to reduce the determination of the basic offspring number to the computation of the spectral radius of a square matrix of size $n_I$.

To apply this result, the conditions of transmission of the infection should be identical in all classes and at all times, in the sense precisely defined by (54). The function $F_+$ being the same for the whole system (independent of the class and time), this condition relates only to the forces of infection $\lambda_q$. Except in very special conditions where the sizes of the classes match exactly, one has to assume that, for any $k \in \{1, \dots, m\}$, $q \in \mathcal{Q}(\mathcal{I}_k)$, the function

$$ x^* \mapsto \partial_{x^I} [F_+(\lambda_q(x))x^{U*}] (x^*) $$

is constant on the cone $\mathbf{X}_{\text{eq}}$. In particular, due to the linearity of the first factor with respect to $x^{U*}$, one should have, for any scalar $\alpha > 0$,

$$ \frac{\partial [F_+(\lambda)\alpha x^{U*}]}{\partial \lambda} \frac{\partial \lambda_q}{\partial x^I}(\alpha x^*) = \frac{\partial [F_+(\lambda)x^{U*}]}{\partial \lambda} \frac{\partial \lambda_q}{\partial x^I}(x^*) = \frac{\partial [F_+(\lambda)\alpha x^{U*}]}{\partial \lambda} \frac{1}{\alpha} \frac{\partial \lambda_q}{\partial x^I}(x^*) $$

A natural condition is therefore that the partial derivative with respect to $x^I$ of the functions $\lambda_q$ defining the forces of infection, are homogeneous of degree $-1$. It is sufficient for this that the functions $\lambda_q$ themselves are homogeneous of degree 0, a characteristic of frequency-dependent transmission.

One shows without difficulty that (54) is verified e.g. when the functions $\lambda_q$ are all identical in every class $q \in \mathcal{Q}(\mathcal{I}_k)$, and constant on the cone $\mathbf{X}_{\text{eq}}$. (This constant function is homogeneous of degree 0.)
Thus, Corollary 5 establishes in particular that, if the conditions of the transmission are identical in all locations (classes), then the value of the basic reproduction number is identical whatever the mixing conditions (that is, whatever the commutation pattern and the number and composition of the classes between two successive switchings). Up to our knowledge, this property of invariance relatively to the mixing is new. See an application to the SIR model in Section 5.2.5.

As a last remark, notice that when Corollary 5 holds, the evolution of the infected compartments of any subpopulation $p \in \mathcal{P}$ in the vicinity of $x_{I,p}^*$ is identical on $\mathcal{I}_k$ to that of the uncoupled autonomous linear equation $\dot{x}_{I,p} = (V^I + G_k)(x_{I,p} - x_{I,p}^*)$.

*Proof of Corollary 5.* Let $v \in \mathbb{R}^{n_I}$, and $k \in \{1, \dots, m\}$. For any $q \in \mathcal{Q}(\mathcal{I}_k)$, one has

$$ (\mathbf{1}_{|q|} \otimes v)^{\mathrm{T}} \mathcal{B}_q = (\mathbf{1}_{|q|}^{\mathrm{T}} \otimes v^{\mathrm{T}} (V^I + F_q)) $$

for $\mathcal{B}_q$ in (48b).

Then, left-multiplying the matrix $e^{\mathcal{J}_{II,k}(\mathcal{X}^*)(\tau_k - \tau_{k-1})}$ by the vector $(\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}}$ yields

$$ (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} e^{\mathcal{J}_{II,k}(\mathcal{X}^*)(\tau_k - \tau_{k-1})} = ((\mathbf{1}_{|q_1|}^{\mathrm{T}} \otimes v^{\mathrm{T}} e^{(V^I + F_{q_1})(\tau_k - \tau_{k-1})}) \dots (\mathbf{1}_{|q_{|\mathcal{Q}(\mathcal{I}_k)|}|}^{\mathrm{T}} \otimes v^{\mathrm{T}} e^{(V^I + F_{q_{|\mathcal{Q}(\mathcal{I}_k)|}})(\tau_k - \tau_{k-1})})). $$

We now use assumption (54), which states that $F_{q_1} = \dots = F_{q_{|\mathcal{Q}(\mathcal{I}_k)|}} = G_k$, to deduce

$$ (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} e^{\mathcal{J}_{II,k}(\mathcal{X}^*)(\tau_k - \tau_{k-1})} = (\mathbf{1}_{|\mathcal{P}|}^{\mathrm{T}} \otimes v^{\mathrm{T}} e^{(V^I + G_k)(\tau_k - \tau_{k-1})}), $$

as the number of elements in all classes is equal to the total number of subpopulations, see (46). One thus obtains

$$ (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} e^{\mathcal{J}_{II,k}(\mathcal{X}^*)(\tau_k - \tau_{k-1})} = (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} (I_{|\mathcal{P}|} \otimes e^{(V^I + G_k)(\tau_k - \tau_{k-1})}). \quad (56) $$

From the previous formula, one deduces the same property for $e^{\hat{\mathcal{J}}_{II,k}(\mathcal{X}^*)}$. In effect, (52b) yields

$$ \begin{align*} (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} e^{\hat{\mathcal{J}}_{II,k}(\mathcal{X}^*)(\tau_k - \tau_{k-1})} &= (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} (\mathcal{M}_k^{\mathrm{T}} \otimes I_{n_I}) e^{\mathcal{J}_{II,k}(\mathcal{X}^*)(\tau_k - \tau_{k-1})} (\mathcal{M}_k \otimes I_{n_I}) \\ &= (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} e^{\mathcal{J}_{II,k}(\mathcal{X}^*)(\tau_k - \tau_{k-1})} (\mathcal{M}_k \otimes I_{n_I}) \\ &\quad \text{(because } \mathbf{1}_{|\mathcal{P}|}^{\mathrm{T}} \mathcal{M}_k = \mathbf{1}_{|\mathcal{P}|}^{\mathrm{T}}, \text{ as } \mathcal{M}_k \text{ is a permutation matrix}) \\ &= (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} (I_{|\mathcal{P}|} \otimes e^{(V^I + G_k)(\tau_k - \tau_{k-1})}) (\mathcal{M}_k \otimes I_{n_I}) \\ &\quad \text{(thanks to (56))} \\ &= (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} (\mathcal{M}_k \otimes e^{(V^I + G_k)(\tau_k - \tau_{k-1})}) \\ &= (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} (I_{|\mathcal{P}|} \otimes e^{(V^I + G_k)(\tau_k - \tau_{k-1})}) \\ &\quad \text{(because } \mathbf{1}_{|\mathcal{P}|}^{\mathrm{T}} \mathcal{M}_k = \mathbf{1}_{|\mathcal{P}|}^{\mathrm{T}}) \\ &= (\mathbf{1}_{|\mathcal{P}|}^{\mathrm{T}} \otimes v^{\mathrm{T}} e^{(V^I + G_k)(\tau_k - \tau_{k-1})}) \\ &= (\mathbf{1}_{|\mathcal{P}|} \otimes e^{(V^I + G_k)^{\mathrm{T}}(\tau_k - \tau_{k-1})} v)^{\mathrm{T}}. \end{align*} $$

In order to finish the proof, one applies this property recursively for $k = m, \dots, 1$:

$$ \begin{align*} (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} e^{\hat{\mathcal{J}}_{II,m}(\mathcal{X}^*)(\tau_m - \tau_{m-1})} \dots e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(\tau_2 - \tau_1)} e^{\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)(\tau_1 - \tau_0)} &= (\mathbf{1}_{|\mathcal{P}|} \otimes e^{(V^I + G_m)^{\mathrm{T}}(\tau_m - \tau_{m-1})} v)^{\mathrm{T}} e^{\hat{\mathcal{J}}_{II,m-1}(\mathcal{X}^*)(\tau_{m-1} - \tau_{m-2})} \dots e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(\tau_2 - \tau_1)} e^{\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)(\tau_1 - \tau_0)} \\ &= (\mathbf{1}_{|\mathcal{P}|} \otimes e^{(V^I + G_1)^{\mathrm{T}}(\tau_1 - \tau_0)} \dots e^{(V^I + G_m)^{\mathrm{T}}(\tau_m - \tau_{m-1})} v)^{\mathrm{T}} \\ &= (\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}} (I_{|\mathcal{P}|} \otimes e^{(V^I + G_m)(\tau_m - \tau_{m-1})} \dots e^{(V^I + G_1)(\tau_1 - \tau_0)}). \end{align*} $$

Choose now for $v \in \mathbb{R}^{n_I}$ the positive eigenvector corresponding to the dominant eigenvalue of the positive matrix $e^{(V^I + G_m)(\tau_m - \tau_{m-1})} \dots e^{(V^I + G_1)(\tau_1 - \tau_0)}$. The vector $(\mathbf{1}_{|\mathcal{P}|} \otimes v)^{\mathrm{T}}$ is then a left-eigenvector of the matrix $e^{\hat{\mathcal{J}}_{II,m}(\mathcal{X}^*)(\tau_m - \tau_{m-1})} \dots e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(\tau_2 - \tau_1)} e^{\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)(\tau_1 - \tau_0)}$. Being positive, it corresponds to its dominant eigenvalue. This shows (55) and achieves the demonstration of Corollary 5. $\square$
### 5.2.4 Computing the exponentials of the blocks $\mathcal{B}_q$ in the Jacobian matrix $\mathcal{J}_{II,k}(\mathcal{X}^*)$

The following result is useful to reduce the dimension of the matrix product that appears in the expression of $\mathcal{R}_0$ in (53).

**Lemma 1.** For $\mathcal{B}_q$ defined in (48b), one has for any $t \ge 0$,

$$e^{\mathcal{B}_q t} = (I_{|q|} \otimes e^{V^I t}) + \int_0^t \left( e^{V^I(t-s)} \partial_{x^I} [F_+(\lambda_q(x))x_{q,1}^{U^*}]|_{x=x_q^*} \right. \\ \left. \qquad \vdots \qquad e^{V^I(t-s)} \partial_{x^I} [F_+(\lambda_q(x))x_{q,|q|}^{U^*}]|_{x=x_q^*} \right) e^{(V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}])|_{x=x_q^*} s} ds \ (\mathbf{1}_{|q|} \otimes I_{n_I})^\top. \quad (57)$$

In particular, for any infection model (27) presenting a unique infected compartment ($n_I = 1$), one has

$$e^{\mathcal{B}_q t} = e^{V^I t} I_{|q|} + e^{V^I t} \left( e^{\partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}]}|_{x=x_q^*} t - 1 \right) \begin{pmatrix} \frac{\partial_{x^I} [F_+(\lambda_q(x))x_{q,1}^{U^*}]}{\partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}]}|_{x=x_q^*} \\ \vdots \\ \frac{\partial_{x^I} [F_+(\lambda_q(x))x_{q,|q|}^{U^*}]}{\partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}]}|_{x=x_q^*} \end{pmatrix} \mathbf{1}_{|q|}^\top. \quad (58)$$

Using (57) to compute $e^{\mathcal{B}_q t}$ amounts essentially to compute the $|q| + 1$ matrices of size $n_I \times n_I$

$$e^{V^I t} \quad \text{and} \quad \int_0^t e^{-V^I s} \partial_{x^I} [F_+(\lambda_q(x))x_{q,i}^{U^*}]|_{x=x_q^*} e^{(V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}])|_{x=x_q^*} s} ds, \quad i = 1, \dots, |q|.$$

Notice also that formula (58), which holds true when $n_I = 1$, may be generalised for $n_I > 1$ when the matrix $V^I$ commutes with every $\partial_{x^I} [F_+(\lambda_q(x))x_{q,i}^{U^*}]|_{x=x_q^*}$, $i = 1, \dots, |q|$.

*Proof of Lemma 1.* By definition, $e^{\mathcal{B}_q t}$ is characterised by the fact that: for any $y(0) \in \mathbb{R}^{n_I|q|}$,

$$y(t) = e^{\mathcal{B}_q t} y(0), \quad t \ge 0 \quad \Leftrightarrow \quad \dot{y} = \mathcal{B}_q y, \quad t \ge 0.$$

Let $\bar{y} := (\mathbf{1}_{|q|} \otimes I_{n_I})^\top y = y_1 + \dots + y_q$. From the structure of $\mathcal{B}_q$ in (48b) and the fact that $x_{q,1}^{U^*} + \dots + x_{q,|q|}^{U^*} = x_q^{U^*}$, one deduces, for any $i = 1, \dots, |q|$,

$$\dot{\bar{y}} = (V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}]|_{x=x_q^*}) \bar{y}, \quad \dot{y}_i = V^I y_i + \partial_{x^I} [F_+(\lambda_q(x))x_{q,i}^{U^*}]|_{x=x_q^*} \bar{y},$$

which writes as the following triangular system

$$\begin{pmatrix} \dot{\bar{y}} \\ \dot{y}_i \end{pmatrix} = \begin{pmatrix} V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}]|_{x=x_q^*} & 0_{n_I \times n_I} \\ \partial_{x^I} [F_+(\lambda_q(x))x_{q,i}^{U^*}]|_{x=x_q^*} & V^I \end{pmatrix} \begin{pmatrix} \bar{y} \\ y_i \end{pmatrix}.$$

By integration of the first equation one gets

$$\bar{y}(t) = e^{(V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}])|_{x=x_q^*} t} \bar{y}(0) = e^{(V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}])|_{x=x_q^*} t} (\mathbf{1}_{|q|} \otimes I_{n_I})^\top y(0);$$

and the second equation thus yields, for any $i = 1, \dots, |q|$,

$$y_i(t) = e^{V^I t} y_i(0) + \int_0^t e^{V^I(t-s)} \partial_{x^I} [F_+(\lambda_q(x))x_{q,i}^{U^*}]|_{x=x_q^*} e^{(V^I + \partial_{x^I} [F_+(\lambda_q(x))x_q^{U^*}])|_{x=x_q^*} s} ds \ (\mathbf{1}_{|q|} \otimes I_{n_I})^\top y(0).$$
This yields

$$
y(t) = (I_{|q|} \otimes e^{V^I t}) y(0) \\
+ \int_0^t (I_{|q|} \otimes e^{V^I (t-s)}) \left( \begin{array}{c} \partial_{x^I} [F_+(\lambda_q(x))x_{q,1}^{U*}] \\ \vdots \\ \partial_{x^I} [F_+(\lambda_q(x))x_{q,|q|}^{U*}] \end{array} \right) \Bigg|_{x=x_q^*} e^{(V^I + \partial_{x^I} [F_+(\lambda_q(x))x_{q}^{U*}])|_{x=x_q^*}} s ds (\mathbf{1}_{|q|} \otimes I_{n_I})^\top y(0),
$$

and finally, as $y(t) = e^{\mathcal{B}_q t} y(0)$, one identifies

$$
e^{\mathcal{B}_q t} = (I_{|q|} \otimes e^{V^I t}) \\
+ \int_0^t (I_{|q|} \otimes e^{V^I (t-s)}) \left( \begin{array}{c} \partial_{x^I} [F_+(\lambda_q(x))x_{q,1}^{U*}] \\ \vdots \\ \partial_{x^I} [F_+(\lambda_q(x))x_{q,|q|}^{U*}] \end{array} \right) \Bigg|_{x=x_q^*} e^{(V^I + \partial_{x^I} [F_+(\lambda_q(x))x_{q}^{U*}])|_{x=x_q^*}} s ds (\mathbf{1}_{|q|} \otimes I_{n_I})^\top.
$$

which yields (57).

Deduction of formula (58) in case where $n_I = 1$ and all matrix components are then scalar, is then straightforward. $\square$

5.2.5 An illustrative example

• In order to illustrate the previous result, we consider the simple example shown in Figures 1 and 2. It consists in five subpopulations, $\mathcal{P} = \{p_1, p_2, p_3, p_4, p_5\}$, with $T$-periodic commutations between two different modes. The reference ordering mentioned in Section 5.2.2 is for example the order of the enumeration of the elements of $\mathcal{P}$ in Figure 1. To fix the ideas, we put $\mathcal{I} := \mathbb{R}_+ \setminus \mathbb{N}$, and $\mathcal{I} = \mathcal{I}_1 \cup \mathcal{I}_2$, with

$$
\mathcal{I}_1 := (0, \tau) \mod T, \quad \mathcal{I}_2 := (\tau, T) \mod T,
$$

for some $0 < \tau < T$. This is a switching pattern from the set $\Lambda_2$ defined in (11). We assume that the corresponding partitions are given by the two partitions $\mathcal{R}_1, \mathcal{R}_2$ presented in Figure 2, that is

$$
\mathcal{R}(\mathcal{I}_1) = \mathcal{R}_1 := \{\{p_1, p_2, p_3\}, \{p_4, p_5\}\}, \quad \mathcal{R}(\mathcal{I}_2) = \mathcal{R}_2 := \{\{p_1\}, \{p_2, p_5\}, \{p_3, p_4\}\}, \quad (59)
$$

in such a way that $m = 2$, $|Q(\mathcal{I}_1)| = 2$, $|Q(\mathcal{I}_2)| = 3$ for the quotient sets $Q(\mathcal{I}_1) = Q_1 = \mathcal{P} \setminus \mathcal{R}_1$, $Q(\mathcal{I}_2) = Q_2 = \mathcal{P} \setminus \mathcal{R}_2$.

Choosing for example to order the classes of subpopulations, and the subpopulations within each class, in the order in which they appear in (59), yields the correspondence

$$
p_1 = p_{1,1}, p_2 = p_{1,2}, p_3 = p_{1,3}, p_4 = p_{2,1}, p_5 = p_{2,2}
$$

during $\mathcal{I}_1$, and

$$
p_1 = p_{1,1}, p_2 = p_{2,1}, p_3 = p_{3,1}, p_4 = p_{3,2}, p_5 = p_{2,2}
$$

during $\mathcal{I}_2$. With this choice, the permutation matrices $\mathcal{M}_k$, $k = 1, 2$, defined in Section 5.2.2 are:

$$
\mathcal{M}_1 := I_5, \quad \mathcal{M}_2 := \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{pmatrix}. \qquad (60)
$$

Theorem 4 provides the exact value of $\mathcal{R}_0$, namely

$$
\mathcal{R}_0 = \rho \left( e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(T-\tau)} e^{\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)\tau} \right) = \rho \left( \mathcal{M}_2^\top e^{\mathcal{J}_{II,2}(\mathcal{X}^*)(T-\tau)} \mathcal{M}_2 e^{\mathcal{J}_{II,1}(\mathcal{X}^*)\tau} \right). \quad (61)
$$

When $\mathcal{R}_0 < 1$, Theorem 4 yields local asymptotic stability of the equilibrium point $\mathcal{X}^*$, due to the fact that the first identity in (33) holds here with an equality, see Remark 3. When $\mathcal{R}_0 > 1$, Theorem 4 yields instability.
• We consider e.g. the SIR model (9), for which $n = 3$, $n_I = 1$ and $n_\lambda = 1$. The global evolution is governed by equation (12), that is here

$$ \dot{I}_i = \beta_1(\mathcal{I}_1) \frac{I_1 + I_2 + I_3}{N_1 + N_2 + N_3} S_i - (\gamma + \mu)I_i, \quad i = 1, 2, 3 \qquad (62a) $$

$$ \dot{S}_i = -\beta_1(\mathcal{I}_1) \frac{I_1 + I_2 + I_3}{N_1 + N_2 + N_3} S_i + \mu(I_i + R_i), \quad i = 1, 2, 3, \qquad (62b) $$

$$ \dot{I}_i = \beta_2(\mathcal{I}_1) \frac{I_4 + I_5}{N_4 + N_5} S_i - (\gamma + \mu)I_i, \quad i = 4, 5, \qquad (62c) $$

$$ \dot{S}_i = -\beta_2(\mathcal{I}_1) \frac{I_4 + I_5}{N_4 + N_5} S_i + \mu(I_i + R_i), \quad i = 4, 5, \qquad (62d) $$

$$ \dot{R}_i = \gamma I_i - \mu R_i, \quad i = 1, 2, 3, 4, 5, \qquad (62e) $$

on $\mathcal{I}_1$, and

$$ \dot{I}_1 = \beta_1(\mathcal{I}_2) \frac{I_1}{N_1} S_1 - (\gamma + \mu)I_1, \qquad (63a) $$

$$ \dot{S}_1 = -\beta_1(\mathcal{I}_2) \frac{I_1}{N_1} S_1 + \mu(I_1 + R_1), \qquad (63b) $$

$$ \dot{I}_i = \beta_2(\mathcal{I}_2) \frac{I_2 + I_5}{N_2 + N_5} S_i - (\gamma + \mu)I_i, \quad i = 2, 5, \qquad (63c) $$

$$ \dot{S}_i = -\beta_2(\mathcal{I}_2) \frac{I_2 + I_5}{N_2 + N_5} S_i + \mu(I_i + R_i), \quad i = 2, 5, \qquad (63d) $$

$$ \dot{I}_i = \beta_3(\mathcal{I}_2) \frac{I_3 + I_4}{N_3 + N_4} S_i - (\gamma + \mu)I_i, \quad i = 3, 4, \qquad (63e) $$

$$ \dot{S}_i = -\beta_3(\mathcal{I}_2) \frac{I_3 + I_4}{N_3 + N_4} S_i + \mu(I_i + R_i), \quad i = 3, 4, \qquad (63f) $$

$$ \dot{R}_i = \gamma I_i - \mu R_i, \quad i = 1, 2, 3, 4, 5, \qquad (63g) $$

on $\mathcal{I}_2$.

The decomposition of the right-hand side $f$ for the SIR model (9) has been given in the formulas (18). In particular,

$$ V^I = -(\gamma + \mu), \quad \lambda(x) = \beta \frac{I}{S + I + R}, \quad F_+(\lambda) = \lambda \begin{pmatrix} 1 & 0 \end{pmatrix}, $$

and for any $x^* \in \mathbf{X}_{\text{eq}}$,

$$ \frac{\partial [F_+(\lambda)x^{U*}]}{\partial \lambda} \frac{\partial \lambda}{\partial x^I}(x^*) = \begin{pmatrix} 1 & 0 \end{pmatrix} x^{U*} \frac{\beta}{\mathbf{1}_{n_V}^\top x^{U*}} = \beta \frac{S^*}{S^*} = \beta. \qquad (64) $$

• With these computations, let us now apply the preceding results to analyze the behaviour of this system, in the vicinity of a given equilibrium point $\mathcal{X}^*$. The latter is uniquely defined by the numbers of susceptible individuals $S_i^*, i = 1, \dots, 5$, in the five subpopulations. At equilibrium, all individuals are susceptible and one has necessarily $S_i^* = |p_i|$. For the subpopulations chosen in Figure 1 this means

$$ S_1^* = 6, \quad S_2^* = 1, \quad S_3^* = 5, \quad S_4^* = 3, \quad S_5^* = 7. \qquad (65) $$

In the general case, applying (61) necessitates first to assess the value of the exponential of the two matrices $\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)\tau$ and $\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(T - \tau)$. The first matrix may be computed from the fact that $\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*) + (\gamma + \mu)I_5 = \mathcal{J}_{II,1}(\mathcal{X}^*) + (\gamma + \mu)I_5$ is equal to

$$ \begin{pmatrix} \beta_1(\mathcal{I}_1) \frac{S_1^*}{S_1^* + S_2^* + S_3^*} & \beta_1(\mathcal{I}_1) \frac{S_1^*}{S_1^* + S_2^* + S_3^*} & \beta_1(\mathcal{I}_1) \frac{S_1^*}{S_1^* + S_2^* + S_3^*} & 0 & 0 \\ \beta_1(\mathcal{I}_1) \frac{S_2^*}{S_1^* + S_2^* + S_3^*} & \beta_1(\mathcal{I}_1) \frac{S_2^*}{S_1^* + S_2^* + S_3^*} & \beta_1(\mathcal{I}_1) \frac{S_2^*}{S_1^* + S_2^* + S_3^*} & 0 & 0 \\ \beta_1(\mathcal{I}_1) \frac{S_3^*}{S_1^* + S_2^* + S_3^*} & \beta_1(\mathcal{I}_1) \frac{S_3^*}{S_1^* + S_2^* + S_3^*} & \beta_1(\mathcal{I}_1) \frac{S_3^*}{S_1^* + S_2^* + S_3^*} & 0 & 0 \\ 0 & 0 & 0 & \beta_2(\mathcal{I}_1) \frac{S_4^*}{S_4^* + S_5^*} & \beta_2(\mathcal{I}_1) \frac{S_4^*}{S_4^* + S_5^*} \\ 0 & 0 & 0 & \beta_2(\mathcal{I}_1) \frac{S_5^*}{S_4^* + S_5^*} & \beta_2(\mathcal{I}_1) \frac{S_5^*}{S_4^* + S_5^*} \end{pmatrix} $$
The ad hoc ordering introduced in Section 5.2.2 provides the second matrix $\mathcal{J}_{II,2}(\mathcal{X}^*) + (\gamma + \mu)I_5$ with a block diagonal structure, for example:

$$ \begin{pmatrix} \beta_1(\mathcal{I}_2) & 0 & 0 & 0 & 0 \\ 0 & \beta_2(\mathcal{I}_2) \frac{S_2^*}{S_2^*+S_5^*} & \beta_2(\mathcal{I}_2) \frac{S_2^*}{S_2^*+S_5^*} & 0 & 0 \\ 0 & \beta_2(\mathcal{I}_2) \frac{S_5^*}{S_2^*+S_5^*} & \beta_2(\mathcal{I}_2) \frac{S_5^*}{S_2^*+S_5^*} & 0 & 0 \\ 0 & 0 & 0 & \beta_3(\mathcal{I}_2) \frac{S_3^*}{S_3^*+S_4^*} & \beta_3(\mathcal{I}_2) \frac{S_3^*}{S_3^*+S_4^*} \\ 0 & 0 & 0 & \beta_3(\mathcal{I}_2) \frac{S_4^*}{S_3^*+S_4^*} & \beta_3(\mathcal{I}_2) \frac{S_4^*}{S_3^*+S_4^*} \end{pmatrix} $$

The value back in the reference ordering is such that $\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*) + (\gamma + \mu)I_5$ equals

$$ \begin{pmatrix} \beta_1(\mathcal{I}_2) & 0 & 0 & 0 & 0 \\ 0 & \beta_2(\mathcal{I}_2) \frac{S_2^*}{S_2^*+S_5^*} & 0 & 0 & 0 \\ 0 & 0 & \beta_3(\mathcal{I}_2) \frac{S_3^*}{S_3^*+S_4^*} & \beta_3(\mathcal{I}_2) \frac{S_3^*}{S_3^*+S_4^*} & 0 \\ 0 & 0 & \beta_3(\mathcal{I}_2) \frac{S_4^*}{S_3^*+S_4^*} & \beta_3(\mathcal{I}_2) \frac{S_4^*}{S_3^*+S_4^*} & 0 \\ 0 & \beta_2(\mathcal{I}_2) \frac{S_5^*}{S_2^*+S_5^*} & 0 & 0 & \beta_2(\mathcal{I}_2) \frac{S_5^*}{S_2^*+S_5^*} \end{pmatrix} $$

Notice that in general the two matrices $\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)$ and $\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)$ do not commute; so that the product of the exponentials is not the exponential of the sum.

* Two degenerate cases that present a constant contact pattern are easy to study, as they actually boil down to the situation of Theorem 3. These are the following.

* When $\tau = 0$, then (61) yields

$$ \mathcal{R}_0 = \rho \left( \mathcal{M}_2^\top e^{\mathcal{J}_{II,2}(\mathcal{X}^*)T} \mathcal{M}_2 \right) = \rho \left( e^{\mathcal{J}_{II,2}(\mathcal{X}^*)T} \right) = e^{-(\gamma+\mu)T} \max \left\{ e^{\beta_1(\mathcal{I}_2)T}; e^{\beta_2(\mathcal{I}_2)T}; e^{\beta_3(\mathcal{I}_2)T} \right\}. $$

* When $\tau = T$, then (61) gives

$$ \mathcal{R}_0 = \rho \left( e^{\mathcal{J}_{II,1}(\mathcal{X}^*)T} \right) = e^{-(\gamma+\mu)T} \max \left\{ e^{\beta_1(\mathcal{I}_1)T}; e^{\beta_2(\mathcal{I}_1)T} \right\}. $$

* Let us compute $e^{\mathcal{J}_{II,1}(\mathcal{X}^*)\tau}$ and $e^{\mathcal{J}_{II,2}(\mathcal{X}^*)(T-\tau)}$ in the general case, thanks to formula (58). In view of (65), one has

$$ \frac{S_1^*}{S_1^* + S_2^* + S_3^*} = \frac{1}{2}, \quad \frac{S_2^*}{S_1^* + S_2^* + S_3^*} = \frac{1}{12}, \quad \frac{S_3^*}{S_1^* + S_2^* + S_3^*} = \frac{5}{12}, \quad \frac{S_4^*}{S_4^* + S_5^*} = \frac{3}{10}, \quad \frac{S_5^*}{S_4^* + S_5^*} = \frac{7}{10}, \\ \frac{S_2^*}{S_2^* + S_5^*} = \frac{1}{8}, \quad \frac{S_5^*}{S_2^* + S_5^*} = \frac{7}{8}, \quad \frac{S_3^*}{S_3^* + S_4^*} = \frac{5}{8}, \quad \frac{S_4^*}{S_3^* + S_4^*} = \frac{3}{8}. $$

One deduces that

$$ e^{\mathcal{J}_{II,1}(\mathcal{X}^*)\tau} = e^{-(\gamma+\mu)\tau} \left( I_5 + \text{diag} \left\{ \frac{1}{12} (e^{\beta_1(\mathcal{I}_1)\tau} - 1) \begin{pmatrix} 6 & 6 & 6 \\ 1 & 1 & 1 \\ 5 & 5 & 5 \end{pmatrix} ; \frac{1}{10} (e^{\beta_2(\mathcal{I}_2)\tau} - 1) \begin{pmatrix} 3 & 3 \\ 7 & 7 \end{pmatrix} \right\} \right) $$

and

$$ e^{\mathcal{J}_{II,2}(\mathcal{X}^*)(T-\tau)} = e^{-(\gamma+\mu)(T-\tau)} \times \left( I_5 + \text{diag} \left\{ (e^{\beta_1(\mathcal{I}_2)(T-\tau)} - 1); \frac{1}{8}(e^{\beta_2(\mathcal{I}_2)(T-\tau)} - 1) \begin{pmatrix} 1 & 1 \\ 7 & 7 \end{pmatrix} ; \frac{1}{8}(e^{\beta_3(\mathcal{I}_2)(T-\tau)} - 1) \begin{pmatrix} 5 & 5 \\ 3 & 3 \end{pmatrix} \right\} \right). $$

Denoting for simplicity

$$ \phi_{11} := \frac{1}{12}(e^{\beta_1(\mathcal{I}_1)\tau} - 1), \quad \phi_{21} := \frac{1}{10}(e^{\beta_2(\mathcal{I}_2)\tau} - 1), $$

$$ \phi_{12} := (e^{\beta_1(\mathcal{I}_2)(T-\tau)} - 1), \quad \phi_{22} := \frac{1}{8}(e^{\beta_2(\mathcal{I}_2)(T-\tau)} - 1), \quad \phi_{32} := \frac{1}{8}(e^{\beta_3(\mathcal{I}_2)(T-\tau)} - 1), $$
with the expression of $M_2$ in (60), one finally obtains that

$$e^{(\gamma+\mu)T} e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(T-\tau)} e^{\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)\tau} = \begin{pmatrix} 1+\phi_{12} & 0 & 0 & 0 & 0 \\ 0 & 1+\phi_{22} & 0 & 0 & \phi_{22} \\ 0 & 0 & 1+5\phi_{32} & 5\phi_{32} & 0 \\ 0 & 0 & 3\phi_{32} & 1+3\phi_{32} & 0 \\ 0 & 7\phi_{22} & 0 & 0 & 1+7\phi_{22} \end{pmatrix} \times \begin{pmatrix} 1+6\phi_{11} & 6\phi_{11} & 6\phi_{11} & 0 & 0 \\ \phi_{11} & 1+\phi_{11} & \phi_{11} & 0 & 0 \\ 5\phi_{11} & 5\phi_{11} & 1+5\phi_{11} & 0 & 0 \\ 0 & 0 & 0 & 1+3\phi_{21} & 3\phi_{21} \\ 0 & 0 & 0 & 7\phi_{21} & 1+7\phi_{21} \end{pmatrix}. \quad (66)$$

• By comparison, the stationary system usually considered when assuming fast commuter movements in metapopulation models [30, Section 7.2.1.4] is obtained by considering the weighted mean of the transmission terms. Accordingly, the blocks of the Jacobian matrix of the corresponding autonomous system are obtained as weighted means, yielding an approximation of $\rho(e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(T-\tau)}e^{\hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)\tau})$ in (61) by the expression

$$\rho(e^{\hat{\mathcal{J}}_{II,2}(\mathcal{X}^*)(T-\tau)} + \hat{\mathcal{J}}_{II,1}(\mathcal{X}^*)\tau).$$

Here,

$$\tau \hat{\mathcal{J}}_{II,1}(\mathcal{X}^*) + (T - \tau) \hat{\mathcal{J}}_{II,2}(\mathcal{X}^*) = -T(\gamma + \mu)I_5 + \tau \begin{pmatrix} \frac{1}{2}\beta_1(\mathcal{I}_1) & \frac{1}{2}\beta_1(\mathcal{I}_1) & \frac{1}{2}\beta_1(\mathcal{I}_1) & 0 & 0 \\ \frac{1}{12}\beta_1(\mathcal{I}_1) & \frac{1}{12}\beta_1(\mathcal{I}_1) & \frac{1}{12}\beta_1(\mathcal{I}_1) & 0 & 0 \\ \frac{5}{12}\beta_1(\mathcal{I}_1) & \frac{5}{12}\beta_1(\mathcal{I}_1) & \frac{5}{12}\beta_1(\mathcal{I}_1) & 0 & 0 \\ 0 & 0 & 0 & \frac{3}{10}\beta_2(\mathcal{I}_1) & \frac{3}{10}\beta_2(\mathcal{I}_1) \\ 0 & 0 & 0 & \frac{7}{10}\beta_2(\mathcal{I}_1) & \frac{7}{10}\beta_2(\mathcal{I}_1) \end{pmatrix} + (T - \tau) \begin{pmatrix} \beta_1(\mathcal{I}_2) & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{8}\beta_2(\mathcal{I}_2) & 0 & 0 & \frac{1}{8}\beta_2(\mathcal{I}_2) \\ 0 & 0 & \frac{5}{8}\beta_3(\mathcal{I}_2) & \frac{5}{8}\beta_3(\mathcal{I}_2) & 0 \\ 0 & 0 & \frac{7}{8}\beta_3(\mathcal{I}_2) & \frac{7}{8}\beta_3(\mathcal{I}_2) & 0 \\ 0 & \frac{7}{8}\beta_2(\mathcal{I}_2) & 0 & 0 & \frac{7}{8}\beta_2(\mathcal{I}_2) \end{pmatrix}. \quad (67)$$

• As a last remark and an illustration of Corollary 5, consider now the situation where the spreading is uniform in the whole population during each of the two phases, that is

$$\beta_1(\mathcal{I}_1) = \beta_2(\mathcal{I}_1) =: \beta(\mathcal{I}_1), \quad \beta_1(\mathcal{I}_2) = \beta_2(\mathcal{I}_2) = \beta_3(\mathcal{I}_2) =: \beta(\mathcal{I}_2).$$

One checks that the vector $1_5^I$ is a left eigenvector of the matrix in (66), associated to the eigenvalue

$$e^{(\gamma+\mu)T} e^{\beta(\mathcal{I}_2)(T-\tau)} e^{\beta(\mathcal{I}_1)\tau}. \quad (68)$$

The values of $G_1, G_2$ provided by (54) are respectively $\beta(\mathcal{I}_1), \beta(\mathcal{I}_2)$, see (64). One then verifies that the value in (68) is identical to the one provided by formula (55), as expected by Corollary 5.

As a side comment, this value is also equal to the one obtained through use of the stationary approximation, namely the spectral radius of the matrix in (67).

# 6 Conclusion

A new framework for the modelling and analysis of the spread of epidemic diseases in complex commuting populations was proposed in this paper. It consists in a general class of compartmental ODE models, where the population is divided into homogeneous subpopulations of different sizes. At any time is defined a partition of the set of these subpopulations, whose classes are comprised of those currently at a common physical location. These partitions are piecewise constant and periodic in time, in order to model repetitive contact patterns. Each subpopulation thus represents here a maximal set of individuals
perpetually and homogeneously in contact from the point of view of epidemic spread. On top of this,
every subpopulation is further subdivided, as is usually the case, into several compartments accounting
for the epidemiological status, and more generally for any heterogeneity trait pertinent to describe the
disease spread with all the desired precision [25].

In order to allow for adequate representation of the mixtures and separations of the subpopulations,
it was necessary to revisit the epidemiological model. To this end, we introduced a new general setting
showcasing the special role of the force of infection. It led to consider compartmental models that are
jointly linear with respect to the state variable, and affine with respect to the force of infection (itself a
function of the state variable).

Based on this formulation, a limited set of assumptions was adopted, which essentially transposes the hypotheses introduced by van den Driessche and Watmough [42] into this larger setting. The behaviour around the disease-free equilibrium points of such systems, and especially the possibility of epidemic bursts, has then been characterised in terms of spectral radius of matrix products. In a didactic purpose, we gradually considered the case of a unique subpopulation (Theorem 2); then systems with several subpopulations in stationary contact pattern (Theorem 3); and finally in full generality so-called ‘epidemiological models in commuting population’, that is, systems with several subpopulations interacting through periodically alternating mixing (Theorem 4). As an interesting consequence, this allowed to show (Corollary 5) that, when the transmission conditions are identical in the whole population, then the value of the basic reproduction number is actually independent of the contact pattern, and is the same than what it would be if all subpopulations were brought together. Up to our knowledge, such result is original. An elementary example was provided to illustrate the computations underlying the results.

The framework proposed in this paper is believed to be original, and may serve as a basis for studying
general model properties. First of all, the proposed class of compartmental models provides a significant
extension and unification of many models previously published in the literature. This larger setting should
help identify the key characteristics of epidemiological models and analyse their behaviour. Second, the
modelling of the spatial aspect in terms of *partitions of the population* rather than in terms of *fixed
physical locations* seems well adapted to the description of the spread of an epidemic in commuting
populations. This option allows to keep track of the identity of each group, while reducing the conceptual
and computational complexity necessary to represent the system evolution compared to the classical
modelling choices. The number of variables is here proportional to the number of subpopulation classes
defined by the partition, which represent the ‘active’ locations, and not of the number of locations or
patches, or a power of this quantity. This is particularly suitable for an exact description of time-
varying conditions of disease propagation, by opposition with the stationary, averaged, settings usually
considered. It allows to analyse theoretically and assess numerically the effects of such approximations.

The proposed modelling framework should allow for extensive computational investigation. It permits
for example to address the issues of robustness with respect to the parametric uncertainties of the model
— be they related to the infection process (as e.g. the contact rates, or the latency, recovery and loss of
immunity rates) or to the commuting process and the mixing times. In this framework, the question of
the stability robustness of a DFE boils down to the study of the extremal value of the spectral radius
over a given set of matrices, or of the generalized spectral radius [46, 29] of a specific set of matrices.
Also, it should permit comprehensive exploration of the final epidemic size for such periodic commuting
systems. This relevant quantity represents the total number of individuals infected during an outbreak
(see e.g. [2] and references therein), and its computation is quite complicated. These topics will be the
subject of forthcoming research.

Acknowledgements

Pierre-Alexandre Bliman is grateful to Carlos Canudas de Wit (GIPSA-lab) for having the opportunity to work together during COVID pandemic in the Inria working group on Healthy mobility control. The discussion and reflection initiated at that time nurtured the inspiration for the present paper. He is also grateful to Abderrahman Iggidr (Inria), Alain Rapaport (Inrae) and Gauthier Sallet (Institut Élie Cartan de Lorraine) for valuable discussions on the topics considered here. Boureima Sangaré is deeply thankful to the Laboratoire de Mathématiques de Besançon for hosting him for three months and supporting
his trip to Sorbonne Université. Assane Savadogo is very much thankful to Prof. Christophe Cambier (Sorbonne Université), coordinator of the PDI-MSC program of IRD which supported this work. He is also thankful to Inria Paris centre for financial support. Last, the funding of ANR is acknowledged, under the project ANR-23-CE48-0004.

## References

[1] R. M. ANDERSON AND R. M. MAY, *Infectious diseases of humans: dynamics and control*, Oxford University Press, 1991.

[2] V. ANDREASEN, *The final size of an epidemic and its relation to the basic reproduction number*, Bulletin of Mathematical Biology, 73 (2011), pp. 2305-2321.

[3] J. ARINO, *Diseases in metapopulations*, in Modeling and dynamics of infectious diseases, World Scientific, 2009, pp. 64-122.

[4] J. ARINO, J. R. DAVIS, D. HARTLEY, R. JORDAN, J. M. MILLER, AND P. VAN DEN DRIESSCHE, *A multi-species epidemic model with spatial dynamics*, Mathematical Medicine and Biology, 22 (2005), pp. 129-142.

[5] J. ARINO AND P. VAN DEN DRIESSCHE, *Disease spread in metapopulations*, Fields Institute Communications, 48 (2006), pp. 1-13.

[6] N. BACAËR, *Approximation of the basic reproduction number R₀ for vector-borne diseases with a periodic vector population*, Bulletin of Mathematical Biology, 69 (2007), pp. 1067-1091.

[7] N. BACAËR AND E. H. AIT DADS, *On the biological interpretation of a definition for the parameter R₀ in periodic population models*, Journal of Mathematical Biology, 65 (2012), pp. 601-621.

[8] N. BACAËR AND M. G. M. GOMES, *On the final size of epidemics with seasonality*, Bulletin of Mathematical Biology, 71 (2009), pp. 1954-1966.

[9] N. BACAËR AND R. OUIFKI, *Growth rate and basic reproduction number for population models with a simple periodic factor*, Mathematical Biosciences, 210 (2007), pp. 647-658.

[10] A. BERMAN AND R. J. PLEMMONS, *Nonnegative matrices in the mathematical sciences*, SIAM, 1994.

[11] S. M. BLOWER, P. M. SMALL, AND P. C. HOPEWELL, *Control strategies for tuberculosis epidemics: new models for old problems*, Science, 273 (1996), pp. 497-500.

[12] B. BONZI, A. FALL, A. IGGIDR, AND G. SALLET, *Stability of differential susceptibility and infectivity epidemic models*, Journal of Mathematical Biology, 62 (2011), pp. 39-64.

[13] F. BRAUER AND C. CASTILLO-CHAVEZ, *Mathematical models for communicable diseases*, SIAM, 2012.

[14] V. CAPASSO, *Mathematical structures of epidemic systems*, vol. 97, Springer, 1993.

[15] V. CAPASSO AND G. SERIO, *A generalization of the Kermack-McKendrick deterministic epidemic model*, Mathematical Biosciences, 42 (1978), pp. 43-61.

[16] C. CASTILLO-CHAVEZ AND Z. FENG, *To treat or not to treat: The case of tuberculosis*, Journal of Mathematical Biology, 35 (1997), pp. 629-656.

[17] S. COLLINSON AND J. M. HEFFERNAN, *Modelling the effects of media during an influenza epidemic*, BMC Public Health, 14 (2014), pp. 1-10.

[18] J. CUI, Y. SUN, AND H. ZHU, *The impact of media on the control of infectious diseases*, Journal of Dynamics and Differential Equations, 20 (2008), pp. 31-53.
[19] O. DIEKMANN AND J. A. P. HEESTERBEEK, *Mathematical epidemiology of infectious diseases: model building, analysis and interpretation*, vol. 5, John Wiley & Sons, 2000.

[20] Z. FENG AND J. X. VELASCO-HERNÁNDEZ, *Competitive exclusion in a vector-host model for the dengue fever*, Journal of Mathematical Biology, 35 (1997), pp. 523–544.

[21] H. GUO, M. Y. LI, AND Z. SHUAI, *Global dynamics of a general class of multistage models for infectious diseases*, SIAM Journal on Applied Mathematics, 72 (2012), pp. 261–279.

[22] I. HANSKI AND M. GILPIN, eds., *Metapopulation Biology. Ecology, Genetics, and Evolution*, Academic Press, 1997.

[23] H. W. HETHCOTE, *An immunization model for a heterogeneous population*, Theoretical Population Biology, 14 (1978), pp. 338–349.

[24] H. W. HETHCOTE, *A thousand and one epidemic models*, in Frontiers in Mathematical Biology, Springer, 1994, pp. 504–515.

[25] H. W. HETHCOTE AND J. W. VAN ARK, *Epidemiological models for heterogeneous populations: Proportionate mixing, parameter estimation, and immunization programs*, Mathematical Biosciences, 84 (1987), pp. 85–118.

[26] J. M. HYMAN, J. LI, AND E. A. STANLEY, *The differential infectivity and staged progression models for the transmission of HIV*, Mathematical Biosciences, 155 (1999), pp. 77–109.

[27] A. IGGIDR, J. MBANG, G. SALLET, AND J.-J. TEWA, *Multi-compartment models*, Discrete and Continuous Dynamical Systems - Series S, 2007 (2007), pp. 506–519.

[28] J. A. JACQUEZ, C. P. SIMON, J. KOOPMAN, L. SATTENSPIEL, AND T. PERRY, *Modeling and analyzing HIV transmission: The effect of contact patterns*, Mathematical Biosciences, 92 (1988), pp. 119–199.

[29] R. JUNGERS, *The joint spectral radius: theory and applications*, vol. 385, Springer Science & Business Media, 2009.

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

[31] D. H. KNIPL, G. RÖST, AND J. WU, *Epidemic spread and variation of peak times in connected regions due to travel-related infections—dynamics of an antigravity-type delay differential model*, SIAM Journal on Applied Dynamical Systems, 12 (2013), pp. 1722–1762.

[32] Y. LI, C. MA, AND J. CUI, *The effect of constant and mixed impulsive vaccination on SIS epidemic models incorporating media coverage*, The Rocky Mountain Journal of Mathematics, (2008), pp. 1437–1455.

[33] W.-M. LIU, H. W. HETHCOTE, AND S. A. LEVIN, *Dynamical behavior of epidemiological models with nonlinear incidence rates*, Journal of Mathematical Biology, 25 (1987), pp. 359–380.

[34] X. LIU AND L. YANG, *Stability analysis of an SEIQV epidemic model with saturated incidence rate*, Nonlinear Analysis: Real World Applications, 13 (2012), pp. 2671–2679.

[35] O. MASON AND R. SHORTEN, *On linear copositive Lyapunov functions and the stability of switched positive linear systems*, IEEE Transactions on Automatic Control, 52 (2007), pp. 1346–1349.

[36] J. C. MILLER, *A note on the derivation of epidemic final sizes*, Bulletin of Mathematical Biology, 74 (2012), pp. 2125–2141.

[37] Y. NAKATA AND G. RÖST, *Global analysis for spread of infectious diseases via transportation networks*, Journal of Mathematical Biology, 70 (2015), pp. 1411–1456.

[38] L. SATTENSPIEL AND K. DIETZ, *A structured epidemic model incorporating geographic mobility among regions*, Mathematical Biosciences, 128 (1995), pp. 71–91.
[39] L. SATTENSPIEL AND C. P. SIMON, *The spread and persistence of infectious diseases in structured populations*, Mathematical Biosciences, 90 (1988), pp. 341–366.

[40] C. SUN, W. YANG, J. ARINO, AND K. KHAN, *Effect of media-induced social distancing on disease transmission in a two patch setting*, Mathematical Biosciences, 230 (2011), pp. 87–95.

[41] H. R. THIEME, *Mathematics in population biology*, vol. 12, Princeton University Press, 2003.

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

[43] W. WANG, *Epidemic models with nonlinear infection forces*, Mathematical Biosciences & Engineering, 3 (2005), pp. 267–279.

[44] W. WANG AND X.-Q. ZHAO, *An epidemic model in a patchy environment*, Mathematical biosciences, 190 (2004), pp. 97–112.

[45] ——, *Threshold dynamics for compartmental epidemic models in periodic environments*, Dynamics and Differential Equations, 20 (2008), pp. 699–717.

[46] F. WIRTH, *The generalized spectral radius and extremal norms*, Linear Algebra and its Applications, 342 (2002), pp. 17–40.

[47] D. XIAO AND S. RUAN, *Global analysis of an epidemic model with nonmonotone incidence rate*, Mathematical Biosciences, 208 (2007), pp. 419–429.