EPIDEMIC MODELS UNDER MOBILITY ON MULTI-LAYER NETWORKS

By

Vishal Abhishek

A THESIS

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Mechanical Engineering – Master of Science

2021

ABSTRACT

EPIDEMIC MODELS UNDER MOBILITY ON MULTI-LAYER NETWORKS

By

Vishal Abhishek

We study epidemic spreading models namely, SIS and SIR models, under mobility on multi-

layer networks. In particular, we consider a patchy environment in which each patch com-

prises individuals belonging the different classes, e.g., individuals in different socio-economic

strata. We model the mobility of individuals of each class across different patches through an

associated Continuous Time Markov Chain (CTMC). The topology of these multiple CTMCs

constitute the multi-layer network of mobility. At each time, individuals in the multi-layer

network of spatially-distributed patches move according to their CTMC and subsequently

interact with the local individuals in the patch according to SIS or SIR models. We es-

tablish the existence of various equilibria under different parameter regimes and establish

their (almost) global asymptotic stability using Lyapunov techniques. We also derive simple

conditions that highlight the influence of the multi-layer network on the stability of these

equilibria. We numerically illustrate that the derived model provides a good approximation

to the stochastic model with a finite population and also demonstrate the influence of the

multi-layer network structure.

Next, we extend some of the results to the case of weakly connected networks. Here, we

use the notion of strongly connected components and input to state stability to study the

stability of equilibria.

Finally, we consider a resource allocation problem to maximize the rate of convergence

to an equilibrium. We show that under certain assumptions the problem can be formulated

as a geometric program. We provide numerical illustrations to corroborate the results.

ACKNOWLEDGEMENTS

I would like to express my sincere thanks to my advisor Prof. Vaibhav Srivastava, without

whom this work would not have been possible. I express my deepest gratitude for his contin-

uous guidance and motivation. I would like to thank my committee members Prof. Ranjan

Mukherjee and Prof. Xiaobo Tan for giving their valuable time and inputs.

I would also like to thank all the members of D-CYPHER Lab for making my work at the

lab fruitful as well as cheerful.

Lastly, I would like to thank my parents for their sacrifices and continuous support.

Vishal Abhishek

iii

TABLE OF CONTENTS

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Epidemic Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Epidemic Models under Mobility . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Resource Allocation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Mathematical Notations
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 2 EPIDEMIC MODELS UNDER MOBILITY ON SINGLE-LAYER

NETWORKS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 SIS Model under Markovian Mobility . . . . . . . . . . . . . . . . . . . . . .
2.2 Analysis of SIS Model under Markovian Mobility . . . . . . . . . . . . . . .
2.3 Numerical Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 3 EPIDEMIC MODELS UNDER MOBILITY ON MULTI-LAYER

NETWORKS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1 SIS Model under Multi-layer Markovian Mobility . . . . . . . . . . . . . . .
3.2 Analysis of SIS Model under Multi-layer Markovian Mobility . . . . . . . . .
3.3 Numerical Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 SIR Model under Multi-layer Markovian Mobility . . . . . . . . . . . . . . .
3.5 Numerical Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 4 SIS EPIDEMIC MODEL UNDER MOBILITY ON MULTI-LAYER

NETWORKS: WEAKLY CONNECTED CASE . . . . . . . . . . . .
4.1 SIS Model under Mobility with Non-strongly Connected Layers
. . . . . . .
4.2 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 5 RESOURCE ALLOCATION . . . . . . . . . . . . . . . . . . . . . .
5.1 Resource allocation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1
Single layer mobility . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 Multilayer mobility . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Numerical Illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 6 CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . .

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

v

1
1
2
3
4
4

5
5
7
15

18
18
20
30
32
37

39
39
41

44
44
44
46
46

49

50

61

iv

LIST OF FIGURES

Figure 2.1: Stochastic simulation of epidemic spread under mobility. Line graph,
, pi(0) = 0.01. Each iteration in stochas-

n = 20, ν(i) = 0.2, qij = ν(i)
Dout
tic model corresponds to time-step 0.01 sec.

. . . . . . . . . . . . . . . .

Figure 2.2: Simulation of deterministic model of epidemic spread under mobility,
with same equilibrium distribution of population over 4 different graph
structure with stable endemic equilibrium. n = 20, pi(0) = 0.01.

. . . . .

Figure 2.3: Stable disease-free equilibrium with curing rates computed as per the λ2
sufficient condition (statement (iv), Corollary 1) for stability of disease-
free equilibrium. Graph: Complete, n = 20, pi(0) = 0.01.

. . . . . . . . .

Figure 3.1: Stochastic simulation of epidemic spread under mobility. Complete-Line
, pi(0) = 0.01. Each iteration in

graphs, n = 20, ν(i) = 0.2, qij = ν(i)
Dout
stochastic model corresponds to time-step 0.01 sec.

. . . . . . . . . . . .

Figure 3.2: Simulation of deterministic model of epidemic spread under 2 layer mo-
bility, over different graph structure with stable endemic equilibrium.
n = 10, pi(0) = 0.01.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.3: Simulation of deterministic model of SIR epidemic spread under 2 layer

mobility, over Line-Ring graph structure. n = 10, pi(0) = 0.01. . . . . . .

Figure 4.1: Non-strongly connected multilayer network with shaded sink components

Figure 5.1: Decay rate λ against maximum allowable cost C. Graph1: Line, Graph

2: Ring, n = 10, νi = 0.2.

. . . . . . . . . . . . . . . . . . . . . . . . . .

15

16

17

30

31

38

39

48

v

CHAPTER 1

INTRODUCTION

In this thesis, we study epidemic models: namely, SIS (Susceptible-Infected-Susceptible) and

SIR (Susceptible-Infected-Recovered) models, under mobility on multi-layer networks. The

following sections introduce the work done.

1.1 Epidemic Models

Contagion dynamics are used to model a variety of phenomenon such as spread of influ-

ence, disease and rumors. Epidemic propagation models are as a class of contagion models

that have been used in the context of disease spread [1, 2], spread of computer viruses [3, 4],

routing in mobile communication networks [5], and spread of rumors [6]. Epidemic prop-

agation in patchy environments refers to the epidemic spread process in an environment

comprised of disjoint spatially distributed regions (patches). In these models, individuals

interact within each patch and also move across different patches according to a CTMC.

In this thesis, we consider a generalized epidemic propagation model in patchy envi-

ronment in which individuals within each patch belong to multiple classes, and individuals

within each class move according to an associated CTMC. This leads to a multi-layer mobility

model and we study its interaction with epidemic propagation. Using Lyapunov techniques,

we characterize the steady state behavior of the model under different parameter regimes

and characterize the influence of mobility on epidemic dynamics.

Epidemic models have been extensively studied in the literature. The two most widely

studied models are SIS (Susceptible-Infected-Susceptible) and SIR (Susceptible-Infected-

Recovered) models, wherein individuals are classified into one of the three categories: sus-

ceptible, infected or recovered. In contrast to classical SIS/SIR models where the dynamics

of the fraction of the population in each category [2] is studied, the Network models consider

patches clustered into different nodes, and the patch-level dynamics is determined by the

1

local SIS/SIR interactions as well as the interactions with neighboring patches in the network

graph [7, 8, 8, 9, 10]. Authors in [11, 12] study network epidemic dynamics in discrete time

setting.

Some common generalizations of the SIR/SIS models include: SEIR model [7, 13], where

an additional classification “exposed" is introduced, SIRS [2, 12], where individuals get tem-

porary immunity after recovery and then become susceptible again, and SIRI [14, 15], where

after recovery, agents become susceptible with a different rate of infection. The network

epidemic dynamics have also been studied for time-varying networks [16, 17, 18].

1.2 Epidemic Models under Mobility

The terms population dispersal and network mobility have been used interchangeably in

the literature. Epidemic spread under mobility has been modeled and analyzed as reaction-

diffusion process in [19, 20]. Epidemic spread in a patchy environment with population

dispersal has been modeled and studied in [21, 22, 23].

In these works, the mobility or

dispersal patterns depend on the state (susceptible or infected) of the individuals, and con-

ditions for global stability of the disease-free equilibrium and an endemic equilibrium are

derived. When the mobility patterns are identical for all individuals, then these models

reduce to a model similar to the single-layer version of the multi-layer model studied in this

thesis.

Epidemic spread with mobility on a multiplex network of patches has been modeled and

studied in [24]. Authors of this work consider a discrete model in which, at each time,

individuals randomly move to another node, participate in epidemic propagation and then

return to their home node. A multi-species SEIR epidemic model with population dispersal

has been analyzed in [25] and conditions for the global stability of a disease-free equilibrium

are derived. Stability results for endemic equilibrium for the single species case are also

derived. The population dispersal model in this work is identical to the multi-layer mobility

model studied in this thesis and it reduces to the single layer mobility model for a single

2

species.

In contrast to this work, we focus on SIS epidemic model and study endemic

equilibria for the general multiplex network.

In this work, we study a coupled epidemic-mobility model comprised of a set of patches

located in a multi-layer network of spatially distributed regions.

Individuals within each

patch (region) can travel across regions according to a Continuous Time Markov Chain

(CTMC) characterising their mobility pattern and upon reaching a new region participate

in the local SIS epidemic process. We extend the results for the deterministic network SIS

and SIR model [8, 9, 10, 11] to the proposed model and characterize its steady state and

stability properties.

We extend the results for the SIS model to the case of weakly connected networks. Here,

we use the notion of strongly connected components to divide the nodes of the multi-layer

networks into several batches and analyze the stability of epidemic model over these nodes

only. We use the notion of input to state stability to incorporate the effects of the remaining

nodes.

1.3 Resource Allocation

We formulate a resource allocation problem for the SIS model under markovian mobility

where we allocate resources so as to have the disease free equilibrium (DFE) as the stable

equilibrium. We consider two types of resources:

i. Preventive resource: This can be used to change the infection rate such that βi ∈

[βi, ¯βi]. This resource is applied to a node with cost function fi(βi).

ii. Corrective resource: This can be used to change the recovery rate such that δi ∈ [δi, ¯δi].

The corresponding cost function is gi(δi).

We show that under certain assumptions the problem can be formulated as a geometric

program.

3

1.4 Thesis Organization

Chapter 2 discusses SIS epidemic model under mobility over single layer networks. Chap-

ter 3 discusses epidemic models: SIS and SIR under mobility on multi-layer networks. Chap-

ter 4 extends the results for weakly connected networks. Chapter 5 discusses resource allo-

cation formulation and its solution as a geometric program. Finally, Chapter 6 concludes

the thesis. First, we give below mathematical notations used in this thesis:

1.5 Mathematical Notations

For any two real vectors x, y ∈ Rn, we denote:
x (cid:29) y, if xi > yi for all i ∈ {1, . . . , n},
x ≥ y, if xi ≥ yi for all i ∈ {1, . . . , n},
x > y, if xi ≥ yi for all i ∈ {1, . . . , n} and x (cid:54)= y.
For a square matrix G, radial abscissa µ : Rn×n → R is defined by

µ(G) = max{Re(λ) | λ is an eigenvalue of G},

where Re(·) denotes the real part of the argument. Spectral radius ρ is defined by

ρ(G) = max{|λ| | λ is an eigenvalue of G},

where |(·)| denotes the absolute value of the argument.
For any vector x = [x1, . . . , xn](cid:62), X = diag(x) is a diagonal matrix with Xii = xi for all
i ∈ {1, . . . , n}.

4

CHAPTER 2

EPIDEMIC MODELS UNDER MOBILITY ON SINGLE-LAYER

NETWORKS

In this chapter, we study the SIS model under mobility on a single-layer network. The sin-

gle layer corresponds to the presence of one mobility pattern with its characteristic Markov

Chain and digraph. We assume that the digraph is strongly connected. We give complete

characterization of the equilibria of the system, including its existence, uniqueness and sta-

bility. We also give numerical illustrations to support our results.

2.1 SIS Model under Markovian Mobility

We consider n sub-populations of individuals that are located in distinct spatial regions.

We assume the individuals within each sub-population can be classified into two categories:
(i) susceptible, and (ii) infected. Let pi ∈ [0, 1] (respectively, 1−pi) be the fraction of infected
(respectively, susceptible) individuals within sub-population i ∈ {1, . . . , n}. We assume that
the individuals within each sub-population can travel to regions associated with other sub-
populations. Let the connectivity of these regions be modeled by a digraph G = (V,E), where
V = {1, . . . , n} is the node set and E ⊂ V × V is the edge set. We model the mobility of
individuals on graph G using a Continuous Time Markov Chain (CTMC) with a stationary
generator matrix Q, whose (i, j)-th entry is qij. The entry qij ≥ 0, i (cid:54)= j, is the instantaneous
transition rate from node i to node j, and −qii = νi is the total rate of transition out of node
j(cid:54)=i qij. Here, qij > 0, if (i, j) ∈ E; and qij = 0, otherwise. Let xi(t) ∈ (0, 1)
be the fraction of the total population that constitutes the sub-population at node i at time

i, i.e., νi =(cid:80)
t. It follows that(cid:80)n

i=1 xi = 1. Define p := [p1, . . . , pn](cid:62) and x := [x1, . . . , xn](cid:62).

We model the interaction of mobility with the epidemic process as follows. At each time t,
individuals at each node move on graph G according to the CTMC with generator matrix Q
and interact with individuals within their current node according to an SIS epidemic process.

5

For the epidemic process at node i, let βi > 0 and δi ≥ 0 be the infection and recovery rate,
respectively. We let B > 0 and D ≥ 0 be the positive and non-negative diagonal matrices
with entries βi and δi, i ∈ {1, . . . , n}, respectively. Similarly we define P as a diagonal matrix
with entries pi. We now derive the continuous time dynamics that captures the interaction
of mobility and the SIS epidemic dynamics.

Proposition 1 (SIS model under mobility) The dynamics of the fractions of the in-

fected sub-population p and the fractions of the total population x that constitute the sub-

population at each node under Markovian mobility model with generator matrix Q, and in-

fection and recovery matrices B and D, respectively, are

where L(x) is a matrix with entries

˙p = (B − D − L(x))p − P Bp
˙x = Q(cid:62)x,
(cid:80)
−qji



j(cid:54)=i qji

if i = j,

xj
xi

,

xj
xi

,

lij(x) =

otherwise.

(2.1a)

(2.1b)

Proof: Consider a small time increment h > 0 at time t. Then the fraction of the total

population present at node i after the evolution of CTMC in time-interval [t, t + h) is

xi(t + h) = xi(t)(1 − νih) +

qjixj(t)h + o(h).

(2.2)

(cid:88)

j(cid:54)=i

Individuals within each node interact according to SIS dynamics. Thus, the fraction of

infected population present at node i is:

(cid:88)
xi(t + h)pi(t + h) = −xi(t)δipi(t)h + xi(t)βipi(t)(1 − pi(t))h
+ xi(t)pi(t)(1 − νih) +

j(cid:54)=i

qjipj(t)xj(t)h + o(h).

(2.3)

The first two terms on the right side of (2.3) correspond to epidemic process within each

node, whereas the last two terms correspond to infected individuals coming from other nodes

6

due to mobility. Using the expression of xi from (2.2) in (2.3) and taking the limit h → 0+
gives

(cid:88)

j(cid:54)=i

qji(pj − pi)

xj
xi

.

(2.4)

˙pi

=

−δipi + βipi(1 − pi) +

Similarly taking limits in (2.2) yields

˙xi = −νixi +

(cid:88)

j(cid:54)=i

qjixj.

Rewriting (2.4) and (2.5) in vector form establishes the proposition.

(2.5)

(cid:3)

Remark 1 (Comparison with other models) The epidemic propagation models in the-

oretical ecology incorporate spatial aspects by using a partial differential equation that is

obtained by adding a spatial diffusion operator to the infected population dynamics [26].

Since, Laplacian matrix is a diffusion operator on a graph, dynamics (2.1) can be interpreted

as a network equivalent of the population models with spatial aspects. The dependence of

the Laplacian matrix on x in (2.1) is more general than the constant diffusion coefficient
(cid:3)

discussed in [26].

2.2 Analysis of SIS Model under Markovian Mobility

In this section, we analyze the SIS model under mobility (2.1) under the following stan-

dard assumption:

Assumption 1 Digraph G is strongly connected which is equivalent to matrix Q being irre-
(cid:3)
ducible [27].
Let v be the right eigenvector of Q(cid:62) associated with eigenvalue at 0. We assume that
v is scaled such that its inner product with the associated left eigenvector 1n is unity,
i.e., 1(cid:62)
n v = 1. We call an equilibrium point (p∗, x∗), an endemic equilibrium point, if at
equilibrium the disease does not die out, i.e., p∗ (cid:54)= 0, otherwise, we call it a disease-free
equilibrium point. Let L∗ := L(x∗) = L(v).

7

Theorem 1 (Existence and Stability of Equilibria) For the SIS model under Marko-

vian mobility (2.1) with Assumption 1, the following statements hold

i. if p(0) ∈ [0, 1]n, then p(t) ∈ [0, 1]n for all t > 0. Also, if p(0) > 0n, then p(t) (cid:29) 0n

for all t > 0;

ii. the model admits a disease-free equilibrium at (p∗, x∗) = (0n, v);
iii. the model admits an endemic equilibrium at (p∗, x∗) = (¯p, v), ¯p (cid:29) 0, if and only if

µ(B − D − L∗) > 0;

iv. the disease-free equilibrium is globally asymptotically stable if and only if µ(B − D −

L∗) ≤ 0 and is unstable otherwise;

v. the endemic equilibrium is almost globally asymptotically stable if µ(B − D − L∗) > 0

with region of attraction p(0) ∈ [0, 1]n such that p(0) (cid:54)= 0n.

Proof: The first part of statement (i) follows from the fact that ˙p is either directed

tangent or inside of the region [0, 1]n at its boundary which are surfaces with pi = 0 or 1 .
For the second part of (i), we rewrite (3.1a) as:

˙p = (B(I − P ) + A(x))p − E(t)p

where L(x) = C(x) − A(x) with C(x) composed of the diagonal terms of L(x), A(x) is
the non-negative matrix corresponding to the off-diagonal terms, and E(t) = C(x(t)) + D

(cid:82) t

is a diagonal matrix. Now, consider a variable change y(t) := e

0 E(t)dtp(t). The rest of the

proof is same as in [8, Theorem 4.2 (i)].

The second statement follows by inspection.

The proof of the third statement is presented in Appendix A.1.

Stability of disease-free equilibria: To prove the fourth statement, we begin by estab-

lishing sufficient conditions for instability. The linearization of (2.1) at (p, x) = (0, v) is

8

 =

 ˙p

˙x

B − D − L∗ 0n×n

Q(cid:62)

0n×n



 .

p

x

(2.6)

Since the system matrix in (2.6) is block-diagonal, its eigenvalues are the eigenvalues of
the block-diagonal sub-matrices. Further, since radial abscissa µ(Q(cid:62)) is zero, a sufficient
condition for instability of the disease-free equilibrium is that µ(B − D − L∗) > 0.

For the case of µ(B − D − L∗) ≤ 0, we now show that the disease-free equilibrium is
a globally asymptotically stable equilibrium. Since (B − D − L∗) is an irreducible Metzler
matrix with µ(B − D − L∗) ≤ 0, there exists a positive diagonal matrix R such that

R(B − D − L∗) + (B − D − L∗)(cid:62)R = −K,

where K is a positive semi-definite matrix [10, Proposition 1 (iv), Lemma A.1]. Define
˜L := L(x) − L∗ and r := (cid:107)R(cid:107), where (cid:107) · (cid:107) denotes the the induced two norm of the matrix.
Since x(0) (cid:29) 0, under Assumption 1, xi(t) is lower bounded by some positive constant
and hence, ˜L is bounded and continuously differentiable. Since x is bounded and exponen-
tially converges to x∗, it follows that (cid:107) ˜L(x)(cid:107) locally exponentially converges to (cid:107) ˜L(x∗)(cid:107) = 0

and(cid:82) t
Consider the Lyapunov-like function V (p, t) = p(cid:62)Rp− 2nr(cid:82) t

0 (cid:107) ˜L(cid:107)dt is bounded for all t > 0.

0 (cid:107) ˜L(cid:107)dt. It follows from the

above arguments that V is bounded. Therefore,

˙V = p(cid:62)R ˙p + ˙p(cid:62)Rp − 2nr(cid:107) ˜L(cid:107)

= p(cid:62)(R(B − D − L∗) + (B − D − L∗)(cid:62)R)p
− 2p(cid:62)R(L(x) − L∗)p − 2p(cid:62)RP Bp − 2nr(cid:107) ˜L(cid:107)
= −p(cid:62)Kp − 2p(cid:62)R ˜L(x)p − 2p(cid:62)RP Bp
− 2nr(cid:107) ˜L(cid:107)
≤ −p(cid:62)Kp + 2nr(cid:107) ˜L(cid:107) − 2nr(cid:107) ˜L(cid:107) − 2p(cid:62)RP Bp
≤ −2p(cid:62)RP Bp ≤ 0.

9

(2.7)

Since all the signals and their derivatives are bounded, it follows that ¨V (t) is bounded and

hence ˙V is uniformly continuous in t. Therefore from Barbalat’s lemma and its application
to Lyapunov-like functions [28, Lemma 4.3, Chapter 4] it follows that ˙V → 0 as t → ∞.
Consequently, from (3.8), p(cid:62)RP Bp → 0. Since R > 0, B > 0 and pi ≥ 0, p(t) → 0 as
t → ∞. This establishes global attractivity of the disease-free equilibrium point. We now
establish its stability.

We note that (cid:107) ˜L(x)(cid:107) is a real analytic function of x, for x (cid:29) 0. Therefore, there exists
a region (cid:107)x − x∗(cid:107) < δ1 in which (cid:107) ˜L(x)(cid:107) ≤ k1(cid:107)x − x∗(cid:107) for some k1 > 0. Also, since
x − x∗ is globally exponentially stable, (cid:107)x(t) − x∗(cid:107) ≤ k2e−αt(cid:107)x(0) − x∗(cid:107) for some k2,
(cid:82) t
, then (cid:107) ˜L(x)(cid:107) ≤ k1k2e−αt(cid:107)x(0) − x∗(cid:107). This implies
α > 0. Thus, if (cid:107)x(0) − x∗(cid:107) <
0 (cid:107) ˜L(cid:107)dt ≤ k

α(cid:107)x(0) − x∗(cid:107), where k := k1k2. Now, since ˙V (p, t) ≤ 0,

δ1
k2

V (p(0), 0) = p(0)(cid:62)Rp(0)

≥ V (p(t), t)
≥ p(t)(cid:62)Rp(t) − 2
≥ Rmin(cid:107)p(t)(cid:107)2 − 2

nrk(cid:107)x(0) − x∗(cid:107)
nrk(cid:107)x(0) − x∗(cid:107)

α

,

α

where Rmin = mini(Ri). Equivalently,

(cid:107)p(t)(cid:107)2 ≤ r
Rmin

(cid:107)p(0)(cid:107)2 + 2

nrk(cid:107)x(0) − x∗(cid:107)

αRmin

.

It follows using stability of x dynamics, that for any 
 > 0, there exists δ > 0 , such that
(cid:107)x(0) − x∗(cid:107)2 + (cid:107)p(0)(cid:107)2 ≤ δ2 ⇒ (cid:107)p(t)(cid:107)2 + (cid:107)x(t) − x∗(cid:107)2 ≤ 
2. This establishes stability.
Together, global attractivity and stability prove the fourth statement.

Stability of endemic equilibria: Finally, we prove the fifth statement. To this end, we

first establish an intermediate result.

Lemma 1 For the dynamics (3.1a), if pi(t) → 0 as t → ∞, for some i ∈ {1, . . . , n}, then
p(t) → 0 as t → ∞.

10

Proof: The dynamics of pi are

˙pi = (βi − δi − lii(x))pi −(cid:88)

j(cid:54)=i

lij(x)pj − βip2
i .

(2.8)

It can be easily seen that ¨pi is bounded and hence ˙pi is uniformly continuous in t. Now
if pi(t) → 0 as t → ∞, it follows from Barbalat’s lemma [28, Lemma 4.2] that ˙pi → 0.
Therefore, from (2.8) and the fact that −lij(x) ≥ 0 and pi ≥ 0, it follows that pj(t) → 0 for
all j such that −lij(x) (cid:54)= 0. Using Assumption 1 and applying the above argument at each
node implies p(t) → 0.
(cid:3)

Define ˜p := p − p∗, P∗ := diag(p∗) and ˜P := diag(˜p). Then

˙˜p = (B − D − L(x) − P B)p
= (B − D − L∗ − P∗B)p∗ + (B − D − L∗ − P∗B)˜p
− ˜L(x)p − ˜P Bp
= (B − D − L∗ − P∗B)˜p − ˜L(x)p − ˜P Bp.

where (B − D − L∗ − P∗B)p∗ = 0, as (p∗, x∗) is an equilibrium point.

Note that (B − D − L∗ − P∗B) is an irreducible Metzler matrix. The Perron-Frobenius
theorem for irreducible Metzler matrices [27] implies µ(B − D − L∗ − P∗B) = 0 and the
associated eigenvector p∗ (cid:29) 0n. Also, this means there exists a positive-diagonal matrix R2
and a positive semi-definite matrix K2 such that

R2(B − D − L∗ − P∗B) + (B − D − L∗ − P∗B)(cid:62)R2 = −K2.

Similar to the proof of the fourth statement, take V2(˜p, t) = ˜p(cid:62)R2˜p − 2nr2

(cid:82) t
0 (cid:107) ˜L(cid:107)dt,

11

where r2 := (cid:107)R2(cid:107). Then,

˙V2 = ˜p(cid:62)R2 ˙˜p + ˙˜p(cid:62)R2˜p − 2nr2(cid:107) ˜L(cid:107)
= ˜p(cid:62)(R2(B − D − L∗ − P∗B)
+ (B − D − L∗ − P∗B)(cid:62)R2)˜p
− 2˜p(cid:62)R2 ˜L(x)p − 2˜p(cid:62)R2 ˜P Bp − 2nr2(cid:107) ˜L(cid:107)
= −˜p(cid:62)K2˜p − 2˜p(cid:62)R2 ˜L(x)p − 2˜p(cid:62)R2 ˜P Bp − 2nr2(cid:107) ˜L(cid:107)
≤ −˜p(cid:62)K2˜p + 2nr2(cid:107) ˜L(cid:107) − 2nr2(cid:107) ˜L(cid:107) − 2˜p(cid:62)R2 ˜P Bp
≤ −2˜p(cid:62)R2 ˜P Bp = −2

(R2)iβi ˜p2

n(cid:88)

i pi ≤ 0.

i=1

It can be easily shown that ¨V2 is bounded implying ˙V2 is uniformly continuous. Applying
Barbalat’s lemma [28, Lemma 4.2] gives ˙V2 → 0 as t → ∞. Now, since R2 and B are positive
diagonal matrices this implies that ˜pipi → 0, for each i. Using Lemma 1, and the fact that
p = 0 is an unstable equilibrium for µ(B − D− L∗) > 0, we have ˜p → 0 as long as p(0) (cid:54)= 0.
Stability can be established similarly to the disease-free equilibrium case. This concludes
(cid:3)

the proof of the theorem.

Corollary 1 (Stability of disease-free equilibria) For the SIS epidemic model under
Markovian mobility (2.1) with Assumption 1 and the disease-free equilibrium (p∗, x∗) =
(0n, v) the following statements hold

i. a necessary condition for stability is δi > βi − νi, for each i ∈ {1, . . . , n};
ii. a necessary condition for stability is that there exists some i ∈ {1, . . . , n} such that

δi ≥ βi;

iii. a sufficient condition for stability is δi ≥ βi, for each i ∈ {1, . . . , n};

iv. a sufficient condition for stability is

(cid:16)

(cid:114)

1 +

1 +

λ2

(cid:0)δi−βi−m(cid:1)(cid:17)2

λ2

(cid:80)

i wi

+ m ≥ 0,

n + 1

12

where w is a positive left eigenvector of L∗ such that w(cid:62)L∗ = 0 with maxi wi = 1,
m = mini(δi − βi), W = diag(w), and λ2 is the the second smallest eigenvalue of
2(W L∗ + L∗(cid:62)W ).

1

Proof: We begin by proving the first two statements. First, we note that L∗
ii = νi. This
can be verified by evaluating L∗ = L(v) and utilising the fact that Q(cid:62)v = 0. The necessary
and sufficient condition for the stability of disease-free equilibrium is µ(B − D − L∗) ≤ 0.
Since, B − D − L∗ is an irreducible Metzler matrix, a necessary condition for µ ≤ 0 is that
its diagonal terms are strictly negative, i.e., βi − δi − νi < 0, for each i ∈ {1, . . . , n}. This
gives the statement (i).

Perron-Frobenius theorem for irreducible Metzler matrices implies that there exists a real
eigenvalue equal to µ with positive eigenvector, i.e., (B − D − L∗)y = µy, where y (cid:29) 0n.
Since, µ ≤ 0, written component-wise for i∗, where yi∗ = min(yi) :

(cid:88)

j(cid:54)=i∗

lij(yj − yi∗)

(βi∗ − δi∗ − νi∗)yi∗ − (cid:88)
⇒ (βi∗ − δi∗)yi∗ ≤ (cid:88)

j(cid:54)=i∗
⇒ (βi∗ − δi∗)yi∗ ≤ (νi∗ +

j(cid:54)=i∗

lijyj ≤ 0
(cid:88)

lij)yi∗ +

j(cid:54)=i∗
lij(yj − yi∗)

⇒ βi∗ − δi∗ ≤ 0.

This proves the statement (ii).

Since, L∗ is a Laplacian matrix, if δi ≥ βi, for each i ∈ {1, . . . , n}, from Gershgorin disks

theorem [27], µ ≤ 0, which proves the third statement.

For the last statement, we use an eigenvalue bound for perturbed irreducible Laplacian

matrix of a digraph [29, Theorem 6], stated below:

Let H = L + ∆, where L is an n × n irreducible Laplacian matrix and ∆ (cid:54)= 0 is a

13

non-negative diagonal matrix, then
Re(λ(H)) ≥

(cid:16)

(cid:114)

1 +

1 +

(cid:17)2

> 0,

n + 1

λ2

λ2(cid:80)

i wi∆i

where, w is a positive left eigenvector of L such that w(cid:62)L = 0 with maxi wi = 1, W =
diag(w), and λ2 is the second smallest eigenvalue of 1
Now, in our case necessary and sufficient condition for stability of disease-free equilibrium

2(W L + L(cid:62)W ).

is:

Re(λ(L∗ + D − B)) = Re(λ(L∗ + ∆ + mI))

= Re(λ(L∗ + ∆)) + m ≥ 0

where, m = mini(δi − βi) and ∆ = D − B − mI. Applying the eigenvalue bound with
H = L∗ + ∆ gives the sufficient condition (iv).
(cid:3)

Remark 2 It can be shown that v is the left eigenvector associated with eigenvalue zero for
both Q and L∗, i.e., v(cid:62)Q = v(cid:62)L∗ = 0 and thus can be re-scaled to compute w =
maxi(vi) v.
(cid:3)

1

Remark 3 For a given graph and the associated mobility transition rates in dynamics (2.1),
let m = mini(δi − βi) and i∗ = argmini(δi − βi). Then, there exist δi’s, i (cid:54)= i∗, that satisfy
statement (iv) of Corollary 1 if m > mlower, where
mlower = − λ2
4n + 1

.

(cid:3)

Remark 4 (Influence of mobility on stability of disease-free equilibrium.) The

statement (iv) of Corollary 1 characterizes the influence of mobility on the stability of disease-

free equilibria. In particular, λ2 is a measure of “intensity" of mobility and m is a measure of
largest deficit in the recovery rate compared with infection rate among nodes. The sufficient

condition in statement (iv) states explicitly how mobility can allow for stability of disease-free
(cid:3)

equilibrium even under deficit in recovery rate at some nodes.

14

2.3 Numerical Illustrations

We start with numerical simulation of epidemic model with mobility in which we treat

epidemic spread as well as mobility as stochastic processes. We take 20 simulations with

same initial conditions and parameters and take the average of the results. The fraction of

infected populations for different cases are shown in Fig. 2.1. We take a line graph and the

mobility transition rates being equal among out going neighbors of a node. The two cases

relate to the stable disease-free equilibrium and stable endemic equilibrium respectively. We

have chosen heterogeneous curing or infection rates to elucidate the influence of mobility.

If the curing rates, infection rates and the initial fraction of infected population is same

for all the nodes, mobility does not play any role. The corresponding simulations of the

deterministic model as per Proposition 1 are also shown for comparison. Figure 2.1 (a)
corresponds to the case δi ≥ βi for each i, whereas Fig. 2.1 (c) corresponds to the case
δi < βi for each i. The results support statements (iii) and (ii) of Corollary 1 and lead to,
respectively, the stable disease-free equilibrium and the stable endemic equilibrium.

(a) Stable disease-free
equilibrium: Stochastic
model

(b) Stable disease-free
equilibrium: Determin-
istic model

Stable

endemic
(c)
equilibrium: Stochastic
model

(d) Stable endemic equi-
librium: Deterministic
model

Figure 2.1: Stochastic simulation of epidemic spread under mobility. Line graph, n = 20,
ν(i) = 0.2, qij = ν(i)
, pi(0) = 0.01. Each iteration in stochastic model corresponds to
Dout

time-step 0.01 sec.

Once we have established the correctness of deterministic model predictions with the

stochastic simulations, we study the simulations of deterministic model only. We study the

effect of mobility over 4 different mobility graph structure - line graph, ring graph, star graph

and a complete graph. First we keep the equilibrium distribution of population same for all

15

the four graphs by using instantaneous transition rates from Metropolis-Hastings algorithm

[30]. This shows the effect of different mobility graph structure on epidemic spread while

the equilibrium population distribution remains the same. Fig. 2.2 shows the fractions of

infected population trajectories for 20 nodes connected with 4 different graph structures. The

nodes have heterogeneous curing rates and these rates are the same across different graph

structures. The values of equilibrium fractions are affected by the presence of mobility and

are different for different graph structures. As seen in Fig. 2.2, star graph has the widest

distribution of equilibrium infected fraction values whereas complete graph has the narrowest

of the four.

(a) Line graph

(b) Ring graph

(c) Star graph

(d) Complete graph

Figure 2.2: Simulation of deterministic model of epidemic spread under mobility, with same
equilibrium distribution of population over 4 different graph structure with stable endemic

equilibrium. n = 20, pi(0) = 0.01.

Next, we verify the statement (iv) of Corollary 1, where one can have some curing rates

δi less than the infection rates βi but still have stable disease-free equilibrium. We take
a complete graph of n = 20 nodes with given mobility transition rates which give us w,
L∗ and λ2. We take a given set of values of βi. Next, we compute mlower = − λ2
4n+1
and take 0.8 times of this value as m in order to compute δi’s that satisfy statement (iv)
of Corollary 1. For our case the values are: βi = 0.3, λ2 = 0.2105, mlower = −0.0026,
m = 0.8 mlower = −0.0021, δ1 = δn = βi + m and the rest δi computed to satisfy the
condition which gives δ1 = δn = 0.2979 and δi = 0.3198 for i ∈ {2, . . . , n− 1}. Fig. 2.3 shows
the trajectories of infected fraction populations. As can be seen the trajectories converge to

the disease-free equilibrium.

16

Figure 2.3: Stable disease-free equilibrium with curing rates computed as per the λ2

sufficient condition (statement (iv), Corollary 1) for stability of disease-free equilibrium.

Graph: Complete, n = 20, pi(0) = 0.01.

17

EPIDEMIC MODELS UNDER MOBILITY ON MULTI-LAYER NETWORKS

CHAPTER 3

In this chapter, we study SIS and SIR models under mobility on multi-layer networks. The

presence of multi-layer arises from the presence of different mobility patterns with their own

Markov chains and digraphs. We give complete characterization for existence, uniqueness

and stability of equilibria. We give numerical illustrations to support our results.

3.1 SIS Model under Multi-layer Markovian Mobility

We consider n sub-population of individuals that are located in distinct spatial regions

(patches). We assume the individuals within each patch can be classified into two categories:

(i) susceptible, and (ii) infected. We assume that the individuals within each patch are further

grouped into m classes which decide how they travel to other patches. Let the connectivity
of these patches corresponding to the mobility pattern of each class α ∈ {1, . . . , m} be
modeled by a digraph Gα = (V,E α), where V = {1, . . . , n} is the node (patch) set and
E α ⊂ V × V is the edge set. We model the mobility of individuals on each graph Gα using a
Continuous Time Markov Chain (CTMC) with generator matrix Qα, whose (i, j)-th entry
ij ≥ 0, i (cid:54)= j, is the instantaneous transition rate from node i to node
is the total rate of transition out of node i, i.e., να
ij. Here,
i (t) be the number of individuals of
i ) be the fraction of infected
n](cid:62),

is qα
j, and −qα
ij > 0, if (i, j) ∈ E α; and qα
qα
class α in patch i at time t. Let pα

i ∈ [0, 1] (respectively, 1 − pα

ij = 0, otherwise. Let xα

i = (cid:80)

(respectively, susceptible) sub-population of class α at patch i. Define pα := [pα
n](cid:62), p := [(p1)(cid:62), . . . , (pm)(cid:62)](cid:62) and x := [(x1)(cid:62), . . . , (xm)(cid:62)](cid:62).

xα := [xα

1 , . . . , xα

1 , . . . , pα

ij. The entry qα

ii = να
i

j(cid:54)=i qα

We model the interaction of mobility with the epidemic process as follows. At each time t,
individuals of each class α within each node move on graph Gα according to the CTMC with
generator matrix Qα and then interact with individuals within their current node according
to an SIS epidemic process. For the epidemic process at node i, let βi > 0 and δi ≥ 0 be

18

the infection and recovery rate, respectively. We let Bα > 0 and Dα ≥ 0 be the positive
and non-negative diagonal matrices with entries βi and δi, i ∈ {1, . . . , n}, respectively. Let
B and D be the positive and non-negative diagonal matrices with block-diagonal entries Bα
and Dα, α ∈ {1, . . . , m}, respectively. Let P α := diag(pα) and P := diag(p). We now
derive the continuous time dynamics that captures the interaction of mobility and the SIS

epidemic dynamics.

Proposition 2 (SIS model under mobility) The dynamics of the fractions of the in-

fected sub-population p and the number of individuals xα under multi-layer Markovian mo-

bility model with generator matrices Qα, and infection and recovery matrices B and D,

respectively, are

˙p = (BF (x) − D − L(x))p − P BF (x)p
˙xα = (Qα)(cid:62)xα,

(3.1a)

(3.1b)

where L is an nm× nm block-diagonal matrix with block-diagonal terms Lα, α ∈ {1, . . . , m},
Lα(x) is a matrix with entries

(cid:80)



lα
ij(x) =

xα
j
xα
i

,

if i = j,

otherwise,

j(cid:54)=i qα
ji
xα
j
xα
i

ji

,

−qα

F (x) := [ ¯F(cid:62)(x), . . . , ¯F(cid:62)(x)](cid:62) be a row-concatenated nm×nm matrix with each n×nm block-
row as ¯F (x) := [F 1(x), . . . , F m(x)], and F α as a diagonal matrix with entries f α
i(cid:80)

, i.e., the fraction of total population at node i contributed by class α.

i (x) :=

xα
α xα
i

Proof: Consider a small time increment h > 0 at time t. Then the number of individuals

of class α present at node i after the evolution of CTMC in time-interval [t, t + h) is

(cid:88)

j(cid:54)=i

xα
i (t + h) = xα

i (t)(1 − να

i h) +

19

qα
jixα

j (t)h + o(h).

(3.2)

After the mobility, individuals within each node interact according to SIS dynamics. Thus,

the fraction of infected population present at node i is:

i (t + h)

i (t + h)pα
xα
= −xα

i (t)δipα
i (t)pα

i (t)h + xα
i (t)(1 − να

(cid:88)

i (t)βi ¯pi(t)(1 − pα
j xα
i h) +

qα
jipα

+ xα

i (t))h

j (t)h + o(h).

where ¯pi is the fraction of infected population at node i and is given as:

j(cid:54)=i

(cid:88)

¯pi :=

i pα
f α
i .

(3.3)

α

The first two terms on the right side of (3.3) correspond to epidemic process within each

node, whereas the last two terms correspond to infected individuals coming from other nodes
from (3.2) in (3.3) and taking the limit h → 0+

due to mobility. Using the expression of xα
i

gives

˙pα
i

=

−δipα

i

+ βi ¯pi(1 − pα

i ) − lα

iipα
i

− (cid:88)

j(cid:54)=i

ijpα
lα
j .

(3.4)

Writing above in vector form gives:

˙pα = (−Dα − Lα(xα))pα + Bα ¯F (x)p − P αBα ¯F (x)p.

Similarly taking limits in (3.2) yields

i = −να
˙xα

i xα

i +

(cid:88)

j(cid:54)=i

qα
jixα
j .

Rewriting (3.4) and (3.6) in vector form establishes the proposition.

(3.5)

(3.6)

(cid:3)

3.2 Analysis of SIS Model under Multi-layer Markovian Mobility

In this section, we analyze the SIS model under multi-layer mobility (3.1) under the

following standard assumption:

Assumption 2 Digraph Gα is strongly connected, for all α ∈ {1, . . . , m}, which is equivalent
(cid:3)
to matrices Qα being irreducible [27].

20

Let vα be the right eigenvector of (Qα)(cid:62) associated with eigenvalue at 0. We assume
that vα is scaled such that its inner product with the associated left eigenvector 1n is
n vα = 1. Define v := [N 1(v1)(cid:62), . . . , N m(vm)(cid:62)](cid:62), where N α is the total
unity, i.e., 1(cid:62)
number of individuals belonging to class α, for α ∈ {1, . . . , m}. We call an equilibrium point
(p∗, x∗), an endemic equilibrium point, if at equilibrium the disease does not die out, i.e.,
p∗ (cid:54)= 0, otherwise, we call it a disease-free equilibrium point. Let F∗ := F (x∗) = F (v) and
L∗ := L(x∗) = L(v). It can be verified that F∗ admits the splitting, F∗ = I − M, where I
is the identity matrix of appropriate dimensions and M is a Laplacian matrix.

Theorem 2 (Existence and Stability of Equilibria) For the SIS model under multi-

layer Markovian mobility (3.1) with Assumption 2, the following statements hold

i. if p(0) ∈ [0, 1]nm, then p(t) ∈ [0, 1]nm for all t > 0. Also, if p(0) > 0nm, then

p(t) (cid:29) 0nm for all t > 0;

ii. the model admits a disease-free equilibrium at (p∗, x∗) = (0nm, v);
iii. the model admits an endemic equilibrium at (p∗, x∗) = (¯p, v), ¯p (cid:29) 0, if and only if

µ(BF∗ − D − L∗) > 0;

iv. the disease-free equilibrium is globally asymptotically stable if and only if µ(BF∗− D−

L∗) ≤ 0 and is unstable otherwise;

v. the endemic equilibrium is almost globally asymptotically stable if µ(BF∗−D−L∗) > 0

with region of attraction p(0) ∈ [0, 1]nm such that p(0) (cid:54)= 0nm.

Proof: The first part of statement (i) follows from the fact that ˙p is either tangent or
i = 0 or 1 .

directed inside of the region [0, 1]nm at its boundary which are surfaces with pα

This can be seen from (3.4). For the second part of (i), we rewrite (3.1a) as:

˙p = ((I − P )BF (x) + A(x))p − E(t)p

21

where L(x) = C(x) − A(x) with C(x) composed of the diagonal terms of L(x), A(x) is
the non-negative matrix corresponding to the off-diagonal terms, and E(t) = C(x(t)) + D is

a diagonal matrix. Now, consider a variable change y(t) := e

0 E(τ )dτ p(t). Differentiating

(cid:82) t

y(t) and using above gives:

(cid:82) t
(cid:82) t
(cid:82) t
0 −E(τ )dτ e
0 E(τ )dτ ((I − P )BF (x) + A(x))e
(cid:82) t
(cid:82) t
0 −E(τ )dτ y
0 E(τ )dτ ((I − P )BF (x) + A(x))e

˙y = e

= e

0 E(τ )dτ p

Now, the coefficient matrix of y above is always non-negative and strongly connected. The

rest of the proof is same as in [8, Theorem 4.2 (i)].

The second statement follows by inspection.

The proof of the third statement is presented in Appendix A.2.

Stability of disease-free equilibria: To prove the fourth statement, we begin by estab-

lishing sufficient conditions for instability. The linearization of (3.1) at (p, x) = (0, v) is

 ˙p

 =

˙x

BF∗ − D − L∗

0



p

 .

x

0
Q(cid:62)

(3.7)

Since the system matrix in (3.7) is block-diagonal, its eigenvalues are the eigenvalues of
the block-diagonal sub-matrices. Further, since radial abscissa µ(Q(cid:62)) is zero, a sufficient
condition for instability of the disease-free equilibrium is that µ(BF∗ − D − L∗) > 0.

For the case of µ(BF∗ − D − L∗) ≤ 0, we now show that the disease-free equilibrium is
a globally asymptotically stable equilibrium. It can be seen from the definitions of matrices
F∗ and L∗, that under Assumption 2, (BF∗ − D − L∗) is an irreducible Metzler matrix.
Together with µ(BF∗ − D − L∗) ≤ 0, implies there exists a positive diagonal matrix R such
that

R(BF∗ − D − L∗) + (BF∗ − D − L∗)(cid:62)R = −K,

where K is a positive semi-definite matrix [10, Proposition 1 (iv), Lemma A.1]. Define
˜L := L(x) − L∗, ˜F := F (x) − F∗ and r := (cid:107)R(cid:107), where (cid:107) · (cid:107) denotes the the induced two
norm of the matrix.

22

Since x(0) (cid:29) 0, under Assumption 2, xα

i (t) is lower bounded by some positive con-
stant and hence, ˜L and ˜F are bounded and continuously differentiable. Since x is bounded
and exponentially converges to x∗, it follows that (cid:107) ˜L(x)(cid:107) and (cid:107) ˜F (x)(cid:107) locally exponentially

converge to 0 and(cid:82) t
Consider the Lyapunov-like function V (p, t) = p(cid:62)Rp − 2nmr(cid:82) t

0 (cid:107) ˜F(cid:107)dt are bounded for all t > 0.

0 (cid:107) ˜L(cid:107)dt and(cid:82) t

0 (B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))dt.

It

follows from the above arguments that V is bounded. Therefore,

˙V = 2p(cid:62)R ˙p − 2nmr(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))

= p(cid:62)(R(BF∗ − D − L∗) + (BF∗ − D − L∗)(cid:62)R)p
+ 2p(cid:62)R(B ˜F − ˜L)p − 2p(cid:62)RP BF p
− 2nmr(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
= −p(cid:62)Kp + 2p(cid:62)R(B ˜F − ˜L)p − 2p(cid:62)RP BF p
− 2nmr(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
≤ −p(cid:62)Kp + 2nmr(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
− 2nmr(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107)) − 2p(cid:62)RP BF p
≤ −2p(cid:62)RP BF p ≤ 0.

(3.8)

Since all the signals and their derivatives are bounded, it follows that ¨V (t) is bounded and

hence ˙V is uniformly continuous in t. Therefore from Barbalat’s lemma and its application
to Lyapunov-like functions [28, Lemma 4.3, Chapter 4] it follows that ˙V → 0 as t → ∞.
Consequently, from (3.8), p(cid:62)RP BF p → 0. Since R > 0, B > 0, F ≥ 0 with Fkk > 0
and pk ≥ 0, p(t) → 0 as t → ∞. This establishes global attractivity of the disease-free
equilibrium point. We now establish its stability.

We note that since, for x (cid:29) 0, (B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107)) is a real analytic function of x, ∃ a
region (cid:107)x − x∗(cid:107) < δ1 in which (B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107)) ≤ k1(cid:107)x − x∗(cid:107) for some k1 > 0. Also, since
x − x∗ is globally exponentially stable, (cid:107)x(t) − x∗(cid:107) ≤ k2e−αt(cid:107)x(0) − x∗(cid:107) for k2, α > 0.
, then (B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107)) ≤ k1k2e−αt(cid:107)x(0) − x∗(cid:107). This implies
Thus, if (cid:107)x(0) − x∗(cid:107) <

δ1
k2

23

(cid:82) t
0 (B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))dt ≤ k

α(cid:107)x(0) − x∗(cid:107), where k := k1k2. Now, since ˙V (p, t) ≤ 0,
V (p(0), 0) = p(0)(cid:62)Rp(0)

≥ V (p(t), t)
≥ p(t)(cid:62)Rp(t) − 2
≥ Rmin(cid:107)p(t)(cid:107)2 − 2

nmrk(cid:107)x(0) − x∗(cid:107)
nmrk(cid:107)x(0) − x∗(cid:107)

α

α

,

where Rmin = mini(Ri). Equivalently,

(cid:107)p(t)(cid:107)2 ≤ r
Rmin

(cid:107)p(0)(cid:107)2 + 2

nmrk(cid:107)x(0) − x∗(cid:107)

αRmin

.

It follows using stability of x dynamics, that for any 
 > 0, there exists δ > 0 , such that
(cid:107)x(0) − x∗(cid:107)2 + (cid:107)p(0)(cid:107)2 ≤ δ2 ⇒ (cid:107)p(t)(cid:107)2 + (cid:107)x(t) − x∗(cid:107)2 ≤ 
2. This establishes stability.
Together, global attractivity and stability prove the fourth statement.

Stability of endemic equilibria: Finally, we prove the fifth statement. To this end, we

first establish an intermediate result.

Lemma 2 For the dynamics (3.1a), if pα
α ∈ {1, . . . , m}, then p(t) → 0 as t → ∞.

i (t) → 0 as t → ∞, for some i ∈ {1, . . . , n} and

Proof:

It can be easily seen from (3.4) that ¨pα
i

continuous in t. Now if pα

is bounded and hence ˙pα
i

is uniformly
i (t) → 0 as t → ∞, it follows from Barbalat’s lemma [28, Lemma
i ≥ 0, it follows
ij(x) (cid:54)= 0. Using Assumption 2 and applying the above
(cid:3)

ij(x) ≥ 0 and pα

4.2] that ˙pα

i → 0. Therefore, from (3.4) and the fact that −lα

j (t) → 0 for all j such that −lα

that pα
argument for each class at each node implies p(t) → 0.

24

Define ˜p := p − p∗, P∗ := diag(p∗) and ˜P := diag(˜p). Then

˙˜p = (BF − D − L − P BF )p
= (BF∗ − D − L∗ − P∗BF∗)p∗
+ (BF∗ − D − L∗ − P∗BF∗)˜p
+ (B ˜F − ˜L)p − P B ˜F p − ˜P BF∗p
= ((I − P∗)BF∗ − D − L∗)˜p + ((I − P )B ˜F − ˜L)p
− ˜P BF∗p.

where we have used (BF∗ − D − L∗ − P∗BF∗)p∗ = 0, as (p∗, x∗) is an equilibrium point.
Note that (BF∗− D− L∗− P∗BF∗) = ((I − P∗)BF∗− D− L∗) is an irreducible Metzler
matrix and p∗ (cid:29) 0 is its positive eigenvector associated with eigenvalue at zero. Therefore,
the Perron-Frobenius theorem for irreducible Metzler matrices [27] implies µ((I − P∗)BF∗−
D − L∗) = 0. Also, this means there exists a positive-diagonal matrix R2 and a positive
semi-definite matrix K2 such that

R2((I − P∗)BF∗ − D − L∗)

+ ((I − P∗)BF∗ − D − L∗)(cid:62)R2 = −K2.

Similar to the proof of the fourth statement, take V2(˜p, t) = ˜p(cid:62)R2˜p− 2nmr2

(cid:82) t
0 (B(cid:107) ˜F(cid:107) +

25

(cid:107) ˜L(cid:107))dt, where r2 := (cid:107)R2(cid:107). Then,

˙V2 = 2˜p(cid:62)R2 ˙˜p − 2nmr2(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
= ˜p(cid:62)(R2((I − P∗)BF∗ − D − L∗)
+ ((I − P∗)BF∗ − D − L∗)(cid:62)R2)˜p
+ 2˜p(cid:62)R2((I − P )B ˜F − ˜L)p − 2˜p(cid:62)R2 ˜P BF∗p
− 2nmr2(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
= −˜p(cid:62)K2˜p + 2˜p(cid:62)R2((I − P )B ˜F − ˜L)p
− 2˜p(cid:62)R2 ˜P BF∗p − 2nmr2(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
≤ −˜p(cid:62)K2˜p + 2nmr2(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
− 2nmr2(B(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107)) − 2˜p(cid:62)R2 ˜P BF∗p
≤ −2˜p(cid:62)R2 ˜P BF∗p
nm(cid:88)
(R2)kβkF∗
≤ −2

kpk ≤ 0.

kk ˜p2

k=1

The last inequality above follows from the fact that R2 > 0, B > 0 and F ≥ 0 with
diagonal terms Fkk > 0. It can be easily shown that ¨V2 is bounded implying ˙V2 is uniformly
continuous. Applying Barbalat’s lemma [28, Lemma 4.2] gives ˙V2 → 0 as t → ∞. This
implies that ˜pkpk → 0. Using Lemma 2, and the fact that p = 0 is an unstable equilibrium
for µ(BF∗ − D − L∗) > 0, we have ˜p → 0 as long as p(0) (cid:54)= 0. Stability can be established
similarly to the disease-free equilibrium case. This concludes the proof of the theorem. (cid:3)

Corollary 2 (Stability of disease-free equilibria) For the SIS epidemic model under

Multi-layer Markovian mobility (3.1) with Assumption 2 and the disease-free equilibrium
(p∗, x∗) = (0, v) the following statements hold

i. a necessary condition for stability is that for each i ∈ {1, . . . , n}, ∃α ∈ {1, . . . , m} such

that δi > βi − να
i ;

26

ii. a necessary condition for stability is that there exists some i ∈ {1, . . . , n} such that

δi ≥ βi;

iii. a sufficient condition for stability is δi ≥ βi, for each i ∈ {1, . . . , n};

iv. a sufficient condition for stability is

(cid:114)

(cid:16)

1 +

(cid:0)δi−βi−s(cid:1)(cid:17)2

λ2
λ2

+ s ≥ 0,

nm + 1

(cid:80)

i wi

1 +

where w is a positive left eigenvector of (BM + L∗) such that w(cid:62)(BM + L∗) = 0 with
maxi wi = 1, s = mini(δi−βi), W = diag(w), and λ2 is the second smallest eigenvalue
of 1
2

(cid:0)W (BM + L∗) + (BM + L∗)(cid:62)W(cid:1).

ii = να

Proof: We begin by proving the first two statements. First, we note that (Lα)∗
i . This
can be verified by evaluating L∗ = L(v) and utilising the fact that Q(cid:62)v = 0. The necessary
and sufficient condition for the stability of disease-free equilibrium is µ(BF∗ − D − L∗) ≤ 0.
Note that BF∗ − D − L∗ is an irreducible Metzler matrix. Perron-Frobenius theorem for
irreducible Metzler matrices implies that there exists a real eigenvalue equal to µ with positive
eigenvector, i.e., (BF∗ − D − L∗)y = µy, where y (cid:29) 0. Rename components of y as
y(nα+i) = yα

i to write y = [(y1)(cid:62), . . . , (ym)(cid:62)](cid:62).

Let for each i ∈ {1, . . . , n}, y

i = min{y1
ki

i , . . . , ym

i }. Since µ ≤ 0, written component-wise

27

for (nki + i)-th component(cid:88)
βif∗α
⇒(cid:88)
−(cid:88)

α

α

i yα

βif∗α
ki
i +
i y
∗ki
j ≤ 0
ki
l
ij y

α

i −(cid:88)

ki

ki
i )y

i − (δi + ν
(cid:88)

j(cid:54)=i
i − y
(yα

βif∗α

i

∗ki
l
ij y

j ≤ 0
ki

i ) − (δi + ν
ki

ki
ki
i )y
i

j(cid:54)=i

⇒ (βi − δi − ν

ki
i )y

i ≤ −(cid:88)

ki

α

βif∗α

i

i − y
(yα

ki
i ) +

(cid:88)

j(cid:54)=i

∗ki
l
ij y

ki
j

∗ki
ij ≤ 0 and that there exists j ∈ {1, . . . , n}

n }. Similar to the proof of the first statement

Here we have used facts: (cid:80)

ki
i < 0

⇒ (βi − δi − νk
i )y
⇒ βi − δi − ν
ki
i < 0.
i = 1, f∗α

α f∗α

Let y

1, . . . , ym

i > 0, l
∗ki
ij < 0. This proves the statement (i).
such that l
i be min{y1
ki
(cid:88)
⇒(cid:88)
−(cid:88)

i −(cid:88)

βif∗α
∗ki
ij (y
l

βif∗α
i y
∗ki
l
ij y

j − y
ki

βif∗α

i yα

ki
i )y

ki

α

α

i

j(cid:54)=i
i − y
(yα
i ) ≤ 0
ki

∗ki
l
ij y

ki
i +

i − (δi + ν
(cid:88)
i −(cid:88)
i ≤ −(cid:88)

j(cid:54)=i

ki

ki

α

α

j(cid:54)=i

⇒ (βi − δi)y

ki
ki
i )y
i

j ≤ 0
ki

i ) − (δi + ν
ki
(cid:88)

ki
i ) +

j(cid:54)=i

βif∗α

i

i − y
(yα

∗ki
ij (y
l

j − y
ki

ki
i )

⇒ (βi − δi)y
i ≤ 0
ki
⇒ βi − δi ≤ 0.

Here we have used an additional fact: ν

(cid:88)

j(cid:54)=i

ki
i +

∗ki
ij = 0. This proves statement (ii).
l

Let F∗ = I − M where M is a Laplacian matrix which can be seen from the definition of
F . Now BF∗ − D− L∗ = B − D− (BM + L∗). Since (BM + L∗) is an irreducible Laplacian

28

matrix, if δi ≥ βi, for each i ∈ {1, . . . , n}, from Gershgorin disks theorem [27], µ ≤ 0, which
proves the third statement.

For the last statement, we use an eigenvalue bound for perturbed irreducible Laplacian

matrix of a digraph [29, Theorem 6], stated below:

Let H = A + ∆, where A is an n × n irreducible Laplacian matrix and ∆ (cid:54)= 0 is a

non-negative diagonal matrix, then
Re(λ(H)) ≥

(cid:16)

(cid:114)

1 +

1 +

(cid:17)2

> 0,

n + 1

λ2

λ2(cid:80)

i wi∆i

where, w is a positive left eigenvector of A such that w(cid:62)A = 0 with maxi wi = 1, W =
diag(w), and λ2 is the second smallest eigenvalue of 1
Now, in our case necessary and sufficient condition for stability of disease-free equilibrium

2(W A + A(cid:62)W ).

is:

Re(λ(BM + L∗ + D − B)) = Re(λ(BM + L∗ + ∆ + sI))
= Re(λ(BM + L∗ + ∆)) + s
≥ 0,

where, s = mini(δi − βi) and ∆ = D − B − sI. Applying the eigenvalue bound with
H = BM + L∗ + ∆ gives the sufficient condition (iv).
(cid:3)

Remark 5 For given graphs and the associated mobility transition rates in dynamics (3.1),
let s = mini(δi − βi) and i∗ = argmini(δi − βi). Then, there exist δi’s, i (cid:54)= i∗, that satisfy
statement (iv) of Corollary 3 if s > slower, where

slower = − λ2

4mn + 1

.

(cid:3)

Remark 6 (Influence of mobility on stability of disease-free equilibrium.) The

statement (iv) of Corollary 3 characterizes the influence of mobility on the stability of disease-

free equilibria. In particular, λ2 is a measure of “intensity" of mobility and s is a measure of

29

largest deficit in the recovery rate compared with infection rate among nodes. The sufficient

condition in statement (iv) states explicitly how mobility can allow for stability of disease-free
(cid:3)

equilibrium even under deficit in recovery rate at some nodes.

3.3 Numerical Illustrations

We start with numerical simulation of epidemic model with multi-layer mobility in which

we treat epidemic spread as well as mobility as stochastic processes. The fraction of infected

populations for different cases are shown in Fig. 3.1. The corresponding simulations of

the deterministic model as per Proposition 2 are also shown for comparison. We take two

mobility network layers: a complete graph and a line graph with the mobility transition

rates being equal among out going neighbors of a node for both the graphs. The two cases

relate to the stable disease-free equilibrium and stable endemic equilibrium respectively. If

the curing rates, infection rates and the initial fraction of infected population are the same

for all the nodes, mobility does not play any role. Therefore, we have chosen heterogeneous

curing or infection rates to elucidate the influence of mobility. Figure 3.1 (a) corresponds
to the case δi ≥ βi for each i, whereas Fig. 3.1 (c) corresponds to the case δi < βi for each
i. The results support statements (iii) and (ii) of Corollary 3 and lead to, respectively, the

stable disease-free equilibrium and the stable endemic equilibrium.

(a) Stable disease-free
equilibrium: Stochastic
model

(b) Stable disease-free
equilibrium: Determin-
istic model

Stable

(c)
endemic
equilibrium: Stochastic
model

(d) Stable endemic equi-
librium: Deterministic
model

Figure 3.1: Stochastic simulation of epidemic spread under mobility. Complete-Line
graphs, n = 20, ν(i) = 0.2, qij = ν(i)
, pi(0) = 0.01. Each iteration in stochastic model
Dout

corresponds to time-step 0.01 sec.

30

Once we have established the correctness of deterministic model predictions with the

stochastic simulations, we study the simulations using only the deterministic model. We

study the effect of multi-layer mobility over different pairs of mobility graph structures -

line-line graph, line-ring graph and line-star graph. We choose different population size for

the two mobility layers and take the mobility transition rates such as to keep the equilibrium

distribution of population the same for both the layers across all pairs (taken as uniform

equilibrium distribution) by using instantaneous transition rates from Metropolis-Hastings

algorithm [30]. This shows the effect of different mobility graph structure on epidemic

spread while the equilibrium population distribution remains the same. Fig. 3.2 shows the

fractions of infected population trajectories for 10 nodes connected with different pairs of

graph structures. The values of equilibrium fractions are affected by the presence of mobility

and are different for different graph structures.

(a) Line-Line
graph 1

graphs;

(b) Line-Line
graph 2

graphs;

(c) Line-Ring
graph 1

graphs;

(d) Line-Ring
graph 2

graphs;

(e) Line-Star
graph 1

graphs;

(f) Line-Star
graph 2

graphs;

Figure 3.2: Simulation of deterministic model of epidemic spread under 2 layer mobility,

over different graph structure with stable endemic equilibrium. n = 10, pi(0) = 0.01.

31

3.4 SIR Model under Multi-layer Markovian Mobility

We consider n sub-populations of individuals that are located in distinct spatial regions

(patches). We assume that the individuals within each patch can be classified into three

categories: (i) susceptible, (ii) infected and (iii) recovered. Additionally, we assume that

these individuals are further grouped into m classes depending on how they travel to other

patches. Let the connectivity of these patches corresponding to the mobility pattern of each
class α ∈ {1, . . . , m} be modeled by a digraph Gα = (V,E α), where V = {1, . . . , n} is the
node (patch) set and E α ⊂ V×V is the edge set. We model the mobility of individuals on each
graph Gα using a Continuous Time Markov Chain (CTMC) with generator matrix Qα, whose
ij ≥ 0, i (cid:54)= j, is the instantaneous transition rate from node
(i, j)-th entry is qα
i to node j, and −qα
ij.
j(cid:54)=i qα
i (t) be the number of individuals
Here, qα
i ∈ [0, 1]) be the fraction of
infected (respectively, susceptible) individuals within individuals of class α at patch i. Define
n](cid:62), p := [(p1)(cid:62), . . . , (pm)(cid:62)](cid:62),

is the total rate of transition out of node i, i.e., να

ij > 0, if (i, j) ∈ E α; and qα

i ∈ [0, 1] (respectively, sα

of class α in patch i at time t. Let pα

ij = 0, otherwise. Let xα

i =(cid:80)

1 , . . . , pα

pα := [pα
s := [(s1)(cid:62), . . . , (sm)(cid:62)](cid:62) and x := [(x1)(cid:62), . . . , (xm)(cid:62)](cid:62).

1 , . . . , sα

1 , . . . , xα

n](cid:62), xα := [xα

n](cid:62), sα := [sα

ij. The entry qα

ii = να
i

For the epidemic process at node i, let βi > 0 and δi ≥ 0 be the infection and recovery
rate, respectively. We let Bα > 0 and Dα ≥ 0 be the positive and non-negative diagonal
matrices with entries βi and δi, i ∈ {1, . . . , n}, respectively. Let B and D be the positive
and non-negative diagonal matrices with block-diagonal entries Bα and Dα, α ∈ {1, . . . , m},
respectively. Let P α := diag(pα), P := diag(p) and S := diag(s). We now derive the

continuous time dynamics that captures the interaction of mobility and the SIR epidemic

dynamics.

32

Proposition 3 (SIR model under mobility) The

dynamic model for SIR epidemic process with multi-layer Markovian mobility is

˙s = −SBF (x)p − L(x)s
˙p = (SBF (x) − D − L(x))p
˙xα = (Qα)(cid:62)xα,

(3.9a)

(3.9b)

(3.9c)

where L is an nm×nm block-diagonal matrix with block-diagonal terms Lα, α ∈ {1, . . . , m},

Lα(x) is a matrix with entries

lα
ij(x) =



(cid:80)

j(cid:54)=i qα
ji
xα
j
xα
i

ji

,

−qα

xα
j
xα
i

,

if i = j,

otherwise,

F (x) := [ ¯F(cid:62)(x), . . . , ¯F(cid:62)(x)](cid:62) be a row-concatenated nm×nm matrix with each n×nm block-
row as ¯F (x) := [F 1(x), . . . , F m(x)], and F α as a diagonal matrix with entries f α
i(cid:80)

, i.e., the fraction of total population at node i contributed by class α.

i (x) :=

xα
α xα
i

Proof: The proof follows similarly to that in Proposition 2.

(cid:3)

We analyze the SIR model under multi-layer mobility (3.9) under the strongly connected

assumption: Assumption 2 as well as following standard assumption:

Assumption 3 There exists a node k such that δk > 0.

(cid:3)

Let vα be the right eigenvector of (Qα)(cid:62) associated with eigenvalue at 0. We assume
that vα is scaled such that its inner product with the associated left eigenvector 1n is unity,
i.e., 1(cid:62)
n vα = 1. Define v := [N 1(v1)(cid:62), . . . , N m(vm)(cid:62)](cid:62), where N α is the total number of
individuals belonging to class α, for α ∈ {1, . . . , m}.

Theorem 3 (Existence and properties of equilibria) For the SIR model with multi-

layer Markovian mobility (3.9) under Assumptions 2 and 3, the following statements hold

33

i. if p(0) and s(0) ∈ [0, 1]nm , then p(t) and s(t) ∈ [0, 1]nm for all t > 0 ;

ii. if p(0) > 0 and sα(0) > 0 for each α, then p(t) (cid:29) 0 and s(t) (cid:29) 0 for all t > 0 ;

iii. the equilibrium points (p∗, s∗, x∗) belong to the set

{(0, [k11(cid:62)

n , k21(cid:62)

n , ..., km1(cid:62)

n ](cid:62), v) | k1, k2, ..., km ∈ R≥0};

iv. the set of equilibria

n , k21(cid:62)

n , ..., km1(cid:62)

n ](cid:62), v) | k1, k2, ..., km ∈ R≥0} is globally asymptotically at-

{(0, [k11(cid:62)
tractive.

(i) and (ii) follow similarly to the proof in Theorem 2. Define L∗ := L(x∗),
Proof:
S∗ := diag(s∗) and F∗ := F (x∗) . To establish statement (iii), premultiply (3.9a) and
(3.9b) with x∗(cid:62) at equilibrium

−x∗(cid:62)S∗BF∗p∗ = 0
x∗(cid:62)S∗BF∗p∗ − x∗(cid:62)Dp∗ = 0.

(3.10a)

(3.10b)

Here, we have used the fact that x∗(cid:62)L∗ = 0, which can be seen from the fact that x(cid:62)L(x) =
x(cid:62)Q. Also, since x∗ (cid:29) 0 and S∗BF∗p∗ ≥ 0, Dp∗ ≥ 0, (3.10) yields

S∗BF∗p∗ = 0
Dp∗ = 0.

(3.11a)

(3.11b)

Using Assumption 3 in (3.11b) implies p∗α
k = 0 for each α at node k with δk > 0. Using (3.11)
in (3.9b) at equilibrium gives L∗p∗ = 0 or equivalently L∗αp∗α = 0 for each α. Therefore
under strong connectivity assumption of each layer (Assumption 2) p∗ = 0. Further using
(3.11a) in (3.9a) at equilibrium yields L∗s∗ = 0, or equivalently L∗αs∗α = 0 which gives:
s∗α = kα1n for each α. This proves statement (iii).
For statement (iv), consider a Lyapunov candidate function V3 = x(cid:62)(2s + p). It follows

34

that

˙V3 = x(cid:62)(−2SBF p − 2Ls + SBF p − Dp − Lp)

+ x(cid:62)Q(2s + p)
= x(cid:62)(−SBF p − Dp)
≤ 0.

Now, using LaSalle’s invariance theorem [31, Theorem 4.4], all trajectory asymptotically goes

to the largest invariant set with ˙V3 = 0. This further implies all trajectory asymptotically
goes to an invariant set with Dp = 0 and SBF p = 0. Using this fact in (3.9b) at equilibrium
expanded for each mobility layer under Assumptions 2 and 3 implies p∗ = 0 is globally
attractive.

Next consider a Lyapunov candidate function V4 = s(cid:62)X∗s−2nmr(cid:82) t

diag(x∗) and r := (cid:107)X∗(cid:107). Then,

0 ((cid:107) ˜L(cid:107))dt with X∗ :=

˙V4 = 2s(cid:62)X∗ ˙s − 2nmr((cid:107) ˜L(cid:107))

= −2s(cid:62)X∗BF Sp + s(cid:62)(X∗(−L) + (−L)(cid:62)X∗)s
− 2nmr((cid:107) ˜L(cid:107))
= −2s(cid:62)X∗BF Sp − s(cid:62)(X∗(L∗) + (L∗)(cid:62)X∗)s
− s(cid:62)(X∗( ˜L) + ( ˜L)(cid:62)X∗)s − 2nmr((cid:107) ˜L(cid:107))
≤ −2s(cid:62)X∗BF Sp − s(cid:62)(X∗(L∗) + (L∗)(cid:62)X∗)s
+ 2nmr((cid:107) ˜L(cid:107)) − 2nmr((cid:107) ˜L(cid:107))
≤ −2s(cid:62)X∗BF Sp − s(cid:62)(X∗(L∗) + (L∗)(cid:62)X∗)s
≤ −s(cid:62)(X∗(L∗) + (L∗)(cid:62)X∗)s ≤ 0.

The last inequality follows as the matrix X∗(L∗) + (L∗)(cid:62)X∗ is a symmetric Laplacian
and hence a symmetric positive semi-definite matrix. To see this, note that (X∗(L∗) +
(L∗)(cid:62)X∗)1 = X∗(L∗)1+(L∗)(cid:62)X∗1 = 0+(L∗)(cid:62)x∗ = 0. Additionally, this matrix is a block

35

diagonal matrix with block elements as strongly connected symmetric Laplacian matrices.
n ](cid:62).
Using Barbalat’s lemma, we get ˙V4 → 0. This in turn leads to s → [k11(cid:62)
(cid:3)
This proves statement (iv).

n , . . . , km1(cid:62)

n , k21(cid:62)

An epidemic outbreak is an event in which the total number of infected individuals in

the system (summed over all the layers and nodes) increase before eventually reaching a

disease-free state. As evident from Theorem 3, the total number of infected individuals

ultimately goes to zero. The epidemic outbreak is characterized by the increase in the size

of the infected population in the early phase of the transient response.

Define smax(t) as the greatest element in s(t) taken over all layers and nodes.

Corollary 3 (Epidemic outbreak) For the SIR epidemic model under multi-layer Marko-

vian mobility (3.9) under Assumption 2, the following statements hold

i. For a single layer network, if smax(0)B − D ≤ 0 then there is no epidemic outbreak

and total infected population monotonically decreases to zero;

ii. If S(0)BF (0) − D > 0, then there is an epidemic outbreak at t = 0.

Proof: Using (3.9b), we first write the expression for the rate of change of total infected

population for the system, NI = x(cid:62)p :

˙NI = x(cid:62) ˙p + ˙x(cid:62)p

= x(cid:62)(SBF (x) − D − L(x))p + x(cid:62)Qp
= x(cid:62)(SBF (x) − D)p

(3.12)

where (3.12) follows using x(cid:62)L(x) = x(cid:62)Q, a consequence of definitions of matrices L(x)
and Q. It can be shown from (3.9a) that smax(t) monotonically decreases with time. This is
a consequence of negative first term and negative Laplacian second term in the right hand

side of (3.9a). Now, for the special case of a single layer network F (x) = I, therefore in the
right hand side of (3.12), we can see that SB−D ≤ smax(t)B−D ≤ smax(0)B−D. Further,

36

since x and p are non-negative, if smax(0)B − D ≤ 0, then for a single layer network the
right hand side of (3.12) is non-positive and hence NI monotonically decreases. This proves
statement (i).

Statement (ii) follows by evaluating SBF − D at t = 0 and making it positive to make
right side of (3.12) positive at t = 0 and hence NI increases at t = 0 giving rise to an initial
(cid:3)
outbreak.

3.5 Numerical Illustrations

In this section, we numerically illustrate our results on multi-layer SIR epidemic model.

We choose different population size for the two mobility layers and select the mobil-

ity transition rates using the Metropolis-Hastings algorithm [30] such that the equilibrium

distribution of population is the same for both the layers (taken as uniform equilibrium

distribution).

Figures 3.3 (a) and (b) show the fractions of infected population whereas Figures 3.3 (c)

and (d) show the fraction of susceptible population trajectories for 10 nodes connected with

2-layers of graph structures. Layer 1 is line graph and layer 2 is ring graph. The initial

population distribution at the 10 nodes for layer 1 and layer 2 are

5 × [700, 500, 300, 100, 500, 700, 800, 900, 600, 500],
5 × [300, 300, 200, 100, 200, 300, 400, 400, 500, 200],

and

respectively. The infection and curing rates for the 10 nodes are

[0.31, 0.32, 0.35, 0.36, 0.5, 0.3, 0.3, 0.1, 0.1, 0.1],

and

[0.3, 0.22, 0.21, 0.25, 0.3, 0.21, 0.23, 0.24, 0.21, 0.22],

respectively. The initial fraction of infected population is taken as 0.01 for all layers and

nodes with no recovered population.

37

(a) Infected population;
graph 1

(b) Infected population;
graph 2

(c) Susceptible popula-
tion; graph 1

(d) Susceptible popula-
tion; graph 2

Figure 3.3: Simulation of deterministic model of SIR epidemic spread under 2 layer

mobility, over Line-Ring graph structure. n = 10, pi(0) = 0.01.

38

CHAPTER 4

SIS EPIDEMIC MODEL UNDER MOBILITY ON MULTI-LAYER

NETWORKS: WEAKLY CONNECTED CASE

In this chapter, we relax the assumption that a layer need to be strongly connected. We

characterize the stability of SIS model under mobility for a weakly connected multi-layer

network.

4.1 SIS Model under Mobility with Non-strongly Connected Layers

Figure 4.1: Non-strongly connected multilayer network with shaded sink components

In this section, we relax the Assumption 1, so that the digraph representing mobility on

a layer need not be strongly connected. For each layer, at equilibrium, the entire population

will be restricted to strongly connected sink components and hence we track infected popu-

lation only among nodes belonging to these sink components. Furthermore, at equilibrium,

the effect of nodes shared between different layers appears in the same fashion as before
through mixing matrix F∗ of the sink component nodes under consideration. Note that the
matrix BF∗ − D − L∗ considered for sink component nodes only, need not be strongly con-
nected (as the inter-layer connections among restricted nodes shared between different layers

need not be strongly connected) in contrast to the earlier case. Therefore, we can partition

the total restricted nodes, comprised of strongly connected sink components on each layer,

into different multi-layer network batches, which are strongly connected in themselves but at

equilibrium disconnected from each other. Figure 4.1 shows the depiction. Note that there

39

Layer 1Layer 2Batch 1Batch 1Batch 3Batch 2may be nodes which do not belong to any of these batches. For any given batch, the effect

of nodes outside that batch will be shown to be vanishing at equilibrium with ISS cascaded

connection to that batch. We develop the formulation as below:

Consider a multi-layered population of n nodes comprising of m class. The dynamics of

fraction of infected population for each class at each node is given by (2) reproduced below:

˙p = ((I − P )BF (x) − D − L(x))p
˙xα = (Qα)(cid:62)xα,

(4.1a)

(4.1b)

Since the strongly connected sink components for all layers can be grouped into different
batches, we consider a general batch comprising of 1 ≤ l ≤ m layers with class α contributing
qα strongly connected sink components. Stacking all the sink components of each layer
one below the other and then such stacks one below the other gives s := [s1, ..., sl] with

sk := [s(k,1), ..., s(k,qk)]. The dynamics is given by writing the dynamics of the cocnerned

class of concerned nodes out of (4.1) and separating the contribution from the rest of the

elements:

˙s = ((I − S) ¯BF1(x) − ¯D − L1(x))s

+ (I − S) ¯BF2(x)u − L2(x)u
= ((I − S) ¯BF1(x) − ¯D − L1(x))s
+ (I − S) ¯Bc1 + c2

(4.2)

where, the elements of the matrix F and the matrix L for the concerned classes on

concerned nodes are written as F1 + F2 and L1 + L2 with F1 and L1 being square matrices
formed by retaining the rows i and columns i of the matrix F and L, F2 and L2 being
rectangular formed by rows i and columns j where, pi is a concerned element of the original
vector p, and pj is not. It is clear from the definition of the matrix F2 whose elements are 0

40

j(cid:80)

xα
α xα
j

xα
j
xα
i

ji

or

and elements of matrix −L2 with elements qα

, that the elements (≥ 0) of these
matrices are either identically zero or has globally asymptotically stable equilibrium at zero.

1 − ¯D − L∗

Hence, c1 and c2 asymptotically go to zero. Note that at equilibrium the block diagonal
matrix L∗
1 comprises of strongly connected laplacian matrices as the diagonal blocks. Also,
the matrix ¯BF∗
1 now determines the existence and nature of equilibrium point as
before. Unlike the strongly connected case, the block diagonal ¯Dα for a strongly connected
sink component on a layer may be zero resulting in s∗α = 1. To avoid this situation here
we assume each δi > 0. We show the stability analysis as below:

4.2 Stability Analysis

First, we use an ISS based approach to prove the stability of the SIS model over strongly
connected multi-layer networks: Since, (Qα)(cid:62) has exactly one eigenvalue with multiplicity
1 at zero and all others in negative half plane, we can have a coordinate transformation
zα = (V α)−1x such that zα = [N α, y(cid:62)](cid:62) with yα = 0 (an n − 1 dimensional zero vector)
and N α is the total number of people in class α. It is easy to see that yα = 0 is globally
asymptotically stable. Defining z = [(z1)(cid:62), ..., (zm)(cid:62)](cid:62) and y = [(y1)(cid:62), ..., (ym)(cid:62)](cid:62), y = 0
(an m(n-1) dimensional zero vector) is globally asymptotically stable. Also, with a given

population of each class, we can write F (y) and L(y) as functions of y alone and it can
further be seen that (cid:107) ˜F (y)(cid:107) and (cid:107) ˜L(y)(cid:107) are positive definite functions of y. Consider the
cascaded system:

˙p = (BF∗ − D − L∗)p − P BF (y)p + (B ˜F (y) − ˜L(y))p

˙y = Ay.

(4.3a)

(4.3b)

41

Consider the Lyapunov function V (p) = p(cid:62)Rp. Therefore,

˙V = 2p(cid:62)R ˙p

= p(cid:62)(R(BF∗ − D − L∗) + (BF∗ − D − L∗)(cid:62)R)p
+ 2p(cid:62)R(B ˜F − ˜L)p − 2p(cid:62)RP B( ˜F + F∗)p
= −p(cid:62)Kp + 2p(cid:62)R(B(I − P ) ˜F − ˜L)p − 2p(cid:62)RP BF∗p
≤ −2p(cid:62)RP BF∗p + 2k1(b(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107)), k1, b > 0
≤ −2(1 − 
)p(cid:62)RP BF∗p − 2
p(cid:62)RP BF∗p
+ 2k1(b(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107)), ∀0 < 
 < 1
≤ −2(1 − 
)p(cid:62)RP BF∗p − k2(cid:107)p(cid:107)3∞ + 2k1(b(cid:107) ˜F(cid:107) + (cid:107) ˜L(cid:107))
≤ −2(1 − 
)p(cid:62)RP BF∗p − k2(cid:107)p(cid:107)3∞ + ρ1((cid:107)y(cid:107)∞)
≤ −2(1 − 
)p(cid:62)RP BF∗p, (cid:107)p(cid:107)∞ ≥ ρ2((cid:107)y(cid:107)∞)

(4.4)

Here, ρ1 is a suitable class κ function as per [31, Lemma 4.3] and define ρ2 = 3(cid:113) ρ1

,
k2
a class κ function of (cid:107)y(cid:107)∞. Since the last term in (4.4) is a continous negative definite
function of p in pi ∈ [0, 1], it follows from [31, Theorem 4.19] that dynamics (4.3a) is ISS
(Input to state stable) with respect to input y. Note that (4.3b) has globally asymptotically

stable equilibrium at origin. Further using [31, Lemma 4.7] over the cascaded system with

y as input to dynamics for p implies that the cascaded system has globally asymptotically

stable equilibrium at origin. Next, we modify this prove for the case of weakly connected

networks as shown below:
Stability of epidemic equilibria: Consider the Lyapunov function V (s) = s(cid:62)Rs. There-

42

fore,

˙V = 2s(cid:62)R ˙s

1) + ( ¯BF∗

1 − ¯D − L∗

1)(cid:62)R)s
1 − ¯D − L∗
= s(cid:62)(R( ¯BF∗
+ 2s(cid:62)R( ¯B ˜F1 − ˜L1)s − 2s(cid:62)RS ¯B( ˜F1 + F∗
1 )s
+ 2s(cid:62)(I − S) ¯Bc1 + 2s(cid:62)c2
= −s(cid:62)Ks + 2s(cid:62)R(B(I − S) ˜F1 − ˜L1)s − 2s(cid:62)RS ¯BF∗
1 s
+ 2s(cid:62)(I − S) ¯Bc1 + 2s(cid:62)c2
≤ −2s(cid:62)RS ¯BF∗

1 s + 2k1(b1(cid:107) ˜F1(cid:107) + b2(cid:107) ˜L1(cid:107) + b3(cid:107)c1(cid:107) + (cid:107)c2(cid:107))

1 s − 2
s(cid:62)RS ¯BF∗
1 s

, k1, b1, b2, b3 > 0
≤ −2(1 − 
)s(cid:62)RS ¯BF∗
+ 2k1(b1(cid:107) ˜F1(cid:107) + b2(cid:107) ˜L1(cid:107) + b3(cid:107)c1(cid:107) + (cid:107)c2(cid:107)), ∀0 < 
 < 1
≤ −2(1 − 
)s(cid:62)RS ¯BF∗
+ 2k1(b1(cid:107) ˜F1(cid:107) + b2(cid:107) ˜L1(cid:107) + b3(cid:107)c1(cid:107) + (cid:107)c2(cid:107))
≤ −2(1 − 
)s(cid:62)RS ¯BF∗
≤ −2(1 − 
)s(cid:62)RS ¯BF∗

1 s − k2(cid:107)s(cid:107)3∞ + ρ1((cid:107)y(cid:107)∞)
1 s, (cid:107)s(cid:107)∞ ≥ ρ2((cid:107)y(cid:107)∞)

1 s − k2(cid:107)s(cid:107)3∞

(4.5)

Here, ρ1 is a suitable class κ function as per [31, Lemma 4.3] and define ρ2 = 3(cid:113) ρ1

, a
class κ function of (cid:107)y(cid:107)∞. Note, that a vector y can still be constructed as in (4.3b) for
a weakly connected graph, where c1(y), c2(y), ˜F1(y) and ˜L1(y) are functions of y and y
asymptotically goes to zero. Since the last term in (4.5) is a continous negative definite
function of s in pi ∈ [0, 1], it follows from [31, Theorem 4.19] that dynamics (4.2) is ISS
(Input to state stable) with respect to input y. The rest of the argument follows same as

k2

the ISS based epidemic equilibrium stability proof presented above.

Stability of endemic equilibria: The stability of endemic equilibrium can be shown by

doing Lyapunov analysis similar to the strongly connected case and the extra terms can be

incorporated into ISS terms similar to the previous proof.

43

CHAPTER 5

RESOURCE ALLOCATION

In this chapter, we formulate a resource allocation problem for the SIS epidemic model under

mobility on single-layer and multi-layer network. We show that under certain assumptions

the problem can be formulated as a geometric program. We give numerical illustrations to

support our results.

5.1 Resource allocation

We formulate a resource allocation problem for the SIS model under markovian mobility

where we allocate resources so as to have the disease free equilibrium (DFE) as the stable

equilibrium. We consider two types of resources:

i. Preventive resource: This can be used to change the infection rate such that βi ∈

[βi, ¯βi]. This resource is applied to a node with cost function fi(βi).

ii. Corrective resource: This can be used to change the recovery rate such that δi ∈ [δi, ¯δi].

The corresponding cost function is gi(δi).

5.1.1 Single layer mobility
It is shown in Appendix A.3 that the matrix B − D − L∗ is similar to B − D + Q(cid:62), hence
the condition for stability of DFE is given by:

µ(B − D + Q(cid:62)) ≤ 0

Since B − D + Q(cid:62) is an irreducible Metzler matrix, µ is a real and simple eigenvalue
referred to as λ1 henceforth. We use the approach similar to the one in [32]. We use a lemma
resulting from Perron-Frobenius theorem for irreducible non-negative matrix:

44

Lemma 3 If M is an irreducible non-negative matrix, with λ1 = µ(M ) being its radial-
abscissa, then λ1 = inf{λ ∈ R : M u ≤ λu for u (cid:29) 0}.

By applying above lemma one can minimize λ1(M ) by minimizing λ such that:

≤ 1,

M ui
λui

ui > 0.

(5.1)
In order to maximize the decay rate of DFE, we minimize λ1(B− D + Q(cid:62)) with a budget
constraint over total resource cost. It can be seen that λ1(B− D + Q(cid:62)) = λ1(B + ˆD + ˆQ(cid:62))−
¯∆ − 1 − ¯ν, where ˆD = diag(ˆδi), ˆδi = ¯∆ + 1 − δi, ¯∆ = max{¯δi}n
i (cid:54)= j,
j(cid:54)=i qij. Note that minimising λ1(B −
ˆQii = ˆqii = ¯ν − νi, ¯ν = max{νi}n
D + Q(cid:62)) is same as minimising λ1(B + ˆD + ˆQ(cid:62)), where latter corresponds to an irreducible
non-negative matrix. To pose the optimization problem as a Geometric Program, we restrict

i=1, where νi =(cid:80)

i=1, ˆQij = Qij

for

ourselves to preventive and corrective resources only. Further, we assume the preventive cost
f (βi) a posynomial and corrective cost g(δi) = ˜g( ¯∆ + 1 − δi) = ˜g(ˆδi), with ˜g a posynomial.
Consider a strongly connected digraph G = (V,E), where V = {1, . . . , n} is the node (patch)
set and E ⊂ V × V is the edge set. The generator matrix for mobility is Q, whose (i, j)-th
entry is qij, where qij > 0 for (i, j) ∈ E and 0 otherwise. Utilising lemma 3, (5.1) on matrix
(B + ˆD+ ˆQ(cid:62)), the resource allocation problem [32] can be stated as a geometric programming
problem:

λ

min

(cid:80)
λ,ui,βi,ˆδi
j=1:n(B + ˆD + ˆQ(cid:62))ijuj
(cid:88)
[fi(βi) + ˜gi(ˆδi)] ≤ C,

λui

≤ 1,

s.t.

where i ∈ [1, n].

i

βi ≤ βi ≤ ¯βi,
¯∆ + 1 − ¯δi ≤ ˆδi ≤ ¯∆ + 1 − δi.

45

(5.2)

(5.3)

(5.4)

(5.5)

(5.6)

5.1.2 Multilayer mobility

The condition for stability of DFE is given by:

µ(BF∗ − D − L∗) ≤ 0

Since BF∗ − D − L∗ is an irreducible Metzler matrix, µ is a real and simple eigenvalue
referred to as λ1 henceforth. We consider a variable transformation similar to the single
layer case: λ1(BF∗ − D − L∗) = λ1(BF∗ + ˆD + ˆA) − ¯∆ − 1 − ¯ν, where ˆD = blkdiag( ˆDα),
ij = −Lα∗
ˆDα = diag(ˆδi), ˆδi = ¯∆+1−δi, ¯∆ = max{¯δi}n
i (cid:54)= j,
ij. Note that minimising
ˆAα
ii = rα
λ1(BF∗ − D − L∗) is same as minimising λ1(BF∗ + ˆD + ˆA), where latter corresponds to
an irreducible non-negative matrix. Similar to the single-layer case, restricting ourselves

i=1, ˆA = blkdiag( ˆAα), ˆAα

i = (cid:80)

i , ¯ν = max{να

i=1,α=1, where να

i = ¯ν − να

i }n,m

for

ij

j(cid:54)=i qα

to posynomial costs f (βi) and ˜g(ˆδi), the resource allocation problem can be stated as a
geometric problem:

λ

min

(cid:80)
λ,ui,βi,ˆδi
j=1:nm(BF∗ + ˆD + ˆA)ijuj
(cid:88)
[fi(βi) + ˜gi(ˆδi)] ≤ C,

λui

s.t.

≤ 1,

for i = 1 to nm,

for i = 1 to n,

i

βi ≤ βi ≤ ¯βi,
for i = 1 to n,
¯∆ + 1 − ¯δi ≤ ˆδi ≤ ¯∆ + 1 − δi,

for i = 1 to n.

(5.7)

(5.8)

(5.9)

(5.10)

(5.11)

5.2 Numerical Illustration

We give numerical illustration for the budget constrained resource allocation problem for

the case of multi-layer network. We take two layers: layer 1 as line graph and layer 2 as a

ring graph. We use mobility transition rates as those obtained with the condition of equal

outgoing rates to neighboring nodes with total outgoing transition rate fixed as 0.2. The

46

total number of individuals on layer 1 is taken as 300 whereas those on layer 2 as 500. The

cost functions are taken as

g(δi) =

1

¯∆ + 1 − δi

−

f (βi) =

1
βi
¯∆ + 1 − δi

1

− 1
¯βi

,

= ˜g(ˆδi) =

−

1
ˆδi

1

¯∆ + 1 − δi

.

The bounds on the infection and curing rates are taken as

βi = 0.1,

¯βi = 0.4,

δi = 0.1,

¯δi = 0.4.

We solve the resulting geometric programming for minimisation of λ1 corresponding to decay
rate of desease-free equilibrium. Note that λ = λ1(BF∗ + ˆD + ˆA) − 1 − ¯∆ − ¯ν, where λ
corresponds to the decay rate of disease-free equilibrium. Figure 5.1 shows the plot of

obtained minimum λ against the maximum allowable cost (budget) C. As can be seen, λ

is positive for very low values of budget implying unstable disease-free equilibrium. As the
budget C increases λ saturates around a value of −0.3. Note that the actual cost used is
less than or equal to the budget and need not be equal to the budget C.

47

Figure 5.1: Decay rate λ against maximum allowable cost C. Graph1: Line, Graph 2:

Ring, n = 10, νi = 0.2.

48

CHAPTER 6

CONCLUSIONS

We derived a continuous-time model for epidemic propagation under Markovian mobility

across a network of sub-populations. The epidemic spread within each node has been modeled

as SIS and SIR population models. The derived models have been analysed to establish the

existence, uniqueness and stability of disease-free equilibrium and an endemic equilibrium

under different conditions. Some necessary and some sufficient conditions for stability of

disease-free equilibrium have been established. We also provided numerical studies to support

our results and elucidated the effect of mobility on epidemic propagation. We extended the

stability results for the SIS model to the case of non-strongly connected layers. Further, we

formulated a budget constrained resource allocation problem as a geometric program and

provided numerical illustrations.

49

APPENDIX

50

A.1 Proof of Theorem 1 (iii), Chapter 2: Existence of an endemic

equilibrium

We show below that in the case of µ(B−D−L∗) > 0 , there exists an endemic equilibrium

p∗, i.e.,

˙p|p=p∗ = (B − D − L∗ − P∗B)p∗ = 0.

(1)

We use Brouwer’s fixed point theorem, similar to the derivation in [9]. Rearranging the
terms and writing the above as an equation in p to be satisfied at non-trivial equilibrium p∗
leads to:

(L∗ + D)((L∗ + D)−1B − I)p = P Bp.

(2)
Define A := (L∗ +D)−1B. Since A−1 = B−1(L∗ +D) is a non-singular M-matrix, its inverse
A is non-negative [33]. Rearranging (2) leads to

p = H(p) = (I + AP )−1Ap.

(3)

Now we show that H(p) as defined above is a monotonic function in the sense that

p2 ≥ p1 implies H(p2) ≥ H(p1). Define ˜p := p2 − p1 and ˜P := diag(˜p). Then,

H(p2) − H(p1)
= (A−1 + P2)−1p2 − (A−1 + P1)−1p1
= (A−1 + P2)−1(p2 − (A−1 + P2)(A−1 + P1)−1p1)
= (A−1 + P2)−1(˜p − ˜P (A−1 + P1)−1p1)
= (A−1 + P2)−1(I − diag((A−1 + P1)−1p1))˜p.

(4)

Since (A−1 + P2) = B−1(L∗ + D) + P2 is an M-matrix its inverse and hence the first term

51

above is non-negative. The second term is shown to be non-negative as below:

(I − diag((A−1 + P1)−1p1))
= (I − diag((I + AP1)−1AP11n))
= diag((I − (I + AP1)−1AP1)1n)
= diag((I + AP1)−11n)
= diag((A−1 + P1)−1A−11n)
≥ 0,

where we have used the identity:

(I + X)−1 = I − (I + X)−1X,

(5)

(6)

in the second line. The last inequality in (5) holds as A−11n = B−1(L∗+D)1n = B−1D1n ≥
0n and (A−1 + P1)−1 ≥ 0 the inverse of an M-matrix. The last term in the last line of (4)
is ˜p ≥ 0n. This implies that H(p) is a monotonic function. Also, result in (5) implies that
H(p) ≤ 1n for all p ∈ [0, 1]n. Therefore H(1n) ≤ 1n.

Convergent splitting property of irreducible M-matrices [33] implies µ(B − D − L∗) > 0
if and only if R0 = ρ(A) = ρ((L∗ + D)−1B) > 1. Here ρ(A) is spectral radius of A. Since
A is an irreducible non-negative matrix, Perron-Frobenius theorem implies ρ(A) is a simple
eigenvalue with right eigenvector u satisfying Au = ρ(A)u = R0u , with u (cid:29) 0n. Define
. Now, we find a value of 
 > 0 such that H(
u) ≥ 
u as below:
U := diag(u) and γ :=

R0−1
R0

H(
u) − 
u = (I + 
AU )−1A
u − 
u

= (I − (I + 
AU )−1
AU )
R0u − 
u

= 
R0(
= 
R0(γu − (I + 
AU )−1
AU u).

R0

(R0 − 1)

u − (I + 
AU )−1
AU u)

(7)

52

Now, the expression in the brackets in the last line is a continous function of 
 and is
equal to γu (cid:29) 0n at 
 = 0. Therefore, there exists an 
 > 0 such that H(
u) − 
u ≥ 0n
or equivalently, H(
u) ≥ 
u. Taking the closed compact set K = [
u, 1n], H : K → K is a
continuous function. Therefore, by Brouwer’s fixed point theorem, there exists a fixed point
in K. This proves the existence of a non-trivial equilibrium p∗ (cid:29) 0n when µ(B−D−L∗) > 0
or equivalently R0 > 1. The uniqueness is further shown in the following proposition.

Proposition 4 If the mapping H has a strictly positive fixed point, then it is unique.

Proof: The proof is similar to the proof of [10, Proposition A.3] and is given below:

Assume there are two strictly positive fixed points: p∗ and q∗. Define

η := max

p∗
q∗

i

i

,

k := arg max

p∗
q∗

i

i

Therefore, p∗ ≤ ηq∗. Lets assume η > 1. First we will show that H(ηq∗) < ηH(q∗) as
follows:

H(ηq∗) − ηH(q∗)
= (I + AηQ∗)−1Aηq∗ − η(I + AQ∗)−1Aq∗
= ((I + AηQ∗)−1A − (I + AQ∗)−1A)ηq∗
= ((A−1 + ηQ∗)−1 − (A−1 + Q∗)−1)ηq∗
= (A−1 + ηQ∗)−1(I − (A−1 + ηQ∗)(A−1 + Q∗)−1)ηq∗
= (A−1 + ηQ∗)−1(−(η − 1)(A−1 + Q∗)−1)ηq∗

< 0

(8)

where the last inequality uses result that inverse of a non-singular M-matrix is non-negative
and non-singular, that η > 1 and, that q∗ (cid:29) 0 by assumption. Consequently

p∗
k = Hk(p∗) ≤ Hk(ηq∗) < ηHk(q∗) = ηq∗
k,

(9)

53

k = p∗

k < p∗

k by definition, if η > 1, we have from above p∗

Since, ηq∗
k, a contradiction. Hence,
η ≤ 1 which implies p∗ ≤ q∗. By switching the roles of p∗ and q∗ and repeating the above
argument we can show q∗ ≤ p∗. Thus p∗ = q∗ and hence there is a unique strictly positive
(cid:3)
fixed point.

A.2 Proof of Theorem 2 (iii), Chapter 3: Existence of an endemic

equilibrium: Multi-layer Case

Here we assume that there exists atleast one node with positive recovery rate, i.e., δi > 0

for atleast one i. The case with no recovery at all nodes is trivial and leads to p∗ = 1.

We first state some properties of M-matrices, which we will use in the proof.

Theorem 4 (Properties of M-matrix, [33]) For a real Z-matrix (i.e., a matrix with all
off-diagonal terms non-positive) A ∈ Rn×n, the following statements are equivalent to A
being a non-singular M-matrix

i. Stability: real part of each eigenvalue of A is positive;

ii. Inverse positivity: A−1 ≥ 0 (for irreducible A, A−1 > 0);

iii. Regular splitting: A has a convergent regular splitting, i.e., A has a representation
A = M −N, where M−1 ≥ 0, N ≥ 0 (called regular splitting), with M−1N convergent,
i.e., ρ(M−1N ) < 1;

iv. Convergent regular splitting: every regular splitting of A is convergent. Further, for

a singular M-matrix (i.e. singular Z-matrix with real part of eigenvalues non-negative)
regular splitting of A gives ρ(M−1N ) = 1;

v. Semi-positivity: there exists x (cid:29) 0 such that Ax (cid:29) 0;

54

vi. Modified semi-positivity: there exists x (cid:29) 0 such that y = Ax > 0 and matrix ˆA

if Aij (cid:54)= 0 or yi (cid:54)= 0,

0 otherwise

defined by

is irreducible.

ˆAij =

1

A consequence of Theorem 4 (vi) is that an irreducible Laplacian matrix perturbed with

a non-negative diagonal matrix with atleast one positive element is a non-singular M-matrix
(take x = 1 (cid:29) 0). This implies that block diagonal submatrices of the matrix L∗ + D are
all non-singular M-matrices (since δi ≥ 0 with strict inequality for atleast one i) and hence
L∗ + D is a non-singular M-matrix. Similar arguments imply B−1(L∗ + D) is a non-singular
M-matrix.

We show below that in the case of µ(BF∗ − D − L∗) > 0 there exists an endemic

equilibrium p∗ (cid:29) 0, i.e.,

˙p|p=p∗ = (BF∗ − D − L∗ − P∗BF∗)p∗ = 0.

(10)

We use Brouwer’s fixed point theorem, similar to the derivation in [9]. We split the
non-negative matrix F∗ as F∗ = I − M, where M is a Laplacian matrix. Rearranging the
terms and writing the above as an equation in p to be satisfied at p∗ leads to

(L∗ + D)((L∗ + D)−1B − I)p = (P B + (I − P )BM )p
= B(P + (I − P )M )p.

(11)

Define A := (L∗ + D)−1B. Since A−1 = B−1(L∗ + D) is a non-singular M-matrix, its

inverse A is non-negative [33]. Rearranging (11) leads to

p = H(p) = (I + A(P + (I − P )M ))−1Ap.

(12)

55

Now we show that H(p) as defined above is a monotonic function in the sense that

p2 ≥ p1 implies H(p2) ≥ H(p1). Define ˜p := p2 − p1 and ˜P := diag(˜p). Then,

=

=

=

(cid:17)−1
(cid:16)
H(p2) − H(p1)
(cid:17)−1
−(cid:16)
A−1 + P2 + (I − P2)M
(cid:17)−1(cid:16)
(cid:16)
A−1 + P1 + (I − P1)M
p1
(cid:16)
(cid:17)(cid:16)
A−1 + P2 + (I − P2)M
p2−
(cid:17)−1(cid:16)
(cid:16)
A−1 + P1 + (I − P1)M
A−1 + P2 + (I − P2)M
A−1 + P2 + (I − P2)M
= (A−1 + P2 + (I − P2)M )−1(cid:16)
− ˜P (I − M )

A−1 + P1 + (I − P1)M

(cid:17)−1

(cid:16)

(cid:17)

p1

− diag

(I − M )(A−1 + P1 + (I − P1)M )−1p1

p2

˜p

I

(cid:17)−1

(cid:17)

p1

(13)

(cid:16)

(cid:17)(cid:17)

˜p

Since (A−1 + P2 + (I − P2)M ) = B−1(L∗ + D) + P2 + (I − P2)M is a non-singular
M-matrix (consider theorem 4 (vi) with x = 1 (cid:29) 0), its inverse and hence the first term

56

above is non-negative. The second term is shown to be non-negative as below

1

(cid:17)

(14)

(cid:16)

I − diag

I − diag

=

(cid:16)

(cid:17)(cid:17)
(cid:16)
(cid:17)(cid:17)
(cid:16)
(I − M )(A−1 + P1 + (I − P1)M )−1p1
(cid:17)
(cid:16)(cid:16)
(I − M )(A−1 + P1 + (I − P1)M )−1P11
(cid:16)(cid:0)I − M − (I − M )(I + AP1 + A(I − P1)M )−1
I − (I − M )(I + AP1 + A(I − P1)M )−1AP1
(AP1 + A(I − P1)M )(cid:1)1
(cid:16)
(I − M )(cid:0)I − (I + AP1 + A(I − P1)M )−1
(AP1 + A(I − P1)M )(cid:1)1
(cid:17)
(cid:16)
(cid:17)
(cid:16)
(I − M ) (I + AP1 + A(I − P1)M )−1 1

F∗(cid:0)A−1 + P1 + (I − P1)M(cid:1)−1A−11

(cid:17)
(cid:17)

= diag

= diag

= diag

= diag

= diag
≥ 0,

where we have used the identity

(I + X)−1 = I − (I + X)−1X,

(15)
and M 1 = 0 , as M is a Laplacian matrix. The last inequality in (14) holds as A−11 =
B−1(L∗ + D)1 = B−1D1 ≥ 0 and (A−1 + P1 + (I − P1)M )−1 ≥ 0, since it is the inverse
of an M-matrix. The last term in the last line of (13) is ˜p ≥ 0. This implies that H(p) is
a monotonic function. Also, argument similar to above can be used to show that H(p) ≤ 1
for all 0 ≤ p ≤ 1. Therefore, H(1) ≤ 1.

Applying the converse of Theorem 4 (iv), with Z-matrix as (L∗ + D) − BF∗, where
(L∗ + D)−1 ≥ 0, BF∗ ≥ 0 implies µ (BF∗ − (D + L∗)) > 0 if and only if R0 = ρ(AF∗) =
ρ(A(I − M )) > 1. Now, A is a block-diagonal matrix with block-diagonal terms as Aα =
(L∗α + Dα)−1Bα, which are inverse of irreducible non-singular M-matrices and hence are
positive. Using the expression for F gives AF∗ = [(A1)(cid:62) ¯F(cid:62)(x∗), . . . , (Am)(cid:62) ¯F(cid:62)(x∗)](cid:62).

57

Since Aα > 0 and ¯F∗ ≥ 0 with no zero column, AF∗ > 0 and hence irreducible. Since AF∗
is an irreducible non-negative matrix, Perron-Frobenius theorem implies ρ(AF∗) is a simple
eigenvalue satisfying AF∗u = ρ(AF∗)u = R0u with u (cid:29) 0. Using F∗ = I − M implies:

Au = R0u + AM u

= (R0 − 1)u + (I + AM )u.
R0−1
R0

 ∈ (0, 
0) implies H(
u) ≥ 
u as below:

Define U := diag(u) and γ :=

(16)

. Putting p = 
u , we show that ∃ 
0 such that

H(
u) − 
u

=(cid:0)I + 
AU + A(I − 
U )M(cid:1)−1A
u − 
u
(cid:16)(cid:0)I + 
AU + A(I − 
U )M(cid:1)−1(R0 − 1)u
+(cid:0)I + 
AU + A(I − 
U )M(cid:1)−1(I + AM )u − u

= 


(cid:17)

≡ 
K(
).

Now we evaluate K(
) at 
 = 0 :

K(0)

= (I + AM )−1(R0 − 1)u + (I + AM )−1(I + AM )u − u
= (I + AM )−1(R0 − 1)u
= (R0 − 1)(A−1 + M )−1A−1u
(A−1 + M )−1F∗u

(R0 − 1)

R0

=
= γ(A−1 + M )−1F∗u
(cid:29) 0.

(17)

(18)

The last inequality follows as γ and u are both positive, and (A−1+M )−1F∗ = (B−1(L+
D) + M )−1F∗ > 0 as B−1(L + D) + M is an irreducible M-matrix and hence its inverse

58

is positive and F∗ ≥ 0 with no zero column. Since K(
) is a continuous function of 
 , ∃

0 such that 
0 > 
 > 0 implies K(
) (cid:29) 0 and therefore, H(
u) ≥ 
u. Therefore there
exists an 
 > 0 such that H(
u) − 
u ≥ 0 or equivalently, H(
u) ≥ 
u. Taking the closed
compact set J = [
u, 1], H(p) : J → J is a continuous function of p. Brouwer’s fixed point
theorem implies there exists a fixed point of H in J. This proves the existence of an endemic
equilibrium p∗ (cid:29) 0 when µ(BF∗ − D − L∗) > 0 or equivalently R0 > 1. The uniqueness is
further shown in the following proposition.

Proposition 5 If the mapping H has a strictly positive fixed point, then it is unique.

Proof: Assume there are two strictly positive fixed points: 0 (cid:28) p∗ (cid:28) 1 and 0 (cid:28) q∗ (cid:28)
1. Strict inequality compared to 1 is assumed which can be easily proved for any equilibrium

point using (2.4) under Assumptions 1 and 3. Define
p∗
q∗

k := arg max

η := max

,

,

p∗
q∗

i

i

zi = min(ηq∗

i , 1)

i

i

Therefore, p∗ ≤ z ≤ ηq∗. Lets assume η > 1, which implies q∗ (cid:28) z. First we will show
that H(z) (cid:28) ηH(q∗) as follows:

H(z) − ηH(q∗)
= (A−1 + Z + (I − Z)M )−1z
− η(A−1 + Q∗ + (I − Q∗)M )−1q∗
≤ (A−1 + Z + (I − Z)M )−1
(I − (A−1 + Z + (I − Z)M )(A−1 + Q∗ + (I − Q∗)M )−1)ηq∗
= W−1(I − I − ((Z − Q∗) − (Z − Q∗)M ))(A−1 + Q∗ + (I − Q∗)M )−1)ηq∗
= −W−1(Z − Q∗)(I − M )(A−1 + Q∗ + (I − Q∗)M )−1ηq∗
= −W−1(Z − Q∗)F∗(A−1 + Q∗ + (I − Q∗)M )−1ηq∗
(cid:28) 0

(19)

59

where the last inequality uses result that, W−1 and (A−1 + Q∗ + (I − Q∗)M )−1 inverse of
non-singular M-matrices are non-negative and non-singular (hence with no zero rows), that
F∗ is non-negative with no zero rows, that Z−Q∗ has all elements positive and, that q∗ (cid:29) 0
by assumption. Consequently

k = Hk(p∗) ≤ Hk(z) < ηHk(q∗) = ηq∗
p∗
k,
k < p∗

k = p∗

k by definition, if η > 1, we have from above p∗

Since, ηq∗
k, a contradiction. Hence,
η ≤ 1 which implies p∗ ≤ q∗. By switching the roles of p∗ and q∗ and repeating the above
argument we can show q∗ ≤ p∗. Thus p∗ = q∗ and hence there is a unique strictly positive
(cid:3)
fixed point.
A.3 Similarity of B − D − L∗ and B − D + Q(cid:62)

(20)

Using definition of L∗ it can be seen that:
x∗
x∗

(cid:88)

L∗
ii =

j

,

i

qji
x∗
x∗

i

i

j(cid:54)=i
= −qii
= −qii

(Q(cid:62)x∗ = 0)

,

Also, Lij = −qji

xj
xi

. This implies:

L∗ = −(X∗)−1Q(cid:62)X∗,

B − D − L∗ = (X∗)−1(B − D + Q(cid:62))X∗

Here we have used the fact that B and D are diagonal matrices.

60

BIBLIOGRAPHY

61

BIBLIOGRAPHY

[1] Roy M Anderson and Robert M May. Infectious Diseases of Humans: Dynamics and

Control. Oxford University Press, 1992.

[2] D. Easley and J. Kleinberg. Networks, Crowds, and Markets: Reasoning About a Highly

Connected World. Cambridge University Press, 2010.

[3] Jon Kleinberg. The wireless epidemic. Nature, 449(7160):287, 2007.

[4] Pu Wang, Marta C González, César A Hidalgo, and Albert-László Barabási. Under-
standing the spreading patterns of mobile phone viruses. Science, 324(5930):1071–1076,
2009.

[5] Xiaolan Zhang, Giovanni Neglia, Jim Kurose, and Don Towsley. Performance modeling

of epidemic routing. Computer Networks, 51(10):2867–2891, 2007.

[6] Fang Jin, Edward Dougherty, Parang Saraf, Yang Cao, and Naren Ramakrishnan. Epi-
demiological modeling of news and rumors on twitter. In Proceedings of the Workshop
on Social Network Mining and Analysis, page 8. ACM, 2013.

[7] C. Nowzari, V. M. Preciado, and G. J. Pappas. Analysis and control of epidemics: A
survey of spreading processes on complex networks. IEEE Control Systems Magazine,
36(1):26–46, Feb 2016.

[8] Wenjun Mei, Shadi Mohagheghi, Sandro Zampieri, and Francesco Bullo. On the dynam-
ics of deterministic epidemic propagation over networks. Annual Reviews in Control,
44:116–128, 2017.

[9] A Fall, Abderrahman Iggidr, Gauthier Sallet, and Jean-Jules Tewa. Epidemiologi-
cal models and Lyapunov functions. Mathematical Modelling of Natural Phenomena,
2(1):62–83, 2007.

[10] Ali Khanafer, Tamer Başar, and Bahman Gharesifard. Stability of epidemic models

over directed graphs: A positive systems approach. Automatica, 74:126–134, 2016.

[11] Hyoung Jun Ahn and Babak Hassibi. Global dynamics of epidemic spread over complex
networks. In 52nd IEEE Conference on Decision and Control, pages 4579–4585. IEEE,
2013.

[12] Navid Azizan Ruhi and Babak Hassibi. SIRS epidemics on complex networks: Concur-
rence of exact Markov chain and approximated models. In 2015 54th IEEE Conference
on Decision and Control (CDC), pages 2919–2926. IEEE, 2015.

[13] Mehran Mesbahi and Magnus Egerstedt. Graph Theoretic Methods in Multiagent Net-

works, volume 33. Princeton University Press, 2010.

62

[14] Jesús Gómez-Gardeñes, Alessandro S de Barros, Suani TR Pinho, and Roberto FS
Andrade. Abrupt transitions from reinfections in social contagions. Europhysics Letters,
110(5):58006, 2015.

[15] Renato Pagliara, Biswadip Dey, and Naomi Ehrich Leonard. Bistability and resurgent

epidemics in reinfection models. IEEE Control Systems Letters, 2(2):290–295, 2018.

[16] VS Bokharaie, O Mason, and F Wirth. Spread of epidemics in time-dependent net-
works. In Proceedings of the 19th International Symposium on Mathematical Theory of
Networks and Systems–MTNS, volume 5, 2010.

[17] M. Ogura and V. M. Preciado. Stability of spreading processes over time-varying large-
scale networks. IEEE Transactions on Network Science and Engineering, 3(1):44–57,
Jan 2016.

[18] P. E. Paré, C. L. Beck, and A. Nedić. Epidemic processes over time-varying networks.

IEEE Transactions on Control of Network Systems, 5(3):1322–1334, Sep. 2018.

[19] Vittoria Colizza and Alessandro Vespignani. Epidemic modeling in metapopulation
systems with heterogeneous coupling pattern: Theory and simulations. Journal of
Theoretical Biology, 251(3):450–467, 2008.

[20] Joan Saldaña. Continuous-time formulation of reaction-diffusion processes on heteroge-

neous metapopulations. Phys. Rev. E, 78:012902, Jul 2008.

[21] Wendi Wang and Xiao-Qiang Zhao. An epidemic model in a patchy environment.

Mathematical biosciences, 190(1):97–112, 2004.

[22] Yu Jin and Wendi Wang. The effect of population dispersal on the spread of a disease.

Journal of mathematical analysis and applications, 308(1):343–364, 2005.

[23] Michael Y Li and Zhisheng Shuai. Global stability of an epidemic model in a patchy

environment. Canadian Applied Mathematics Quarterly, 17(1):175–187, 2009.

[24] D Soriano-Paños, L Lotero, A Arenas, and J Gómez-Gardeñes. Spreading processes in
multiplex metapopulations containing different mobility networks. Physical Review X,
8(3):031039, 2018.

[25] Julien Arino, Jonathan R Davis, David Hartley, Richard Jordan, Joy M Miller, and
P Van Den Driessche. A multi-species epidemic model with spatial dynamics. Mathe-
matical Medicine and Biology, 22(2):129–142, 2005.

[26] Lewi Stone, Frank Hilker, and Guy Katriel. SIR models. In Encyclopedia of Theoretical

Ecology, pages 648–658. University of California Press, 2012.

[27] F. Bullo. Lectures on Network Systems. CreateSpace, 1 edition, 2018. With contributions

by J. Cortes, F. Dorfler, and S. Martinez.

[28] Jean-Jacques E Slotine and Weiping Li. Applied Nonlinear Control, volume 199. Prentice

Hall Englewood Cliffs, NJ, 1991.

63

[29] Chai Wah Wu. On bounds of extremal eigenvalues of irreducible and m-reducible ma-

trices. Linear Algebra and its Applications, 402:29–45, 2005.

[30] W. K. Hastings. Monte carlo sampling methods using markov chains and their appli-

cations. Biometrika, 57(1):97–109, 1970.

[31] H. K. Khalil. Nonlinear Systems. Prentice Hall, third edition, 2002.

[32] Cameron Nowzari, Victor M Preciado, and George J Pappas. Optimal resource al-
location for control of networked epidemic models. IEEE Transactions on Control of
Network Systems, 4(2):159–169, 2015.

[33] Abraham Berman and Robert J Plemmons. Nonnegative Matrices in the Mathematical

Sciences, volume 9. Siam, 1994.

64