*Article* **Stability of Equilibria of Rumor Spreading Model under Stochastic Perturbations**

#### **Leonid Shaikhet**

Department of Mathematics, Ariel University, Ariel 40700, Israel; leonid.shaikhet@usa.net

Received: 2 February 2020; Accepted: 11 February 2020; Published: 18 February 2020

**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.

**Keywords:** rumor spreading model; white noise; stochastic differential equations; asymptotic mean square stability; stability in probability; linear matrix inequality

#### **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–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–29]). Following [26], we will consider the rumor spreading model (the so-called I2SR-model) in the form

$$\begin{aligned} I(t) &= p - \lambda\_1 I(t) S\_1(t) - \lambda\_2 I(t) S\_2(t) - q I(t), \\ \dot{S}\_1(t) &= \lambda\_1 I(t) S\_1(t) + a S\_2(t) - \delta\_1 S\_1(t) R(t) - q S\_1(t), \\ \dot{S}\_2(t) &= \lambda\_2 I(t) S\_2(t) - a S\_2(t) - \delta\_2 S\_2(t) R(t) - q S\_2(t), \\ \dot{R}(t) &= \delta\_1 S\_1(t) R(t) + \delta\_2 S\_2(t) R(t) - q R(t), \end{aligned} \tag{1}$$

where *I*(*t*), *S*1(*t*), *S*2(*t*), *R*(*t*) 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*, *α*, *δ*1, *δ*2, *λ*1, *λ*<sup>2</sup> are positive parameters.

Please note that the sense of the parameters *p*, *q*, *α*, *δ*1, *δ*2, *λ*1, *λ*<sup>2</sup> 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 *λ<sup>i</sup>* and *δ<sup>i</sup>* unlike from [26], where these parameters are considered in the form of the product of two parameters: *λik* and *δik*, *i* = 1, 2. We will not suppose in the general case as it is made in [26] that *p* = *q* and *δ*<sup>1</sup> = *δ*2. We will correct also some errors and inaccuracies which there are in [26]. For example, in [26] it is supposed that *λ*<sup>2</sup> > *λ*<sup>1</sup> (p. 856) but in the numerical examples the following values are used: *λ*<sup>1</sup> = 0.05 and *λ*<sup>2</sup> = 0.007 or *λ*<sup>2</sup> = 0.003 (p. 862), all equilibria and stability conditions are obtained

under the assumption *δ*<sup>1</sup> = *δ*<sup>2</sup> = *δ* (p. 857) but in the numerical examples one can see *δ*<sup>1</sup> = 0.007 and *δ*<sup>2</sup> = 0.59 (p. 862) or *δ*<sup>2</sup> = 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*∗ <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup>, *R*∗) of the model (1) are defined by the system of algebraic equations

$$\begin{cases} (\lambda\_1 S\_1 + \lambda\_2 S\_2 + q)I = p, \\ (\delta\_1 R - \lambda\_1 I + q)S\_1 = aS\_{2\prime} \\ (\delta\_2 R - \lambda\_2 I + a + q)S\_2 = 0, \\ (\delta\_1 S\_1 + \delta\_2 S\_2 - q)R = 0, \end{cases} \tag{2}$$

that follows from (1) by the condition that *I*(*t*), *S*1(*t*), *S*2(*t*), *R*(*t*) are constants.

Please note that the solution of the system (2) is not unique. Solving the system (2) gives the following five equilibria *Ei* = (*I*<sup>∗</sup> *<sup>i</sup>* , *S*<sup>∗</sup> 1*i* , *S*∗ 2*i* , *R*∗ *<sup>i</sup>* ), *i* = 0, ..., 4, where (see Appendix A.1)

*E*<sup>0</sup> =(*I* ∗ <sup>0</sup> , 0, 0, 0), *I* ∗ <sup>0</sup> <sup>=</sup> *<sup>p</sup> q* ; *E*<sup>1</sup> =(*I* ∗ <sup>1</sup> , *S*<sup>∗</sup> 11, 0, 0), *I* ∗ <sup>1</sup> <sup>=</sup> *<sup>q</sup> λ*1 , *S*∗ <sup>11</sup> <sup>=</sup> *<sup>p</sup> <sup>q</sup>* <sup>−</sup> *<sup>q</sup> λ*1 ; *E*<sup>2</sup> =(*I* ∗ <sup>2</sup> , *S*<sup>∗</sup> 12, 0, *R*<sup>∗</sup> <sup>2</sup> ), *I* ∗ <sup>2</sup> <sup>=</sup> *<sup>p</sup>δ*<sup>1</sup> *q*(*δ*<sup>1</sup> + *λ*1) , *S*∗ <sup>12</sup> <sup>=</sup> *<sup>q</sup> δ*1 , *R*∗ <sup>2</sup> <sup>=</sup> *<sup>p</sup>λ*<sup>1</sup> *<sup>q</sup>*(*δ*<sup>1</sup> <sup>+</sup> *<sup>λ</sup>*1) <sup>−</sup> *<sup>q</sup> δ*1 ; *E*<sup>3</sup> =(*I* ∗ <sup>3</sup> , *S*<sup>∗</sup> 13, *S*<sup>∗</sup> 23, 0), *I* ∗ <sup>3</sup> <sup>=</sup> *<sup>α</sup>* <sup>+</sup> *<sup>q</sup> λ*2 , *S*∗ <sup>13</sup> <sup>=</sup> *<sup>α</sup>*(*pλ*<sup>2</sup> <sup>−</sup> *<sup>q</sup>*(*<sup>α</sup>* <sup>+</sup> *<sup>q</sup>*)) *q*(*λ*<sup>2</sup> − *λ*1)(*α* + *q*) , *S*∗ <sup>23</sup> <sup>=</sup> (*q*(*λ*<sup>2</sup> <sup>−</sup> *<sup>λ</sup>*1) <sup>−</sup> *αλ*1)(*pλ*<sup>2</sup> <sup>−</sup> *<sup>q</sup>*(*<sup>α</sup>* <sup>+</sup> *<sup>q</sup>*)) *<sup>λ</sup>*2*q*(*λ*<sup>2</sup> <sup>−</sup> *<sup>λ</sup>*1)(*<sup>α</sup>* <sup>+</sup> *<sup>q</sup>*) ; *E*<sup>4</sup> =(*I* ∗ <sup>4</sup> , *S*<sup>∗</sup> 14, *S*<sup>∗</sup> 24, *R*<sup>∗</sup> 4 ), if (*δ*<sup>2</sup> − *δ*1)(*λ*2*δ*<sup>1</sup> − *λ*1*δ*2) = 0 then *S*<sup>∗</sup> <sup>14</sup> is a positive root of the quadratic equation *S*2 <sup>1</sup> <sup>−</sup> *<sup>ν</sup>*1*S*<sup>1</sup> <sup>+</sup> *<sup>ν</sup>*<sup>2</sup> <sup>=</sup> 0, *<sup>ν</sup>*<sup>1</sup> <sup>=</sup> *<sup>q</sup><sup>α</sup>* <sup>+</sup> *<sup>p</sup>δ*<sup>2</sup> *<sup>q</sup>*(*δ*<sup>2</sup> <sup>−</sup> *<sup>δ</sup>*1) <sup>+</sup> *<sup>q</sup>*(*λ*<sup>2</sup> <sup>+</sup> *<sup>δ</sup>*2) *λ*2*δ*<sup>1</sup> − *λ*1*δ*<sup>2</sup> , *<sup>ν</sup>*<sup>2</sup> <sup>=</sup> *<sup>α</sup>q*(*λ*<sup>2</sup> <sup>+</sup> *<sup>δ</sup>*2) (*δ*<sup>2</sup> − *δ*1)(*λ*2*δ*<sup>1</sup> − *λ*1*δ*2) , *S*∗ <sup>14</sup> = ⎧ ⎪⎨ ⎪⎩ *αq*2(*δ* + *λ*2) *<sup>δ</sup>*(*λ*<sup>2</sup> <sup>−</sup> *<sup>λ</sup>*1)(*q<sup>α</sup>* <sup>+</sup> *<sup>p</sup>δ*) if *<sup>δ</sup>*<sup>2</sup> <sup>=</sup> *<sup>δ</sup>*<sup>1</sup> <sup>=</sup> *<sup>δ</sup>*, *<sup>λ</sup>*<sup>2</sup> <sup>&</sup>gt; *<sup>λ</sup>*1, *α δ*<sup>2</sup> − *δ*<sup>1</sup> if *λ*2*δ*<sup>1</sup> = *λ*1*δ*2, *δ*<sup>2</sup> > *δ*1, (3)

$$S\_{24}^{\*} = \frac{1}{\delta\_2} (\eta - \delta\_1 S\_{14}^{\*}), \quad I\_4^{\*} = \frac{p}{\lambda\_1 S\_{14}^{\*} + \lambda\_2 S\_{24}^{\*} + q}, \quad R\_4^{\*} = \frac{\lambda\_2 I\_4^{\*} - \kappa - q}{\delta\_2}, \quad S\_{14}^{\*} < \frac{q}{\delta\_1}, \quad I\_4^{\*} > \frac{\kappa + q}{\lambda\_2}.$$

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

Putting *N*(*t*) = *I*(*t*) + *S*1(*T*) + *S*2(*t*) + *R*(*t*) and summing all equations of the system (1), we obtain

$$\dot{N}(t) = p - qN(t), \qquad N(t) = \left(N(0) - \frac{p}{q}\right)e^{-qt} + \frac{p}{q}, \qquad \lim\_{t \to \infty} N(t) = \frac{p}{q}. \tag{4}$$

In accordance with (4) for all equilibria we have

$$N^\* = I\_i^\* + S\_{1i}^\* + S\_{2i}^\* + R\_i^\* = \frac{p}{q}, \qquad i = 0, \ldots, 4. \tag{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*(*t*), *S*1(*t*), *S*2(*t*), *R*(*t*)) from the equilibrium (*I*∗, *S*∗ <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup>, *R*∗) and are of the type of white noise (*w*˙ <sup>0</sup>(*t*), *w*˙ <sup>1</sup>(*t*), *w*˙ <sup>2</sup>(*t*), *w*˙ <sup>3</sup>(*t*)), where (*w*0(*t*), *w*1(*t*), *w*2(*t*), *w*3(*t*)) are mutually independent standard Wiener processes. Therefore, we obtain the following system of the Ito stochastic differential equations [33]

$$\begin{aligned} \dot{I}(t) &= p - \lambda\_1 I(t) S\_1(t) - \lambda\_2 I(t) S\_2(t) - q I(t) + \sigma\_0 (I(t) - I^\*) \dot{w}\_0(t), \\ S\_1(t) &= \lambda\_1 I(t) S\_1(t) + a S\_2(t) - \delta\_1 S\_1(t) R(t) - q S\_1(t) + \sigma\_1 (S\_1(t) - S\_1^\*) \dot{w}\_1(t), \\ \dot{S}\_2(t) &= \lambda\_2 I(t) S\_2(t) - a S\_2(t) - \delta\_2 S\_2(t) R(t) - q S\_2(t) + \sigma\_2 (S\_2(t) - S\_2^\*) \dot{w}\_2(t), \\ \dot{R}(t) &= \delta\_1 S\_1(t) R(t) + \delta\_2 S\_2(t) R(t) - q R(t) + \sigma\_3 (R(t) - R^\*) \dot{w}\_3(t). \end{aligned} \tag{6}$$

Please note that the equilibrium (*I*∗, *S*∗ <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup> , *R*∗) of the deterministic system (1) is also a solution of the system with stochastic perturbations (6).

Let (*I*∗, *S*∗ <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup> , *R*∗) be one of the equilibria of the system (1). Putting in (6) *I*(*t*) = *y*0(*t*) + *I*∗, *S*1(*t*) = *y*1(*t*) + *S*<sup>∗</sup> <sup>1</sup> , *S*2(*t*) = *y*2(*t*) + *S*<sup>∗</sup> <sup>2</sup> , *R*(*t*) = *y*3(*t*) + *R*∗, we obtain

$$\begin{aligned} \dot{y}\_0(t) &= p - (y\_0(t) + I^\*) [\lambda\_1 (y\_1(t) + S\_1^\*) + \lambda\_2 (y\_2(t) + S\_2^\*) + q] + \sigma\_0 y\_0(t) \dot{w}\_0(t), \\ \dot{y}\_1(t) &= (y\_1(t) + S\_1^\*) [\lambda\_1 (y\_0(t) + I^\*) - \delta\_1 (y\_3(t) + R^\*) - q] + a (y\_2(t) + S\_2^\*) + \sigma\_1 y\_1(t) \dot{w}\_1(t), \\ \dot{y}\_2(t) &= (y\_2(t) + S\_2^\*) [\lambda\_2 (y\_0(t) + I^\*) - \delta\_2 (y\_3(t) + R^\*) - a - q] + \sigma\_2 y\_2(t) \dot{w}\_2(t), \\ \dot{y}\_3(t) &= (y\_3(t) + R^\*) [\delta\_1 (y\_1(t) + S\_1^\*) + \delta\_2 (y\_2(t) + S\_2^\*) - q] + \sigma\_3 y\_3(t) \dot{w}\_3(t). \end{aligned} \tag{7}$$

It is clear that stability of the zero solution of the system (7) is equivalent to stability of the equilibrium (*I*∗, *S*∗ <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup>, *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{aligned} \dot{z}\_{0}(t) &= -p(I^{\*})^{-1}z\_{0}(t) - \lambda\_{1}I^{\*}z\_{1}(t) - \lambda\_{2}I^{\*}z\_{2}(t) + \sigma\_{0}z\_{0}(t)\dot{w}\_{0}(t),\\ \dot{z}\_{1}(t) &= \lambda\_{1}S\_{1}^{\*}z\_{0}(t) - (q + \delta\_{1}R^{\*} - \lambda\_{1}I^{\*})z\_{1}(t) + \alpha z\_{2}(t) - \delta\_{1}S\_{1}^{\*}z\_{3}(t) + \sigma\_{1}z\_{1}(t)\dot{w}\_{1}(t),\\ \dot{z}\_{2}(t) &= \lambda\_{2}S\_{2}^{\*}z\_{0}(t) - (a + q - \lambda\_{2}I^{\*} + \delta\_{2}R^{\*})z\_{2}(t) - \delta\_{2}S\_{2}^{\*}z\_{3}(t) + \sigma\_{2}z\_{2}(t)\dot{w}\_{2}(t),\\ \dot{z}\_{3}(t) &= \delta\_{1}R^{\*}z\_{1}(t) + \delta\_{2}R^{\*}z\_{2}(t) - (q - \delta\_{1}S\_{1}^{\*} - \delta\_{2}S\_{2}^{\*})z\_{3}(t) + \sigma\_{3}z\_{3}(t)\dot{w}\_{3}(t). \end{aligned} \tag{8}$$

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

$$dz(t) = Az(t)dt + \mathbb{C}(z(t))dw(t),\tag{9}$$

where *z*(*t*)=(*z*0(*t*), *z*1(*t*), *z*2(*t*), *z*3(*t*)) , *w*(*t*)=(*w*0(*t*), *w*1(*t*), *w*2(*t*), *w*3(*t*)) , *C*(*z*(*t*)) = *diag*(*σ*0*z*0(*t*), ..., *σ*3*z*3(*t*)),

$$A = \begin{bmatrix} -p(I^\*)^{-1} & -\lambda\_1 I^\* & -\lambda\_2 I^\* & 0\\ \lambda\_1 S\_1^\* & -(q + \delta\_1 R^\* - \lambda\_1 I^\*) & a & -\delta\_1 S\_1^\*\\ \lambda\_2 S\_2^\* & 0 & -(a + q - \lambda\_2 I^\* + \delta\_2 R^\*) & -\delta\_2 S\_2^\*\\ 0 & \delta\_1 R^\* & \delta\_2 R^\* & -(q - \delta\_1 S\_1^\* - \delta\_2 S\_2^\*) \end{bmatrix}. \tag{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*∗ <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup>, *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 ε*<sup>1</sup> > 0 *and ε*<sup>2</sup> > 0 *there exists δ* > 0 *such that the solution y*(*t*)=(*y*0(*t*), *y*1(*t*), *y*2(*t*), *y*3(*t*)) *of the system* (7) *satisfies the condition* **<sup>P</sup>**{sup*t*≥<sup>0</sup> <sup>|</sup>*y*(*t*)<sup>|</sup> <sup>&</sup>gt; *<sup>ε</sup>*1} <sup>&</sup>lt; *<sup>ε</sup>*<sup>2</sup> *provided that* **<sup>P</sup>**{|*y*(0)<sup>|</sup> <sup>&</sup>lt; *<sup>δ</sup>*} <sup>=</sup> <sup>1</sup>*.*

**Definition 2.** *The zero solution of the system* (9) *is called:*


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* (*Vt*), two continuous derivatives (∇*V* and <sup>∇</sup>2*V*) with respect to *<sup>z</sup>* and has the form [30,33]

$$LV(t, z(t)) = V\_l(t, z(t)) + \nabla V'(t, z(t))Az(t) + \frac{1}{2}Tr[\mathbb{C}(z(t))\nabla^2 V(t, z(t))\mathbb{C}(z(t))].\tag{11}$$

**Theorem 1.** *Let there exist a function <sup>V</sup>*(*t*, *<sup>z</sup>*) *with continuous derivatives Vt,* <sup>∇</sup>*V,* <sup>∇</sup>2*V, positive constants c*1*, c*2*, c*3*, such that the following conditions hold:*

$$\mathbb{E}V(t, z(t)) \ge c\_1 \mathbb{E}|z(t)|^2, \quad \mathbb{E}V(0, z(0)) \le c\_2 \mathbb{E}|z(0)|^2, \quad \mathbb{E}LV(t, z(t)) \le -c\_3 \mathbb{E}|z(t)|^2.$$

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

**Lemma 1.** *Let there exist a positive definite matrix P* = *pij (i*, *j* = 1, 2, 3, 4*) such that the matrix* (10) *with the equilibrium (I*∗, *S*∗ <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup> , *R*∗*) satisfies the linear matrix inequality (LMI)*

$$PA + A'P + P\_{\sigma} < 0, \qquad P\_{\sigma} = \text{diag}\{p\_{11}\sigma\_{0}^{2}, \dots, p\_{44}\sigma\_{3}^{2}\}. \tag{12}$$

*Then the equilibrium (I*∗, *S*∗ <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup>, *R*∗*) of the system* (6) *is stable in probability.*

**Proof.** For the function *V*(*t*, *z*) = *z Pz* from (11) and LMI (12) for some *c* > 0 we have

$$\begin{aligned} LV(t, z(t)) &= 2z'(t)PAz(t) + Tr[C(z(t))PC(z(t))] \\ &= z'(t)(PA + A'P + P\_\sigma)z(t) \leq -c|z(t)|^2. \end{aligned}$$

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*∗ <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup>, *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* × *n-matrix A* = *aij is defined as follows:*

$$T\_k = \sum\_{1 \le i\_1 < \dots < i\_k \le n} \left| \begin{array}{ccccc} a\_{i\_1 i\_1} & \dots & a\_{i\_1 i\_k} \\ \dots & \dots & \dots \\ a\_{i\_k i\_1} & \dots & a\_{i\_k i\_k} \end{array} \right|, \qquad k = 1, \dots, n.$$

*Here, in particular, <sup>T</sup>*<sup>1</sup> <sup>=</sup> *Tr*(*A*)*, Tn* <sup>=</sup> det(*A*)*, Tn*−<sup>1</sup> <sup>=</sup> *<sup>n</sup>* ∑ *i*=1 *Aii, where Aii is the algebraic complement of the diagonal element aii of the matrix A.*

**Lemma 2.** *[30,31] Let Tk, k* = 1, 2, 3, 4*, be the trace of the k-th order of a* 4 × 4*-matrix A. The matrix A is the Hurwitz matrix if and only if*

$$T\_1 < 0, \quad T\_1 T\_2 < T\_3 < 0, \quad 0 < T\_1^2 T\_4 < \left(T\_1 T\_2 - T\_3\right) T\_3. \tag{13}$$

*A* 3 × 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*<sup>0</sup> = ( *<sup>p</sup> <sup>q</sup>* , 0, 0, 0)

**Theorem 2.** *If*

$$\frac{1}{\lambda\_1} > \frac{p}{q^{2'}}, \qquad \frac{1}{\lambda\_2} \left(1 + \frac{a}{q}\right) > \frac{p}{q^{2'}} \tag{14}$$

*and*

$$
\sigma\_0^2 < 2\eta, \quad \sigma\_1^2 < 2\left(q - \lambda\_1 \frac{p}{q}\right), \quad \sigma\_2^2 < 2\left(a + q - \lambda\_2 \frac{p}{q}\right), \quad \sigma\_3^2 < 2\eta. \tag{15}
$$

*then the equilibrium E*<sup>0</sup> *is stable in probability.*

$$\textbf{Proof.}\text{ For the equilibrium } E\_0 = \left(\frac{p}{q}, 0, 0, 0\right) \text{ the system (8) takes the form}$$

$$\begin{aligned} \dot{z}\_0(t) &= -qz\_0(t) - \lambda\_1 pq^{-1}z\_1(t) - \lambda\_2 pq^{-1}z\_2(t) + \sigma\_0 z\_0(t)\dot{w}\_0(t), \\ \dot{z}\_1(t) &= -(q - \lambda\_1 pq^{-1})z\_1(t) + az\_2(t) + \sigma\_1 z\_1(t)\dot{w}\_1(t), \\ \dot{z}\_2(t) &= -(\mathfrak{a} + q - \lambda\_2 pq^{-1})z\_2(t) + \sigma\_2 z\_2(t)\dot{w}\_2(t), \\ \dot{z}\_3(t) &= -qz\_3(t) + \sigma\_3 z\_3(t)\dot{w}\_3(t). \end{aligned} \tag{16}$$

The conditions (14) provide negativity of the coefficients before *z*1(*t*) and *z*2(*t*) 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 *<sup>z</sup>*0(*t*) and *<sup>z</sup>*1(*t*) and can be considered separately. Since lim*t*→<sup>∞</sup> **<sup>E</sup>***z*<sup>2</sup> <sup>2</sup>(*t*) = 0 then the system of first two Equation (16) for *z*0(*t*) and *z*1(*t*) can be considered without the process *z*2(*t*), i.e.,

$$\begin{cases} \dot{z}\_0(t) = -qz\_0(t) - \lambda\_1 pq^{-1} z\_1(t) + \sigma\_0 z\_0(t) \dot{w}\_0(t), \\ \dot{z}\_1(t) = -(q - \lambda\_1 pq^{-1}) z\_1(t) + \sigma\_1 z\_1(t) \dot{w}\_1(t). \end{cases} \tag{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*<sup>0</sup> of the system (6). The proof is completed.

**Remark 2.** *One can check that by the conditions* (14) *and* (15) *the matrix*

$$A = \begin{bmatrix} -q & -\lambda\_1 pq^{-1} & -\lambda\_2 pq^{-1} & 0\\ 0 & -(q - \lambda\_1 pq^{-1}) & a & 0\\ 0 & 0 & -(a + q - \lambda\_2 pq^{-1}) & 0\\ 0 & 0 & 0 & -q \end{bmatrix} \tag{18}$$

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

**Example 1.** *Put*

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

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

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

*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{aligned} T\_1 &= -1.8286 < 0, \quad T\_2 = 1.1286 > 0, \quad T\_3 = -0.2640 < 0, \quad T\_4 = 0.0189 > 0, \\ T\_3 - T\_1 T\_2 &= 1.7997 > 0, \quad (T\_1 T\_2 - T\_3) T\_3 - T\_1^2 T\_4 = 0.4119 > 0, \end{aligned}$$

*hold too. Therefore, the equilibrium E*<sup>0</sup> *is stable in probability.*

*In Figure 1 one can see 30 trajectories of the system* (6) *solution for the equilibrium E*<sup>0</sup> *with the initial condition I*(0) = 1.7*, S*1(0) = 0.9*, S*2(0) = 0.7*, R*(0) = 0.5*: all trajectories I*(*t*) *(yellow), S*1(*t*) *(green), S*2(*t*) *(blue), R*(*t*) *(red) converge to the equilibrium E*<sup>0</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup>, *R*∗)=(1.1429, 0, 0, 0)*.*

**Figure 1.** 30 trajectories of the system (6) solution with the initial condition *I*(0) = 1.7, *S*1(0) = 0.9, *S*2(0) = 0.7, *R*(0) = 0.5: all trajectories *I*(*t*) (yellow), *S*1(*t*) (green), *S*2(*t*) (blue), *R*(*t*) (red) converge to the equilibrium *E*<sup>0</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup> , *R*∗)=(1.1429, 0, 0, 0).

*4.2. Stability of the Equilibrium E*<sup>1</sup> = *q λ*1 , *p <sup>q</sup>* <sup>−</sup> *<sup>q</sup> λ*1 , 0, 0 

**Theorem 3.** *If*

$$1 \frac{1}{\lambda\_1} + \frac{1}{\delta\_1} > \frac{p}{q^2} > \frac{1}{\lambda\_1}, \qquad \frac{1}{\lambda\_2} \left( 1 + \frac{a}{q} \right) > \frac{1}{\lambda\_1},\tag{20}$$

*and*

$$
\sigma\_0^2 < \frac{2p\lambda\_1}{q}, \quad \sigma\_1^2 < \frac{2p\lambda\_1}{q + \left[\frac{q}{p\lambda\_1}\left(1 - \frac{q^2}{p\lambda\_1}\right)\right]^{-1}}, \quad \sigma\_2^2 < 2\left(a + q - q\frac{\lambda\_2}{\lambda\_1}\right), \quad \sigma\_3^2 < 2\left(q - \delta\_1\left(\frac{p}{q} - \frac{q}{\lambda\_1}\right)\right), \tag{21}
$$

*then the equilibrium E*<sup>1</sup> *is stable in probability.*

**Proof.** For the equilibrium *E*<sup>1</sup> = *q λ*1 , *p <sup>q</sup>* <sup>−</sup> *<sup>q</sup> λ*1 , 0, 0 the system (8) takes the form

$$\begin{aligned} \dot{z}\_{0}(t) &= -pq^{-1}\lambda\_{1}z\_{0}(t) - qz\_{1}(t) - q\lambda\_{2}\lambda\_{1}^{-1}z\_{2}(t) + \sigma\_{0}z\_{0}(t)\dot{w}\_{0}(t),\\ \dot{z}\_{1}(t) &= \lambda\_{1}\left(\frac{p}{q} - \frac{q}{\lambda\_{1}}\right)z\_{0}(t) + az\_{2}(t) - \delta\_{1}\left(\frac{p}{q} - \frac{q}{\lambda\_{1}}\right)z\_{3}(t) + \sigma\_{1}z\_{1}(t)\dot{w}\_{1}(t),\\ \dot{z}\_{2}(t) &= -\left(a + q - q\lambda\_{2}\lambda\_{1}^{-1}\right)z\_{2}(t) + \sigma\_{2}z\_{2}(t)\dot{w}\_{2}(t),\\ \dot{z}\_{3}(t) &= -\left(q - \delta\_{1}\left(\frac{p}{q} - \frac{q}{\lambda\_{1}}\right)\right)z\_{3}(t) + \sigma\_{3}z\_{3}(t)\dot{w}\_{3}(t). \end{aligned} \tag{22}$$

The conditions (20) provide positivity of the nonzero component of the equilibrium *E*<sup>1</sup> and negativity of the coefficients before *z*2(*t*) and *z*3(*t*) 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(*t*) and *z*1(*t*) and can be considered separately. Since lim*t*→<sup>∞</sup> **<sup>E</sup>***z*<sup>2</sup> <sup>2</sup>(*t*) = 0 and lim*t*→<sup>∞</sup> **<sup>E</sup>***z*<sup>2</sup> <sup>3</sup>(*t*) = 0 then the system of first two Equation (22) for *z*0(*t*) and *z*1(*t*) can be considered without the processes *z*2(*t*), *z*3(*t*), i.e.,

$$\begin{aligned} \dot{z}\_0(t) &= -pq^{-1}\lambda\_1 z\_0(t) - qz\_1(t) + \sigma v z\_0(t)\dot{w}\_0(t),\\ \dot{z}\_1(t) &= \lambda\_1 \left(\frac{p}{q} - \frac{q}{\lambda\_1}\right) z\_0(t) + \sigma\_1 z\_1(t)\dot{w}\_1(t). \end{aligned} \tag{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*<sup>1</sup> of the system (6). The proof is completed.

**Remark 3.** *One can check that by the conditions* (20) *and* (21) *the matrix*

$$A = \begin{bmatrix} -pq^{-1}\lambda\_1 & -q & -q\lambda\_2\lambda\_1^{-1} & 0\\ \lambda\_1\left(\frac{p}{q} - \frac{q}{\lambda\_1}\right) & 0 & a & -\delta\_1\left(\frac{p}{q} - \frac{q}{\lambda\_1}\right) \\ 0 & 0 & -(a + q - q\lambda\_2\lambda\_1^{-1}) & 0 \\ 0 & 0 & 0 & -\left(q - \delta\_1\left(\frac{p}{q} - \frac{q}{\lambda\_1}\right)\right) \end{bmatrix} \tag{24}$$

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

**Example 2.** *Put*

$$\begin{aligned} \alpha &= 0.4, \quad \lambda\_1 = 0.65, \quad \lambda\_2 = 0.75, \quad \delta\_1 = \delta\_2 = 0.2, \quad p = 0.9, \quad q = 0.7, \\ \sigma\_0 &= 1.01, \quad \sigma\_1 = 0.41, \quad \sigma\_2 = 0.76, \quad \sigma\_3 = 1.14. \end{aligned} \tag{25}$$

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

$$\begin{aligned} \frac{1}{\lambda\_1} + \frac{1}{\delta\_1} &= 6.538 > \frac{p}{q^2} = 1.837 > \frac{1}{\lambda\_1} = 1.538, \qquad \frac{1}{\lambda\_2} \left( 1 + \frac{a}{q} \right) = 2.095 > \frac{1}{\lambda\_1} = 1.538, \\\sigma\_0^2 &= 1.0201 < \frac{2p\lambda\_1}{q} < 1.671, \quad \sigma\_1^2 = 0.1681 < \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( a + q - q \frac{\lambda\_2}{\lambda\_1} \right) = 0.5846, \quad \sigma\_3^2 = 1.2996 < 2 \left( q - \delta\_1 \left( \frac{p}{q} - \frac{q}{\lambda\_1} \right) \right) = 1.3165. \end{aligned}$$

*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{aligned} T\_1 &= -1.7863 < 0, \quad T\_2 = 1.0818 > 0, \quad T\_3 = -0.2511 < 0, \quad T\_4 = 0.0183 > 0, \\ T\_3 - T\_1 T\_2 &= 1.6813 > 0, \quad (T\_1 T\_2 - T\_3) T\_3 - T\_1^2 T\_4 = 0.3638 > 0, \end{aligned}$$

*hold too. Therefore, the equilibrium E*<sup>1</sup> *is stable in probability.*

*In Figure 2 one can see 30 trajectories of the system* (6) *solution for the equilibrium E*<sup>1</sup> *with the initial condition I*(0) = 1.7*, S*1(0) = 0.9*, S*2(0) = 0.7*, R*(0) = 0.5*: all trajectories I*(*t*) *(yellow), S*1(*t*) *(green), S*2(*t*) *(blue), R*(*t*) *(red) converge to the equilibrium E*<sup>1</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup>, *R*∗)=(1.0769, 0.2088, 0, 0)*.*

**Figure 2.** 30 trajectories of the system (6) solution with the initial condition *I*(0) = 1.7, *S*1(0) = 0.9, *S*2(0) = 0.7, *R*(0) = 0.5: all trajectories *I*(*t*) (yellow), *S*1(*t*) (green), *S*2(*t*) (blue), *R*(*t*) (red) converge to the equilibrium *E*<sup>1</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup> , *R*∗)=(1.0769, 0.2088, 0, 0).

*4.3. Stability of the Equilibrium E*<sup>2</sup> = (*I*<sup>∗</sup> <sup>2</sup> , *S*<sup>∗</sup> <sup>12</sup>, 0, *R*<sup>∗</sup> 2 )

For the equilibrium *E*<sup>2</sup> the system (8) takes the form

$$\begin{aligned} \dot{z}\_{0}(t) &= -q(1 + \lambda\_{1}\delta\_{1}^{-1})z\_{0}(t) - \lambda\_{1}l\_{2}^{\*}z\_{1}(t) - \lambda\_{2}l\_{2}^{\*}z\_{2}(t) + \sigma\_{0}z\_{0}(t)\dot{w}\_{0}(t),\\ \dot{z}\_{1}(t) &= q\lambda\_{1}\delta\_{1}^{-1}z\_{0}(t) - (q + \delta\_{1}R\_{2}^{\*} - \lambda\_{1}l\_{2}^{\*})z\_{1}(t) + az\_{2}(t) - qz\_{3}(t) + \sigma\_{1}z\_{1}(t)\dot{w}\_{1}(t),\\ \dot{z}\_{2}(t) &= -(a + q - \lambda\_{2}l\_{2}^{\*} + \delta\_{2}R\_{2}^{\*})z\_{2}(t) + \sigma\_{2}z\_{2}(t)\dot{w}\_{2}(t),\\ \dot{z}\_{3}(t) &= \delta\_{1}R\_{2}^{\*}z\_{1}(t) + \delta\_{2}R\_{2}^{\*}z\_{2}(t) + \sigma\_{3}z\_{3}(t)\dot{w}\_{3}(t), \end{aligned} \tag{26}$$

where *I*∗ <sup>2</sup> and *R*<sup>∗</sup> <sup>2</sup> are defined in (3).

**Lemma 3.** *If*

$$\frac{p}{q^2} > \frac{1}{\lambda\_1} + \frac{1}{\delta\_1}, \qquad 1 + \frac{a}{q} > \frac{p(\lambda\_2 \delta\_1 - \lambda\_1 \delta\_2)}{q^2(\delta\_1 + \lambda\_1)} + \frac{\delta\_2}{\delta\_1}, \tag{27}$$

*then the matrix*

$$A = \begin{bmatrix} -q(1 + \lambda\_1 \delta\_1^{-1}) & -\lambda\_1 I\_2^\* & -\lambda\_2 I\_2^\* & 0\\ q\lambda\_1 \delta\_1^{-1} & 0 & a & -q\\ 0 & 0 & -(a + q - \lambda\_2 I\_2^\* + \delta\_2 R\_2^\*) & 0\\ 0 & \delta\_1 R\_2^\* & \delta\_2 R\_2^\* & 0 \end{bmatrix} \tag{28}$$

*of the system* (26) *is the Hurwitz matrix.*

**Proof.** The first and the second conditions (27) provide respectively a positivity of *R*∗ <sup>2</sup> and a negativity of the coefficient before *z*2(*t*) in the third equation of the system (26). Please note that the inequality

$$
\sigma\_2^2 < 2q \left( 1 + \frac{a}{q} - \frac{p(\lambda\_2 \delta\_1 - \lambda\_1 \delta\_2)}{q^2(\delta\_1 + \lambda\_1)} - \frac{\delta\_2}{\delta\_1} \right) \tag{29}
$$

is the necessary and sufficient condition for asymptotic mean square stability of the zero solution of the equation for *<sup>z</sup>*2(*t*) of the system (26). Therefore, by this condition lim*t*→<sup>∞</sup> **<sup>E</sup>***z*<sup>2</sup> <sup>2</sup>(*t*) = 0 it is enough to show that the matrix

$$A = \begin{bmatrix} -q(1 + \lambda\_1 \delta\_1^{-1}) & -\lambda\_1 I\_2^\* & 0\\ q\lambda\_1 \delta\_1^{-1} & 0 & -q\\ 0 & \delta\_1 R\_2^\* & 0 \end{bmatrix} \tag{30}$$

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

$$T\_1 = -q(1 + \lambda\_1 \delta\_1^{-1}) < 0, \quad T\_2 = q\lambda\_1^2 \delta\_1^{-1} I\_2^\* + q\delta\_1 R\_2^\* > 0, \quad T\_3 = -q^2 \delta\_1 R\_2^\* (1 + \lambda\_1 \delta\_1^{-1}) < 0,$$

and

$$\begin{aligned} T\_1 T\_2 &= \left( -q \left( 1 + \lambda\_1 \delta\_1^{-1} \right) \right) \left( q \lambda\_1^2 \delta\_1^{-1} I\_2^\* + q \delta\_1 \mathcal{R}\_2^\* \right) \\ &< -q^2 \delta\_1 \mathcal{R}\_2^\* (1 + \lambda\_1 \delta\_1^{-1}) = T\_3. \end{aligned}$$

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 σ*<sup>2</sup> <sup>0</sup> *, <sup>σ</sup>*<sup>2</sup> <sup>1</sup> *, <sup>σ</sup>*<sup>2</sup> <sup>3</sup> *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*<sup>2</sup> *is stable in probability.*

#### **Example 3.** *Put*

$$\begin{aligned} &a = 0.4, \quad \lambda\_1 = 1, \quad \lambda\_2 = 1.3, \quad \delta\_1 = 0.5, \quad \delta\_2 = 0.7, \quad p = 0.9, \quad q = 0.5, \\ &\sigma\_0 = 0.55, \quad \sigma\_1 = 0.30, \quad \sigma\_2 = 0.72, \quad \sigma\_3 = 0.44. \end{aligned} \tag{31}$$

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

$$\begin{array}{ll} \frac{p}{q^2} = 3.6 > \frac{1}{\lambda\_1} + \frac{1}{\delta\_1} = 3, & 1 + \frac{a}{q} = 1.8 > \frac{p(\lambda\_2 \delta\_1 - \lambda\_1 \delta\_2)}{q^2(\delta\_1 + \lambda\_1)} + \frac{\delta\_2}{\delta\_1} = 1.28, \\\sigma\_2^2 = 0.5184 < 2q \left( 1 + \frac{a}{q} - \frac{p(\lambda\_2 \delta\_1 - \lambda\_1 \delta\_2)}{q^2(\delta\_1 + \lambda\_1)} - \frac{\delta\_2}{\delta\_1} \right) = 0.52. \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*<sup>2</sup> *is stable in probability.*

*In Figure 3 one can see 30 trajectories of the system* (6) *solution for the equilibrium E*<sup>2</sup> *with the initial condition I*(0) = 1.7*, S*1(0) = 0.9*, S*2(0) = 0.7*, R*(0) = 0.5*: all trajectories I*(*t*) *(yellow), S*1(*t*) *(green), S*2(*t*) *(blue), R*(*t*) *(red) converge to the equilibrium E*<sup>2</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup> , *R*∗)=(0.6, 1, 0, 0.2)*. In accordance with* (5) *I*∗ + *S*∗ <sup>1</sup> + *S*<sup>∗</sup> <sup>2</sup> + *<sup>R</sup>*<sup>∗</sup> = *pq*−<sup>1</sup> = 1.8*.*

**Figure 3.** 30 trajectories of the system (6) solution with the initial condition *I*(0) = 1.7, *S*1(0) = 0.9, *S*2(0) = 0.7, *R*(0) = 0.5: all trajectories *I*(*t*) (yellow), *S*1(*t*) (green), *S*2(*t*) (blue), *R*(*t*) (red) converge to the equilibrium *E*<sup>2</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup> , *R*∗)=(0.6, 1, 0, 0.2).

#### *4.4. Stability of the Equilibrium E*<sup>3</sup> = (*I*<sup>∗</sup> <sup>3</sup> , *S*<sup>∗</sup> <sup>13</sup>, *S*<sup>∗</sup> <sup>23</sup>, 0)

For the equilibrium *E*<sup>3</sup> the system (8) takes the form

$$\begin{aligned} \dot{z}\_{0}(t) &= -p\lambda\_{2}(a+q)^{-1}z\_{0}(t) - \lambda\_{1}\lambda\_{2}^{-1}(a+q)z\_{1}(t) - (a+q)z\_{2}(t) + \sigma\_{0}z\_{0}(t)\dot{w}\_{0}(t),\\ \dot{z}\_{1}(t) &= \lambda\_{1}S\_{13}^{\*}z\_{0}(t) - (q-\lambda\_{1}\lambda\_{2}^{-1}(a+q))z\_{1}(t) + az\_{2}(t) - \delta\_{1}S\_{13}^{\*}z\_{3}(t) + \sigma\_{1}z\_{1}(t)\dot{w}\_{1}(t),\\ \dot{z}\_{2}(t) &= \lambda\_{2}S\_{23}^{\*}z\_{0}(t) - \delta\_{2}S\_{23}^{\*}z\_{3}(t) + \sigma\_{2}z\_{2}(t)\dot{w}\_{2}(t),\\ \dot{z}\_{3}(t) &= -(q-\delta\_{1}S\_{13}^{\*} - \delta\_{2}S\_{23}^{\*})z\_{3}(t) + \sigma\_{3}z\_{3}(t)\dot{w}\_{3}(t),\end{aligned} \tag{32}$$

where *S*∗ <sup>13</sup>, *S*<sup>∗</sup> <sup>23</sup> are defined in (3).

**Lemma 4.** *If*

$$\frac{p}{q^2} > \frac{1}{\lambda\_2} \left( 1 + \frac{a}{q} \right), \qquad \frac{1}{\lambda\_1} > \frac{1}{\lambda\_2} \left( 1 + \frac{a}{q} \right), \qquad q > \delta\_1 S\_{13}^\* + \delta\_2 S\_{23}^\* \tag{33}$$

*then the matrix*

$$A = \begin{bmatrix} -p\lambda\_2(a+q)^{-1} & -\lambda\_1\lambda\_2^{-1}(a+q) & -(a+q) & 0\\ \lambda\_1 S\_{13}^\* & -(q-\lambda\_1\lambda\_2^{-1}(a+q)) & a & -\delta\_1 S\_{13}^\*\\ \lambda\_2 S\_{23}^\* & 0 & 0 & -\delta\_2 S\_{23}^\*\\ 0 & 0 & 0 & -(q-\delta\_1 S\_{13}^\* - \delta\_2 S\_{23}^\*) \end{bmatrix} \tag{34}$$

*of the system* (26) *is the Hurwitz matrix.*

**Proof.** The conditions (33) provide a positivity of *S*∗ <sup>13</sup> and *S*<sup>∗</sup> <sup>23</sup> and a negativity of the diagonal elements of the matrix (34). Please note that the inequality

$$
\sigma\_3^2 < 2(q - \delta\_1 S\_{13}^\* - \delta\_2 S\_{23}^\*) \tag{35}
$$

is the necessary and sufficient condition for asymptotic mean square stability of the zero solution of the equation for *<sup>z</sup>*3(*t*) of the system (32). Therefore, by this condition lim*t*→<sup>∞</sup> **<sup>E</sup>***z*<sup>2</sup> <sup>3</sup>(*t*) = 0 and it is enough to show that the matrix

$$A = \begin{bmatrix} -p\lambda\_2(a+q)^{-1} & -\lambda\_1\lambda\_2^{-1}(a+q) & -(a+q) \\ \lambda\_1 S\_{13}^\* & -(q-\lambda\_1\lambda\_2^{-1}(a+q)) & a \\ \lambda\_2 S\_{23}^\* & 0 & 0 \end{bmatrix} \tag{36}$$

with

$$\begin{aligned} T\_1 &= -p\lambda\_2(\mathfrak{a}+q)^{-1} - (q - \lambda\_1\lambda\_2^{-1}(\mathfrak{a}+q)) < 0, \\ T\_2 &= p(\mathfrak{a}+q)^{-1}(\lambda\_2q - \lambda\_1(\mathfrak{a}+q)) + (\mathfrak{a}+q)[\lambda\_1^2\lambda\_2^{-1}S\_{13}^\* + \lambda\_2S\_{23}^\*] > 0, \\ T\_3 &= -q(\mathfrak{a}+q)(\lambda\_2 - \lambda\_1)S\_{23}^\* < 0, \end{aligned}$$

is the Hurwitz matrix, i.e., *T*1*T*<sup>2</sup> < *T*3. The proof is completed.

**Corollary 2.** *If the conditions* (33) *and* (35) *hold then for small enough σ*<sup>2</sup> <sup>0</sup> *, <sup>σ</sup>*<sup>2</sup> <sup>1</sup> *, <sup>σ</sup>*<sup>2</sup> <sup>2</sup> *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*<sup>3</sup> *is stable in probability.*

#### **Example 4.** *Put*

$$\begin{aligned} \alpha &= 0.8, \quad \lambda\_1 = 0.3, \quad \lambda\_2 = 0.9, \quad \delta\_1 = 0.8, \quad \delta\_2 = 0.7, \quad p = 1.2, \quad q = 0.6, \\\sigma\_0 &= 0.91, \quad \sigma\_1 = 0.50, \quad \sigma\_2 = 0.40, \quad \sigma\_3 = 0.70. \end{aligned} \tag{37}$$

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

$$\begin{array}{llll} \frac{p}{q^2} = 3.33 > \frac{1}{\lambda\_2} \left( 1 + \frac{a}{q} \right) = 2.59, & \frac{1}{\lambda\_1} = 3.33 > \frac{1}{\lambda\_2} \left( 1 + \frac{a}{q} \right) = 2.59, \\\ 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. \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*<sup>1</sup> = −0.9048 < 0*, T*<sup>2</sup> = 0.2362 > 0*, T*<sup>3</sup> = −0.0320 < 0*, T*<sup>3</sup> − *T*1*T*<sup>2</sup> = 0.1817 > 0*, the conditions* (13) *hold too. Therefore, the equilibrium E*<sup>3</sup> *is stable in probability.*

*In Figure 4 one can see 30 trajectories of the system* (6) *solution for the equilibrium E*<sup>3</sup> *with the initial condition I*(0) = 1.9*, S*1(0) = 0.8*, S*2(0) = 0.4*, R*(0) = 0.4*: all trajectories I*(*t*) *(yellow), S*1(*t*) *(green), S*2(*t*) *(blue), R*(*t*) *(red) converge to the equilibrium E*<sup>3</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup>, *R*∗)=(1.5556, 0.3810, 0.0634, 0)*. In accordance with* (5) *I*∗ + *S*∗ <sup>1</sup> + *S*<sup>∗</sup> <sup>2</sup> + *<sup>R</sup>*<sup>∗</sup> = *pq*−<sup>1</sup> = <sup>2</sup>*.*

**Figure 4.** 30 trajectories of the system (6) solution with the initial condition *I*(0) = 1.9, *S*1(0) = 0.8, *S*2(0) = 0.4, *R*(0) = 0.4: all trajectories *I*(*t*) (yellow), *S*1(*t*) (green), *S*2(*t*) (blue), *R*(*t*) (red) converge to the equilibrium *E*<sup>3</sup> = (*I*∗, *S*<sup>∗</sup> , *S*<sup>∗</sup> , *R*∗)=(1.5556, 0.3810, 0.0634, 0).

*4.5. Stability of the Equilibrium E*<sup>4</sup> = (*I*<sup>∗</sup> , *S*<sup>∗</sup> , *S*<sup>∗</sup> , *R*<sup>∗</sup> )

For the equilibrium *E*<sup>4</sup> the system (8) by virtue of (2) takes the form

$$\begin{aligned} \dot{z}\_{0}(t) &= -p(I\_{4}^{\*})^{-1}z\_{0}(t) - \lambda\_{1}I\_{4}^{\*}z\_{1}(t) - \lambda\_{2}I\_{4}^{\*}z\_{2}(t) + \sigma\_{0}z\_{0}(t)\dot{w}\_{0}(t), \\ \dot{z}\_{1}(t) &= \lambda\_{1}S\_{14}^{\*}z\_{0}(t) - aS\_{24}^{\*}(S\_{14}^{\*})^{-1}z\_{1}(t) + az\_{2}(t) - \delta\_{1}S\_{14}^{\*}z\_{3}(t) + \sigma\_{1}z\_{1}(t)\dot{w}\_{1}(t), \\ \dot{z}\_{2}(t) &= \lambda\_{2}S\_{24}^{\*}z\_{0}(t) - \delta\_{2}S\_{24}^{\*}z\_{3}(t) + \sigma\_{2}z\_{2}(t)\dot{w}\_{2}(t), \\ \dot{z}\_{3}(t) &= \delta\_{1}R\_{4}^{\*}z\_{1}(t) + \delta\_{2}R\_{4}^{\*}z\_{2}(t) + \sigma\_{3}z\_{3}(t)\dot{w}\_{3}(t), \end{aligned} \tag{38}$$

where *I*∗ , *S*<sup>∗</sup> , *S*<sup>∗</sup> , *R*<sup>∗</sup> are defined in (3).

Let us show that the matrix

$$A = \begin{bmatrix} -p(I\_4^\*)^{-1} & -\lambda\_1 I\_4^\* & -\lambda\_2 I\_4^\* & 0\\ \lambda\_1 S\_{14}^\* & -aS\_{24}^\*(S\_{14}^\*)^{-1} & a & -\delta\_1 S\_{14}^\*\\ \lambda\_2 S\_{24}^\* & 0 & 0 & -\delta\_2 S\_{24}^\*\\ 0 & \delta\_1 R\_4^\* & \delta\_2 R\_4^\* & 0 \end{bmatrix} \tag{39}$$

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

$$\begin{cases} T\_1 = -p(I\_4^\*)^{-1} - a S\_{24}^\* (S\_{14}^\*)^{-1} < 0, \\ T\_2 = ap S\_{24}^\* (I\_4^\* S\_{14}^\*)^{-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 = -a S\_{24}^\* (\lambda\_1 \lambda\_2 I\_4^\* + \delta\_1 \delta\_2 R\_4^\*) - p(I\_4^\*)^{-1} R\_4^\* (\delta\_1^2 S\_{14}^\* + \delta\_2^2 S\_{24}^\*) - a (S\_{24}^\*)^2 (S\_{14}^\*)^{-1} (\lambda\_2^2 I\_4^\* + \delta\_2^2 R\_4^\*) < 0, \\ T\_4 = ap R\_4^\* (I\_4^\* S\_{14}^\*)^{-1} (\delta\_1^2 (S\_{14}^\*)^2 + \delta\_2^2 (S\_{24}^\*)^2) + (\lambda\_2 \delta\_1 - \lambda\_1 \delta\_2)^2 I\_4^\* S\_{14}^\* S\_{24}^\* R\_4^\* > 0. \end{cases}$$

**Example 5.** *Put*

$$\begin{aligned} \alpha &= 0.4, \quad \lambda\_1 = 0.3, \quad \lambda\_2 = 0.9, \quad \delta\_1 = \delta\_2 = 0.8, \quad p = 1.2, \quad q = 0.6, \\\sigma\_0 &= 0.51, \quad \sigma\_1 = 0.51, \quad \sigma\_2 = 0.55, \quad \sigma\_3 = 0.34. \end{aligned} \tag{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*<sup>1</sup> = −1.3259 < 0*, T*<sup>2</sup> = 0.7020 > 0*, T*<sup>3</sup> = −0.1828 < 0*, T*<sup>4</sup> = 0.0187 > 0*,*

*<sup>T</sup>*<sup>3</sup> <sup>−</sup> *<sup>T</sup>*1*T*<sup>2</sup> <sup>=</sup> 0.7479 <sup>&</sup>gt; <sup>0</sup>*,* (*T*1*T*<sup>2</sup> <sup>−</sup> *<sup>T</sup>*3)*T*<sup>3</sup> <sup>−</sup> *<sup>T</sup>*<sup>2</sup> <sup>1</sup> *T*<sup>4</sup> = 0.1039 > 0*. Therefore, the equilibrium E*<sup>4</sup> *is stable in probability.*

*In Figure 5 one can see 30 trajectories of the system* (6) *solution for the equilibrium E*<sup>4</sup> *with the initial condition I*(0) = 1.9*, S*1(0) = 0.8*, S*2(0) = 0.7*, R*(0) = 0.4*: all trajectories I*(*t*) *(yellow), S*1(*t*) *(green), S*2(*t*) *(blue), R*(*t*) *(red) converge to the equilibrium E*<sup>4</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup>, *S*<sup>∗</sup> <sup>2</sup>, *R*∗)=(1.1765, 0.4250, 0.3250, 0.0735)*. In accordance with* (5) *I*∗ + *S*∗ <sup>1</sup> + *S*<sup>∗</sup> <sup>2</sup> + *<sup>R</sup>*<sup>∗</sup> = *pq*−<sup>1</sup> = <sup>2</sup>*.*

*Please note that decreasing δ*<sup>2</sup> *from δ*<sup>2</sup> = 0.8 *to δ*<sup>2</sup> = 0.7*, we obtain that S*<sup>∗</sup> <sup>1</sup> *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*<sup>4</sup> *a bit changed E*<sup>4</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup> , *R*∗)=(1.1534, 0.4543, 0.3379, 0.0544)*, I*<sup>∗</sup> + *S*<sup>∗</sup> <sup>1</sup> + *S*<sup>∗</sup> <sup>2</sup> + *<sup>R</sup>*<sup>∗</sup> = *pq*−<sup>1</sup> = 2*, but remains stable in probability and the conditions* (13) *hold with T*<sup>1</sup> = −1.3379 < 0*, T*<sup>2</sup> = 0.6971 > 0*, <sup>T</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>0.1686 <sup>&</sup>lt; <sup>0</sup>*, T*<sup>4</sup> <sup>=</sup> 0.0143 <sup>&</sup>gt; <sup>0</sup>*, T*<sup>3</sup> <sup>−</sup> *<sup>T</sup>*1*T*<sup>2</sup> <sup>=</sup> 0.7641 <sup>&</sup>gt; <sup>0</sup>*,* (*T*1*T*<sup>2</sup> <sup>−</sup> *<sup>T</sup>*3)*T*<sup>3</sup> <sup>−</sup> *<sup>T</sup>*<sup>2</sup> <sup>1</sup> *T*<sup>4</sup> = 0.1033 > 0*.*

**Figure 5.** 30 trajectories of the system (6) solution with the initial condition *I*(0) = 1.9, *S*1(0) = 0.8, *S*2(0) = 0.7, *R*(0) = 0.4: all trajectories *I*(*t*) (yellow), *S*1(*t*) (green), *S*2(*t*) (blue), *R*(*t*) (red) converge to the equilibrium *E*<sup>4</sup> = (*I*∗, *S*<sup>∗</sup> <sup>1</sup> , *S*<sup>∗</sup> <sup>2</sup> , *R*∗)=(1.1765, 0.4250, 0.3250, 0.0735).

#### **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–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.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Appendix A**

#### *Appendix A.1. Equilibria of the System* (1)

The equilibria *E*0, ..., *E*<sup>3</sup> of the system (1) are obtained from the system (2) quite simply (see (3)). To get the equilibrium *E*<sup>4</sup> note that from the second and the third equations of the system (2) we obtain

$$R = \frac{1}{\delta\_1} (\alpha S\_2 S\_1^{-1} - q + \lambda\_1 I) = \frac{1}{\delta\_2} (\lambda\_2 I - (\alpha + q)).$$

From this and the first equation of the system (2) we have

$$I = \frac{(\alpha + q)\delta\_1 + (\alpha S\_2 S\_1^{-1} - q)\delta\_2}{\lambda\_2 \delta\_1 - \lambda\_1 \delta\_2} = \frac{p}{\lambda\_1 S\_1 + \lambda\_2 S\_2 + q}.$$

and therefore

$$p\left(\left(a+q\right)\delta\_1 + \left(aS\_2S\_1^{-1} - q\right)\delta\_2\right)\left(\lambda\_1S\_1 + \lambda\_2S\_2 + q\right) = p\left(\lambda\_2\delta\_1 - \lambda\_1\delta\_2\right).\tag{A1}$$

From the last equation of the system (2) and *R*<sup>∗</sup> = 0 it follows that

$$S\_2 S\_1^{-1} = \frac{1}{\delta\_2} \left( \frac{q}{S\_1} - \delta\_1 \right). \tag{A2}$$

.

Substituting (A2) into (A1) we obtain the equation for *S*<sup>1</sup>

$$\begin{pmatrix} (a+q)\delta\_1 + a\left(\frac{q}{S\_1} - \delta\_1\right) - q\delta\_2\\ q(a - (\delta\_2 - \delta\_1)S\_1)\left(q(\lambda\_2 + \delta\_2) - (\lambda\_2\delta\_1 - \lambda\_1\delta\_2)S\_1\right) = p\delta\_2(\lambda\_2\delta\_1 - \lambda\_1\delta\_2)S\_1\\ q(a - (\delta\_2 - \delta\_1)S\_1)\left(q(\lambda\_2 + \delta\_2) - (\lambda\_2\delta\_1 - \lambda\_1\delta\_2)S\_1\right) = p\delta\_2(\lambda\_2\delta\_1 - \lambda\_1\delta\_2)S\_1 \end{pmatrix}$$
  $\begin{array}{c} \lambda\_1 = q\delta\_2 + q\delta\_1 - q\delta\_2 \end{array}$   $\begin{array}{c} \lambda\_1 = q\delta\_2 - q\delta\_1 + q\delta\_2 - q\delta\_1 \end{array}$ 

If (*δ*<sup>2</sup> <sup>−</sup> *<sup>δ</sup>*1)(*λ*2*δ*<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*1*δ*2) <sup>=</sup> 0 then *<sup>S</sup>*<sup>1</sup> is a positive root of the quadratic equation *<sup>S</sup>*<sup>2</sup> <sup>1</sup> − *ν*1*S*<sup>1</sup> + *ν*<sup>2</sup> = 0, where

$$\nu\_1 = \frac{q\alpha + p\delta\_2}{q(\delta\_2 - \delta\_1)} + \frac{q(\lambda\_2 + \delta\_2)}{\lambda\_2\delta\_1 - \lambda\_1\delta\_2}, \quad \nu\_2 = \frac{\alpha q(\lambda\_2 + \delta\_2)}{(\delta\_2 - \delta\_1)(\lambda\_2\delta\_1 - \lambda\_1\delta\_2)}.$$

**Remark A1.** *Please note that a positive root of the quadratic equation S*<sup>2</sup> <sup>1</sup> − *ν*1*S*<sup>1</sup> + *ν*<sup>2</sup> = 0 *may not exist, for instance, if ν*<sup>1</sup> < 0 *and ν*<sup>2</sup> > 0*. In this case a positive equilibrium E*<sup>4</sup> *does not exist too. On the other hand for some values of the parameters the quadratic equation S*<sup>2</sup> <sup>1</sup> − *ν*1*S*<sup>1</sup> + *ν*<sup>2</sup> = 0 *may have two positive roots, for instance, if ν*<sup>1</sup> > 0 *and* 0 < 4*ν*<sup>2</sup> < *ν*<sup>2</sup> <sup>1</sup> *: S*<sup>∗</sup> <sup>1</sup> = <sup>1</sup> <sup>2</sup> (*ν*<sup>1</sup> ± 3 *ν*2 <sup>1</sup> − 4*ν*2)*. In this case there are two equilibria of the type of E*4*.*

$$\text{If } \delta\_1 = \delta\_2 = \delta\_\prime \lambda\_2 > \lambda\_1 \text{ then } S\_1^\* = \frac{aq^2(\delta + \lambda\_2)}{\delta(\lambda\_2 - \lambda\_1)(qa + p\delta)}. \text{ If } \lambda\_2 \delta\_1 = \lambda\_1 \delta\_2, \delta\_2 > \delta\_1 \text{ then } S\_1^\* = \frac{a}{\delta\_2 - \delta\_1}.$$
 $\text{If } S\_1^\* \text{ is defined then via (2)}$ 

$$S\_2^\* = \frac{1}{\delta\_2} (q - \delta\_1 S\_1^\*), \qquad I^\* = \frac{p}{\lambda\_1 S\_1^\* + \lambda\_2 S\_2^\* + q}, \qquad R^\* = \frac{\lambda\_2 I^\* - a - q}{\delta\_2}.$$

For positivity of the equilibrium *E*<sup>4</sup> must be *S*∗ <sup>1</sup> <sup>&</sup>lt; *<sup>q</sup> δ*1 and *<sup>I</sup>*<sup>∗</sup> <sup>&</sup>gt; *<sup>α</sup>* <sup>+</sup> *<sup>q</sup> λ*2 . *Appendix A.2. Stability of the System of Two Stochastic Differential Equations*

Consider the system of two stochastic differential equations

$$\begin{aligned} \dot{x}\_1(t) &= a\_{11}x\_1(t) + a\_{12}x\_2(t) + \sigma\_1 x\_1(t)\dot{w}\_1(t),\\ \dot{x}\_2(t) &= a\_{21}x\_1(t) + a\_{22}x\_2(t) + \sigma\_1 x\_2(t)\dot{w}\_2(t), \end{aligned} \tag{A3}$$

where *aij*, *σi*, *i*, *j* = 1, 2, are constants, *w*1(*t*) and *w*2(*t*) are mutually independent standard Wiener processes [30,33].

**Lemma A1.** *[30] Put <sup>A</sup>* <sup>=</sup> *aij , <sup>i</sup>*, *<sup>j</sup>* <sup>=</sup> 1, 2*, Ai* <sup>=</sup> det(*A*) + *<sup>a</sup>*<sup>2</sup> *ii, <sup>μ</sup><sup>i</sup>* <sup>=</sup> <sup>1</sup> 2*σ*<sup>2</sup> *<sup>i</sup> , i* = 1, 2*, and suppose that the following conditions hold*

$$\begin{array}{ll} \operatorname{Tr}(A) = a\_{11} + a\_{22} < 0, & \operatorname{det}(A) = a\_{11}a\_{22} - a\_{12}a\_{21} > 0, \\ \mu\_{1} < \frac{|\operatorname{Tr}(A)| \operatorname{det}(A)}{A\_{2}}, & \mu\_{2} < \frac{|\operatorname{Tr}(A)| \operatorname{det}(A) - A\_{2}\mu\_{1}}{A\_{1} - |\operatorname{Tr}(A)|\mu\_{1}}. \end{array} \tag{A4}$$

*Then the zero solution of the system* (A3) *is asymptotically mean square stable.*

**Remark A2.** *Please note that if a*12*a*<sup>21</sup> = 0 *then the last two conditions in* (A4) *take the form μ*<sup>1</sup> < −*a*11*, μ*<sup>2</sup> < −*a*22*.*

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Dynamics of HIV-TB Co-Infection Model**

#### **Nita H Shah 1,\*, Nisha Sheoran 1,\* and Yash Shah <sup>2</sup>**


Received: 3 February 2020; Accepted: 5 March 2020; Published: 11 March 2020

**Abstract:** According to World Health Organization (WHO), the population suffering from human immunodeficiency virus (HIV) infection over a period of time may suffer from TB infection which increases the death rate. There is no cure for acquired immunodeficiency syndrome (AIDS) to date but antiretrovirals (ARVs) can slow down the progression of disease as well as prevent secondary infections or complications. This is considered as a medication in this paper. This scenario of HIV-TB co-infection is modeled using a system of non-linear differential equations. This model considers HIV-infected individual as the initial stage. Four equilibrium points are found. Reproduction number *R*<sup>0</sup> is calculated. If *R*<sup>0</sup> >1 disease persists uniformly, with reference to the reproduction number, backward bifurcation is computed for pre-AIDS (latent) stage. Global stability is established for the equilibrium points where there is no Pre-AIDS TB class, point without co-infection and for the endemic point. Numerical simulation is carried out to validate the data. Sensitivity analysis is carried out to determine the importance of model parameters in the disease dynamics.

**Keywords:** Co-infection of HIV-TB; equilibrium point; reproduction number; stability analysis; backward bifurcation

#### **1. Introduction**

In the public health sector, human immunodeficiency virus (HIV) continues to be the major health threat globally, having claimed more than 32 million lives to date [1]. There were approximately 37.9 million people living with HIV at the end of 2018 [1]. The human immunodeficiency virus (HIV) is a virus that spread through certain body fluids, attacking the body- s immune system, specifically the CD4 cells. The immune function is typically measured by CD4 cell count. Over time, HIV can destroy so many of these cells that the body can't fight against infections and diseases, which paves the way for many opportunistic diseases. One such disease is tuberculosis (TB). It is a contagious disease caused by bacteria called Mycobacterium tuberculosis. The bacteria mostly attack the lungs, but can also damage other parts of the body. The population living with HIV are 15–22 times more likely to develop TB [2]. It is the most commonly occurring illness among HIV-infected individuals, including among those taking antiretroviral treatment (ART). This interaction explains the fact that HIV and TB co-infection is a deadly human syndemic, where syndemic refers to the convergence of two or more diseases that exacerbate the burden of the disease [3]. For the treatment of HIV, HIV drugs called antiretrovirals (ARVs) are advised. ART reduces the risk of TB infection in people living with HIV by 65% [4]. It plays a significant role in preventing TB.

Mathematical modelling has enhanced understanding of disease dynamics. The first compartmental model was given by Kermack and McKendrick [5]. Some basic papers like [6,7] have constructed mathematical models by formulating non-linear differential equation for their respective models and have worked out the critical point/equilibrium points of the respective system and various related properties. In some the related research, many authors have worked out various types of HIV-TB co-infection model. Kirschner et al. [8] developed a model for HIV-1 and TB coinfection inside a host. This was the

first attempt to understand how TB affects the dynamics of HIV-infected individuals. TB is known to be the common serious opportunistic infection occurring in HIV individuals and it occurs in more than 50% of the acquired immunodeficiency syndrome (AIDS) cases in developing countries. Naresh et al. [9] developed a simple nonlinear mathematical model dividing the population into four sub-classes, namely the susceptible, TB-infective, HIV-infective and AIDS patients. The treatment class in the HIV-AIDS co-infection model was first introduced by Huo et al. [10], however, Bhunu et al. [11] in his co-infection model considered all aspects of TB and HIV transmission dynamics with both HIV and TB treatment. This paper incorporated ARTs for AIDS cases and studied its implication on TB. However, the author did not consider the case where individual co-infected with HIV-TB can effectively recover from TB infection. Another HIV-TB co-infection model was formulated by Roeger et al. [12], assuming TB-infected individuals in the active stage of disease to be sexually inactive. Singh et al. [13] studied the transmission dynamics of the HIV/AIDS epidemic model considering three different latent stages based on treatment. Torres et al. [14], in his model, incorporates both TB and AIDS treatment for individuals suffering with either or both disease.

The model formulated in this paper considers the susceptible class to be HIV-infected. The paper is organized as follows. The model is formulated and its description is given in Section 2. Calculation of reproduction number and uniform persistence of the disease is shown in Section 2.3. In Section 2.4, global stability for all the equilibrium points is done. In Section 3, backward bifurcation is established. The sensitivity of reproduction number is done in Section 3.1. Section 3.2 presents a numerical simulation. The paper concludes in Section 4.

#### **2. Mathematical Model**

We begin with seven mutually exclusive compartmental models showing HIV-TB co-infection. In this model, the human population is divided into sub-populations as follows: acute HIV-infected individuals (*H*), co-infected with HIV-TB (*HTB*), Pre-AIDS stage(*PA*), infected individuals undergoing any type of treatment say ARV's and any TB treatment (*M*), Pre-AIDS stage with TB disease (*PATB*), HIV-infected individuals showing clinical AIDS symptoms (*A*), HIV-infected individuals with AIDS symptoms coinfected with TB disease (*ATB*).

The notations and parametric values assumed in the paper for the study of dynamical system of HIV-TB co-infection model is tabulated in Table 1.


**Table 1.** Parametric definitions and its values.

In this paper, the susceptible class is considered to be HIV-infected (acute HIV infection). This class is increased by recruitments of newly HIV-infected individuals at the rate *B*. All the individuals in their respective compartments suffer from natural death at the constant rate μ. Individuals undergoing medication (treatment) through ARTs lower the rate of progression from HIV disease to AIDS, as HIV can never be cured.

Here, the individuals infected with HIV develop a very weak immune system, which means they are likely to get infected by many opportunistic diseases. As TB is considered to be one of the most commonly occurring disease among HIV patients [15], the individual infected with HIV gets TB disease moving towards *HTB* by rate β1. The HIV-infected individuals are also assumed to progress to the asymptotic pre-AIDS class (*PA*) at the rate β2. The HIV-infected and co-infected individuals undergoing treatment move to class *M* with the rates β<sup>3</sup> and β4, respectively. Similarly, individuals with a co-infection of HIV-TB move towards *PATB* with the rate β5. Individuals showing symptoms of AIDS (*PA*) suffer from full-blown AIDS, joining *A*, at the rate β9, and they are more likely to develop TB, progressing to class *PATB* with the rate β8. *PA* class individuals undergoing ARTs treatment (anti-retroviral therapy) join *M* at the rate β6. Individuals in *PATB* are treated for TB at the constant rate β7, joining *M*, and some of them can also develop full-blown AIDS, moving to *ATB* class with the constant rate β10. Treated individuals, recovered from TB but still with HIV infection (as it cannot be cured) move to full-blown AIDS (*A*) with the constant rate β11. Individuals suffering with AIDS have such a badly damaged immune system that they get an increasing number of severe illnesses (here, TB) and hence move towards full-blown AIDS-TB class with the constant rate β12. The death rates μ*<sup>D</sup>* and μ*DTB* are considered as deaths due to individuals infected with AIDS and AIDS-TB, respectively.

#### *2.1. HIV-TB Co-infection Model*

Considering the aforementioned assumptions and Figure 1 gives rise to the following set of non-linear differential equations for the HIV-TB co-infection model:

$$\begin{aligned} \frac{dM}{dt} &= B - \beta\_1 HH\_{TB} - (\beta\_2 + \beta\_3 + \mu)H\\ \frac{dH\_{TB}}{dt} &= \beta\_1 HH\_{TB} - \beta\_5 H\_{TB}P\_{ATB} - (\mu + \beta\_4)H\_{TB} \\ \frac{dP\_A}{dt} &= \beta\_2 H - \beta\_8 P\_A P\_{ATB} - (\mu + \beta\_6 + \beta\_9)P\_A\\ \frac{dH}{dt} &= \beta\_3 H + \beta\_4 H\_{TB} + \beta\_6 P\_A + \beta\_7 P\_{ATB} - (\mu + \beta\_{11})M \\ \frac{dP\_{ATB}}{dt} &= \beta\_5 H\_{TB}P\_{ATB} + \beta\_8 P\_A P\_{ATB} - (\mu + \beta\_7 + \beta\_{10})P\_{ATB} \\ \frac{dA}{dt} &= \beta\_9 P\_A + \beta\_{11}M - (\mu + \mu\_D + \beta\_{12})A \\ \frac{dA\_{TB}}{dt} &= \beta\_{10} P\_{ATB} + \beta\_{12}A - (\mu + \mu\_{DTB})A\_{TB} \end{aligned} \tag{1}$$

where *N*(*t*) = *H*(*t*) + *HTB*(*t*) + *PA*(*t*) + *M*(*t*) + *PATB*(*t*) + *A*(*t*) + *ATB*(*t*). The system satisfies the conditions:

$$H(t) \ge 0,\\ H\_{\rm TB}(t) \ge 0,\\ P\_A(t) \ge 0,\\ M(t) \ge 0,\\ P\_{\rm ATB}(t) \ge 0,\\ A(t) \ge 0,\\ A\_{\rm TB}(t) \ge 0$$

Adding the above set of differential equations, we get,

$$\begin{aligned} \frac{dN(t)}{dt} &= B - \mu(H + H\_{TB} + P\_A + M + P\_{ATB} + A + A\_{TB}) - \mu\_D A - \mu\_{DTB} A\_{TB} \\ &\le B - \mu(H + H\_{TB} + P\_A + M + P\_{ATB} + A + A\_{TB}) \end{aligned}$$

Hence, *dN*(*t*) *dt* <sup>≤</sup> *<sup>B</sup>* <sup>−</sup> <sup>μ</sup>*N*, so that lim*t*→∞sup*<sup>N</sup>* <sup>≤</sup> *<sup>B</sup>* μ

The feasible region for the system is defined as

$$\Lambda = \left\{ (H\_\prime H\_{\rm TB}, P\_{\rm A\prime} M\_\prime P\_{\rm ATB}, A\_\prime A\_{\rm TB}) : 0 \le H + H\_{\rm TB} + P\_A + M + P\_{\rm ATB} + A + A\_{\rm TB} \le \frac{B}{\mu} \right\}$$

We assume *L*<sup>1</sup> = β<sup>2</sup> + β3, *L*<sup>2</sup> = β<sup>6</sup> + β9, *L*<sup>3</sup> = β<sup>7</sup> + β10. The modified system is

$$\begin{aligned} \frac{d\mathcal{H}}{dt} &= B - \beta\_1 H H \tau\_B - (L\_1 + \mu) H \\ \frac{d\mathcal{H}\_{TB}}{dt} &= \beta\_1 H H\_{TB} - \beta\_5 H\_{TB} P\_{ATB} - (\mu + \beta\_4) H\_{TB} \\ \frac{d\mathcal{P}\_A}{dt} &= \beta\_2 H - \beta\_8 P\_A P\_{ATB} - (L\_2 + \mu) P\_A \\ \frac{d\mathcal{M}}{dt} &= \beta\_3 H + \beta\_4 H\_{TB} + \beta\_6 P\_A + \beta\_7 P\_{ATB} - (\mu + \beta\_{11}) M \\ \frac{d\mathcal{P}\_{ATB}}{dt} &= \beta\_5 H\_{TB} P\_{ATB} + \beta\_8 P\_A P\_{ATB} - (L\_3 + \mu) P\_{ATB} \\ \frac{d\mathcal{A}}{dt} &= \beta\_9 P\_A + \beta\_{11} M - (\mu + \mu\_D + \beta\_{12}) A \\ \frac{d\mathcal{A}\_{TB}}{dt} &= \beta\_{10} P\_{ATB} + \beta\_{12} A - (\mu + \mu\_{DTB}) A\_{TB} \end{aligned} \tag{2}$$

System (1) and (2) are equivalent, hence Λ is also the feasible region for system (2).

**Figure 1.** Transmission of individuals in different compartments.

#### *2.2. Equilibrium Solutions*

Equating *dH dt* <sup>=</sup> *dHTB dt* <sup>=</sup> *dPA dt* <sup>=</sup> *dM dt* <sup>=</sup> *dPATB dt* <sup>=</sup> *dA dt* <sup>=</sup> *dATB dt* = 0 and solving for the compartments following are the equilibria:

1. *E*1(*H*1, 0, *PA*<sup>1</sup> , *M*1, 0, *A*1, *ATB*<sup>1</sup> )

$$\begin{array}{c} H\_{1} = \frac{\mathcal{B}}{L\_{1} + \mu}, H\_{TB\_{1}} = 0, P\_{A1} = \frac{\mathcal{B}\beta\_{2}}{(L\_{1} + \mu)(L\_{2} + \mu)}, M\_{1} = \frac{\mathcal{B}(L\_{2}\beta\_{3} + \beta\_{2}\beta\_{6} + \beta\_{3}\mu)}{(L\_{1} + \mu)(L\_{2} + \mu)(\beta\_{11} + \mu)}, \\ P\_{ATB\_{1}} = 0, A\_{1} = \frac{\mathcal{B}(\beta\_{11}(L\_{2}\beta\_{3} + \beta\_{2}\beta\_{6} + \beta\_{3}\mu) + \beta\_{2}\beta\_{7}(\beta\_{11} + \mu))}{(L\_{1} + \mu)(L\_{2} + \mu)(\beta\_{11} + \mu)(\mu + \mu\_{D} + \beta\_{12})}, \\ A\_{TB\_{1}} = \frac{\mathcal{B}\beta\_{12}(\beta\_{11}(L\_{2}\beta\_{3} + \beta\_{2}\beta\_{6} + \beta\_{3}\mu) + \beta\_{2}\beta\_{9}(\beta\_{11} + \mu))}{(L\_{1} + \mu)(L\_{2} + \mu)(\beta\_{11} + \mu)(\mu + \mu\_{D} + \beta\_{12})(\mu + \mu\_{D})}. \end{array}$$

2. *E*2(*H*2, 0, *PA*<sup>2</sup> , *M*2, *PATB*<sup>2</sup> , *A*2, *ATB*<sup>2</sup> )

*H*<sup>2</sup> = *<sup>B</sup> <sup>L</sup>*1+<sup>μ</sup> , *HTB*<sup>2</sup> <sup>=</sup> 0, *PA*<sup>2</sup> <sup>=</sup> *<sup>L</sup>*3+<sup>μ</sup> <sup>β</sup><sup>8</sup> , *PATB*<sup>2</sup> <sup>=</sup> *<sup>B</sup>*β2β8−(*L*1+μ)(*L*2+μ)(*L*3+μ) <sup>β</sup>8(*L*1+μ)(*L*3+μ) , *<sup>M</sup>*<sup>2</sup> <sup>=</sup> *<sup>B</sup>*β8(*L*3β3+β2β7+β3μ)+(*L*1+μ)(*L*3+μ)(β6(*L*3+μ)−β7(*L*2+μ)) <sup>β</sup>8(*L*1+μ)(*L*3+μ)(β11+μ) , *A*<sup>2</sup> = *B*β8β11(*L*3β<sup>3</sup> + β2β<sup>7</sup> + β3μ)+(*L*<sup>3</sup> + μ) 2 (*L*<sup>1</sup> + μ) (β9(β<sup>11</sup> + μ) + β6β11) − (*L*<sup>1</sup> + μ)(*L*<sup>2</sup> + μ)(*L*<sup>3</sup> + μ)β7β11 <sup>β</sup>8(*L*1+μ)(*L*3+μ)(μ+μ*D*+β12) , *ATB*<sup>2</sup> = *B*β8β11β12(β2β<sup>7</sup> + *L*3β<sup>3</sup> + μβ3)+(β<sup>11</sup> + μ)μ*DB*β2β8β<sup>10</sup> − (*L*<sup>1</sup> + μ)(*L*<sup>2</sup> + μ) (*L*<sup>3</sup> + μ)(β7β11β<sup>12</sup> + (β<sup>11</sup> + μ)μ*D*β10)+(β<sup>12</sup> + μ)(β<sup>11</sup> + μ)*B*β2β8β<sup>10</sup> + (*L*<sup>1</sup> + μ) (*L*<sup>3</sup> + μ)(β6β<sup>11</sup> + (β<sup>11</sup> + μ)β9β12) − (*L*<sup>1</sup> + μ)(*L*<sup>2</sup> + μ)(*L*<sup>3</sup> + μ)(β<sup>11</sup> + μ)(β<sup>12</sup> + μ) β8(*L*1+μ)(*L*3+μ)(β11+μ)(μ+μ*DTB*)(μ+μ*D*+β12)

#### 3. *E*3(*H*3, *HTB*<sup>3</sup> , *PA*<sup>3</sup> , *M*3, 0, *A*3, *ATB*<sup>3</sup> )

*<sup>H</sup>*<sup>3</sup> <sup>=</sup> <sup>β</sup>4+<sup>μ</sup> <sup>β</sup><sup>1</sup> , *HTB*<sup>3</sup> <sup>=</sup> *<sup>B</sup>*β1−(*L*1+μ)(β4+μ) <sup>β</sup>1(β4+μ) , *PA*<sup>3</sup> <sup>=</sup> <sup>β</sup>2(β4+μ) <sup>β</sup>1(*L*2+μ), *<sup>M</sup>*<sup>3</sup> <sup>=</sup> (*L*2+μ)(*B*β1β4−(β4+μ)(*L*1β4−β3β4−β3μ−β4μ))+β2β6(β4+μ) 2 <sup>β</sup>1(*L*2+μ)(β4+μ)(β11+μ) , *PATB*<sup>3</sup> = 0, *A*<sup>3</sup> = β4β11(*B*β<sup>1</sup> − *L*1(β<sup>4</sup> + μ))(*L*<sup>2</sup> + μ)+(β<sup>4</sup> + μ)(*L*2β11(β3β<sup>4</sup> + β3μ − β4μ) <sup>−</sup>β4β11μ2)+(β<sup>4</sup> + <sup>μ</sup>) 2 (β2β9(β<sup>11</sup> + μ) + β2β6β<sup>11</sup> + β3β11μ) β1(*L*2+μ)(β4+μ)(β11+μ)(μ+μ*D*+β12) *ATB*<sup>3</sup> = β12(β4β11(*B*β<sup>1</sup> − *L*1(β<sup>4</sup> + μ))(*L*<sup>2</sup> + μ)+(β<sup>4</sup> + μ)(*L*2β11(β3β<sup>4</sup> + β3μ − β4μ) <sup>−</sup>β4β11μ2)) + (β<sup>4</sup> + <sup>μ</sup>) 2 (β2β9(β<sup>11</sup> + μ) + β2β6β<sup>11</sup> + β3β11μ) β1(*L*2+μ)(β4+μ)(β11+μ)(μ+μ*D*+β12)(μ+μ*DTB*)

4. Endemic Equilibrium point *E*∗ (*H*∗ , *H*∗ *TB*, *P*<sup>∗</sup> *<sup>A</sup>*, *M*<sup>∗</sup> , *P*∗ *ATB*, *A*<sup>∗</sup> , *A*∗ *TB*)

$$\begin{split} H^\* &= \frac{r(\beta\_8(\beta\_4+\mu) - \beta\_5(L\_2+\mu)) + \beta\beta\_5}{\beta\varepsilon(L\_1+\mu) + \beta\_1(L\_3+\mu) - \beta\_2\beta\_5}, H^\*\_{TB} = \frac{-\beta\_8r + L\_3 + \mu}{\beta\_5}, P^\*\_A = r, \\ &\qquad r\beta\_1(\beta\_8(\beta\_4+\mu) - \beta\_5(L\_2+\mu)) + \beta\_5(B\beta\_1} \\ P^\*\_{ATB} &= \frac{-(L\_1+\mu)(\beta\_4+\mu) + \beta\_2(\beta\_4+\mu)) - \beta\_1(L\_3+\mu)(\beta\_4+\mu)}{\beta\varepsilon(L\_1\beta\mathfrak{s}+L\_3\mathfrak{s}\_1 + \beta\_1\mathfrak{s}\_1 - \beta\_2\beta\_5 + \beta\_5\mu)}. \end{split}$$

*M*∗ = *r*[(β5β<sup>6</sup> − β4β8)[β5(*L*<sup>1</sup> + μ) + β1(*L*<sup>3</sup> + μ) − β2β5]+(β3β<sup>5</sup> + β7β11)(β8(β<sup>4</sup> + μ) −β5(*L*<sup>2</sup> + μ))] + β5β7[(*B*β<sup>1</sup> − (*L*<sup>1</sup> + μ)(β<sup>4</sup> + μ)) + β2(β<sup>4</sup> + μ)] + β4β5(*L*<sup>1</sup> + μ) (*L*<sup>3</sup> + <sup>μ</sup>)+(*L*<sup>3</sup> + <sup>μ</sup>)β1(β4(*L*<sup>3</sup> + <sup>μ</sup>) <sup>−</sup> <sup>β</sup>7(β<sup>4</sup> + <sup>μ</sup>)) <sup>−</sup> <sup>β</sup>2β4β5(*L*<sup>3</sup> + <sup>μ</sup>) + *<sup>B</sup>*β3β<sup>2</sup> 5 <sup>β</sup>5(β11+μ)(β5(*L*1+μ)+β1(*L*3+μ)−β2β5) , *A*∗ = *r*(β11(β1β<sup>7</sup> + β3β5)(β8(β<sup>4</sup> + μ) − β5(*L*<sup>2</sup> + μ)) + (β1(*L*<sup>3</sup> + μ) + β5(*L*<sup>1</sup> + μ) − β2β5) ((β<sup>11</sup> + μ)β5β<sup>9</sup> + (β5β<sup>6</sup> − β4β8)β11)) + β11(β5β7(*B*β<sup>1</sup> − (*L*<sup>1</sup> + μ)(β<sup>4</sup> + μ)) + β2(β<sup>4</sup> + μ)) +β4β5(*L*<sup>1</sup> + <sup>μ</sup>)(*L*<sup>3</sup> + <sup>μ</sup>) + <sup>β</sup>1(*L*<sup>3</sup> + <sup>μ</sup>)(β4(*L*<sup>3</sup> + <sup>μ</sup>) <sup>−</sup> <sup>β</sup>7(β<sup>4</sup> + <sup>μ</sup>)) <sup>−</sup> <sup>β</sup>2β4β5(*L*<sup>3</sup> + <sup>μ</sup>) + *<sup>B</sup>*β3β<sup>2</sup> 5 <sup>β</sup>5(μ+β11)(μ+μ*D*+β12)(β5(*L*1+μ)+β1(*L*3+μ)−β2β5) , *A*∗ *TB* = *r*[(β8(β<sup>4</sup> + μ) − β5(*L*<sup>2</sup> + μ))((β<sup>11</sup> + μ)β1β10(μ + μ*<sup>D</sup>* + β12) + β11β12(β3β<sup>5</sup> + β1β7)) +(β5(*L*<sup>1</sup> + μ) + β1(*L*<sup>3</sup> + μ) − β2β5)β12((β<sup>4</sup> + μ)β5β<sup>9</sup> + (β5β<sup>6</sup> − β4β8)β11)] + μ*D*β<sup>10</sup> (β<sup>11</sup> + μ)(β<sup>4</sup> + μ)(β2β<sup>5</sup> − β1(*L*<sup>3</sup> + μ) − β5(*L*<sup>1</sup> + μ)) + μ*DB*β1β5β10(β<sup>11</sup> + μ) + β7β11β<sup>12</sup> (β<sup>4</sup> + μ)(β2β<sup>5</sup> − β5μ − *L*3β1) − β10(β<sup>12</sup> + μ)(β<sup>11</sup> + μ) 2 (β5(*L*<sup>1</sup> + μ) + β1(*L*<sup>3</sup> + μ)) + β<sup>4</sup> β11β12(*L*<sup>3</sup> + μ)(β5(*L*<sup>1</sup> + μ) + β1(*L*<sup>3</sup> + μ) − β2β5) + β5β10(β<sup>11</sup> + μ)(β<sup>12</sup> + μ)(β2(β<sup>4</sup> + μ) +*B*β1) − β1β7β11β12(β<sup>5</sup> + μ)(β<sup>4</sup> + μ) + *B*β5β11β12(β1β<sup>7</sup> + β3β5) β5(μ+β11)(μ+μ*DTB*)(μ+μ*D*+β12)(β5(*L*1+μ)+β1(*L*3+μ)−β2β5)

$$\text{where,}\begin{aligned} r &= \text{root}(( (\beta\_1 \beta\_8 (\beta\_5 (L\_2 + \mu) - \beta\_8 (\beta\_4 + \mu)) ) ) \mathbf{Z}^2 + (-\beta\_5 \beta\_8 (\mathbf{B} \beta\_1 - (L\_1 + \mu)(\beta\_4 + \mu)) \\ -\beta\_5^2 (L\_1 + \mu)(L\_2 + \mu) + \beta\_1 (L\_3 + \mu)(\beta\_8 (\beta\_4 + \mu) - \beta\_5 (L\_2 + \mu)) ) \mathbf{Z} + \mathbf{B} \beta\_2 \beta\_5^2 \end{aligned}$$

#### *2.3. Reproduction Number*

The reproduction number measures the expected number of secondary infected individuals produced due to an infected individual during the entire death period in an uninfected population.

In this paper, reproduction number *R*<sup>0</sup> is defined as the number of infected individuals due to an AIDS- or TB-infected individual in the HIV infected-population. It is calculated using next-generation matrix method [16] and is defined as the spectral radius of *FV*−<sup>1</sup> at *E*1.


$$V = \begin{bmatrix} \beta\_4 + \beta\_5 P\_{ATB} + \mu & 0 & 0 & \beta\_5 H\_{TB} & 0 & 0 & 0 \\ 0 & \beta \mathfrak{s} P\_{ATB} + L\_2 + \mu & 0 & \beta \mathfrak{s} P\_A & 0 & 0 & -\beta\_2 \\ -\beta & -\beta\_6 & \beta\_{11} + \mu & -\beta\tau & 0 & 0 & -\beta\_3 \\ 0 & 0 & 0 & L\_3 + \mu & 0 & 0 & 0 \\ 0 & -\beta\_9 & -\beta\_{11} & 0 & \beta\_{12} + \mu\_D + \mu & 0 & 0 \\ 0 & 0 & 0 & -\beta\_{10} & -\beta\_{12} & \mu + \mu\_{DTB} & 0 \\ \beta\_1 H & 0 & 0 & 0 & 0 & 0 & \beta\_1 H\_{TB} + L\_1 + \mu \end{bmatrix}$$

The dominant eigenvalue of *FV*−<sup>1</sup> at *<sup>E</sup>*<sup>1</sup> is *<sup>R</sup>*<sup>0</sup> <sup>=</sup> *<sup>B</sup>*β2β<sup>8</sup> (*L*1+μ)(*L*2+μ)(*L*3+μ) <sup>+</sup> <sup>β</sup>1*<sup>B</sup>* (*L*1+μ)(β4+μ).

#### *2.4. Persistence of Disease*

Now, uniform persistence for the system (1) is constructed. The model system (1) is said to be uniformly persistent if there is a constant *f*, such that any solution (*H*(*t*), *HTB*(*t*), *PA*(*t*), *M*(*t*), *PATB*(*t*), *A*(*t*), *ATB*(*t*)) satisfies [17,18].

$$\lim\_{t \to \infty} \inf H(t) > f\_{\prime} \lim\_{t \to \infty} \inf H\_{TB}(t) > f\_{\prime} \lim\_{t \to \infty} \inf P\_{A}(t) > f\_{\prime} \lim\_{t \to \infty} \inf M(t) > f\_{\prime}.$$
 
$$\lim\_{t \to \infty} \inf P\_{ATB}(t) > f\_{\prime} \lim\_{t \to \infty} \inf A(t) > f\_{\prime} \lim\_{t \to \infty} \inf A\_{TB}(t) > f.$$

Provided that (*H*(0), *HTB*(0), *PA*(0), *M*(0), *PATB*(0), *A*(0), *ATB*(0)) ∈ Λ

**Theorem 1.** *The model (1) is uniformly persistent in* Λ *only if R*<sup>0</sup> > 1*.*

#### *2.5. Stability Analysis*

In this section, global stability is studied for all the equilibrium points obtained.

**Theorem 2.** *Global Stability of E*1(*H*1, 0, *PA*<sup>1</sup> , *M*1, 0, *A*1, *ATB*<sup>1</sup> )

The system (2) of the model can be written as

$$\frac{dX\_1}{dt} = F\_1(X\_1, Z\_1) \tag{3}$$

$$\frac{dZ\_1}{dt} = G\_1(X\_1, Z\_1), G\_1(X\_1, 0) = 0\tag{4}$$

where *X*<sup>1</sup> = (*H*, *PA*, *M*, *A*, *ATB*) and *Z*<sup>1</sup> = (*HTB*, *PATB*). According to this notation, equilibrium point is denoted by *E*<sup>1</sup> = (*X*- <sup>1</sup>, 0), *where X*- <sup>1</sup> = (*H*1, 0, *PA*<sup>1</sup> , *M*1, 0, *A*1, *ATB*<sup>1</sup> ).

By the Castillo Chavez method, the following two condition ensure the global stability of the given equilibrium point:

P.1  $\text{For } \frac{dX\_1}{dt} = F\_1(X\_1, 0)$ ,  $E\_1$  is globally asymptotically stable.

P.2  $G\_1(X\_1, Z\_1) = AZ\_1 - \hat{G}\_1(X\_1, Z\_1)$ , where  $\hat{G}\_1(X\_1, Z\_1) \ge 0$  for  $(X\_1, Z\_1) \in \Lambda$ .

where *A* = *DZ*1*G*1(*X*1, 0) is a M-matrix (matrix with non-negative off diagonal elements) and Λ is the region defined above. We have,

$$F\_1(X\_1, 0) = \begin{bmatrix} B - (\beta\_2 + \beta\_3 + \mu)H \\ \beta\_2 H - (\mu + \beta\_6 + \beta\_7)P\_A \\ \beta\_3 H + \beta\_6 P\_A - (\mu + \beta\_{11})M \\ \beta\_9 P\_A + \beta\_{11}M - (\mu + \mu\_D + \beta\_{12})A \\ \beta\_{12}A - (\mu + \mu\_{DTB})A\_{TB} \end{bmatrix}$$

$$\mathbf{G}\_{1}(\mathbf{X}\_{1},\mathbf{Z}\_{1}) = \begin{bmatrix} \beta\_{1}H\mathbf{H}\_{\mathrm{TB}} - \beta\_{5}\mathbf{H}\_{\mathrm{TB}}\mathbf{P}\_{\mathrm{AT}} - (\mu + \beta\_{4})H\_{\mathrm{TB}} \\ \beta\_{5}H\_{\mathrm{TB}}\mathbf{P}\_{\mathrm{AT}} + \beta\_{8}P\_{\mathrm{A}}\mathbf{P}\_{\mathrm{AT}} - (\mu + \beta\_{7} + \beta\_{10})P\_{\mathrm{AT}} \\ \end{bmatrix} \text{and } \mathbf{G}\_{1}(\mathbf{X}\_{1},0) = \text{0, thus}$$

$$A = D\_{\mathrm{Z}\_{1}}\mathbf{G}\_{1}(\mathbf{X}\_{1}^{\prime},0) = \begin{bmatrix} \beta\_{1}H - (\mu + \beta\_{4}) & 0 \\ 0 & \beta\_{8}P\_{\mathrm{A}} - (\mu + \beta\_{7} + \beta\_{10}) \end{bmatrix}$$

$$\mathbf{G}(\mathbf{X}\_{1},\mathbf{Z}\_{1}) = \begin{bmatrix} \beta\_{5}H\_{\mathrm{TB}}P\_{\mathrm{AT}} \\ -\beta\_{5}H\_{\mathrm{TB}}P\_{\mathrm{AT}} \end{bmatrix} \tag{5}$$

From Equation (5), the condition P.2 is not satisfied, since *<sup>G</sup>*<sup>ˆ</sup> <sup>1</sup>(*X*1,*Z*1) <sup>≥</sup> 0 is not true. Therefore, the equilibrium point *E*<sup>1</sup> may not be globally stable. Here, since disease (HIV-AIDS) persists at this point, it will not be globally stable. Following [19], the backward bifurcation occurs at *R*<sup>0</sup> = 1.

**Theorem 3.** *Global Stability of E*2(*H*2, 0, *PA*<sup>2</sup> , *M*2, *PATB*<sup>2</sup> , *A*2, *ATB*<sup>2</sup> )

The system (2) of the model can be written as

$$\frac{dX\_2}{dt} = F\_2(X\_2, Z\_2) \tag{6}$$

$$\frac{dZ\_2}{dt} = G\_2(X\_2, Z\_2), \ G\_2(X\_2, 0) = 0 \tag{7}$$

where *X*<sup>2</sup> = (*H*, *PA*, *M*, *PATB*, *A*, *ATB*) and *Z*<sup>2</sup> = (*HTB*). According to this notation, the equilibrium point is denoted by *E*<sup>2</sup> = (*X*- <sup>2</sup>, 0), *where X*- <sup>2</sup> = (*H*2, 0, *PA*<sup>2</sup> , *M*2, *PATB*<sup>2</sup> , *A*2, *ATB*<sup>2</sup> ).

Using the Castillo Chavez method [20], the following two condition ensure the global stability of the given equilibrium point:

P.3 For *dX*<sup>2</sup> *dt* = *F*2(*X*2, 0), *E*<sup>2</sup> is globally asymptotically stable. P.4 *<sup>G</sup>*2(*X*2,*Z*2) = *BZ*<sup>2</sup> <sup>−</sup> *<sup>G</sup>*<sup>ˆ</sup> <sup>2</sup>(*X*2,*Z*2), *where <sup>G</sup>*ˆ(*X*2,*Z*2) <sup>≥</sup> 0 for (*X*2,*Z*2) <sup>∈</sup> <sup>Λ</sup>.

where *B* = *DZ*2*G*2(*X*2, 0) is an M-matrix (matrix with non-negative off diagonal elements) and Λ is the region defined above.

The equilibrium point *E*2(*H*2, 0, *PA*<sup>2</sup> , *M*2, *PATB*<sup>2</sup> , *A*2, *ATB*<sup>2</sup> ) is the globally asymptotically stable equilibrium of the system (P.3)–(P.4)

$$\text{we have } F\_2(X\_2, 0) = \begin{bmatrix} B - (\beta\_2 + \beta\_3 + \mu)H \\ \beta\_2 H - \beta\_8 P\_A P\_{ATB} - (\mu + \beta\_6 + \beta\_9)P\_A \\ \beta\_3 H + \beta\_6 P\_A + \beta\_7 P\_{ATB} - (\mu + \beta\_{11})M \\ \beta\_8 P\_A P\_{ATB} - (\mu + \beta\_7 + \beta\_{10})P\_{ATB} \\ \beta\_9 P\_A + \beta\_{11}M - (\mu + \mu\_D + \beta\_{12})A \\ \beta\_{10}P\_{ATB} + \beta\_{12}A - (\mu + \mu\_{DTB})A\_{TB} \end{bmatrix}$$

The eigenvalues of the characteristic polynomial of its Jacobian matrix are given as

$$\begin{aligned} \lambda\_1 &= - (\mu + \beta\_2 + \beta\_3), \lambda\_2 = - (\mu + \beta\_{11}), \lambda\_3 = - (\mu + \mu\_{\rm DTB}), \lambda\_4 = - (\mu + \mu\_{\rm D} + \beta\_{12}), \\ \lambda\_5 &= -\frac{1}{2} ( (\beta\_8 P\_{\rm ATB} + \beta\_6 + \beta\_7 + \beta\_9 + \beta\_{10} + 2\mu - \beta\_8 P\_A) \\ - \sqrt{\beta\_8^2 (P\_A - P\_{\rm ATB})^2 + \left(\beta\_6 - \beta\_7 + \beta\_9 - \beta\_{10}\right)^2 + 2\beta\_8 (P\_A + P\_{\rm ATB}) \left(\beta\_6 - \beta\_7 + \beta\_9 - \beta\_{10}\right)}, \\ \lambda\_6 &= -\frac{1}{2} ( (\beta\_8 P\_{\rm ATB} + \beta\_6 + \beta\_7 + \beta\_9 + \beta\_{10} + 2\mu - \beta\_8 P\_A) \\ + \sqrt{\beta\_8^2 (P\_A - P\_{\rm ATB})^2 + \left(\beta\_6 - \beta\_7 + \beta\_9 - \beta\_{10}\right)^2 + 2\beta\_8 (P\_A + P\_{\rm ATB}) \left(\beta\_6 - \beta\_7 + \beta\_9 - \beta\_{10}\right)} \end{aligned}$$

Here, λ5, λ<sup>6</sup> have a negative real part if (β<sup>6</sup> + β<sup>9</sup> + μ)β8*PA* < (β<sup>7</sup> + β<sup>10</sup> + μ)(β<sup>6</sup> + β<sup>9</sup> + μ + β8*PATB*). Hence, by Routh–Hurwitz criterion, the system is globally asymptotically stable. Next,

$$\begin{aligned} \mathsf{G}\_2(\mathsf{X}\_2, \mathsf{Z}\_2) &= (\beta\_1 H\_2 + \beta\_5 P\_{ATB\_2} - (\mu + \beta\_4)) H\_{TB} - (\beta\_1 (H\_2 - H) H\_{TB} + \beta\_5 (P\_{ATB\_2} - P\_{ATB}) H\_{TB}) \\ &= B H\_{TB} - \mathsf{G}\_2(\mathsf{X}\_2, \mathsf{Z}\_2) \end{aligned}$$

Here, *<sup>G</sup>*<sup>ˆ</sup> <sup>2</sup>(*X*2,*Z*2) <sup>≥</sup> 0, hence the conditions of P.3 and P.4 are satisfied. Hence, by Castillo Chavez the system is globally stable.

**Theorem 4.** *Global stability of E*3(*H*3, *HTB*<sup>3</sup> , *PA*<sup>3</sup> , *M*3, 0, *A*3, *ATB*<sup>3</sup> )

The system (1) of the model can be written as

$$\frac{dX\_3}{dt} = F\_3(X\_3, Z\_3) \tag{8}$$

$$\frac{dZ\_3}{dt} = G\_3(X\_3, Z\_3), \ G\_3(X\_3, 0) = 0 \tag{9}$$

where *X*<sup>3</sup> = (*H*, *HTB*, *PA*, *M*, *A*, *ATB*) and *Z*<sup>3</sup> = (*PATB*). According to this notation, the equilibrium point is denoted by *E*<sup>3</sup> = (*X*- <sup>3</sup>, 0), *where X*- <sup>3</sup> = (*H*3, *HTB*<sup>3</sup> , *PA*<sup>3</sup> , *M*3, 0, *A*3, *ATB*<sup>3</sup> ).

The following two conditions ensure the global stability of this equilibrium point

P.5 For  $\frac{dX\_3}{dt} = F\_3(X\_3, 0)$ ,  $E\_3$  is globally asymptotically stable.

P.6 G $\_3(X\_3, Z\_3) = CZ\_3 - \hat{G}\_3(X\_3, Z\_3)$ , where  $\hat{G}\_3(X\_3, Z\_3) \ge 0$  for  $(X\_3, Z\_3) \in \Lambda$ 

where *C* = *DZ*3*G*3(*X*3, 0) is an M-matrix (matrix with non-negative off diagonal elements) and Λ is the region defined above.

The equilibrium point *E*3(*H*3, *HTB*<sup>3</sup> , *PA*<sup>3</sup> , *M*3, 0, *A*3, *ATB*<sup>3</sup> ) is the globally asymptotically stable equilibrium of the system (P.5)–(P.6)

$$\text{we have } F\_3(X\_3, 0) = \begin{bmatrix} B - \beta\_1 HH\_{TB} - (\beta\_2 + \beta\_3 + \mu)H \\ \beta\_1 HH\_{TB} - (\mu + \beta\_4)H\_{TB} \\ \beta\_2 H - (\mu + \beta\_6 + \beta\_9)P\_A \\ \beta\_3 H + \beta\_4 H\_{TB} + \beta\_6 P\_A - (\mu + \beta\_{11})M \\ \beta\_3 P\_A + \beta\_{11}M - (\mu + \mu\_D + \beta\_{12})A \\ \beta\_{12}A - (\mu + \mu\_{DTB})A\_{TB} \end{bmatrix}$$

The eigenvalues of the characteristic polynomial of its Jacobian matrix are given as

$$\begin{aligned} \lambda\_1 &= - (\mu + \beta\_6 + \beta\_9), \lambda\_2 = - (\mu + \beta\_{11}), \lambda\_3 = - (\mu + \mu\_{\rm{DTB}}), \\ \lambda\_4 &= - (\mu + \mu\_D + \beta\_{12}), \\ \lambda\_5 &= -\frac{1}{2} ( (\beta\_1 H \tau \mu + \beta\_2 + \beta\_3 + \beta\_4 + 2\mu - \beta\_1 H) \\ &- \sqrt{\beta\_1^2 (H - H\_{TB})^2 + (\beta\_2 + \beta\_3 - \beta\_4)^2 + 2\beta\_1 (H + H\_{TB}) (\beta\_2 + \beta\_3 - \beta\_4)} \\ \lambda\_6 &= -\frac{1}{2} ( (\beta\_1 H\_{TB} + \beta\_2 + \beta\_3 + \beta\_4 + 2\mu - \beta\_1 H) \\ &+ \sqrt{\beta\_1^2 (H - H\_{TB})^2 + \left( \beta\_2 + \beta\_3 - \beta\_4 \right)^2 + 2\beta\_1 (H + H\tau B) (\beta\_2 + \beta\_3 - \beta\_4)} ). \end{aligned}$$

here, λ5, λ<sup>6</sup> have negative real part if (β<sup>2</sup> + β<sup>3</sup> + μ)β1*H* < (β<sup>4</sup> + μ)(β<sup>2</sup> + β<sup>3</sup> + μ + *HTB*β1). Hence, by Routh–Hurwitz criterion, the system is globally stable. Next,

$$\begin{aligned} \mathsf{G}\_{3}(X\_{3}, Z\_{3}) &= \left(\beta\_{5}H\_{\mathsf{TB}} + \beta\_{8}P\_{A} - (\mu + \beta\_{7} + \beta\_{10})\right)P\_{ATB} - \left(\beta\_{5}(H\_{TB} - H\_{TB})P\_{ATB} + \beta\_{8}(P\_{A3} - P\_{A})P\_{ATB}\right) \\ &= \mathsf{C}P\_{ATB} - \hat{\mathsf{G}}\_{3}(X\_{3}, Z\_{3}) \end{aligned}$$

Here, *<sup>G</sup>*<sup>ˆ</sup> <sup>3</sup>(*X*3,*Z*3) <sup>≥</sup> 0, hence the conditions P.5 and P.6 are satisfied. Therefore, by Castillo Chavez the system is globally stable.

#### **Theorem 5.** *The endemic equilibrium point E*∗ (*H*∗ ,*H*∗ *TB*,*P*<sup>∗</sup> *A*,*M*<sup>∗</sup> ,*P*∗ *ATB*, *A*<sup>∗</sup> , *A*∗ *TB*) *is globally asymptotically stable.*

**Proof.** Let us assume Lyapunov function

*L*∗ (*t*) = <sup>1</sup> <sup>2</sup> [(*H* − *H*<sup>∗</sup> )+(*HTB* − *H*<sup>∗</sup> *TB*)+(*PA* − *P*<sup>∗</sup> *<sup>A</sup>*)+(*M* − *M*<sup>∗</sup> )+(*PATB* − *P*<sup>∗</sup> *ATB*)+(*A* − *A*<sup>∗</sup> )+(*ATB* − *A*<sup>∗</sup> *TB*)]<sup>2</sup> *dL*∗ *dt* = [(*H* − *H*<sup>∗</sup> )+(*HTB* − *H*<sup>∗</sup> *TB*)+(*PA* − *P*<sup>∗</sup> *<sup>A</sup>*)+(*M* − *M*<sup>∗</sup> )+(*PATB* − *P*<sup>∗</sup> *ATB*)+(*A* − *A*<sup>∗</sup> )+(*ATB* − *A*<sup>∗</sup> *TB*)] [*H*- + *H*- *TB* + *P*- *<sup>A</sup>* + *M*- + *P*- *ATB* + *A*- + *A*- *ATB*] = [(*H* − *H*<sup>∗</sup> )+(*HTB* − *H*<sup>∗</sup> *TB*)+(*PA* − *P*<sup>∗</sup> *<sup>A</sup>*)+(*M* − *M*<sup>∗</sup> )+(*PATB* − *P*<sup>∗</sup> *ATB*)+(*A* − *A*<sup>∗</sup> )+(*ATB* − *A*<sup>∗</sup> *TB*)] [*B* − μ(*H* + *HTB* + *PA* + *M* + *PATB* + *A* + *ATB*) − μ*DA* − μ*DTBATB*] = −[(*H* − *H*<sup>∗</sup> )+(*HTB* − *H*<sup>∗</sup> *TB*)+(*PA* − *P*<sup>∗</sup> *<sup>A</sup>*)+(*M* − *M*<sup>∗</sup> )+(*PATB* − *P*<sup>∗</sup> *ATB*)+(*A* − *A*<sup>∗</sup> )+(*ATB* − *A*<sup>∗</sup> *TB*)] [μ((*H* − *H*<sup>∗</sup> )+(*HTB* − *H*<sup>∗</sup> *TB*)+(*PA* − *P*<sup>∗</sup> *<sup>A</sup>*)+(*M* − *M*<sup>∗</sup> )+(*PATB* − *P*<sup>∗</sup> *ATB*)+(*A* − *A*<sup>∗</sup> )+(*ATB* − *A*<sup>∗</sup> *TB*))] = −μ[(*H* − *H*<sup>∗</sup> )+(*HTB* − *H*<sup>∗</sup> *TB*)+(*PA* − *P*<sup>∗</sup> *<sup>A</sup>*)+(*M* − *M*<sup>∗</sup> )+(*PATB* − *P*<sup>∗</sup> *ATB*)+(*A* − *A*<sup>∗</sup> )+(*ATB* − *A*<sup>∗</sup> *TB*)]<sup>2</sup> ≤ 0

where *B* = μ(*H*∗ + *H*∗ *TB* + *P*<sup>∗</sup> *<sup>A</sup>* + *M*<sup>∗</sup> + *P*<sup>∗</sup> *ATB* + *A*<sup>∗</sup> + *A*<sup>∗</sup> *TB*) + μ*DA* + μ*DTBATB* -

Here, *dL*<sup>∗</sup> *dt* ≤ 0. Hence, by LaSalle Invariance principle [21] the endemic equilibrium point is globally asymptotically stable.

#### **3. Backward Bifurcation**

If the reproduction number *R*<sup>0</sup> > 1, then *PA* > 0, the system (1) exhibits a unique positive solution *E*∗ . Now, on solving system (2), we have

$$F(P\_A^\*) = b\_2 P\_A^{\*2} + b\_1 P\_A^\* + b\_0 \tag{10}$$

where,

$$\begin{aligned} b\_2 &= \beta\_1 \wp\_8(\beta\_5(L\_2 + \mu) - \wp\_8(\beta\_4 + \mu)) \\ b\_1 &= \beta\_8(\mu + \beta\_4) \left[ \beta\_1 (L\_1 + \mu) + \beta\_5 (L\_3 + \mu) \right] - \beta\_1 \wp\_5[(L\_2 + \mu)(L\_3 + \mu) + B\beta\_8] - \beta\_5^2 (L\_1 + \mu)(L\_2 + \mu) \\ b\_0 &= B\beta\_1 \beta\_5^2 \end{aligned}$$

Here, the coefficient *b*<sup>2</sup> < 0, and *b*<sup>0</sup> depends on the value of *R*0. If *R*<sup>0</sup> < 1, then *b*<sup>0</sup> is positive and if *R*<sup>0</sup> > 1, then *b*<sup>0</sup> is negative. For *R*<sup>0</sup> > 1, Equation (10) has two roots, positive and negative.

For *b*<sup>1</sup> > 0, the system has endemic equilibria continuously depending on *R*0; this shows that there exists an interval for *R*0, which has two positive equilibria as follows:

$$I\_1 = \frac{-b\_1 - \sqrt{b\_1^2 - 4b\_2b\_0}}{2b\_2},\\ I\_2 = \frac{-b\_1 + \sqrt{b\_1^2 - 4b\_2b\_0}}{2b\_2}$$

For Backward Bifurcation, setting *b*<sup>1</sup> <sup>2</sup> <sup>−</sup> <sup>4</sup>*b*2*b*<sup>0</sup> = 0 and solving for critical points of *<sup>R</sup>*<sup>0</sup> gives

$$R\_{\mathbb{C}} = 1 - \frac{\frac{[\beta\_{8}(\mu + \beta\_{4})]\beta\_{1}(L\_{1} + \mu) + \beta\_{5}(L\_{3} + \mu)] - \beta\_{1}\beta\_{5}[(L\_{2} + \mu)(L\_{3} + \mu)}{+B\beta\_{8}] - \beta\_{5}^{2}(L\_{1} + \mu)(L\_{2} + \mu)]^{2}((L\_{1} + \mu)(L\_{2} + \mu)(L\_{3} + \mu) - B\beta\_{1}\beta\_{8})}$$
 
$$R\_{\mathbb{C}} = 1 - \frac{+B\beta\_{8}[-\beta\_{5}^{2}(L\_{1} + \mu)(L\_{2} + \mu)]^{2}((L\_{1} + \mu)(L\_{2} + \mu)(L\_{3} + \mu) - B\beta\_{1}\beta\_{8})}{4(\beta\_{1}\beta s(\beta\varepsilon(L\_{2} + \mu) - \beta s(\beta\_{4} + \mu)))B\beta\_{2}\beta\_{5}^{2}(L\_{1} + \mu)(L\_{2} + \mu)(L\_{3} + \mu)}$$

If *RC* < *R*0**,** then, equivalently, *b*<sup>1</sup> <sup>2</sup> <sup>−</sup> <sup>4</sup>*b*2*b*<sup>0</sup> <sup>&</sup>gt; 0 and backward bifurcation occur for the points of *<sup>R</sup>*0, such that *RC* < *R*<sup>0</sup> < 1 [22], as shown in the above Figure 2. Here, *RC*= 0.95 is the critical value after which co-infection attains stability.

**Figure 2.** Pre-AIDS class at equilibria versus *R*0.

#### *3.1. Sensitivity Analysis of R*<sup>0</sup>

In this section, sensitivity indices of *R*<sup>0</sup> with respect to different parameters are calculated as shown in Table 2, using the formula γ*R*<sup>0</sup> <sup>α</sup> <sup>=</sup> <sup>∂</sup>*R*<sup>0</sup> ∂α . <sup>α</sup> *R*0 , where α is the model parameter. These indices show how crucial each parameter is to disease transmission.



The other parameters β5, β11, β12, μ*D*, μ*DTB* do not have any impact on the sensitivity of reproduction number.

#### *3.2. Numerical Simulation*

From Figure 3 it can be observed that about 34% of the total HIV-infected population gets TB infection within 15 months. Approximately 30% of individuals infected with HIV go for treatment in 27 months. Co-infected individuals undergo treatment for TB in 11 months. Within 26 months, approximately 31% of HIV-infectives proceed to next stage, i.e., AIDS. About 22% of pre-AIDS infectives get TB infection and join pre-AIDS TB in 20 months. Individuals in the pre-AIDS class initially undergo medication then, due to ignorance or any other social reason, individuals leave the compartment and, after some time, joins the medication class again. Between approximately 1.7 and 4 months, individuals in the pre-AIDS class (not taking any kind of medication) suffer from AIDS, whereas those undergoing treatment get infected by AIDS within 28 months. This shows that medication is helpful. Even though it does not help the complete elimination of disease, the rate of disease spread can be controlled.

**Figure 3.** Transmission of HIV–TB co-infection.

From Figure 4, we can conclude that individuals in Pre-AIDS class for a longer duration get AIDS at faster rate than the individuals who have just joined the pre-AIDS class.

**Figure 4.** Intensity of pre-AIDS class versus AIDS class.

Figure 5 indicates that individuals suffering from HIV suffer from TB also, and both the compartments stabilize after some time. Figure 6 shows that individuals in the pre-AIDS class also suffer from TB. The trajectory is stable.

**Figure 5.** Behavior of *H* v/s *H*TB.

**Figure 6.** Trajectory of *PATB* and *PA*.

**Figure 7.** Phase transition plot of *M* with *HTB* and *PATB*.

Figure 8 shows that the newly HIV-TB infected individuals and individuals in Pre-AIDS TB class will oscillate cyclically. Here, neither compartment will die out completely nor they will grow indefinitely.

**Figure 8.** Phase transition of co-infection and pre-stage of co-infection.

From Figure 9, it can be observed that, of the total population, 21% are HIV infected and 13% are HIV-TB infected, whereas the percentage of individuals in pre-AIDS and pre-AIDS TB stage is 16% and 6%, respectively. A total of 15% of the population undergoes treatment for both diseases. Since HIV-AIDS is not curable, even after taking treatment, 17% of cases lead to AIDS infections and 12% are infected by TB, moving towards the AIDS TB class.

**Figure 9.** Percentage wise distribution of different compartments.

#### **4. Conclusions**

In this paper, a mathematical model of HIV-TB considering the HIV-infected population is studied. Using data tabulated in Table 1, we have *R*<sup>0</sup> = 2.262 > 1, which shows the persistence of the disease in the society. HIV-AIDS cannot be eradicated completely from the infected population. Next, global stability is shown for the equilibrium points where there is no co-infection, and instances when there is no individual in the pre-AIDS TB class are shown using Castillo Chavez method. The equilibrium points where there is no co-infection and no individual in the pre-AIDS TB class proved to be globally unstable and is said to exhibit bifurcation. The endemic point is proven to be globally stable using Lyapunov function. Backward bifurcation analysis is studied, which indicates that a minimum of 95% of individuals join the pre-AIDS class. Numerical simulation is done to validate the model, which concludes that the medication plays a vital role in controlling disease spread. Here, we can observe that if treatment is provided at the initial stage of disease, its further progression can be prevented, and survival of individuals can be extended. The value of the reproduction number is highly affected by the rate at which individuals join the AIDS TB class. The pie-chart exhibits the distribution of the population in various compartments in the model.

**Author Contributions:** Conceptualization, N.H.S.; Formal analysis, N.S.; Writing—original draft, Y.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors thank DST-FIST file #MSI-097 for technical support to the department (only technical support to department).

**Acknowledgments:** The paper is prepared under the guidance of Nita H. shah.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
