Global stability of multi-group SEIQR epidemic models with stochastic perturbation in computer network

1.   Introduction

As is well known, the internet world has brought great changes in the society. In reality, we know that cyber world is being threatened by the attack of malicious objects. Malicious object is a code that infects computer systems. There are different kinds of malicious objects such as: Worm, Virus, Trojan horse, etc., which differ according to the way they attack computer systems and the malicious actions they perform (see ^[1,2,3]). With the development of the computer network, malicious objects be widely spread through a network, through an online service, through shared computer software or through a mobile storage tool, and so on. Because of the similarity between the transmission of human infectious diseases and transmission of malicious objects in the computer network, some authors employ the epidemic models to describe the transmission of malicious objects in the cyber world (see ^{[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]}).

Considering different contact patterns, different anti-virus software, or distinct number of contacts etc., it is more appropriate to divide individual hosts into groups in modeling epidemic disease. Therefore, it is reasonable to propose multi-group models to describe the transmission dynamics of malicious objects in heterogeneous host populations on computer network. At present, many scholars have focused their study on various forms of multi-group epidemic models (see ^{[18,19,20,21,22,23]}). They have also proved the global stability of the unique endemic equilibrium through Lyapunov function, which is one of the main mathematical challenges in analyzing multi-group models. Particularly, Wang et al. ^[23] proposed the following multi-group SEIQR epidemic model for describing the transmission of malicious objects in computer network

where the total network nodes are divided into $n$ groups of nodes, $n\geq 2$ is an integer. $S_k(t)$, $E_k(t)$, $I_k(t)$, $Q_k(t)$ and $R_k(t)$ express the numbers of susceptible nodes, exposed (infected but not yet infectious) nodes, infectious nodes, quarantined nodes and recovered nodes at time $t$ in the $k$-th group $(1\leq k\leq n)$, respectively. The definitions of all parameters in model (1.1) are listed in Table 1. We assume that the parameters $d^{S}_{k}$, $d^{E}_{k}$, $d^{I}_{k}$, $d^{Q}_{k}$, $d^{R}_{k}$ and $\Lambda_{k}$ are positive and the rest of parameters in model (1.1) ia nonnegative for all $k$. In particular, $\beta_{kj} = 0$ if there is no transmission of the disease between compartments $S_{k}$ and $I_{j}$. In model (1.1), the basic reproduction number $R_{0} = \rho(M_{0})$, the spectral radius of matrix $M_{0} = (\frac{\beta_{kj}\epsilon_{k}\frac{\Lambda_k}{d_k}}{(d_{k}^{E} +\epsilon_{k})(d_{k}^{I}+\alpha_{k}+\delta_{k}+\gamma_{k})})_{n\times n}$, is a threshold which completely determines the persistence or extinction of the disease. It is shown that, if $R_{0}\leq1$, the disease-free equilibrium $E _{0}$ is globally stable in the feasible region and the disease always dies out, and if $R_{0} > 1$, a unique endemic equilibrium $E^{ *}$ exists and is globally stable in the interior of the feasible region, and once the disease appears, it eventually persists at the unique endemic equilibrium level.

On the other hand, there exist uncertainties and random phenomena everywhere in nature ^{[23,24,25,26,27]}. Environmental noises are usually considered to be harmful, which will lead to the disorder of the dynamics ^[20,21]. Nevertheless, the noises also play a positive role in the dynamics of complex nonlinear systems, especially in interdisciplinary physical models and biomathematics models, such as noise induced resonances, noise enhanced stability (NES) and so on ^{[22,23,24,28,29,30]}. According to the noise source, the noises can be divided into the additive noise and the multiplicative noise. The former is not controlled by the system and can be directly introduced to the system, while the latter is related to system parameters and variables. The multiplicative noises can always ensure the nonnegativity of the solution. The two main peculiarities of the presence of the multiplicative noise are the presence of the absorbing barrier in zero population density and the phenomenon of the anomalous fluctuations ^[25,31]. The noise existing in biological systems is caused by environmental fluctuations, which is usually considered as the multiplicative white noise. For example, Caruso et al. ^[26] described the dynamic behavior of an ecosystem of two competing species by a stochastic Lotka-Volterra model with the multiplicative white noise. The multiplicative noise models the interaction between the environment and the species.

For human disease related epidemics, the nature of epidemic growth and spread is random due to the unpredictability in person to person contacts. Because of environmental noises, the deterministic approach has some limitations in the mathematical modeling transmission of an infectious disease, several authors began to consider the effect of white noise on the computer network systems (see ^{[23,24,25,26,27]}).

There are different approaches used in the literature to introduce random perturbations into population models, both from a mathematical and biological perspective (see ^{[23,24,25,26,27,28,29,31]}). One is to perturb the positive equilibria in order for making robust the equilibria of deterministic models. In this situation, the essence of the investigation using the approach is to check if the asymptotic stability of the positive equilibria of deterministic models can be preserved. For example, Wang et al. ^[23] investigated a multi-group SEIQR model with random perturbation around the positive equilibrium of corresponding deterministic model, which revealed that the stochastic stability of endemic equilibrium depends on the magnitude of the intensity of noise as well as the parameters involved within the model. The other important approach is with parameters perturbation. We find that there are many literatures on this approach, see ^[25,26,27] and the references cited therein. In epidemic models, the natural death rate and the disease transmission rate are two of the key parameters to disease transmission. And in the real situation, the natural death rate and the disease transmission rate always fluctuate around some average value due to continuous fluctuation in the environment. For example, El Ansari et al. ^[25] considered a stochastic version of model (1.1) with noises introduced in the rate at which nodes are crashed due to reasons other than the attacks of viruses and the transmission rate, and they proved the various conditions that control the extinction and stability of a nonlinear mathematical spread model with stochastic perturbations.

We now turn to a continuous time SEIQRS model which takes random effects into account. In SEIQRS model (1.1), the natural death rate $d_{k}^{X_{i}}$, where $1\leq k\leq n$ and $(X_1, X_2, X_3, X_4, X_5) = (S, E, I, Q, R)$, is one of the key parameters to disease transmission. May ^[30] pointed out that all the parameters involved in the population model exhibit random fluctuation as the factors controlling them are not constant. And in the real situation, the natural death rate d always fluctuate around some average value due to continuous fluctuation in the environment. In this sense, $d_{k}^{X_{i}}$ can seem as a random variable $\tilde{d}_{k}^{X_{i}}$. More precisely, in $[t, t + dt)$,

where $B_{ik}(t)\; (1\leq k\leq n, i = 1, 2, 3, 4, 5)$ are the independent standard Brownian motion defined on the complete probability space $(\Omega, \{\mathcal{F}_{t}\}_{t\geq0}, P)$ with a filtration $\{\mathcal{F}_{t}\}_{t\geq0}$ satisfying the usual conditions, and $\sigma_{ik}^2$ is the intensity of $B_{ik}(t)$. The reason of adopting $\sigma^2_{ik}\; (1\leq k\leq n, i = 1, 2, 3, 4, 5)$ as the intensity of the noise for the group $S_k$, $E_k$, $I_k$, $Q_k$ and $R_k$, respectively, is considering the difference between the group mobility response to infection risks. And then, in $[t, t + dt)$, $-\tilde{d}_{k}^{X_{i}}dt$ is normally distributed with mean $\mathbb{E}(-\tilde{d}_{k}^{X_{i}}dt) = -d_{k}^{X_{i}}dt$ and variance ${\bf Var}(-\tilde{d}_{k}^{X_{i}}dt) = \sigma_{i}^{2}dt$. Due to ${\bf Var}(-\tilde{d}_{k}^{X_{i}}dt) = \sigma_{i}^{2}dt\rightarrow0$ as $dt\rightarrow0$, this is a biologically reasonable assumption. Indeed this is a well-established way of introducing stochastic environmental noise into biologically realistic population dynamic models.

Therefore, replace $-d_{k}^{X_{i}}dt$ in model (1.1) with $-\tilde{d}_{k}^{X_{i}}dt = -d_{k}^{X_{i}}dt+ \sigma_{ik}dB_{ik}(t) \; (1\leq k\leq n, i = 1, 2, 3, 4, 5)$, and for simplicity, we replace $-\tilde{d}_{k}^{X_{i}}$ with $d_{k}^{X_{i}}$ again, then we can obtain the same SDE epidemic model as the following model (1.2) that is analog to its deterministic version model (1.1) by introducing stochastic perturbation terms to the growth equations of susceptible, infectious, recovered individuals to incorporate the effect of randomly fluctuating environments:

Throughout this paper, we always assume that model (1.2) is defined on a complete probability space $(\Omega, \{\mathcal{F}_{t}\}_{t\geq0}, P)$ with a filtration $\{\mathcal{F}_{t}\}_{t\geq0}$ satisfying the usual conditions (i.e., it is right continuous and $\mathcal{F}_{0}$ contain all P-null sets). Furthermore, we also always assume that the infection rate matrix $B = (\beta_{kj})_{n\times n}$ in model (1.2) is irreducible.

In this paper, we will study the asymptotic behavior of positive solutions of model (1.2) around the disease-free and endemic equilibria of corresponding deterministic model (1.1) in probability meaning by using the theory of graphs, Lyapunov functions method, It$\hat{o}$'s formula and the theory of stochastic analysis. Then by using the theory of stationary distributions of stochastic process we will study the existence of stationary distribution of model (1.2).

The paper is organized as follows. In Section 2, the criterion on the asymptotic behavior of positive solutions of model (1.2) around the disease-free equilibrium of the corresponding deterministic model is stated and proved. In Section 3, the sufficient condition the asymptotic behavior of positive solutions of model (1.2) around the endemic equilibrium of corresponding deterministic model and the existence of stationary distribution are stated and proved. In Section 4, we make some numerical simulations to illustrate our analytical results. Finally, in Section 5, we give a brief conclusion.

2.   Asymptotic behavior around disease-free equilibrium of model (1.1)

We first give a lemma to show that for any positive initial value model (1.2) has a unique positive solution defined on $[0, \infty)$.

Lemma 1. For any initial value in $R_{+}^{5n}$ model (1.2) has a unique positive solution defined for all $t\geq 0$ and the solution remain in $R_{+}^{5n}$ with probability one.

This lemma can be easily proved by using the standard arguments as in ^[14,18] and with the help of Lyapunov function

where positive constant $a$ satisfies $a\leq\min\{\frac{d_{k}^{I}+\alpha_{k}}{\sum_{j = 1}^{n}\beta_{jk}}, k = 1, 2, \cdots, n\}$.

For deterministic model (1.1), in ^[23] the authors have obtained that there is a disease-free equilibrium $E_{0} = (S_{1}^{0}, 0, 0, 0, 0, S_{2}^{0}, 0, 0, 0, 0, \cdots, S_{n}^{0}, 0, 0, 0, 0)$, where $S_{k}^{0} = \frac{\Lambda_{k}}{d_{k}^{S}}$, and if $R_{0}\leq 1$, then $E_{0}$ is globally asymptotically stable, which means the disease will die out. Therefore, it is interesting to study the stability of disease-free equilibrium for controlling the spread of infectious disease. However, for stochastic model (1.2) there is not any disease-free equilibrium. Therefore, it is natural to ask how we can consider the disease will be extinct. In this section we mainly through estimating the asymptotic oscillation around equilibrium $E_{0}$ of any positive solutions of stochastic model (1.2) to reflect whether the disease in stochastic model (1.2) will die out. We have the following result.

Theorem 1. Assume that $R_{0}\leq1$ and the following conditions hold

Then for any positive solution $(S_{k}(t), E_{k}(t), I_{k}(t), Q_{k}(t), R_{k}(t), 1\leq k\leq n)$ of model (1.2) one has

where $A_{k} = (d_{k}^{S}-\sigma_{1k}^{2})$, $B_{k} = \frac{1}{4}(d_{k}^{E}+\epsilon_{k}-\frac{1}{2}\sigma_{2k}^{2})$ and

and positive constants $a_k, d_k, c_k, e_k\; (1\leq k\leq n)$ and $b$ will be confirmed in the proof of the theorem.

Proof. Let $u_{k} = S_{k}-\frac{\Lambda_{k}}{d_{k}^{S}}, v_{k} = E_{k}, w_{k} = I_{k}, y_{k} = Q_{k}, z_{k} = R_{k}\, (1\leq k\leq n)$, then model (1.2) becomes into

Since $B = (\beta_{kj})_{n\times n}$ is irreducible, then $M_{0}$ is also nonnegative and irreducible. Hence, by Lemma A.1 in ^[3], $M_{0}$ has a positive left eigenvector $\eta = (\eta_{1}, \eta_{2}, \cdots, \eta_{n})$ such that

Define a Lyapunov function as follows.

with $V_1 = \frac{1}{2}\sum_{k = 1}^{n}(u_{k}+v_{k})^{2}$, $V_2 = \frac{1}{2}\sum_{k = 1}^{n}a_{k}u_{k}^{2}$, $V_3 = \sum_{k = 1}^{n}\frac{\epsilon_{k}\eta_{k}}{(d_{k}^{E}+\epsilon_{k})(d_{k}^{I}+\alpha_{k}+\delta_{k}+\gamma_{k})} (v_{k}+\frac{d_{k}^{E}+\epsilon_{k}}{\epsilon_{k}}w_{k})$, $V_4 = \sum_{k = 1}^{n}c_{k}w_{k}^{2}$, $V_5 = \sum_{k = 1}^{n}d_{k}y_{k}^{2}$ and $V_6 = \sum_{k = 1}^{n}e_{k}z_{k}^{2},$ where positive constants $a_k, c_k, d_k, e_k\, (1\leq k\leq n)$ and $b$ will be determined later. By It$\hat{o}$'s formula, we get

with $LV = LV_1+b(LV_2+LV_3)+LV_4+LV_5+LV_6,$ where

and

Note from (2.2) that

where $w = (w_1, w_2, \cdots, w_n)^T$. If $R_{0}\leq 1$, then

Furthermore, we also have

and

Choosing $a_{k} = \frac{d_{k}^{S}\eta_{k}\epsilon_{k}}{(d_{k}^{E}+\epsilon_{k})(d_{k}^{I}+\alpha_{k}+ \delta_{k}+\gamma_{k})\Lambda_{k}}\; (1\leq k\leq n)$ and $b = \max_{1\leq k\leq n}\{\frac{(d_{k}^{S}+d_{k}^{E} +\epsilon_{k})^{2}}{2a_{k}(d_{k}^{E} +\epsilon_{k}-\frac{1}{2}\sigma_{2k}^{2})}\},$ then from (2.4)–(2.8) we finally obtain

where $A_k, B_k, C_k, D_k$ and $F_k$ are given in the above.

If (2.1) holds, then $A_{k} > 0$, $B_{k} > 0$ and $D_k > 0$. Further, we can choose $c_{k}$, $d_{k}$ and $e_{k}$ such that

Particularly, we can take

where $\eta_k = \frac{\delta_{k}^{2}}{d_{k}^{Q}+\alpha_{k}+\mu_{k} -\frac{1}{2}\sigma_{4k}^{2}}+\frac{\gamma_{k}^{2}(d_{k}^{Q}+\alpha_{k}+\mu_{k} -\frac{1}{2}\sigma_{4k}^{2})}{\mu_{k}^{2}} > 0.$ Thus, we have

and $F_{k} = e_{k}(d_{k}^{R}-\frac{1}{2}\sigma_{5k}^{2}-\frac{e_{k}\mu_{k}^{2}}{d_{k}(d_{k}^{Q}+\alpha_{k}+\mu_{k} -\frac{1}{2}\sigma_{4k}^{2})}) > 0.$ By integration and taking expectation of both sides of (2.3), from (2.9) we obtain

Therefore,

Consequently,

This completes the proof. □

Remark 1. From Theorem 1, we see that under some conditions the solution of model (1.2) will oscillates around the disease-free equilibrium of deterministic model (1.1), and the intensity of fluctuation is only relation to the intensity of the white noise $B_{1k}(t)$, but do not relation to the intensities of the other white noises. In a biological interpretation, as the intensity of stochastic perturbations is small, the solution of model (1.2) will be close to the disease-free equilibrium of model (1.1) most of the time.

As a special case of model (1.2), when $\sigma_{1k} = 0$, then model (1.2) becomes into

Obviously, $E_{0}$ is also the disease-free equilibrium of model (2.10). From the proof of Theorem 2, we get

which is negative definite if for each $1\leq k\leq n$

Therefore, as a consequence of Theorem 1 we have the following result.

Corollary 1. Assume that $R_{0}\leq 1$ and condition (2.11) holds. Then disease-free equilibrium $E_0$ of model (2.9) is globally stochastically asymptotically stable.

3.   Asymptotic behavior around endemic equilibrium of model (1.1)

Firstly, we introduce some concepts and conclusions of graph theory (see ^[10]). A directed graph $g = (V, E)$ contains a set $V = \{1, 2, \cdots, n\}$ of vertices and a set $E$ of arcs $(k, j)$ leading from initial vertex $k$ to terminal vertex $j$. $A$ subgraph $H$ of $g$ is said to be spanning if $H$ and $g$ have the same vertex set. $A$ directed digraph $g$ is weighted if each arc $(k, j)$ is assigned a positive weight $a_{kj}$. Given a weighted digraph $g$ with $n$ vertices, define the weight matrix $A = (a_{kj})_{n\times n}$ whose entry $a_{kj}$ equals the weight of arc $(k, j)$ if it exists, and 0 otherwise. $A$ weighted digraph is denoted by $(g, A)$. $A$ digraph $g$ is strongly connected if for any pair of distinct vertices, there exists a directed path from one to the other and it is well known that a weighted digraph $(g, A)$ is stronly connected if and only if the weight matrix $A$ is irreducible (see ^[32]).

The Laplacian matrix of graph $(g, A)$ is defined by

Let $c_{k}\, (1\leq k\leq n)$ denote the cofactor of the $k$-th diagonal element of $L_{A}$. The following lemmas are the classical results of graph theory (see ^[21,33]) which will be used in this paper.

Lemma 2. Assume that $A$ is a irreducible matrix and $n\geq 2$. Then $c_{k} > 0$ for all $1\leq k\leq n$.

Lemma 3. Assume that $A$ is a irreducible matrix and $n\geq 2$. Then the following equality holds

where $G_{k}(x_{k})\, (1\leq k\leq n)$ are arbitrary functions.

For model (1.2), we see that there is not any endemic equilibrium. Therefore, in order to study the persistence of disease in model (1.2), we need to study the asymptotic behavior of the endemic equilibrium of model (1.2) which is surrounding the deterministic model (1.1), we obtain the following result.

Theorem 2. Assume that $R_{0} > 1$ and the following conditions hold

Then for any positive solution $(S_{k}(t), E_{k}(t), I_{k}(t), Q_{k}(t), R_{k}(t), 1\leq k\leq n)$ of model (1.2) one has

where $E^* = (S_k^*, E_k^*, I_k^*, Q_k^*, R_k^*, 1\leq k\leq n)$ be the endemic equilibrium of model (1.1), and

and positive constants $r$, $a_k, b_k, c_k$ and $D_k\; (1\leq k\leq n)$ will be confirmed in the proof of the theorem.

Proof. When $R_{0} > 1$, from ^[23] there exits an endemic equilibrium $E^{*}$ of model (1.1), then

Let matrix $A = (a_{kj})_{n\times n}$ with $a_{kj} = \beta_{kj}S_k^*I_j^*, \; k, j = 1, 2, \cdots, n.$ Since $B = (\beta_{kj})_{n\times n}$ is irreducible, then $A$ also is irreducible.

Firstly, define the $C^{2}$-function $V_{1}:R_{+}^{3n}\rightarrow R_{+}$ by

where $c_k\, (1\leq k\leq n)$ are the cofactor of the $k$-th diagonal element of $L_A$. $V_{1}$ is positive definite. From It$\hat{o}$'s formula, by calculating we can get

By Lemma 2, we obtain

Similarly, we also get

Hence

where the last equality is derived from Lemma 3. Substituting (3.3) and (3.4) into (3.2), we have

Secondly, define the $C^{2}$-function $V_{2}:R_{+}^{2n}\rightarrow R_{+}$ as follows.

where $c_k\, (1\leq k\leq n)$ are given as in $V_1$. $V_{2}$ is positive definite. It follows from It$\hat{o}$'s formula that

We have

where the last equality is derived from Lemma 3 such that

We further get

Substituting (3.7) and (3.8) into (3.6), we have

Thirdly, define the $C^{2}$-function $V_{3}:R_{+}^{n}\rightarrow R_{+}$ by

where $c_k\, (1\leq k\leq n)$ are given as in $V_1$. We obtain

Choose $K = \sum_{j = 1}^{n}\beta_{kj}\frac{I_{k}^{*}}{d_{k}^{S}}$, then (3.5) together with (3.9) and (3.10) implies

where $A_{k} = \frac{1}{2} \sum_{k = 1}^{n}c_{k}[(K+2)\sigma_{1k}^{2}S^{*}_{k}+(K+1)\sigma_{2k}^{2}E^{*}_{k} +(K+1)\frac{d_{k}^{E}+\epsilon_{k}}{\epsilon_{k}}\sigma_{3k}^{2}I^{*}_{k}].$

Next, define the $C^{2}$-function $V_{4}:R_{+}^{3n}\rightarrow R_{+}$ by

where $a_{k}\, (1\leq k\leq n)$ are positive constants to be determined later. By calculating, we can get

Since $2(d_{k}^{S}+d_{k}^{E})(S_{k}-S_{k}^{*})(E_{k}-E_{k}^{*})\leq \frac{(d_{k}^{S}+d_{k}^{E})^{2}}{d_{k}^{E}}(S_{k}-S_{k}^{*})^{2}+d_{k}^{E}(E_{k}-E_{k}^{*})^{2}$ and

we further obtain

where $D_{k} = d_{k}^{S}+d_{k}^{E}+\frac{(d_{k}^{S})^{2}}{d_{k}^{E}}+\frac{(d_{k}^{S}+d_{k}^{I}+\alpha_{k}+\delta_{k} +\gamma_{k})^{2}}{d_{k}^{I}+\alpha_{k}+\delta_{k}+\gamma_{k}} +\sigma_{1k}^{2}.$

Further, define the $C^{2}$-function $V_{5}:R_{+}^{n}\rightarrow R_{+}$ by

We obtain

Finally, define the $C^{2}$ functions $V_{6}$ and $V_{7}:R_{+}^{n}\rightarrow R_{+}$ as follows.

where $b_{k}$, $d_{k}\, (1\leq k\leq n)$ are positive constants to be determined later. We get

and

where the last equality is derived by the inequality $2ab\leq a^{2}+b^{2}$.

From (3.12)–(3.15) we obtain

where

From condition (3.1), we can choose positive numbers $r$, $a_k$, $b_k$ and $d_k$ for $k = 1, 2, \cdots, n$ satisfying $d_k < d_{k}^{R}-2\sigma_{5k}^{2}$ and

such that for each $1\leq k\leq n$

Lastly, define a Lyapunov function as follows

By It$\hat{o}$'s formula, we obtain

where (3.11) together with (3.16) implies

By integration and taking expectation of both sides of (3.17), we obtain

Therefore,

This completes the proof. □

As a consequence of Theorem 2, we have the following result on the existence and uniqueness of stationary distribution for model (1.2).

Theorem 3. Assume that all conditions in Theorem 2 hold. Then model (1.2) has a unique stationary distribution $\mu(\cdot)$ in $R_{+}^{5n}$.

Proof. Choose region $\Omega$ in (^[34], Lemma 2.5) by $\Omega = R_{+}^{5n}$. Consider the following inequality

Let region $U_1$ denote all points $(S_k, E_k, I_k, Q_k, R_k, 1\leq k\leq n)$ which satisfy the above inequality with $H = 2\sum_{k = 1}^n\rho_k$ and region $U_2$ denote all points $(S_k, E_k, I_k, Q_k, R_k, 1\leq k\leq n)$ which satisfy the above inequality with $H = 3\sum_{k = 1}^n\rho_k$. Obviously, $U_2$ is a neighborhood of $U_1$ and the closure $\bar{U}_2\subset \Omega$. Then from (3.18), for any $x\in \Omega\setminus U_1$,

which implies condition (ⅱ) in (^[35], Lemma 2.5) is satisfied.

For model (1.2), the diffusion matrix is

Choose a positive constant $M\geq\inf_{\bar{U}_2}\{\sigma_{1i}^{2}S_{i}^{2}, \sigma_{2i}^{2}E_{i}^{2}, \sigma_{3i}^{2}I_{i}^{2}, \sigma_{4i}^{2}Q_{i}^{2}, \sigma_{5i}^{2}R_{i}^{2}, 1\leq i\leq n\}.$ Then,

for all $(S_{i}, E_{i}, I_{i}, Q_{i}, R_{i}, \; 1\leq i\leq n)\in\bar{U}_2$ and $\xi\in R^{5n}$. This implies condition (ⅰ) in (^[34], Lemma 2.5) is also satisfied. Therefore, by (^[34], Lemma 2.5), model (1.2) has a unique stationary distribution $\mu$ in $R_{+}^{5n}$. This completes the proof. □

4.   Numerical simulation

In this section, we analyse the stochastic behaviour of model (1.2) by means of the numerical simulations in order to make readers understand our results more better. The numerical simulation method can be found in ^[36]. The corresponding discretization system of

where time increment $\Delta t > 0$, and $\varepsilon_{1k, i}$, $\varepsilon_{2k, i}$, $\varepsilon_{3k, i}$, $\varepsilon_{4k, i}$, $\varepsilon_{5k, i}$ for $1\leq k\leq n$ are $N(0, 1)$-distributed independent random variables which be generated numerically by pseudo-random number generators.

Example 1. In model (1.2), we choose $n = 2$ and the parameters $\Lambda_{1} = 3.2$, $\epsilon_{1} = 0.1$, $\alpha_{1} = 0.1$, $\beta_{11} = 0.409$, $d_{1}^{S} = 0.9$, $d_{1}^{E} = 0.7$, $d_{1}^{I} = 0.81$, $d_{1}^{Q} = 0.2$, $d_{1}^{R} = 0.65$, $\mu_{1} = 0.3$, $\gamma_{1} = 0.04$, $\beta_{12} = 0.02$, $\delta_{1} = 0.1$, $\sigma_{11} = 0.15$, $\sigma_{21} = 0.1$, $\sigma_{31} = 0.41$, $\sigma_{41} = 0.2$, $\sigma_{51} = 0.3$, $\Lambda_{2} = 7.5$, $\epsilon_{2} = 2.4$, $\alpha_{2} = 0.2$, $\beta_{21} = 0.05$, $d_{2}^{S} = 0.49$, $d_{2}^{E} = 0.25$, $d_{2}^{I} = 0.15$, $d_{2}^{Q} = 0.25$, $d_{2}^{R} = 0.39$, $\mu_{2} = 0.5$, $\gamma_{2} = 0.15$, $\beta_{22} = 0.0014$, $\delta_{2} = 0.43$, $\sigma_{12} = 0.2$, $\sigma_{22} = 0.6$, $\sigma_{32} = 0.5$, $\sigma_{42} = 0.8$ and $\sigma_{52} = 0.8$.

By computing, we have $R_{0}\doteq0.8675 < 1$ and disease-free equilibrium $E_0 = (3.56, 0, 0, 0, 0, 15.31, 0, 0, 0, 0)$ for corresponding deterministic model (1.1), and the conditions in Theorem 1 are satisfied. Therefore, according to the conclusion in Theorem 1 by numerical calculation we can obtain that for the solution $(S_k(t), E_k(t), I_k(t), Q_k(t), R_k(t), k = 1, 2)$ satisfying the initial values $(S_1(0), E_1(0), I_1(0), Q_1(0), R_1(0)) = (0.75, 0.8, 0.8, 0.2, 0.2)$ and $(S_2(0), E_2(0), I_2(0), Q_2(0), R_2(0)) = (1.7, 4.5, 2.7, 4.3, 5)$ one has

where $S_1^0 = 3.56$, $S_2^0 = 15.31$, $A_1 = 0.8775$, $A_2 = 0.48$, $B_{1} = 0.1988$, $B_2 = 0.6175$, $C_1 = 18.21$, $C_2 = 0.04295$, $D_1 = 0.1482$, $D_2 = 0.0077$, $F_1 = 0.1507$ and $F_2 = 3.7551\times10^{5}$.

From the numerical simulations given in Figure 1 we easily see that the above formula (4.1) holds. That is, the solution of stochastic model (1.2) asymptotically oscillates in probability around disease-free equilibrium $E_0$.

In addition, from Figure 1 we also easily see that the mean of susceptible $S_k (t)\; (k = 1, 2)$ tend to $S_k^0$ and all exposed $E_k$, infectious $I_k$, quarantined $Q_k$ and recovered $R_k$ for $k = 1, 2$ tend to zero in probability as $t\to\infty$.

Example 2. In model (1.2), we choose $n = 2$ and the parameters $\Lambda_{1} = 0.8$, $\epsilon_{1} = 0.1$, $\alpha_{1} = 0.1$, $\beta_{11} = 0.109$, $d_{1}^{S} = 0.19$, $d_{1}^{E} = 1.107$, $d_{1}^{I} = 0.081$, $d_{1}^{Q} = 0.2$, $d_{1}^{R} = 0.65$, $\mu_{1} = 0.3$, $\gamma_{1} = 0.04$, $\beta_{12} = 0.02$, $\delta_{1} = 0.01$, $\sigma_{11} = 1.15$, $\sigma_{21} = 1.1$, $\sigma_{31} = 1.41$, $\sigma_{41} = 01.2$, $\sigma_{51} = 1.3$, $\Lambda_{2} = 1.5$, $\epsilon_{2} = 2.4$, $\alpha_{2} = 0.2$, $\beta_{21} = 0.05$, $d_{2}^{S} = 0.49$, $d_{2}^{E} = 0.25$, $d_{2}^{I} = 0.15$, $d_{2}^{Q} = 0.25$, $d_{2}^{R} = 0.39$, $\mu_{2} = 0.5$, $\gamma_{2} = 0.15$, $\beta_{22} = 0.0014$, $\delta_{2} = 0.043$, $\sigma_{12} = 1.2$, $\sigma_{22} = 1.6$, $\sigma_{32} = 0.5$, $\sigma_{42} = 0.8$ and $\sigma_{52} = 0.8$.

By computing, we have $R_{0}\doteq 0.5174\leq1$. Since $d_{1}^{S}-\sigma_{11}^{2} = -1.13 < 0$, $d_{2}^{S}-\sigma_{12}^{2} = -0.33 < 0$, $d_{1}^{R}-\frac{1}{2}\sigma_{51}^{2} = -0.2 < 0$ and $d_{2}^{R}-\frac{1}{2}\sigma_{52}^{2} = -0.46 < 0$, the condition (2.1) in Theorem 1 does not hold. However, from the numerical simulations given in Figure 2, we can see that the solution $(S_k(t), E_k(t), I_k(t), Q_k(t), R_k(t), k = 1, 2)$ of stochastic model (1.2) with initial values $(S_1(0), E_1(0), I_1(0), Q_1(0), R_1(0)) = (0.75, 0.8, 0.8, 0.2, 0.2)$ and $(S_2(0), E_2(0), I_2(0), Q_2(0), R_2(0)) = (1.7, 4.5, 2.7, 4.3, 5)$ asymptotically oscillates in probability around the disease-free equilibrium $E_0 = (4.21, 0, 0, 0, 0, 3.06, 0, 0, 0, 0)$ of corresponding deterministic model (1.1). This example seems to indicate that the condition (2.1) in Theorem 1 can be weakened or taken out.

Example 3. In model (1.2), we choose $n = 2$ and the parameters $\Lambda_{1} = 4.5$, $\epsilon_{1} = 1$, $\alpha_{1} = 0.1$, $\beta_{11} = 1.55$, $d_{1}^{S} = 0.5$, $d_{1}^{E} = 0.15$, $d_{1}^{I} = 0.1$, $d_{1}^{Q} = 0.2$, $d_{1}^{R} = 0.65$, $\mu_{1} = 0.3$, $\gamma_{1} = 0.4$, $\beta_{12} = 1.35$, $\delta_{1} = 0.6$, $\sigma_{11} = 0.3$, $\sigma_{21} = 0.5$, $\sigma_{31} = 0.4$, $\sigma_{41} = 0.2$, $\sigma_{51} = 0.4$, $\Lambda_{2} = 7.5$, $\epsilon_{2} = 2.4$, $\alpha_{2} = 0.2$, $\beta_{21} = 1.5$, $d_{2}^{S} = 0.49$, $d_{2}^{E} = 0.25$, $d_{2}^{I} = 0.15$, $d_{2}^{Q} = 0.25$, $d_{2}^{R} = 0.39$, $\mu_{2} = 0.5$, $\gamma_{2} = 0.15$, $\beta_{22} = 1.24$, $\delta_{2} = 0.43$, $\sigma_{12} = 0.2$, $\sigma_{22} = 0.6$, $\sigma_{32} = 0.5$, $\sigma_{42} = 0.8$ and $\sigma_{52} = 0.3$.

By computing, we have $R_{0}\doteq1.1032 > 1$ and the conditions in Theorem 2 are satisfied. The numerical simulations are given in Figures 3 and 4. Figure 3 shows that the solution $(S_k(t), E_k(t), I_k(t), Q_k(t), R_k(t), k = 1, 2)$ of stochastic model (1.2) satisfying the initial values $(S_1(0), E_1(0), I_1(0), Q_1(0), R_1(0)) = (0.75, 0.8, 0.8, 0.2, 0.2)$ and $(S_2(0), E_2(0), I_2(0), Q_2(0), R_2(0)) = (1.7, 4.5, 2.7, 4.3, 5)$ asymptotically oscillates in probability around the endemic equilibrium $E^* = (0.37, 3.35, 0.27, 0.38, 0.19, 0.79, 2.68, 6.93, 3.14, 6.68)$ of corresponding deterministic model (1.1). Figure 4 shows that the solution has a unique stationary distribution. Therefore, the conclusions of Theorem 3 are validated by the numerical example.

In addition, from Figure 3 we also easily see that the mean value of the solution for stochastic model (1.2) asymptotically oscillates in probability around the endemic equilibrium $E^*$ of corresponding deterministic model (1.1). From Figure 5 we can find the relationship between variances of the solution $(S_k(t), E_k(t), I_k(t), Q_k(t), R_k(t), k = 1, 2)$ and the intensities of noises $(\sigma_{ik}^2, \sigma_{2k}^2, \sigma_{3k}^2, \sigma_{4k}^2, \sigma_{5k}^2, k = 1, 2)$ as time $t$ is enough large.

Example 4. In model (1.2), we choose $n = 2$ and the parameters $\Lambda_{1} = 4.5$, $\epsilon_{1} = 0.1$, $\alpha_{1} = 0.1$, $\beta_{11} = 1.55$, $d_{1}^{S} = 2.05$, $d_{1}^{E} = 1.015$, $d_{1}^{I} = 0.51$, $d_{1}^{Q} = 0.02$, $d_{1}^{R} = 0.65$, $\mu_{1} = 0.3$, $\gamma_{1} = 0.04$, $\beta_{12} = 1.35$, $\delta_{1} = 0.6$, $\sigma_{11} = 2.3$, $\sigma_{21} = 1.5$, $\sigma_{31} = 0.5$, $\sigma_{41} = 0.4$, $\sigma_{51} = 0.4$, $\Lambda_{2} = 7.5$, $\epsilon_{2} = 2.4$, $\alpha_{2} = 0.2$, $\beta_{21} = 1.5$, $d_{2}^{S} = 0.49$, $d_{2}^{E} = 0.25$, $d_{2}^{I} = 0.15$, $d_{2}^{Q} = 0.25$, $d_{2}^{R} = 0.39$, $\mu_{2} = 0.5$, $\gamma_{2} = 0.15$, $\beta_{22} = 0.24$, $\delta_{2} = 0.43$, $\sigma_{12} = 1.2$, $\sigma_{22} = 0.6$, $\sigma_{32} = 0.5$, $\sigma_{42} = 0.8$ and $\sigma_{52} = 0.3$.

By computing, we have $R_{0}\doteq 1.09013 > 1$. Since $d_{1}^{S}-\sigma_{11}^{2} = -10.12 < 0$, $d_{2}^{S}-\sigma_{21}^{2} = -0.43 < 0$ and $d_{1}^{E}-\frac{1}{2}\sigma_{21}^{2} = -0.11 < 0$, the condition (3.1) in Theorem 2 does not hold. However, from the numerical simulations are given in Figures 6 we can see that the solution $(S_k(t), E_k(t), I_k(t), Q_k(t), R_k(t), k = 1, 2)$ of stochastic model (1.2) with initial values $(S_1(0), E_1(0), I_1(0), Q_1(0), R_1(0)) = (0.75, 0.8, 0.8, 0.2, 0.2)$ and $(S_2(0), E_2(0), I_2(0), Q_2(0), R_2(0)) = (1.7, 4.5, 2.7, 4.3, 5)$ asymptotically oscillates in probability around the endemic equilibrium $E^* = (0.43, 3.24, 0.26, 0.37, 0.19, 3.34, 2.22, 5.7, 2.59, 5.51)$ of corresponding deterministic model (1.1). This example seems to indicate that the condition (3.1) in Theorem 2 can be weakened or taken out.

5.   Conclusions

In this research we consider a class of stochastic multi-group SEIQR (susceptible, exposed, infectious, quarantined and recovered) models in computer network. For the deterministic system, if the reproduction number $R_{0} > 1$, the system has unique endemic equilibrium which is globally stable, this means that the disease will persist at the endemic equilibrium level if it is initially present. It is clear that when the disease is endemic, the recovery nodes increases with the increasing quarantine nodes, and finally both reach the steady state values. Thus, it will be of great importance for one to run anti-malicious software to quarantine infected nodes. In order to study the asymptotic behavior of model (1.2), we first introduce the global existence of a positive solution. Then by using the theory of graphs, stochastic Lyapunov functions method, It$\hat{o}$'s formula and the theory of stochastic analysis, we carry out a detailed analysis on the asymptotic behavior of model (1.2). If $R_{0}\leq 1$, the solution of model (1.2) oscillates around the disease-free equilibrium, while if $R_{0} > 1$, the solution of model (1.2) fluctuates around the endemic equilibrium. The investigation of this stochastic model revealed that the stochastic stability of $E^{*}$ depends on the magnitude of the intensity of noise as well as the parameters involved within the model system. finally, numerical methods are employed to illustrate the dynamic behavior of the model. The effect of quarantine on recovered nodes is also analyzed in the stochastic model.

Some interesting topics deserve further consideration. On the one hand, we can solve the corresponding probability density function of various stochastic epidemic models. On the other hand, we need to establish a more complete and systematic theory to obtain more accurate conditions and density function. The reader is referred to ^{[37,38,39,40,41,42,43,44,45]}. These problems are expected to be studied and solved as planned future work.

Acknowledgments

This research is supported by the Natural Science Foundation of Xinjiang of China (Grant Nos. 2020D01C178) and the National Natural Science Foundation of China (Grant Nos. 12101529, 12061079, 72163033, 72174175, 11961071).

Conflict of interest

The authors declare there is no conflicts of interest.