1 Introduction

Nowadays, several infectious diseases are still targeting huge populations. They are considered amongst the principal causes of mortality, especially in many developing countries. Accordingly, mathematical modelling in epidemiology occupies more and more an increasingly preponderant place in public health research. This research discipline contributes indeed to well understand the studied epidemiological phenomenon and apprehend the different factors that can lead to a severe epidemic or even to a dangerous pandemic worldwide. The classical susceptible-infected-recovered (SIR) epidemic model was first introduced in [1]. Nevertheless, in many cases, the infection incubation period may take a long time interval. In this period, an incubated individual remains latent but not yet infectious. Therefore, another class of exposed individuals should be added to SIR and the new epidemic model will have SEIR abbreviation. Furthermore, epidemiological studies have revealed that the phenomenon of mutation causes more and more resistant viruses giving appearance of many new harmful epidemics or even new dangerous pandemics. Indeed, H1N1 flu virus is considered as mutation of the seasonal influenza [2, 3]. Also, the late coronavirus disease COVID-19 caused by the severe acute respiratory syndrome-related coronavirus SARS-Cov-2 is classified as a strain of SARS-CoV-1 [4]. Other processes of mutation were observed in many infections such as tuberculosis, human immunodeficiency virus and dengue fever [5,6,7]. For this reason, the multi-strain SEIR epidemic models present an important tool to study several infectious diseases that include a long incubation period and also various infection strains. The relevance of studying multi-strain models is to find out the different conditions permitting the coexistence of all acting strains. The global dynamics of one-strain SEIR model have been the subject of many investigations by considering either bilinear or nonlinear incidence rates [8,9,10]. The global stability of two-strain SEIR model have been tackled in [11], the authors include to the studied model two incidence functions, the first one is bilinear, while the second function is non-monotonic. Recently, the same problem was studied in [12] by assuming that the two incidence functions are non-monotonic. It is worthy to mention that the incidence rate gives more information of the disease transmission. Hence, the general incidence function has as goal to represent a large set of infection incidence rates. Accordingly, the purpose of this work is to generalize the previous models by taking into account a multi-strain SEIR model with two general incidence rates. Hence, our study will be carried out on the following two- strain generalized epidemic model:

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\mathrm{d}S}{\mathrm{d}t} = \varLambda -f(S,I_1) I_1-g(S,I_2) I_2- \delta S , \\&\frac{\mathrm{d}E_{1}}{\mathrm{d}t} = f(S,I_1) I_1-(\gamma _{1}+\delta )E_{1}, \\&\frac{\mathrm{d}E_{2}}{\mathrm{d}t} = \;g(S,I_2) I_2-(\gamma _{2}+\delta )E_{2}, \\&\frac{\mathrm{d}I_{1}}{\mathrm{d}t} = \gamma _{1} E_{1}-(\mu _{1}+\delta )I_{1}, \\&\frac{\mathrm{d}I_{2}}{\mathrm{d}t} = \gamma _{2} E_{2}-(\mu _{2}+\delta )I_{2},\\&\frac{\mathrm{d}R}{\mathrm{d}t} = \mu _{1}I_{1}+\mu _{2}I_{2}-\delta R, \end{aligned} \right. \end{aligned}$$
(1.1)

with

$$\begin{aligned}&S(0)\ge 0 ,\quad E_{1}(0)\ge 0 ,\quad E_{2}(0)\ge 0 , \\&I_{1}(0)\ge 0 ,\quad I_{2}(0)\ge 0 ,\quad R(0)\ge 0. \end{aligned}$$

where (S) is the number of susceptible individuals, \((E_{1})\) and \((E_{2})\) are, respectively, the numbers of each latent individuals class, \((I_{1})\) and \((I_{2})\) are, respectively, the numbers of each infectious individuals class and (R) is the number of removed individuals. The parameter \(\varLambda \) is the recruitment rate, \(\delta \) is the death rate of the population, \(\gamma _1\) and \( \gamma _2\) are, respectively, the latency rates of strain 1 and strain 2, \(\mu _1\) and \( \mu _2\) are, respectively, the two-strain transfer rates from infected to recovered. The general incidence functions \(f(S,I_1)\) and \(g(S,I_2)\) stand for the infection transmission rates for strain 1 and strain 2, respectively. The incidence functions \(f(S,I_1)\) and \(g(S,I_2)\) are assumed to be continuously differentiable in the interior of \({\mathbb {R}}^{2}_{+}\) and satisfy the same properties as in [13, 14]:

The properties \((H_1)\), \((H_2)\) and \((H_3)\), for the both functions f and g, are easily verified by several classical biological incidence rates such as the bilinear incidence function \(\beta S\)  [1, 15, 16], the saturated incidence function \(\displaystyle \frac{\beta S}{1+ \alpha _{1} S}\) or \(\displaystyle \frac{\beta S}{1+ \alpha _{2} I}\) [17, 18], Beddington–DeAngelis incidence function \(\displaystyle \frac{\beta S}{1+ \alpha _{1} S+\alpha _{2} I}\) [19,20,21], Crowley–Martin incidence function \(\displaystyle \frac{\beta S}{1+ \alpha _{1} S+\alpha _{2} I+ \alpha _{1}\alpha _{2} S I}\) [22,23,24], the specific nonlinear incidence function \(\displaystyle \frac{\beta S}{1+ \alpha _{1} S+\alpha _{2} I+ \alpha _{3} S I}\) [25,26,27,28,29] and non-monotone incidence function \( \displaystyle {\frac{\beta S}{1+ \alpha I^{2}}}\) [30,31,32,33,34]. The flowchart of the two-strain epidemiological SEIR model is illustrated in Fig. 1. Our main contribution centres around the global stability of multi-strain SEIR epidemic model with general incidence rates. In addition, a numerical comparison between our two-strain epidemic model results and COVID-19 clinical data will be conducted. It will be worthy to notice that the dynamics of SEIR COVID-19 epidemic model with two bilinear incidence functions was tackled in [35], and the authors introduce seasonality and stochasticity in order to describe the infection rate parameters. Taking into account non-monotone incidence function, the technique of sliding mode control was used to study an SEIR epidemic model describing COVID-19 disease [36]. The COVID-19 SEIR epidemiological model with Crowley–Martin incidence rate was studied [37], and the authors study the effect of different parameters on the disease spread. In our numerical comparison with COVID-19 clinical data, we will take into account the three latter incidence functions along with Beddington–DeAngelis incidence rate. A brief analysis to a more complex compartmental epidemic model will be established in Appendix of this paper.

The rest of this paper is outlined as follows. In the next section, we will establish the wellposedness of the suggested model by proving the existence, positivity and boundedness of solutions. In Sect. 3, we present an analysis of the model, we calculate the basic reproduction number of our epidemic model and we prove the global stability of the equilibria. Numerical simulations are given in Sect. 4 by using different specific incidence functions, and the comparison between the numerical results and COVID-19 clinical data is conducted in the same section. The generalization of the problem to a more complex compartmental model as well as concluding remarks is given at the end of this paper.

Fig. 1
figure 1

Flowchart of two-strain SEIR model

Full size image

2 The problem wellposedness and steady states

2.1 Positivity and boundedness of solutions

For the problems dealing with population dynamics, all the variables must be positive and bounded. We will assume first that all the model parameters are positive.

Proposition 1

For all non-negative initial data, the solutions of the problem (1.1) exist, remain bounded and non-negative.

Moreover, we have \(N(t)\le {\frac{\varLambda }{\delta }}+ N(0)\).

Proof

By the fundamental theory of differential equations functional framework (see for instance [38] and the references therein), we confirm that there exists a unique local solution to the problem (1.1).

In order to prove the positivity result, we will show that any solution starting from non-negative orthant \({\mathbb {R}}^{6}_{+}=\{(S,E_1,E_2,I_1,I_2,R)\in {\mathbb {R}}^{6} : S \ge 0, \; E_{1} \ge 0, \; E_{2} \ge 0, \; I_{1}\ge 0, \; I_{2}\ge 0\; , R \ge 0\ \}\) remains there forever.

First, let

$$\begin{aligned} T =&\,\hbox {sup}\{\displaystyle {\tau \ge 0 \;|\; \forall t \in [0, \tau ]} \;\hbox {such that}\; S(t) \ge 0,\nonumber \\&E_{1}(t) \ge 0, \; E_{2}(t) \ge 0, \; I_{1}(t)\ge 0, \nonumber \\&I_{2}(t)\ge 0\; , R(t) \ge 0\} \end{aligned}$$
(2.1)

Let us now prove that \(T = +\infty \). Suppose that T is finite; by continuity of solutions, we have

$$\begin{aligned}&S(T) = 0\quad \hbox {or}\quad E_{1}(T) = 0 \quad \hbox {or}\quad E_{2}(T) = 0 \quad \hbox {or}\\&I_{1}(T) = 0 \quad \hbox {or}\quad I_{2}(T) = 0\quad \hbox {or}\quad R(T) = 0. \end{aligned}$$

If \(S(T) = 0\) before the other variables \(E_{1}\), \(E_{2}\), \(I_{1}\), \(I_{2}\), R, become zero. Therefore,

$$\begin{aligned} \displaystyle { \frac{\mathrm{d}S(T)}{\mathrm{d}t}=\lim _{t\rightarrow T^{-}} \frac{S(T)-S(t)}{T-t}=\lim _{t\rightarrow T^{-}}} \frac{-S(t)}{T-t} \le 0.\nonumber \\ \end{aligned}$$
(2.2)

From the first equation of system (1.1), we have

$$\begin{aligned} {\dot{S}}(T)= & {} \varLambda -f(S(T),I_1(T)) \;I_1(T)\nonumber \\&-\,g(S(T),I_2(T)) \;I_2(T)- \delta S(T), \end{aligned}$$
(2.3)

then,

$$\begin{aligned} {\dot{S}}(T) {=} \varLambda -f(0,I_1(T)) \;I_1(T)-g(0,I_2(T))\; I_2(T),\nonumber \\ \end{aligned}$$
(2.4)

However, from \((H_1)\) we have

$$\begin{aligned} {\dot{S}}(T) = \varLambda > 0 \end{aligned}$$
(2.5)

which presents a contradiction.

If \(E_{1}(T)= 0\) before S, \(E_{2}\), \(I_{1}\), \(I_{2}\) and R. Then,

$$\begin{aligned} \displaystyle { \frac{\mathrm{d}E_{1}(T)}{\mathrm{d}t}{=}\lim _{t\rightarrow T^{-}} \frac{E_{1}(T)-E_{1}(t)}{T-t} {=}\lim _{t\rightarrow T^{-}}} \frac{-E_{1}(t)}{T-t} \le 0.\nonumber \\ \end{aligned}$$
(2.6)

Again, from the second equation of the system (1.1) with the fact \(E_{1}(T) = 0\), we will have

$$\begin{aligned} \displaystyle { \frac{\mathrm{d}E_{1}(T)}{\mathrm{d}t}= f(S,I_1) I_1}. \end{aligned}$$
(2.7)

However, from \((H_1)\) and \((H_2)\), \(f(S,I_1)I_1\) is positive, then we will have

$$\begin{aligned} \displaystyle { \frac{\mathrm{d}E_{1}(T)}{\mathrm{d}t}} > 0. \end{aligned}$$
(2.8)

Also, if \( I_{1} = 0 \) before S, \(E_{1}\), \(E_{2}\), \(I_{2}\), R become zero then

$$\begin{aligned} \displaystyle { \frac{\mathrm{d}I_{1}(T)}{\mathrm{d}t}=\lim _{t\rightarrow T^{-}} \frac{I_{1}(T)-I_{1}(t)}{T-t}=\lim _{t\rightarrow T^{-}}} \frac{-I_{1}(t)}{T-t} \le 0.\nonumber \\ \end{aligned}$$
(2.9)

But from the fourth equation of the system (1.1) with \(I_{1}(T)=0 \), we will have

$$\begin{aligned} \displaystyle { \frac{\mathrm{d}I_{1}(T)}{\mathrm{d}t}}= \gamma _{1} E_{1}. \end{aligned}$$
(2.10)

Since \(\gamma _{1}>0\), we have

$$\begin{aligned} \displaystyle { \frac{\mathrm{d}I_{1}(T)}{\mathrm{d}t}}= \gamma _{1} E_{1} > 0. \end{aligned}$$
(2.11)

This leads to contradiction.

Similar proofs for \( E_{2}(t)\), \(I_{2}(t)\) and R(t).

We conclude that T could not be finite, which implies that \(S(t) \ge 0, \; E_{1}(t) \ge 0, \; E_{2}(t) \ge 0, \; I_{1}(t)\ge 0, \; I_{2}(t)\ge 0\; , R(t) \ge 0\) for all positive times. This proves the positivity of solutions.

About boundedness, let the total population

$$\begin{aligned} N(t)= S(t) {+} E_{1}(t) {+} E_{2}(t) {+} I_{1}(t) {+} I_{2}(t) + R(t).\nonumber \\ \end{aligned}$$
(2.12)

From the system (1.1), we have

$$\begin{aligned} \displaystyle {\frac{\mathrm{d}N(t)}{\mathrm{d}t} = \varLambda - \delta N(t)}, \end{aligned}$$
(2.13)

and therefore,

$$\begin{aligned} N(t)= \displaystyle {\frac{\varLambda }{\delta }} + K \mathrm{e}^{-\delta t}, \end{aligned}$$
(2.14)

at \(t=0\), we have

$$\begin{aligned} N(0)= \displaystyle {\frac{\varLambda }{\delta }} + K, \end{aligned}$$
(2.15)

then

$$\begin{aligned} N(t)= \displaystyle {\frac{\varLambda }{\delta }} + (N(0) - \displaystyle {\frac{\varLambda }{\delta }})\mathrm{e}^{-\delta t}, \end{aligned}$$
(2.16)

consequently

$$\begin{aligned} N(t)\le \displaystyle {\frac{\varLambda }{\delta }} + N(0)\mathrm{e}^{-\delta t}, \end{aligned}$$
(2.17)

since \( 0<\mathrm{e}^{-\delta t}\le 1\), for all \( t\ge 0\), we conclude that

$$\begin{aligned} N(t)\le \displaystyle {\frac{\varLambda }{\delta }}+ N(0) \; . \end{aligned}$$
(2.18)

This implies that N(t) is bounded, and so are \(S(t), E_{1}(t), E_{2}(t), I_{1}(t), I_{2}(t) \;\text {and}\; R(t).\) Thus, the local solution can be prolonged to any positive time, which means that the unique solution exists globally. \(\square \)

2.2 The steady states

In this subsection, we show that there exist a disease-free equilibrium and three endemic equilibria. First, since the first five equations of the system (1.1) are independent of R and knowing that the number of the total population verifies the Eq. (2.14), so we can omit the sixth equation of the system (1.1). Therefore, the problem can be reduced to:

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\mathrm{d}S}{\mathrm{d}t} = \varLambda -f(S,I_1) I_1-g(S,I_2) I_2- \delta S , \\&\frac{\mathrm{d}E_{1}}{\mathrm{d}t} = f(S,I_1) I_1-(\gamma _{1}+\delta )E_{1}, \\&\frac{\mathrm{d}E_{2}}{\mathrm{d}t} = g(S,I_2) I_2-(\gamma _{2}+\delta )E_{2}, \\&\frac{\mathrm{d}I_{1}}{\mathrm{d}t} = \gamma _{1} E_{1}-(\mu _{1}+\delta )I_{1}, \\&\frac{\mathrm{d}I_{2}}{\mathrm{d}t} = \gamma _{2} E_{2}-(\mu _{2}+\delta )I_{2},\\ \end{aligned} \right. \end{aligned}$$
(2.19)

with

$$\begin{aligned} R=N-S-E_{1}-E_{2}-I_{1}-I_{2}. \end{aligned}$$
(2.20)

As usual, the basic reproduction number can be defined as the average number of new cases of an infection caused by one infected individual when all the population individuals are susceptibles. We will use the next generation matrix \(F V^{-1}\) to calculate the basic reproduction number \(R_{0}\). The formula that gives us the basic reproduction number is: \(R_{0} = \rho (F V^{-1}),\) where \(\rho \) stands for the spectral radius, F is the non-negative matrix of new infection cases and V is the matrix of the transition of infections associated with the model (2.19). We have

$$\begin{aligned} F&= \left( \begin{matrix} 0 &{} 0 &{}f(\frac{\varLambda }{\delta },0) &{} 0 \\ 0 &{} 0 &{} 0 &{}g(\frac{\varLambda }{\delta },0)\\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{matrix} \right) ,\\ V&= \left( \begin{matrix} \gamma _{1}+\delta &{} 0 &{} 0 &{} 0 \\ 0 &{} \gamma _{2}+\delta &{} 0 &{} 0 \\ -\gamma _{1} &{} 0 &{} \mu _{1}+\delta &{} 0\\ 0 &{} -\gamma _{2} &{} 0 &{} \mu _{2}+\delta \end{matrix} \right) . \end{aligned}$$

So,

$$\begin{aligned} F V^{-1}= \left( \begin{matrix} \displaystyle {\frac{ f(\frac{\varLambda }{\delta },0)\gamma _{1}}{(\gamma _{1}+\delta ) (\mu _{1}+\delta )}} &{} 0 &{} \displaystyle {\frac{f(\frac{\varLambda }{\delta },0)}{(\mu _{1}+\delta )}} &{} 0 \\ 0 &{} \displaystyle {\frac{g(\frac{\varLambda }{\delta },0) \gamma _{2}}{(\gamma _{2}+\delta )(\mu _{2}+\delta )}} &{} 0 &{} \displaystyle {\frac{g(\frac{\varLambda }{\delta },0)}{(\mu _{1}+\delta )}}\\ 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 \end{matrix} \right) . \end{aligned}$$

The basic reproduction number is the spectral radius of the matrix \(F V^{-1}\). This fact implies that

$$\begin{aligned} R_{0}=\max \{R^{1}_{0}, R^{2}_{0}\}, \end{aligned}$$
(2.21)

with

$$\begin{aligned} R^{1}_{0}=\displaystyle {\frac{ f(\frac{\varLambda }{\delta },0)\gamma _{1}}{(\gamma _{1}+\delta ) (\mu _{1}+\delta )}} \end{aligned}$$
(2.22)

and

$$\begin{aligned} R^{2}_{0}=\displaystyle {\frac{g(\frac{\varLambda }{\delta },0) \gamma _{2}}{(\gamma _{2}+\delta )(\mu _{2}+\delta )}}. \end{aligned}$$
(2.23)

We denote

$$\begin{aligned} a = \gamma _{1}+\delta , \; b = \gamma _{2}+\delta , \; c = \mu _{1}+\delta , \; e= \mu _{2}+\delta , \end{aligned}$$

then

$$\begin{aligned} R^{1}_{0}=\displaystyle {\frac{ f(\frac{\varLambda }{\delta },0)\gamma _{1}}{ac}} \end{aligned}$$
(2.24)

and

$$\begin{aligned} R^{2}_{0}=\displaystyle {\frac{g(\frac{\varLambda }{\delta },0) \gamma _{2}}{be}}. \end{aligned}$$
(2.25)

We call \(R^{1}_{0}\) the strain 1 reproduction number and \(R^{2}_{0}\) the strain 2 reproduction number.

Theorem 1

The problem (2.19) have the disease-free equilibrium \(E_f\) and three endemic equilibria \({\mathcal {E}}_{s_{1}}\), \({\mathcal {E}}_{s_{2}}\) and \({\mathcal {E}}_{s_{t}}\). Moreover, we have

  • The strain 1 endemic equilibrium \({\mathcal {E}}_{s_{1}}\) exists when \(R^{1}_{0}>1\).

  • The strain 2 endemic equilibrium \({\mathcal {E}}_{s_{2}}\) exists when \(R^{2}_{0}>1\).

  • The third endemic equilibrium \({\mathcal {E}}_{s_{t}}\) exists when \(R^{1}_{0}>1\) and \(R^{1}_{0}>1\).

Proof

In order to find the steady states of the system (2.19), we solve the following equations

$$\begin{aligned}&\varLambda -f(S,I_1) I_1-g(S,I_2) I_2- \delta S =0, \end{aligned}$$
(2.26)
$$\begin{aligned}&f(S,I_1) I_1-(\gamma _{1}+\delta )E_{1}=0, \end{aligned}$$
(2.27)
$$\begin{aligned}&g(S,I_2) I_2-(\gamma _{2}+\delta )E_{2}=0, \end{aligned}$$
(2.28)
$$\begin{aligned}&\gamma _{1} E_{1}-(\mu _{1}+\delta )I_{1}=0, \end{aligned}$$
(2.29)
$$\begin{aligned}&\gamma _{2} E_{2}-(\mu _{2}+\delta )I_{2}=0. \end{aligned}$$
(2.30)

From where, we obtain

  • When \(I_1=0 \;\text {and}\; I_2=0\), we find the disease-free equilibrium

    $$\begin{aligned} {\mathcal {E}}_{f}=\left( \displaystyle {\frac{\varLambda }{\delta }},0,0,0,0\right) . \end{aligned}$$

  • When \(I_1\ne 0 \;\text {and}\; I_2=0\), we find the strain 1 endemic equilibrium defined as follows

    $$\begin{aligned}&{\mathcal {E}}_{s_{1}}=\left( S^{*}_{1},\frac{1}{a}(\varLambda -\delta S^{*}_{1}),0, \frac{\gamma _{1}}{ac}(\varLambda -\delta S^{*}_{1}),0\right) ,\quad \\&\hbox {where } \; S^{*}_{1} \in \left[ 0,\displaystyle \frac{\varLambda }{\delta }\right] . \end{aligned}$$

    Define now a function \(\varPsi \; \text{ on } \; [0, +\infty [ \) as follows

    $$\begin{aligned} \varPsi (S)= f(S,\displaystyle {\frac{\gamma _{1}}{ac}} (\varLambda -\delta S)) - \frac{ac}{\gamma _{1}}. \end{aligned}$$
    (2.31)

    We have

    $$\begin{aligned} \displaystyle {\frac{\partial \varPsi (S)}{\partial S}}= \frac{\partial f(S,I_1)}{\partial S}+\frac{\partial f(S,I_1)}{\partial I_1} \left( \displaystyle \frac{-\delta \gamma _{1}}{ac}\right) ,\nonumber \\ \end{aligned}$$
    (2.32)

    using the conditions \((H_2)\) and \((H_3)\), we deduce that

    $$\begin{aligned} \displaystyle {\frac{\partial \varPsi (S)}{\partial S}} \ge 0. \end{aligned}$$
    (2.33)

    However, \(\varPsi (0)= f(0,I^{*}_{1,s_{1}})- \displaystyle \frac{ac}{\gamma _1}= - \displaystyle \frac{ac}{\gamma _1} < 0.\) Therefore, for \(R^{1}_{0} > 1\), we he have

    $$\begin{aligned} \displaystyle {\varPsi \left( \frac{\varLambda }{\delta }\right) }= f\left( \frac{\varLambda }{\delta },0\right) - \displaystyle \frac{ac}{\gamma _1}= \displaystyle \frac{ac}{\gamma _1}(R^{1}_{0}-1) > 0.\nonumber \\ \end{aligned}$$
    (2.34)

    Hence, there exists a unique strain 1 endemic equilibrium

    $$\begin{aligned} \displaystyle {{\mathcal {E}}_{s_{1}}=\left( S^{*}_{1},E^{*}_{1,s_{1}},E^{*}_{2,s_{1}},I^{*}_{1,s_{1}},I^{*}_{2,s_{1}}\right) }, \end{aligned}$$
    (2.35)

    with \(S^{*}_{1} \in \big ]0,\displaystyle \frac{\varLambda }{\delta }\big [\), \(E^{*}_{1,s_{1}}>0\), \(I^{*}_{1,s_{1}}>0\) and \(E^{*}_{2,s_{1}}=I^{*}_{2,s_{1}}=0\).

  • When \(I_2\ne 0 \;\text {and}\; I_1=0\), we find the strain 2 endemic equilibrium defined as follows

    $$\begin{aligned}&{\mathcal {E}}_{s_{2}}=\left( S^{*}_{2},0, \frac{1}{b}(\varLambda -\delta S^{*}_{2}),0,\frac{\gamma _{2}}{be}(\varLambda -\delta S^{*}_{2})\right) , \quad \\&\hbox {where } S^{*}_{2} \in \left[ 0,\displaystyle \frac{\varLambda }{\delta }\right] . \end{aligned}$$

    Define also a function \(\varPhi \; \text{ on } \; [0, +\infty [ \) as follows

    $$\begin{aligned} \varPhi (S)= g(S,\displaystyle {\frac{\gamma _{2}}{be}(\varLambda -\delta S)}) - \frac{be}{\gamma _{2}}. \end{aligned}$$
    (2.36)

    We have

    $$\begin{aligned} \displaystyle {\frac{\partial \varPhi (S)}{\partial S}}= \frac{\partial g(S,I_2)}{\partial S}+\frac{\partial g(S,I_2)}{\partial I_2} \left( \displaystyle \frac{-\delta \gamma _{2}}{be}\right) ,\nonumber \\ \end{aligned}$$
    (2.37)

    using the conditions \((H_2)\) and \((H_3)\), we conclude that

    $$\begin{aligned} \displaystyle {\frac{\partial \varPhi (S)}{\partial S}} \ge 0. \end{aligned}$$
    (2.38)

    However, \(\varPhi (0)= g(0,I^{*}_{2,s_{2}})- \displaystyle \frac{be}{\gamma _2}= - \displaystyle \frac{be}{\gamma _2} < 0\). So, for \(R^{2}_{0} > 1\), we have

    $$\begin{aligned} \displaystyle {\varPhi \left( \frac{\varLambda }{\delta }\right) }= g\left( \frac{\varLambda }{\delta },0\right) - \displaystyle \frac{be}{\gamma _2}= \displaystyle \frac{be}{\gamma _2}(R^{2}_{0}-1) > 0.\nonumber \\ \end{aligned}$$
    (2.39)

    Hence, there exists a unique strain 2 endemic equilibrium

    $$\begin{aligned} \displaystyle {{\mathcal {E}}_{s_{2}}=\left( S^{*}_{2},E^{*}_{1,s_{2}},E^{*}_{2,s_{2}},I^{*}_{1,s_{2}},I^{*}_{2,s_{2}}\right) }, \end{aligned}$$
    (2.40)

    with \(S^{*}_{2} \in \left]0,\displaystyle \frac{\varLambda }{\delta }\right[\), \(E^{*}_{2,s_{2}}>0\), \(I^{*}_{2,s_{2}}>0\) and \(E^{*}_{1,s_{2}}=I^{*}_{1,s_{2}}=0\).

  • When \(I_1\ne 0 \;\text {and}\; I_2 \ne 0\), we find the third endemic equilibrium defined as follows

    $$\begin{aligned} \displaystyle {{\mathcal {E}}_{t}=\left( S^{*}_{t},E^{*}_{1,t},E^{*}_{2,t},I^{*}_{1,t},I^{*}_{2,t}\right) }, \end{aligned}$$
    (2.41)

    where

    $$\begin{aligned} E^{*}_{1,t}= & {} \displaystyle {\frac{c}{\gamma _{1}}I^{*}_{1,t}} \; , \; E^{*}_{2,t} = \displaystyle {\frac{e}{\gamma _{2}}I^{*}_{2,t}}, \end{aligned}$$
    (2.42)
    $$\begin{aligned} S^{*}_{t}= & {} \displaystyle \frac{1}{\delta }\left[ \varLambda -\displaystyle \frac{f(\frac{\varLambda }{\delta },0)}{R^{1}_{0}} I^{*}_{1,t}- \displaystyle \frac{g(\frac{\varLambda }{\delta },0)}{R^{2}_{0}} I^{*}_{2,t}, \right] \nonumber \\ \end{aligned}$$
    (2.43)

    with \( \varLambda \ge \displaystyle \frac{f(\frac{\varLambda }{\delta },0)}{R^{1}_{0}} I^{*}_{1,t}+\displaystyle \frac{g(\frac{\varLambda }{\delta },0)}{R^{2}_{0}} I^{*}_{2,t}\), \(R^{1}_{0}>1\) and \(R^{2}_{0}>1\). \(\square \)

3 Global stability of equilibria

3.1 Global stability of disease-free equilibrium

Theorem 2

If \(R_{0}\le 1\), then the disease-free equilibrium \({\mathcal {E}}_{f}\) is globally asymptotically stable.

Proof

First, we consider the following Lyapunov function in \({\mathbb {R}}^{5}_{+}\):

$$\begin{aligned} {\mathcal {L}}_f(S,E_{1},E_{2},I_{1},I_{2})= & {} S-S_{0}^{*}-\int ^{S}_{S_{0}^{*}}\frac{f(S_{0}^{*},0)}{f(X,0)}\,\mathrm{d}X\nonumber \\&\,+E_{1}+E_{2}+\frac{a}{\gamma _{1}} I_{1}+\frac{b}{\gamma _{2}} I_{2}.\nonumber \\ \end{aligned}$$
(3.1)

The time derivative is given by

$$\begin{aligned} \dot{{\mathcal {L}}}_f(S,E_{1},E_{2},I_{1},I_{2})&= {\dot{S}}-\frac{f(S_{0}^{*},0)}{f(S,0)} {\dot{S}}+\dot{E_{1}}+\dot{E_{2}}\nonumber \\&\quad +\,\frac{a}{\gamma _{1}}\dot{I_{1}}+ \frac{b}{\gamma _{2}}\dot{I_{2}}, \nonumber \\&=\delta S_{0}^{*}\left( 1-\frac{S}{S_{0}^{*}} \right) \left( 1- \frac{f(S_{0}^{*},0)}{f(S,0)}\right) \nonumber \\&\quad +\,\frac{ac}{\gamma _{1}}I_1 \left( \frac{f(S,I_1)}{f(S,0)} R_{0}^{1}-1\right) \end{aligned}$$
(3.2)
$$\begin{aligned}&\quad +\,\frac{be}{\gamma _{2}}I_2 \left( \frac{f(S_{0}^{*},0)}{f(S,0)}\frac{g(S,I_2)}{g(S_{0}^{*},0)} R_{0}^{2}-1 \right) \end{aligned}$$
(3.3)
$$\begin{aligned}&\nonumber \le \delta S_{0}^{*}\left( 1-\frac{S}{S_{0}^{*}} \right) \left( 1- \frac{f(S_{0}^{*},0)}{f(S,0)}\right) \\&\quad +\,\frac{ac}{\gamma _{1}}I_1 \left( R_{0}^{1}-1\right) \nonumber \\&\quad +\,\frac{be}{\gamma _{2}}I_2 \left( \displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)}\displaystyle \frac{g(S,I_2)}{g(S_{0}^{*},0)} R_{0}^{2}-1 \right) \end{aligned}$$
(3.4)

We will discuss two cases:

  • If \( S \le S_{0}^{*}\), using \((H_2),\) we will have \(\displaystyle \frac{g(S,I_2)}{g(S_{0}^{*},0)} \le 1\), then,

    $$\begin{aligned} \dot{{\mathcal {L}}}_f(S,E_{1},E_{2},I_{1},I_{2})&\le \delta S_{0}^{*}\left( 1-\frac{S}{S_{0}^{*}} \right) \left( 1- \frac{f(S_{0}^{*},0)}{f(S,0)}\right) \nonumber \\&\quad +\,\frac{ac}{\gamma _{1}}I_1\left( R_{0}^{1}-1\right) \nonumber \\&\quad +\,\frac{be}{\gamma _{2}}I_2 \left( \displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)} R_{0}^{2}-1 \right) . \end{aligned}$$
    (3.5)

    Since \(R_{0}^{2} \le \displaystyle \frac{f(S,0)}{f(S_{0}^{*},0)} \le 1\), we obtain

    $$\begin{aligned} \displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)} R_{0}^{2}-1 \le 0. \end{aligned}$$
    (3.6)

    Otherwise, \(1- \displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)} \le 0,\) therefore

    $$\begin{aligned} \delta S_{0}^{*}\left( 1-\displaystyle \frac{S}{S_{0}^{*}} \right) \left( 1- \displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)}\right) \le 0. \end{aligned}$$
    (3.7)

  • If \( S_{0}^{*} < S \), using \((H_2),\) we will have \(\displaystyle \frac{g(S,I_2)}{g(S_{0}^{*},0)} > 1\) and \(\displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)} < 1\) then,

    $$\begin{aligned} \dot{{\mathcal {L}}}_f(S,E_{1},E_{2},I_{1},I_{2})&\le \delta S_{0}^{*}\left( 1-\frac{S}{S_{0}^{*}} \right) \left( 1- \frac{f(S_{0}^{*},0)}{f(S,0)}\right) \nonumber \\&\quad +\,\frac{ac}{\gamma _{1}}I_1\left( R_{0}^{1}-1\right) \nonumber \\&\quad +\, \frac{be}{\gamma _{2}}I_2 \left( \displaystyle \frac{g(S,I_2)}{g(S_{0}^{*},0)} R_{0}^{2}-1 \right) . \end{aligned}$$
    (3.8)

    Since \(R_{0}^{2}< \displaystyle \frac{g(S_{0}^{*},0)}{g(S,I_2)} < 1\), we obtain

    $$\begin{aligned} \displaystyle \frac{g(S,I_2)}{g(S_{0}^{*},0)} R_{0}^{2}-1 < 0. \end{aligned}$$
    (3.9)

    From \(\displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)} < 1\), we have

    $$\begin{aligned} \delta S_{0}^{*}\left( 1-\displaystyle \frac{S}{S_{0}^{*}} \right) \left( 1- \displaystyle \frac{f(S_{0}^{*},0)}{f(S,0)}\right) \le 0. \end{aligned}$$
    (3.10)

By the above discussion, we deduce that, if \(R_{0}^{2}\le 1\) and \(R_{0}^{1}\le 1\) (which means that \(R_{0} \le 1\)), then

$$\begin{aligned} \dot{{\mathcal {L}}}_f(S,E_{1},E_{2},I_{1},I_{2}) \le 0. \end{aligned}$$
(3.11)

Thus, the disease-free equilibrium point \({\mathcal {E}}_{f}\) is globally asymptotically stable when \(R_{0} \le 1\). \(\square \)

3.2 Global stability of strain 1 endemic equilibrium

For the global stability of \({\mathcal {E}}_{s_{1}}\), we assume that the function f satisfies the following condition:

$$\begin{aligned}&\left( 1- \frac{f(S,I_1)}{f(S,I^{*}_{1,s_{1}})}\right) \left( \frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}-\frac{I_1}{I^{*}_{1,s_{1}}} \right) \le 0,\nonumber \\&\qquad \qquad \forall \; S, I_1 > 0.\nonumber \\ \end{aligned}$$
(3.12)

Theorem 3

The strain 1 endemic equilibrium \({\mathcal {E}}_{s_{1}}\) is globally asymptotically stable when \(R^{2}_0 \le 1< R^{1}_0\).

Proof

First, we consider the Lyapunov function \({\mathcal {L}}_1\) defined by:

$$\begin{aligned}&{\mathcal {L}}_1(S,E_{1},E_{2},I_{1},I_{2}) = S-S_{1}^{*}-\int ^{S}_{S_{1}^{*}}\frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(X,I^{*}_{1,s_{1}})}\,\mathrm{d}X +\,E^{*}_{1}\left( \frac{E_{1}}{E^{*}_{1,s_{1}}}-\ln \left( \frac{E_{1}}{E^{*}_{1,s_{1}}}\right) -1\right) \nonumber \\&+\,E_{2}+\frac{a}{\gamma _{1}}I^{*}_{1,s_{1}}\left( \frac{I_{1}}{I^{*}_{1,s_{1}}}-\ln \left( \frac{I_{1}}{I^{*}_{1,s_{1}}}\right) -1\right) +\frac{b}{\gamma _{2}}I_{2}. \end{aligned}$$
(3.13)

The time derivative is given by

$$\begin{aligned} \dot{{\mathcal {L}}}_1(S,E_{1},E_{2},I_{1},I_{2}) =&\left( \varLambda -f(S,I_1)I_1-g(S,I_2)I_2- \delta S\right) \nonumber \\&\quad \displaystyle {\left( 1-\frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})}\right) }\nonumber \\&\nonumber +\,\left( f(S,I_1)I_1-aE_{1}\right) \displaystyle {\left( 1-\frac{E^{*}_{1,s_{1}}}{E_{1}}\right) }\\&+\,\displaystyle {\left( g(S,I_2)I_2-bE_{2}\right) }\nonumber \\&+\,\displaystyle {\frac{a}{\gamma _{1}} \left( \gamma _{1}E_{1}-cI_{1}\right) \left( 1-\frac{I^{*}_{1,s_{1}}}{I_{1}}\right) }\nonumber \\&+\,\displaystyle {\frac{b}{\gamma _{2}}}\left( \gamma _{2}E_{2}-eI_{2}\right) . \end{aligned}$$
(3.14)

We have

$$\begin{aligned} \left\{ \begin{aligned}&\varLambda = \delta S^{*}_{1}+f(S_{1}^{*},I^{*}_{1,s_{1}})I^{*}_{1,s_{1}},\\&f(S_{1}^{*},I^{*}_{1,s_{1}})I^{*}_{1,s_{1}}= \displaystyle \frac{ac}{\gamma _{1}}I^{*}_{1,s_{1}} =aE^{*}_{1,s_{1}},\\&\displaystyle \frac{E^{*}_{1,s_{1}}}{I^{*}_{1,s_{1}}}=\displaystyle \frac{c}{\gamma _{1}}. \end{aligned}\right. \end{aligned}$$
(3.15)

Therefore,

$$\begin{aligned} \dot{{\mathcal {L}}}_1(S,E_{1},E_{2},I_{1},I_{2})\; =&\nonumber \;\displaystyle {\delta S^{*}_{1}}\left( 1-\frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})} \right) \\&-\,\displaystyle {\delta S \left( 1-\frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})}\right) }\nonumber \\&\nonumber +\,f(S_{1}^{*},I^{*}_{1,s_{1}})I^{*}_{1,s_{1}}\nonumber \\&-\,f(S_{1}^{*},I^{*}_{1,s_{1}})I^{*}_{1,s_{1}} \frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})}\nonumber \\&\nonumber +\,\frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})}f(S,I_1)I_1\\&-\,f(S,I_1)\frac{I_1E^{*}_{1,s_{1}}}{E_1} +aE^{*}_{1,s_{1}}\nonumber \\&-\displaystyle \frac{ac}{\gamma _{1}}I_1-\,\frac{aI^{*}_{1,s_{1}}}{I_1}E_1 +\displaystyle \frac{ac}{\gamma _{1}}I^{*}_{1,s_{1}}\nonumber \\&+\,\frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})} \; g(S,I_2)I_2-\frac{be}{\gamma _{2}}I_2. \end{aligned}$$
(3.16)

Then,

$$\begin{aligned}&\dot{{\mathcal {L}}}_1(S,E_{1},E_{2},I_{1},I_{2})= \;aE^{*}_{1,s_{1}} \left( 4-\frac{aE^{*}_{1,s_{1}}}{f(S,I^{*}_{1,s_{1}})I^{*}_{1,s_{1}}}-\frac{f(S,I_1)I_1}{aE_1}-\frac{I^{*}_{1,s_{1}}E_1}{I_1E^{*}_{1,s_{1}}} -\frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}\right) \nonumber \\&\qquad +\,aE^{*}_{1,s_{1}}\left( \frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}+\frac{f(S,I_1)}{f(S,I^{*}_{1,s_{1}})}\frac{I_1}{I^{*}_{1,s_{1}}}-\frac{I_1}{I^{*}_{1,s_{1}}}-1 \right) \nonumber \\&\qquad +\,\frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})} \; g(S,I_2)I_2-\frac{be}{\gamma _{2}}I_2 \end{aligned}$$
(3.17)
$$\begin{aligned}&\dot{{\mathcal {L}}}_1(S,E_{1},E_{2},I_{1},I_{2})= \;aE^{*}_{1,s_{1}} \left( 4-\frac{aE^{*}_{1,s_{1}}}{f(S,I^{*}_{1,s_{1}})I^{*}_{1,s_{1}}}-\frac{f(S,I_1)I_1}{aE_1}-\frac{I^{*}_{1,s_{1}}E_1}{I_1E^{*}_{1,s_{1}}} -\frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}\right) \nonumber \\&\qquad +\,aE^{*}_{1,s_{1}}\left( \frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}+\frac{f(S,I_1)}{f(S,I^{*}_{1,s_{1}})}\frac{I_1}{I^{*}_{1,s_{1}}}-\frac{I_1}{I^{*}_{1,s_{1}}}-1 \right) \nonumber \\&\qquad +\, \frac{be}{\gamma _2}I_2 \left( \displaystyle \frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})} \; \displaystyle \frac{g(S,I_2)}{g(S_{0}^{*},0)}\; R_{0}^{2}-1\right) . \end{aligned}$$
(3.18)

From \((H_2)\) and \((H_3)\), we have \( \displaystyle \frac{g(S,I_2)}{g(S_{0}^{*},0)} \le 1,\) then,

$$\begin{aligned}&\dot{{\mathcal {L}}}_1(S,E_{1},E_{2},I_{1},I_{2})\le \;aE^{*}_{1,s_{1}} \left( 4-\frac{aE^{*}_{1,s_{1}}}{f(S,I^{*}_{1,s_{1}})I^{*}_{1,s_{1}}}-\frac{f(S,I_1)I_1}{aE_1}-\frac{I^{*}_{1,s_{1}}E_1}{I_1E^{*}_{1,s_{1}}} -\frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}\right) \nonumber \\&\qquad +\,aE^{*}_{1,s_{1}}\left( \frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}+\frac{f(S,I_1)}{f(S,I^{*}_{1,s_{1}})}\frac{I_1}{I^{*}_{1,s_{1}}}-\frac{I_1}{I^{*}_{1,s_{1}}}-1 \right) +\, \frac{be}{\gamma _2}I_2 \left( \displaystyle \frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})} {R_{0}^{2}}-1\right) , \end{aligned}$$
(3.19)

From (3.12), we have

$$\begin{aligned}&\left( \frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}+\frac{f(S,I_1)}{f(S,I^{*}_{1,s_{1}})}\frac{I_1}{I^{*}_{1,s_{1}}}-\frac{I_1}{I^{*}_{1,s_{1}}}-1 \right) \nonumber \\&\quad = \left( 1- \frac{f(S,I_1)}{f(S,I^{*}_{1,s_{1}})}\right) \left( \frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}-\frac{I_1}{I^{*}_{1,s_{1}}} \right) \le \;\;0, \end{aligned}$$
(3.20)

by the relation between arithmetic and geometric means, we have

$$\begin{aligned} \left( 4-\frac{aE^{*}_{1,s_{1}}}{f(S,I^{*}_{1,s_{1}})I^{*}_{1,s_{1}}}-\frac{f(S,I_1)I_1}{aE_1}-\frac{I^{*}_{1,s_{1}}E_1}{I_1E^{*}_{1,s_{1}}} -\frac{f(S,I^{*}_{1,s_{1}})}{f(S,I_1)}\right) \;\le \; 0,\nonumber \\ \end{aligned}$$
(3.21)

We discuss two cases:

  • If \(S_{1}^{*}\le S\), from \((H_2)\), we will have \(\displaystyle \frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})}\le 1,\) since \(\left( \displaystyle \frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})} R^{2}_0 -1\right) \le 0 \), we obtain, for \(R^{2}_0 \le 1,\) the following \(\dot{{\mathcal {L}}}_1\le 0.\)

  • If \(S \le S_{1}^{*},\) from \((H_2)\), we will have \(\displaystyle \frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})}\ge 1,\) since \(R^{2}_0 \le \displaystyle \frac{f(S,I^{*}_{1,s_{1}})}{f(S_{1}^{*},I^{*}_{1,s_{1}})} \le 1,\) we obtain \(\displaystyle \frac{f(S_{1}^{*},I^{*}_{1,s_{1}})}{f(S,I^{*}_{1,s_{1}})} R^{2}_0 -1 \le 0,\) which implies, \(\dot{{\mathcal {L}}}_1\le 0.\) By the above discussion, we deduce that if \(R^{2}_0 \le 1,\) we will have \( \dot{{\mathcal {L}}}_1\le 0.\)

We conclude that the steady state \({\mathcal {E}}_{s_{1}}\) is globally asymptotically stable when \( R^{2}_0 \le 1\) and \(1 < R^{1}_0\). \(\square \)

3.3 Global stability of strain 2 endemic equilibrium

For the global stability of \({\mathcal {E}}_{s_{2}}\), we assume that the function g satisfies the following condition:

$$\begin{aligned} \left( 1- \frac{g(S,I_2)}{g(S,I^{*}_{2,s_{2}})}\right) \left( \frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}-\frac{I_2}{I^{*}_{2,s_{2}}} \right) \le 0, \quad \forall \; S, I_2 > 0.\nonumber \\ \end{aligned}$$
(3.22)

Theorem 4

The strain 2 endemic equilibrium point \({\mathcal {E}}_{s_{2}}\) is globally asymptotically stable when \( R^{1}_0 \le 1 < R^{2}_0\).

Proof

First, we consider the Lyapunov function \({\mathcal {L}}_2\) defined by:

$$\begin{aligned} {\mathcal {L}}_2(S,E_{1},E_{2},I_{1},I_{2}) =&\nonumber S-S_{2}^{*}-\int ^{S}_{S_{2}^{*}}\frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(X,I^{*}_{2,s_{2}})}\,\mathrm{d}X\\&+\,E^{*}_{2,s_{2}}\left( \frac{E_{2}}{E^{*}_{2,s_{2}}}-\ln \left( \frac{E_{2}}{E^{*}_{2,s_{2}}}\right) -1\right) \nonumber \\&+\,E_{1}+\frac{a}{\gamma _{1}} I_{1}\nonumber \\&+\,\displaystyle {\frac{b}{\gamma _{2}}I^{*}_{2,s_{2}}}\left( \frac{I_{2}}{I^{*}_{2,s_{2}}} -\ln \left( \frac{I_{2}}{I^{*}_{2,s_{2}}}\right) -1\right) . \end{aligned}$$
(3.23)

The time derivative is given by

$$\begin{aligned} \dot{{\mathcal {L}}}_2(S,E_{1},E_{2},I_{1},I_{2}) =&\nonumber \left( \varLambda -f(S,I_1)I_1-g(S,I_2)I_2- \delta S\right) \\&\quad \displaystyle {\left( 1-\frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})}\right) }\nonumber \\&\nonumber +\,\left( g(S,I_2)I_2-bE_{2}\right) \displaystyle {\left( 1-\frac{E^{*}_{2,s_{2}}}{E_{2}}\right) }\\&+\,\displaystyle {\left( f(S,I_1)I_1 -aE_{1}\right) }\nonumber \\&+\,\displaystyle {\frac{b}{\gamma _{2}} \left( \gamma _{2}E_{2}-eI_{2}\right) \left( 1-\frac{I^{*}_{2,s_{2}}}{I_{2}}\right) }\nonumber \\&+\,\displaystyle {\frac{a}{\gamma _{1}}}\left( \gamma _{1}E_{1}-cI_{1}\right) . \end{aligned}$$
(3.24)

We have

$$\begin{aligned} \left\{ \begin{aligned}&\varLambda = \delta S^{*}_{2}+g(S_{2}^{*},I^{*}_{2,s_{2}})I^{*}_{2,s_{2}},\\&g(S_{2}^{*},I^{*}_{2,s_{2}})I^{*}_{2,s_{2}}= \displaystyle \frac{be}{\gamma _{2}}I^{*}_{2,s_{2}} =bE^{*}_{2,s_{2}},\\&\displaystyle \frac{E^{*}_{2,s_{2}}}{I^{*}_{2,s_{2}}}=\displaystyle \frac{e}{\gamma _{2}}. \end{aligned} \right. \end{aligned}$$
(3.25)

Therefore,

$$\begin{aligned} \dot{{\mathcal {L}}}_2(S,E_{1},E_{2},I_{1},I_{2}) =&\nonumber \;\displaystyle {\delta S^{*}_{2}}\left( 1-\frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})} \right) \\&-\,\displaystyle {\delta S \left( 1-\frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})}\right) }\nonumber \\&\nonumber +\,g(S_{2}^{*},I^{*}_{2,s_{2}})I^{*}_{2,s_{2}}\\&-\,g(S_{2}^{*},I^{*}_{2,s_{2}})I^{*}_{2,s_{2}} \frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})}\nonumber \\&\nonumber +\,\frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})}g(S,I_2)I_2\\&-g(S,I_2)\frac{I_2 E^{*}_{2,s_{2}}}{E_2}+\,bE^{*}_{2,s_{2}}-\displaystyle \frac{be}{\gamma _{2}}I_2\nonumber \\&-\,\frac{bI^{*}_{2,s_{2}}}{I_2}E_2 +\displaystyle \frac{be}{\gamma _{2}}I^{*}_{2,s_{2}}\nonumber \\&+\,\frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})} \; f(S,I_1)I_1-\frac{ac}{\gamma _{1}}I_1, \end{aligned}$$
(3.26)

then

$$\begin{aligned}&\dot{{\mathcal {L}}}_2(S,E_{1},E_{2},I_{1},I_{2})= \;bE^{*}_{2,s_{2}} \left( 4-\frac{bE^{*}_{2,s_{2}}}{g(S,I^{*}_{2,s_{2}})I^{*}_{2,s_{2}}}-\frac{g(S,I_2)I_2}{bE_2}-\frac{I^{*}_{2,s_{2}}E_2}{I_2E^{*}_{2,s_{2}}} -\frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}\right) \nonumber \\&\qquad +\,bE^{*}_{2,s_{2}}\left( \frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}+\frac{g(S,I_2)}{g(S,I^{*}_{2,s_{2}})}\frac{I_2}{I^{*}_{2,s_{2}}}-\frac{I_2}{I^{*}_{2,s_{2}}}-1 \right) +\,\frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{12,s_{2}})} \; f(S,I_1)I_1-\frac{ac}{\gamma _{1}}I_1 \end{aligned}$$
(3.27)
$$\begin{aligned}&\dot{{\mathcal {L}}}_2(S,E_{1},E_{2},I_{1},I_{2})= \nonumber \;bE^{*}_{2,s_{2}} \left( 4-\frac{bE^{*}_{2,s_{2}}}{g(S,I^{*}_{2,s_{2}})I^{*}_{2,s_{2}}}-\frac{g(S,I_2)I_2}{bE_2}-\frac{I^{*}_{2,s_{2}}E_2}{I_2E^{*}_{2,s_{2}}} -\frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}\right) \nonumber \\&\qquad +\,bE^{*}_{2,s_{2}}\left( \frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}+\frac{g(S,I_2)}{g(S,I^{*}_{2,s_{2}})}\frac{I_2}{I^{*}_{2,s_{2}}}-\frac{I_2}{I^{*}_{2,s_{2}}}-1 \right) \nonumber \\ {}&\qquad +\, \frac{ac}{\gamma _1}I_1\left( \displaystyle \frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})} \; \displaystyle \frac{f(S,I_1)}{f(S_{0}^{*},0)}\; R_{0}^{1}-1\right) . \end{aligned}$$
(3.28)

From \((H_2)\) and \((H_3)\), we have \( \displaystyle \frac{f(S,I_1)}{f(S_{0}^{*},0)} \le 1,\) then,

$$\begin{aligned}&\dot{{\mathcal {L}}}_2(S,E_{1},E_{2},I_{1},I_{2})\le \;bE^{*}_{2,s_{2}} \left( 4-\frac{bE^{*}_{2,s_{2}}}{g(S,I^{*}_{2,s_{2}})I^{*}_{2,s_{2}}}-\frac{g(S,I_2)I_2}{bE_2}-\frac{I^{*}_{2,s_{2}}E_2}{I_2E^{*}_{2,s_{2}}} -\frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}\right) \nonumber \\&\qquad +\,bE^{*}_{2,s_{2}}\left( \frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}+\frac{g(S,I_2)}{g(S,I^{*}_{2,s_{2}})}\frac{I_2}{I^{*}_{2,s_{2}}}-\frac{I_2}{I^{*}_{2,s_{2}}}-1 \right) +\, \frac{ac}{\gamma _1}I_1 \left( \displaystyle \frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})} R_{0}^{1}-1\right) , \end{aligned}$$
(3.29)

From (3.22), we have

$$\begin{aligned}&\left( \frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}+\frac{g(S,I_2)}{g(S,I^{*}_{2,s_{2}})}\frac{I_2}{I^{*}_{2,s_{2}}}-\frac{I_2}{I^{*}_{2,s_{2}}}-1 \right) \nonumber \\&\quad = \left( 1- \frac{g(S,I_2)}{g(S,I^{*}_{2,s_{2}})}\right) \left( \frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}-\frac{I_2}{I^{*}_{2,s_{2}}} \right) \le \;\;0. \end{aligned}$$
(3.30)

By the relation between arithmetic and geometric means, we have

$$\begin{aligned} \left( 4-\frac{bE^{*}_{2,s_{2}}}{g(S,I^{*}_{2,s_{2}})I^{*}_{2,s_{2}}}-\frac{g(S,I_2)I_2}{bE_2}-\frac{I^{*}_{2,s_{2}}E_2}{I_2E^{*}_{2,s_{2}}} -\frac{g(S,I^{*}_{2,s_{2}})}{g(S,I_2)}\right) \;\le \; 0\nonumber \\ \end{aligned}$$
(3.31)

We discuss two cases:

  • If \(S_{2}^{*}\le S,\) then from \((H_2)\), we will have \(\displaystyle \frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})}\le 1,\) since \(\left( \displaystyle \frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})} R^{1}_0 -1\right) \le 0\), from \(R^{1}_0 \le 1,\) we will have \(\dot{{\mathcal {L}}}_2\le 0.\)

  • If \(S \le S_{2}^{*},\) then from \((H_2)\) we will obtain \(\displaystyle \frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})}\ge 1,\) from \(R^{1}_0 \le \displaystyle \frac{g(S,I^{*}_{2,s_{2}})}{g(S_{2}^{*},I^{*}_{2,s_{2}})} \le 1,\) we will get \(\displaystyle \frac{g(S_{2}^{*},I^{*}_{2,s_{2}})}{g(S,I^{*}_{2,s_{2}})} R^{1}_0 -1 \le 0,\) which implies, \(\dot{{\mathcal {L}}}_2\le 0.\) By the above discussion, we deduce that if \(R^{1}_0 \le 1,\) then \( \dot{{\mathcal {L}}}_2\le 0.\)

We conclude that the steady state \({\mathcal {E}}_{s_{2}}\) is globally asymptotically stable when \( R^{1}_0 \le 1\) and \(1 < R^{2}_0\).\(\square \)

3.4 Global stability of the third endemic equilibrium

For the global stability of \({\mathcal {E}}_{t}\), we assume that the functions f and g satisfy the following condition:

$$\begin{aligned}&\left( 1- \frac{g(S,I_2)}{g(S^{*}_{t},I^{*}_{2,t})} \frac{f(S^{*}_{t},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}\right) \left( \frac{g(S^{*}_{t},I^{*}_{2,t})}{g(S,I_2)} \frac{f(S,I^{*}_{1,t})}{f(S^{*}_{t},I^{*}_{1,t})}- \frac{I_2}{I^{*}_{2,s_{2}}} \right) \;\nonumber \\&\le 0, \;\; \forall \; S, I_1, I_2 > 0. \end{aligned}$$
(3.32)

Theorem 5

The endemic equilibrium \({\mathcal {E}}_{s_{t}}\) is globally asymptotically stable when \(R^{1}_0 > 1\) and \(R^{2}_0 > 1\).

Proof

We consider the Lyapunov function \({\mathcal {L}}_3\) defined by:

$$\begin{aligned} {\mathcal {L}}_3(S,E_{1},E_{2},I_{1},I_{2})=&\nonumber \;S-S_{t}^{*}-\int ^{S}_{S_{t}^{*}}\frac{f(S_{t}^{*},I^{*}_{1,t})}{f(X,I^{*}_{1,t})}\,\mathrm{d}X\\&+\,E^{*}_{1,t}\left( \frac{E_{1}}{E^{*}_{1,t}}-\ln \left( \frac{E_{1}}{E^{*}_{1,t}}\right) -1\right) \nonumber \\&+E^{*}_{2,t}\left( \frac{E_{2}}{E^{*}_{2,t}}-\ln \left( \frac{E_{2}}{E^{*}_{2,t}}\right) -1\right) \nonumber \\&+\,\frac{a}{\gamma _{1}}I^{*}_{1,t}\left( \frac{I_{1}}{I^{*}_{1,t}}-\ln \left( \frac{I_{1}}{I^{*}_{1,t}}\right) -1\right) \nonumber \\&+\,\frac{b}{\gamma _{2}}I^{*}_{2,t}\left( \frac{I_{2}}{I^{*}_{2,t}}-\ln \left( \frac{I_{2}}{I^{*}_{2,t}}\right) -1\right) , \end{aligned}$$
(3.33)

then

$$\begin{aligned} \dot{{\mathcal {L}}}_3(S,E_{1},E_{2},I_{1},I_{2})=&\nonumber \displaystyle {\left( 1-\frac{f(S_{t}^{*},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}\right) {\dot{S}}}\\&+\,\displaystyle {\left( 1-\frac{E^{*}_{1,t}}{E_{1}}\right) }\dot{E_{1}}+\displaystyle {\left( 1-\frac{E^{*}_{2,t}}{E_{2}}\right) }\dot{E_{2}}\nonumber \\&+\,\displaystyle {\frac{a}{\gamma _{1}}\left( 1-\frac{I^{*}_{1,t}}{I_{1}}\right) }\dot{I_{1}}+ \displaystyle {\frac{b}{\gamma _{2}}\left( 1-\frac{I^{*}_{2,t}}{I_{2}}\right) }\dot{I_{2}}. \end{aligned}$$
(3.34)

It is easy to verify that

$$\begin{aligned} \left\{ \begin{aligned}&\varLambda = \delta S^{*}_{t}+f(S_{t}^{*},I^{*}_{1,t})I^{*}_{1,t}+g(S_{t}^{*},I^{*}_{2,t})I^{*}_{2,t},\\&f(S_{t}^{*},I^{*}_{1,t})I^{*}_{1,t}= aE^{*}_{1,t} , \; g(S_{t}^{*},I^{*}_{2,t})I^{*}_{2,t}= bE^{*}_{2,t}, \\&\displaystyle \frac{E^{*}_{1,t}}{I^{*}_{1,t}}=\displaystyle \frac{c}{\gamma _{1}} \; , \; \displaystyle \frac{E^{*}_{2,t}}{I^{*}_{2,t}}=\displaystyle \frac{e}{\gamma _{2}}. \end{aligned} \right. \end{aligned}$$
(3.35)

As a result,

$$\begin{aligned}&\dot{{\mathcal {L}}}_3(S,E_{1},E_{2},I_{1},I_{2})= \nonumber \delta S^{*}_{t}\left( 1-\frac{S}{S_{t}^{*}} \right) \left( 1-\frac{f(S_{t}^{*},I^{*}_{1,t})}{f(S,I^{*}_{1,t})} \right) \nonumber \\&\quad +\, aE^{*}_{1,t}\left( 4-\frac{aE^{*}_{1,t}}{f(S,I^{*}_{1,t})I^{*}_{1,t}}- \frac{f(S,I_1)I_1}{aE_1}-\frac{I^{*}_{1,t}E_1}{I_1E^{*}_{1,t}} -\frac{f(S,I^{*}_{1,t})}{f(S,I_1)}\right) \nonumber \\&\qquad +\, \nonumber bE^{*}_{2,t} \left( 4-\frac{f(S^{*}_{t},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}- \frac{g(S,I_2)I_2}{bE_2}-\frac{I^{*}_{2,t} E_2}{ I_2 E^{*}_{2,t}} -\frac{b E^{*}_{2,t} f(S,I^{*}_{1,t})}{g(S,I_2)f(S^{*}_{t},I^{*}_{1,t})I^{*}_{2,t}}\right) \nonumber \\&\qquad +\, \nonumber aE^{*}_{1,t}\left( \frac{f(S,I^{*}_{1,t})}{f(S,I_1)}+\frac{f(S,I_1)}{f(S,I^{*}_{1,t})}\frac{I_1}{I^{*}_{1,t}}- \frac{I_1}{I^{*}_{1,t}}-1\right) \nonumber \\&\qquad +\, bE^{*}_{2,s_{2}}\left( \frac{g(S^{*}_{t},I^{*}_{2,t})}{g(S,I_2)} \frac{f(S,I^{*}_{1,t})}{f(S^{*}_{t},I^{*}_{1,t})}+ \frac{g(S,I_2)}{g(S^{*}_{t},I^{*}_{2,t})} \frac{f(S^{*}_{t},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}\frac{I_2}{I^{*}_{2,t}}-\frac{I_2}{I^{*}_{2,t}}-1\right) . \end{aligned}$$
(3.36)

Using the following trivial inequalities

$$\begin{aligned}&1- \frac{f(S_{t}^{*},I^{*}_{1,t})}{f(S,I^{*}_{1,t})} \ge 0\;\; \text {for} \;\; S \ge S_{t}^{*} \end{aligned}$$
(3.37)
$$\begin{aligned}&1- \frac{f(S_{t}^{*},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}< 0\;\; \text {for} \;\; S < S_{t}^{*}, \end{aligned}$$
(3.38)

then

$$\begin{aligned} \left( 1-\frac{S}{S_{t}^{*}} \right) \left( 1- \frac{f(S_{t}^{*},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}\right) \le 0, \end{aligned}$$
(3.39)

by the relation between arithmetic and geometric means, we have

$$\begin{aligned}&\left( 4-\frac{aE^{*}_{1,t}}{f(S,I^{*}_{1,t})I^{*}_{1,t}}- \frac{f(S,I_1)I_1}{aE_1}-\frac{I^{*}_{1,t}E_1}{I_1 E^{*}_{1,t}} -\frac{f(S,I^{*}_{1,t})}{f(S,I_1)}\right) \le 0 \end{aligned}$$
(3.40)
$$\begin{aligned}&\left( 4-\frac{f(S^{*}_{t},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}- \frac{g(S,I_2)I_2}{bE_2}-\frac{I^{*}_{2,t} E_2}{ I_2 E^{*}_{2,t}} -\frac{b E^{*}_{2,t} f(S,I^{*}_{1,t})}{g(S,I_2)f(S^{*}_{t},I^{*}_{1,t})I^{*}_{2,t}}\right) \le 0, \end{aligned}$$
(3.41)

from (3.12) we have

$$\begin{aligned}&\left( \frac{f(S,I^{*}_{1,t})}{f(S,I_1)}+\frac{f(S,I_1)}{f(S,I^{*}_{1,t})}\frac{I_1}{I^{*}_{1,t}}-\frac{I_1}{I^{*}_{1,t}}-1 \right) \nonumber \\&\quad = \left( 1- \frac{f(S,I_1)}{f(S,I^{*}_{1,t})}\right) \left( \frac{f(S,I^{*}_{1,t})}{f(S,I_1)}-\frac{I_1}{I^{*}_{1,t}} \right) \le \;\;0, \end{aligned}$$
(3.42)

also from (3.32) we have

$$\begin{aligned}&\left( \frac{g(S^{*}_{t},I^{*}_{2,t})}{g(S,I_2)} \frac{f(S,I^{*}_{1,t})}{f(S^{*}_{t},I^{*}_{1,t})}{+} \frac{g(S,I_2)}{g(S^{*}_{t},I^{*}_{2,t})} \frac{f(S^{*}_{t},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}\frac{I_2}{I^{*}_{2,t}}-1 -\frac{I_2}{I^{*}_{2,t}}\right) \nonumber \\&\quad =\left( 1-\varGamma \right) \left( \frac{1}{\varGamma } -\frac{I_2}{I^{*}_{2,t}} \right) \;\le 0, \end{aligned}$$
(3.43)

with \(\varGamma = \displaystyle \frac{g(S,I_2)}{g(S^{*}_{t},I^{*}_{2,t})} \displaystyle \frac{f(S^{*}_{t},I^{*}_{1,t})}{f(S,I^{*}_{1,t})}\)

We conclude that the steady state \({\mathcal {E}}_{t}\) is globally asymptotically stable when \( 1 < R_0^{1} \) and \( 1 < R^{2}_0 \). \(\square \)

4 Numerical simulations

In this section, we will perform some numerical simulations in order to examine numerically the SEIR infection dynamics under various incidence functions and also validating our theoretical results. Indeed, we will restrict ourselves to four cases, the first one is to consider the model (1.1) along with the simplest two bilinear incidence functions \(f(S,I_1)=\alpha S\) and \(g(S,I_2)=\beta S\), the second case deals with two Beddington–DeAngelis incidence functions, \(f(S,I_1)=\displaystyle \frac{\alpha S}{1+ \omega _{1} S+\omega _{2} I_1}\) and \(g(S,I_2)= \displaystyle \frac{\beta S}{1+ \omega _{3} S+\omega _{4} I_2}\), while the third case is devoted to check the impact of choosing both incidence functions under Crowley–Martin incidence form, \(f(S,I_1)=\displaystyle \frac{\alpha S}{1+ \chi _{1} S+\chi _{2} I_1+ \chi _{1} \chi _{2} S I_1}\) and \(g(S,I_2)= \displaystyle \frac{\beta S}{1+ \chi _{3} S+\chi _{4} I_2+ \chi _{3}\chi _{4} S I_2}\). The last case consists of incorporating into the model two non-monotonic incidence functions \(f(S,I_1)=\displaystyle {\frac{\alpha S}{1+ \alpha _1 I_1^{2}}}\) and \(g(S,I_2)=\displaystyle {\frac{\beta S}{1+ \alpha _2 I_2^{2}}}\).

4.1 The stability of the disease-free equilibrium

The disease-free equilibrium is characterized by the extinction of the infection, and the disease cannot invade the population. In this subsection, attention will be focused to the numerical stability of such disease-free equilibrium. Indeed, Theorem 2 points out that we can expect the stability of disease-free equilibrium when the basic reproduction number is less than unity. This permits us to look for the right model parameters in order to check numerically the stability of the first steady sate.

Fig. 2
figure 2

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the disease-free equilibrium \({\mathcal {E}}_{f}\). The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

Full size image

Figure 2 shows the behaviour of the infection for the following parameter values: \(\varLambda = 1\), \(\alpha = 0.17\), \(\beta = 0.15\), \(\gamma _1 = 0.3\), \(\gamma _2 = 0.4\), \(\mu _1 = 0.65\), \(\mu _2 = 0.75\), \(\delta = 0.2\), \(\omega _1=0.4\), \(\omega _2=0.6\), \(\omega _3=4.5\), \(\omega _4=5.8\), \(\chi _1=2.5\), \(\chi _2=3\), \(\chi _3=6.5\), \(\chi _4=5\), \(\alpha _1= 2.5\) and \(\alpha _2= 3\). We clearly see that the solutions of the system, under the various suggested incidence functions (bilinear, Beddington–DeAngelis, Crowley–Martin and non-monotonic functions), converge towards the same disease-free equilibrium point \({\mathcal {E}}_{f}=(5,0,0,0,0,0)\). In this situation, the nature of the incidence function has no effect on the steady-state value. Consequently, it was remarked that the disease-free equilibrium first nonzero component depends only on the birth and death rates of the susceptible individuals and does not depend, in any way, on the incidence functions parameters. Here, the disease dies out, the susceptible reaches their maximum value, and the other variables vanish. With the chosen parameters, we can easily calculate the two- strain basic reproduction numbers; in our case, we will have \(R_0^{1}=0.6\) and \(R_0^{2}=0.5263\) for the bilinear incidence functions case; we will have also \(R_0^{1}=0.1\) and \(R_0^{2}=0.6579\) for the Beddington–DeAngelis incidences case; \(R_0^{1}=0.1154\) and \(R_0^{2}=0.1645\) are calculated for Crowley–Martin incidence functions case, and finally, we have \(R_0^{1}=0.6\) and \(R_0^{2}=0.5263\) for the non-monotonic incidence functions case. In all these cases, we have the basic reproduction number \(R_0\) is less than unity which is consistent with the theoretical result concerning the stability of the disease-free equilibrium \({\mathcal {E}}_{f}\).

Fig. 3
figure 3

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the strain 1 endemic equilibrium \({\mathcal {E}}_{s_{1}}\). The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

Full size image

4.2 The stability of the endemic equilibria

The dynamics behaviour concerning the stability of the strain 1 endemic equilibrium \({\mathcal {E}}_{s_{1}}\) is shown in Fig. 3. Indeed, we have chosen the following parameter values: \(\varLambda = 1\), \(\alpha = 0.5\), \(\beta = 0.12\), \(\gamma _1 = 0.4\), \(\gamma _2 = 0.3\), \(\mu _1 = 0.4\), \(\mu _2 = 0.75\), \(\delta = 0.2\), \(\omega _1=0.4\), \(\omega _2=0.85\), \(\omega _3=1.5\), \(\omega _4=0.05\), \(\chi _1=4\), \(\chi _2=3.5\), \(\chi _3=0.3\), \(\chi _4=0.05\), \(\alpha _1=2\) and \(\alpha _2=1.5\). For the bilinear incidence functions case, the solution converges towards the strain 1 endemic equilibrium (1.8, 1.0667, 0, 0.7111, 0, 1.4222) and with the chosen parameters we have \(R_0^{1}=2.7778\) and \(R_0^{2}=0.3789\). For the model with Beddington–DeAngelis, we observe the convergence towards \((3.1788, 0.6071, 0, 0.4047, 0, 0.8094)\) and we have \(R_0^{1}=11.1111\) and \(R_0^{2}=0.3609\); for the third case, Crowley–Martin incidence type, we observe the convergence towards the equilibrium (2.3664, 0.8779, 0, 0.5852, 0, 1.1705) and we have \(R_0^{1}=11.1111\) and \(R_0^{2}=0.1024\). For the last non-monotonic incidence functions case, we see in the figure the convergence towards the steady state (2.7223, 0.7592, 0, 0.5062, 0, 1.0123) and we have \(R_0^{1}=2.7778\) and \(R_0^{2}=0.3789\). We remark that, for the all four cases, the strain 1 basic reproduction number \(R_0^{1}\) is greater than unity while the other strain 2 basic reproduction number \(R_0^{2}\) is less than one which confirms our theoretical findings concerning the stability of the strain 1 endemic equilibrium \({\mathcal {E}}_{s_{1}}\). This endemic equilibrium is characterized by the vanishing the strain 2 latent and infectious individuals.

Fig. 4
figure 4

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the strain 2 endemic equilibrium \({\mathcal {E}}_{s_{2}}\). The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

Full size image
Fig. 5
figure 5

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the endemic equilibrium \({\mathcal {E}}_{t}\). The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

Full size image

The strain 2 endemic equilibrium \({\mathcal {E}}_{s_{2}}\) stability is illustrated in Fig. 4. In this figure, the following parameters are chosen: \(\varLambda = 1\), \(\alpha = 0.2\), \(\beta = 0.4\), \(\gamma _1 = 0.4\), \(\gamma _2 = 0.3\), \(\mu _1 = 1\), \(\mu _2 = 0.5\), \(\delta = 0.2\), \(\omega _1=0.4\), \(\omega _2=0.85\), \(\omega _3=1.5\), \(\omega _4=0.05\), \(\chi _1=4\), \(\chi _2=3.5\), \(\chi _3=0.3\), \(\chi _4=0.05\), \(\alpha _1=2\) and \(\alpha _2=1.5\). We observe the convergence of the solution towards the steady state (2.9167, 0, 0.8333, 0, 0.3571, 0.8928) for the bilinear incidence functions case and with the adopted parameters we have \(R_0^{1}=0.5556\) and \(R_0^{2}=1.7143\). For the model with Beddington–DeAngelis, we observe the convergence towards \((4.1559, 0, 0.3376, 0, 0.1447, 0.3617)\) and we have \(R_0^{1}=0.5291\) and \(R_0^{2}=6.8571\); for the third case, Crowley–Martin incidence type, we observe the convergence towards the equilibrium (3.6876, 0, 0.5250, 0, 0.2250, 0.5624) and we have \(R_0^{1}=0.1502\) and \(R_0^{2}=6.8571\). For the last non-monotonic incidence functions case, we see in the figure the convergence towards the steady state (3.2918, 0, 0.6833, 0, 0.2928, 0.7321) and we have \(R_0^{1}=0.5556\) and \(R_0^{2}=1.7143\). We remark that, for all the four cases, the strain 1 basic reproduction number \(R_0^{1}\) is less than unity while the other strain 2 basic reproduction number \(R_0^{2}\) is greater than one which is in good agreement with our theoretical findings concerning the stability of the strain 2 endemic equilibrium \({\mathcal {E}}_{s_{2}}\). From the components of this endemic equilibrium, we remark the clearance of the strain 1 latent and infectious individuals.

The last endemic equilibrium \({\mathcal {E}}_{t}\) stability behaviour is depicted in Fig. 5 for \(\varLambda = 1\), \(\alpha = 0.6\), \(\beta = 0.6\), \(\gamma _1 = 0.5\), \(\gamma _2 = 0.5\), \(\mu _1 = 0.15\), \(\mu _2 = 0.15\), \(\delta = 0.2\), \(\omega _1 =1.2\), \(\omega _2=0.06\), \(\omega _3=1.2\), \(\omega _4=0.06\), \(\chi _1=0.4\), \(\chi _2=0.05\), \(\chi _3=0.4\), \(\chi _4=0.05\), \(\alpha _1=0.14\) and \(\alpha _2=0.145\). As the previous figures, the dynamics for the different suggested incidence functions are illustrated. For the first one, the bilinear incidence function, we observe the convergence towards (0.8167, 0.5196, 0.6757, 0.7422, 0.9652, 1.2806) and with the adopted parameters we have \(R_0^{1}=6.1224\) and \(R_0^{2}=6.1224\). For Beddington–DeAngelis case, we observe the convergence towards \((1.5783, 0.4888, 0.4888, 0.6983, 0.6983, 1.0474)\) and we have \(R_0^{1}=23.5479\) and \(R_0^{2}=23.5479\), for the third case, Crowley–Martin incidence type, we observe the convergence towards the equilibrium \((1.1353, 0.5519, 0.5523, 0.7884, 0.7890, 1.1831)\) and we have \(R_0^{1}=24.4898\) and \(R_0^{2}=24.4898\). About the last non-monotonic incidence functions case, we see the convergence towards the steady state \((0.8982, 0.5904, 0.5815, 0.8433, 0.8309, 1.2557)\) and we have \(R_0^{1}=6.1224\) and \(R_0^{2}=6.1224\). We remark that both the two-strain basic reproduction numbers \(R_0^{1}\) and \(R_0^{2}\) are greater than unity which confirms our theoretical findings concerning the stability of the last endemic equilibrium \({\mathcal {E}}_{t}\). This last endemic equilibrium is characterized by the persistence of all strains latent and infectious individuals. The numerical simulations confirm that the model with general incidence functions encompasses a large number of classical well-known incidence functions. Therefore, the generalized mathematical model can give a wide view about the stability of the different problem equilibria.

4.3 Comparison with COVID-19 clinical data

As we have mentioned in the introduction, the recent pandemic COVID-19 is a multi-strain infection. Therefore, the main interest of this subsection is to compare the numerical simulations resulting from our multi-strain epidemic model with COVID-19 clinical data. We have chosen to make our comparison the Moroccan clinical data during the year 2020 in the period between March 31 and June 20 [39, 40].

Figure 6 shows the time evolution of infected cases for the following parameter values: \(\varLambda = 1.5\), \(\alpha = 1.67\), \(\beta = 0.88\), \(\delta = 0.2\), \(\gamma _1 = 1.05\), \(\gamma _ 2 = 6.8\), \(\mu _1 = 0.005\), \(\mu _2 = 0.087\), \(\omega _1=0.3\), \(\omega _2=0.5\), \(\omega _3=0.1\), \(\omega _4=0.3\), \(\chi _1=0.01\), \(\chi _2=0.01\), \(\chi _3= 2.5\), \(\chi _4= 2.8\), \(\alpha _1=0.0005\) and \(\alpha _2=0.007\). We observe that a significant relationship exists between the curve representing the COVID-19 clinical data and the numerical simulations resulting from our mathematical model for the different incidence functions. Indeed, a good fit between the infected cases given by the mathematical model and the clinical data is observed. Each numerical result with the incidence function (bilinear, Beddington–DeAngelis, Crowley–Martin or non-monotonic function) fits well the clinical data, especially for a certain period of observation. Hence, the multi-strain mathematical model with generalized incidence functions is more appropriate to represent the studied disease.

Fig. 6
figure 6

Time evolution of infected cases with the bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). The clinical infected cases are illustrated by magenta circles. (Color figure online)

Full size image
Fig. 7
figure 7

Time evolution of infected cases with the bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). The clinical infected cases are illustrated by magenta circles. (Color figure online)

Full size image

Figure 7 (left-hand side) shows the time evolution of infected cases for the following parameter values: \(\varLambda = 2\), \(\alpha = 0.01\), \(\beta = 0.01\), \(\gamma _1 = 0.3\), \(\gamma _2 = 0.4\), \(\mu _1 = 0.07\), \(\mu _2 = 0.5\), \(\delta = 0.2\), \(\omega _1=0.5\), \(\omega _2=0.6\), \(\omega _3=0.1\), \(\omega _4=0.3\), \(\chi _1=1.5\), \(\chi _2=2\), \(\chi _3=0.7\), \(\chi _4=3.7\), \(\alpha _1= 2\) and \(\alpha _2= 1.5\). We see that the solution of our model, under the various suggested incidence functions, converges towards the same disease-free equilibrium point \({\mathcal {E}}_{f}\). In this situation, the disease dies out, the susceptible reaches their maximum value and the other variables vanish. Within the chosen parameters, we can easily calculate the basic reproduction number; in our case we will have \(R_0 = 0.23\) for the bilinear incidence functions case; we will have also \(R_0 = 0.28\) for the Beddington–DeAngelis incidences case; \(R_0 = 0.2\) are calculated for Crowley–Martin incidence functions case, and finally, we have \(R_0 = 0.23\) for the non-monotonic incidence functions case.

The COVID-19 disease persistence case is illustrated in Fig. 7 (right-hand side). In this figure, the following parameters are chosen: \(\varLambda = 1.2\), \(\alpha = 0.1\), \(\beta = 0.52\), \(\gamma _1 = 0.4\), \(\gamma _2 = 0.1\), \(\mu _1 = 0.23\), \(\mu _2 = 0.2\), \(\delta = 0.03\), \(\omega _1=0.5\), \(\omega _2=0.6\), \(\omega _3=0.1\), \(\omega _4=0.3\), \(\chi _1=1.5\), \(\chi _2=4.5\), \(\chi _3=0.7\), \(\chi _4=0.12\), \(\alpha _1=2\) and \(\alpha _2=1.5\). We remark the convergence of the solution towards the endemic equilibrium for all the taken incidence functions. Indeed, for the case with bilinear incidence functions, we have the basic reproduction number is greater than unity \(R_0=14.31\). For other cases, the basic reproduction number is also greater than one; indeed for Beddington–DeAngelis case, we have \(R_0= 36.6956\); for the third case, Crowley–Martin incidence type, we have \(R_0= 14.3113\). For the last non-monotonic incidence functions case, we have \(R_0=14.3113\) which means that the disease persists. We can conclude that the numerical simulations fit well COVID-19 clinical data. Our numerical simulations reveal two scenarios of evolution for this pandemic, and within the first scenario the disease will die out. The other scenario happens when the basic reproduction number is greater than unity; in this case the disease will persist. In this situation, it will be important to eventually undertake some strategies like quarantine, isolation, wearing of masks, disinfection and if it will be possible vaccination.

Fig. 8
figure 8

Effect of quarantine strategy on the SEIR model dynamics

Full size image

4.4 The effect of quarantine strategy

In this subsection, we will study the effect of quarantine strategy in controlling the infection spread. To illustrate this effect and for simplicity, we will restrict ourselves to the case of the bilinear incidence function. More precisely, we will take the following incidence forms:

$$\begin{aligned} f(S,I_1) {=} \alpha (1-u_1) S \; \; \text {and} \; \; g(S,I_2)= \beta (1-u_2) S,\nonumber \\ \end{aligned}$$
(4.1)

where \(u_1\) stands for the efficiency of the quarantine strategy concerning the first strain infection rate, while \(u_2\) represents the efficacy of the quarantine measures concerning the second-strain infection rate. Our numerical simulations will be performed in order to check the impact of the quarantine strategy for the case of the last endemic equilibrium.

Figure 8 shows the time evolution of the susceptible, the two- strain latent individuals, two-strain infected individuals and the recovered for the following parameter values: \(\varLambda = 1\), \(\alpha = 0.6\), \(\beta = 0.6\), \(\gamma _1 = 0.5\), \(\gamma _2 = 0.5\), \(\mu _1 = 0.15\), \(\mu _2 = 0.15\), \(\delta = 0.2\) and for different values of the strategy controls \(u_1\) and \(u_2\). When no strategy is applied, we find the same result as in Fig. 5. By increasing the efficiency of the quarantine measures, we observe an interesting result. Indeed, the number of the susceptible individuals increases when \(u_1\) and \(u_2\) increase. For higher values of these two control parameters, the number of two-strain latent and infected individuals is reduced considerably, which means that the quarantine strategy can reduce the infection in efficient manner.

5 Conclusion

In this paper, we have studied the global stability of two-strain epidemic model with two general incidence functions. The model included six compartments, namely the susceptible, two categories of the exposed, two categories of the infected and the removed individuals, this kind of model takes the abbreviation SEIR. We have established the existence, positivity and boundedness of solutions results which guarantee the wellposedness of our SEIR model. The disease-free equilibrium, the endemic equilibrium with respect to strain 1, the endemic equilibrium with respect to strain 2 and the endemic equilibrium with respect to both strains are given. By using an appropriate Lyapunov functionals, the global stability of the equilibria is established depending on the basic reproduction number \(R_0\), the strain 1 reproduction number \(R^{1}_0\) and the strain 2 reproduction number \(R^{2}_0\). Numerical simulations are performed in order to confirm our different theoretical results. It was observed that the model with a generalized incidence function encompasses a large number of classical incidence functions and it can give more clear view about the equilibria stability. In addition, a numerical comparison between our model results and COVID-19 clinical data is conducted. We have observed a good fit between our numerical simulations and the clinical data which indicates that our multi-strain mathematical model can fit and predict the evolution of the infection. We have observed that strict quarantine measures can reduce significantly the infection spread. It is worthy to notice that a generalization form of this problem to a more complex compartmental model is proposed in Appendix of this paper. For this complex model, we have given its basic reproduction number and each strain reproduction number. Furthermore, the forms of the Lyapunov functionals for this complex compartmental model are also formulated.