Stability of Spreading Processes over TimeVarying LargeScale Networks
Abstract
In this paper, we analyze the dynamics of spreading processes taking place over timevarying networks. A common approach to model timevarying networks is via Markovian random graph processes. This modeling approach presents the following limitation: Markovian random graphs can only replicate switching patterns with exponential interswitching 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 interswitching 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.
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 timeinvariant (TI) networks have been made during the last decade (see [4] and [5] for recent surveys). Under the timeinvariance assumption, we find several recent breakthroughs in the literature, such as the rigorous analysis of meanfield 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 timevarying (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 wellsuited 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 activitydriven 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 timevarying 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 edgeMarkovian graphs, was proposed in [27] and analyzed in [28]. In this model, edges appear and disappear independently of each other according to Markov processes. EdgeMarkovian graphs have been used to model, for example, intermittentlyconnected mobile networks [29, 30] and information spread thereon [31]. Taylor et al. derived in [32] the value of the epidemic threshold in edgeMarkovian 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 edgeswitching signals are identical, these assumptions are not satisfied in realworld 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 radiofrequency 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 nonidentical edgeswitching signals able to replicate nonexponential switching patterns.
In this paper, we study spreading processes over TV networks with heterogeneous rates and nonidentical edgeswitching 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 edgeMarkovian graphs [27]. We furthermore observe that, due to the modeling flexibility of aggregated Markov processes, the proposed class of TV networks can replicate humancontact durations following nonexponential 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 timevarying version of the popular intertwined SIS model [9]. Third, we derive conditions for these stochastic differential equations to be globally exponentially stable around the infectionfree 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 largescale 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 aggregatedMarkovian 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 entrywise application. The Euclidean norm of is denoted by . Let denote the identity matrix. A square matrix is said to be Metzler if its offdiagonal 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 nodes^{1}^{1}1For 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 timehomogeneous 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 timehomogeneous. 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 continuoustime 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 upperbounding linear models. Let denote the adjacency matrix of , and define , , , and . Then, it is known that the solutions (, , ) of the linear differential equation
(1) 
upperbounds 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 continuoustime Markov process described above, we can also model the dynamics of a spread using a discretetime networked Markov chain with an exponential number of states [8]. As in the continuoustime case, we can show [8] that the solution of the linear difference equation
(3) 
where , upperbounds the probabilities (, , ). As in the continuoustime case, it turns out that if the discretetime 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 AggregatedMarkovian Random Graph Processes
In this subsection, we introduce two new models of TV networks: the aggregatedMarkovian edgeindependent (AMEI) and the aggregatedMarkovian arcindependent (AMAI) models. As we mentioned above, these models generalize the class of edgeMarkovian timevarying 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 aggregatedMarkovian edgeindependent (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 uppertriangular entries of the timevarying 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 selfloops. AMEI graphs extend the class of edgeMarkovian graphs [27] and allow us to model a wider class of edge processes. For example, in an edgeMarkovian 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 statetransition 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@
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 aggregatedMarkovian arcindependent (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 aggregatedMarkovian 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 aggregatedMarkovian link dynamics die out exponentially fast? In other words, we consider the problem of analyzing the stability of the upperbounding linear model (1) when the graph is an aggregatedMarkovian 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 diseasefree 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 diseasefree 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 largescale 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 continuoustime Markov processes to spreading processes modeled as discretetime 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 wellsuited for largescale 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 (timevarying) 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 aggregatedMarkov 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 onetoone 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 entrywise, 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 largescale 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 bettersuited 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 diseasefree equilibrium of when the underlying timevarying 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 componentwise, as follows
(6) 
Let , , and
Define
(7) 
If
(8) 
then the diseasefree 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 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 graphtheoretical 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 aggregatedMarkovian arcindependent (AMAI) networks. In the rest of the subsection, we derive stability conditions for spreading processes taking place in aggregatedMarkovian edgeindependent (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 diseasefree 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 timevarying 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 diseasefree 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 entrywise for every , the solution of is always upperbounded 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 diseasefree 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 timevarying networks in terms of the spectral abscissa of , which can be interpreted as an aggregated static network based on longtime 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 timevarying graph models found in the literature:
Example 3.7 (Dynamic smallworld networks)
We consider a timevarying version of the smallworld 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 (EdgeMarkovian graph)
In this example, we consider the edgeMarkovian graph model proposed in [27]. In this model, all the undirected edges () are timevarying and modeled by independent valued Markov processes sharing the same activation rate and deactivation rate . In this case, we have
and therefore . The stability condition (15) hence reads .
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 aggregatedMarkovian 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 discretetime stochastic process.
3.3 DiscreteTime Epidemic Dynamics
In this subsection, we briefly present a discretetime version of the stability analysis provided above. We define a discretetime dynamic random graph as a stochastic process taking values in the set of directed graphs with nodes. In the discretetime case, aggregatedMarkovian edge and arcindependent dynamic random graphs are defined in the same way as in the continuoustime case. Let us consider the discretetime 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 diseasefree equilibrium of in the same way as we did for the continuoustime case:
Definition 3.9
Let be a discretetime edgeindependent dynamic random graph. We say that the diseasefree 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 edgeindependent case (AMEI graph process), since the analysis is identical for the arcindependent case. The next theorem provides a sufficient condition for almost sure exponential stability of the diseasefree equilibrium of .
Theorem 3.10
Consider an irreducible and aperiodic^{2}^{2}2We say that an AMEI random graph process in discretetime 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 discretetime. Let and . Define
(17) 
where is defined in (9). If
(18) 
then the diseasefree 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
In this section, we illustrate the effectiveness of our results with numerical simulations. In our illustrations, we consider a discretetime edgeindependent dynamic random graph with nodes. This random graph process is specified by the two following rules:

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

If , then nodes and are connected in via a stochastic link described by a valued timehomogeneous 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ősRé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ősRé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 discretetime spreading process with homogeneous infection and recovery rates, i.e., and for all . Using Theorem 3.10, we deduce that the diseasefree 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 reinfecting a randomly chosen node immediately after the infection process dies (i.e., when all the nodes become susceptible). A similar reinfection 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 metastable number of infected nodes when reinfection 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 reinfection. 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 infectionfree 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 timehomogeneous 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 arcindependent. 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