1 Introduction

Mathematical modeling has been extensively used in last few decades for studying the dynamics of the infectious diseases. After the pioneering work of Kermack and McKendrick (1927), different kinds of mathematical models have been proposed in order to achieve a better understanding of the transmission pattern of the infectious diseases (see Dubey et al. (2016a), Dubey et al. (2016b), Hattaf et al. (2012), Khatua and Kar (2020a), Khatua and Kar (2020b), Khatua et al. (2021), Mandal et al. (2020), Mandal et al. (2020) ). The proper choice of the incidence function is of great practical importance in studying the disease dynamics. There are several nonlinear incidence functions in literature for describing the disease transmission rate, for example bi-linear incidence rate (Martcheva (2015)), nonlinear incidence rate (Korobeinikov and Maini (2004)), saturated incidence rate (Zhang and Liu (2008)), non-monotone incidence rate (Xiao and Ruan (2007)), and general incidence rate (Ávila-Vales et al. (2017)). However, from a practical point of view, out of these, the non-monotone incidence is most important as it is capable to incorporate the behavioral changes of the susceptible population when the infected population is very large. Xiao and Ruan (2007) studied the global behavior of an epidemic system with non-monotonic incidence rate, while Muroya et al. (2011) investigated a delayed SIRS type epidemic model with a non-monotonic incidence rate. By using the monotone iterative techniques, they established the sufficient conditions for the global asymptotic stability of endemic equilibrium under the condition \(R_0>1\). Baba and Hincal (2017) investigated the dynamics of a two strain epidemic model with bi-linear and non-monotone incidence rate. They have shown that any strain with highest basic reproduction ratio will automatically outperform the other strain. Meskaf et al. (2020) studied the two-strain SEIR type model with non-monotone incidence rates. Very recently, Kumar (2021) studied the effects of non-monotonic functional response on a disease transmission model. Some other model-based research works using non-monotone incidence rate can be found in Feng and Qiu (2018), Li and Teng (2018), Kar and Batabyal (2010) and the references cited therein.

Also, there are many diseases in which more than one pathogen strains are observed due to the process of virus mutation, for example influenza, dengue, COVID-19, and some sexually transmitted diseases. So, it is better to use the multi-strain models as they can describe those diseases more accurately in comparison with the single strain models. Baba et al. (2018) studied an epidemic model consisting of two strains with vaccine for each strain. They concluded that if both the basic reproduction numbers are less than unity, then the disease dies out, and if one of the ratios is less than unity, then epidemic occurs with respect to the other. While Baba and Hincal (2018) proposed and analyzed an epidemic model for influenza with three strains and they studied the effects of vaccination and awareness, Nandi et al. (2019) studied the complex dynamical behavior of an epidemic model with two diseases, and also they pointed out the optimal control strategy. Similar kind of research articles on multi-strain models can be found in Li et al. (2013), Liu et al. (2018), Meskaf et al. (2020), and the references cited in those articles.

Moreover, spatial heterogeneity of the populations plays an important role in the spread of infectious diseases. Ignoring spatial diffusion in model is to ignore the reality. In this regard, reaction–diffusion equations are used to include the spatial movement of the populations into the models. Wang et al. (2012) investigated the complex dynamics of a reaction–diffusion epidemic model. Their results indicate that the spread of the epidemic is strongly influenced by the diffusion. Li and Bie (2019) studied the dynamics of a diffusive SIRS type epidemic system in spatially heterogeneous environment. They concluded that the spontaneous infection enhances the persistence of disease. Han and Lei (2019) studied the global dynamics of a diffusive SEIR epidemic model with nonlinear incidence. There are many research works investigating the role of spatial diffusion on the dynamics of the diseases, for example Holmes et al. (1994), Lotfi et al. (2014), Kuniya and Wang (2018), and the references cited therein. However, to our best knowledge no attempt has been made to study the dynamics of a two-strain reaction–diffusion epidemic model with non-monotone incidence rate for both the strains and virus mutation. In this present study, we will study the dynamical behavior of a diffusive two-strain epidemic model incorporating non-monotonic incidence rate.

The rest of this article is organized as follows: In Sect. 2, we formulate the two-strain epidemic model considering the spatial diffusion and virus mutation. Then, the existence, positivity and uniform boundedness of the solutions of the model system are studied in Sect. 3. All the possible equilibrium points of the system and their respective global stability analysis are investigated in the Sect. 4. Finally, in Sect. 5, some conclusions are provided.

2 Model Formulation

Here, we develop a two-strain diffusive epidemic model with virus mutation and non-monotone incidence rate. We assume that the total population is divided into four classes, namely susceptible class S(t), infected class with strain-1 \(I_{1}(t)\), infected class with strain-2 \(I_{2}(t)\) and recovered class R(t). Further we consider the following assumptions to construct the model.

\(\bullet\) The recruitment rate to the susceptible class is taken as \(\varLambda .\)

\(\bullet\) The susceptible individuals can become infected by strain-1 or strain-2 and move to the infected compartments \(I_{1}\) and \(I_{2}\) at a rate \(\alpha _{1}\) and \(\alpha _{2}\) respectively.

\(\bullet\) Following Xiao and Ruan (2007), non-monotonic incidence rate \(\frac{kSI}{1+\alpha I^{2}}\) is considered. The incidence function is increasing when I is small and decreasing when I is large. This type of incidence is considered to interpret the psychological effect: for a very large number of infective individuals the infection force may decrease as the number of infective individuals increases, because in the presence of large number of infectives the population may tend to reduce the number of contacts per unit time.

  • The infected individuals with one of the strain cannot be infected with the other, i.e., there is no super-infection. However, it is considered that the strain-1 mutate into strain-2 at a rate \(\xi\).

  • The infected individuals of \(I_{1}\) and \(I_{2}\) compartments recover at a rate \(\rho _{1}\) and \(\rho _{2}\) respectively.

  • The natural mortality rate for each classes is taken as d and the disease induced mortality rate for due to strain-1 and strain-2 is taken as \(\mu _{1}\) and \(\mu _{2}\) respectively.

  • Moreover, the spatial diffusion is taken into consideration, i.e., the population move randomly described as Brownian random motion.

Based on the above assumptions, our reaction–diffusion model takes the following form:

$$\begin{aligned} \begin{aligned}\frac{\partial S}{\partial t}&= \varLambda - \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}\\&\quad -\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-dS+d_{1}\nabla ^{2}S,\;\;\;x\in \varOmega , t>0,\\ \frac{\partial I_{1}}{\partial t}&=\frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\left( d+\mu _{1}+\rho _{1} \right) I_{1}\\ {}&\quad -\xi I_{1}+d_{2}\nabla ^{2}I_{1} ,\;\;\;x\in \varOmega , t>0,\\&\frac{\partial I_{2}}{ \partial t}=\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}\\ {}&\quad -\left( d+\mu _{2}+\rho _{2} \right) I_{2}+\xi I_{1}+d_{3}\nabla ^{2}I_{2},\;\;\;x\in \varOmega , t>0,\\\frac{\partial R}{\partial t} &=\rho _{1}I_{1}+\rho _{2}I_{2}-dR+d_{4}\nabla ^{2}R,\;\;\;x\in \varOmega , t>0, \end{aligned} \end{aligned}$$
(1)

with homogeneous Neumann boundary conditions

$$\begin{aligned} \frac{\partial S}{\partial \nu }=0,\frac{\partial I_{1}}{\partial \nu }=0,\frac{\partial I_{2}}{\partial \nu }=0,\frac{\partial R}{\partial \nu }=0,\;\;\;x\in \partial \varOmega , t>0 \end{aligned}$$
(2)

and initial conditions

$$\begin{aligned} S(x,0)= & {} S_{0}(x)\ge 0,I_{1}(x,0)=I_{10}(x)\ge \nonumber \\ 0,I_{2}(x,0)= & {} I_{20}(x)\ge 0,R(x,0)=R_{0}(x)\ge 0. \end{aligned}$$
(3)

Here, \(\varOmega\) is a bounded domain in \(R^n\) with smooth boundary \(\partial \varOmega\), and \(\frac{\partial S}{\partial \nu }\), \(\frac{\partial I_{1}}{\partial \nu }\), \(\frac{\partial I_{2}}{\partial \nu }\), and \(\frac{\partial R}{\partial \nu }\) are, respectively, the normal derivatives of \(S, I_1, I_2\) and R on \(\partial \varOmega\).

In the above, \(S(x, t), I_{1}(x, t), I_{2}(x, t), R(x, t)\), respectively, stand for the density of susceptible, infected with strain-1, infected with strain-2 and recovered populations at location x and time t. The positive constants \(d_{1}, d_{2}, d_{3}, d_{4}\) are used to measure their respective migration (movement) mobility.

3 Global Existence, Positivity, and Boundedness

Proposition 1

The system (1) has a unique solution defined on \([0, \; \infty )\) and this solution remains non-negative for all \(t\ge 0\).

Proof

The system (1) can be written abstractly in the Banach space \(X=C(\bar{\varOmega })\times C(\bar{\varOmega })\) of the following form:

$$\begin{aligned} u^{'}(t)= & Au(t)+F(u(t)),\;t>0,\\ u(0)= & x_{0}\in X. \end{aligned}$$

where \(u=Col(S,I_{1},I_{2},R)\), \(A(u(t))=Col(d_{1}\nabla ^{2}S,d_{2}\nabla ^{2}I_{1},d_{3}\nabla ^{2}I_{2},d_{4}\nabla ^{2}R)\) and

$$\begin{aligned} F(u(t))=\begin{pmatrix} \varLambda - \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}} -dS\\ \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}} -\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\xi I_{1}\\ \frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-\left( d+\mu _{2}+\rho _{2} \right) I_{2}+\xi I_{1}\\ \rho _{1}I_{1}+\rho _{2}I_{2}-dR \end{pmatrix}. \end{aligned}$$

However, F is locally Lipschitz in X. Therefore, it can be concluded from Pazy (1983) that the system (1) admits a local solution on \([0, \; T_{max}]\), where \(T_{max}\) is the maximum time for which the solution of the system (1) exist.

In addition, the system (1) can be written as:

$$\begin{aligned} \begin{aligned}\frac{\partial S}{\partial t}-d_{1}\nabla ^{2}S&= \varLambda - \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}} -dS\\&=F_{1}(S,I_{1},I_{2},R),\\&\frac{\partial I_{1}}{\partial t}-d_{2}\nabla ^{2}I_{1}\\&=\frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}} -\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\xi I_{1}=F_{2}(S,I_{1},I_{2},R),\\\frac{\partial I_{2}}{ \partial t}-d_{3}\nabla ^{2}I_{2}&=\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}} -\left( d+\mu _{2}+\rho _{2} \right) I_{2}+\xi I_{1}\\& =F_{3}(S,I_{1},I_{2},R),\\&\frac{\partial R}{\partial t}-d_{4}\nabla ^{2}R =\rho _{1}I_{1}+\rho _{2}I_{2}-dR\\&=F_{4}(S,I_{1},I_{2},R) . \end{aligned} \end{aligned}$$
(4)

Here the functions \(F_{i}(S,I_{1},I_{2},R),i=1,2,3,4\) are continuously differentiable functions and also

$$\begin{aligned} F_{1}(0,I_{1},I_{2},R)= & {} \varLambda >0,\\ F_{2}(S,0,I_{2},R)= & {} 0\ge 0,\\ F_{3}(S,I_{1},0,R)= & {} \xi _{1}I_{1}\ge 0,\\ F_{4}(S,I_{1},I_{2},0)= & {} \rho _{1}I_{1}+\rho _{2}I_{2}\ge 0. \end{aligned}$$

Since the initial data of the system (1) are non-negative, all the solutions of the system (1) remain positive. \(\square\)

Proposition 2

There exists a positive constant \(K_{1}\) independent of initial data such that the solution \((S,I_{1},I_{2},R)\) of the system (1) satisfies

$$\begin{aligned}&\mid \mid S(.,t)\mid \mid _{L^{\infty }(\varOmega )}+\mid \mid I_{1}(.,t)\mid \mid _{L^{\infty }(\varOmega )} +\mid \mid I_{2}(.,t)\mid \mid _{L^{\infty }(\varOmega )}\\&\quad + \mid \mid R(.,t)\mid \mid _{L^{\infty }(\varOmega )}\le K_{1}, \forall t\ge T_0, \end{aligned}$$

for some large time \(T_0 > 0\).

Proof

First adding the equations of the system (1), we obtain

$$\begin{aligned} \begin{aligned}&\frac{\partial }{\partial t}\left( S+I_{1}+I_{2}+R\right) \\ {}&\quad -d_{1} \nabla ^{2}S- d_{2}\nabla ^{2}I_{1}-d_{3}\nabla ^{2}I_{2}-d_{4}\nabla ^{2}R\\&= \varLambda -d\left( S+I_{1}+I_{2}+R\right) -\mu _{1}I_{1} \\ {}&\quad - \mu _{2}I_{2}, \;\;\forall t\in (0,T_{max}) \end{aligned} \end{aligned}$$
(5)

Now integrating both sides of Eq. (5) over \(\varOmega\) with respect to x by parts and using the zero Neumann boundary conditions, we then obtain

$$\begin{aligned} \begin{aligned}&\frac{d}{dt}\int _{\varOmega }\left( S+I_{1}+I_{2}+R\right) dx\\&\quad =\int _{\varOmega } \varLambda dx-\int _{\varOmega }d \left( S+I_{1}+I_{2}+R\right) dx\\&\qquad -\int _{\varOmega } \mu _{1}I_{1} dx-\int _{\varOmega }\mu _{2}I_{2} dx\\&\quad \le \int _{\varOmega }\varLambda dx-\int _{\varOmega }d \left( S+I_{1}+I_{2}+R\right) dx \end{aligned} \end{aligned}$$
(6)
$$\begin{aligned} \text {Let}\;\;V(t)=\int _{\varOmega }\left( S+I_{1}+I_{2}+R\right) dx. \end{aligned}$$

Then, from inequality (6), we have

$$\begin{aligned} \frac{dV}{dt}\le \mid \varOmega \mid \varLambda -d V\left( t\right) \end{aligned}$$

Now integrating the above inequality from 0 to t, we get

$$\begin{aligned} V(t)\le & {} V(0)e^{-dt}+\frac{\mid \varOmega \mid \varLambda }{d}\left( 1-e^{-dt}\right) \\&\text {i.e.}\int _{\varOmega }\left( S+I_{1}+I_{2}+R\right) dx\\\le & {} e^{-dt}\int _{\varOmega } \left( S_{0}(x)+I_{10}(x)+I_{20}(x)+R(x)\right) dx\\&+\frac{\mid \varOmega \mid \varLambda }{d}\left( 1-e^{-dt}\right) \forall t\ge 0 \end{aligned}$$

Thus, we obtain

$$\begin{aligned} \underset{t \rightarrow \infty }{lim\; sup}\int _{\varOmega } \left( S+I_{1}+I_{2}+R\right) dx \le \frac{\mid \varOmega \mid \varLambda }{d}. \end{aligned}$$
(7)

Now in view of [Theorem 3.3 of Dung (1998) ], there exists a positive constant \(K_{1}\), which does not depend on initial data such that

$$\begin{aligned}&\mid \mid S(.,t)\mid \mid _{L^{\infty }(\varOmega )}+\mid \mid I_{1}(.,t)\mid \mid _{L^{\infty }(\varOmega )}\\&\quad +\mid \mid I_{2}(.,t)\mid \mid _{L^{\infty }(\varOmega )}+\mid \mid R(.,t)\mid \mid _{L^{\infty }(\varOmega )}\le K_{1}, \forall t\ge T_0. \end{aligned}$$

for some large time \(T_0 > 0\). Hence, the result follows. \(\square\)

Note that R does not appear in the first three equations, so it is equivalent to study the following sub-system:

$$\begin{aligned} \begin{aligned}&\frac{\partial S}{\partial t}= \varLambda - \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-dS+d_{1}\nabla ^{2}S,\\&\frac{\partial I_{1}}{\partial t}=\frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}} -\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\xi I_{1}+d_{2}\nabla ^{2}I_{1} ,\\&\frac{\partial I_{2}}{ \partial t}=\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-\left( d+\mu _{2}+\rho _{2} \right) I_{2}+\xi I_{1}+d_{3}\nabla ^{2}I_{2}.\\ \end{aligned} \end{aligned}$$
(8)

4 Model Analysis

In the absence of spatial movement, we find the basic reproduction number of the system (1) as:

$$\begin{aligned} R_{0}=max\left\{ R_{0}^1,R_{0}^2\right\} , \end{aligned}$$

where

$$\begin{aligned} R_{0}^1=\frac{\alpha _{1} \varLambda }{(d+\mu _{1}+\rho _{1}+\xi )d}\;\text {and}\; R_{0}^2=\frac{\alpha _{2} \varLambda }{(d+\mu _{2}+\rho _{2})d}. \end{aligned}$$

Here, \(R_{0}^1(R_{0}^2)\) is the number of secondary infection cases generated by infection of the first (second) strain, called the reproduction number of the first (second) strain.

Also, we define two invasion reproduction numbers. By definition, the invasion reproduction number of the strain-1 \(R_{1}^{2}\) is given by

$$\begin{aligned} R_{1}^{2}=\frac{R_{0}^1}{R_{0}^2}=\frac{\alpha _{1}(d+\mu _{2}+\rho _{2})}{\alpha _{2}(d+\mu _{1}+\rho _{1}+\xi )}, \end{aligned}$$

which gives the number of secondary infections that an individual infected by the strain-1 will generate in a population in which the strain-2 is at equilibrium.

Similar to \(R_{1}^{2}\), the invasion reproduction number of the strain-2 \(R_{2}^{1}\) is given by

$$\begin{aligned} R_{2}^{1}=\frac{R_{0}^2}{R_{0}^1}=\frac{\alpha _{2}(d+\mu _{1}+\rho _{1}+\xi )}{\alpha _{1}(d+\mu _{2}+\rho _{2})}. \end{aligned}$$

4.1 Existence of Equilibria

There may be at most four equilibria for the non-spatial system corresponding to the spatial system (8), and they are given below

(i) Disease-free equilibrium:

The disease-free equilibrium point is \(E_{0}=\left( S^0, I_{1,0}, I_{2,0}\right)\), where \(S^0=\frac{\varLambda }{d}, I_{1,0}=I_{2,0}=0\)

(ii) Strain-1 endemic equilibrium:

To derive the strain-1 endemic equilibrium point \(E_{1}=\left( S_1, I_{1,1}, I_{2,1}\right)\), we consider \(I_{2,0}=0\) and we note that \(I_{2,0}=0 \implies I_{1,0}=0\) and \(\varLambda -d S_1=0\). Thus, no strain-1 endemic equilibrium can exist.

(iii) Strain-2 endemic equilibrium:

The strain-2 endemic equilibrium is \(E_{2}=\left( S_2, I_{1,2}, I_{2,2}\right)\), where

\(I_{1,2}=0\), \(S_{2}=\frac{\varLambda }{d+\frac{\alpha _{2}I_{2,2}}{1+\gamma _{2}I_{2,2}^{2}}}=\frac{ \left( 1+\gamma _{2}I_{2,2}^{2}\right) \varLambda }{d \left( 1+\gamma _{2}I_{2,2}^{2}\right) +\alpha _{2}I_{2,2}}\) and \(I_{2,2}\) is the positive root of the following equation

$$\begin{aligned} d \gamma _{2}I_{2,2}^{2}+ \alpha _{2} I_{2,2}-\left( \dfrac{\varLambda \alpha _{2}}{C}-d\right) =0 \end{aligned}$$
(9)

and it is given by \(I_{2,2}=\dfrac{-\alpha _{2}+\sqrt{\alpha _{2}^2+4 d \gamma _{2} \left( \dfrac{\varLambda \alpha _{2}}{C}-d\right) }}{2d \gamma _{2}}\).

Again, for the positivity of \(I_{2,2}\), we have \(\varLambda \alpha _{2}-Cd>0 \; \text {i.e.}\; \varLambda \alpha _{2} -d(d+\mu _{2}+\rho _{2})>0.\) i.e. \(R_0^2>1\).

(iv) Strain-1 and strain-2 endemic equilibrium:

Let us choose \(A=d+\mu _{1}+\rho _{1}+\xi _{1}\), \(B=d+\mu _{1}+\rho _{1}\), \(C=d+\mu _{2}+\rho _{2}.\)

Then, the endemic equilibrium point \(E_3=(S_3, I_{1,3}, I_{2,3})\) is given by

$$\begin{aligned} I_{1,3}= & {} \dfrac{\left( \frac{\alpha _{2}S_3}{1+\gamma _{2}I_{2,3}^{2}} -C\right) I_{2,3}}{-\xi _{1}}\\ S_3= & {} \dfrac{\varLambda \xi _{1}-B C I_{2,3}-C \xi _{1} I_{2,3}}{d \xi _{1}-B \frac{\alpha _{2} I_{2,3}}{1+\gamma _{2}I_{2,3}^{2}}}\\= & {} \dfrac{\left( \varLambda \xi _{1} - A C I_{2,3}\right) (1+\gamma _{2}I_{2,3}^{2}) }{\xi _1 ( d + \alpha _{2} I_{2,3} + d\gamma _{2}I_{2,3}^{2})-A \alpha _{2} I_{2,3}} \end{aligned}$$

and \(I_{2,3}\) is obtained from the following equation:

$$\begin{aligned}&a_{0} I_{2,3}^6+a_{1} I_{2,3}^5+a_{2} I_{2,3}^4\nonumber \\&\quad +a_{3} I_{2,3}^3+a_{4} I_{2,3}^2+a_{5} I_{2,3}+a_{6}=0 \end{aligned}$$
(10)

where

$$\begin{aligned} a_{0}= & {} A C^2 d^2 \gamma _{1} \gamma _{2}^2>0,\\ a_{1}= & {} A C d \gamma _2 (2 C \alpha _2 \gamma _1 + \alpha _1 \gamma _2 \xi _{1}) >0,\\ a_{2}= & {} -A^2 C \alpha _1 \alpha _2 \gamma _2 - d \alpha _1 \gamma _2^2 \varLambda \xi _{1}^2 + A (C^2 \gamma _1 (\alpha _2^2 + 2 d^2 \gamma _2) \\&+ d^2 \gamma _2^2 \xi _{1}^2 + C \alpha _2 \gamma _2 (-2 d \gamma _1 \varLambda + \alpha _1 \xi _{1})),\\ a_{3}= & {} -2 A^2 d \alpha _2 \gamma _2 \xi _{1} - \alpha _1 \alpha _2 \gamma _2 \varLambda \xi _{1}^2 \\&+ A (2 C^2 d \alpha _2 \gamma _1 + \alpha _2 \gamma _2 \xi _{1} ( \alpha _1 \varLambda + 2 d \xi _{1}) \\&+ C (-2 \alpha _2^2 \gamma _1 \varLambda + 2 d \alpha _1 \gamma _2 \xi _{1})),\\ a_{4}= & {} A^3 \alpha _2^2 - 2 d \alpha _1 \gamma _2 \varLambda \xi _{1}^2 - A^2 \alpha _2 (C \alpha _1 + 2 \alpha _2 \xi _{1})\\&+ A (C^2 d^2 \gamma _1 + 2 d^2 \gamma _2 \xi _{1}^2 + C (-2 d \alpha _2 \gamma _1 \varLambda + \alpha _1 \alpha _2 \xi _{1})\\&+ \alpha _2^2 (\gamma _1 \varLambda ^2 + \xi _{1}^2)),\\ a_{5}= & {} \xi _{1} (-2 A^2 d \alpha _2 - \alpha _1 \alpha _2 \varLambda \xi _{1}\\&+ A (C d \alpha _1 + \alpha _1 \alpha _2 \varLambda + 2 d \alpha _2 \xi _{1})),\\ a_{6}= & {} d (A d - \alpha _1 \varLambda ) \xi _{1}^2=d^2\xi _{1}^2(d+\mu _{1}\\&+\rho _{1}+\xi _{1})(1-R_0^1). \end{aligned}$$

Now, we note that under the conditions \(R_0^1>1\), \(R_0^2>1\), \(B\alpha _{1} \gamma _{2}>C\alpha _{2}\gamma _{1}\), \(B \alpha _{2}>C \alpha _1\), \(Ad>B \alpha _{2}\) and \(AB \alpha _2 (B \alpha _2-C \alpha _1)+ 2 d \gamma _2 \xi _{1}^2 (Ad-\varLambda \alpha _1 ) + A \gamma _1(\varLambda \alpha _2-Cd)^2<0\), we have \(a_0>0,a_1>0, a_2<0, a_3<0, a_4<0, a_5<0, a_6<0\). Hence, using Descartes' rule of sign we conclude that Eq. (10) has a unique root and we write the following proposition.

Proposition 3

If the conditions \(R_0^1>1\), \(R_0^2>1\), \(B\alpha _{1} \gamma _{2}>C\alpha _{2}\gamma _{1}\), \(B \alpha _{2}>C \alpha _1\), \(Ad>B \alpha _{2}\) and \(AB \alpha _2 (B \alpha _2-C \alpha _1)+ 2 d \gamma _2 \xi _{1}^2 (Ad-\varLambda \alpha _1 ) + A \gamma _1(\varLambda \alpha _2-Cd)^2<0\), are satisfied, then there is a unique positive solution of Eq. (10).

So, for this unique \(I_{2,3}\), we have a \(I_{1,3}\) and hence a \(S_{3}\). Thus, the endemic \(E_3=(S_3, I_{1,3}, I_{2,3})\) is unique under Proposition 3.

4.2 Global Stability Analysis

In this section, we investigate the global asymptotic stability of the model system (8) around all the possible equilibrium points. Here, we follow the technique as described by Hattaf and Yousfi (2013) to establish the global stability of the diffusive system. The Lyapunov functional are constructed based upon the Proposition 2.1. of Hattaf and Yousfi (2013). In fact, first we formulate the Lyapunov functional for the non-diffusive system, i.e., for the following ordinary differential equations system

$$\begin{aligned} \begin{aligned}&\frac{dS}{ dt}= \varLambda - \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-dS,\\&\frac{d I_{1}}{ dt}=\frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}} -\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\xi I_{1},\\&\frac{d I_{2}}{ dt}=\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-\left( d+\mu _{2}+\rho _{2} \right) I_{2}+\xi I_{1}.\\ \end{aligned} \end{aligned}$$
(11)

Then, we construct the Lyapunov functional for the diffusive system (8).

Theorem 1

The infection-free equilibrium point of the system (8) is globally asymptotically stable if \(R_0^2<1\) and \(\varLambda \alpha _{1} <d(d+\mu _{1}+\rho _{1})\).

Proof

Following Hattaf and Yousfi (2013), we first construct the Lyapunov functional for system (11) as:

$$\begin{aligned} L_1=S-S^0-S^0ln\left( \frac{S}{S^0} \right) +I_1+I_2. \end{aligned}$$
(12)

Now calculating the time derivative of Eq. (12) along the solutions of the system (11), we get

$$\begin{aligned} \frac{dL_1}{dt}= & \left( 1- \frac{S^0}{S}\right) \frac{dS}{dt} +\frac{dI_1}{dt} +\frac{dI_2}{dt}\nonumber \\= & \left( 1- \frac{S^0}{S}\right) \left( \varLambda - \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-dS \right) \nonumber \\&+ \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}} -\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\xi I_{1}+\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}\nonumber \\&-\left( d+\mu _{2}+\rho _{2} \right) I_{2}+\xi I_{1} \nonumber \\= & dS^0\left( 2- \frac{S}{S^0}-\frac{S^0}{S}\right) + \frac{S^0}{S}\left( \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}+ \frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}} \right) \nonumber \\&-\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\left( d+\mu _{2}+\rho _{2} \right) I_{2}\nonumber \\&\le dS^0\left( 2- \frac{S}{S^0}-\frac{S^0}{S}\right) +\frac{S^0}{S}\left( \alpha _{1}SI_{1}+\alpha _{2}SI_{2} \right) \nonumber \\&-\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\left( d+\mu _{2}+\rho _{2} \right) I_{2}\nonumber \\= & dS^0\left( 2- \frac{S}{S^0}-\frac{S^0}{S}\right) +\left( \frac{\varLambda \alpha _{1}}{d}-(d+\mu _{1}+\rho _{1}) \right) I_1\nonumber \\&+\left( \frac{\varLambda \alpha _{2}}{d}-(d+\mu _{2}+\rho _{2}) \right) I_2\nonumber \\= & dS^0\left( 2- \frac{S}{S^0}-\frac{S^0}{S}\right) + \left( \frac{\varLambda \alpha _{1}}{d}-(d+\mu _{1}+\rho _{1}) \right) I_1\nonumber \\&+(R_0^2-1)(d+\mu _{2}+\rho _{2})I_2 \end{aligned}$$
(13)

Now, by applying A.M.\(\ge\)G.M., we have \(\left( 2- \frac{S}{S^0}-\frac{S^0}{S}\right) \le 0\). In addition, under the conditions \(R_0^2<1\) and \(\varLambda \alpha _{1} <d(d+\mu _{1}+\rho _{1})\), we observe that \(\frac{dL_1}{dt}\le 0\).

Now we construct the Lyapunov functional for the diffusive system (8), at the disease-free state \(E_0\) as follows:

$$\begin{aligned} W_1=\int _{\varOmega } L_1(S(x,t),I(x,t),I_2(x,t))dx \end{aligned}$$

Calculating the time derivative of \(W_1\) along the solutions of the system (8), we get

$$\begin{aligned} \frac{dW_1}{dt}\le & {} \int _{\varOmega }\left\{ dS^0\left( 2- \frac{S}{S^0}-\frac{S^0}{S}\right) +\left( \frac{\varLambda \alpha _{1}}{d}-(d+\mu _{1}+\rho _{1}) \right) I_1\right. \nonumber \\&\left. +(R_0^2-1)(d+\mu _{2}+\rho _{2})I_2\right\} \nonumber \\&-d_1S^0\int _{\varOmega }\frac{|\nabla S|^2}{S^2} dx \end{aligned}$$
(14)

Now, under the conditions \(R_0^2<1\) and \(\varLambda \alpha _{1} <d(d+\mu _{1}+\rho _{1})\), we have \(\frac{dW_1}{dt}\le 0\). Also, the largest set in \(\Theta _1=\left\{ (S,I_1,I_2):\frac{dW_1}{dt}=0 \right\}\) is the singleton \(E_0\). Hence, from the LaSalle invariance principle (La Salle 1976), we conclude that the infection-free equilibrium point \(E_0\) is globally asymptotically stable under the stated conditions. Hence the theorem is proved. \(\square\)

Theorem 2

The strain-2 endemic equilibrium point \(E_2(S_2,0,I_{2,2})\) is globally asymptotically stable if \(R_0^2>1\), \(S_2<min\left\{ \frac{d+\mu _{1}+\rho _{1}}{\alpha _{1}},\frac{d+\mu _{2}+\rho _{2}}{\alpha _{2}} \right\}\) and \(\xi I_1+ \frac{\alpha _{1} S I_2}{1+\gamma _{2} I_2^2}>2(d+\mu _{2}+\rho _{2})I_2\).

Proof

As earlier, following Hattaf and Yousfi (2013), we first construct the Lyapunov functional for system (11) as

$$\begin{aligned} L_2=S-S_2-S_2ln\left( \frac{S}{S_2} \right) +I_1+I_2-I_{2,2}-I_{2,2}ln\left( \frac{I_2}{I_{2,2}} \right) . \end{aligned}$$
(15)

Now calculating the time derivative of the equation (15) along the solutions of the system (11), we get

$$\begin{aligned} \frac{dL_2}{dt}= \left( 1- \frac{S_2}{S}\right) \frac{dS}{dt} +\frac{dI_1}{dt} +\left( 1-\frac{I_{2,2}}{I_2} \right) \frac{dI_2}{dt}\nonumber \\= \left( 1- \frac{S_2}{S}\right) \left( \varLambda - \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}}-dS \right) \nonumber \\+ \frac{\alpha _{1}SI_{1}}{1+\gamma _{1}I_{1}^{2}}-\left( d+\mu _{1}+\rho _{1} \right) I_{1}-\xi I_{1}\nonumber \\+ \left( 1-\frac{I_{2,2}}{I_2} \right) \left( \frac{\alpha _{2}SI_{2}}{1+\gamma _{2}I_{2}^{2}} -\left( d+\mu _{2}+\rho _{2} \right) I_{2}+\xi I_{1}\right) \nonumber \\= dS_2\left( 2- \frac{S}{S_2}-\frac{S_2}{S}\right) +\frac{\alpha _{2} S_2I_{2,2}}{1+\gamma _{2} I_{2,2}^2}\nonumber \\-\frac{S_2}{S}\left( \frac{\alpha _{2} S_2I_{2,2}}{1+\gamma _{2} I_{2,2}^2}-\frac{\alpha _{1} SI_{1}}{1+\gamma _{1} I_{1}^2}-\frac{\alpha _{2} SI_{2}}{1+\gamma _{2} I_{2}^2} \right) \nonumber \\-(d+\mu _{1}+\rho _{1})I_1\nonumber \\-(d+\mu _{2}+\rho _{2})I_2-\frac{\alpha _{2} SI_{2,2}}{1+\gamma _{2} I_{2}^2}\nonumber \\+(d+\mu _{2}+\rho _{2})I_{2,2}-\frac{\xi I_1I_{2,2}}{I_2}\nonumber \\\le dS_2\left( 2- \frac{S}{S_2}-\frac{S_2}{S}\right) +\frac{\alpha _{2} S_2I_{2,2}}{1+\gamma _{2} I_{2,2}^2}\nonumber \\+\frac{\alpha _{1} S_2I_{1}}{1+\gamma _{1} I_{1}^2}+\frac{\alpha _{2} S_2I_{2}}{1+\gamma _{2} I_{2}^2}\nonumber \\-(d+\mu _{1}+\rho _{1})I_1-(d+\mu _{2}+\rho _{2})I_2\nonumber \\-\frac{\alpha _{2} SI_{2,2}}{1+\gamma _{2} I_{2}^2}+(d+\mu _{2}+\rho _{2})I_{2,2}-\frac{\xi I_1I_{2,2}}{I_2}\nonumber \\\le dS_2\left( 2- \frac{S}{S_2}-\frac{S_2}{S}\right) +\frac{\alpha _{2} S_2I_{2,2}}{1+\gamma _{2} I_{2,2}^2}\nonumber \\+\alpha _{1} S_2I_{1}+\alpha _{2} S_2I_{2}-(d+\mu _{1}+\rho _{1})I_1-(d+\mu _{2}+\rho _{2})I_2\nonumber \\-\frac{\alpha _{2} SI_{2,2}}{1+\gamma _{2} I_2^2}+(d+\mu _{2}+\rho _{2})I_{2,2}-\frac{\xi I_1I_{2,2}}{I_2}\nonumber \\= {} dS_2\left( 2- \frac{S}{S_2}-\frac{S_2}{S}\right) +\left( \alpha _{1} S_2-(d+\mu _{1}+\rho _{1})\right) I_1\nonumber \\+\left( \alpha _{2} S_2-(d+\mu _{2}+\rho _{2})\right) I_2\nonumber \\-\frac{I_{2,2}}{I_2}\left( \frac{\alpha _{2} S I_2}{1+\gamma _{2} I_2^2}+\xi I_1 -2(d+\mu _{2}+\rho _{2})I_2\right) \end{aligned}$$
(16)

Now, by applying A.M.\(\ge\)G.M., we have \(\left( 2- \frac{S}{S_2}-\frac{S_2}{S}\right) \le 0\). Further, under the conditions \(S_2<min\left\{ \frac{d+\mu _{1}+\rho _{1}}{\alpha _{1}}, \frac{d+\mu _{2}+\rho _{2}}{\alpha _{2}} \right\}\) and \(\xi I_1+ \frac{\alpha _{1} S I_2}{1+\gamma _{2} I_2^2}\,>\,2(d+\mu _{2}+\rho _{2})I_2\), we note that \(\frac{dL_2}{dt}\le 0\).

Now we construct the Lyapunov functional for the diffusive system (8), at the strain-2 endemic steady state \(E_2\) as follows:

$$\begin{aligned} W_2=\int _{\varOmega } L_2(S(x,t),I_1(x,t),I_2(x,t))dx \end{aligned}$$

Calculating the time derivative of \(W_2\) along the solutions of the system (8), we get

$$\begin{aligned} \frac{dW_2}{dt}\le \int _{\varOmega } dS_2\left( 2- \frac{S}{S_2} -\frac{S_2}{S}\right) +\left( \alpha _{1} S_2-(d+\mu _{1}+\rho _{1})\right) I_1\\+\left( \alpha _{2} S_2-(d+\mu _{2}+\rho _{2})\right) I_2\\-\frac{I_{2,2}}{I_2}\left( \frac{\alpha _{2} S I_2}{1+\gamma _{2} I_2^2}+\xi I_1 -2(d+\mu _{2}+\rho _{2})I_2\right) dx\\-d_1S_2\int _{\varOmega }\frac{|\nabla S|^2}{S^2} dx\\-d_2I_{2,2}\int _{\varOmega }\frac{|\nabla I_2|^2}{I_2^2} dx \end{aligned}$$

Now, under the conditions as stated earlier, we have \(\frac{dW_2}{dt}\le 0\). Also, the largest set in \(\Theta _2=\left\{ (S,I_1,I_2):\frac{dW_2}{dt}=0 \right\}\) is the singleton \(E_2\). Hence, from the LaSalle invariance principle (La Salle 1976), we conclude that the strain-2 endemic equilibrium point \(E_2\) is globally asymptotically stable under the stated conditions. Hence, the theorem follows. \(\square\)

Theorem 3

The strain-1 and strain-2 endemic equilibrium point \(E_3(S_3,I_{1,3},I_{2,3})\) is globally asymptotically stable under the conditions \(R_0^1>1\), \(R_0^2>1\) and either \(\frac{I_2}{I_{2,3}}\le 1\le \frac{I_1}{I_{1,3}}\), \(\frac{I_2}{I_{2,3}}\le \frac{S}{S_3}\le \frac{I_1}{I_{1,3}}\) or \(\frac{I_1}{I_{1,3}}\le 1\le \frac{I_2}{I_{2,3}}\), \(\frac{I_1}{I_{1,3}} \le \frac{S}{S_3}\le \frac{I_2}{I_{2,3}}\).

Proof

As earlier, following Hattaf and Yousfi (2013), we first construct the Lyapunov functional for system (11) as:

$$\begin{aligned} L_3= & S-S_3-S_3ln\left( \frac{S}{S_3} \right) +I_1-I_{1,3}-I_{1,3}ln \left( \frac{I_1}{I_{1,3}}\right) \nonumber \\&+I_2-I_{2,3}-I_{2,3}ln\left( \frac{I_2}{I_{2,3}} \right) . \end{aligned}$$
(17)

Now calculating the time derivative of the equation (17) along the solutions of the system (11), we get

$$\begin{aligned} \frac{dL_3}{dt}= & \left( 1- \frac{S_3}{S}\right) \frac{dS}{dt} +\left( 1-\frac{I_{1,3}}{I_1} \right) \frac{dI_1}{dt}\nonumber \\&+\left( 1-\frac{I_{2,3}}{I_2} \right) \frac{dI_2}{dt}\nonumber \\= & \left( S-S_3\right) \frac{\dot{S}}{S}+\left( I_1-I_{1,3}\right) \frac{\dot{I_1}}{I_1}+\left( I_2-I_{2,3}\right) \frac{\dot{I_2}}{I_2} \end{aligned}$$
(18)

Now for clarity, we calculate the terms \(\frac{\dot{S}}{S}\), \(\frac{\dot{I_1}}{I_1}\) and \(\frac{\dot{I_2}}{I_2}\) separately.

$$\begin{aligned} \frac{\dot{S}}{S}= & \frac{\varLambda }{S}-\frac{\alpha _{1}I_{1}}{1+\gamma _{1}I_{1}^{2}}\\&-\frac{\alpha _{2}I_{2}}{1+\gamma _{2}I_{2}^{2}}-d\\= & {} \left( \frac{\varLambda }{S}-\frac{\alpha _{1}I_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{2}I_{2}}{1+\gamma _{2}I_{2}^{2}}-d\right) \\&-\left( \frac{\varLambda }{S_3}- \frac{\alpha _{1}I_{1,3}}{1+\gamma _{1}I_{1,3}^{2}}-\frac{\alpha _{2}I_{2,3}}{1+\gamma _{2}I_{2,3}^{2}}-d \right) \\= & \frac{\varLambda \left( S_3-S \right) }{SS_3}-\alpha _{1} \left( \frac{I_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{I_{1,3}}{1+\gamma _{1}I_{1,3}^{2}} \right) \\&-\alpha _{2}\left( \frac{I_{2}}{1+\gamma _{2}I_{2}^{2}}-\frac{I_{2,3}}{1+\gamma _{2}I_{2,3}^{2}} \right) \\= & \frac{\varLambda \left( S_3-S \right) }{SS_3}-\alpha _{1} \left( \frac{I_{1}}{1+\gamma _{1}I_{1}^{2}}-\frac{I_{1,3}}{1+\gamma _{1}I_{1}^{2}}\right. \\&\left. +\frac{I_{1,3}}{1+\gamma _{1}I_{1}^{2}}-\frac{I_{1,3}}{1+\gamma _{1}I_{1,3}^{2}} \right) \\&-\alpha _{2}\left( \frac{I_{2}}{1+\gamma _{2}I_{2}^{2}}-\frac{I_{2,3}}{1+\gamma _{2}I_{2}^{2}}+\frac{I_{2,3}}{1+\gamma _{2}I_{2}^{2}}-\frac{I_{2,3}}{1+\gamma _{2}I_{2,3}^{2}}\right) \\&=\frac{\varLambda \left( S_3-S \right) }{SS_3}-\alpha _{1}\left( \frac{I_1-I_{1,3}}{(1+\gamma _{1}I_{1}^{2})}-\frac{\gamma _{1} I_{1,3}(I_1^2-I_{1,3}^2)}{(1+\gamma _{1}I_{1}^{2})(1+\gamma _{1}I_{1,3}^{2})} \right) \\&-\alpha _{2} \left( \frac{I_2-I_{2,3}}{(1+\gamma _{2}I_{2}^{2})}-\frac{\gamma _{2} I_{2,3}(I_2^2-I_{2,3}^2)}{(1+\gamma _{2}I_{2}^{2})(1+\gamma _{2}I_{2,3}^{2})} \right) \end{aligned}$$
$$\begin{aligned} \begin{aligned} \frac{\dot{I_1}}{I_1}&=\frac{\alpha _{1}S}{1+\gamma _{1}I_{1}^{2}}-\left( d+\mu _{1}+\rho _{1} +\xi \right) \\&=\frac{\alpha _{1}S}{1+\gamma _{1}I_{1}^{2}}-\frac{\alpha _{1}S_3}{1+\gamma _{1}I_{1,3}^{2}}\\&=\alpha _{1} \left( \frac{S}{1+\gamma _{1}I_{1}^{2}}-\frac{S_3}{1+\gamma _{1}I_{1}^{2}}\right. \\&\left. \quad+\frac{S_3}{1+\gamma _{1}I_{1}^{2}}-\frac{S_3}{1+\gamma _{1}I_{1,3}^{2}}\right) \\&= \alpha _{1} \left( \frac{S-S_3}{(1+\gamma _{1}I_{1}^{2})}\right. \\&\left.\quad -\frac{\gamma _{1} S_3 (I_1^2-I_{1,3}^2)}{(1+\gamma _{1}I_{1}^{2})(1+\gamma _{1}I_{1,3}^{2})} \right) \end{aligned} \end{aligned}$$
$$\begin{aligned} \begin{aligned} \frac{\dot{I_2}}{I_2}&=\frac{\alpha _{2}S}{1+\gamma _{2}I_{2}^{2}}-\left( d+\mu _{2}+\rho _{2} \right) +\frac{\xi I_{1}}{I_{2}}\\&= \frac{\alpha _{2}S}{1+\gamma _{2}I_{2}^{2}}-\left( d+\mu _{2}+\rho _{2} \right) \\&+\frac{\xi I_{1}}{I_{2}}-\left( \frac{\alpha _{2}S_3}{1+\gamma _{2}I_{2,3}^{2}}-\left( d+\mu _{2}+\rho _{2} \right) +\frac{\xi I_{1,3}}{I_{2,3}} \right) \\&= \alpha _2 \left( \frac{S}{1+\gamma _{2}I_{2}^{2}}-\frac{S_3}{1+\gamma _{2}I_{2,3}^{2}} \right) \\&\quad+\xi \left( \frac{ I_{1}}{I_{2}}-\frac{ I_{1,3}}{I_{2,3}}\right) \\&= \alpha _2 \left( \frac{S}{1+\gamma _{2}I_{2}^{2}}-\frac{S_3}{1+\gamma _{2}I_{2}^{2}}+\frac{S_3}{1+\gamma _{2}I_{2}^{2}}-\frac{S_3}{1+\gamma _{2}I_{2,3}^{2}} \right) \\&\quad+\xi \left( \frac{ I_{1}}{I_{2}}-\frac{ I_{1,3}}{I_{2,3}}-\frac{ I_{1,3}}{I_{2}}+\frac{ I_{1,3}}{I_{2}}\right) \\&=\alpha _{2} \left( \frac{S-S_3}{(1+\gamma _{2}I_{2}^{2}) }+\frac{\gamma _{2} S_3(I_{2,3}^2-I_2^2)}{(1+\gamma _{2}I_{2}^{2})(1+\gamma _{2}I_{2,3}^{2})}\right) \\&\quad+\xi \left( \frac{I_1-I_{1,3}}{I_2}-\frac{I_{1,3}(I_2-I_{2,3})}{I_2I_{2,3}} \right) \end{aligned} \end{aligned}$$

Now substituting the values of \(\frac{\dot{S}}{S}\), \(\frac{\dot{I_1}}{I_1}\) and \(\frac{\dot{I_2}}{I_2}\) in the equation (18), we get

$$\begin{aligned}&\frac{dL_3}{dt}=(S-S_3) \left( \frac{\varLambda \left( S_3-S \right) }{SS_3}-\alpha _{1}\left( \frac{I_1-I_{1,3}}{(1+\gamma _{1}I_{1}^{2})}\right. \right. \\&\left. \left. -\frac{\gamma _{1}I_{1,3}(I_1^2-I_{1,3}^2)}{(1+\gamma _{1}I_{1}^{2})(1+\gamma _{1}I_{1,3}^{2})} \right) \right. \\&\left. -\alpha _{2} \left( \frac{I_2-I_{2,3}}{(1+\gamma _{2}I_{2}^{2})}-\frac{\gamma _{2} I_{2,3}(I_2^2-I_{2,3}^2)}{(1+\gamma _{2}I_{2}^{2})(1+\gamma _{2}I_{2,3}^{2})} \right) \right) \\&+ (I_1-I_{1,3})\left( \alpha _{1} \left( \frac{S-S_3}{(1+\gamma _{1}I_{1}^{2})}-\frac{\gamma _{1} S_3 (I_1^2-I_{1,3}^2)}{(1+\gamma _{1}I_{1}^{2})(1+\gamma _{1}I_{1,3}^{2})} \right) \right) \\&+(I_2-I_{2,3}) \left( \alpha _{2} \left( \frac{S-S_3}{(1+\gamma _{2}I_{2}^{2}) }+\frac{\gamma _{2} S_3(I_{2,3}^2-I_2^2)}{(1+\gamma _{2}I_{2}^{2})(1+\gamma _{2}I_{2,3}^{2})}\right) +\right. \\&\left. \xi \left( \frac{I_1-I_{1,3}}{I_2}-\frac{I_{1,3}(I_2-I_{2,3})}{I_2I_{2,3}} \right) \right) \end{aligned}$$

After some simplifications, we finally obtain

$$\begin{aligned} \begin{aligned} \frac{dL_3}{dt}&=-\frac{\varLambda (S-S_3)^2}{SS_3}-\frac{\xi I_{1,3} (I_2-I_{2,3})^2}{I_2I_{2,3}}\\&\quad+\frac{\xi (I_1-I_{1,3})(I_2-I_{2,3}) }{I_2}\\&\quad+\frac{\alpha _{1} \gamma _{1} (I_1-I_{1,3})(I_1+I_{1,3})(SI_{1,3}-S_3I_1)}{(1+\gamma _{1}I_{1}^{2}) (1+\gamma _{1}I_{1,3}^{2})}\\&\quad+\frac{\alpha _{2} \gamma _{2} (I_2-I_{2,3})(I_2+I_{2,3})(SI_{2,3}-S_3I_2)}{(1+\gamma _{2}I_{2}^{2}) (1+\gamma _{2}I_{2,3}^{2})} \end{aligned} \end{aligned}$$
(19)

Now we note that if either \(\frac{I_2}{I_{2,3}}\le 1\le \frac{I_1}{I_{1,3}}\), \(\frac{I_2}{I_{2,3}}\le \frac{S}{S_3}\le \frac{I_1}{I_{1,3}}\) or \(\frac{I_1}{I_{1,3}}\le 1\le \frac{I_2}{I_{2,3}}\), \(\frac{I_1}{I_{1,3}} \le \frac{S}{S_3}\le \frac{I_2}{I_{2,3}}\), then \(\frac{dL_3}{dt}\le 0\).

Now we construct the Lyapunov functional for the diffusive system (8), at the strain-1 and strain-2 endemic steady state \(E_3\) as follows:

$$\begin{aligned} W_3=\int _{\varOmega } L_3(S(x,t),I_1(x,t),I_2(x,t))dx \end{aligned}$$

Calculating the time derivative of \(W_3\) along the solutions of the system (8), we get

$$\begin{aligned} \begin{aligned} \frac{dW_3}{dt}&=\int _{\varOmega } -\frac{\varLambda (S-S_3)^2}{SS_3}-\frac{\xi I_{1,3} (I_2-I_{2,3})^2}{I_2I_{2,3}}\\&\quad+\frac{\xi (I_1-I_{1,3})(I_2-I_{2,3}) }{I_2}\\&\quad+\frac{\alpha _{1} \gamma _{1} (I_1-I_{1,3})(I_1+I_{1,3})(SI_{1,3}-S_3I_1)}{(1+\gamma _{1}I_{1}^{2})(1+\gamma _{1}I_{1,3}^{2})}\\&\quad+\frac{\alpha _{2} \gamma _{2} (I_2-I_{2,3})(I_2+I_{2,3})(SI_{2,3}-S_3I_2)}{(1+\gamma _{2}I_{2}^{2})(1+\gamma _{2}I_{2,3}^{2})} dx\\&\quad-d_1S_3\int _{\varOmega }\frac{|\nabla S|^2}{S^2} dx-d_2I_{1,3}\int _{\varOmega }\frac{|\nabla I_2|^2}{I_2^2} dx\\&\quad-d_3I_{2,3}\int _{\varOmega }\frac{|\nabla I_3|^2}{I_3^2} dx \end{aligned} \end{aligned}$$

Therefore, under the conditions as stated earlier, we get \(\frac{dW_3}{dt}\le 0\). Also, the largest set in \(\Theta _3=\left\{ (S,I_1,I_2):\frac{dW_3}{dt}=0 \right\}\) is the singleton \(E_3\). Hence, from the LaSalle invariance principle (La Salle 1976), it follows that the strain-1 and strain-2 endemic equilibrium point \(E_3\) is globally asymptotically stable under the stated conditions. Hence, the theorem follows. \(\square\)

5 Conclusions

In this paper, we have proposed and analyzed a diffusive epidemic model with two strains and non-monotonic incidence rates for both the strains. These types of incidence rate are considered to interpret the psychological effects. The recent pandemic outbreak due to SARS-CoV-2 has such psychological effects on the common people; also several strains of the virus are found. So, this work can be applied to study transmission mechanisms of the corona virus disease. The biological feasibility conditions of the solutions of the system are studied. All the possible equilibrium points of the system are obtained and their global asymptotic stability are investigated. The parametric conditions are derived for the long-term occurrence of the disease-free steady state and the strain-1 and strain-2 endemic equilibrium point. It is observed that the model system has no strain-1 endemic steady state. This scenario occurs possibly due to the presence of the virus mutation. It indicates that the disease due to strain-1 do not stay permanently in the community.

An epidemic model becomes more useful from a practical point of view when it is calibrated with the real world data. Also it is quite important to visualize the obtained theoretical results in graphical way. We would like to address these issues in our future research which requires more advanced computational skills. However, this theoretical study may be very useful to explain the dynamical characteristics of a diffusive two strain epidemic model with virus mutation and non-monotonic incidence rate.