Stability of Equilibria of Rumor Spreading Model under Stochastic Perturbations

Abstract

The known mathematical model of rumor spreading, which is described by a system of four nonlinear differential equations and is very popular in research, is considered. It is supposed that the considered model is influenced by stochastic perturbations that are of the type of white noise and are proportional to the deviation of the system state from its equilibrium point. Sufficient conditions of stability in probability for each from the five equilibria of the considered model are obtained by virtue of the Routh–Hurwitz criterion and the method of linear matrix inequalities (LMIs). The obtained results are illustrated by numerical analysis of appropriate LMIs and numerical simulations of solutions of the considered system of stochastic differential equations. The research method can also be used in other applications for similar nonlinear models with the order of nonlinearity higher than one.

1. Introduction

There are two classes of mathematical models of the type of epidemics: medical epidemics (see, for instance, the so-called SIR-epidemic model [1,2,3]) and different social epidemics (see, for instance, the alcohol consumption model [4] or the model of obesity epidemic [5]). During the last two decades, the rumor spreading model, that is an epidemic of the social type too, is extremely popular in research (see, [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]). Following [26], we will consider the rumor spreading model (the so-called I2SR-model) in the form

$$\begin{array}{l}\dot{I}\left(t\right)=p-{\lambda }_{1}I\left(t\right)S_1\left(t\right)-{\lambda }_{2}I\left(t\right)S_2\left(t\right)-qI\left(t\right),\\ {\dot{S}}_1\left(t\right)={\lambda }_{1}I\left(t\right)S_1\left(t\right)+\alpha S_2\left(t\right)-{\delta }_1{S}_1\left(t\right)R\left(t\right)-q{S}_1\left(t\right),\hfill \\ {\dot{S}}_2\left(t\right)={\lambda }_{2}I\left(t\right)S_2\left(t\right)-\alpha S_2\left(t\right)-{\delta }_2{S}_2\left(t\right)R\left(t\right)-q{S}_2\left(t\right),\hfill \\ \dot{R}\left(t\right)={\delta }_1{S}_1\left(t\right)R\left(t\right)+{\delta }_2{S}_2\left(t\right)R\left(t\right)-qR\left(t\right),\hfill \end{array}$$

(1)

where $I\left(t\right),S_1\left(t\right),S_2\left(t\right),R\left(t\right)$ are respectively the density of ignorants, the low rate of active spreaders, the high rate of active spreaders and stiflers at time t, $p,q,\alpha ,{\delta }_1,{\delta }_2,{\lambda }_1,{\lambda }_2$ are positive parameters.

Please note that the sense of the parameters $p,q,\alpha ,{\delta }_1,{\delta }_2,{\lambda }_1,{\lambda }_2$ that are used in the rumor spreading model (1) are described in [26]. We will consider the system (1) as a mathematical object and show how stability of nonlinear mathematical models of the similar type can be investigated under influence of stochastic perturbations. In particular, we will consider here the simple parameters ${\lambda }_i$ and ${\delta }_i$ unlike from [26], where these parameters are considered in the form of the product of two parameters: ${\lambda }_i\overline{k}$ and ${\delta }_i\overline{k}$, $i=1,2$. We will not suppose in the general case as it is made in [26] that $p=q$ and ${\delta }_1={\delta }_2$. We will correct also some errors and inaccuracies which there are in [26]. For example, in [26] it is supposed that ${\lambda }_2>{\lambda }_1$ (p. 856) but in the numerical examples the following values are used: ${\lambda }_1=0.05$ and ${\lambda }_2=0.007$ or ${\lambda }_2=0.003$ (p. 862), all equilibria and stability conditions are obtained under the assumption ${\delta }_1={\delta }_2=\delta$ (p. 857) but in the numerical examples one can see ${\delta }_1=0.007$ and ${\delta }_2=0.59$ (p. 862) or ${\delta }_2=0.009$ (p. 863) and so on.

The purpose of the proposed research is to calculate of equilibria of the system (1) and to obtain stability conditions for each from these equilibria under assumption that the system is exposed to stochastic perturbations. Sufficient conditions of stability in probability for each from the five equilibria of the considered model are obtained by virtue of the Routh–Hurwitz criterion [30] and the method of linear matrix inequalities (LMIs) [31,32]. The proposed research method can be used for a lot of other similar nonlinear models in different applications.

2. Equilibria of the Model

Equilibria $E=(I^{*},S_1^{*},S_2^{*},R^{*})$ of the model (1) are defined by the system of algebraic equations

$$\begin{array}{c}({\lambda }_1{S}_1+{\lambda }_2{S}_2+q)I=p,\hfill \\ ({\delta }_{1}R-{\lambda }_{1}I+q)S_1=\alpha S_2,\hfill \\ ({\delta }_{2}R-{\lambda }_{2}I+\alpha +q)S_2=0,\hfill \\ ({\delta }_1{S}_1+{\delta }_2{S}_2-q)R=0,\hfill \end{array}$$

(2)

that follows from (1) by the condition that $I\left(t\right),S_1\left(t\right),S_2\left(t\right),R\left(t\right)$ are constants.

Please note that the solution of the system (2) is not unique. Solving the system (2) gives the following five equilibria $E_i=(I_i^{*},S_{1i}^{*},S_{2i}^{*},R_i^{*})$, $i=0,\dots ,4$, where (see Appendix A.1)

$$\begin{array}{ll}{E}_0=& (I_0^{*},0,0,0),I_0^{*}={\displaystyle \frac{p}{q}};\\ E_1=& (I_1^{*},S_{11}^{*},0,0),I_1^{*}={\displaystyle \frac{q}{{\lambda }_1}},S_{11}^{*}={\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}};\\ E_2=& (I_2^{*},S_{12}^{*},0,R_2^{*}),I_2^{*}={\displaystyle \frac{p{\delta }_1}{q({\delta }_1+{\lambda }_1)}},S_{12}^{*}={\displaystyle \frac{q}{{\delta }_1}},R_2^{*}={\displaystyle \frac{p{\lambda }_1}{q({\delta }_1+{\lambda }_1)}}-{\displaystyle \frac{q}{{\delta }_1}};\\ E_3=& (I_3^{*},S_{13}^{*},S_{23}^{*},0),\\ & I_3^{*}={\displaystyle \frac{\alpha +q}{{\lambda }_2}},S_{13}^{*}={\displaystyle \frac{\alpha (p{\lambda }_2-q(\alpha +q))}{q({\lambda }_2-{\lambda }_1)(\alpha +q)}},S_{23}^{*}={\displaystyle \frac{(q({\lambda }_2-{\lambda }_1)-\alpha {\lambda }_1)(p{\lambda }_2-q(\alpha +q))}{{\lambda }_{2}q({\lambda }_2-{\lambda }_1)(\alpha +q)}};\\ E_4=& (I_4^{*},S_{14}^{*},S_{24}^{*},R_4^{*}),\\ & \mathrm{if}({\delta }_2-{\delta }_1)({\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2)\ne 0\mathrm{then}{S}_{14}^{*}\mathrm{is}\mathrm{a}\mathrm{positive}\mathrm{root}\mathrm{of}\mathrm{the}\mathrm{quadratic}\mathrm{equation}\\ & S_1^2-{\nu }_1{S}_1+{\nu }_2=0,{\nu }_1={\displaystyle \frac{q\alpha +p{\delta }_2}{q({\delta }_2-{\delta }_1)}}+{\displaystyle \frac{q({\lambda }_2+{\delta }_2)}{{\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2}},{\nu }_2={\displaystyle \frac{\alpha q({\lambda }_2+{\delta }_2)}{({\delta }_2-{\delta }_1)({\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2)}},\\ & S_{14}^{*}=\left\{\begin{array}{ll}{\displaystyle \frac{\alpha q^2(\delta +{\lambda }_2)}{\delta ({\lambda }_2-{\lambda }_1)(q\alpha +p\delta )}}& \mathrm{if}{\delta }_2={\delta }_1=\delta ,{\lambda }_2>{\lambda }_1,\\ {\displaystyle \frac{\alpha }{{\delta }_2-{\delta }_1}}& \mathrm{if}{\lambda }_2{\delta }_1={\lambda }_1{\delta }_2,{\delta }_2>{\delta }_1,\end{array}\right.\\ & S_{24}^{*}={\displaystyle \frac{1}{{\delta }_2}}(q-{\delta }_1{S}_{14}^{*}),I_4^{*}={\displaystyle \frac{p}{{\lambda }_1{S}_{14}^{*}+{\lambda }_2{S}_{24}^{*}+q}},R_4^{*}={\displaystyle \frac{{\lambda }_2{I}_4^{*}-\alpha -q}{{\delta }_2}},S_{14}^{*}<{\displaystyle \frac{q}{{\delta }_1}},I_4^{*}>{\displaystyle \frac{\alpha +q}{{\lambda }_2}}.\end{array}$$

(3)

It is supposed that all nonzero elements of all equilibria are positive.

Putting $N\left(t\right)=I\left(t\right)+S_1\left(T\right)+S_2\left(t\right)+R\left(t\right)$ and summing all equations of the system (1), we obtain

$$\begin{array}{c}\dot{N}\left(t\right)=p-qN\left(t\right),N\left(t\right)=\left(N\left(0\right)-{\displaystyle \frac{p}{q}}\right)e^{-qt}+{\displaystyle \frac{p}{q}},\underset{t\to \infty }{lim}N\left(t\right)={\displaystyle \frac{p}{q}}.\hfill \end{array}$$

(4)

In accordance with (4) for all equilibria we have

$$\begin{array}{c}{N}^{*}=I_i^{*}+S_{1i}^{*}+S_{2i}^{*}+R_i^{*}={\displaystyle \frac{p}{q}},i=0,\dots ,4.\hfill \end{array}$$

(5)

3. Stochastic Perturbations, Centralization, and Linearization

Let us suppose that the system (1) is exposed to stochastic perturbations which are directly proportional to the deviation of the system (1) state ($I\left(t\right),S_1\left(t\right),S_2\left(t\right),R\left(t\right)$) from the equilibrium ($I^{*},S_1^{*},S_2^{*},R^{*}$) and are of the type of white noise (${\dot{w}}_0\left(t\right),{\dot{w}}_1\left(t\right),{\dot{w}}_2\left(t\right),{\dot{w}}_3\left(t\right)$), where ($w_0\left(t\right),w_1\left(t\right),w_2\left(t\right),w_3\left(t\right)$) are mutually independent standard Wiener processes. Therefore, we obtain the following system of the Ito stochastic differential equations [33]

$$\begin{array}{c}\dot{I}\left(t\right)=p-{\lambda }_{1}I\left(t\right)S_1\left(t\right)-{\lambda }_{2}I\left(t\right)S_2\left(t\right)-qI\left(t\right)+{\sigma }_0(I\left(t\right)-I^{*}){\dot{w}}_0\left(t\right),\hfill \\ {\dot{S}}_1\left(t\right)={\lambda }_{1}I\left(t\right)S_1\left(t\right)+\alpha S_2\left(t\right)-{\delta }_1{S}_1\left(t\right)R\left(t\right)-q{S}_1\left(t\right)+{\sigma }_1(S_1\left(t\right)-S_1^{*}){\dot{w}}_1\left(t\right),\hfill \\ {\dot{S}}_2\left(t\right)={\lambda }_{2}I\left(t\right)S_2\left(t\right)-\alpha S_2\left(t\right)-{\delta }_2{S}_2\left(t\right)R\left(t\right)-q{S}_2\left(t\right)+{\sigma }_2(S_2\left(t\right)-S_2^{*}){\dot{w}}_2\left(t\right),\hfill \\ \dot{R}\left(t\right)={\delta }_1{S}_1\left(t\right)R\left(t\right)+{\delta }_2{S}_2\left(t\right)R\left(t\right)-qR\left(t\right)+{\sigma }_3(R\left(t\right)-R^{*}){\dot{w}}_3\left(t\right).\hfill \end{array}$$

(6)

Please note that the equilibrium ($I^{*},S_1^{*},S_2^{*},R^{*}$) of the deterministic system (1) is also a solution of the system with stochastic perturbations (6).

Let ($I^{*},S_1^{*},S_2^{*},R^{*}$) be one of the equilibria of the system (1). Putting in (6) $I\left(t\right)=y_0\left(t\right)+I^{*}$, $S_1\left(t\right)=y_1\left(t\right)+S_1^{*}$, $S_2\left(t\right)=y_2\left(t\right)+S_2^{*}$, $R\left(t\right)=y_3\left(t\right)+R^{*}$, we obtain

$$\begin{array}{c}{\dot{y}}_0\left(t\right)=p-(y_0\left(t\right)+I^{*})[{\lambda }_1(y_1\left(t\right)+S_1^{*})+{\lambda }_2(y_2\left(t\right)+S_2^{*})+q]+{\sigma }_0{y}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{y}}_1\left(t\right)=(y_1\left(t\right)+S_1^{*})[{\lambda }_1(y_0\left(t\right)+I^{*})-{\delta }_1(y_3\left(t\right)+R^{*})-q]+\alpha (y_2\left(t\right)+S_2^{*})+{\sigma }_1{y}_1\left(t\right){\dot{w}}_1\left(t\right),\hfill \\ {\dot{y}}_2\left(t\right)=(y_2\left(t\right)+S_2^{*})[{\lambda }_2(y_0\left(t\right)+I^{*})-{\delta }_2(y_3\left(t\right)+R^{*})-\alpha -q]+{\sigma }_2{y}_2\left(t\right){\dot{w}}_2\left(t\right),\hfill \\ {\dot{y}}_3\left(t\right)=(y_3\left(t\right)+R^{*})[{\delta }_1(y_1\left(t\right)+S_1^{*})+{\delta }_2(y_2\left(t\right)+S_2^{*})-q]+{\sigma }_3{y}_3\left(t\right){\dot{w}}_3\left(t\right).\hfill \end{array}$$

(7)

It is clear that stability of the zero solution of the system (7) is equivalent to stability of the equilibrium ($I^{*},S_1^{*},S_2^{*},R^{*}$) of the system (6).

Removing from the system (7) nonlinear terms and using the system for equilibria (2) we obtain the linear part of the system (7)

$$\begin{array}{c}{\dot{z}}_0\left(t\right)=-p{\left(I^{*}\right)}^{-1}{z}_0\left(t\right)-{\lambda }_1{I}^{*}{z}_1\left(t\right)-{\lambda }_2{I}^{*}{z}_2\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)={\lambda }_1{S}_1^{*}{z}_0\left(t\right)-(q+{\delta }_1{R}^{*}-{\lambda }_1{I}^{*})z_1\left(t\right)+\alpha z_2\left(t\right)-{\delta }_1{S}_1^{*}{z}_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\hfill \\ {\dot{z}}_2\left(t\right)={\lambda }_2{S}_2^{*}{z}_0\left(t\right)-(\alpha +q-{\lambda }_2{I}^{*}+{\delta }_2{R}^{*})z_2\left(t\right)-{\delta }_2{S}_2^{*}{z}_3\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\hfill \\ {\dot{z}}_3\left(t\right)={\delta }_1{R}^{*}{z}_1\left(t\right)+{\delta }_2{R}^{*}{z}_2\left(t\right)-(q-{\delta }_1{S}_1^{*}-{\delta }_2{S}_2^{*})z_3\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right).\hfill \end{array}$$

(8)

Let us present the system (8) in the matrix form

$$\begin{array}{c}dz\left(t\right)=Az\left(t\right)dt+C\left(z\right(t\left)\right)dw\left(t\right),\hfill \end{array}$$

(9)

where $z\left(t\right)={(z_0\left(t\right),z_1\left(t\right),z_2\left(t\right),z_3\left(t\right))}^{\prime }$, $w\left(t\right)={(w_0\left(t\right),w_1\left(t\right),w_2\left(t\right),w_3\left(t\right))}^{\prime }$, $C\left(z\left(t\right)\right)=diag({\sigma }_0{z}_0\left(t\right),\dots ,{\sigma }_3{z}_3\left(t\right))$,

$$A=\left[\begin{array}{cccc}-p{\left(I^{*}\right)}^{-1}& -{\lambda }_1{I}^{*}& -{\lambda }_2{I}^{*}& 0\\ {\lambda }_1{S}_1^{*}\hfill & -(q+{\delta }_1{R}^{*}-{\lambda }_1{I}^{*})& \alpha & -{\delta }_1{S}_1^{*}\\ {\lambda }_2{S}_2^{*}\hfill & 0& -(\alpha +q-{\lambda }_2{I}^{*}+{\delta }_2{R}^{*})& -{\delta }_2{S}_2^{*}\\ 0\hfill & {\delta }_1{R}^{*}& {\delta }_2{R}^{*}& -(q-{\delta }_1{S}_1^{*}-{\delta }_2{S}_2^{*})\end{array}\right].$$

(10)

Remark 1.

The order of nonlinearity of the nonlinear system (7) is higher than one. For systems of such type a sufficient condition for asymptotic mean square stability of the zero solution of its linear part (9) provides stability in probability of the zero solution of the initial nonlinear system (7) [30]. Therefore, a sufficient condition for asymptotic mean square stability of the zero solution of the linear Equation (9) provides stability in probability of the equilibrium ($I^{*},S_1^{*},S_2^{*},R^{*}$) of the initial system (6).

Following Remark 1, below we will have sufficient conditions for asymptotic mean square stability of the zero solution of the linear Equation (9) for each from the equilibria (3).

4. Stability of the Equilibria

Consider some definitions and statements that will be used below [30].

Definition 1.

The zero solution of the system (7) is called stable in probability if for any ${\epsilon }_1>0$ and ${\epsilon }_2>0$ there exists $\delta >0$ such that the solution $y\left(t\right)={(y_0\left(t\right),y_1\left(t\right),y_2\left(t\right),y_3\left(t\right))}^{\prime }$ of the system (7) satisfies the condition $\mathbf{P}\{{sup}_{t\ge 0}\left|y\left(t\right)\right|>{\epsilon }_1\}<{\epsilon }_2$ provided that $\mathbf{P}\left\{\right|y\left(0\right)|<\delta \}=1$.

Definition 2.

The zero solution of the system (9) is called:

-

mean square stable if for each $\epsilon >0$ there exists a $\delta >0$ such that ${\mathbf{E}\left|z\left(t\right)\right|}^2<\epsilon$, $t\ge 0$, provided that ${\mathbf{E}\left|z\left(0\right)\right|}^2<\delta$;

-

asymptotically mean square stable if it is mean square stable and the solution $z\left(t\right)$ of Equation (9) satisfies the condition ${\lim }_{t\to \infty }\mathbf{E}{\left|z\left(t\right)\right|}^2=0$ provided that ${\mathbf{E}\left|z\left(0\right)\right|}^2<\infty$.

The generator L of the Ito stochastic differential Equation (9) is defined on the functions $V(t,z)$ which have one continuous derivative with respect to t ($V_t$), two continuous derivatives ($\nabla V$ and ${\nabla }^{2}V$) with respect to z and has the form [30,33]

$$\begin{array}{c}LV(t,z\left(t\right))=V_t(t,z\left(t\right))+\nabla V^{\prime }(t,z\left(t\right))Az\left(t\right)+{\displaystyle \frac{1}{2}}Tr\left[C\left(z\left(t\right)\right){\nabla }^{2}V(t,z\left(t\right))C\left(z\left(t\right)\right)\right].\hfill \end{array}$$

(11)

Theorem 1.

Let there exist a function $V(t,z)$ with continuous derivatives $V_t$, $\nabla V$, ${\nabla }^{2}V$, positive constants $c_1$, $c_2$, $c_3$, such that the following conditions hold:

$$\begin{array}{c}\mathbf{E}V(t,z\left(t\right))\ge c_1{\mathbf{E}\left|z\left(t\right)\right|}^2,\mathbf{E}V(0,z\left(0\right))\le c_2{\mathbf{E}\left|z\left(0\right)\right|}^2,\mathbf{E}LV(t,z\left(t\right))\le -c_3\mathbf{E}{\left|z\left(t\right)\right|}^2.\hfill \end{array}$$

Then the zero solution of Equation (9) is asymptotically mean square stable.

Lemma 1.

Let there exist a positive definite matrix $P=\parallel p_{ij}\parallel$ ($i,j=1,2,3,4$) such that the matrix (10) with the equilibrium ($I^{*},S_1^{*},S_2^{*},R^{*}$) satisfies the linear matrix inequality (LMI)

$$\begin{array}{c}PA+A^{\prime }P+P_{\sigma }<0,P_{\sigma }=diag\{p_{11}{\sigma }_0^2,\dots ,p_{44}{\sigma }_3^2\}.\hfill \end{array}$$

(12)

Then the equilibrium ($I^{*},S_1^{*},S_2^{*},R^{*}$) of the system (6) is stable in probability.

Proof. 

For the function $V(t,z)=z^{\prime }Pz$ from (11) and LMI (12) for some $c>0$ we have

$$\begin{array}{ll}LV(t,z(t\left)\right)& =2z^{\prime }\left(t\right)PAz\left(t\right)+Tr\left[C\left(z\left(t\right)\right)PC\left(z\left(t\right)\right)\right]\\ & =z^{\prime }\left(t\right)(PA+A^{\prime }P+P_{\sigma })z\left(t\right)\le -c{\left|z\left(t\right)\right|}^2.\end{array}$$

From Theorem 1 it follows that the zero solution of Equation (9) is asymptotically mean square stable. Via Remark 1 one can conclude that the equilibrium ($I^{*},S_1^{*},S_2^{*},R^{*}$) of the system (6) is stable in probability. The proof is completed. □

Note to satisfy the LMI (12) the matrix A must be the Hurwitz matrix [30,31].

Definition 3.

The trace of the $k-$th order of a $n\times n$-matrix $A=\parallel a_{ij}\parallel$ is defined as follows:

$$\begin{array}{c}{T}_k={\displaystyle \sum _{1\le i_1<\dots <i_k\le n}}\left|\begin{array}{ccc}{a}_{i_1{i}_1}\hfill & \dots \hfill & a_{i_1{i}_k}\hfill \\ \dots \hfill & \dots \hfill & \dots \hfill \\ a_{i_k{i}_1}\hfill & \dots \hfill & a_{i_k{i}_k}\hfill \end{array}\right|,k=1,\dots ,n.\hfill \end{array}$$

Here, in particular, $T_1=Tr\left(A\right)$, $T_n=det\left(A\right)$, $T_{n-1}={\displaystyle \sum _{i=1}^n}{A}_{ii}$, where $A_{ii}$ is the algebraic complement of the diagonal element $a_{ii}$ of the matrix A.

Lemma 2.

[30,31] Let $T_k$, $k=1,2,3,4$, be the trace of the k-th order of a $4\times 4$-matrix A. The matrix A is the Hurwitz matrix if and only if

$$T_1<0,T_1{T}_2<T_3<0,0<T_1^2{T}_4<(T_1{T}_2-T_3)T_3.$$

(13)

A $3\times 3$-matrix A is the Hurwitz matrix if and only if first two conditions (13) hold.

In general, the LMI (12) for each equilibrium (3) must be numerically investigated via MATLAB. However, in some particular cases this process can be simplified and analytical conditions can be obtained. Below it is shown in investigation of stability of the equilibria (3).

4.1. Stability of the Equilibrium $E_0=(\frac{p}{q},0,0,0)$

Theorem 2.

If

$$\begin{array}{c}{\displaystyle \frac{1}{{\lambda }_1}}>{\displaystyle \frac{p}{q^2}},{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right)>{\displaystyle \frac{p}{q^2}},\hfill \end{array}$$

(14)

and

$$\begin{array}{c}{\sigma }_0^2<2q,{\sigma }_1^2<2\left(q-{\lambda }_1{\displaystyle \frac{p}{q}}\right),{\sigma }_2^2<2\left(\alpha +q-{\lambda }_2{\displaystyle \frac{p}{q}}\right),{\sigma }_3^2<2q,\hfill \end{array}$$

(15)

then the equilibrium $E_0$ is stable in probability.

Proof. 

For the equilibrium $E_0=\left({\displaystyle \frac{p}{q}},0,0,0\right)$ the system (8) takes the form

$$\begin{array}{c}{\dot{z}}_0\left(t\right)=-q{z}_0\left(t\right)-{\lambda }_{1}p{q}^{-1}{z}_1\left(t\right)-{\lambda }_{2}p{q}^{-1}{z}_2\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)=-(q-{\lambda }_{1}p{q}^{-1})z_1\left(t\right)+\alpha z_2\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\hfill \\ {\dot{z}}_2\left(t\right)=-(\alpha +q-{\lambda }_{2}p{q}^{-1})z_2\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\hfill \\ {\dot{z}}_3\left(t\right)=-q{z}_3\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right).\hfill \end{array}$$

(16)

The conditions (14) provide negativity of the coefficients before $z_1\left(t\right)$ and $z_2\left(t\right)$ in the second and the third equations (16). It is known [30] that the last two inequalities (15) are the necessary and sufficient conditions for asymptotic mean square stability of the zero solutions of the last two equations in (16) which do not depend on $z_0\left(t\right)$ and $z_1\left(t\right)$ and can be considered separately. Since ${\lim }_{t\to \infty }\mathbf{E}{z}_2^2\left(t\right)=0$ then the system of first two Equation (16) for $z_0\left(t\right)$ and $z_1\left(t\right)$ can be considered without the process $z_2\left(t\right)$, i.e.,

$$\begin{array}{c}{\dot{z}}_0\left(t\right)=-q{z}_0\left(t\right)-{\lambda }_{1}p{q}^{-1}{z}_1\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)=-(q-{\lambda }_{1}p{q}^{-1})z_1\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right).\hfill \end{array}$$

(17)

Via Remark A2 (see Appendix A.2) the first two inequalities (15) are sufficient for asymptotic mean square stability of the zero solution of the system (17). Therefore, the conditions (14), (15) provide asymptotic mean square stability of the zero solution of the system (16) and via Remark 1 stability in probability of the equilibrium $E_0$ of the system (6). The proof is completed. □

Remark 2.

One can check that by the conditions (14) and (15) the matrix

$$A=\left[\begin{array}{cccc}-q\hfill & -{\lambda }_{1}p{q}^{-1}& -{\lambda }_{2}p{q}^{-1}& 0\\ 0& -(q-{\lambda }_{1}p{q}^{-1})& \alpha & 0\\ 0& 0& -(\alpha +q-{\lambda }_{2}p{q}^{-1})& 0\\ 0& 0& 0& -q\end{array}\right]$$

(18)

of the system (16) satisfies the conditions (13).

Example 1.

Put

$$\begin{array}{c}\alpha =0.4,{\lambda }_1=0.5,{\lambda }_2=0.7,{\delta }_1={\delta }_2=0.2,p=0.8,q=0.7,\hfill \\ {\sigma }_0=1.18,{\sigma }_1=0.50,{\sigma }_2=0.77,{\sigma }_3=1.18.\hfill \end{array}$$

(19)

By these values of the parameters the conditions (14) and (15) hold:

$$\begin{array}{c}{\displaystyle \frac{1}{{\lambda }_1}}=2>{\displaystyle \frac{p}{q^2}}=1.63,{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right)=2.24>{\displaystyle \frac{p}{q^2}}=1.63,\hfill \\ {\sigma }_0^2=1.3924<2q=1.4,{\sigma }_1^2=0.25<2\left(q-{\lambda }_1{\displaystyle \frac{p}{q}}\right)=0.257,\hfill \\ {\sigma }_2^2=0.5929<2\left(\alpha +q-{\lambda }_2{\displaystyle \frac{p}{q}}\right)=0.6,{\sigma }_3^2=1.3924<2q=1.4.\hfill \end{array}$$

Using MATLAB it was shown that by the values of the parameters (19) the matrix (18) satisfies the LMI (12). The conditions (13) with

$$\begin{array}{c}{T}_1=-1.8286<0,T_2=1.1286>0,T_3=-0.2640<0,T_4=0.0189>0,\hfill \\ T_3-T_1{T}_2=1.7997>0,(T_1{T}_2-T_3)T_3-T_1^2{T}_4=0.4119>0,\hfill \end{array}$$

hold too. Therefore, the equilibrium $E_0$ is stable in probability.

In Figure 1 one can see 30 trajectories of the system (6) solution for the equilibrium $E_0$ with the initial condition $I\left(0\right)=1.7$, $S_1\left(0\right)=0.9$, $S_2\left(0\right)=0.7$, $R\left(0\right)=0.5$: all trajectories $I\left(t\right)$ (yellow), $S_1\left(t\right)$ (green), $S_2\left(t\right)$ (blue), $R\left(t\right)$ (red) converge to the equilibrium $E_0=(I^{*},S_1^{*},S_2^{*},R^{*})=(1.1429,0,0,0)$.

4.2. Stability of the Equilibrium $E_1=\left(\frac{q}{{\lambda }_1},\frac{p}{q}-\frac{q}{{\lambda }_1},0,0\right)$

Theorem 3.

If

$$\begin{array}{c}{\displaystyle \frac{1}{{\lambda }_1}}+{\displaystyle \frac{1}{{\delta }_1}}>{\displaystyle \frac{p}{q^2}}>{\displaystyle \frac{1}{{\lambda }_1}},{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right)>{\displaystyle \frac{1}{{\lambda }_1}},\hfill \end{array}$$

(20)

and

$$\begin{array}{c}{\sigma }_0^2<{\displaystyle \frac{2p{\lambda }_1}{q}},{\sigma }_1^2<{\displaystyle \frac{2p{\lambda }_1}{q+{\left[\frac{q}{p{\lambda }_1}\left(1-\frac{q^2}{p{\lambda }_1}\right)\right]}^{-1}}},{\sigma }_2^2<2\left(\alpha +q-q{\displaystyle \frac{{\lambda }_2}{{\lambda }_1}}\right),{\sigma }_3^2<2\left(q-{\delta }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)\right),\hfill \end{array}$$

(21)

then the equilibrium $E_1$ is stable in probability.

Proof. 

For the equilibrium $E_1=\left({\displaystyle \frac{q}{{\lambda }_1}},{\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}},0,0\right)$ the system (8) takes the form

$$\begin{array}{c}{\dot{z}}_0\left(t\right)=-p{q}^{-1}{\lambda }_1{z}_0\left(t\right)-q{z}_1\left(t\right)-q{\lambda }_2{\lambda }_1^{-1}{z}_2\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)={\lambda }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)z_0\left(t\right)+\alpha z_2\left(t\right)-{\delta }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)z_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\hfill \\ {\dot{z}}_2\left(t\right)=-(\alpha +q-q{\lambda }_2{\lambda }_1^{-1})z_2\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\hfill \\ {\dot{z}}_3\left(t\right)=-\left(q-{\delta }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)\right)z_3\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right).\hfill \end{array}$$

(22)

The conditions (20) provide positivity of the nonzero component of the equilibrium $E_1$ and negativity of the coefficients before $z_2\left(t\right)$ and $z_3\left(t\right)$ in the last two equations (22). The last two inequalities (21) are the necessary and sufficient conditions for asymptotic mean square stability of the zero solutions of last two equations in (22) [30] which do not depend on $z_0\left(t\right)$ and $z_1\left(t\right)$ and can be considered separately. Since ${\lim }_{t\to \infty }\mathbf{E}{z}_2^2\left(t\right)=0$ and ${\lim }_{t\to \infty }\mathbf{E}{z}_3^2\left(t\right)=0$ then the system of first two Equation (22) for $z_0\left(t\right)$ and $z_1\left(t\right)$ can be considered without the processes $z_2\left(t\right)$, $z_3\left(t\right)$, i.e.,

$$\begin{array}{c}{\dot{z}}_0\left(t\right)=-p{q}^{-1}{\lambda }_1{z}_0\left(t\right)-q{z}_1\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)={\lambda }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)z_0\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right).\hfill \end{array}$$

(23)

Via Remark A2 (see Appendix A.2) first two inequalities (21) are sufficient for asymptotic mean square stability of the zero solution of the system (23). Therefore, the conditions (20) and (21) provide asymptotic mean square stability of the zero solution of the system (22) and via Remark 1 stability in probability of the equilibrium $E_1$ of the system (6). The proof is completed. □

Remark 3.

One can check that by the conditions (20) and (21) the matrix

$$A=\left[\begin{array}{cccc}-p{q}^{-1}{\lambda }_1& -q& -q{\lambda }_2{\lambda }_1^{-1}& 0\\ {\lambda }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)& 0& \alpha & -{\delta }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)\\ 0& 0& -(\alpha +q-q{\lambda }_2{\lambda }_1^{-1})& 0\\ 0& 0& 0& -\left(q-{\delta }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)\right)\end{array}\right]$$

(24)

of the system (22) satisfies the conditions (13).

Example 2.

Put

$$\begin{array}{c}\alpha =0.4,{\lambda }_1=0.65,{\lambda }_2=0.75,{\delta }_1={\delta }_2=0.2,p=0.9,q=0.7,\hfill \\ {\sigma }_0=1.01,{\sigma }_1=0.41,{\sigma }_2=0.76,{\sigma }_3=1.14.\hfill \end{array}$$

(25)

By these values of the parameters the conditions (20) and (21) hold:

$$\begin{array}{l}{\displaystyle \frac{1}{{\lambda }_1}}+{\displaystyle \frac{1}{{\delta }_1}}=6.538>{\displaystyle \frac{p}{q^2}}=1.837>{\displaystyle \frac{1}{{\lambda }_1}}=1.538,{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right)=2.095>{\displaystyle \frac{1}{{\lambda }_1}}=1.538,\\ {\sigma }_0^2=1.0201<{\displaystyle \frac{2p{\lambda }_1}{q}}<1.671,{\sigma }_1^2=0.1681<{\displaystyle \frac{2p{\lambda }_1}{q+{\left[\frac{q}{p{\lambda }_1}\left(1-\frac{q^2}{p{\lambda }_1}\right)\right]}^{-1}}}=0.2001,\\ {\sigma }_2^2=0.5776<2\left(\alpha +q-q{\displaystyle \frac{{\lambda }_2}{{\lambda }_1}}\right)=0.5846,{\sigma }_3^2=1.2996<2\left(q-{\delta }_1\left({\displaystyle \frac{p}{q}}-{\displaystyle \frac{q}{{\lambda }_1}}\right)\right)=1.3165.\end{array}$$

Using MATLAB it was shown that by the values of the parameters (25) the matrix (24) satisfies the LMI (12), the conditions (13) with

$$\begin{array}{c}{T}_1=-1.7863<0,T_2=1.0818>0,T_3=-0.2511<0,T_4=0.0183>0,\hfill \\ T_3-T_1{T}_2=1.6813>0,(T_1{T}_2-T_3)T_3-T_1^2{T}_4=0.3638>0,\hfill \end{array}$$

hold too. Therefore, the equilibrium $E_1$ is stable in probability.

In Figure 2 one can see 30 trajectories of the system (6) solution for the equilibrium $E_1$ with the initial condition $I\left(0\right)=1.7$, $S_1\left(0\right)=0.9$, $S_2\left(0\right)=0.7$, $R\left(0\right)=0.5$: all trajectories $I\left(t\right)$ (yellow), $S_1\left(t\right)$ (green), $S_2\left(t\right)$ (blue), $R\left(t\right)$ (red) converge to the equilibrium $E_1=(I^{*},S_1^{*},S_2^{*},R^{*})=(1.0769,0.2088,0,0)$.

4.3. Stability of the Equilibrium $E_2=(I_2^{*},S_{12}^{*},0,R_2^{*})$

For the equilibrium $E_2$ the system (8) takes the form

$$\begin{array}{c}{\dot{z}}_0\left(t\right)=-q(1+{\lambda }_1{\delta }_1^{-1})z_0\left(t\right)-{\lambda }_1{I}_2^{*}{z}_1\left(t\right)-{\lambda }_2{I}_2^{*}{z}_2\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)=q{\lambda }_1{\delta }_1^{-1}{z}_0\left(t\right)-(q+{\delta }_1{R}_2^{*}-{\lambda }_1{I}_2^{*})z_1\left(t\right)+\alpha z_2\left(t\right)-q{z}_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\hfill \\ {\dot{z}}_2\left(t\right)=-(\alpha +q-{\lambda }_2{I}_2^{*}+{\delta }_2{R}_2^{*})z_2\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\hfill \\ {\dot{z}}_3\left(t\right)={\delta }_1{R}_2^{*}{z}_1\left(t\right)+{\delta }_2{R}_2^{*}{z}_2\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right),\hfill \end{array}$$

(26)

where $I_2^{*}$ and $R_2^{*}$ are defined in (3).

Lemma 3.

If

$$\begin{array}{c}{\displaystyle \frac{p}{q^2}}>{\displaystyle \frac{1}{{\lambda }_1}}+{\displaystyle \frac{1}{{\delta }_1}},1+{\displaystyle \frac{\alpha }{q}}>{\displaystyle \frac{p({\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2)}{q^2({\delta }_1+{\lambda }_1)}}+{\displaystyle \frac{{\delta }_2}{{\delta }_1}},\hfill \end{array}$$

(27)

then the matrix

$$A=\left[\begin{array}{cccc}-q(1+{\lambda }_1{\delta }_1^{-1})\hfill & -{\lambda }_1{I}_2^{*}& -{\lambda }_2{I}_2^{*}& 0\\ q{\lambda }_1{\delta }_1^{-1}\hfill & 0& \alpha & -q\\ 0\hfill & 0& -(\alpha +q-{\lambda }_2{I}_2^{*}+{\delta }_2{R}_2^{*})& 0\\ 0\hfill & {\delta }_1{R}_2^{*}& {\delta }_2{R}_2^{*}& 0\end{array}\right]$$

(28)

of the system (26) is the Hurwitz matrix.

Proof. 

The first and the second conditions (27) provide respectively a positivity of $R_2^{*}$ and a negativity of the coefficient before $z_2\left(t\right)$ in the third equation of the system (26). Please note that the inequality

$$\begin{array}{c}{\sigma }_2^2<2q\left(1+{\displaystyle \frac{\alpha }{q}}-{\displaystyle \frac{p({\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2)}{q^2({\delta }_1+{\lambda }_1)}}-{\displaystyle \frac{{\delta }_2}{{\delta }_1}}\right)\hfill \end{array}$$

(29)

is the necessary and sufficient condition for asymptotic mean square stability of the zero solution of the equation for $z_2\left(t\right)$ of the system (26). Therefore, by this condition ${\lim }_{t\to \infty }\mathbf{E}{z}_2^2\left(t\right)=0$ it is enough to show that the matrix

$$A=\left[\begin{array}{ccc}-q(1+{\lambda }_1{\delta }_1^{-1})\hfill & -{\lambda }_1{I}_2^{*}& 0\\ q{\lambda }_1{\delta }_1^{-1}\hfill & 0& -q\\ 0\hfill & {\delta }_1{R}_2^{*}& 0\end{array}\right]$$

(30)

is the Hurwitz matrix. Really, for the matrix (30) we have

$$\begin{array}{c}{T}_1=-q(1+{\lambda }_1{\delta }_1^{-1})<0,T_2=q{\lambda }_1^2{\delta }_1^{-1}{I}_2^{*}+q{\delta }_1{R}_2^{*}>0,T_3=-q^2{\delta }_1{R}_2^{*}(1+{\lambda }_1{\delta }_1^{-1})<0,\hfill \end{array}$$

and

$$\begin{array}{ll}{T}_1{T}_2& =(-q(1+{\lambda }_1{\delta }_1^{-1}))(q{\lambda }_1^2{\delta }_1^{-1}{I}_2^{*}+q{\delta }_1{R}_2^{*})\\ & <-q^2{\delta }_1{R}_2^{*}(1+{\lambda }_1{\delta }_1^{-1})=T_3.\end{array}$$

Therefore, the matrix (30) is the Hurwitz matrix. Therefore the matrix (28) is the Hurwitz matrix too. The proof is completed. □

Corollary 1.

If the conditions (27) and (29) hold then for small enough ${\sigma }_0^2$, ${\sigma }_1^2$, ${\sigma }_3^2$ the LMI (12) holds. It means that the zero solution of the linear system (26) is asymptotically mean square stable and therefore (Remark 1) the equilibrium $E_2$ is stable in probability.

Example 3.

Put

$$\begin{array}{c}\alpha =0.4,{\lambda }_1=1,{\lambda }_2=1.3,{\delta }_1=0.5,{\delta }_2=0.7,p=0.9,q=0.5,\hfill \\ {\sigma }_0=0.55,{\sigma }_1=0.30,{\sigma }_2=0.72,{\sigma }_3=0.44.\hfill \end{array}$$

(31)

By these values of the parameters the conditions (27) and (29) hold:

$$\begin{array}{c}{\displaystyle \frac{p}{q^2}}=3.6>{\displaystyle \frac{1}{{\lambda }_1}}+{\displaystyle \frac{1}{{\delta }_1}}=3,1+{\displaystyle \frac{\alpha }{q}}=1.8>{\displaystyle \frac{p({\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2)}{q^2({\delta }_1+{\lambda }_1)}}+{\displaystyle \frac{{\delta }_2}{{\delta }_1}}=1.28,\hfill \\ {\sigma }_2^2=0.5184<2q\left(1+{\displaystyle \frac{\alpha }{q}}-{\displaystyle \frac{p({\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2)}{q^2({\delta }_1+{\lambda }_1)}}-{\displaystyle \frac{{\delta }_2}{{\delta }_1}}\right)=0.52.\hfill \end{array}$$

Using MATLAB it was shown that by the values of the parameters (31) the matrix (28) satisfies the LMI (12), via Lemma 3 the conditions (13) hold too. Therefore, the equilibrium $E_2$ is stable in probability.

In Figure 3 one can see 30 trajectories of the system (6) solution for the equilibrium $E_2$ with the initial condition $I\left(0\right)=1.7$, $S_1\left(0\right)=0.9$, $S_2\left(0\right)=0.7$, $R\left(0\right)=0.5$: all trajectories $I\left(t\right)$ (yellow), $S_1\left(t\right)$ (green), $S_2\left(t\right)$ (blue), $R\left(t\right)$ (red) converge to the equilibrium $E_2=(I^{*},S_1^{*},S_2^{*},R^{*})=(0.6,1,0,0.2)$. In accordance with (5) $I^{*}+S_1^{*}+S_2^{*}+R^{*}=p{q}^{-1}=1.8$.

4.4. Stability of the Equilibrium $E_3=(I_3^{*},S_{13}^{*},S_{23}^{*},0)$

For the equilibrium $E_3$ the system (8) takes the form

$$\begin{array}{l}{\dot{z}}_0\left(t\right)=-p{\lambda }_2{(\alpha +q)}^{-1}{z}_0\left(t\right)-{\lambda }_1{\lambda }_2^{-1}(\alpha +q)z_1\left(t\right)-(\alpha +q)z_2\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)={\lambda }_1{S}_{13}^{*}{z}_0\left(t\right)-(q-{\lambda }_1{\lambda }_2^{-1}(\alpha +q))z_1\left(t\right)+\alpha z_2\left(t\right)-{\delta }_1{S}_{13}^{*}{z}_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\hfill \\ {\dot{z}}_2\left(t\right)={\lambda }_2{S}_{23}^{*}{z}_0\left(t\right)-{\delta }_2{S}_{23}^{*}{z}_3\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\hfill \\ {\dot{z}}_3\left(t\right)=-(q-{\delta }_1{S}_{13}^{*}-{\delta }_2{S}_{23}^{*})z_3\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right),\hfill \end{array}$$

(32)

where $S_{13}^{*}$, $S_{23}^{*}$ are defined in (3).

Lemma 4.

If

$$\begin{array}{c}{\displaystyle \frac{p}{q^2}}>{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right),{\displaystyle \frac{1}{{\lambda }_1}}>{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right),q>{\delta }_1{S}_{13}^{*}+{\delta }_2{S}_{23}^{*},\hfill \end{array}$$

(33)

then the matrix

$$A=\left[\begin{array}{cccc}-p{\lambda }_2{(\alpha +q)}^{-1}\hfill & -{\lambda }_1{\lambda }_2^{-1}(\alpha +q)& -(\alpha +q)& 0\\ {\lambda }_1{S}_{13}^{*}\hfill & -(q-{\lambda }_1{\lambda }_2^{-1}(\alpha +q))& \alpha & -{\delta }_1{S}_{13}^{*}\\ {\lambda }_2{S}_{23}^{*}\hfill & 0& 0& -{\delta }_2{S}_{23}^{*}\\ 0\hfill & 0& 0& -(q-{\delta }_1{S}_{13}^{*}-{\delta }_2{S}_{23}^{*})\end{array}\right]$$

(34)

of the system (26) is the Hurwitz matrix.

Proof. 

The conditions (33) provide a positivity of $S_{13}^{*}$ and $S_{23}^{*}$ and a negativity of the diagonal elements of the matrix (34). Please note that the inequality

$$\begin{array}{c}{\sigma }_3^2<2(q-{\delta }_1{S}_{13}^{*}-{\delta }_2{S}_{23}^{*})\hfill \end{array}$$

(35)

is the necessary and sufficient condition for asymptotic mean square stability of the zero solution of the equation for $z_3\left(t\right)$ of the system (32). Therefore, by this condition ${\lim }_{t\to \infty }\mathbf{E}{z}_3^2\left(t\right)=0$ and it is enough to show that the matrix

$$A=\left[\begin{array}{ccc}-p{\lambda }_2{(\alpha +q)}^{-1}\hfill & -{\lambda }_1{\lambda }_2^{-1}(\alpha +q)& -(\alpha +q)\\ {\lambda }_1{S}_{13}^{*}\hfill & -(q-{\lambda }_1{\lambda }_2^{-1}(\alpha +q))& \alpha \\ {\lambda }_2{S}_{23}^{*}\hfill & 0& 0\end{array}\right]$$

(36)

with

$$\begin{array}{c}{T}_1=-p{\lambda }_2{(\alpha +q)}^{-1}-(q-{\lambda }_1{\lambda }_2^{-1}(\alpha +q))<0,\hfill \\ T_2=p{(\alpha +q)}^{-1}({\lambda }_{2}q-{\lambda }_1(\alpha +q))+(\alpha +q)[{\lambda }_1^2{\lambda }_2^{-1}{S}_{13}^{*}+{\lambda }_2{S}_{23}^{*}]>0,\hfill \\ T_3=-q(\alpha +q)({\lambda }_2-{\lambda }_1)S_{23}^{*}<0,\hfill \end{array}$$

is the Hurwitz matrix, i.e., $T_1{T}_2<T_3$. The proof is completed. □

Corollary 2.

If the conditions (33) and (35) hold then for small enough ${\sigma }_0^2$, ${\sigma }_1^2$, ${\sigma }_2^2$ the LMI (12) holds. It means that the zero solution of the linear system (32) is asymptotically mean square stable and therefore (Remark 1) the equilibrium $E_3$ is stable in probability.

Example 4.

Put

$$\begin{array}{c}\alpha =0.8,{\lambda }_1=0.3,{\lambda }_2=0.9,{\delta }_1=0.8,{\delta }_2=0.7,p=1.2,q=0.6,\hfill \\ {\sigma }_0=0.91,{\sigma }_1=0.50,{\sigma }_2=0.40,{\sigma }_3=0.70.\hfill \end{array}$$

(37)

By these values of the parameters the conditions (33) and (35) hold:

$$\begin{array}{c}{\displaystyle \frac{p}{q^2}}=3.33>{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right)=2.59,{\displaystyle \frac{1}{{\lambda }_1}}=3.33>{\displaystyle \frac{1}{{\lambda }_2}}\left(1+{\displaystyle \frac{\alpha }{q}}\right)=2.59,\hfill \\ q=0.6>{\delta }_1{S}_{13}^{*}+{\delta }_2{S}_{23}^{*}=0.1715,{\sigma }_3^2=0.49<2(q-{\delta }_1{S}_{13}^{*}-{\delta }_2{S}_{23}^{*})=0.5015.\hfill \end{array}$$

Using MATLAB it was shown that by the values of the parameters (37) the matrix (34) satisfies the LMI (12), for the matrix (36) $T_1=-0.9048<0$, $T_2=0.2362>0$, $T_3=-0.0320<0$, $T_3-T_1{T}_2=0.1817>0$, the conditions (13) hold too. Therefore, the equilibrium $E_3$ is stable in probability.

In Figure 4 one can see 30 trajectories of the system (6) solution for the equilibrium $E_3$ with the initial condition $I\left(0\right)=1.9$, $S_1\left(0\right)=0.8$, $S_2\left(0\right)=0.4$, $R\left(0\right)=0.4$: all trajectories $I\left(t\right)$ (yellow), $S_1\left(t\right)$ (green), $S_2\left(t\right)$ (blue), $R\left(t\right)$ (red) converge to the equilibrium $E_3=(I^{*},S_1^{*},S_2^{*},R^{*})=(1.5556,0.3810,0.0634,0)$. In accordance with (5) $I^{*}+S_1^{*}+S_2^{*}+R^{*}=p{q}^{-1}=2$.

4.5. Stability of the Equilibrium $E_4=(I_4^{*},S_{14}^{*},S_{24}^{*},R_4^{*})$

For the equilibrium $E_4$ the system (8) by virtue of (2) takes the form

$$\begin{array}{c}{\dot{z}}_0\left(t\right)=-p{\left(I_4^{*}\right)}^{-1}{z}_0\left(t\right)-{\lambda }_1{I}_4^{*}{z}_1\left(t\right)-{\lambda }_2{I}_4^{*}{z}_2\left(t\right)+{\sigma }_0{z}_0\left(t\right){\dot{w}}_0\left(t\right),\hfill \\ {\dot{z}}_1\left(t\right)={\lambda }_1{S}_{14}^{*}{z}_0\left(t\right)-\alpha S_{24}^{*}{\left(S_{14}^{*}\right)}^{-1}{z}_1\left(t\right)+\alpha z_2\left(t\right)-{\delta }_1{S}_{14}^{*}{z}_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\hfill \\ {\dot{z}}_2\left(t\right)={\lambda }_2{S}_{24}^{*}{z}_0\left(t\right)-{\delta }_2{S}_{24}^{*}{z}_3\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\hfill \\ {\dot{z}}_3\left(t\right)={\delta }_1{R}_4^{*}{z}_1\left(t\right)+{\delta }_2{R}_4^{*}{z}_2\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right),\hfill \end{array}$$

(38)

where $I_4^{*}$, $S_{14}^{*}$, $S_{24}^{*}$, $R_4^{*}$ are defined in (3).

Let us show that the matrix

$$A=\left[\begin{array}{cccc}-p{\left(I_4^{*}\right)}^{-1}\hfill & -{\lambda }_1{I}_4^{*}& -{\lambda }_2{I}_4^{*}& 0\\ {\lambda }_1{S}_{14}^{*}\hfill & -\alpha S_{24}^{*}{\left(S_{14}^{*}\right)}^{-1}& \alpha & -{\delta }_1{S}_{14}^{*}\\ {\lambda }_2{S}_{24}^{*}\hfill & 0& 0& -{\delta }_2{S}_{24}^{*}\\ 0\hfill & {\delta }_1{R}_4^{*}& {\delta }_2{R}_4^{*}& 0\end{array}\right]$$

(39)

of the system (38) is the Hurwitz matrix. Really, the conditions (13) for the matrix (39) hold with

$$\begin{array}{l}{T}_1=-p{\left(I_4^{*}\right)}^{-1}-\alpha S_{24}^{*}{\left(S_{14}^{*}\right)}^{-1}<0,\hfill \\ T_2=\alpha p{S}_{24}^{*}{\left(I_4^{*}{S}_{14}^{*}\right)}^{-1}+I_4^{*}({\lambda }_1^2{S}_{14}^{*}+{\lambda }_2^2{S}_{24}^{*})+R_4^{*}({\delta }_1^2{S}_{14}^{*}+{\delta }_2^2{S}_{24}^{*})>0,\\ T_3=-\alpha S_{24}^{*}({\lambda }_1{\lambda }_2{I}_4^{*}+{\delta }_1{\delta }_2{R}_4^{*})-p{\left(I_4^{*}\right)}^{-1}{R}_4^{*}({\delta }_1^2{S}_{14}^{*}+{\delta }_2^2{S}_{24}^{*})-\alpha {\left(S_{24}^{*}\right)}^2{\left(S_{14}^{*}\right)}^{-1}({\lambda }_2^2{I}_4^{*}+{\delta }_2^2{R}_4^{*})<0,\\ T_4=\alpha p{R}_4^{*}{\left(I_4^{*}{S}_{14}^{*}\right)}^{-1}({\delta }_1^2{\left(S_{14}^{*}\right)}^2+{\delta }_2^2{\left(S_{24}^{*}\right)}^2)+{({\lambda }_2{\delta }_1-{\lambda }_1{\delta }_2)}^2{I}_4^{*}{S}_{14}^{*}{S}_{24}^{*}{R}_4^{*}>0.\end{array}$$

Example 5.

Put

$$\begin{array}{c}\alpha =0.4,{\lambda }_1=0.3,{\lambda }_2=0.9,{\delta }_1={\delta }_2=0.8,p=1.2,q=0.6,\hfill \\ {\sigma }_0=0.51,{\sigma }_1=0.51,{\sigma }_2=0.55,{\sigma }_3=0.34.\hfill \end{array}$$

(40)

Using MATLAB it was shown that by the values of the parameters (40) the matrix (39) satisfies the LMI (12), the conditions (13) hold too: $T_1=-1.3259<0$, $T_2=0.7020>0$, $T_3=-0.1828<0$, $T_4=0.0187>0$, $T_3-T_1{T}_2=0.7479>0$, $(T_1{T}_2-T_3)T_3-T_1^2{T}_4=0.1039>0$. Therefore, the equilibrium $E_4$ is stable in probability.

In Figure 5 one can see 30 trajectories of the system (6) solution for the equilibrium $E_4$ with the initial condition $I\left(0\right)=1.9$, $S_1\left(0\right)=0.8$, $S_2\left(0\right)=0.7$, $R\left(0\right)=0.4$: all trajectories $I\left(t\right)$ (yellow), $S_1\left(t\right)$ (green), $S_2\left(t\right)$ (blue), $R\left(t\right)$ (red) converge to the equilibrium $E_4=(I^{*},S_1^{*},S_2^{*},R^{*})=(1.1765,0.4250,0.3250,0.0735)$. In accordance with (5) $I^{*}+S_1^{*}+S_2^{*}+R^{*}=p{q}^{-1}=2$.

Please note that decreasing ${\delta }_2$ from ${\delta }_2=0.8$ to ${\delta }_2=0.7$, we obtain that $S_1^{*}$ unlike from the previous case is calculated via quadratic equation (see (3)). By that with the same values of all other parameters the equilibrium $E_4$ a bit changed $E_4=(I^{*},S_1^{*},S_2^{*},R^{*})=(1.1534,0.4543,0.3379,0.0544)$, $I^{*}+S_1^{*}+S_2^{*}+R^{*}=p{q}^{-1}=2$, but remains stable in probability and the conditions (13) hold with $T_1=-1.3379<0$, $T_2=0.6971>0$, $T_3=-0.1686<0$, $T_4=0.0143>0$, $T_3-T_1{T}_2=0.7641>0$, $(T_1{T}_2-T_3)T_3-T_1^2{T}_4=0.1033>0$.

5. Conclusions

In this paper, it is shown how the dynamics of the very popular I2SR rumor spreading model can be investigated under stochastic perturbations. It is shown that for some equilibria of the considered model it is possible to get conditions for stability in probability in an analytical form, for other equilibria stability condition can be obtained numerically by an appropriate linear matrix inequality via MATLAB.

The proposed way of research can be used for more detail investigation both the considered I2SR rumor spreading model and also all other known type of rumor spreading models [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].

Besides, this research method can be used for a detailed investigation of many other nonlinear mathematical models (with the order of nonlinearity higher than one) in different other applications. In particular, the proposed method can be used for systems with exponential nonlinearity [34,35], together with stochastic perturbations of the type of white noise other types of perturbations can be used, for instance, perturbations of the type of Poisson’s jumps [35], the method does not depend on the dimension of the considered system and can be used for systems of more than four equations.