In-homogeneous Virus Spread in Networks

TL;DR

This paper extends the NIMFA virus spread model to heterogeneous networks, characterizing steady states and thresholds using generalized Laplacians and eigenvalues, and analyzing the effects of curing rates on infection probabilities.

Contribution

It introduces a full heterogeneous extension of NIMFA, defines a generalized Laplacian for metastable states, and characterizes the epidemic threshold via eigenvalues of a modified adjacency matrix.

Findings

Steady-state infection probabilities are convex in own curing rates.

Critical threshold is determined by the largest eigenvalue of a modified adjacency matrix.

The model applies to any network with N nodes, generalizing previous homogeneous models.

Abstract

Our $N$-intertwined model (now called NIMFA) for virus spread in any network with $N$ nodes is extended to a full heterogeneous setting. The metastable steady-state nodal infection probabilities are specified in terms of a generalized Laplacian, that possesses analogous properties as the classical Laplacian in graph theory. The critical threshold that separates global network infection from global network health is characterized via an $N$ dimensional vector that makes the largest eigenvalue of a modified adjacency matrix equal to unity. Finally, the steady-state infection probability of node $i$ is convex in the own curing rate $\delta_{i}$, but concave in the curing rates $\delta_{j}$ of the other nodes $1\leq j\neq i\leq N$ in the network.

Equations

$$Q_{i}\left(t\right)=\left[\begin{array}[c]{cc}-q_{1;i}&q_{1;i}\\ q_{2;i}&-q_{2;i}\end{array}\right]$$

$$q_{1;i}=\sum _{j=1}^N{\beta }_j{a}_{ij}{1}_{\left\{X_j(t)=1\right\}}$$

$$\mathrm{Pr}[X_i(t+\Delta t)=1|X_i\left(t\right)=0]=q_{1;i}\Delta t+o(\Delta t)$$

$$E\left[q_{1;i}\right]=E\left[\sum _{j=1}^N{\beta }_j{a}_{ij}{1}_{\left\{X_j(t)=1\right\}}\right]$$

$$E\left[q_{1;i}\right]=\sum _{j=1}^N{\beta }_j{a}_{ij}\mathrm{Pr}[X_j(t)=1]$$

$$\overline{Q_{i}(t)}=\left[\begin{array}[c]{cc}-E\left[q_{1;i}\right]&E\left[q_{1;i}\right]\\ \delta_{i}&-\delta_{i}\end{array}\right]$$

$$\frac{d{v}_i\left(t\right)}{dt}=\sum _{j=1}^N{\beta }_j{a}_{ij}{v}_j\left(t\right)-v_i\left(t\right)\left(\sum _{j=1}^N{\beta }_j{a}_{ij}{v}_j\left(t\right)+{\delta }_i\right)$$

$$\left\{\begin{array}[c]{c}\frac{dv_{1}\left(t\right)}{dt}=\sum_{j=1}^{N}\beta_{j}a_{1j}v_{j}\!\left(t\right)\!-v_{1}\!\left(t\right)\!\left(\sum_{j=1}^{N}\beta_{j}a_{1j}v_{j}\!\left(t\right)+\delta_{1}\right)\\ \frac{dv_{2}\left(t\right)}{dt}=\sum_{j=1}^{N}\beta_{j}a_{2j}v_{j}\!\left(t\right)\!-v_{2}\!\left(t\right)\!\left(\sum_{j=1}^{N}\beta_{j}a_{2j}v_{j}\!\left(t\right)+\delta_{2}\right)\\ \vdots\\ \frac{dv_{N}\left(t\right)}{dt}=\sum_{j=1}^{N}\beta_{j}a_{Nj}v_{j}\!\left(t\right)\!-v_{N}\!\left(t\right)\!\left(\sum_{j=1}^{N}\beta_{j}a_{Nj}v_{j}\!\left(t\right)+\delta_{N}\right)\end{array}\right.$$

$$V\left(t\right)=\left[\begin{array}[c]{cccc}v_{1}\left(t\right)&v_{2}\left(t\right)&\cdots&v_{N}\left(t\right)\end{array}\right]^{T}$$

$$\frac{dV\left(t\right)}{dt}=A\text{diag}\left({\beta }_j\right)V\left(t\right)-\text{diag}\left(v_i\left(t\right)\right)\left(A\text{diag}\left({\beta }_j\right)V\left(t\right)+C\right)$$

$$A\text{diag}\left({\beta }_i\right)V_{\infty }-\text{diag}\left(v_{i\infty }\right)\left(A\text{diag}\left({\beta }_i\right)V_{\infty }+C\right)=0$$

$$w=A\text{diag}\left({\beta }_i\right)V_{\infty }+C$$

$$w-C=\text{diag}\left(v_{i\infty }\right)w$$

$$\left(I-\text{diag}\left(v_{i\infty }\right)\right)w=C$$

$$w=\text{diag}\left(\frac{1}{1-v_{i\infty }}\right)C$$

$$A\text{diag}\left({\beta }_i\right)V_{\infty }$$

$$=\text{diag}\left(\frac{{\delta }_i}{1-v_{i\infty }}\right)V_{\infty }$$

$$\sum _{j=1}^N{a}_{ij}{\beta }_j{v}_{j\infty }=\frac{v_{i\infty }{\delta }_i}{1-v_{i\infty }}$$

$$\mathcal{Q}\left(\frac{1}{{\tau }_i\left(1-v_{i\infty }\right)}\right){\tilde{V}}_{\infty }=0$$

$$\mathcal{Q}\left(q_i\right)$$

$$=\text{diag}\left(q_i-d_i\right)+Q$$

$$\sum _{j=1}^N{\beta }_j{v}_{j\infty }{y}_j=0$$

$$0={\lambda }_N\left(\mathcal{Q}\right)\le {\lambda }_{N-1}\left(\mathcal{Q}\right)\le \dots \le {\lambda }_1\left(\mathcal{Q}\right)$$

$$q_k-d_k\le \lambda \le d_k+q_k$$

$$0<v_{j\infty }\le 1-\frac{1}{{\tau }_j{d}_j}$$

$$\sum _{k=1}^{N-1}{\lambda }_k\left(\mathcal{Q}\right)=\sum _{i=1}^N\frac{1}{{\tau }_i\left(1-v_{i\infty }\right)}$$

$$\text{trace}\left({\mathcal{Q}}^2\left(q_i\right)\right)$$

$$=\sum _{i=1}^N\frac{1}{{\tau }_i^2{\left(1-v_{i\infty }\right)}^2}+2L$$

$$\sum _{k=1}^{N-1}{\lambda }_k^2\left(\mathcal{Q}\right)=\sum _{i=1}^N\frac{1}{{\tau }_i^2{\left(1-v_{i\infty }\right)}^2}+2L$$

$$u^{T}A\text{diag}\left({\beta }_i\right)V_{\infty }=u^T\text{diag}\left(\frac{{\delta }_i}{1-v_{i\infty }}\right)V_{\infty }$$

$$\left(u^T\text{diag}\left(\frac{{\delta }_i}{1-v_{i\infty }}\right)-D^T\text{diag}\left({\beta }_i\right)\right)V_{\infty }=0$$

$$\mathcal{Q}\left(q_i\right){\tilde{V}}_{\infty }=\text{diag}\left(q_i-d_i\right){\tilde{V}}_{\infty }+Q{\tilde{V}}_{\infty }$$

Full text

I Introduction

This paper generalizes our $N$-Intertwined Mean-Field Approximation (NIMFA) for virus spread in networks, presented in [12] and [11, Chapter 17], to a heterogeneous setting. Heterogeneity rather than homogeneity abounds in real networks. For example, in data communications networks, the transmission capacity, age, performance, installed software, security level and other properties of networked computers are generally different. Social and biological networks are very diverse: a population often consists of a mix of weak and strong, or old and young species or of completely different types of species. The network topology for transport by airplane, car, train, ship is different. Many more examples can be added illustrating that homogeneous networks are the exception rather than the rule. This diversity in the “nodes” and “links” of real networks will thus likely affect the spreading pattern of viruses, that are here understood as malicious challenges of a network.

NIMFA approximates the continuous-time Markov susceptible-infected-susceptible (SIS) epidemic process on a network with $N$ nodes, that was earlier considered by Ganesh *et al. *[3] and by Wang et al. [14] in discrete-time. Each node in the network is either infected or healthy. In a heterogeneous setting, an infected node $i$ can infect its neighbors with an infection rate $\beta_{i}$, but it is cured with curing rate $\delta_{i}$. Once cured and healthy, the node is again prone to the virus. Both infection and curing processes are independent.

Previously in [12], only a homogeneous virus spread was investigated, where all infection rates $\beta_{i}=\beta$ and all curing rates $\delta_{i}=\delta$ were the same for each node. We believe that the extension to a full heterogeneous setting is, perhaps, the best SIS model that we can achieve. The exact Markovian model, described and analyzed in [12], has $2^{N}$ states, which makes it infeasible to compute for realistic sizes of networks. Moreover, the exact Markovian model possesses as steady-state the overall healthy state, which is an absorbing state, that is, unfortunately, only reached after an extreme and unrealistically long time. The heterogeneous NIMFA makes one approximation, a mean field approximation as shown in Section II and in [12], that results in a set of $N$ non-linear equations. Hence, NIMFA trades computational feasibility, a reduction of $2^{N}$ linear equations to $N$ non-linear ones, at the expense of exactness. The last point, the accuracy of NIMFA is shown in [12] (and further in [6]) to be overall remarkably good, with a worst case performance near the critical threshold, which is a realistic and observable artifact of the metastable steady-state that does not exist in the exact Markovian steady-state. Below the critical epidemic threshold, infection vanishes exponentially fast in time and above the critical threshold the network stays infected to a degree determined by the effective infection vector $\tau$, with components $\tau_{i}=\frac{\beta_{i}}{\delta_{i}}$.

A major new insight is that the metastable steady-state can be written in terms of a generalized Laplacian matrix that bears similar deep properties as the Laplacian matrix of a graph (see e.g. [1], [2] and [10]). In a heterogeneous setting, the critical threshold is characterized by an effective infection vector, instead of one scalar in the homogeneous case equal to $\tau_{\text{hom;}c}=\frac{1}{\lambda_{\max}\left(A\right)}$, where $\lambda_{\max}\left(A\right)$ is the largest eigenvalue of the adjacency matrix $A$ of the graph. This critical vector determines a critical surface in the $N$-dimensional space spanned by the vector components $\tau_{1},\ldots,\tau_{N}$. We also prove that the steady-state infection probability $v_{i\infty}$ of node $i$ is convex in the curing rate $\delta_{i}$, given all other curing rates $\delta_{j}$ are the same.** **This convexity result is applied in a virus protection game played by the individual and selfish nodes in a network [7].

II $N$-intertwined continuous Markov chains with $2$ states

This section extends the homogeneous NIMFA in [12] to a heterogeneous setting. Although analogous to the corresponding section in [12], its inclusion makes this paper self-contained.

By separately observing each node, we will model the virus spread in a bi-directional network specified by a symmetric adjacency matrix $A$. Every node $i$ at time $t$ in the network has two states: infected with probability $\Pr[X_{i}\left(t\right)=1]$ and healthy with probability $\Pr[X_{i}\left(t\right)=0]$. At each moment $t$, a node can only be in one of two states, thus $\Pr[X_{i}(t)=1]+\Pr[X_{i}(t)=0]=1$. If we apply Markov theory, the infinitesimal generator $Q_{i}\left(t\right)$ of this two-state continuous Markov chain is,

[TABLE]

with $q_{2;i}=\delta_{i}$ and

[TABLE]

where the indicator function $1_{x}=1$ if the event $x$ is true else it is zero. The coupling of node $i$ to the rest of the network is described by an infection rate $q_{1;i}$ that is a random variable, which essentially makes the process doubly stochastic. This observation is crucial. For, using the definition of the infinitesimal generator [8, p. 181],

[TABLE]

the continuity and differentiability shows that this process is not Markovian anymore. The random nature of $q_{1;i}$ is removed by an additional conditioning to all possible combinations of rates, which is equivalent to conditioning to all possible combinations of the states $X_{j}(t)=1$ (and their complements $X_{j}(t)=0$) of the neighbors of node $i$. Hence, the number of basic states dramatically increases. Eventually, after conditioning each node in such a way, we end up with a $2^{N}$– state Markov chain, studied in [12].

Instead of conditioning, we replace the actual, random infection rate by an effective or average infection rate, which is basically a mean field approximation,

[TABLE]

In general, we may take the expectation over the rates $\beta_{i}$, the network topology via the matrix $A$ and the states $X_{j}(t)$. Since we assume that both the infection rates $\beta_{i}$ and the network are constant and given, we only average over the states. Using $E\left[1_{x}\right]=\Pr\left[x\right]$ (see e.g. [8]), we replace $q_{1;i}$ by

[TABLE]

which results in an effective infinitesimal generator,

[TABLE]

The effective $\overline{Q_{i}(t)}$ allows us to proceed with Markov theory. Denoting $v_{i}\left(t\right)=\Pr[X_{i}(t)=1]$ and recalling that $\Pr[X_{i}(t)=0]=1-v_{i}\left(t\right)$, the Markov differential equation [11, (10.11) on p. 208] for state $X_{i}(t)=1$ turns out to be non-linear

[TABLE]

Each node obeys a differential equation as (2),

[TABLE]

Written in matrix form, with

[TABLE]

we arrive at

[TABLE]

where diag$\left(v_{i}\left(t\right)\right)$ is the diagonal matrix with elements $v_{1}\left(t\right),v_{2}\left(t\right),\ldots\,,v_{N}\left(t\right)$ and the curing rate vector is $C=\left(\delta_{1},\delta_{2},\ldots,\delta_{N}\right)$.

We note that $A$diag$\left(\beta_{i}\right)$ is, in general and opposed to the homogeneous setting, not symmetric anymore, unless $A$ and diag$\left(\beta_{i}\right)$ commute, in which case the eigenvalue $\lambda_{i}\left(A\text{diag}\left(\beta_{i}\right)\right)=\lambda_{i}\left(A\right)\beta_{i}$ and both $\beta_{i}$ and $\lambda_{i}\left(A\right)$ have a same eigenvector $x_{i}$.

III General in-homogenous steady-state

III-A The steady-state equation

The metastable steady-state follows from (3) as

[TABLE]

where $V_{\infty}=\lim_{t\rightarrow\infty}V\left(t\right)$. We define the vector

[TABLE]

and write the stead-state equation as

[TABLE]

or

[TABLE]

Ignoring extreme virus spread conditions (the absence of curing ($\delta_{i}=0$) and an infinitely strong infection rate $\beta_{i}\rightarrow\infty$), then the infection probabilities $v_{i\infty}$ cannot be one such that the matrix $\left(I-\text{diag}\left(v_{i\infty}\right)\right)=$ diag$\left(1-v_{i\infty}\right)$ is invertible. Hence,

[TABLE]

Invoking the definition (4) of $w$, we obtain

[TABLE]

The $i$-th row of (5) yields the nodal steady state equation,

[TABLE]

Let $\widetilde{V}_{\infty}=$ diag$\left(\beta_{i}\right)V_{\infty}$ and the effective spreading rate for node $i$, $\tau_{i}=\frac{\beta_{i}}{\delta_{i}}$, then we arrive at

[TABLE]

where the symmetric matrix

[TABLE]

can be interpreted as a generalized Laplacian111All eigenvalues of the Laplacian $Q=\Delta-A$ in a connected graph are positive, except for the smallest one that is zero. Hence, $Q$ is positive semi-definite. Much more properties of the Laplacian $Q$ are found e.g. in [1] and [2]., because $\mathcal{Q}\left(d_{i}\right)=Q=\Delta-A$, where $\Delta=$ diag$\left(d_{i}\right)$. The observation that the non-linear set of steady-state equations can be written in terms of the generalized Laplacian $\mathcal{Q}\left(q_{i}\right)$ is fortunate, because, as will be shown in Section III-B, the powerful theory of the “normal” Laplacian $Q$ applies.

The modified steady-state vector $\widetilde{V}_{\infty}$ is orthogonal to each row (or, by symmetry, each column) vector of $\mathcal{Q}\left(\frac{1}{\tau_{i}\left(1-v_{i\infty}\right)}\right)$. A non-zero modified steady-state vector $\widetilde{V}_{\infty}$ is thus only possible provided $\det\mathcal{Q}\left(\frac{1}{\tau_{i}\left(1-v_{i\infty}\right)}\right)=0$. In other words, the generalized Laplacian $\mathcal{Q}\left(\frac{1}{\tau_{i}\left(1-v_{i\infty}\right)}\right)$ should have a zero eigenvalue with the modified steady-state vector $\widetilde{V}_{\infty}$ as corresponding eigenvector. Since the vectors $B=\left(\beta_{1},\beta_{2},\ldots,\beta_{N}\right)$ and $C=\left(\delta_{1},\delta_{2},\ldots,\delta_{N}\right)$ are given, the non-linear eigenvector problem (7) has, in general, a solution that cannot simply be recast to the homogeneous case where $B=\beta u$ and $C=\delta u$ (or $\beta_{i}=\beta$ and $\delta_{i}=\delta$ for all $1\leq i\leq N$) in which the all-one vector $u=\left(1,1,\ldots,1\right)$.

III-B The generalized Laplacian $\mathcal{Q}\left(q_{i}\right)$

Since $\mathcal{Q}\left(q_{i}\right)$ is symmetric, all eigenvectors are orthogonal such that, with $\widetilde{V}_{\infty}=$ diag$\left(\beta_{i}\right)V_{\infty}$,

[TABLE]

where $y$ is the eigenvector belonging to eigenvalue $\lambda\left(\mathcal{Q}\left(q_{i}\right)\right)\neq 0$.

Theorem 1

If the network $G$ is connected, all eigenvalues of $\mathcal{Q}\left(q_{i}\right)$ are positive, except for the smallest one $\lambda_{N}\left(\mathcal{Q}\right)=0$.

Proof: The theorem is a consequence of the Perron-Frobenius Theorem (see e.g. [4]) for a non-negative, irreducible matrix. Indeed, consider the non-negative matrix $q_{\max}I-\mathcal{Q}\left(q_{i}\right)$, where $q_{\max}=\max_{1\leq i\leq N}q_{i}$, whose eigenvalues are $\xi_{k}=q_{\max}-\lambda_{k}\left(\mathcal{Q}\right)$ for $1\leq k\leq N$. Since $G$ is connected, then $q_{\max}I-\mathcal{Q}\left(q_{i}\right)$ is irreducible and the Perron-Frobenius Theorem states that the largest eigenvalue $r=\max_{1\leq k\leq N}\xi_{k}$ of $q_{\max}I-\mathcal{Q}\left(q_{i}\right)$ is positive and simple and the corresponding eigenvector $x_{r}$ has positive components. Hence, $\mathcal{Q}\left(q_{i}\right)x_{r}=\left(q_{\max}-r\right)x_{r}$. Since eigenvectors of a symmetric matrix are orthogonal while $\widetilde{V}_{\infty}^{T}x_{r}>0$, $x_{r}$ must be proportional to $\widetilde{V}_{\infty}$, and thus $q_{\max}=r$. Since there is only one such eigenvector $x_{r}$ and since the eigenvalue $r>q_{\max}-\lambda_{k}\left(\mathcal{Q}\right)$ for all $k$ (except that $k$ for which $\lambda_{k}\left(\mathcal{Q}\right)=0$, which is thus the smallest eigenvalue), all other eigenvalues of $\mathcal{Q}\left(q_{i}\right)$ must exceed zero. $\square\vskip 6.0pt plus 2.0pt minus 2.0pt$

If the graph $G$ is disconnected which means that $A$ is reducible [8], the Theorem 1 still applies (see e.g. [4]), however, under the slightly weakened form that $x_{r}$ has non-negative components (instead of positive, hence, zero components can occur) and that the largest eigenvalue $r$ is non-zero (not necessarily strict positive). The consequence is that more than one zero eigenvalue can occur. From the point of virus spread, we may ignore disconnected graphs, because the theory can be applied to each connected component (cluster) of the network $G$. The symmetry of $\mathcal{Q}\left(q_{i}\right)$ implies that all eigenvalues are real and can be ordered. By Theorem 1, we have

[TABLE]

Gerschgorin’s theorem [15, p. 71-75] indicates that the eigenvalues of $\mathcal{Q}\left(q_{i}\right)$ are centered around $q_{i}$ with radius equal to the degree $d_{i}$, i.e. an eigenvalue $\lambda$ of $\mathcal{Q}\left(q_{i}\right)$ lies in an interval $\left|\lambda-q_{k}\right|\leq d_{k}$ for some $1\leq k\leq N$. Thus, there is an eigenvalue $\lambda$ of $\mathcal{Q}\left(q_{i}\right)$ that obeys

[TABLE]

A solution of (7) requires that at least one eigenvalue of $\mathcal{Q}\left(q_{i}\right)$ is zero, while Theorem 1 states that there is only one zero eigenvalue. Hence, precisely one, say the $j$-th, of the Gerschgorin line segments that contain the eigenvalue $\lambda_{N}\left(\mathcal{Q}\right)=0$, must obey $q_{j}\leq d_{j}$ to have a non-zero solution of (7). However, more Gerschgorin segments may obey $q_{k}-d_{k}\leq 0$. This couples $\frac{1}{\tau_{j}\left(1-v_{j\infty}\right)}\leq d_{j}$ for at least one $j$ component and shows that, when $v_{j\infty}\rightarrow 1$, there must hold that $\tau_{j}\rightarrow\infty$. Hence, for at least one component $j$, there holds that

[TABLE]

where the lower bound follows, by the Perron-Frobenius Theorem, from the fact that the network