Stability of Spreading Processes over Time-Varying Large-Scale Networks

Stability of Spreading Processes over Time-Varying Large-Scale Networks

Masaki Ogura,  and Victor M. Preciado,  The authors are with the Department of Electrical and Systems Engineering at the University of Pennsylvania, Philadelphia, PA 19104.
Abstract

In this paper, we analyze the dynamics of spreading processes taking place over time-varying networks. A common approach to model time-varying networks is via Markovian random graph processes. This modeling approach presents the following limitation: Markovian random graphs can only replicate switching patterns with exponential inter-switching times, while in real applications these times are usually far from exponential. In this paper, we introduce a flexible and tractable extended family of processes able to replicate, with arbitrary accuracy, any distribution of inter-switching times. We then study the stability of spreading processes in this extended family. We first show that a direct analysis based on Itô’s formula provides stability conditions in terms of the eigenvalues of a matrix whose size grows exponentially with the number of edges. To overcome this limitation, we derive alternative stability conditions involving the eigenvalues of a matrix whose size grows linearly with the number of nodes. Based on our results, we also show that heuristics based on aggregated static networks approximate the epidemic threshold more accurately as the number of nodes grows, or the temporal volatility of the random graph process is reduced. Finally, we illustrate our findings via numerical simulations.

Dynamic random graphs, complex networks, epidemics, stochastic processes, random matrix theory.

1 Introduction

Understanding spreading processes in complex networks is a central question in the field of Network Science with applications in rumor propagation [1], malware spreading [2], epidemiology, and public health [3]. In this direction, important advances in the analysis of spreading processes over time-invariant (TI) networks have been made during the last decade (see [4] and [5] for recent surveys). Under the time-invariance assumption, we find several recent breakthroughs in the literature, such as the rigorous analysis of mean-field approximations [6], the connection between the eigenvalues of a network and epidemic thresholds [7, 8, 9], new modeling tools for analysis of multilayer networks [10], and the use of control and optimization tools to contain epidemic outbreaks [11, 12, 13, 14].

In practice, most spreading processes of practical interest take place in networks with time-varying (TV) topologies, such as human contact networks, online social networks, biological, and ecological networks [15]. A naive, but common, approach to analyze dynamic processes in TV networks is to build an aggregated static topology based on time averages. When the time scales of the network evolution are comparable to those of the spreading process, the static aggregated graph is not well-suited to study the dynamics of TV networks, as pointed out in [16, 17]. Using extensive numerical simulations, it has been observed that the speed of spreading of a disease in a TV network can be substantially slower than in its aggregated static representation [18]. This observation is also supported by the study of spreading processes in a real TV network constructed from a mobile phone dataset [19]. More recent studies point out the key role played by the addition and removal of links [20], as well as the distribution of contact durations between nodes [21] in the dynamics of the spread. The works mentioned above provide empirical evidence about the nontrivial effect that the dynamics of the network has on the behavior of spreading processes.

Apart from these empirical studies, we also find several works providing theoretical support to these numerical observations. The authors in [22, 23] derived the value of the epidemic threshold in particular types of TV networks, assuming that all nodes in the network present a homogeneous infection and recovery rates. In [24], Perra et al. proposed a model of temporal network, called the activity-driven model, able to replicate burstiness in node/edge activity and analyzed contagion processes taking place in this model in [25]. Karsai et al. proposed in [26] a time-varying network using a reinforcement process able to replicate the emergence of strong and weak ties observed in a mobile call dataset. A wide and flexible class of TV network model, called edge-Markovian graphs, was proposed in [27] and analyzed in [28]. In this model, edges appear and disappear independently of each other according to Markov processes. Edge-Markovian graphs have been used to model, for example, intermittently-connected mobile networks [29, 30] and information spread thereon [31]. Taylor et al. derived in [32] the value of the epidemic threshold in edge-Markovian graphs and proposed control strategies to contain an epidemic outbreak, assuming homogeneous spreading and recovery rates.

It is worth remarking that, although most existing analyses of spreading processes over dynamic networks rely on the assumption that nodes present homogeneous rates and the Markov processes used to generate edge-switching signals are identical, these assumptions are not satisfied in real-world networks. Furthermore, most theoretical analyses mentioned above assume that edges appear and disappear in the network according to exponential distributions. However, as was found in several experimental studies [33, 34, 35] utilizing radio-frequency identification devices (RFID), the probability distribution of contact durations in human proximity networks is far from exponential. Therefore, it is of practical interest to develop tools to analyze spreading processes in TV networks with heterogenous rates and non-identical edge-switching signals able to replicate non-exponential switching patterns.

In this paper, we study spreading processes over TV networks with heterogeneous rates and non-identical edge-switching signals. First, we propose a wide class of TV networks where edges appear and disappear following aggregated Markov processes [36]. Since the class of aggregated Markov processes contains all Markov processes, the proposed class includes all edge-Markovian graphs [27]. We furthermore observe that, due to the modeling flexibility of aggregated Markov processes, the proposed class of TV networks can replicate human-contact durations following non-exponential distributions with an arbitrary accuracy. Second, we theoretically analyze the dynamics of spreading process over the proposed class of dynamic graphs. We model the dynamics of spreading processes on this class using a set of stochastic differential equations, which is a time-varying version of the popular -intertwined SIS model [9]. Third, we derive conditions for these stochastic differential equations to be globally exponentially stable around the infection-free equilibrium, i.e., the infection ‘dies out’ exponentially fast over time. One of the main challenges in this analysis is to derive computationally tractable stability conditions. For example, as will be illustrated later, if we simply apply a recently proposed stability condition [37] based on the direct application of Itô’s formula for jump processes (see, e.g., [38]), we obtain a stability condition that requires the computation of the largest eigenvalue of a matrix whose size grows exponentially with the number of edges in the network. Hence, this condition is hard (if not impossible) to verify for large-scale networks. In this paper, we use spectral graph theory [39] to derive stability conditions for the set of stochastic difference equations in terms of the largest eigenvalues of a matrix whose size grows linearly with the number of nodes in the network.

This paper is organized as follows. In Section 2, we introduce the notation and preliminary facts used in the paper and, then, introduce a flexible model of TV networks model called aggregated-Markovian random graph process. In Section 3, we state computationally tractable stability conditions for spreading processes taking place in this dynamic network model. In Section 4, we illustrate our results with some numerical examples. The proofs of the theorems are presented in Section 5.

2 Preliminaries and Problem Statement

This section starts by introducing notation in Subsection 2.1. Then, in Subsection 2.2, we review preliminary facts about spreading processes over static networks. In Subsection 2.3, we introduce the class of TV network models under consideration. We finalize this section by stating the problem to be studied in Subsection 2.4.

2.1 Notation

For a positive integer , define the set . For , define . Let denote the sign function over , which can be extended to matrices by entry-wise application. The Euclidean norm of is denoted by . Let denote the identity matrix. A square matrix is said to be Metzler if its off-diagonal entries are nonnegative. The spectral abscissa of a square matrix , denoted by , is defined as the maximum real part of its eigenvalues. We say that is Hurwitz stable if . Also, we define the matrix measure [40] of by . For a square random matrix , its expectation is denoted by and its variance is defined by . The diagonal matrix with diagonal elements , , is denoted by .

A directed graph is defined as the pair , where is a set of nodes111For simplicity in notation, we sometimes refer to node as . and is a set of directed edges (also called diedges or arcs), defined as ordered pairs of nodes. By convention, we say that is a directed edge from pointing towards . The adjacency matrix of a directed graph is defined as the matrix such that if , and otherwise.

A stochastic process taking values in a set is said to be an aggregated Markov process (see, e.g., [36]) if there exists a time-homogeneous Markov process taking values in a set  and a function such that . Throughout the paper, it is assumed that is finite, is surjective, and all the Markov processes are time-homogeneous. We say that a Markov process is irreducible if all the states in can be reached from any other state in (see, e.g., [41] for more details). We say that an aggregated Markov process  is irreducible if is irreducible. We will use the next lemma in our proofs.

Lemma 2.1 ([41])

Let be an irreducible Markov process taking values in a finite set . Then, has a unique stationary distribution . Moreover,

for every and .

2.2 Spreading Process over Static Networks

In this subsection, we review a model of spreading processes over static networks called the Heterogeneous Networked SIS (HeNeSIS) model, which is an extension of the popular -intertwined SIS model [9] to the case of nodes with heterogeneous rates. This model can be described using a continuous-time Markov process, as follows. Let be a directed graph, where nodes in represent individuals and diedges represent interactions between them (which we consider to be directed). At a given time , each node can be in one of two possible states: susceptible or infected. We define the variable as if node  is infected at time , and if is susceptible. In the HeNeSIS model, when a node  is infected, it can randomly transition to the susceptible state according to a Poisson process with rate , called the recovery rate of node . On the other hand, if node is in the susceptible state, it can randomly transition to the infected state according to a Poisson process with rate , where is called the infection rate of node . In other words, the rate of transition of node from susceptible to infected is proportional to the number of its infected neighbors.

Since the above Markov process has a total of possible states (two states per node), its analysis is very hard for arbitrary contact networks of large size. A popular approach to simplify the analysis of this type of Markov process is to consider upper-bounding linear models. Let denote the adjacency matrix of , and define , , , and . Then, it is known that the solutions (, , ) of the linear differential equation

(1)

upper-bounds the evolution of (, , ) from the exact Markov process with states [9]. Thus, if the dynamics in (1) is globally exponentially stable, the infection dies out exponentially fast in the exact Markov process [13, 6]. In the case of homogeneous infection and recovery rates, i.e., and for all , the dynamics (1) is globally exponentially stable if and only if [9]

(2)

Apart from the continuous-time Markov process described above, we can also model the dynamics of a spread using a discrete-time networked Markov chain with an exponential number of states [8]. As in the continuous-time case, we can show [8] that the solution of the linear difference equation

(3)

where , upper-bounds the probabilities (, , ). As in the continuous-time case, it turns out that if the discrete-time linear system in (3) is globally exponentially stable, the infection dies out exponentially fast in the exact Markov process [6]. In the case of homogeneous infection and recovery rates, the system in (3) is stable if and only if (2) holds [8].

2.3 Aggregated-Markovian Random Graph Processes

In this subsection, we introduce two new models of TV networks: the aggregated-Markovian edge-independent (AMEI) and the aggregated-Markovian arc-independent (AMAI) models. As we mentioned above, these models generalize the class of edge-Markovian time-varying networks [27]. Let be a stochastic graph process taking values in the set of directed graphs with nodes. Let denote the adjacency matrix of for each . The class of dynamic random graphs studied in this paper is defined as follows.

Definition 2.2

Consider a collection of aggregated Markov processes where are -valued functions and are stochastically independent Markov processes for . An aggregated-Markovian edge-independent (AMEI) network is a random graph process in which for and for all . If is irreducible for all possible pairs , then we say that is irreducible.

In other words, the upper-triangular entries of the time-varying adjacency matrix of an AMEI graph are independent aggregated Markov processes. The lower triangular entries are constructed via symmetry and the diagonal entries are zero. Hence, the graph is undirected by construction and does not contain self-loops. AMEI graphs extend the class of edge-Markovian graphs [27] and allow us to model a wider class of edge processes. For example, in an edge-Markovian graph, the time it takes for an edge to switch from connected to disconnected (or vice versa) follows an exponential distribution. In contrast, in an AMEI graph, we can design the function and the Markov process to fit any desired distribution for the contact durations with arbitrary precision, as illustrated in the following example:

Example 2.3

Let be the Markov process over the state space  presenting the state-transition diagram shown in Fig. 2.3. Consider an edge modeled by the aggregated Markov process , where the function is defined by and for all and . Therefore, the edge is active if and only if is in one of the states , , . From the diagram in Fig. 2.3, we see that the edge is active from the time enters until the time arrives to . It can be proved that the time elapsed between these two events follows a Coxian probability distribution, which forms a dense subset in the set of nonnegative random variables [42]. Therefore, by appropriately choosing the transition rates , , , , , , it is possible to tune the distribution followed by the time the edge is active in order to follow any distribution with an arbitrary accuracy (when is large enough). From the symmetry of the diagram in Fig. 2.3, we can also tune the parameters , , , , , so that the time the edge is inactive follows any given probability distribution with an arbitrary accuracy.

\xymatrix@!C *+=[Fo]c_1 \xy@@ix@

Fig. 1: A Markov process described in Example 2.3.

In the following definition, we introduce a directed version of AMEI graphs.

Definition 2.4

Consider a collection of aggregated Markov processes where are -valued functions and are stochastically independent Markov processes for and . An aggregated-Markovian arc-independent (AMAI) network is a random graph process in which for and for all . If is irreducible for all possible pairs , then we say that is irreducible.

Once the class of aggregated-Markovian networks is introduced, we analyze spreading processes taking on these random graph processes.

2.4 Problem Statement

In this paper, we address the following question: Under what conditions does a spreading process taking place in a TV network with aggregated-Markovian link dynamics die out exponentially fast? In other words, we consider the problem of analyzing the stability of the upper-bounding linear model (1) when the graph is an aggregated-Markovian network. In this case, the entries of the adjacency matrix are random processes and, therefore, (1) is the following stochastic linear differential equation:

where is a random matrix process. In the following section, we will derive conditions under which the disease-free equilibrium of is stable, i.e., the infection probabilities converge to zero as exponentially fast, in the following sense:

Definition 2.5

Consider a random graph process defined by either an AMEI or AMAI graph. We say that the disease-free equilibrium of is almost surely exponentially stable if there exists such that

for all initial states and . The supremum of satisfying the above condition is called the decay rate.

3 Stability of Spreading Processes in Random Graph Processes

In this section, we state the main results in this paper. In particular, we derive conditions for the spreading model to be stable when the network structure varies according to an AMEI (or AMAI) random graph process. In Subsection 3.1, we apply a stability condition [37] based on the direct application of the Itô formula for jump processes to derive stability conditions in terms of the largest eigenvalue of a matrix whose size depends exponentially on the number of edges in the graph. Due to the exponential size of this matrix, this condition is hard (if not impossible) to apply in the analysis of large-scale networks. Motivated by this limitation, we propose in Subsection 3.2 an alternative approach using tools from the theory of random matrices. In particular, we derive stability conditions in terms of the largest eigenvalue of a matrix whose size depends linearly on the number of nodes. In Subsection 3.3, we extend our results from continuous-time Markov processes to spreading processes modeled as discrete-time Markov chains.

3.1 Stability Conditions: Exponential Matrix Size

In this subsection, we show that the result in [37, Theorem 7.1] based on the direct application of the Itô formula for jump processes results in stability conditions that are not well-suited for large-scale graphs. We illustrate this idea through the analysis of the stability of the spread when is an AMAI random graph process. For simplicity in our exposition, we here temporarily assume that all the processes are Markovian; i.e., the mappings are the identities.

To state our claims, we need to introduce the following notations. Let be an AMAI random graph process (Definition 2.4). Let be the set of (time-varying) directed edges in (i.e., those directed edges such that is not the zero stochastic process). Let be the number of diedges in  and label the diedges using integers . Recall that these edges are not always present in the graph, but they appear (link is ‘on’) and disappear (link is ‘off’) according to a collection of aggregated-Markov processes. Therefore, at a particular time instant, only a subset of the edges in are present in the graph. The set of present edges can be represented by an element of , since we have possible edges that can be either ‘on’ of ‘off’. In other words, at a particular time, the contact network is one of the possible subgraphs of .

To analyze the resulting dynamics, we consider the set of possible subgraphs of and label them using integers . We define the one-to-one function , as the function that maps a particular subgraph of into the binary string in that indicates the subset of present edges in the subgraph. For the subgraph labeled , we denote by the corresponding binary sequence and by the -th entry in this sequence. In other words, means that the edge labeled in is present in the subgraph labeled . Furthermore, we denote by the adjacency matrix representing the structure of the subgraph labeled .

Given two labels , define the edit distance  to be the minimum number of links that must be added and/or removed from the subgraph labeled to construct the subgraph labeled , i.e., . If (i.e., the subgraphs differ in a single diedge), we define as the (only) index such that (i.e., the label of the single edge that must added/removed). Finally, assume that the transition probabilities of the Markov process () are given by

for some nonnegative constants and . We define and for each .

Under the above notations, the following stability condition readily follows from [37, Theorem 7.1]:

Proposition 3.1

Let be an AMAI random graph process. Define entry-wise, as follows. For ,

For , we let . Then, for every , the probability that is infected at time in the HeNeSIS model converges to zero exponentially fast as if the matrix

(4)

where denotes the Kronecker product of matrices, is Hurwitz stable.

As stated before, this stability condition is not applicable to large-scale networks, since it requires the computation of the eigenvalue with the largest real part of a matrix whose size increases exponentially with the number of edges in the graph . In the following subsection, we present alternative stability conditions that are better-suited for stability analysis of large networks.

3.2 Stability Analysis: Linear Matrix Size

In this subsection, we present sufficient conditions for almost sure exponential stability of the disease-free equilibrium of when the underlying time-varying network is described by an AMEI (or AMAI) random graph process. For clarity in our exposition, we state and discuss our main results in this section, leaving the details of the proofs for Subsections 5.1 and 5.2.

We now state our main results. Let us consider two positive constants and , and define the decreasing function by

(5)

where is the number of the nodes in the network. The following theorem provides conditions for a spreading process taking place in an AMAI random graph process to be almost surely exponentially stable.

Theorem 3.2

Let be an irreducible AMAI random graph process and define component-wise, as follows

(6)

Let , , and

Define

(7)

If

(8)

then the disease-free equilibrium of is almost surely exponentially stable with a decay rate greater than or equal to

where denotes the optimal value of in the maximization problem (7).

Proof:

See Subsection 5.1. \qed

Remark 3.3

The existence of the limit in (6) is guaranteed by the irreducibility of as will be discussed in Subsection 5.1, where we will also give an explicit representation of the limit.

The major computational cost for checking the stability condition in (8) comes from that of finding the matrix measure . Since the matrix has size , we can compute in operations using Lanczos’s algorithm [43], where and are the number of nodes and diedges in . In contrast, computing the dominant eigenvalue of the -dimensional matrix in (4) requires operations, since this matrix has nonzero entries. Therefore, the result in Theorem 3.2 (based on random graph-theoretical results) represents a major computational improvement with respect to the result in Proposition 3.1 (based on a direct application of the Itô formula for jump processes [37]).

Above, we have analyzed the stability of spreading processes in aggregated-Markovian arc-independent (AMAI) networks. In the rest of the subsection, we derive stability conditions for spreading processes taking place in aggregated-Markovian edge-independent (AMEI) networks. The next theorem is the AMEI counterpart of Theorem 3.2.

Theorem 3.4

Let be an irreducible AMEI random graph process. Let

(9)

Define

(10)

If

(11)

then the disease-free equilibrium of is almost surely exponentially stable with a decay rate greater than or equal to

(12)

where denotes the optimal value of in the maximization problem (10).

Proof:

See Subsection 5.2. \qed

Theorems 3.2 and 3.4 provide threshold conditions for stability of spreading processes in time-varying networks when agents present heterogeneous infection and recovery rates, and , respectively. In the particular case of networks with homogeneous rates, i.e., and for all , our results have an appealing interpretation that can be related with existing results in the literature. The following theorem states a threshold stability condition in the homogeneous case:

Theorem 3.5

Let be an irreducible AMEI random graph process. Assume that and for all . Let

(13)

Define

(14)

If

(15)

then the disease-free equilibrium of is almost surely exponentially stable with decay rate greater than or equal to

where denotes the optimal value of in the optimization problem (14). Moreover, if

(16)

then .

Proof:

See Subsection 5.2. \qed

Remark 3.6

Since entry-wise for every , the solution of is always upper-bounded by that of the deterministic differential equation , whose solution converges to zero as if and only if . Therefore, if the inequality (16) is false, then the disease-free equilibrium of is almost surely exponentially stable and, hence, we do not need to check condition (15). In other words, condition (16) guarantees that stability analysis of is not a trivial problem.

As shown in Subsection 2.2, the epidemic threshold in a static network of agents with homogeneous rates is given by . Condition (15) provides a similar epidemic threshold for time-varying networks in terms of the spectral abscissa of , which can be interpreted as an aggregated static network based on long-time averages. Furthermore, is a multiplicative factor that modifies the epidemic threshold corresponding to the aggregated static network. Below, we derive explicit expressions of for a few particular examples of time-varying graph models found in the literature:

\includegraphics

[width=3.5cm]SW_cropped.pdf

Fig. 2: Dynamic small-world network. Solid lines represent static edges, while dashed lines represent temporarily active edges modeled by aggregated-Markov processes.
Example 3.7 (Dynamic small-world networks)

We consider a time-varying version of the small-world network model studied in [44]. The network consists of nodes and static edges in the set  (see Fig. 2). Apart from these static edges, we also consider a collection of dynamic edges that switch on and off over time. In particular, we allow any other pair of nodes () to be dynamically connected according to an irreducible aggregated Markov process . For simplicity in our analysis, we assume that the stationary probability of an edge being active, defined by , equals a constant independent of and . In this case, the matrix  takes the form

The spectral abscissa of  is given by . Therefore, the stability condition (15) reads .

Example 3.8 (Edge-Markovian graph)

In this example, we consider the edge-Markovian graph model proposed in [27]. In this model, all the undirected edges  () are time-varying and modeled by independent -valued Markov processes  sharing the same activation rate  and de-activation rate . In this case, we have

and therefore . The stability condition (15) hence reads .

\includegraphics

[width=0.32]varyRatioAndDelta_n100

(a)      
\includegraphics

[width=0.32]varyRatioAndDelta_n1000

(b)        
\includegraphics

[width=0.32]varyRatioAndDelta_n10000

(c)          
Fig. 3: Plots of as and vary for the following cases: (a) and , (b) and , and (c) and .

In what follows, we elaborate on the behavior of the factor . In particular, we will analyze its dependence on defined in (13). Notice that  is zero when is either or  for all pairs , which corresponds to the case of a static, deterministic graph . In contrast, the maximum value of is achieved when , i.e., edges are either ‘on’ or ‘off’ with equal probability (asymptotically). Therefore, the term can be interpreted as a measure of structural variability, in the long run. In Fig. 3, we illustrate the behavior of the factor  as we modify the network size , the epidemic ratio , and the uncertainty measure . We observe that the smaller the value of , the larger the value of . In other words, as we reduce the structural variability (i.e., we reduce ), we obtain a better approximation using an aggregated static model (i.e., approaches ). Furthermore, notice that the value of tends to as increases. This indicates that the aggregated static network approximates the epidemic threshold more accurately as the network grows.

In this section, we have introduced our main theoretical results, namely, we have studied the epidemic threshold in the wide class of aggregated-Markovian random graph processes. Before we illustrate our results with some numerical simulations in Section 4, we will discuss the case when the epidemics is modeled as a discrete-time stochastic process.

3.3 Discrete-Time Epidemic Dynamics

In this subsection, we briefly present a discrete-time version of the stability analysis provided above. We define a discrete-time dynamic random graph as a stochastic process taking values in the set of directed graphs with nodes. In the discrete-time case, aggregated-Markovian edge- and arc-independent dynamic random graphs are defined in the same way as in the continuous-time case. Let us consider the discrete-time epidemic model in (3) taking place in the dynamics random graph . The resulting dynamics is described in the following stochastic difference equation:

where denotes the adjacency matrix of .

We define almost sure exponential stability of the disease-free equilibrium of in the same way as we did for the continuous-time case:

Definition 3.9

Let be a discrete-time edge-independent dynamic random graph. We say that the disease-free equilibrium of is almost surely exponentially stable if there exists such that

for all and (). We call the supremum of satisfying the above condition the decay rate.

For simplicity in our exposition, we focus only on the edge-independent case (AMEI graph process), since the analysis is identical for the arc-independent case. The next theorem provides a sufficient condition for almost sure exponential stability of the disease-free equilibrium of .

Theorem 3.10

Consider an irreducible and aperiodic222We say that an AMEI random graph process in discrete-time is aperiodic if all the aggregated Markov processes in the adjacency matrix of the graph are the images of aperiodic Markov chains. AMEI random graph process in discrete-time. Let and . Define

(17)

where is defined in (9). If

(18)

then the disease-free equilibrium of is almost surely exponentially stable with decay rate greater than or equal to

(19)

where denotes the value of giving the solution to the optimization problem (17).

Proof:

See Subsection 5.3. \qed

4 Numerical Simulations

Fig. 5: The averaged numbers of the infected nodes at time versus . Red dashed line: Epidemic threshold predicted by Theorem 3.10. Blue dotted line: Epidemic threshold predicted by using aggregated static network .
\includegraphics

[width=0.3]00060.pdf

(a)
\includegraphics

[width=0.3]00075.pdf

(b)
\includegraphics

[width=0.3]00090.pdf

(c)
Fig. 6: Sample paths of the number of infected nodes for different values of .
Fig. 4: The upper bound of decay rates. Red dashed line shows the epidemic threshold given by Theorem 3.10.  

In this section, we illustrate the effectiveness of our results with numerical simulations. In our illustrations, we consider a discrete-time edge-independent dynamic random graph with nodes. This random graph process is specified by the two following rules:

  1. Construct an undirected Erdős-Rényi random graph over the nodes of [45]. Denote the resulting adjacency matrix by .

  2. If , then nodes and are connected in via a stochastic link described by a -valued time-homogeneous Markov chain . The transition probability from to in this Markov process is given by ; while the transition probability from to is given by . If , then nodes and are not connected in at any time.

One can easily deduce that, according to the second rule, if . In the following simulations, we let and consider a realization of an the Erdős-Rényi graph with an edge probability of (in step 1 above). The values of the transition probabilities  (in step 2 above) are heterogeneous over the links of the Erdős-Rényi graph. In particular, we have chosen for a collection of random values from the Gaussian , where we have truncated those values to , and to . Finally, we let .

We consider a discrete-time spreading process with homogeneous infection and recovery rates, i.e., and for all . Using Theorem 3.10, we deduce that the disease-free equilibrium of is almost surely exponentially stable if . The upper bound of the decay rate versus is shown in Fig. 6. In order to check the accuracy of this epidemic threshold, we compare this threshold with the one obtained from the following simulation. We simulate the spreading process over the dynamic random graph described above for various values of . For each value of , we generate sample paths of the spreading process. In this simulation, we prevent a spreading process from dying out by re-infecting a randomly chosen node immediately after the infection process dies (i.e., when all the nodes become susceptible). A similar re-infection mechanism is employed for simulating spreading processes over static networks [46]. After the simulation, we compute the averaged number of infected nodes at time . This is the meta-stable number of infected nodes when re-infection of nodes is allowed. We then determine the metastable number of infected nodes in the original spreading process without reinfection as , where the subtraction of one compensates the effect of re-infection. Finally, we define the empirical epidemic threshold as the maximum value of such that .

Fig. 6 shows the value of as varies. We see that the above defined empirical threshold  lies between and , confirming that our condition indeed gives a sufficient condition for stability. On the other hand, the epidemic threshold based on the aggregated static network  (i.e., the supremum of such that ) is and, therefore, overestimates the actual epidemic threshold. In Fig. 6, we show multiple realizations of the sample paths used to obtain Fig. 6 for , , and .

5 Proofs of Main Results

This section presents the proofs of the theorems presented in Section 3.

5.1 Proof of Theorem 3.2

We begin by recalling a basic result on the stability analysis of a class of stochastic differential equations called switched linear systems. Let be an aggregated Markov process defined by the mapping and a Markov process  with state space . For each , there is an associated state matrix . An aggregated Markov jump linear system is defined by the following stochastic differential equation

(20)

where and . We remark that, if is the identity mapping, then the system in (20) is a Markov jump linear system [47]. We say that the system in (20) is almost surely exponentially stable if there exists such that for all and . The supremum of satisfying the above condition is called the decay rate. In order to prove Theorem 3.2, we will need the following criterion for almost sure exponential stability of aggregated Markov jump linear systems.

Lemma 5.1

Assume that is an irreducible Markov process and let be its unique stationary distribution. Assume that

(21)

Then, the aggregated Markov jump linear system in (20) is almost surely exponentially stable with a decay rate greater than or equal to .

Proof:

The above statement is known to be true for Markov jump linear systems [48, Theorem 4.2], i.e., when is the identity mapping. Let us now consider an arbitrary function . Define for each . Then, we see that the system in (20) is equivalent to the Markov jump linear system , for which we can guarantee almost sure exponential stability with decay rate grater than or equal to via the condition [48, Theorem 4.2]. This condition is equivalent to (21) by the definition of . \qed

Using Lemma 5.1, we can prove the following preliminary result, which will be useful in the proof of Theorem 3.2.

Proposition 5.2

Consider an irreducible AMAI random graph process . For all distinct and ordered pairs , let be an independent Bernoulli random variable with mean . Define the random matrix

(22)

where denotes the -matrix whose entries are all zero except its -entry. If , then, the infection-free equilibrium of is almost surely exponentially stable with decay rate greater than or equal to .

Proof:

Define the direct product , which is a stochastic process with state space . Define by . Also, for , define the matrix

(23)

Then, from Definition 2.4, we see that . Therefore, we can rewrite as the aggregated Markov jump linear system .

By the irreducibility of , each time-homogeneous Markov process has a unique stationary distribution on . Then, the unique stationary distribution of  is . Therefore, by Lemma 5.1, if , then is almost surely exponentially stable with decay rate greater than or equal to . Therefore, to complete the proof of the proposition, we need to show that

(24)

in the sense that the random matrices appearing in the both sides of the equation share the same probability distribution. From (23) and the definition of the matrix , we have

(25)

Since maps into , the random variable is the Bernoulli random variable with mean equal to

where we used Lemma 2.1. Moreover, all the random variables (, ) are independent of each other because is arc-independent. From these observations and the equation (25), we obtain (24), as desired. \qed

In order to complete the proof of Theorem 3.2, we need the following result from the theory of random graphs concerning the maximum eigenvalue of the sum of random symmetric matrices.

Proposition 5.3 ([49])

Let , , be independent random symmetric matrices. Let be a nonnegative constant such that for every with probability one. Also let . Then the sum satisfies

(26)

for every (where was defined in (5)).

Proof:

Let be arbitrary. Under the assumptions stated in the proposition, it is shown in [49] that, for every ,

where is defined by for . It is easy to show that the right hand side of the inequality takes its minimum value with respect to when . Substituting this particular value of into the inequality, we readily obtain (26). \qed

In order to complete the proof of Theorem 3.2, we also need the next lemma concerning Metzler matrices.

Lemma 5.4 ([50, Lemma 2])

Let and be Metzler matrices. If , then we have and .

We can now provide a proof of Theorem 3.2.

Proof:

Let and define . Notice that

(27)

because by (24). Also, from the definition of in (22), we have

where ( and ) are random, independent symmetric matrices. We then apply Proposition 5.3 to the random matrix , for which we choose . Clearly, the first constant term of has zero variance, while a simple computation shows that . Therefore

(28)

Combining (26), (27), and (28), we obtain the estimate

(29)

for .

Let