Statistical Properties of SIS Processes with Heterogeneous Nodal Recovery Rates in Networks

Abstract

The modeling and analysis of epidemic processes in networks have attracted much attention over the past few decades. A major underlying assumption is that the recovery process and infection process are homogeneous, allowing the associated theoretical studies to be conducted in a convenient manner. However, the recovery and infection processes usually exhibit heterogeneous rates in the real world, which makes it difficult to characterize the general relations between the dynamics and the underlying network structure. In this work, we focus on the susceptible–infected–susceptible (SIS) epidemic process with heterogeneous recovery rates in a finite-size complete graph. Specifically, we study the metastable-state statistical properties of SIS epidemic dynamics with two different nodal recovery rates in complete graphs. We propose approximate solutions to the metastable-state expectation and the variance in the number of infected nodes within the framework of the mean-field approximation method. We also derive several upper and lower bounds of the steady-state probability that a node is in the infected state. We verify the proposed approximate solutions of the mean and variance via simulations. This work provides insights into the fluctuations in the statistical properties of epidemic processes with complex dynamical behaviors in networks.

1. Introduction

The modeling and analysis of epidemic-like dynamical processes in complex networks have gained much attention among researchers over the past few decades [1,2,3,4,5,6]. From the view point of dynamics, two basic models are commonly employed to describe epidemic spread in networks. They are named the susceptible–infectious–susceptible (SIS) model and the susceptible–infectious–recovery (SIR) model. It is convenient to characterize the change in dynamics in populations by means of the notation named the system state, which can be defined as the collection of all nodal infection states as well as the network structure. Many efforts have been made to investigate how the system parameters impact the system state. Specifically, it is of great interest to investigate how many infected nodes there when the system state changes over a long time period. It is also intriguing to gain insights into the relations between the number of infected nodes and the system parameters such as the infection rate and the network structure. Pastor-Satorras et al. [2] reviewed related works about the epidemic spread in systems with static or dynamical single-layer network structures and provided insights into the impact of network structure on the steady-state system state. Bianconi [3] provided a thorough overview of multilayer networks, which denotes the structure of a network of networks, and reported how epidemic spread is affected by the structure of a multilayer network. Usually, it is possible to determine the epidemic threshold by the eigenstructure of the topology of a finite-size network [2,3]. It is usually assumed that the infection process and the recovery process that constitutes the epidemic processes are Markovian processes, as studied in most previous works. Tomovski et al. [7] and Basnarkov et al. [8] investigated non-Markovian SIS and SIR epidemic processes in networks, respectively. Furthermore, Han et al. [9] extended the results for the Markovian case to non-Markovian epidemics in networks with time-changing topology.

From the view point of methodology, several theoretical frameworks have been proposed to investigate the dynamical statistical properties of the processes and obtain physical insights into the impact of the network’s structure on the dynamics. Typical methodologies and frameworks include the mean-field approximation [2,10,11], message passing [12], master equation [13] and Markov theory [14]. Markov theory is capable of exactly capturing the time-evolution change of each system state but is limited to small-size networks because of its computational complexity [12]. The mean-field approximation framework appears in various forms with respect to the abstraction of the system. One merit of the mean-field approximation framework is that it is possible to strike a balance between the simplicity and accuracy of the theory.

A common underlying assumption made with the basic SIS model or the SIR model is that all nodal infection rates are equal and that the same applies to the nodal recovery rates. This simplicity of the model makes it possible to theoretically and experimentally study the epidemic process and obtain several intriguing insights, such as the relationship between the network’s structure and the epidemic threshold, in a convenient manner. However, a real-world epidemic process usually shows heterogeneous nodal recovery and nodal infection rates, which makes it difficult to propose a theory suited to any general scenario that would enable us to obtain physical insights into its general dynamics. Compared with studies on the basic SIS processes in complex networks, studies on SIS processes with heterogeneous nodal recovery and infection rates are still limited.

In this work, we focus on the SIS process with heterogeneous recovery rates in a finite-size network that was previously used to characterize some real-world scenarios of interest, such as the impact of a mixture of mitigation methods on interest flooding attacks in named data networking (NDN) [15]. Specifically, we employ an individual-based mean-field approximation method to model and analyze the SIS epidemic process with two types of nodal recovery rates and an identical infection rate. We derive the approximate solutions to the metastable-state mean value and the variance value of the infected nodes. These approximations are capable of characterizing the metastable state of the dynamics in networks as well as the impact of the network structure on epidemic dynamics. We also derive some bounds of the metastable-state probability of a node being in the infected state. Then, we experimentally verify the proposed theoretical models.This work is a first for predicting variance for the scenario of interest and could keep a balance between theoretical complexity and computational complexity, though it has its scope and accuracy limitations.

This paper is organized as follows. Section 2 provides a literature review, and the methodology is described in Section 3. Section 4 presents the simulation results and corresponding theoretical results. We discuss the work and present our conclusions in Section 5 and Section 6, respectively.

2. Literature Review

The SIS processes with homogeneous and fixed nodal infection rates and recovery rates in a static and single-layer network have been analyzed thoroughly in many previous studies. Pastor-Satorras et al. [1] modified the basic compartmental model and proposed a heterogeneous mean-field approximation model to investigate the impact of the statistical properties of large networks on the dynamical behaviors of the SIS process in a coarse-grained manner. Van Mieghem [16] proposed the N-intertwined mean-field approximation (NIMFA), which is a time-continuous individual-based mean-field model for finite-size networks. The metastable state of the homogeneous SIS process in finite-size networks was claimed to be approximated by the steady state of the NIMFA model.

Several steady-state statistical properties of the dynamics were also verified under the theoretical framework of NIMFA. For instance, the epidemic threshold can be approximated by the inverse of the maximum eigenvalue of the adjacency matrix of the underlying network. The metastable-state mean fraction of infected nodes can also be calculated via NIMFA. Chakrabarti et al. [17] proposed a discrete-time model to study the statistical properties of the SIS process in networks. Omic and Van Mieghem [18] applied the NIMFA framework to investigate the metastable-state variance of the number of infected nodes. Several general approximations were derived to calculate the variance for the SIS process with homogeneous infection and recovery rates in complete and bipartite graphs.

It is of great interest to investigate the heterogeneous SIS process, which is defined as the SIS process where at least one nodal infection rate or recovery rate is not the same as that of the others. However, the modeling of SIS processes with heterogeneous recovery rates or infection rates has attracted limited attention. Van Mieghem and Omic [19] investigated the SIS process with heterogeneous virus-spread behaviors in networks with finite sizes and extended the NIMFA model to a heterogeneous NIMFA (HNIMFA) model to analyze this complex type of scenario in a general manner. Similarly, Jiao et al. [20] extended the NIMFA framework to accommodate the SIS process with two different infection rates and the same recovery rates in finite-size bipartite and star graphs. Jiao et al. [20] derived the expressions to estimate the metastable-state mean value of the infection fraction and experimentally verified the accuracy of the approximations. Cui et al. [21] proposed a susceptible–infected–water–susceptible model with different infection rates, which describes the disease spread in different mediums. Pagliara and Leonard [22] proposed a susceptible–infected–recovered–infected model where the nodal susceptibility could change after recovery.

In other research regarding SIS process modeling, Yuan et al. [23] investigated the impact of time-varying infection rates on infection behaviors. Wang et al. [15] applied the HNIMFA framework to characterize the scenario of interest flooding attacks and mitigation methods. The SIS model with various recovery rates was employed to estimate the metastable-state mean value of the infection fraction and assess the impact of a mixture of mitigation methods on the infection fraction. Jing et al. [24] proposed an approximate SIS pairwise model to investigate the variances of system state variables and the temporal variance of the infection. The authors also found that the variances of the pairwise model and the mean-field SIS model are nearly equal for dense networks. Finally, Guo et al. [25] employed the HNIMFA framework to estimate the metastable-state average variance of the number of infected nodes for bipartite and star graphs and verified its accuracy.

It is of great interest to investigate the statistical properties of the SIS epidemic dynamics with complex heterogeneous infection and recovery behaviors. The obtained insights into the complex epidemic dynamics will improve understanding and control of real-world epidemics.

3. Methodology

3.1. Notation and Model

We first introduce notation before describing the general HNIMFA model. The dynamics of interest are those of the SIS process with heterogeneous nodal infection and recovery rates in a network represented by the undirected graph $G(N,L)$, composed of L links and N nodes. The topology of the network (graph) is denoted by a symmetric $N\times N$ adjacency matrix $A:=\left\{a_{ij}\right|0,1\}$, where $a_{ij}$ is equal to 1 if there is a link between nodes i and j; and 0 otherwise. We denote by $d_i$ the degree of the node i, such that $d_i={\sum }_{j=1}^N{a}_{ij}={\sum }_{i=1}^N{a}_{ij}$ for the undirected graph.

Regarding the SIS process in a network as a system, the state of the system can be specified by the viral state of each node and the topology. At any time, the viral state of a node can only be one of two states: infected or susceptible. In order to specify the viral state of a node i, the Bernoulli random variable $X_i\left(t\right)\in \{0,1\}$ is introduced. Specifically, the random variable $X_i\left(t\right)=1$ means node i is in the infection state at time t. The random variable $X_i\left(t\right)=0$ means node i is in the susceptible state at time t. The probability of node i being in the infected or susceptible state is represented by $v_i\left(t\right)=Pr[X_i\left(t\right)=1]$ or $1-v_i\left(t\right)=Pr\left[X_i\left(t\right)\right]=0$, respectively. Using the properties of the Bernoulli random variable provides the expectation and variance of the nodal viral state $X_i$, which are

(1)$\begin{array}{ccc}& E\left[X_i\right]& =1·Pr[X_i=1]+0·Pr[X_i=0]=v_i,\hfill \end{array}$

(2)$\begin{array}{ccc}& \mathrm{Var}\left[X_i\right]& =(1-v_i)v_i.\hfill \end{array}$

Figure 1 provides instances of two Poisson processes constituting the SIS process in a network of a pairwise type. Any process between a nodal pair is independent of others. We consider a scenario where node i in the infected state is connected with node j in the susceptible state, as shown in Figure 1a. There exists an infection process between node i and node j, where the infected node i tries to infect the susceptible node j. This infection is a Poisson process with a rate of ${\beta }_i$, such that node j will probably become infected soon, changing its state to susceptible. On the other hand, any infected node will recover and change to a susceptible state, which is also characterized by a Poisson process, as shown in Figure 1b. Node i in the infected state will recover following a Poisson process at a rate of ${\delta }_i$. All the infection and recovery processes simultaneously existing in a finite-size network of interest constitute a heterogeneous SIS process. Additionally, we introduce the notation ${\tau }_i={\beta }_i/{\delta }_i$, which is named the effective infection rate and used to simplify the differential equations of the dynamics in the following.

The general HNIMFA model is elaborated in the following. One can depict the SIS process with heterogeneous infection and recovery rates in a network using a Markov chain. The exact Markov differential equation for the viral state of node i is expressed as follows:

(3)$\begin{array}{cc}\hfill {\displaystyle \frac{dE\left[X_i\left(t\right)\right]}{dt}}& =E\left[-{\delta }_i{X}_i+\sum _{j=1}^N{a}_{ij}{\beta }_j{X}_j(1-X_i)\right].\hfill \end{array}$

The first-order NIMFA model [16] employs the assumption of independence approximation:

(4)$\begin{array}{c}\hfill Pr[X_i=1,X_j=1]=Pr[X_i=1]·Pr[X_j=1].\end{array}$

In other words, the approximation means that the random variables $X_i$ and $X_j$ are assumed to be independent. Using the property of the Bernoulli random variable (1) and assuming independence, we arrive at

(5)$\begin{array}{c}\hfill E[X_i=1,X_j=1]=E[X_i=1]·E[X_j=1]=v_i{v}_j.\end{array}$

Using the approximation (5), the Markov differential Equation (3) turns out to be

(6)$\begin{array}{cc}\hfill {\displaystyle \frac{d{v}_i\left(t\right)}{dt}}& =-{\delta }_i{v}_i+\sum _{j=1}^N{a}_{ij}{\beta }_j{v}_j(1-v_i)\hfill \\ \hfill & =\sum _{j=1}^N{a}_{ij}{\beta }_j{v}_j-v_i\left(\sum _{j=1}^N{a}_{ij}{\beta }_j{v}_j+{\delta }_i\right).\hfill \end{array}$

The Markov differential Equations (3) and (6) were proposed and derived in ([19], (2) on p. 2). Employing the assumption of independence results in the overestimation of $d{v}_i/dt$ [16]. Similarly, we can obtain the differential equations for all the nodal viral states:

(7)$\begin{array}{cc}\hfill {\displaystyle \frac{d{v}_1\left(t\right)}{dt}}& =\sum _{j=1}^N{\beta }_j{a}_{j1}{v}_j\left(t\right)-v_1\left(t\right)\left(\sum _{j=1}^N{\beta }_j{a}_{j1}{v}_j\left(t\right)+{\delta }_1\right),\hfill \\ \hfill & \cdots \hfill \\ \hfill {\displaystyle \frac{d{v}_N\left(t\right)}{dt}}& =\sum _{j=1}^N{\beta }_j{a}_{jN}{v}_j\left(t\right)-v_N\left(t\right)\left(\sum _{j=1}^N{\beta }_j{a}_{jN}{v}_j\left(t\right)+{\delta }_N\right).\hfill \end{array}$

The steady state of a system composed of the SIS process in a given finite-size network actually corresponds to the metastable state [16]. The steady state is denoted by $v_{i\infty }$:

(8)$\begin{array}{c}\hfill {\left.{\displaystyle \frac{d{v}_i\left(t\right)}{dt}}\right|}_{t\to \infty }=0.\end{array}$

For the steady-state regime, the differential Equation (6) for node i’s infection probability can be rewritten as

(9)$\begin{array}{c}\hfill {\delta }_i{v}_i=\sum _{j=1}^N{a}_{ij}{\beta }_j{v}_j(1-v_i).\end{array}$

For the steady-state regime, applying some manipulations, the nodal infection probability (9) can be expressed as

(10)$\begin{array}{c}\hfill v_i=1-{\displaystyle \frac{1}{1+{\tau }_i{\sum }_{j=1}^N{a}_{ji}{v}_j}}.\end{array}$

We also denote by y or $y_{\infty }$ and I or $I_{\infty }$ the metastable-state fraction of infected nodes and the metastable-state number of infected nodes, respectively.

The mean value $E\left[I\right]$ of the metastable-state number of infected nodes can be expressed as the sum of the expectation of the probability of all the nodes, according to the statistical property of the random variable following a binomial distribution composed of Bernoulli random variables. Thus, we arrive at

(11)$E\left[I\left(t\right)\right]=\sum _{i=1}^{N}E\left[X_i\right]=\sum _{i=1}^N{v}_i.$

One exact expression of the variance $\mathrm{Var}\left[I\right(t\left)\right]$ of the metastable-state number I of infected nodes can be expressed as

(12)$\mathrm{Var}\left[I\left(t\right)\right]=\sum _{i=1}^N\mathrm{Var}\left[X_i\right]+2\sum _{i=1}^N\sum _{j=1}^{i-1}Cov\left[X_i,X_j\right],$

which was previously obtained in ([26], p. 30).

Applying the assumption (5) of independence to the exact expression (12) of the metastable-state variance and substituting (2) into (12) yields

(13)$\mathrm{Var}\left[I\left(t\right)\right]=\sum _{i=1}^N\mathrm{Var}\left[X_i\right]=\sum _{i=1}^N(1-v_i)v_i.$

In other words, employing the independence approximation makes it possible to conveniently calculate the variance, but this method results in some deviations.

3.2. Mean Value and Variance of Steady-State Infection Fraction

In this work, we confined ourselves to the scenario demonstrated in Figure 2, where the heterogeneous SIS process spreads in the network with the topology of a complete graph. A complete graph $K_N$ of size N is defined as a graph where there exists a link between any nodal pair. Specifically, the infection rate of each node is identical, such that ${\beta }_i=\beta$ holds for any i. The nodal recovery rates are different, meaning that ${\delta }_i\ne {\delta }_j$ if $i\ne j$. We divided the node set of the whole network into two groups, represented by $G_a$ and $G_b$. The recovery rates of the nodes in groups $G_a$ and $G_b$ are denoted by ${\delta }_a$ and ${\delta }_b$, respectively. In other words, the infection rate as well as the recovery rate of any node i is either the setup $C_a=(\beta ,{\delta }_a)$ or $C_b=(\beta ,{\delta }_b)$. This scenario was previously proposed to model a flooding attack in an NDN network [15].

For the considered scenario in this work, we first prove that one node is the same as the other nodes with respect to the expectation of the steady-state nodal probability of being in the infected state. The proof is depicted as Lemma 1.

Lemma 1. 

All nodes in a complete graph are divided into two groups, $G_a$ and $G_b$, with the same infection rate but heterogeneous recovery rates. A node i is configured either with the parameters $C_a=(\beta ,{\delta }_a)$ or $C_b=(\beta ,{\delta }_b)$. The steady-state infection probability $v_i$ of node i of group $G_a$ is equal to the probability for any other node j of the same group $G_a$.

Proof of Lemma 1. 

A proof by contradiction is used, which is detailed in Appendix A. □

Regarding the steady-state infection probability $v_i$ as a function of the nodal recovery rate, we denote it by $v_i:=v_i\left({\delta }_i\right)$. Before proceeding, we introduce $N_a$ and $N_b$ to represent the numbers of nodes in $G_a$ and $G_b$, respectively.

Lemma 2. 

For the considered scenario shown in Figure 2, the steady-state infection probability $v_i$ of group $G_a$ is not equal to the probability for any other node j of the same group $G_b$. Specifically, we denote by $v_a=v_i\left({\delta }_a\right)$ and $v_b=v_i\left({\delta }_b\right)$ the probabilities for the groups $G_a$ and $G_b$, respectively. Then, we arrive at $v_a\le v_b$ if ${\delta }_a\ge {\delta }_b$.

Proof of Lemma 2. 

Using Lemma 1, one can derive this lemma. The details of the proof are shown in Appendix B. □

Lemma 3. 

We denote $y_a=v_{i\infty }\left({\delta }_a\right)$, $y_b=v_{i\infty }\left({\delta }_b\right)$, $N_a=pN$, and $N_b=(1-p)N$. For the scenario of interest as shown in Figure 2, the steady-state expectation of the fraction $y_{\infty }$ of infected nodes is

(14) $y_{\infty }=p{y}_a+(1-p)y_b,$

where $y_a$ and $y_b$ can be derived from the equations

(15) $\begin{array}{c}\hfill \beta (N_a-1)y_a^2+({\delta }_a-\beta (N_a-1)+\beta N_b{y}_b)y_a-\beta N_b{y}_b=0\end{array}$

(16) $\begin{array}{c}\hfill \beta (N_b-1)y_b^2+({\delta }_b-\beta (N_b-1)+\beta N_a{y}_a)y_b-\beta N_a{y}_a=0\end{array}$

Proof of Lemma 3. 

Using Lemma 1, one can derive this lemma. The details of the proof are shown in Appendix C. A similar brief proof was previously provided ([15], (24–26) on p. 1867). □

Furthermore, we can obtain the expression of the steady-state variance of the number of infected nodes as follows.

Proposition 1. 

The variance $\mathrm{Var}\left[I\right]$ of the metastable-state average number of infected nodes for the SIS process with heterogeneous recovery rates in a complete network for the scenario of interest, as shown in Figure 2, can be estimated using

(17) $\mathrm{Var}\left[I\right]=N_a\left(1-y_a\right)y_a+N_b\left(1-y_b\right)y_b,$

where $y_a$ and $y_b$ are the solutions of Equations (15) and (16).

Proof of Proposition 1. 

Applying the equality Lemma 1 and substituting the terms $y_a$ and $y_b$ into the approximate expression (13) yields the approximation (17). □

3.3. Bounds on Steady-State Infection

From (10), we can derive several bounds of the steady-state infection probability. Let $v_{max}={max}_j\left\{v_j\right\}$ and ${\tau }_{max}=max\{{\tau }_a,{\tau }_b\}$; then, we write (10) for the concerned scenario as

(18)$\begin{array}{ccc}\hfill v_i& \le \hfill & \hfill {\displaystyle \frac{1}{1+{\tau }_{max}{\sum }_{j=1}^N{a}_{ji}{v}_{max}}}\\ \hfill & =\hfill & \hfill {\displaystyle \frac{1}{1+(N-1){\tau }_{max}{v}_{max}}}.\end{array}$

Additionally, we can rewrite (10) as

(19)$\begin{array}{ccc}\hfill v_i& =& 1-{\displaystyle \frac{1}{1+{\tau }_i{\sum }_{j=1}^N{a}_{ji}-{\tau }_i{\sum }_{j=1}^N{a}_{ji}(1-v_j)}}\hfill \\ & =& 1-{\displaystyle \frac{1}{1+{\tau }_i{d}_i-{\tau }_i{\sum }_{j=1}^N{a}_{ji}(1-v_j)}}\hfill \\ & \le & {\displaystyle \frac{1}{1+{\tau }_i{d}_i}}.\hfill \end{array}$

Combining this with $d_i=N-1$ for the concerned scenario leads to

(20)$\begin{array}{c}\hfill v_i\le 1-{\displaystyle \frac{1}{1+(N-1){\tau }_i}}.\end{array}$

We denote $v_{min}={min}_j\left\{v_j\right\}$. Based on Lemma 1, we obtain $v_{min}=min\left\{v_a,v_b\right\}$. Substituting $v_{min}$ into (10) yields

(21)$\begin{array}{ccc}\hfill v_i& \ge & 1-{\displaystyle \frac{1}{1+{\tau }_i{\sum }_{j=1}^N{a}_{ij}{v}_{min}}}\hfill \end{array}$

(22)$\begin{array}{ccc}& =& 1-{\displaystyle \frac{1}{1+{\tau }_i(N-1)v_{min}}}.\hfill \end{array}$

Upon combining the bounds (20) and (22) of the metastable-state probability of a node being infected, we obtain the following relation:

(23)$\begin{array}{c}\hfill 1-{\displaystyle \frac{1}{1+{\tau }_i(N-1)v_{min}}}\le v_i\le 1-{\displaystyle \frac{1}{1+(N-1){\tau }_i}}.\end{array}$

3.4. Simulation Method

The simulation method employed in this work is called the $\epsilon$-SIS model. Before introducing the $\epsilon$-SIS model, we first depict the simulation of the SIS process in a fixed network. One can denote by a vector $X=\left[X_1,\cdots ,X_N\right]$ the state of the system of the SIS process in a complete graph. For a fixed N, the number $\left|\right|X\left|\right|$ of possible states is equal to $2^N$, which is large but still limited. It is possible to characterize the state change using a Markov chain. Notably, one absorbing state exists and can be represented by $X=[0,\cdots ,0]$. Once trapped in the absorbing state, the state of the system does not change to any other state. For a finite-size state space, the probability that the system traverses from any other state to the absorbing state is non-zero. Thus, the steady state of the SIS process in a finite-size network is the absorbing state.

On the other hand, the steady state of the NIMFA model actually corresponds to the metastable state of the SIS process in a finite-size network [16]. It is hard to precisely define and detect the metastable state. Van Mieghem et al. [27] came up with the so-called $\epsilon$-SIS model, which can be obtained by assigning each node an independent self-infection process such that the absorbing state of the pure SIS process is eliminated. The $\epsilon$-SIS model with a properly chosen self-infection rate provides an approximation to the pure SIS process, and it can avoid the detection of the metastable state of the SIS process. Thus, it was previously employed as the simulation method to obtain empirical results with reliable accuracy for several types of networks [28].

The self-infection process in the $\epsilon$-SIS model is an independent Poisson process assigned to each node i with a non-negative rate ${\epsilon }_i$. Even if the state of any node becomes susceptible and no infected node exists, each node may be infected by itself with a non-zero probability. In this way, the absorbing state is eliminated by merely introducing the self-infection process into the exact and pure SIS process. As a result, the steady state of the $\epsilon$-SIS model approximately corresponds to the metastable state of the SIS process. It is better to set the critical parameter $\epsilon$ as $\epsilon <1/N$ [27].

We employed the $\epsilon$-SIS model as the simulation method and implemented a discrete-time $\epsilon$-SIS simulator to produce empirical results in this work. We initiated a simulation instance by choosing a small number of the nodes as the infected nodes, with the others being susceptible. The simulation instance continued to run for a relatively long time. The changes in the system state were recorded and stored during the simulation. At the end of the simulation instance, we computed several statistical properties of the system based on the records. We set the simulation time step and the end time to ${10}^{-3}$ and 200 time units, respectively. We also set the self-infection rate to $\epsilon ={10}^{-3}$ for each node. For the scenario of interest, we set the infection and recovery rates to $\beta ,{\delta }_a=0.5,{\delta }_b=20,\epsilon ={10}^{-3}$. For a fixed infection setup ($\beta$, ${\delta }_a$, ${\delta }_b$, $\epsilon$), we only needed to run one simulation instance, record the state changes, and finally, calculate the steady-state statistics, such as the average number of infected nodes.

We show the time evolution of a single run of the simulation, e.g., the time-varying change in the number of infected nodes recorded in a simulation instance. We changed the setup of the infection rate $\beta$ for several different mixed proportions, $p=(0.2,0.5,0.8)$, of nodes with a recovery rate of ${\delta }_a$. Then, we calculated and demonstrated the metastable-state mean number of infected nodes as well as its corresponding variance.

4. Results

4.1. Time Evolution and Steady State

We developed a simulator based on the $\epsilon$-SIS model, depicted in detail in Section 3.4. The time evolution of the SIS process with heterogeneous recovery rates for one specific parameter set is demonstrated in Figure 3. At the beginning of the simulation instance, a small number of nodes were selected to be infected, while the others were susceptible. The fraction $y\left(t\right)$ of infected nodes was recorded as the simulation proceeded step by step, shown as the dotted curve. The steady-state fraction of infected nodes was estimated using (14), shown as the straight line. Additionally, the standard deviation was calculated using the root of the proposed approximation (17) of the variance.

The steady-state fraction (14) of infected nodes as a function of the infection rate for several different mixed proportions of nodes—$p=0.2$, $p=0.5$, and $p=0.8$—with a recovery rate of ${\delta }_a$ was verified via the simulations, as shown in Figure 4. The markers and curves represent the simulation and theoretical results, respectively. For any marker, the corresponding data were obtained by running a single simulation instance for the associated parameter setup and calculating the time-averaged value of the recorded number of infected nodes.

4.2. Variance of Steady-State Infection Fraction

Comparisons of the simulation results for the variance and the proposed approximation (17) of the variance for various portions $p=0.2$, $p=0.5$, and $p=0.8$ are shown in Figure 5, Figure 6 and Figure 7, respectively. The setups of these comparisons correspond to those in Figure 4. Each data point was calculated from the records via one single run of the simulation. First, we conducted the simulation, recorded the number of infected nodes, and calculated the time-averaged value of the number of infected nodes. Then, we could calculate the time-averaged value of the variance based on the time-averaged mean. The simulation results demonstrate that the variability in the number of infected nodes of the SIS process with heterogeneous infection rates is relatively small and stable for the regime of a low or very high infection rate.

5. Discussion

Figure 3 demonstrates the time evolution of the fraction of infected nodes for one instance of the simulation as well as a comparison between the simulation and the theory. For the setup of the simulation, the SIS process with heterogeneous recovery rates had a metastable state. Also, the theoretical framework, i.e., HNIMFA with heterogeneous recovery rates, can provide an approximation of the metastable-state mean value of the fraction of infected nodes of the SIS process with promising accuracy. Although the framework introduces an overestimation of $d{v}_i/dt$, the dense graph rather than the sparse graph is more resilient to the impact of the independence assumption [28].

On the other hand, it appears that the HNIMFA can provide approximations to the SIS process with heterogeneous recovery rates in finite-size complete graphs for various setups, as demonstrated in Figure 4. For the theory employing the independence assumption (4), there are observable deviations near the phase transition point, which is also referred to as the epidemic threshold.

According to the simulation results shown in Figure 5, Figure 6 and Figure 7, the variability in the number of infected nodes is tiny, while the infection rate is low and followed by a sharp increase to a maximum. Then, the variability decreases to a non-zero value. The maximum fluctuation of the variability occurs near the phase transition point, where the discrepancy between the theoretical mean infection fraction (14) and the corresponding simulation results also reaches its maximum.

The proposed variance expression (17) can provide rough approximations of the simulation results for different parameter setups with varying accuracy. For all the setups, the simulations show that the underlying trend for the regime of a high infection rate is stable. This trend is also predicted by the proposed theory. The proposed theory underestimates the simulations for $p=0.2$, as shown in Figure 5, while it overestimates the simulations for $p=0.8$, as shown in Figure 7. Notably, the approximation accuracy is sound for the regime with a relatively high infection rate for the case where $p=0.5$.

The reasons that there are discrepancies between the proposed variance expression and the simulations are two-fold. Firstly, the approximation (17) omits the co-variance part of the exact expression (12). This kind of manipulation makes it possible to estimate the variance easily, but it introduces theoretical deviations. Secondly, the network topology plays an important role in determining the range of deviations if deviations are inevitable. For a typical topology, such as a complete graph where each pair of nodes is connected, the nodal-pair impact is large such that the co-variance may become large.

6. Conclusions

The heterogeneity of epidemics’ dynamics poses challenges for their theoretical modeling and analysis and makes it hard to exactly predict the statistical properties of the topological impact on the dynamics. In this work, we proposed employing the HNIMFA framework to analyze the SIS process with heterogeneous recovery rates in finite-size networks. Specifically, we focused on the scenario where the SIS process in a finite-size complete graph shows two nodal recovery rates. We first proposed several lemmas to prove the equality of nodes with the same recovery rate. We also proposed expressions for approximating the metastable-state mean value and variance of the number of infected nodes for the SIS process with heterogeneous recovery rates. Then, we employed the $\epsilon$-SIS model as a simulation method and conducted numerous simulations to verify the proposed theories. We finally analyzed the reasons for the presence of deviations.

In our future work, we will employ methods employing high-order approximations to improve the accuracy of the approximate computations [13,29,30]. The basic NIMFA framework, as a node-based approximation, employs low-order mean-field approximation and neglects co-variance in the expression of the approximation of the variance, which leads to estimation inaccuracy for the variance. Employing high-order approximations may help to propose theories with higher accuracy. In addition to NIMFA, other theory frameworks [13,29,30] which could predict the mean value may also be capable of predicting variance more accurately. On the other hand, it would be best to maintain a balance between the order of approximation and the simplicity of the model [10,13]. Markov theory is accurate but computationally complex for relatively large finite-size networks. The NIMFA framework has a simple expression but decreases the prediction accuracy. High-order pairwise approximations or improved master equations [10,13] have the potential to predict the variance accurately while having simple expressions. Additionally, we will extend HNIMFA to accommodate the discrete-time SIS process, which shows statistical differences from a continuous-time type of model [11].

Author Contributions

Conceptualization, D.G. and L.J.; methodology, D.G.; software, D.G. and L.J.; validation, D.G., L.J. and W.F.; formal analysis, D.G.; investigation, D.G.; resources, D.G. and L.J.; data curation, D.G. and L.J.; writing—original draft preparation, D.G., L.J. and W.F.; writing—review and editing, D.G., L.J. and W.F.; visualization, D.G. and L.J.; supervision, D.G. and L.J.; project administration, D.G.; funding acquisition, D.G. and L.J. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

We would like to thank the anonymous reviewers for their insightful comments and suggestions.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

SIS	Susceptible–infected–susceptible
NIMFA	N-intertwined mean-field approximation
HNIMFA	Heterogeneous N-intertwined mean-field approximation
NDN	Named data networking

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Demonstration of Poisson processes constituting the SIS process in a pairwise network. (a) The infection process between nodes i and j, depicting the infected node i infecting its neighbor j, which is in a susceptible state, at a rate of [Forumla omitted. See PDF.]. (b) The recovery process shows an infected node i recovering at a rate of [Forumla omitted. See PDF.].

Figure 2. Schematic diagram of the scenario of interest in this work. The topology of the network of interest is a complete graph. Any node is configured with either the parameters [Forumla omitted. See PDF.] or with [Forumla omitted. See PDF.], as shown by the circles with horizontal lines and the circles with dots, respectively.

Figure 3. The fraction [Forumla omitted. See PDF.] of infected nodes versus the simulation time for a single run of the SIS process in a complete graph with [Forumla omitted. See PDF.] nodes. The blue dotted curve represents the simulation data obtained by running the [Forumla omitted. See PDF.]-SIS model. The thick red line and the dashed green lines represent the expected fraction (14) of infected nodes and the corresponding standard deviations calculated based on the variance (17), respectively.

Figure 4. The steady-state mean value of fraction [Forumla omitted. See PDF.] of infected nodes versus the infection rate. Scattered markers and curves represent the simulation results and theoretical results obtained from (14), respectively.

Figure 5. The steady-state variance [Forumla omitted. See PDF.] of the number of infected nodes for a mixture portion [Forumla omitted. See PDF.]. Markers and curves represent the simulation results and theoretical results calculated with (17), respectively.

Figure 6. The steady-state variance [Forumla omitted. See PDF.] of the number of infected nodes for mixture portion [Forumla omitted. See PDF.]. Markers and curves represent the simulation results and the theoretical results calculated by (17).

Figure 7. The steady-state variance [Forumla omitted. See PDF.] of the number of infected nodes for a mixture portion [Forumla omitted. See PDF.]. Markers and curves represent the simulation results and theoretical results calculated with (17), respectively.

Appendix A. Proof of Lemma 1

We consider a scenario where the SIS process in a complete graph has an identical nodal infection rate $\beta$ and one of two nodal recovery rates: ${\delta }_a$ or ${\delta }_b$. We can divide the network into two groups, where all the nodes within a group have the same infection rate and one of the possible recovery rates: ${\delta }_a$ or ${\delta }_b$. We assume that the first group includes the first k nodes, with the nodal number ranging from 1 to k, and the second group includes the last $N-k+1$ nodes, with the nodal number ranging from $k+1$ to N. Without loss of generality, we focus on nodes 1 and 2, which are in the same group. The steady-state infection probability (10) can be rewritten as (A1)$$\begin{array}{c}\hfill v_1=1-{\displaystyle \frac{1}{1+\left({\tau }_a{v}_2+{\tau }_a{v}_3+{\tau }_a{v}_4+\cdots +{\tau }_a{v}_k\right)+\left({\tau }_a{v}_{k+1}+\cdots +{\tau }_a{v}_N\right)}},\end{array}$$ (A2)$$\begin{array}{c}\hfill v_2=1-{\displaystyle \frac{1}{1+\left({\tau }_a{v}_1+{\tau }_a{v}_3+{\tau }_a{v}_4+\cdots +{\tau }_a{v}_k\right)+\left({\tau }_a{v}_{k+1}+\cdots +{\tau }_a{v}_N\right)}},\end{array}$$ where ${\tau }_a=\beta /{\delta }_a$. Using $C={\tau }_a{v}_3+\cdots +{\tau }_a{v}_k+\cdots +{\tau }_a{v}_N$ and substituting C into (A1) and (A2) leads to (A3)$$\begin{array}{c}\hfill v_1=1-{\displaystyle \frac{1}{1+{\tau }_a{v}_2+C}},\end{array}$$ (A4)$$\begin{array}{c}\hfill v_2=1-{\displaystyle \frac{1}{1+{\tau }_a{v}_1+C}}.\end{array}$$

A proof by contradiction is used to prove that $v_1$ is equal to $v_2$. One can assume that $v_1>v_2$. Then, upon using $c_1=v_1$ and $c_2=v_2$, it is apparent that (A5)$$\begin{array}{c}\hfill 1-{\displaystyle \frac{1}{{\tau }_a{c}_1+1+C}}<1-{\displaystyle \frac{1}{{\tau }_a{c}_2+1+C}}.\end{array}$$ Combining Equations (A3) and (A4) and inequality (A5) leads to $v_1<v_2$, which apparently contradicts the assumption $v_1>v_2$.

In summary, the probability $v_1$ is equal to $v_2$ if the nodes belong to the same group.

Appendix B. Proof of Lemma 2

Without loss of generality, we assume that the numbers of nodes of $G_a$ and $G_b$ are $N_a$ and $N_b$, respectively. Let the nodal infection probabilities of the nodes of $G_a$ and $G_b$ be $v_a$ and $v_b$, respectively; then, we recast (9) as (A6)$$\begin{array}{c}\hfill v_a=1-{\displaystyle \frac{1}{1+{\tau }_a((N_a-1)v_a+N_b{v}_b)}},\end{array}$$ for any node i of $G_a$ and (A7)$$\begin{array}{c}\hfill v_b=1-{\displaystyle \frac{1}{1+{\tau }_b(N_a{v}_a+(N_b-1)v_b)}},\end{array}$$ for any node j of $G_b$.

We take as an instance ${\delta }_a\ge {\delta }_b$ such that ${\tau }_a\le {\tau }_b$ without loss of generality before proceeding. Using these assumptions, we arrive at (A8)$$\begin{array}{ccc}\hfill v_a& \le & 1-{\displaystyle \frac{1}{1+{\tau }_b(N_a-1)v_a+{\tau }_b{N}_b{v}_b}},\hfill \end{array}$$ where the right-hand side is obtained by replacing ${\tau }_a$ with the ${\tau }_b$ of (A6).

Performing some manipulations of (A8) leads to (A9)$$\begin{array}{c}\hfill v_a\le {\displaystyle \frac{1}{1+{\tau }_b{N}_a{v}_a+{\tau }_b(N_b-1)v_b-{\tau }_b(v_a-v_b)}}.\end{array}$$

By employing a proof by contradiction and assuming $v_a\ge v_b$, we apply this assumption to the right-hand side of (A9), and then, obtain the following relation:(A10)$$\begin{array}{c}\hfill v_a\le {\displaystyle \frac{1}{1+{\tau }_b{N}_a{v}_a+{\tau }_b(N_b-1)v_b}}=v_b,\end{array}$$ which apparently contradicts the assumption $v_a\ge v_b$.

To sum up, we arrive at the inequality $v_a\le v_b$ if ${\delta }_a\ge {\delta }_b$.

Appendix C. Proof of Lemma 3

As shown in Figure 2, the nodes in a complete graph are divided into two sets with respect to the associated nodal recovery rate. Within each set, the nodes are identical, as proved in Lemma 1. For the steady-state regime, the relation that $v_{i\infty }\left({\delta }_a\right)$ is equal to $v_{j\infty }\left({\delta }_a\right)$ holds. Without loss of generality, we denote the numbers of nodes with nodal recovery rates of ${\delta }_a$ and ${\delta }_b$ by $N_a=pN$ and $N_b=(1-p)N$, respectively. The steady-state mean value of the fraction $y_{\infty }$ of infected nodes is thus expressed as (A11)$$\begin{array}{c}\hfill y_{\infty }={\displaystyle \frac{N_a{v}_i\left({\delta }_a\right)}{N}}+{\displaystyle \frac{N_b{v}_j\left({\delta }_b\right)}{N}}=p{y}_a+(1-p)y_b,\end{array}$$ where $y_a=v_{i\infty }\left({\delta }_a\right)$ for any node i with a nodal recovery rate of ${\delta }_a$.

Substituting the notation $N_a$ and $N_b$ into (9) and applying Lemma 1 yields (A12)$$\begin{array}{c}\hfill {\delta }_a{y}_a=\beta (N_a-1)y_a(1-y_a)+\beta N_b{y}_b(1-y_a),\end{array}$$ for any node i with a nodal recovery rate of ${\delta }_i={\delta }_a$. After some manipulations, we rewrite (A12) as (A13)$$\begin{array}{c}\hfill \beta (N_a-1)y_a^2+\left({\delta }_a-\beta (N_a-1)+\beta N_b{y}_b\right)y_a-\beta N_b{y}_b=0.\end{array}$$ In a similar manner, we can obtain the equation for any node j with a nodal recovery rate of ${\delta }_j={\delta }_b$ as (A14)$$\begin{array}{c}\hfill \beta (N_b-1)y_b^2+\left({\delta }_b-\beta (N_b-1)+\beta N_a{y}_a\right)y_b-\beta N_a{y}_a=0.\end{array}$$

A similar proof can be found in ([15], (24–26) on p. 1867), which is presented briefly without extensive detail.

References

1. Pastor-Satorras, R.; Vespignani, A. Epidemic dynamics and endemic states in complex networks. Phys. Rev. E; 2001; 63, 066117. [DOI: https://dx.doi.org/10.1103/PhysRevE.63.066117] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11415183]

2. Pastor-Satorras, R.; Castellano, C.; Van Mieghem, P.; Vespignani, A. Epidemic processes in complex networks. Rev. Mod. Phys.; 2015; 87, pp. 925-979. [DOI: https://dx.doi.org/10.1103/RevModPhys.87.925]

3. Bianconi, G. Multilayer Networks: Structure and Function; Oxford University Press: Oxford, UK, 2018; [DOI: https://dx.doi.org/10.1093/oso/9780198753919.001.0001]

4. de Arruda, G.F.; Rodrigues, F.A.; Moreno, Y. Fundamentals of spreading processes in single and multilayer complex networks. Phys. Rep.; 2018; 756, pp. 1-59. [DOI: https://dx.doi.org/10.1016/j.physrep.2018.06.007]

5. Majhi, S.; Perc, M.; Ghosh, D. Dynamics on higher-order networks: A review. J. R. Soc. Interface; 2022; 19, 20220043. [DOI: https://dx.doi.org/10.1098/rsif.2022.0043]

6. de Arruda, G.F.; Aleta, A.; Moreno, Y. Contagion dynamics on higher-order networks. Nat. Rev. Phys.; 2024; 6, pp. 468-482. [DOI: https://dx.doi.org/10.1038/s42254-024-00733-0]

7. Tomovski, I.; Basnarkov, L.; Abazi, A. Discrete-Time Non-Markovian SEIS Model on Complex Networks. IEEE Trans. Netw. Sci. Eng.; 2022; 9, pp. 552-563. [DOI: https://dx.doi.org/10.1109/TNSE.2021.3125191]

8. Basnarkov, L.; Tomovski, I.; Sandev, T.; Kocarev, L. Non-Markovian SIR epidemic spreading model of COVID-19. Chaos Solitons Fractals; 2022; 160, 112286. [DOI: https://dx.doi.org/10.1016/j.chaos.2022.112286]

9. Han, L.; Lin, Z.; Yin, Q.; Tang, M.; Guan, S.; Boguñá, M. Non-Markovian epidemic spreading on temporal networks. Chaos Solitons Fractals; 2023; 173, 113664. [DOI: https://dx.doi.org/10.1016/j.chaos.2023.113664]

10. Ward, J.A.; Tapper, A.; Simon, P.L.; Mann, R.P. Micro-scale foundation with error quantification for the approximation of dynamics on networks. Commun. Phys.; 2022; 5, 71. [DOI: https://dx.doi.org/10.1038/s42005-022-00834-1]

11. Silva, D.H.; Rodrigues, F.A.; Ferreira, S.C. Accuracy of discrete- and continuous-time mean-field theories for epidemic processes on complex networks. Phys. Rev. E; 2024; 110, 014302. [DOI: https://dx.doi.org/10.1103/PhysRevE.110.014302]

12. Newman, M.E.J. Message passing methods on complex networks. Proc. R. Soc. A Math. Phys. Eng. Sci.; 2023; 479, 20220774. [DOI: https://dx.doi.org/10.1098/rspa.2022.0774]

13. Ortega, E.; Machado, D.; Lage-Castellanos, A. Dynamics of epidemics from cavity master equations: Susceptible-infectious-susceptible models. Phys. Rev. E; 2022; 105, 024308. [DOI: https://dx.doi.org/10.1103/PhysRevE.105.024308] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35291082]

14. Cator, E.; Van Mieghem, P. Susceptible-infected-susceptible epidemics on the complete graph and the star graph: Exact analysis. Phys. Rev. E; 2013; 87, 012811. [DOI: https://dx.doi.org/10.1103/PhysRevE.87.012811] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23410392]

15. Wang, K.; Guo, D.; Quan, W. Analyzing NDN NACK on Interest Flooding Attack via SIS Epidemic Model. IEEE Syst. J.; 2020; 14, pp. 1862-1873. [DOI: https://dx.doi.org/10.1109/JSYST.2019.2923841]

16. Van Mieghem, P.; Omic, J.; Kooij, R. Virus Spread in Networks. IEEE/ACM Trans. Netw.; 2009; 17, pp. 1-14. [DOI: https://dx.doi.org/10.1109/TNET.2008.925623]

17. Chakrabarti, D.; Wang, Y.; Wang, C.; Leskovec, J.; Faloutsos, C. Epidemic thresholds in real networks. ACM Trans. Inf. Syst. Secur.; 2008; 10, pp. 1:1-1:26. [DOI: https://dx.doi.org/10.1145/1284680.1284681]

18. Omic, J.; Van Mieghem, P. Epidemic Spreading in Networks–Variance of the Number of Infected Nodes; Technical Report report20090707 Quantum and Computer Engineering Department, Delft University of Technology: Delft, The Netherlands, 2009.

19. Van Mieghem, P.; Omic, J. In-homogeneous Virus Spread in Networks. arXiv; 2014; arXiv: 1306.2588

20. Jiao, L.; Yin, H.; Guo, D.; Lyu, Y. Heterogeneous Malware Spread Process in Star Network. Proceedings of the 2017 IEEE 37th International Conference on Distributed Computing Systems Workshops (ICDCSW); Atlanta, GA, USA, 5–8 June 2017; pp. 265-269. [DOI: https://dx.doi.org/10.1109/ICDCSW.2017.35]

21. Cui, S.; Liu, F.; Jardón-Kojakhmetov, H.; Cao, M. Discrete-time layered-network epidemics model with time-varying transition rates and multiple resources. Automatica; 2024; 159, 111303. [DOI: https://dx.doi.org/10.1016/j.automatica.2023.111303]

22. Pagliara, R.; Leonard, N.E. Adaptive Susceptibility and Heterogeneity in Contagion Models on Networks. IEEE Trans. Autom. Control; 2021; 66, pp. 581-594. [DOI: https://dx.doi.org/10.1109/TAC.2020.2985300]

23. Yuan, X.; Yao, Y.; Li, X.; Feng, M. Impact of time-dependent infection rate and self-isolation awareness on networked epidemic propagation. Nonlinear Dyn.; 2024; 112, pp. 15653-15669. [DOI: https://dx.doi.org/10.1007/s11071-024-09832-0]

24. Jing, X.; Liu, G.; Jin, Z. Stochastic dynamics of an SIS epidemic on networks. J. Math. Biol.; 2022; 84, 50. [DOI: https://dx.doi.org/10.1007/s00285-022-01754-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35513730]

25. Guo, D.; Jiao, L.; Jiao, J.; Meng, K. Variance of the Infection Number of Heterogeneous Malware Spread in Network. Appl. Sci.; 2024; 14, 3972. [DOI: https://dx.doi.org/10.3390/app14103972]

26. Van Mieghem, P. Performance Analysis of Communications Networks and Systems; Cambridge University Press: Cambridge, MA, USA, 2006; [DOI: https://dx.doi.org/10.1017/CBO9780511616488]

27. Van Mieghem, P.; Cator, E. Epidemics in networks with nodal self-infection and the epidemic threshold. Phys. Rev. E; 2012; 86, 016116. [DOI: https://dx.doi.org/10.1103/PhysRevE.86.016116] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23005500]

28. Li, C.; van de Bovenkamp, R.; Van Mieghem, P. Susceptible-infected-susceptible model: A comparison of N-intertwined and heterogeneous mean-field approximations. Phys. Rev. E; 2012; 86, 026116. [DOI: https://dx.doi.org/10.1103/PhysRevE.86.026116]

29. Cator, E.; Van Mieghem, P. Second-order mean-field susceptible-infected-susceptible epidemic threshold. Phys. Rev. E; 2012; 85, 056111. [DOI: https://dx.doi.org/10.1103/PhysRevE.85.056111]

30. Wang, K.; Liu, X.F.; Guo, D. High-Order Mean-Field Approximations for Adaptive Susceptible-Infected-Susceptible Model in Finite-Size Networks. Complexity; 2021; 2021, 6637761. [DOI: https://dx.doi.org/10.1155/2021/6637761]