Autocorrelation of the susceptible-infected-susceptible process on networks.

In this paper, we focus on the autocorrelation of the susceptible-infected-susceptible (SIS) process on networks. The N-intertwined mean-field approximation (NIMFA) is applied to calculate the autocorrelation properties of the exact SIS process. We derive the autocorrelation of the infection state of each node and the fraction of infected nodes both in the steady and transient states as functions of the infection probabilities of nodes. Moreover, we show that the autocorrelation can be used to estimate the infection and curing rates of the SIS process. The theoretical results are compared with the simulation of the exact SIS process. Our work fully utilizes the potential of the mean-field method and shows that NIMFA can indeed capture the autocorrelation properties of the exact SIS process.

INTRODUCTION

The susceptible-infected-susceptible (SIS) process [1] is a basic epidemic model which models the spread of viruses, information, opinions, and computer malware on networks. In the SIS model, each node in the network can be either infected or susceptible (healthy). The infection state of node for at time is denoted by a Bernoulli random variable : infected or susceptible (healthy) . The SIS model has simple local rules that nodes can be infected by their infected neighbors and be cured by themselves. The infection and curing processes are independent and both Poisson processes with infection rate and curing rate , respectively. By tuning the effective infection rate , the phase transition of infection persistence emerges at an epidemic threshold determined by the network [2,3]. If the effective infection rate is below the epidemic threshold, then the virus dies out quickly and every node becomes healthy. Above the threshold, the infection can persist in the network for a very long time [4].

In this paper, we focus on the autocorrelation of the SIS process. Locally, an individual node in the network can be infected and cured repeatedly so that the infection state at two different time points can be autocorrelated. The autocorrelation of the infection state of a node between time and is The numerator on the right-hand side in (1) is the covariance of the infection state and , and the denominator normalizes the covariance. If time , then the infection state is fully correlated with itself, and the autocorrelation is . If and are independent, then the autocorrelation is . The autocorrelation is symmetric: . Given the initial infection state of the network , the infection states and are positively correlated [5], Corollary 1 such that and the autocorrelation . The autocorrelation contains one second-moment term , but the rest of the terms can be calculated given the first-moment infection probabilities and .

The autocorrelation contains information about the change of the infection state of each node. A large autocorrelation implies that the change of the infection state is slow, and the infection state is more likely to be identical between time and . While a smaller autocorrelation indicates that the infection state between time and is more independent. Globally, the fluctuating fraction of infected nodes is also autocorrelated, and the autocorrelation and its spectral analysis of in real epidemics can be traced back to Anderson et al. [6]. By analyzing the autocorrelation and its spectrum of the incidence data of pertussis, mumps, and measles, Anderson et al. [6] indicate statistically significant seasonal and the longer-term resurgence of those diseases and find that vaccination increases the periods of the longer-term oscillations of the incidence data. However, in the basic networked SIS model, the autocorrelation of the infection state is infeasible to be calculated, because the SIS model is a -state Markov process [3,7] and the computational complexity is exponentially high regarding of the network size . Previously, Meier et al. [8], Supplementary Information E analyzed the correlation of the infection state of the SIS model for small time intervals, but the calculation involves higher-order moments. In this paper, we apply the -intertwined mean-field approximation (NIMFA) [3] to study the autocorrelation of the infection state and the fraction of infected nodes both in the transient and steady states. Particularly in the steady state, we derive the explicit formula of the autocorrelation of the infection state, which is an exponentially decreasing function of time delay. The accuracy of the NIMFA autocorrelation is evaluated by simulating the exact SIS process. The result indicates that NIMFA, as an approximate stochastic process, well captures the autocorrelation properties of the exact SIS process. Moreover, the autocorrelation can also be used to estimate the infection and curing rates of the SIS process.

THE SIS PROCESS ON NETWORKS

The undirected and unweighted network with nodes is denoted by its symmetric adjacency matrix : If node and are connected, then ; otherwise, . The infection probability of node at time is just the expectation of its infection state , and the prevalence is defined as the expectation of the fraction of infected nodes .

The exact SIS process

The SIS process can be described by a Markov process, and there are states in total, including one all-healthy absorbing state [3,7]. The state transition of each node in the -state Markov process can be described as Since the all-healthy state is absorbing and the network is finite, the SIS process will enter the absorbing state when . However, the SIS process can also stay in the metastable state for a long time where the infection probability of every node is almost constant.

The infection probability of node for follows the governing equation [9], Equation (3) describes the exact Markovian SIS process, but higher-order moments of the infection states are involved in Eq. (3). In total, equations are needed to solve the process [9, p. 452] and the complexity increases exponentially with network size . Furthermore, the analysis of the SIS process is not tractable without approximation, not even for the complete graph [10].

The -intertwined mean-field approximation

NIMFA [11] approximates the exact Markovian SIS process by assuming independence , which is equivalent to approximating the infection rate due to all neighbors in (2) by its mean . For Bernoulli random variables, uncorrelation and independence are equivalent [12, footnote 5]. Under NIMFA, the governing equation is where is the NIMFA infection probability of node at time and approximates the exact infection probability . The NIMFA epidemic threshold is , where is the largest eigenvalue of the adjacency matrix . If the effective infection rate , then the infection can persist on the network and the steady-state infection probability is constant [3,13]. The steady state of NIMFA corresponds to the metastable state of the exact SIS process. If , then the NIMFA SIS process will eventually enter the all-healthy state . NIMFA has been successfully applied to analyze the first-order moments of the SIS process [3]. For example, the NIMFA infection probability and the prevalence well approximate the expectations of the infection state and the fraction of infected nodes , respectively. However, NIMFA has not yet been applied to approximate the autocorrelation properties. Since NIMFA omits the correlation between neighbors, the autocorrelation is the only second-moment property that is possibly captured by NIMFA.

To avoid ambiguity, we denote the NIMFA infection state of node at time by another Bernoulli random variable : infected and susceptible . Thus, we actually approximate the statistical properties of the infection state by those of in NIMFA. In the steady state and for , we denote the infection state of node by .

Under NIMFA, the transition of the infection state of node following Eq. (4) can be denoted by a two-state Markov process [14], and the transition rate of becomes a determined function of time. The whole system is composed of intertwined 2-state Markov processes instead of being a -state Markov process. Corresponding to (2), the transition of the NIMFA infection state is The infinitesimal generator of the Markov process (5) is

AUTOCORRELATION IN THE STEADY STATE

In the steady state, the NIMFA autocorrelation of the infection state of node with time lag is defined by where since is a Bernoulli random variable. By further derivation (see Appendix A), we obtain the autocorrelation as a function of the steady-state infection probability and the curing rate , where we assume that the time lag is positive without loss of generality. Since the autocorrelation is symmetric , for . The NIMFA infection probability in (8) can be obtained by solving the NIMFA Eq. (4) numerically.

With a fixed and time lag , the autocorrelation in (8) decreases with the infection rate because the infection probability increases correspondingly. A larger infection rate implies a faster state transition from healthy to infected, and the autocorrelation of the infection state is smaller consequently. A larger leads to a faster transition from infected to healthy, but, simultaneously, the infection probability of each neighbor becomes smaller. Therefore, the state transition of each node is slower from the healthy state to the infected state, and the effect of the curing rate is unclear. Only in special networks can the effect of the curing rate be determined. For example, the infection probabilities of all nodes are equal to in a -regular graph [3], and then the autocorrelation function becomes Formula (9) indicates that the autocorrelation of the infection state does not depend on the curing rate in regular graphs, which enables us to adjust the autocorrelation while keeping the effective infection rate unchanged. In regular graphs, the effect of the decrease (increase) of is exactly compensated by the increase (decrease) of in (8). The autocorrelation under other mean-field approximations can also be derived with the same procedure. For example, the heterogeneous mean-field approximation (HMF) assumes statistical equivalence among the nodes with the same degree [2], and the autocorrelation under HMF has the same form as the NIMFA autocorrelation (see Appendix B). In the case of regular graphs, HMF is equivalent [15] to NIMFA and then their approximate autocorrelations are identical.

Generally, the NIMFA infection probability of node with degree for is bounded by [3] in a connected network with minimum degree , and the NIMFA autocorrelation (8) is thus bounded by

The largest eigenvalue of the adjacency matrix follows , and then the effective infection rate can either be larger or smaller than when is above the threshold . Equation (8) indicates that the autocorrelation has another upper bound when (i.e., above the threshold). If , then and the upper bound (11) is tighter. If , then the upper bound in (10) is tighter, and we can rewrite (10) as In (12), the upper bound is just the product of the lower bound and the term . In a network with large degree deviation , the bound (12) is loose. In the regular graph, , and the upper bound achieves the exact NIMFA autocorrelation (9) while the lower bound does not.

In a heterogeneous network, e.g., the scale-free network, the degree can diverge in the thermodynamic limit . Thus, if and , then both the upper and lower bound in (12) converge to zero, and the autocorrelation . If and , then the lower bound converges to zero. Consequently, the autocorrelation is loosely bounded by .

From a global point of view, the fraction of infected nodes in the steady state can be approximated by . The autocorrelation of is just a linear combination of the autocorrelation of each node (see Appendix A),

AUTOCORRELATION IN THE TRANSIENT STATE

In this section, we consider the NIMFA autocorrelation of the SIS process at two arbitrary time points and , respectively. Different from that in the steady state in Sec. III, the infinitesimal generator (6) is a determined function of time given the initial state. The two-state Markov process (5) of each node is thus a time-inhomogeneous process. Calculating the process (5) allows us to analyze the autocorrelation of the epidemic process in the transient regime before the metastable state or the regime before the all-healthy steady state when the effective infection rate .

We denote the NIMFA autocorrelation of node between time and as Following a similar derivation as Eq. (13) in the steady state, the autocorrelation of the fraction of infected nodes is also a linear combination of the autocorrelation of each node,

Similarly to the steady-state autocorrelation in Sec. III, we only use the infection probabilities in the calculation, and the joint expectation in (14) becomes a crucial term. The calculation of the joint expectation involves the time-dependent transition matrix of which the element is . The computation of the autocorrelation functions (14) and (15), requires us first to calculate the matrix .

The matrix follows the time-inhomogeneous Kolmogorov forward equation where is the NIMFA infinitesimal generator (6). We can apply the Magnus expansion [16,17] to analyze the NIMFA transition matrix in Eq. (16). A brief introduction of the Magnus expansion can be found in Appendix C. Although the calculation of the exact NIMFA transition probability is not possible, approximations of allowing a fair comparison between NIMFA and the exact SIS process can be made with restricted error. First, there exists a matrix such that the solution of Eq. (16) is . Second, if (see the derivation of (C5) for details in Appendix C) then the exponent matrix can be expanded into a convergent Magnus series . Specifically, by only preserving the first term, i.e., , in the convergent Magnus series of , we can achieve a third-order accuracy (see Appendix C) for the time length , i.e., Equation (18) holds because holds for a matrix as can be verified by evaluating their power series. Using the Taylor expansion of the infinitesimal generator at time , the solution (18) becomes Only the first two terms of the Taylor expansion of the infinitesimal generator are preserved in (19) since the error is in (18). The first term on the right-hand side of (19) can be calculated by matrix diagonalization described in Appendix A. The derivative of the infinitesimal generator involves from Eq. (6), which is where denotes the neighbors of node . The calculation in Eq. (20) involves the infection probabilities of two-hop neighbors of node . Specifically, the transition probability that node remains infected after time units is Different from that in the steady state [see Eq. (8)], the infection probabilities of neighbors of node always appear in the calculation of the transition matrix in the transient state as indicated in (20). Higher-order accuracy is also possible by preserving more terms of the Magnus series, and higher-order derivative , which can be calculated by the infection probabilities of all nodes within hops from node , is involved. For example, if we preserve the second term in the Magnus expansion of , which can be calculated by the Taylor expansion as then we can achieve an accuracy of because (see Appendix C) and the calculation involves the infection probabilities of neighbors within three hops. For NIMFA, preserving more terms is not always reasonable, because the infection probability of each node can only be solved numerically. When more Magnus terms are preserved, the inaccuracy is mainly caused by the numerical method which solves the nonlinear NIMFA Eq. (4). For example, using the fourth-order Runge-Kutta method [17, p. 200], the error of the infection probabilities is of order .

For a time interval , the Magnus expansion of the exponent may not converge. The time interval can be divided into subintervals with length in which the Magnus series converges. The NIMFA transition matrix between time and can be written as by the Chapman-Kolmogorov equation [see Eq. (C2)]. Equation (21) is also applicable to a small time interval to obtain a more accurate result. An -order accuracy regarding the time delay is achieved for the transition matrix using Eq. (21) if the accuracy is for each .

The analysis in this section allows us to calculate and compare the NIMFA autocorrelation with the exact SIS process since the error can be controlled, even though the exact NIMFA autocorrelation is not feasible in the transient state.

NUMERICAL AND SIMULATION RESULTS

In this section, we compare the NIMFA autocorrelation with the autocorrelation of the exact SIS process from the simulation. The simulation of the exact SIS process is implemented by the Gillespie algorithm (Monte Carlo method) [18-20] and the theoretical results are obtained by solving the NIMFA Eq. (4) numerically (fourth-order Runge-Kutta method [17, p. 200]). In the steady state, we run the simulation for time units with the curing rate and sample the infection state of each node every 0.001 time unit. In other words, we obtain the infection state for from simulation. We only use the state sequence sampled after to ensure that the SIS process is in the metastable state. Moreover, the time series of the fraction of infected nodes can be calculated as . In the transient state, realizations of the infection states and are obtained to calculate the autocorrelation between two arbitrary time and .

Steady state

Figures 1 to 3 show the NIMFA autocorrelation and the simulated autocorrelation of the infection state of randomly selected nodes in an Erdős-Rényi (ER) graph, a regular graph with degree 26, and a star graph, respectively. The NIMFA autocorrelation is a very accurate approximation on those graphs. Figure 1 shows that the autocorrelation of the infection state is not sensitive to the value of the curing rate , which is reasonable because the deviation of the degree is small and the result is similar to that of the regular graph in Fig. 2. In Fig. 2, the autocorrelation of the infection state is identical to formula (9) that the autocorrelation is invariant to the curing rate in regular graphs. Figure 3 shows the autocorrelation of the infection state in a star graph. The autocorrelation of the hub node is much smaller than that of the leaf nodes since the infection probability of the hub node is larger. The cross correlation of the infection states between neighbors shown in Fig. 1 to 3 is approximately 0, which leads to the effectiveness of NIMFA since NIMFA omits the cross correlation between neighbors.

FIG. 1.

A randomly selected node is evaluated in an Erdős-Rényi (ER) network with the link connecting probability 0.4 and . The autocorrelation is approximately constant for different value of . The cross correlation between the node with one of its neighbors is also plotted, which is almost zero.

FIG. 3.

The autocorrelation of the infection state of the hub and a leaf node in a star graph with . The NIMFA autocorrelation shows a very good approximation and the cross correlation between hub and leaf nodes is approximately 0.

FIG. 2.

The autocorrelation of the infection state of a 26-regular graph with . The results are similar to those of the ER graph. The autocorrelation is invariant to .

Figure 4 shows the autocorrelation of the infection state of a node in a cycle graph and NIMFA fails to capture the autocorrelation. Actually, NIMFA also fails to approximate the prevalence as shown in Fig. 4. In the situation of the cycle graph, the cross correlation of the infection states between neighbors is much larger than zero and NIMFA itself is a bad approximation. The accuracy of mean-field methods has been studied in Refs. [21-23], which is beyond the scope of this paper.

FIG. 4.

The correlation of the infection state of a node and the prevalence in a cycle graph with . Initially all nodes are infected to prevent the inaccuracy caused by early die-out [24]. The NIMFA autocorrelation is much smaller than the exact one, and the cross correlations between neighbors and second-hop neighbors are very large.

We also calculate the autocorrelation of the fraction of infected nodes . Figure 5 shows that NIMFA can also approximate the autocorrelation of the fraction of infected nodes in the star graph corresponding to Fig. 3.

FIG. 5.

The autocorrelation of the fraction of infected nodes in the metastable state.

Transient state

In the transient state, we validate the NIMFA autocorrelation on the star graph where the NIMFA infection probabilities are accurate while nodes have very different degrees.

Figure 6 shows the joint expectation of the infection states and the corresponding NIMFA approximation of the leaf and hub nodes. For the leaf node and the hub node, the convergent time delay of the Magnus series of are and from (17), respectively. Figure 6 indicates that the NIMFA joint expectation (the blue lower curve) is accurate comparing with the exact joint expectation for a small time delay , i.e., for the leaf node. For a large time delay, the inaccuracy is due to either the omission of term in (19) or that the NIMFA transition probability matrix itself is a bad approximation, but we can eliminate the possibility of the latter using Eq. (21). As the black middle curve in Fig. 6 indicated, the NIMFA joint expectation is indeed a good approximation using Eq. (21) with subinterval length 0.01.

FIG. 6.

The joint expectation of the infection state and the corresponding NIMFA approximation of the SIS process on the star graph.

From a global point of view of the network, Fig. 7 presents the autocorrelation of the fraction of infected nodes and the corresponding NIMFA approximation , which are in the transient state of the SIS process before the metastable state. The exact autocorrelation is well fitted by NIMFA. Interestingly, the decay of the autocorrelation in the transient state is also exponential as shown in Fig. 7, but we cannot demonstrate exponential decay as opposed to the steady state.

FIG. 7.

The autocorrelation of the fraction of infected nodes and the corresponding NIMFA approximation of the SIS process on the star graph.

In this section, we have tested our method on different networks with size 50, but for larger networks, the results are similar. In a conclusion, NIMFA captures the autocorrelation properties of the exact SIS process except in the cases that NIMFA is not applicable even for approximating the first-moment properties, i.e., the infection probabilities and the prevalence .

ESTIMATING THE CURING RATE AND THE INFECTION RATE : AN APPLICATION

In real epidemics, a disease agency may have the infection-state data by monitoring individuals periodically but no information about the rates. We consider the reverse problem of estimating the curing rate and the infection rate , given the sequence of the infection state of node in the metastable state. From Eq. (8), the curing rate is Formula (22) can be used to estimate the curing rate of the SIS process. In formula (22), we can approximate the infection probability as , while the autocorrelation , which approximates the exact autocorrelation in (1), is just the autocorrelation of the binary infection sequence . Furthermore, using the NIMFA equation in the metastable state , we can eliminate and (22) becomes Under NIMFA, the curing rate can be estimated by (22) without knowing the underlying network. However, to estimate the infection rate , formula (23) involves the network information. We rewrite (23) as and sum over all nodes where is the degree of node . After rearrangement of the above equation, we obtain Thus, the estimation of the infection rate requires either the degree of every node for all as in (24) or the local topology information about node , i.e., for all as in (23).

Using the binary infection-state sequence obtained by simulation, we estimate the curing rate and the infection rate by (22) and (23), respectively. In Fig. 8, the value of the estimated rates times the time lag is plotted for a leaf node of the star graph corresponding to Fig. 3. The slopes of the linear fitting functions (red curves in Fig. 8) are the estimated rates, and both the estimated infection rate and the curing rate are 1.00 while both the real rates equal to 1.

FIG. 8.

The estimation of the infection rate and the curing rate using (23) and (22) for the star graph corresponding to Fig. 3. The curves are and of a leaf node versus . Both the estimated and are 1.00, while the real values of rates are both 1.

CONCLUSION

In this paper, we study the autocorrelation, the only second-moment property captured by NIMFA, of the SIS process. We obtained the explicit formula of the autocorrelation, i.e., Eq. (8), under NIMFA in the steady state, and the steady-state autocorrelation follows an exponential decay with the time lag. Interestingly, the steady-state autocorrelation is independent of the curing rate in regular graphs. Moreover, using the Magnus expansion, we are able to calculate the autocorrelation in the transient state of the SIS process. Our analysis of the transient state not only allows the study of the SIS process above or below the epidemic threshold but also opens an avenue for the study of the critical behavior [25].

We also evaluated our results by simulation. Although NIMFA assumes that there is no correlation between the infection states of neighbors, i.e., for , we show that the NIMFA autocorrelation () is generally accurate by simulation, and the accuracy depends on the accuracy of the NIMFA infection probabilities. If NIMFA can capture the first-order moments, i.e., the infection probability of each node and the prevalence, under certain SIS parameters and networks, then NIMFA can also be applied to approximate the autocorrelation properties. Finally, we show that our results can be used to estimate the infection and curing rate of the SIS process.