Global Stability of Multi-Strain SEIR Epidemic Model with Vaccination Strategy

Abstract

A three-strain SEIR epidemic model with a vaccination strategy is suggested and studied in this work. This model is represented by a system of nine nonlinear ordinary differential equations that describe the interaction between susceptible individuals, strain-1-vaccinated individuals, strain-1-exposed individuals, strain-2-exposed individuals, strain-3-exposed individuals, strain-1-infected individuals, strain-2-infected individuals, strain-3-infected individuals, and recovered individuals. We start our analysis of this model by establishing the existence, positivity, and boundedness of all the solutions. In order to show global stability, the model has five equilibrium points: The first one stands for the disease-free equilibrium, the second stands for the strain-1 endemic equilibrium, the third one describes the strain-2 equilibrium, the fourth one represents the strain-3 equilibrium point, and the last one is called the total endemic equilibrium. We establish the global stability of each equilibrium point using some suitable Lyapunov function. This stability depends on the strain-1 reproduction number $R_0^1$, the strain-2 basic reproduction number $R_0^2$, and the strain-3 reproduction number $R_0^3$. Numerical simulations are given to confirm our theoretical results. It is shown that in order to eradicate the infection, the basic reproduction numbers of all the strains must be less than unity.

1. Introduction

Multi-strain models present a very important part in mathematical modeling in order to well understand infectious disease spread. Indeed, many infectious diseases such as human immunodeficiency virus (VIH), tuberculosis, and coronavirus disease (COVID-19) can be analyzed by using different multi-strain epidemic models because these diseases contain usually two or more strains [1,2,3,4,5].

To describe the infection, many works use the classical $SIR$ epidemic model, with S representing the susceptible individuals, I the infected individuals, and R the removed individuals. The first $SIR$ epidemic model was proposed by Kermack and Mc Kendrick in [6]. When infection takes a specific time to appear in infected individuals, another class describing the exposed individuals is added to the $SIR$ epidemic model for a good description of the infection dynamics. These new categories of epidemic models are abbreviated as $SEIR$. Many works have used this model to describe the infection dynamic of infectious diseases [7,8,9,10,11,12,13]. The infection rate of a disease can be defined as the number of newly infected individuals in a specific time [14]. The famous one is a bilinear incidence under the form $\zeta SI$ or $\zeta SI/N$, with zeta as the infection rate and N as the population size. Some mathematical models have used these incidence functions [15,16,17,18,19]. Since mutation is among the characteristics of viruses, a virus can experience several strains. In the case of two strains, the exposed class of the individual for $SEIR$ is divided into two sub-classes, $E_1$ and $E_2$; the first one stands for strain-1-exposed individuals, and $E_2$ stands for strain-2-exposed individuals. The same process applies to the infected population $I_1$ and $I_2$, and they are divided into two sub-populations; the first refers to strain-1-infected individuals, and $I_2$ stands for strain-2-infected individuals. Some multi-strain $SEIR$ epidemic models have used bilinear or non-monotonic incidence rates [20,21,22,23]. Likewise and recently, Bentaleb and Amine [24] suggested a two-strain epidemic model with bilinear and non-monotonic incidence rates. The authors began the analysis of the model by giving the different theorems of existence, positivity, and boundedness of the model solutions, which they demonstrated after the global stability of the equilibrium points, and they gave some numerical simulations in the last part of their work. This last proposed model was improved by Meskaf et al. in [25] by proposing a two-strain epidemic model with non-monotonic incidence rates. More recently, in [26], Yaagoub et al. suggested a two-strain epidemic model with treatment. The authors started their work by proving the existence, positivity, and boundedness of the suggested model solution, and they gave different theorems of the global stability of the equilibria in order to perform numerical simulations for confirming the theoretical results and showing the effect of treatment on infection.

Vaccination is a very effective way to fight most infectious diseases such as COVID-19 [27]. Therefore, developing safe and effective vaccines significantly reduces morbidity and mortality rates. Some mathematical models considered this vaccination strategy in their proposed model [28,29,30,31,32,33,34,35,36]. The $SVEIR$ model is inspired by $SEIR$ models, which take into consideration the vaccination strategy. In the literature, some authors use these $SVEIR$ models to describe the infection transmission of some diseases [37,38,39,40,41,42]. Recently, in [43], Baba et al. suggested a two-strain $SVEIR$ model with a bilinear incidence rate and vaccination strategy. They gave the different theorems of existence, positivity, and boundedness of solutions and also showed the global stability of the equilibria; they finished their work with some numerical simulations and discussions. In this context, and motivated by the previous works, we suggest an $SVEIR$ three-strain epidemic model with a vaccination strategy. More precisely, in our model, we investigate the vaccine effect only on the first strain because several studies have shown that vaccination of only one strain reduces the total infection.

$$\left\{\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =\Lambda -{\alpha }_{1}S{I}_1-{\alpha }_{2}S{I}_2-{\alpha }_{3}S{I}_3-(\sigma +\lambda )S,\hfill & \\ \hfill {\displaystyle \frac{dV}{dt}}& =\sigma S-{\eta }_{2}V{I}_2-{\eta }_{3}V{I}_3-\lambda V,\hfill & \\ \hfill {\displaystyle \frac{d{E}_1}{dt}}& ={\alpha }_{1}S{I}_1-({\gamma }_1+\lambda )E_1,\hfill & \\ \hfill {\displaystyle \frac{d{E}_2}{dt}}& ={\eta }_{2}V{I}_2+{\alpha }_{2}S{I}_2-({\gamma }_2+\lambda )E_2,\hfill & \\ \hfill {\displaystyle \frac{d{E}_3}{dt}}& ={\eta }_{3}V{I}_3+{\alpha }_{3}S{I}_3-({\gamma }_3+\lambda )E_3,\hfill & \\ \hfill {\displaystyle \frac{d{I}_1}{dt}}& ={\gamma }_1{E}_1-({\mu }_1+\delta )I_1,\hfill & \\ \hfill {\displaystyle \frac{d{I}_2}{dt}}& ={\gamma }_2{E}_2-({\mu }_2+\lambda )I_2,\hfill & \\ \hfill {\displaystyle \frac{d{I}_3}{dt}}& ={\gamma }_3{E}_3-({\mu }_3+\lambda )I_3,\hfill & \\ \hfill {\displaystyle \frac{dR}{dt}}& ={\mu }_1{I}_1+{\mu }_2{I}_2+{\mu }_3{I}_3-\lambda R.\hfill \end{array}\right.$$

(1)

In the above, $S,V,E_1,E_2,E_3,I_1,I_2,I_3$, and R represent, respectively, the compartment of susceptible individuals, vaccinated individuals, strain-1-exposed individuals, stain-2-exposed individuals, stain-3-exposed individuals, strain-1-infected individuals, stain-2-infected individuals, stain-3-infected individuals, and recovered individuals. The parameters of this model (1) are given in

Table 1. Description of the parameters in the model.

Parameters	Description
$\Lambda$	The recruitment rate of the population
$\lambda$	The natural mortality rate
${\alpha }_1$	The infection rate of strain-1
${\alpha }_2$	The infection rate of strain-2
${\alpha }_3$	The infection rate of strain-3
$\sigma$	The vaccination rate of the strain-1 individuals
${\eta }_2$	The transmission rate of vaccinated individuals to strain-2
${\eta }_3$	The transmission rate of vaccinated individuals to strain-3
$1/{\gamma }_1$	The average latency period of the strain-1
$1/{\gamma }_2$	The average latency period of the strain-2
$1/{\gamma }_3$	The average latency period of the strain-3
$1/{\mu }_1$	The average infection period of the strain-1
$1/{\mu }_2$	The average infection period of the strain-2
$1/{\mu }_3$	The average infection period of the strain-3

, and the description of all model (1) elements is represented in Figure 1. We assume that there is no reinfection of recovered individuals.

This work is divided as follows: In Section 2, we give some results of the existence, positivity, and boundedness of model (1) solutions. In Section 3, we prove the global stability of equilibrium points. Numerical simulations are presented in Section 4 to validate the different results found in the theoretical part. The last section concludes this work.

2. Existence, Positivity, and Boundedness of Solutions

In this section, we prove that model (1) has a unique, non-negative, and bounded solution for all $t\ge 0$.

Proposition 1.

For any non-negative initial condition, model (1) has a unique solution. In addition, this solution remains non-negative and bounded for all $t\ge 0$.

Proof. 

First, we prove that system (1) has a unique solution. We can reformulate the model (1) as follows:

$$\dot{U}=f\left(U\right),$$

(2)

with

$$U={(S,V,E_1,E_2,E_3,I_1,I_2,I_3,R)}^T$$

(3)

and

$$F\left(U\right)=\left(\begin{array}{c}\Lambda -{\alpha }_{1}S{I}_1-{\alpha }_{2}S{I}_2-{\alpha }_{3}S{I}_3+(\sigma +\lambda )S\\ \sigma S-{\eta }_{2}V{I}_2-{\eta }_{3}V{I}_3+\lambda V\\ {\alpha }_{1}S{I}_1-({\gamma }_1+\lambda )E_1\\ {\eta }_{2}V{I}_2+{\alpha }_{2}S{I}_2-({\gamma }_2+\lambda )E_2\\ {\eta }_{3}V{I}_3+{\alpha }_{3}S{I}_3-({\gamma }_3+\lambda )E_3\\ {\gamma }_1{E}_1-({\mu }_1+\delta )I_1\\ {\gamma }_2{E}_2-({\mu }_2+\lambda )I_2\\ {\gamma }_3{E}_3-({\mu }_3+\lambda )I_3\\ {\mu }_1{I}_1+{\mu }_2{I}_2+{\mu }_3{I}_3-\lambda R\end{array}\right).$$

(4)

We remark that F is a Lipschitz function; moreover, we have

$$\Vert F\left(U_1\right)-F\left(U_2\right)\Vert \le k{\Vert U_1-U_2\Vert }_1,\forall U_1,U_2\in {\mathbb{R}}_{+}^9$$

(5)

with

$$k=max\{A,B,C,D,E,G,H,L,M\}$$

(6)

and $A={\alpha }_1{I}_1+{\alpha }_2{I}_2+{\alpha }_3{I}_3+(\sigma +\lambda ),B=\sigma +\lambda ,C={\gamma }_1+\lambda ,D={\gamma }_2+\lambda ,E={\gamma }_3+\lambda ,G={\mu }_1+\lambda ,H={\mu }_2+\lambda ,L={\mu }_3+\lambda ,M=\lambda .$

Thus, model (1) has a unique solution in ${\mathbb{R}}_{+}^9$.

Now, we prove that this solution remains non-negative.

$$\left\{\begin{array}{ccc}\hfill \dot{S}{\mid }_{S=0}& =\lambda >0,\hfill & \\ \hfill \dot{V}{\mid }_{V=0}& =\sigma S\ge 0,\hfill & \\ \hfill \dot{E_1}{\mid }_{E_1=0}& ={\alpha }_{1}S{I}_1\ge 0,\hfill & \\ \hfill \dot{E_2}{\mid }_{E_2=0}& ={\eta }_{2}V{I}_2+{\alpha }_{2}S{I}_2\ge 0,\hfill & \\ \hfill \dot{E_3}{\mid }_{E_3=0}& ={\eta }_{3}V{I}_3+{\alpha }_{3}2S{I}_3\ge 0,\hfill & \\ \hfill \dot{I_1}{\mid }_{I_1=0}& ={\gamma }_1{E}_1\ge 0,\hfill & \\ \hfill \dot{I_2}{\mid }_{I_2=0}& ={\gamma }_2{E}_2\ge 0,\hfill & \\ \hfill \dot{I_3}{\mid }_{I_3=0}& ={\gamma }_3{E}_3\ge 0,\hfill & \\ \hfill \dot{R}{\mid }_{R=0}& ={\mu }_1{I}_1+{\mu }_2{I}_2+{\mu }_3{I}_3\ge 0.\hfill & \end{array}\right.$$

(7)

Thus, this solution remains non-negative for all $t\ge 0$.

Finally, for the boundedness of this solution, we verify that the biologically feasible region

$$P=\{(S,V,E_1,E_2,E_3,I_1,I_2,I_3,R)\in {\mathbb{R}}_{+}^{9}suchthatS+V+E_1+E_2+E_3+I_1+I_2+I_3+R\le {\displaystyle \frac{\Lambda }{\lambda }}\}$$

(8)

is positively invariant.

Let the total population

$$N\left(t\right)=S\left(t\right)+V\left(t\right)+E_1\left(t\right)+E_2\left(t\right)+E_3\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+I_3\left(t\right)+R\left(t\right),$$

(9)

By adding all equations of the system (1), we have

$$\dot{N}=\Lambda -\lambda N\left(t\right),$$

(10)

then,

$$N\left(t\right)\le {\displaystyle \frac{\Lambda }{\lambda }}\left(N\left(0\right)-{\displaystyle \frac{\Lambda }{\lambda }}\right)e^{-\lambda t},$$

(11)

therefore,

$$li{m}_{\begin{array}{c}t\to +\infty \end{array}}N\left(t\right)={\displaystyle \frac{\Lambda }{\lambda }},$$

(12)

since $\exists t_0>0$ such as $\forall t\ge t_0$, we will have

$$N\left(t\right)\le {\displaystyle \frac{\Lambda }{\lambda }}.$$

(13)

Then, we conclude that P is positively invariant. Therefore, we can conclude that model (1) has a unique, positive, and bounded solution $\in {\mathbb{R}}_{+}^9$. □

3. Analysis of the Model

In this section, we prove that there exists a disease-free equilibrium point and four equilibrium points. The global stability of these equilibrium points using the Lyapunov functional method is proved. As the first eight equations of system (1) are independent of the ninth equation, and also the total number of population N is determined using Equation (15), system (1) can be reduced to the following system:

$$\left\{\begin{array}{ccc}\hfill {\displaystyle \frac{ds}{dt}}& =\Lambda -{\alpha }_{1}S{I}_1-{\alpha }_{2}S{I}_2-{\alpha }_{3}S{I}_3-(\sigma +\lambda )S,\hfill & \\ \hfill {\displaystyle \frac{dV}{dt}}& =\sigma S-{\eta }_{2}V{I}_2-{\eta }_{3}V{I}_3-\lambda V,\hfill & \\ \hfill {\displaystyle \frac{d{E}_1}{dt}}& ={\alpha }_{1}S{I}_1-({\gamma }_1+\lambda )E_1,\hfill & \\ \hfill {\displaystyle \frac{d{E}_2}{dt}}& ={\eta }_{2}V{I}_2+{\alpha }_{2}S{I}_2-({\gamma }_2+\lambda )E_2,\hfill & \\ \hfill {\displaystyle \frac{d{E}_3}{dt}}& ={\eta }_{3}V{I}_3+{\alpha }_{3}S{I}_3-({\gamma }_3+\lambda )E_3,\hfill & \\ \hfill {\displaystyle \frac{d{I}_1}{dt}}& ={\gamma }_1{E}_1-({\mu }_1+\delta )I_1,\hfill & \\ \hfill {\displaystyle \frac{d{I}_2}{dt}}& ={\gamma }_2{E}_2-({\mu }_2+\lambda )I_2,\hfill & \\ \hfill {\displaystyle \frac{d{I}_3}{dt}}& ={\gamma }_3{E}_3-({\mu }_3+\lambda )I_3.\hfill \end{array}\right.$$

(14)

with

$$R=N-S-V-E_1-E_2-E_3-I_1-I_2-I_3.$$

(15)

3.1. The Basic Reproduction Number Calculation

The basic reproduction number $R_0$ is the number of secondary infection cases caused by one infected individual in a population constituted only by susceptible individuals, and mathematically, the basic reproduction number is the spectral radius of the matrix called the next-generation matrix $F{V}^{-1}$, with F as the positive matrix of new infection cases and V as the matrix of the transition of the infections.

Let

$$a=\sigma +\lambda ,b={\gamma }_1+\lambda ,c={\gamma }_2+\lambda ,d={\gamma }_3+\lambda ,e={\mu }_1+\lambda ,f={\mu }_2+\lambda ,g={\mu }_3+\lambda ,$$

(16)

$$F=\left(\begin{array}{cccccc}0& 0& 0& {\displaystyle \frac{{\alpha }_1\sigma }{a}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{{\alpha }_1\sigma }{a}}+{\displaystyle \frac{{\eta }_2\sigma \Lambda }{\lambda a}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{{\alpha }_3\sigma \Lambda }{a}}+{\displaystyle \frac{{\eta }_3\sigma \Lambda }{\lambda a}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$$

(17)

and

$$V=\left(\begin{array}{cccccc}b& 0& 0& 0& 0& 0\\ 0& c& 0& 0& 0& 0\\ 0& 0& d& 0& 0& 0\\ -{\gamma }_1& 0& 0& e& 0& 0\\ 0& -{\gamma }_2& 0& 0& f& 0\\ 0& 0& -{\gamma }_3& 0& 0& g\end{array}\right).$$

(18)

Thus, we have

$$F{V}^{-1}=\left(\begin{array}{cccccc}{\displaystyle \frac{{\alpha }_1\Lambda {\gamma }_1}{abe}}& 0& 0& {\displaystyle \frac{{\alpha }_1\sigma }{ae}}& 0& 0\\ 0& {\displaystyle \frac{{\alpha }_2\Lambda {\gamma }_2}{acf}}+{\displaystyle \frac{{\eta }_2\Lambda {\gamma }_2}{\lambda acf}}& 0& 0& {\displaystyle \frac{{\alpha }_2\lambda }{af}}+{\displaystyle \frac{{\eta }_2\sigma \Lambda }{\lambda af}}& 0\\ 0& 0& {\displaystyle \frac{{\alpha }_3\Lambda {\gamma }_3}{adg}}+{\displaystyle \frac{{\eta }_3\Lambda {\gamma }_3}{\lambda adg}}& 0& 0& {\displaystyle \frac{{\alpha }_3\sigma \Lambda }{ag}}+{\displaystyle \frac{{\alpha }_3\sigma \Lambda }{\lambda ag}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right).$$

(19)

Therefore, the basic reproduction number of model (14) is

$$R_0=max\{R_0^1,R_0^2,R_0^3\},$$

(20)

with

$$R_0^1={\displaystyle \frac{{\alpha }_1\Lambda {\gamma }_1}{abe}},$$

(21)

$$R_0^2={\displaystyle \frac{{\alpha }_2\Lambda {\gamma }_2}{acf}}+{\displaystyle \frac{{\eta }_2\Lambda {\gamma }_2}{\lambda acf}}$$

(22)

and

$$R_0^3={\displaystyle \frac{{\alpha }_3\sigma \Lambda }{ag}}+{\displaystyle \frac{{\alpha }_3\sigma \Lambda }{\lambda ag}}.$$

(23)

3.2. Steady States

Model (14) has one disease-free equilibrium point, and the other four endemic equilibrium points are given by

• The disease-free equilibrium ${\mathcal{E}}_0=(S_0,V_0,0,0,0,0,0,0)$, where

  $$S_0={\displaystyle \frac{\Lambda }{\sigma +\lambda }},V_0={\displaystyle \frac{\sigma \Lambda }{\lambda a}}.$$

  (24)
• The strain-1 endemic equilibrium ${\mathcal{E}}_{s_1}=(S_{s_1}^{*},V_{s_1}^{*},E_{1,s_1}^{*},E_{2,s_1}^{*},E_{3,s_1}^{*},I_{1,s_1}^{*},I_{2,s_1}^{*},I_{3,s_1}^{*})$, where

  $$S_{s_1}^{*}={\displaystyle \frac{be}{{\alpha }_1{\gamma }_1}},V_{s_1}^{*}={\displaystyle \frac{\sigma }{\lambda }}{S}_{s_1}^{*}={\displaystyle \frac{\sigma be}{\lambda {\alpha }_1{\gamma }_1}},I_{1,s_1}^{*}={\displaystyle \frac{a}{{\alpha }_1}}(R_0^1-1),$$

  (25)

  $$E_{1,s_1}^{*}=e{I}_{1,s_1}^{*}={\displaystyle \frac{ea}{{\gamma }_1{\alpha }_1}}(R_0^1-1),E_{2,s_1}^{*}=I_{2,s_1}^{*}=E_{3,s_1}^{*}=I_{3,s_1}^{*}=0.$$

  (26)
• The strain-2 endemic equilibrium ${\mathcal{E}}_{s_2}=(S_{s_2}^{*},V_{s_2}^{*},E_{1,s_2}^{*},E_{2,s_2}^{*},E_{3,s_2}^{*},I_{1,s_2}^{*},I_{2,s_2}^{*},I_{3,s_2}^{*})$, where

  $$S_{s_2}^{*}={\displaystyle \frac{\Lambda }{{\alpha }_2{I}_{2,s_2}^{*}+a}},V_{s_2}^{*}={\displaystyle \frac{\sigma \Lambda }{({\alpha }_2{I}_{2,s_2}^{*}+a)({\alpha }_2{I}_{2,s_2}^{*}+\lambda )}},E_{2,s_2}^{*}={\displaystyle \frac{f}{{\gamma }_2}}{I}_{2,s_2}^{*},$$

  (27)

  $$I_{2,s_2}^{*}={\displaystyle \frac{-B_2-\sqrt{B_2^2-4A_2{C}_2}}{2A_2}},E_{1,s_2}^{*}=0,E_{3,s_2}^{*}=0,I_{1,s_2}^{*}=0,I_{3,s_2}^{*}=0.$$

  (28)

  where $I_{2,s_2}^{*}$ is solution of the equation $A_2^2{x}^2+B_{2}x+C_2=0,$ with $A_2={\alpha }_2^2\Lambda -{\displaystyle \frac{fc{\alpha }_2^2}{{\gamma }_2}},B_2={\eta }_2\sigma \Lambda +{\alpha }_2\Lambda \lambda -{\displaystyle \frac{fc}{{\gamma }_2({\alpha }_2\lambda +{\alpha }_{2}a)}},C_2=\lambda a$.
• The strain-3 endemic equilibrium ${\mathcal{E}}_{s_3}=(S_{s_3}^{*},V_{s_3}^{*},E_{1,s_3}^{*},E_{2,s_3}^{*},E_{3,s_3}^{*},I_{1,s_3}^{*},I_{2,s_3}^{*},I_{3,s_3}^{*})$, where

  $$S_{s_3}^{*}={\displaystyle \frac{\Lambda }{{\alpha }_3{I}_{3,s_3}^{*}+a}},V_{s_3}^{*}={\displaystyle \frac{\sigma \Lambda }{({\alpha }_3{I}_{3,s_3}^{*}+a)({\alpha }_3{I}_{3,s_3}^{*}+\lambda )}},E_{3,s_3}^{*}={\displaystyle \frac{g}{{\gamma }_3}}{I}_{3,s_3}^{*},$$

  (29)

  $$I_{3,s_3}^{*}={\displaystyle \frac{-B_3-\sqrt{B_3^2-4A_3{C}_3}}{2A_3}},E_{1,s_3}^{*}=0,E_{2,s_3}^{*}=0,I_{1,s_3}^{*}=0,I_{2,s_3}^{*}=0.$$

  (30)

  where $I_{3,s_3}^{*}$ is solution of the equation $A_3^2{x}^2+B_{3}x+C_3=0,$ with $A_3={\alpha }_3^2\Lambda -{\displaystyle \frac{dg{\alpha }_3^2}{{\gamma }_3}},B_3={\eta }_3\sigma \Lambda +{\alpha }_3\Lambda \lambda -{\displaystyle \frac{dg}{{\gamma }_3({\alpha }_3\lambda +{\alpha }_{3}a)}},C_3=\lambda a$.
• The total strain endemic equilibrium ${\mathcal{E}}_t=(S_t^{*},V_t^{*},E_{1,t}^{*},E_{2,t}^{*},E_{3,t}^{*},I_{1,t}^{*},I_{2,t}^{*},I_{3,s_t}^{*})$, where

  $$S_t^{*}={\displaystyle \frac{be}{{\alpha }_1{\gamma }_1}},V_t^{*}={\displaystyle \frac{{\alpha }_{2}be}{{\eta }_2{\alpha }_1{\gamma }_1}}\left({\displaystyle \frac{R_0^1}{R_0^2}}-1\right)={\displaystyle \frac{{\alpha }_{2}be}{{\eta }_3{\alpha }_1{\gamma }_3}}\left({\displaystyle \frac{R_0^1}{R_0^3}}-1\right),$$

  (31)

  $$E_{1,t}^{*}={\displaystyle \frac{e}{{\gamma }_1}}{I}_{1,t}^{*},E_{2,t}^{*}={\displaystyle \frac{f}{{\gamma }_2}}{E}_{2,t}^{*},E_{3,t}^{*}={\displaystyle \frac{g}{{\gamma }_3}}{I}_{3,t}^{*}.$$

  (32)

  where $I_{1,t}^{*},I_{2,t}^{*}$ and $I_{3,t}^{*}$ are the roots of the following equations:

  $$a_{1}x+b_{1}y+c_{1}z+d_1=0,$$

  (33)

  $$a_{2}y+b_{2}z+d_2=0$$

  (34)

  and

  $$a_{3}x+b_{3}y+c_{3}z+d_3=0.$$

  (35)

  where $a_1={\alpha }_1{S}_t^{*},b_1={\alpha }_2{S}_t^{*},c_1={\alpha }_3{S}_t^{*},d_1=-\lambda -a{S}_t^{*},a_2={\eta }_2{V}_t^{*},b_2={\eta }_3{V}_t^{*},d_2=\sigma S_t^{*}-\lambda V_t^{*},a_3={\displaystyle \frac{be}{{\gamma }_1}},b_3={\displaystyle \frac{fc}{{\gamma }_2}},c_3={\displaystyle \frac{dg}{{\gamma }_3}},d_3=a{S}_t^{*}+\lambda V_t^{*}-\Lambda -\sigma S_t^{*}.$

Remark 1.

From the components of the equilibrium points, we conclude that points exist when$R_0^1>1,R_0^2>1$, and$R_0^3>1$.

3.3. Global Stability

In this section, we give the different theorems concerning the global stability of the different equilibrium points:

Theorem 1.

If $R_0^1\le 1,R_0^2\le 1$, and $R_0^3\le 1$, then the disease-free equilibrium ${\mathcal{E}}_f$ is globally asymptotically stable.

Proof. 

We consider the following Lyapunov function $L_f$ in ${\mathbb{R}}_{+}^8$:

$$L_f=S_0\left({\displaystyle \frac{S}{S_0}}-ln\left({\displaystyle \frac{S}{S_0}}\right)-1\right)+V_0\left({\displaystyle \frac{V}{V_0}}-ln\left({\displaystyle \frac{V}{V_0}}\right)-1\right)+E_1+E_2+E_3+{\displaystyle \frac{b}{{\gamma }_1}}{I}_1+{\displaystyle \frac{c}{{\gamma }_2}}{I}_2+{\displaystyle \frac{d}{{\gamma }_3}}{I}_3.$$

(36)

The time derivative of $L_f$ is given by

$$\begin{array}{cc}\hfill {\dot{L}}_f& =\dot{S}+\dot{V}+{\dot{E}}_1+{\dot{E}}_2+{\dot{E}}_3+{\displaystyle \frac{b}{{\gamma }_1}}{\dot{I}}_1+{\displaystyle \frac{c}{{\gamma }_2}}{\dot{I}}_2+{\displaystyle \frac{d}{{\gamma }_3}}{\dot{I}}_3-{\displaystyle \frac{S_0}{S}}\dot{S}-{\displaystyle \frac{V_0}{V}}\dot{V}\hfill \\ & +\Lambda -\lambda S-\lambda V-{\displaystyle \frac{be}{{\gamma }_1}}{I}_1-{\displaystyle \frac{fc}{{\gamma }_2}}{I}_2-{\displaystyle \frac{dg}{{\gamma }_3}}{I}_3-\Lambda {\displaystyle \frac{S_0}{S}}+\alpha S_0{I}_1\hfill \\ & +{\alpha }_2{S}_0{I}_2+{\alpha }_3{S}_0{I}_3+\sigma S_0+\lambda S_0-\sigma S{\displaystyle \frac{V_0}{V}}+{\eta }_2{V}_0{I}_2+{\eta }_3{V}_0{I}_3+\lambda V_0.\hfill \end{array}$$

(37)

As $S_0={\displaystyle \frac{\Lambda }{\sigma +\lambda }}$ and $V_0={\displaystyle \frac{\sigma \Lambda }{\lambda (\sigma +\lambda )}}$, we will have

$$\begin{array}{cc}\hfill {\dot{L}}_f& =\lambda S_0\left(2-{\displaystyle \frac{S_0}{S}}-{\displaystyle \frac{S}{S_0}}\right)+\sigma S_0\left(3-{\displaystyle \frac{S_0}{S}}-{\displaystyle \frac{S{V}_0}{S_{0}V}}-{\displaystyle \frac{V}{V_0}}\right)+{\displaystyle \frac{be}{{\gamma }_1}}{I}_1\left(R_0^1-1\right)\hfill \\ & +{\displaystyle \frac{fc}{{\gamma }_2}}{I}_2\left(R_0^2-1\right)+{\displaystyle \frac{dg}{{\gamma }_3}}{I}_3\left(R_0^3-1\right).\hfill \end{array}$$

(38)

Thus, when $R_0^1\le 1,R_0^2\le 1$ and $R_0^3\le 1$, we will have ${\dot{L}}_f\le 0$. Thus, the disease-free equilibrium point ${\mathcal{E}}_0$ is globally asymptotically stable. □

Theorem 2.

If $R_0^2\le 1,R_0^3\le 1$, and $R_0^1>1$, then the strain-1 endemic equilibrium point ${\mathcal{E}}_{s_1}$ is globally asymptotically stable.

Proof. 

We consider the following Lyapunov function $L_1$ in ${\mathbb{R}}_{+}^8$:

$$\begin{array}{cc}\hfill L_1& =S_{s_1}^{*}\left({\displaystyle \frac{S}{S_{s_1}^{*}}}-ln\left({\displaystyle \frac{S}{S_{s_1}^{*}}}\right)-1\right)+V_{s_1}^{*}\left({\displaystyle \frac{V}{V_{s_1}^{*}}}-ln\left({\displaystyle \frac{V}{V_{s_1}^{*}}}\right)-1\right)+E_{1,s_1}^{*}\left({\displaystyle \frac{E_1}{E_{1,s_1}^{*}}}-ln\left({\displaystyle \frac{E_1}{E_{1,s_1}^{*}}}\right)-1\right)\hfill \\ & +E_2+E_3+{\displaystyle \frac{b}{{\gamma }_1}}{I}_{1,s_1}^{*}\left({\displaystyle \frac{I_1}{I_{1,s_1}^{*}}}-ln\left({\displaystyle \frac{I_1}{I_{1,s_1}^{*}}}\right)-1\right)+{\displaystyle \frac{c}{{\gamma }_2}}{I}_2+{\displaystyle \frac{d}{{\gamma }_3}}{I}_3.\hfill \end{array}$$

(39)

The time derivative of $L_1$ is given by

$$\begin{array}{cc}\hfill {\dot{L}}_1& =\left(1-{\displaystyle \frac{S_{s_1}^{*}}{S}}\right)\dot{S}+\left(1-{\displaystyle \frac{V_{s_1}^{*}}{S}}\right)\dot{V}+\left(1-{\displaystyle \frac{E_{1,s_1}^{*}}{E_1}}\right)\dot{E_1}+\dot{E_2}+\dot{E_3}+{\displaystyle \frac{b}{{\gamma }_1}}\left(1-{\displaystyle \frac{I_{1,s_1}^{*}}{I_1}}\right)\dot{I_1}+{\displaystyle \frac{c}{{\gamma }_2}}\dot{I_2}+{\displaystyle \frac{d}{{\gamma }_3}}\dot{I_3}\hfill \\ & =\Lambda -\lambda S-\lambda V-{\displaystyle \frac{be}{{\gamma }_1}}-{\displaystyle \frac{fc}{{\gamma }_2}}-{\displaystyle \frac{dg}{{\gamma }_3}}{I}_3-\Lambda {\displaystyle \frac{S_{s_1}^{*}}{S}}+{\alpha }_1{S}_{s_1}^{*}{I}_1+{\alpha }_2{S}_{s_1}^{*}{I}_2+{\alpha }_3{S}_{s_1}^{*}{I}_3+(\sigma +\lambda )S_{s_1}^{*}\hfill \\ & -\sigma S{\displaystyle \frac{V_{s_1}^{*}}{V}}+{\eta }_2{V}_{s_1}^{*}{I}_2+{\eta }_3{V}_{s_1}^{*}{I}_3+\lambda V_{s_1}^{*}-{\alpha }_{1}S{I}_1{\displaystyle \frac{E_{1,s_1}^{*}}{E_1}}+b{E}_{1,s_1}^{*}-b{E}_1{\displaystyle \frac{I_{1,s_1}^{*}}{I_1}}+{\displaystyle \frac{be}{{\gamma }_1}}{I}_{1,s_1}^{*}.\hfill \end{array}$$

(40)

As $E_1$ is an equilibrium point of system (14), we will have

$$\left\{\begin{array}{ccc}\hfill \Lambda & ={\alpha }_1{S}_{s_1}^{*}{I}_{1,s_1}^{*}+(\sigma +\lambda )S_{s_1}^{*},\hfill & \\ \hfill \sigma S_{s_1}^{*}& =\lambda V_{s_1}^{*},\hfill & \\ \hfill b{E}_{1,s_1}^{*}& ={\alpha }_1{S}_{s_1}^{*}{I}_{1,s_1}^{*},\hfill & \\ \hfill E_{1,s_1}^{*}& ={\displaystyle \frac{e}{{\gamma }_1}}{I}_{1,s_1}^{*}.\hfill & \end{array}\right.$$

(41)

Thus, after some simplifications and factorizations, we will have

$$\begin{array}{cc}\hfill {\dot{L}}_1& =\sigma S_{s_1}^{*}\left(3-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_1}^{*}}{S_{s_1}^{*}V}}-{\displaystyle \frac{V}{V_{s_1}^{*}}}\right)+\lambda S_{s_1}^{*}\left(2-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_1}^{*}}}\right)\hfill \\ & +b{E}_{1,s_1}^{*}\left(3-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{E_1{I}_{1,s_1}^{*}}{E_{1,s_1}^{*}{I}_1}}-{\displaystyle \frac{S{I}_1{E}_{1,s_1}^{*}}{S_{s_1}^{*}{I}_{1,s_1}^{*}{E}_1}}\right)+I_1\left({\alpha }_1{S}_{s_1}^{*}-{\displaystyle \frac{be}{{\gamma }_1}}\right)\hfill \\ & +I_2\left({\alpha }_2{S}_{s_1}^{*}+{\eta }_2{V}_{s_1}^{*}-{\displaystyle \frac{fc}{{\gamma }_2}}\right)+I_3\left({\alpha }_3{S}_{s_1}^{*}+{\eta }_3{V}_{s_1}^{*}-{\displaystyle \frac{dg}{{\gamma }_3}}\right)\hfill \end{array}$$

(42)

$$\begin{array}{cc}& \le \sigma S_{s_1}^{*}\left(3-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_1}^{*}}{S_{s_1}^{*}V}}-{\displaystyle \frac{V}{V_{s_1}^{*}}}\right)+\lambda S_{s_1}^{*}\left(2-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_1}^{*}}}\right)\hfill \\ & +b{E}_{1,s_1}^{*}\left(3-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{E_1{I}_{1,s_1}^{*}}{E_{1,s_1}^{*}{I}_1}}-{\displaystyle \frac{S{I}_1{E}_{1,s_1}^{*}}{S_{s_1}^{*}{I}_{1,s_1}^{*}{E}_1}}\right)+I_1\left(\alpha S_{s_1}^{*}-{\displaystyle \frac{be}{{\gamma }_1}}\right)\hfill \\ & +{\displaystyle \frac{fc}{{\gamma }_2}}\left(R_0^2-1\right)+{\displaystyle \frac{dg}{{\gamma }_3}}\left(R_0^3-1\right).\hfill \end{array}$$

(43)

As the arithmetic mean is greater than or equal to the geometric mean, we will have

$$3-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_1}^{*}}{S_{s_1}^{*}V}}-{\displaystyle \frac{V}{V_{s_1}^{*}}}\le 0,$$

(44)

$$3-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{E_1{I}_{1,s_1}^{*}}{E_{1,s_1}^{*}{I}_1}}-{\displaystyle \frac{S{I}_1{E}_{1,s_1}^{*}}{S_{s_1}^{*}{I}_{1,s_1}^{*}{E}_1}}\le 0$$

(45)

and

$$2-{\displaystyle \frac{S_{s_1}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_1}^{*}}}\le 0.$$

(46)

Moreover,

$$S_{s_1}^{*}-{\displaystyle \frac{be}{{\gamma }_1}}=0.$$

(47)

Then, when $R_0^2\le 1$ and $R_0^3\le 1$, we will have ${\dot{L}}_1\le 0$. Thus, the strain-1 endemic equilibrium point ${\mathcal{E}}_{s_1}$ is globally asymptotically stable. □

For the global stability of the equilibrium point ${\mathcal{E}}_{s_2}$, we assume that this point verifies the following condition:

$${\alpha }_2{S}_{s_2}^{*}+{\eta }_2{V}_{s_2}^{*}-{\displaystyle \frac{fc}{{\gamma }_2}}\le 0.\left(H_1\right)$$

(48)

Theorem 3.

If $R_0^1\le 1,R_0^3\le 1$, and $R_0^2>1$, then the strain-2 endemic equilibrium point ${\mathcal{E}}_{s_2}$ is globally asymptotically stable.

Proof. 

We consider the following Lyapunov function $L_2$ in ${\mathbb{R}}_{+}^8$:

$$\begin{array}{cc}\hfill L_2& =S_{s_2}^{*}\left({\displaystyle \frac{S}{S_{s_2}^{*}}}-ln\left({\displaystyle \frac{S}{S_{s_2}^{*}}}\right)-1\right)+V_{s_2}^{*}\left({\displaystyle \frac{V}{V_{s_2}^{*}}}-ln\left({\displaystyle \frac{V}{V_{s_2}^{*}}}\right)-1\right)+E_1\hfill \\ & +E_{2,s_2}^{*}\left({\displaystyle \frac{E_2}{E_{2,s_2}^{*}}}-ln\left({\displaystyle \frac{E_2}{E_{2,s_2}^{*}}}\right)-1\right)+E_3+{\displaystyle \frac{b}{{\gamma }_1}}{I}_1+{\displaystyle \frac{c}{{\gamma }_2}}{I}_{2,s_2}^{*}\left({\displaystyle \frac{I_2}{I_{2,s_2}^{*}}}-ln\left({\displaystyle \frac{I_2}{I_{2,s_2}^{*}}}\right)-1\right)+{\displaystyle \frac{d}{{\gamma }_3}}{I}_3.\hfill \end{array}$$

(49)

The time derivative of $L_2$ is given by

$$\begin{array}{cc}\hfill {\dot{L}}_2& =\left(1-{\displaystyle \frac{S_{s_2}^{*}}{S}}\right)\dot{S}+\left(1-{\displaystyle \frac{V_{s_2}^{*}}{S}}\right)\dot{V}+\dot{E_1}+\left(1-{\displaystyle \frac{E_{2,s_2}^{*}}{E_2}}\right)\dot{E_2}+\dot{E_3}+{\displaystyle \frac{b}{{\gamma }_1}}\dot{I_1}+{\displaystyle \frac{c}{{\gamma }_2}}\left(1-{\displaystyle \frac{I_{2,s_2}^{*}}{I_2}}\right)\dot{I_2}\hfill \\ & +{\displaystyle \frac{c}{{\gamma }_2}}\dot{I_2}+{\displaystyle \frac{d}{{\gamma }_3}}\dot{I_3}\hfill \end{array}$$

(50)

$$\begin{array}{cc}& =\Lambda -\lambda S-\lambda V-{\displaystyle \frac{be}{{\gamma }_1}}-{\displaystyle \frac{fc}{{\gamma }_2}}-{\displaystyle \frac{dg}{{\gamma }_3}}{I}_3-\Lambda {\displaystyle \frac{S_{s_2}^{*}}{S}}+{\alpha }_1{S}_{s_2}^{*}{I}_1+{\alpha }_2{S}_{s_2}^{*}{I}_2+{\alpha }_3{S}_{s_2}^{*}{I}_3+(\sigma +\lambda )S_{s_2}^{*}\hfill \\ & -\sigma S{\displaystyle \frac{V_{s_2}^{*}}{V}}+{\eta }_2{V}_{s_2}^{*}{I}_2+{\eta }_3{V}_{s_2}^{*}{I}_3+\lambda V_{s_2}^{*}-{\alpha }_{2}S{I}_2{\displaystyle \frac{E_{2,s_2}^{*}}{E_2}}-{\eta }_{2}V{I}_2{\displaystyle \frac{E_{2,s_2}^{*}}{E_2}}+c{E}_{2,s_2}^{*}-c{E}_2{\displaystyle \frac{I_{2,s_2}^{*}}{I_2}}+{\displaystyle \frac{fc}{{\gamma }_2}}{I}_{2,s_2}^{*}.\hfill \end{array}$$

(51)

As ${\mathcal{E}}_{s_2}$ is an equilibrium point of system (14), we will have

$$\left\{\begin{array}{ccc}\hfill \Lambda & ={\alpha }_2{S}_{s_2}^{*}{I}_{2,s_2}^{*}+(\sigma +\lambda )S_{s_2}^{*},\hfill & \\ \hfill \sigma S_{s_2}^{*}& =\lambda V_{s_2}^{*}+{\eta }_2{V}_{s_2}^{*}{I}_{2,s_2}^{*},\hfill & \\ \hfill c{E}_{2,s_2}^{*}& ={\alpha }_2{S}_{s_2}^{*}{I}_{2,s_2}^{*}+{\eta }_2{V}_{s_2}^{*}{I}_{2,s_2}^{*},\hfill & \\ \hfill E_{2,s_2}^{*}& ={\displaystyle \frac{f}{{\gamma }_2}}{I}_{1,s_2}^{*}.\hfill & \end{array}\right.$$

(52)

Thus, after some simplifications and factorizations, we will have

$$\begin{array}{cc}\hfill {\dot{L}}_2& =\sigma S_{s_2}^{*}\left(3-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_2}^{*}}{S_{s_2}^{*}V}}-{\displaystyle \frac{V}{V_{s_2}^{*}}}\right)+\lambda S_{s_2}^{*}\left(2-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_2}^{*}}}\right)\hfill \\ & +(c{E}_{2,s_2}^{*}-{\eta }_2{V}_{s_2}^{*}{I}_{2,s_2}^{*})\left(3-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{E_2{I}_{2,s_2}^{*}}{E_{2,s_2}^{*}{I}_2}}-{\displaystyle \frac{S{I}_2{E}_{2,s_2}^{*}}{S_{s_2}^{*}{I}_{2,s_2}^{*}{E}_2}}\right)\hfill \\ & +I_1\left({\alpha }_1{S}_{s_2}^{*}-{\displaystyle \frac{be}{{\gamma }_1}}\right)+I_2\left({\alpha }_2{S}_{s_2}^{*}+{\eta }_2{V}_{s_2}^{*}-{\displaystyle \frac{fc}{{\gamma }_2}}\right)+I_3\left({\alpha }_3{S}_{s_2}^{*}+{\eta }_3{V}_{s_2}^{*}-{\displaystyle \frac{dg}{{\gamma }_3}}\right)-{\eta }_{2}V{I}_2{\displaystyle \frac{E_{2,s_2}^{*}}{E_2}}\hfill \end{array}$$

(53)

$$\begin{array}{cc}& \le \sigma S_{s_2}^{*}\left(3-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_2}^{*}}{S_{s_2}^{*}V}}-{\displaystyle \frac{V}{V_{s_2}^{*}}}\right)+\lambda S_{s_2}^{*}\left(2-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_2}^{*}}}\right)\hfill \\ & +(c{E}_{2,s_2}^{*}-{\eta }_2{V}_{s_2}^{*}{I}_{2,s_2}^{*})\left(3-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{E_2{I}_{2,s_2}^{*}}{E_{2,s_2}^{*}{I}_2}}-{\displaystyle \frac{S{I}_2{E}_{2,s_2}^{*}}{S_{s_2}^{*}{I}_{2,s_2}^{*}{E}_2}}\right)\hfill \\ & +{\displaystyle \frac{be}{{\gamma }_1}}{I}_1\left(R_0^1-1\right)+\left({\alpha }_2{S}_{s_2}^{*}+{\eta }_2{V}_{s_2}^{*}-{\displaystyle \frac{fc}{{\gamma }_2}}\right)+{\displaystyle \frac{dg}{{\gamma }_3}}{I}_3\left(R_0^3-1\right).\hfill \end{array}$$

(54)

As the arithmetic mean is greater than or equal to the geometric mean, we will have

$$\begin{array}{c}\hfill 3-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_2}^{*}}{S_{s_2}^{*}V}}-{\displaystyle \frac{V}{V_{s_2}^{*}}}\le 0,\end{array}$$

(55)

$$3-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{E_2{I}_{2,s_2}^{*}}{E_{2,s_2}^{*}{I}_2}}-{\displaystyle \frac{S{I}_2{E}_{2,s_2}^{*}}{S_{s_2}^{*}{I}_{2,s_2}^{*}{E}_2}}\le 0$$

(56)

and

$$2-{\displaystyle \frac{S_{s_2}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_2}^{*}}}\le 0.$$

(57)

As ${\mathcal{E}}_{s_2}$ verifies condition $H_1$, when $R_0^1\le 1$ and $R_0^3\le 1$, we will have ${\dot{L}}_2\le 0$. Thus, the strain-2 endemic equilibrium point is globally asymptotically stable. □

For the global stability of the equilibrium point ${\mathcal{E}}_{s_3}$, we assume that this point verifies the following condition:

$${\alpha }_3{S}_{s_3}^{*}+{\eta }_3{V}_{s_3}^{*}-{\displaystyle \frac{dg}{{\gamma }_3}}\le 0.\left(H_2\right)$$

(58)

Theorem 4.

If $R_0^1\le 1,R_0^2\le 1$ and $R_0^3>1$, then the strain-3 endemic equilibrium point ${\mathcal{E}}_{s_3}$ is globally asymptotically stable.

Proof. 

We consider the following Lyapunov function $L_3$ in ${\mathbb{R}}_{+}^8$:

$$\begin{array}{cc}\hfill L_3& =S_{s_3}^{*}\left({\displaystyle \frac{S}{S_{s_3}^{*}}}-ln\left({\displaystyle \frac{S}{S_{s_3}^{*}}}\right)-1\right)+V_{s_3}^{*}\left({\displaystyle \frac{V}{V_{s_3}^{*}}}-ln\left({\displaystyle \frac{V}{V_{s_3}^{*}}}\right)-1\right)+E_1+E_2\hfill \\ & +E_{3,s_3}^{*}\left({\displaystyle \frac{E_3}{E_{3,s_3}^{*}}}-ln\left({\displaystyle \frac{E_3}{E_{3,s_3}^{*}}}\right)-1\right)+{\displaystyle \frac{b}{{\gamma }_1}}{I}_1+{\displaystyle \frac{c}{{\gamma }_2}}{I}_2+{\displaystyle \frac{d}{{\gamma }_3}}{I}_{3,s_3}^{*}\left({\displaystyle \frac{I_3}{I_{3,s_3}^{*}}}-ln\left({\displaystyle \frac{I_3}{I_{3,s_3}^{*}}}\right)-1\right).\hfill \end{array}$$

(59)

The time derivative of $L_3$ is given by

$$\begin{array}{cc}\hfill {\dot{L}}_3& =\left(1-{\displaystyle \frac{S_{s_3}^{*}}{S}}\right)\dot{S}+\left(1-{\displaystyle \frac{V_{s_3}^{*}}{V}}\right)\dot{V}+\dot{E_1}+\dot{E_2}+\left(1-{\displaystyle \frac{E_{3,s_3}^{*}}{E_3}}\right)\dot{E_3}+{\displaystyle \frac{b}{{\gamma }_1}}\dot{I_1}+{\displaystyle \frac{c}{{\gamma }_2}}{I}_2\hfill \\ & +{\displaystyle \frac{d}{{\gamma }_3}}\left(1-{\displaystyle \frac{I_{3,s_3}^{*}}{I_3}}\right)\dot{I_3}\hfill \\ & =\Lambda -\lambda S-\lambda V-{\displaystyle \frac{be}{{\gamma }_1}}{I}_1-{\displaystyle \frac{fc}{{\gamma }_2}}{I}_2-{\displaystyle \frac{dg}{{\gamma }_3}}{I}_3-\Lambda {\displaystyle \frac{S_{s_3}^{*}}{S}}+{\alpha }_1{S}_{s_3}^{*}{I}_1+{\alpha }_2{S}_{s_3}^{*}{I}_2+{\alpha }_3{S}_{s_3}^{*}{I}_3+(\sigma +\lambda )S_{s_3}^{*}\hfill \end{array}$$

(60)

$$\begin{array}{cc}& -\sigma S{\displaystyle \frac{V_{s_3}^{*}}{V}}+{\eta }_2{V}_{s_3}^{*}{I}_2+{\eta }_3{V}_{s_3}^{*}{I}_3+\lambda V_{s_3}^{*}-{\eta }_3{V}_{s_3}^{*}{I}_3{\displaystyle \frac{E_{3,s_3}^{*}}{E_3}}-{\alpha }_{3}S{I}_3{\displaystyle \frac{E_{3,s_3}^{*}}{E_3}}+d{E}_{3,s_3}^{*}-d{E}_3{\displaystyle \frac{I_{3,s_3}^{*}}{I_3}}\hfill \\ & +{\displaystyle \frac{dg}{{\gamma }_3}}{I}_{3,s_3}^{*}.\hfill \end{array}$$

(61)

As ${\mathcal{E}}_{s_3}$ is an equilibrium point of system (14), we will have

$$\left\{\begin{array}{ccc}\hfill \Lambda & ={\alpha }_3{S}_{s_3}^{*}{I}_{3,s_3}^{*}+(\sigma +\lambda )S_{s_3}^{*},\hfill & \\ \hfill \sigma S_{s_3}^{*}& =\lambda V_{s_3}^{*}+{\eta }_3{V}_{s_3}^{*}{I}_{3,s_3}^{*},\hfill & \\ \hfill d{E}_{3,s_3}^{*}& ={\alpha }_3{S}_{s_3}^{*}{I}_{3,s_3}^{*}+{\eta }_3{V}_{s_3}^{*}{I}_{3,s_3}^{*},\hfill & \\ \hfill E_{3,s_3}^{*}& ={\displaystyle \frac{g}{{\gamma }_3}}{I}_{3,s_3}^{*}.\hfill & \end{array}\right.$$

(62)

Thus, after some simplifications and factorizations, we will have

$$\begin{array}{cc}\hfill {\dot{L}}_3& =\sigma S_{s_3}^{*}\left(3-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_3}^{*}}{S_{s_3}^{*}V}}-{\displaystyle \frac{V}{V_{s_3}^{*}}}\right)+\lambda S_{s_3}^{*}\left(2-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_3}^{*}}}\right)\hfill \\ & +(d{E}_{3,s_3}^{*}-{\eta }_3{V}_{s_3}^{*}{I}_{3,s_3}^{*})\left(3-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{E_3{I}_{3,s_3}^{*}}{E_{3,s_3}^{*}{I}_3}}-{\displaystyle \frac{S{I}_3{E}_{3,s_3}^{*}}{S_{s_3}^{*}{I}_{3,s_3}^{*}{E}_3}}\right)\hfill \\ & +I_1\left({\alpha }_1{S}_{s_3}^{*}-{\displaystyle \frac{be}{{\gamma }_1}}\right)+I_2\left({\alpha }_2{S}_{s_3}^{*}-{\displaystyle \frac{fc}{{\gamma }_2}}\right)+I_3\left({\alpha }_3{S}_{s_3}^{*}+{\eta }_3{V}_{s_3}^{*}-{\displaystyle \frac{dg}{{\gamma }_3}}\right)-{\eta }_{3}V{I}_2{\displaystyle \frac{E_{3,s_3}^{*}}{E_3}}\hfill \end{array}$$

(63)

$$\begin{array}{cc}& \le \sigma S_{s_3}^{*}\left(3-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_3}^{*}}{S_{s_3}^{*}V}}-{\displaystyle \frac{V}{V_{s_3}^{*}}}\right)+\lambda S_{s_3}^{*}\left(2-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_3}^{*}}}\right)\hfill \\ & +(d{E}_{3,s_3}^{*}-{\eta }_3{V}_{s_3}^{*}{I}_{3,s_3}^{*})\left(3-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{E_3{I}_{3,s_3}^{*}}{E_{3,s_3}^{*}{I}_3}}-{\displaystyle \frac{S{I}_3{E}_{3,s_3}^{*}}{S_{s_3}^{*}{I}_{3,s_3}^{*}{E}_3}}\right)\hfill \\ & +{\displaystyle \frac{be}{{\gamma }_1}}{I}_1\left(R_0^1-1\right)+{\displaystyle \frac{fc}{{\gamma }_2}}{I}_2\left(R_0^2-1\right)+I_2\left({\alpha }_3{S}_{s_3}^{*}+{\eta }_3{V}_{s_3}^{*}-{\displaystyle \frac{dg}{{\gamma }_3}}\right).\hfill \end{array}$$

(64)

As the arithmetic mean is greater than or equal to the geometric mean, we will have

$$\begin{array}{c}\hfill 3-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{S{V}_{s_3}^{*}}{S_{s_3}^{*}V}}-{\displaystyle \frac{V}{V_{s_3}^{*}}}\le 0,\end{array}$$

(65)

$$3-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{E_3{I}_{2,s_3}^{*}}{E_{3,s_3}^{*}{I}_3}}-{\displaystyle \frac{S{I}_3{E}_{3,s_3}^{*}}{S_{s_3}^{*}{I}_{3,s_3}^{*}{E}_3}}\le 0$$

(66)

and

$$2-{\displaystyle \frac{S_{s_3}^{*}}{S}}-{\displaystyle \frac{S}{S_{s_3}^{*}}}\le 0.$$

(67)

As ${\mathcal{E}}_{s_3}$ verifies condition $H_2$, when $R_0^1\le 1$ and $R_0^2\le 1$, we will have ${\dot{L}}_3\le 0$. Thus, the strain-3 endemic equilibrium is globally asymptotically stable. □

Theorem 5.

If $1<R_0^1<max\{R_0^2,R_0^3\}$, then the total endemic equilibrium point ${\mathcal{E}}_t$ is globally asymptotically stable.

Proof. 

We consider the following Lyapunov function $L_4$ in ${\mathbb{R}}_{+}^8$:

$$\begin{array}{cc}\hfill L_4& =S_t^{*}\left({\displaystyle \frac{S}{S_t^{*}}}-ln\left({\displaystyle \frac{S}{S_t^{*}}}\right)-1\right)+V_t^{*}\left({\displaystyle \frac{V}{V_t^{*}}}-ln\left({\displaystyle \frac{V}{V_t^{*}}}\right)-1\right)+E_{1,t}^{*}\left({\displaystyle \frac{E_1}{E_{1,t}^{*}}}-ln\left({\displaystyle \frac{E_1}{E_{1,t}^{*}}}\right)-1\right)\hfill \\ & +E_{2,t}^{*}\left({\displaystyle \frac{E_2}{E_{2,t}^{*}}}-ln\left({\displaystyle \frac{E_2}{E_{2,t}^{*}}}\right)-1\right)+E_{3,t}^{*}\left({\displaystyle \frac{E_3}{E_{3,t}^{*}}}-ln\left({\displaystyle \frac{E_3}{E_{3,t}^{*}}}\right)-1\right)+{\displaystyle \frac{b}{{\gamma }_1}}{I}_{1,t}^{*}\left({\displaystyle \frac{I_1}{I_{1,t}^{*}}}-ln\left({\displaystyle \frac{I_1}{I_{1,t}^{*}}}\right)-1\right)\hfill \\ & +{\displaystyle \frac{c}{{\gamma }_2}}{I}_{2,t}^{*}\left({\displaystyle \frac{I_2}{I_{2,t}^{*}}}-ln\left({\displaystyle \frac{I_2}{I_{2,t}^{*}}}\right)-1\right)+{\displaystyle \frac{d}{{\gamma }_3}}{I}_{3,t}^{*}\left({\displaystyle \frac{I_3}{I_{3,t}^{*}}}-ln\left({\displaystyle \frac{I_3}{I_{3,t}^{*}}}\right)-1\right).\hfill \end{array}$$

(68)

The time derivative of $L_4$ is given by

$$\begin{array}{cc}\hfill {\dot{L}}_4& =\left(1-{\displaystyle \frac{S_t^{*}}{S}}\right)\dot{S}+\left(1-{\displaystyle \frac{V_t^{*}}{V}}\right)\dot{V}+\left(1-{\displaystyle \frac{E_{1,t}^{*}}{E_1}}\right)\dot{E_1}+\left(1-{\displaystyle \frac{E_{2,t}^{*}}{E_2}}\right)\dot{E_2}+\left(1-{\displaystyle \frac{E_{3,t}^{*}}{E_3}}\right)\dot{E_3}\hfill \\ & +{\displaystyle \frac{b}{{\gamma }_1}}\left(1-{\displaystyle \frac{I_{1,t}^{*}}{I_1}}\right)\dot{I_1}+{\displaystyle \frac{c}{{\gamma }_2}}\left(1-{\displaystyle \frac{I_{2,t}^{*}}{I_2}}\right)\dot{I_2}+{\displaystyle \frac{d}{{\gamma }_3}}\left(1-{\displaystyle \frac{I_{3,t}^{*}}{I_3}}\right)\dot{I_3}\hfill \\ & =\Lambda -\lambda S-\lambda V-{\displaystyle \frac{be}{{\gamma }_1}}{I}_1-{\displaystyle \frac{fc}{{\gamma }_2}}{I}_2-{\displaystyle \frac{dg}{{\gamma }_3}}{I}_3-\Lambda {\displaystyle \frac{S_t^{*}}{S}}+{\alpha }_1{S}_t^{*}{I}_1+{\alpha }_2{S}_t^{*}{I}_2+{\alpha }_3{S}_t^{*}{I}_3+(\sigma +\lambda )S_t^{*}\hfill \end{array}$$

(69)

$$\begin{array}{cc}& -\sigma S{\displaystyle \frac{V_t^{*}}{V}}+{\eta }_2{V}_t^{*}{I}_2+{\eta }_3{V}_t^{*}{I}_3+\eta V_t^{*}-{\alpha }_{1}S{I}_1{\displaystyle \frac{E_{1,t}^{*}}{E_1}}+b{E}_{1,t}^{*}-{\alpha }_{2}S{I}_2{\displaystyle \frac{E_{2,t}^{*}}{E_2}}-{\eta }_{2}V{I}_2{\displaystyle \frac{E_{2,t}^{*}}{E_2}}+c{E}_{2,t}^{*}\hfill \\ & -{\alpha }_{3}S{I}_3{\displaystyle \frac{E_{3,t}^{*}}{E_3}}-{\eta }_{3}V{I}_3{\displaystyle \frac{E_{3,t}^{*}}{E_3}}+d{E}_{3,t}^{*}-b{E}_1{\displaystyle \frac{I_{1,t}^{*}}{I_1}}+{\displaystyle \frac{be}{{\gamma }_1}}{I}_{1,t}^{*}-c{E}_2{\displaystyle \frac{I_{2,t}^{*}}{I_2}}+{\displaystyle \frac{fc}{{\gamma }_2}}{I}_{2,t}^{*}-d{E}_3{\displaystyle \frac{I_{3,t}^{*}}{I_3}}+{\displaystyle \frac{dg}{{\gamma }_3}}{I}_{3,t}^{*}.\hfill \end{array}$$

(70)

As ${\mathcal{E}}_{s_4}$ is an equilibrium point of system (14), we will have

$$\left\{\begin{array}{ccc}\hfill \Lambda & ={\alpha }_1{S}_t^{*}{I}_{1,t}^{*}+{\alpha }_2{S}_t^{*}{I}_{2,t}^{*}+{\alpha }_3{S}_t^{*}{I}_{3,t}^{*}+(\sigma +\lambda )S_t^{*},\hfill & \\ \hfill \sigma S_t^{*}& =\lambda V_t^{*}+{\eta }_2{V}_t^{*}{I}_{2,t}^{*}+{\eta }_3{V}_t^{*}{I}_{3,t}^{*},\hfill & \\ \hfill b{E}_{1,t}^{*}& ={\alpha }_1{S}_t^{*}{I}_{1,t}^{*},c{E}_{2,t}^{*}={\alpha }_2{S}_t^{*}{I}_{2,t}^{*}+{\eta }_2{V}_t^{*}{I}_{2,t}^{*},d{E}_{3,t}^{*}={\alpha }_3{S}_t^{*}{I}_{3,t}^{*}+{\eta }_3{V}_t^{*}{I}_{3,t}^{*},\hfill & \\ \hfill E_{1,t}^{*}& ={\displaystyle \frac{f}{{\gamma }_1}}{I}_{1,t}^{*},E_{2,t}^{*}={\displaystyle \frac{c}{{\gamma }_2}}{I}_{2,t}^{*},E_{3,t}^{*}={\displaystyle \frac{d}{{\gamma }_3}}{I}_{3,t}^{*}.\hfill & \end{array}\right.$$

(71)

Thus, after some simplifications and factorizations, we will have

$$\begin{array}{cc}\hfill {\dot{L}}_4& =\sigma S_t^{*}\left(3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{S{V}_{s_t}^{*}}{S_t^{*}V}}-{\displaystyle \frac{V}{V_t^{*}}}\right)+\lambda S_t^{*}\left(2-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{S}{S_t^{*}}}\right)+b{E}_{1,t}^{*}\left(3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{E_1{I}_{1,t}^{*}}{E_{1,s_t}^{*}{I}_1}}-{\displaystyle \frac{S{I}_1{E}_{1,s_t}^{*}}{S_t^{*}{I}_{1,t}^{*}{E}_1}}\right)\hfill \\ & +(c{E}_{2,t}^{*}-{\eta }_2{V}_t^{*}{I}_{2,t}^{*})\left(3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{E_2{I}_{2,t}^{*}}{E_{2,s_t}^{*}{I}_2}}-{\displaystyle \frac{S{I}_2{E}_{2,t}^{*}}{S_t^{*}{I}_{2,t}^{*}{E}_3}}\right)\hfill \\ & +(d{E}_{3,t}^{*}-{\eta }_3{V}_t^{*}{I}_{3,t}^{*})\left(3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{E_3{I}_{3,t}^{*}}{E_{3,t}^{*}{I}_3}}-{\displaystyle \frac{S{I}_3{E}_{3,t}^{*}}{S_t^{*}{I}_{3,t}^{*}{E}_3}}\right)\hfill \\ & +I_1\left({\alpha }_1{S}_t^{*}-{\displaystyle \frac{be}{{\gamma }_1}}\right)+I_2\left({\alpha }_2{S}_t^{*}+{\eta }_2{V}_t^{*}-{\displaystyle \frac{fc}{{\gamma }_2}}\right)+I_3\left({\alpha }_3{S}_t^{*}+{\eta }_3{V}_t^{*}-{\displaystyle \frac{dg}{{\gamma }_3}}\right)\hfill \\ & -{\eta }_{2}V{I}_2{\displaystyle \frac{E_{2,t}^{*}}{E_2}}-{\eta }_{3}V{I}_3{\displaystyle \frac{E_{3,t}^{*}}{E_3}}.\hfill \end{array}$$

(72)

As the arithmetic mean is greater than or equal to the geometric mean, we will have

$$3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{S{V}_{s_t}^{*}}{S_t^{*}V}}-{\displaystyle \frac{V}{V_t^{*}}}\le 0,$$

(73)

$$2-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{S}{S_t^{*}}}\le 0,$$

(74)

$$3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{E_1{I}_{1,t}^{*}}{E_{1,s_t}^{*}{I}_1}}-{\displaystyle \frac{S{I}_1{E}_{1,s_t}^{*}}{S_t^{*}{I}_{1,t}^{*}{E}_1}}\le 0,$$

(75)

$$3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{E_2{I}_{2,t}^{*}}{E_{2,s_t}^{*}{I}_2}}-{\displaystyle \frac{S{I}_2{E}_{2,t}^{*}}{S_t^{*}{I}_{2,t}^{*}{E}_3}}\le 0,$$

(76)

and

$$3-{\displaystyle \frac{S_t^{*}}{S}}-{\displaystyle \frac{E_3{I}_{3,t}^{*}}{E_{3,t}^{*}{I}_3}}-{\displaystyle \frac{S{I}_3{E}_{3,t}^{*}}{S_t^{*}{I}_{3,t}^{*}{E}_3}}\le 0.$$

(77)

Moreover, if $R_0^1<R_0^2$ and $R_0^1<R_0^3$, we will have

$${\alpha }_1{S}_t^{*}-{\displaystyle \frac{be}{{\gamma }_1}}\le 0,$$

(78)

$${\alpha }_2{S}_t^{*}+{\eta }_2{V}_t^{*}-{\displaystyle \frac{fc}{{\gamma }_2}}\le 0$$

(79)

and

$${\alpha }_3{S}_t^{*}+{\eta }_3{V}_t^{*}-{\displaystyle \frac{dg}{{\gamma }_3}}\le 0.$$

(80)

Then, when $1<R_0^1<max\{R_0^2,R_0^3\}$, we will have ${\dot{L}}_4\le 0$.Thus, the total endemic equilibrium is globally asymptotically stable. □

4. Numerical Simulations

In order to confirm our theoretical results, some numerical simulations and discussions are presented in this section by using the value of the parameter given in

Table 2. The parameter values of system (1).

Parameters	Figure 2	Figure 3	Figure 4	Figure 5	Figure 6
$\Lambda$	1	1	1	1	1
${\alpha }_1$	$0.2$	$0.5$	$0.2$	$0.2$	$0.2$
${\alpha }_2$	$0.2$	$0.2$	$0.6$	$0.2$	$0.6$
${\alpha }_3$	$0.2$	$0.2$	$0.2$	$0.7$	$0.7$
$\sigma$	$0.2$	$0.2$	$0.2$	$0.2$	$0.2$
$\lambda$	$0.2$	$0.2$	$0.2$	$0.2$	$0.2$
${\eta }_2$	$0.01$	$0.01$	$0.01$	$0.01$	$0.01$
${\eta }_3$	$0.005$	$0.005$	$0.005$	$0.005$	$0.005$
${\gamma }_1$	$0.1$	$0.1$	$0.1$	$0.1$	$0.1$
${\gamma }_2$	$0.1$	$0.1$	$0.1$	$0.1$	$0.1$
${\gamma }_3$	$0.1$	$0.1$	$0.1$	$0.1$	$0.1$
${\mu }_1$	$0.1$	$0.1$	$0.1$	$0.1$	$0.1$
${\mu }_2$	$0.1$	$0.1$	$0.1$	$0.1$	$0.1$
${\mu }_3$	$0.1$	$0.1$	$0.1$	$0.1$	$0.1$

. We performed our numerical simulations with the Runge–Kutta method [44] to test the effect of the vaccine in reducing the infection; in reality, the vaccination against a strain of an infectious disease such as COVID-19 leads to a reduction in the total number of infections.

Figure 2 shows the dynamics of the infection of the three strains SVEIR model. In this figure, we can see that the curves $E_1,E_2,E_3,I_1,I_2,I_3$ drop to zero, except the curves representing the susceptible and vaccinated individuals. This behavior is clearly observed in the zoomed part of the same figure. The obtained numerical result perfectly coincides with our theoretical result given in Theorem 1 concerning the global stability of the disease-free equilibrium ($R_0^1=0.55<1,R_0^2=0.82<1$ and $R_0^3=0.68<1$). Figure 3 describes the evolution of the different three strains of $SVEIR$ model components. In this figure, we notice that the first strain persists in contrast to the others strains that die out. The zoomed part of the same figure more clearly shows the persistence of the first strain. The strain-1 basic reproduction number is greater than 1 ($R_0^1=1.38>1$), and the basic reproduction numbers of the other strains are less than 1 ($R_0^2=0.82<1,R_0^3=0.68<1$). This numerical result verifies our theoretical result given in Theorem 2 concerning the global stability of the strain-1 endemic equilibrium. Figure 4 shows the global stability of the strain-2 endemic equilibrium, it is easy to observe that strain-2 persists while the other strains die out. The zoomed part of the same figure more clearly shows the persistence of the second strain. The basic reproduction number of this strain is greater than 1 ($R_0^2=1.93>1$), while the basic reproduction numbers of the other strains are less than 1 ($R_0^1=0.55<1,R_0^3=0.68<1$), which confirms our theoretical result given in Theorem 3 concerning the global stability of the strain-2 endemic equilibrium. Figure 5 describes the global stability of the strain-3 endemic equilibrium. This figure shows that strain-3 persists, while the other strains die out. The zoomed part of the same figure more clearly shows the persistence of the third strain. The basic reproduction number of this strain is greater than 1 ($R_0^3=2.21>1$), while the basic reproduction numbers of the other strains are less than 1 ($R_0^1=0.55<1$, $R_0^2=0.82<1$), which confirms our theoretical result given in Theorem 4 concerning the global stability of this equilibrium point. For the case of the global stability of the total endemic equilibrium, Figure 6 shows that all the strains persist. The zoomed part of the same figure more clearly shows the persistence of the all acting strains. Indeed, the basic reproduction number of every strain is greater than 1 ($R_0^1=1.38>1,R_0^2=1.93>1$ and $R_0^3=2.21>1$); this numerical result is in a good argument with our theoretical result given in Theorem 5 concerning the global stability of this equilibrium point.

5. Conclusions

In this paper, we analyzed a three-strain epidemic model with a vaccination strategy. The model contains nine compartments, namely susceptible individuals, vaccinated individuals, the three categories of exposed individuals, the three categories of infected individuals, and the recovered individuals. We started the analysis of this model by giving the different results of the existence, positivity, and boundedness of the model solutions. The suggested model has five steady states, namely the disease-free equilibrium, the strain-1 endemic equilibrium, the strain-2 endemic equilibrium, the strain-3 endemic equilibrium, and the total endemic equilibrium. By using the new generation method, we obtained the three basic reproduction numbers $R_0^1,R_0^2$, and $R_0^3$. Next, by using some suitable Lyapunov functions, we gave the global stability of the different steady states. This stability depends on different values of basic reproduction numbers. More precisely, if the three basic reproduction numbers are less than 1, the free equilibrium point is globally asymptotically stable. In addition, if $R_0^1>1,R_0^2\le 1$, and $R_0^3\le 1$, the strain-1 endemic equilibrium point is globally asymptotically stable, while the strain-2 endemic equilibrium point and strain-3 endemic equilibrium point are, respectively, globally asymptotically stable if $R_0^2>1,R_0^1\le 1$, and $R_0^3\le 1$; and $R_0^3>1,R_0^2\le 1$, and $R_0^2\le 1$. Finally, the total endemic equilibrium is globally asymptotically stable if all reproduction numbers are greater than 1. We showed that any strain with a higher reproduction number value outperforms the other strains. Numerical simulations were given in order to confirm and validate our theoretical results. It was observed that, in order to eradicate an infection, the basic reproduction numbers of all the strains must be less than unity.