An SIS epidemic model with time delay and stochastic perturbation on heterogeneous networks

1.   Introduction

Transmission dynamics studies problems arising in the real world, for example, spread of diseases in population, virus propagation in computer networks and diffusion of information. It strongly depends on properties of the contact network. Scale-free network by Barabási and Albert ^[1] can well depict complex connectivity patterns in nature and human society such as, e.g., social network, computer network and World Wide Web. Therefore, compared with classical epidemic models on homogeneous networks, it is more significant to study spreading dynamics on heterogeneous networks, i.e., scale-free networks.

The spreading dynamics on complex networks has attracted increasing attention. A lot of epidemiological models on complex networks, including SI ^[2,3], SIS ^[4,5], SIR ^[6], SIQS ^[7], and so on ^[8,9], have been established successively. The pioneering work of Pastor-Satorras and Vespignani ^[4,5] introduced the SIS model on scale-free networks by mean-field approximation. They showed that the epidemic threshold is infinitesimal with network size increasing, which makes the spread of infections tremendously strengthened. Wang and Dai ^[10] analyzed the network-based SIS model theoretically for the first time. They derived that the epidemic threshold is the critical parameter for global stability of the disease-free and endemic equilibria. In the meantime, d'Onofrio ^[11] conducted further research on this aspect.

In recent years, as for modification and extension of the network-based SIS model, two directions have always been of great concern. One is introducing time delays to simulate incubation periods, infection periods, immunity periods and so on. Models expressed by functional differential equations (FDE) have been formulated. Liu et al. ^[12,13] proposed and discussed SIS epidemic models on scale-free networks with discrete and distributed delays, respectively, in which the time delays represent infection periods. Kang et al. ^[14] established an SIS model with time delay denoting the incubation period of disease in a vector's body and analyzed the global stability of equilibria. Furthermore, Kang and Fu ^[15] considered two transmitting ways (by human and by vector) and discussed another new delayed SIS model on heterogeneous networks. Another way is to enter interactions of uncertain environments in the models. Stochastic differential equation (SDE) SIS models on complex networks have been developed. Bonaccorsi et al. ^[16] proposed an SIS model with stochastic infection rates on networks and proved the conditions for extinction and stochastic persistence of epidemics. Krause et al. ^[17] analyzed dynamical behaviors of a stochastic SIS epidemic model in metapopulation by numerical simulation. Some control approaches were designed to control epidemic spreading. Yang and Jin ^[18] introduced a stochastic SIS model driven by Lévy noise on networks. The stability of the disease-free equilibrium and the sufficient condition for persistence were proved.

However, to the best of our knowledge, there has been little research about the interplay of noise and delay on spreading dynamics in complex networks. We will present an epidemic model with time delay and stochastic perturbation. It is based on the following general version of the delayed SIS model with an infective vector ^[14],

where $\lambda(k)$ is the degree-dependent infection rate which is bounded ^[19], such as $\lambda k$ ^[14], $\lambda c(k)$ ^[20], and so on. $S_k(t)$ and $I_k(t)$, $k = 1, \cdots, n$, represent the relative densities of healthy and infected nodes with degree $k$ at time $t$, respectively. Here, $\mu$ is the recovery rate, and $\gamma$ is the mortality rate which is equal to the birth rate $\Lambda$, i.e., the network size is time invariant. $d_m$ is the nature death of the vector and $\tau$ is the incubation. $e^{-d_m\tau}$ represents the probability of infected vectors who were infected at time $t-\tau$ but did not die during the time period $\tau$. The parameters aforesaid are all positive. The vector's density is simply proportional to $\Theta(t-\tau)$ expressed as

where $P(j)$ is the degree distribution, $\langle k\rangle = \sum_{k = 1}^nkP(k)$ is the average degree of the network, $\varphi (j)$ denotes an infected node, with degree $j$, occupied edges which can transmit the disease ^[21], and $n$ is the maximum degree of nodes in this network. We also define $\langle f(k)\rangle = \sum_{k = 1}^nf(k)P(k)$.

By analogy with the results of Ref. ^[14], the following equivalent system of (1.1)

always has a disease-free equilibrium $I^0 = (0, 0, \cdots, 0)$. If $R_0 < 1$, $I^0$ is globally asymptotically stable. Whereas, $I^0$ is unstable and there exists a globally stable endemic equilibrium $I^* = (I_1^*, I_2^*, \cdots, I_n^*)$ when $R_0 > 1$. And here,

Now, we introduce interactions of random environments to system (1.1) by replacing the infection rate $\lambda (k)$ with

where $\sigma_k > 0, k = 1, \cdots, n$ represent the noise intensities, and $dB_k(t)\; (k = 1, \cdots, n)$ is an n-dimensional standard white noise, i.e., $B_k(t)\; (k = 1, \cdots, n)$ is an n-dimensional standard Brownian motion with $B_k(0) = 0$. It is essential to assume that the diffusion coefficient depends on relative densities $I_k$. The infection rate varies around a mean value, and the variance gets smaller with the relative densities of infected nodes decreasing. It guarantees that the solution has physical meaning, namely, it remains above zero. Moreover, the dependence of noise intensities on the solution is typical in previous articles concerning population dynamics ^[22] and spread dynamics ^[16] with environmental noise. Then the stochastic system we study takes the following form

Let $\mathcal{D} = (0, 1)^n$ be the $n$-th Cartesian power of the interval $(0, 1)$. And denote by $C = C([-\tau, 0], \mathcal{D})$ the Banach space of continuous functions mapping the interval $[-\tau, 0]$ into $\mathcal{D}$ with norm

where $|\phi_k(\theta)|_\tau = \sup_{-\tau\leq\theta\leq 0}|\phi_k(\theta)|$. For a practical consideration, the initial condition of system (1.3) can be considered as

where $\phi = (\phi_1, \phi_2, \cdots, \phi_n)\in C$.

In this paper, we focus on the dynamical behaviors of system (1.3). The existence, uniqueness and boundedness of the solution to system (1.3) are discussed in Section 2. In Section 3, we investigate the dynamics of system (1.3) and present sufficient conditions for the exponential extinction and permanence of the disease, respectively. Numerical simulations are given to demonstrate the theoretical results in Section 4. Section 5 draws the conclusion.

2.   Existence, uniqueness and boundedness of positive solution

Considering the practical meaning, the first concern is whether there exists a global positive and bounded solution to system (1.3). That is to verify the well-posedness of system (1.3). In this section, we shall discuss this issue by the means of Lyapunov analysis method ^[23]. Denote

Theorem 2.1. For any initial condition (1.4), there exists a unique global solution to system (1.3) on $t > 0$. Moreover the solution remains in $\Omega$ almost surely (a.s.) for all times.

Proof. Because of the locally Lipschitzian continuity of the coefficients in system (1.3), for any initial condition (1.4) there exists a unique local solution on $t\in[-\tau, \tau_e)$, where $\tau_e$ is the explosion time ^[24]. By summing the equations of system (1.3), we get

Noting that $A = \gamma$ and initial condition (1.4), it follows that

To verify that the solution is global, we only need to show that $\tau_e = +\infty$. Now, we show this by proving a stronger property of the solution, i.e., it always remains in $\Omega$ a.s.. Define the stopping time as

where we let $\inf\emptyset = +\infty$. According to initial condition (1.4) and property (2.1), $\bar{\tau}$ is the first leaving time of the solution from $\Omega$. Clearly, $\bar{\tau}\leq \tau_e$. Thus, we only need to prove $\bar{\tau} = +\infty$, which implies that the solution remains in $\Omega$ for all times. Next, we will prove it by contradiction. If $\bar{\tau} < \infty$ a.s., there would exist a pair of constants $T > 0$ and $\epsilon\in (0, 1)$ such that $P(\bar{\tau}\leq T) > \epsilon$. Let's define a $C^2$-function $V:\Omega\rightarrow \mathbb{R}_+$ as

From the definition of $\bar{\tau}$, it is easily obtained that

for almost all $\omega\in \{\bar{\tau}\leq T\}$. On the other hand, using Itô's formula on $V$ for $t\in[0, \bar{\tau})$ and $\omega\in \{\bar{\tau}\leq T\}$, one has

where

Integrating both sides of (2.4) from $0$ to $t$ and then letting $t\rightarrow \bar{\tau}$, for almost all $\omega\in \{\bar{\tau}\leq T\}$, we get

which contradicts with (2.3). Thus, $\bar{\tau} = \tau_e = +\infty.$ This completes the proof.

By Theorem 2.1, if initial functions $(S_1(\theta), I_1(\theta), \cdots, S_n(\theta), I_n(\theta))\in\Omega$ for all $\theta\in [-\tau, 0]$, then

That is to say, the bounded region $\Omega$ is the almost surely positive invariant set of system (1.3). From now on, we always assume that $(S_1(t), I_1(t), \cdots, S_n(t), I_n(t))\in\Omega$.

3.   The dynamics of system

Since $S_k+I_k = 1$, as deterministic system (1.1), we analyze the following equivalent system of (1.3) instead

Denote $I = (I_1, \cdots, I_n)^\top$. System (1.3) can be represented by the vector-valued stochastic differential delay equation

where the $k$-th component of $f$ is $\lambda(k)(1-I_k)\Theta(t-\tau)e^{-d_m\tau}-\gamma I_k-\mu I_k$, $g$ is a diagonal matrix with entries $\sigma_k I_k(1-I_k)\Theta(t-\tau)e^{-d_m\tau}$, $k = 1, 2, \cdots, n$, and $B(t)$ is an n-dimensional standard Brownian motion with $B(0) = 0$. Define the differential generator associated with system (3.2) as

By Itô's formula, we have

Obviously, the disease-free equilibrium $I^0$ of system (1.2) is also that of stochastic system (3.1). In the following, the exponential stability of $I^0$ for system (3.1) will be deduced by Razumikhin-Mao type theorem ^[25].

Theorem 3.1. If $R_0 < 1$, then for any initial condition (1.4), the disease-free equilibrium $I^0$ is $p$-th $(p > 0)$ moment exponentially stable in $\mathcal{D}$ for system (3.1), and it is also almost surely exponentially stable in $\mathcal{D}$. Moreover,

where $q\in(1, \frac{1}{R_0})$ is the unique root of equation $(\gamma+\mu)(1-qR_0)\tau = \ln q$ and $|\cdot|$ represents the Euclidean norm of a vector or the trace norm of a matrix.

Proof. First, we prove the first moment exponential stability and a.s. exponential stability. Define the function $\Theta: \mathcal{D}\rightarrow \mathbb{R}_+$ as

Obviously,

where $\underline{\varrho} = \min\{\varphi(k)P(k), k = 1, \cdots, n\}$ and $\overline{\varrho} = \max\{\varphi(k)P(k), k = 1, \cdots, n\}.$ Using Cauchy-Schwarz inequality, we get

Noting that $0 < I_k < 1$, a direct calculation yields

By Theorem 6.4 of Mao ^[25], if $R_0 < 1$, the trivial solution $I^0$ of system (3.1) is the first moment exponentially stable.

For system (3.1), by triangle inequality, it can be easily obtained that

where $\overline{\lambda} = \max\{\lambda(k), k = 1, 2\cdots, n\}.$ Since $0 < I_k < 1$, one has $|1-I| < \sqrt{n}$. Combining with (3.6), it follows from (3.8) that

It holds that $I_k(1-I_k)\leq \frac{1}{4}$ for $0 < I_k < 1$, which together with (3.6) yields

By Theorem 6.4 of Mao ^[25], the trivial solution $I^0$ of system (3.1) is a.s. exponentially stable and the estimate (3.3) of its sample Lyapunov exponent holds.

Next, we prove the $p$-th moment $(p > 0)$ exponential stability. By Jensen's inequality, we get

Since $0 < I_k < 1$, then $\ln|I|^p/t < p\ln \sqrt n$ holds on $t\geq 1$. Using inverse Fatou's lemma, it follows that

which together with (3.3) and (3.9) implies the $p$-th moment exponential stability and the estimate (3.4). This completes the proof.

The condition for extinction of the disease has been derived in Theorem 3.1. Another problem we are interested in is when the disease will prevail. For deterministic system (1.2), as mentioned in the Introduction, when $R_0 > 1$ there is a globally stable endemic equilibrium $I^*$, which implies the prevalence. However, there is no endemic equilibrium for system (3.1). Next, we will study dynamics of system (3.1) around $I^*$ to reveal the persistence of the disease.

Theorem 3.2. If $R_0 > 1$, then for any initial condition (1.4), the solution of (3.1) has the property

where $\sigma_{max}^2 = max\{\sigma_k^2, k = 1, \cdots, n\}$, and $I^* = (I_1^*, \cdots, I_n^*)$ is the endemic equilibrium of corresponding deterministic system (1.2).

Proof. When $R_0 > 1$, the unique endemic equilibrium $I^*$ of deterministic system (1.2) satisfies

where $S_k^* = 1-I_k^*, \Theta^* = \frac{1}{\langle k\rangle}\sum_{j = 1}^n\varphi(j)P(j)I_j^*.$ Define a $C^2$-function $V:\mathcal{D}\rightarrow \mathbb{R}_+$ as

where $S_k = 1-I_k$ and

in which $x = S_k, \Theta$. Clearly, $V(I)$ is positive definite, i.e. $V(I^*) = 0$ and $V(I) > 0, \; I\neq I^*$.

Applying Itô's formula to $V$ along the solution of system (3.1), together with (3.11) and the property $S_k+I_k = 1$, we have

where

Since (3.11) holds, one has

Using Itô's formula on $V_\Theta$ and then substituting (3.13) into it, it follows that

where

From formulas (3.12) and (3.14), it can be obtained that

where $H(x) = x-1-\ln x.$ Since $0 < S_k < 1$, then

Using the fact that $H(x)\geq 0$ for $x > 0$ and the property $0 < S_k, I_k < 1$, together with (3.15), it follows that

Thus,

Integrating both sides of (3.16) from 0 to $t$ yields

Obviously, $M_\Theta$ is a continuous local martingale with $M(0) = 0$. By formula (3.15) and the property $0 < I_k < 1$, we obtain that

By Strong Law of Large Numbers ^[25], we get

Similarly,

Because of the positivity of $V$, it follows from (3.17)-(3.19) that

which implies the property (3.10). Thus, Theorem 3.2 is proved.

Remark 3.1. Theorem 3.2 shows that, when $R_0 > 1$, the Euclidian distance between the solution $I(t)$ and the endemic equilibrium $I^*$ of system (1.2) in time average takes the following form

where $C$ is a positive constant. It indicates that the solution of system (3.1) oscillates around $I^*$. The amplitude is correlated positively with the noise intensities but negatively with the time delay. It is reasonable that the solution is approximately stable, given that the disturbance intensities are sufficiently small. Under this assumption, we infer that the disease will prevail.

We denote by $\bar{I}(t)$ the average relative density of infected nodes at time t, then

From the result of Theorem 3.2, we conclude that system (3.1) is persistent, which also reflects that the disease is prevalent.

Definition 3.1. ^[26] System (3.1) is said to be persistent in the mean, if

Corollary 3.1. If $R_0 > 1$ and

then system (3.1) is persistent in the mean.

Proof. Obviously,

which together with (3.10) and (3.21) yields

Consequently,

where $\vartheta = \max\{\varphi(k)I_k^*, k = 1, \cdots, n\}$. This completes the proof.

4.   Simulation

In this section, numerical simulations of system (3.1) are shown to illustrate the theoretical results aforesaid. Moreover, we numerically simulate the solution of corresponding deterministic system (1.2) for comparison.

Consider a finite scale-free network whose maximum connectivity of any node $n = 100$. The degree distribution of the network is $P(k) = Ck^{-r}$ with $r = 2.3$, $\sum_{k = 1}^n P(k) = 1$. We fix the parameters $\gamma = 0.06, \mu = 0.05, d_m = 0.2.$ Let $\lambda(k) = \lambda k$ and $\varphi(k) = a k^\alpha/(1+bk^\alpha)$ ^[27] with $\alpha = 0.75, a = 0.8, b = 0.01$. The initial functions are $I_k(t) = 0.05, k = 1, \cdots, n$ for $t\in[-\tau, 0].$ By the Milstein's method ^[28], the discretized difference equations of system (3.1) are represented by

where $m = \tau/\Delta t$, $\xi_{k, i}\; (k = 1, \cdots, n, \; i = 1, \cdots, N)$ are independent standard normal random variables, and

First, we consider the values of $\lambda$ and $\tau$ such that $R_0 < 1$. Fig. 1 (imaginary lines) shows that the solution of deterministic system (1.2) converges to zero. Besides, Figure 1 (solid lines) depicts the dynamical behaviors of infected nodes with degree $k = 10, \; 35, \; 80$ by computing one sample path of the solution to system (3.1). The noise intensities $\sigma_k = 10$ and $\sigma_k = 30, k = 1, \cdots, 100$ in (a) and (b), respectively. It can be seen that the noise intensities do not affect the ultimate trend of the solution. It confirms the stability result in Theorem 3.1.

Second, we choose the values of $\lambda$ and $\tau$ such that $R_0 > 1$. Then there is an endemic equilibrium $I^*$ for system (1.2) and it is globally asymptotically stable, which is illustrated in Figure 2 (imaginary lines). Figure 2 (solid lines) shows one sample path of the solution to system (3.1) with the same $\tau = 3$ but different noise intensities. Specifically, the noise intensities $\sigma_k = 0.5$ in (a) and $\sigma_k = 1$ in (b) for $k = 1, \cdots, 100$, which all meet condition (3.21). As expected, we can recognize the behavior in Theorem 3.2, that the solution of system (3.1) fluctuates around $I^*$. Moreover, the oscillation amplitude is positively correlated with the noise intensities.

Finally, we simulate dynamical behaviors of system (3.1) in the case where $R_0 > 1$ and the condition (3.21) is unsatisfied. In Figure 3(a), the average relative density $\bar{I}(t)$ of system (3.1) does not fluctuate around that of deterministic system (1.2). Its fluctuation center seems less than the infection level of system (1.2). Figure 3(b) shows the density $\bar{I}(t)$ of system (3.1) averaged over 100 sample paths. It is not in excess of the infection level of system (1.2). In other words, the infection level of system (1.2) provides an upper bound for it. Therefore, we surmise that the strong noise intensities may reduce the infection size in statistical average meaning.

In addition, we study the influence of the time delay. In Figure 4 (solid lines), we depict one sample path of the solution with the same noise intensities but different time delays. By comparing Figure 2(b), Figure 4(a) and Figure 4(b), it is found out that when $R_0 > 1$, the oscillation amplitude reduces with the time delay increasing. Moreover, as it continues to increase such that $R_0 < 1$, the disease will die out.

5.   Conclusion

An SIS epidemic model with time delay and stochastic perturbation on scale-free networks is established in this paper. Here, we introduce the incubation period of the disease in a vector's body as time delay and enter stochastic perturbation in the infection rate. We prove that it possesses a unique global solution that remains within $\Omega$ whenever it starts from this region. Then we mainly discuss the dynamics of the stochastic system. The basic reproduction number $R_0$ of the corresponding deterministic system is a critical parameter. The disease-free equilibrium $I^0$ is exponentially stable (a.s. and $p$-th moment) when $R_0 < 1$, which means rapid extinction of the disease. Compared with the deterministic system, stochastic perturbation does not change the final extinction of the disease. Conversely, we analyze permanence of the stochastic system. Under the conditions $R_0 > 1$ and sufficiently weak noise intensities, the solution of the stochastic system ultimately fluctuates around endemic equilibrium $I^*$, which implies permanence in the mean. Moreover, the oscillation amplitude gets smaller with the white noise intensities decreasing, while it becomes greater with the time delay decreasing.

These results reveal the combined influences of time delay and stochastic perturbation on the dynamics of SIS epidemic model on networks. One may analyze the same issue about SIR ^[6], SIQS ^[7] model on scale-free networks and so on. In particular, because of little research about the SIRA model ^[29,30] in heterogeneous networks, this model itself and influences of time delay and stochastic perturbation on it are all worth studying. These will be considered in our future work.

Acknowledgments

This research was supported by the Foundation for Basic Disciplines of Army Engineering University under Grant no. KYSZJQZL2011.

Conflict of interest

The authors declare that they have no conflict of interest.