A mathematical model for the impact of disinfectants on the control of bacterial diseases

Abstract

Here, we investigate a mathematical model to assess the impact of disinfectants in controlling diseases that spread in the population via direct contacts with the infected persons and also due to bacteria present in the environment. We find that the disease-free and endemic equilibria of the system are related via a transcritical bifurcation whose direction is forward. Our numerical results show that controlling the transmissions of disease through direct contacts and bacteria present in the environment can help in reducing the disease prevalence. Moreover, fostering the recovery rate and the death rate of bacteria play significant roles in disease eradication. Our numerical observations convey that reducing the bacterial density at the source discharged by the infected population through the use of chemicals has prominent effect in disease control. Overall, our findings manifest that the disinfectants of high quality can completely control the bacterial density and the disease outbreak.

Keywords:

• Epidemic model
• awareness
• disinfectants
• stability
• sensitivity

1. Introduction

The utmost challenging issue for the survival of human lives apart from war, starvation and natural disasters is the spread of infectious diseases around the world. The outbreak of these diseases not only enhances the mortality rate but also requires an enormous amount of wealth to prevent them [14]. Therefore, it is essential to find the means to control the prevalence of diseases. Infectious diseases can be transmitted through the air (air-borne transmission), water (water-borne transmission), and the direct contact with the infected persons and also indirectly via the intermediate objects known as vectors/carriers. Bacterial diseases such as cholera and typhoid still remain as a major public health problem in developing countries, where outbreaks continue to occur due to poor sanitation and personal hygiene [8, 15]. Pathogenes of bacterial diseases, like cholera (V. cholerae), typhoid (Salmonella typhi), diarrhoea and dysentery (Shigella) spread through contaminated food and water or from an infected person to a healthy one as a result of poor hygiene. Some individuals produce growth substrates surrounding the bacterial cells, thereby provide a very conducive environment to the bacteria, resulting in the fast spread of the diseases. The clinically infected animals/humans harbour and shed bacteria in the environment, and act as a source of infection to susceptible humans through direct contact or indirectly through fecal-contaminated environments [13]. Globally, more than 600 million people living in urban areas are still facing a lack of basic sanitation, and the total urban population is expanding by around 80 million per year [49], making it a challenge to keep up with sanitation infrastructure requirements. Consumption of polluted drinking water is estimated to cause more than five hundred thousand diarrhoeal deaths every year [48]. According to WHO, every year there are 1.3–4.0 million cases of cholera and 21,000–143,000 deaths worldwide [3]. Therefore, it has become desirable to wipe out pathogenic bacteria from the environment by some means to control the prevalence of diseases. Pathogenic bacteria, such as E. coli pathogenic strains, Campylobacter, V. cholera , Salmonella, and Shigella are treated by using disinfectants (chlorine, chloramine, ozone, hydrogen peroxide with UV and chemical coagulation), thereby preventing waterborne diseases.

The study of various infectious diseases has been carried out in the past few decades by considering the direct contact of susceptible with infectives, without taking into account the role of bacteria present in the environment [4, 14, 16, 17]. However, there are certain diseases such as tuberculosis and typhoid fever which spread via direct contacts with infectives as well as by the bacteria present in the environment [9, 12, 34–36, 51]. Therefore, it is salient to study the dynamics of such diseases so that proper control strategies can be applied. Mathematical modelling is a useful methodology for evaluating the hypotheses that underlie disease processes, which can be applied to assess the impact of control efforts and predict future trends. Several kinds of mathematical models have been proposed to understand the bacterial transmission and their control. Singh et al. [41] have assumed that disease spreads in the population due to the direct contact between susceptible and infected individuals as well as due to the bacteria transported by the carriers, like flies and ticks from the environment to the water and edibles of human population. They have considered that the density of carrier population follows logistic growth and human-related activities make the environment conducive for the growth of carrier population. They found that the disease can be controlled by controlling the growth rate of density of carrier population; however the immigration of population to susceptible class and human-related activities make the disease endemic. This model has been generalized by Misra et al. [27, 28] for cholera and typhoid diseases by incorporating the fact that the density of carrier population in the environment can be controlled using the chemicals like DDT and chlorine. Authors have concluded that the use of chemicals not only control the carrier population but also reduces the number of infected individuals in the community.

Notably, in the aforementioned studies, the density of bacteria has not been considered explicitly, and the authors have assumed that the bacteria of the disease are present in the environment. Shukla et al. [39] have incorporated the density of bacteria explicitly as a dynamic variable in the modelling phenomenon. Here, authors have assumed that the susceptible individuals contract the infection due to direct contacts with infected individuals and indirectly either through the bacteria transported by the carriers to the edibles or due to the movement of susceptible individuals in the bacteria filled environment. The findings showed that the growth rate of bacteria in the environment due to discharge of excreta of infected individuals causes challenges to control the spread of the disease. A mathematical model for the control of cholera disease using chemicals that reduce the density of bacteria in the aquatic environment was analysed by Misra and Singh [32]. Authors have assumed here a bilinear interaction between susceptible individuals to the density of bacteria, and that the rate of use of chemicals in the aquatic environment is proportional to the density of bacteria in the aquatic environment with some time lag. Shukla et al. [40] have studied the effect of sanitation effort on the spread of bacterial diseases, and found that if sanitation effort is suitably applied to the human habitat, the density of the bacteria declines which in turn reduces the number of infected population. Musa et al. [33] incorporated public health education programmes (that include house-to-house enlightenment, media coverage, reports, etc.) to assess its impact and provide effective strategies to combat typhoid fever in the endemic areas. They obtained wavelet results to show the temporal patterns of typhoid fever outbreaks in Taiwan.

Media coverage is also a key effective approach for controlling the spread of infectious diseases because it raises public awareness of disease transmission modes, which leads to behavioural changes in people's attitudes about disease transmission. Some studies have been conducted to control the spread of infectious diseases with the impact of media [10, 11, 20, 21, 29–31, 37, 38, 44]. In particular, Dubey et al. [10] have studied the joint effects of awareness programmes through media and treatment of disease in infected population. They have concluded that the optimal treatment together with awareness can help in the total eradication of infectious disease from the population. Misra et al. [29] have developed a mathematical model to study the effect of TV and social media advertisements to control the infectious diseases. They observed that more TV and social media advertisements cause destabilization of the system and induce oscillations through Hopf-bifurcation. Recently, Lata et al. [19] have studied a mathematical model for the water-borne diseases, like cholera by incorporating TV and social media advertisements as a dynamic variable. They have assumed that aware individuals adapt some precautionary measures not to contract the infection and they also make some effectual efforts, like chlorination of water, etc., to reduce the density of bacteria in the aquatic environment. The outcomes of this study showed that awareness among the people due to TV and social media advertisements, and chlorination of water made by aware individuals have substantial impact on the control of water-borne diseases. Majumder et al. [23] investigated the impact of saturated treatments on HIV-TB dual epidemic as a consequence of COVID-19. They have explored the joint effects of awareness and treatment by analysing an optimal control problem. Massard et al. [26] investigated a multi-strain epidemic model for COVID-19 by considering the infected as well as asymptomatic cases. The study reveals that the beta/gamma variants could have led to a larger number of infections in both symptomatic and asymptomatic people. The findings of Sun et al. [42] suggest that stamping out the infected individuals and blocking the epidemic zones play significant roles in the prevention and control of the foot-and-mouth diseases. Ma et al. [22] observed that by enhancing the effective coverage of mask use in public, the cumulative confirmed cases can be reduced to a great extent. A 10 year systematic review of mathematical models for dengue fever epidemiology by Aguiar et al. [2] provide insights on general features to be considered to model the aspects of real-world public health problems, such as the current epidemiological scenario (COVID-19 pandemic) we are living in. The study of Asamoah et al. [5] demonstrated that in the absence of vaccination, practising physical or social distancing protocols is the most cost-saving and most effective control intervention in Saudi Arabia.

Besides the aforementioned studies, it is crucial to reduce the density of bacteria in the environment for the eradication of the disease from the community. To achieve this goal, i.e. eradication of disease, in this paper, we formulate and analyse a mathematical model to control the density of bacteria in the environment using chemicals. Here, we assume that a part of the chemical is used to treat the excreta of infected human population, i.e. reducing the density of bacteria at source and the remaining part of the chemical is used to kill the bacteria in the environment.

We arrange the remaining portions of this paper in the following way. In the next section, we propose our epidemic model. In the following section, the model is analysed; positivity and boundedness of the system's solutions are shown; equilibria are obtained and their stability properties are discussed; an explicit expression for the basic reproduction number is derived by using the next-generation matrix method; the existence of transcritical bifurcation between the disease-free and endemic equilibria is shown. The simulation results are carried out in Section 4 to extrapolate the analytical findings and explore the role of disinfectants in controlling the disease. The paper ends with important observations of this study in Section 5.

2. The mathematical model

Here, we propose a mathematical model for the diseases which spread among the human population due to direct contact of the susceptible persons with the infected ones and also through the bacteria of these diseases present in the environment. To this end, we classify the total human population of the considered region into three disjoint groups: susceptible individuals (S), infected individuals (I) and the recovered individuals (R). Let the density of bacteria present in the environment of the considered domain be B. Also, we assume that $C_h$ amount of disinfectants is being used to control the excessive bacterial growth at the hot spots of the infected individuals.

Table 1 summarizes all the considered dynamical variables for the model to be investigated here. In the formulation of the model, the following ecological/epidemiological assumptions are made.

Table 1. Descriptions of variables used in the model system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ).

Variables	Descriptions
S	Number of susceptible human population
I	Number of infected human population
R	Number of recovered human population
B	Density of bacteria present in the environment
$C_h$	Chemical concentration of disinfectants/sanitizer

1. All the new immigrants are susceptible and join the group at a constant rate Λ. The population is homogeneously mixed and disease spreads via the direct contact of susceptible with infected individuals, and also through bacteria present in the environment.

2. All the susceptible individuals living in the habitat are affected by the bacterial populations. Besides its self-growth, the density of bacteria in the environment increases as the number of infected individuals increases.

3. Disinfectants/sanitizers potentially suppress the density of bacteria at source (medical hospitals, infected hot spots) discharged by infected individuals, and are introduced at a rate proportional to density of bacteria present in the environment.

4. The bacterial populations at source discharged by infected individuals could be minimized to a limited extent by the use of disinfectants/sanitizer due to their limited effectiveness and also the cost involved for the purchase.

5. The concentration of chemicals reduces due to its uptake by bacteria present in the environment apart from its natural depletion.

6. The concentration of disinfectants/sanitizers reduces due to their use to control the density of bacteria at source discharged by infected individuals.

7. The infected individuals recover from infection through proper treatment and subsequently join the recovered class.

8. The individuals in the recovered class lose their immunity with passage of time and become susceptible again.

9. There is natural mortality in each human cohort, whereas infected individuals experience disease induced death additionally.

Considering the above assumptions all together, a schematic diagram representing the problem is depicted in Figure 1, and the model equations are obtained as follows: (1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) Let $S\left(0\right)>0$,  $I\left(0\right)\ge 0$,  $R\left(0\right)\ge 0$,  $B\left(0\right)\ge 0$ and $C_h\left(0\right)\ge 0$ be the initial conditions for solving the above set of differential equations.

Figure 1. Schematic diagram of model system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ).

In model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ), we have assumed that the density of bacteria increases in the environment due to the self-growth and discharge by the infected individuals (mainly excreta of infectives). The chemical disinfectants are used to control the bacteria in the environment. A part of the chemical disinfectants (i.e. $k{C}_h$) is assumed to be used for treating the excreta of infectives (i.e. the control of bacteria at source) and the remaining is used to reduce the bacteria in the environment. As the chemical disinfectants have the limited effect on the treatment of excreta (i.e. all the bacteria present in the excreta cannot be killed in the single application of chemical disinfectants), so we have considered its limited impact $\left(i.e.\frac{k{C}_h}{M+k{C}_h}\right)$. Further, it is worthy to note that the coefficient of decay in the density of bacteria present in the environment due to uptake of chemicals (i.e. ${\pi }_1$) should be greater than one as a small amount of chemical disinfectants may kill a large amount of bacteria. Further, the per-capita natural death rate of bacteria will be greater than its self growth rate, i.e. $s_0-s=s_3>0$. Apparently, if $s>s_0$, then there is outburst of the bacterial population in the environment. Thus, to limit the net growth rate of bacteria in the environment, we must impose the condition $s_0>s$ in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ). Also, it may be noted that $s_1>s_2$. For epidemiological reasons, all the parameters involved in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) are considered to be positive, and are described in Table 2.

Table 2. Descriptions of parameters involved in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ).

Parameters	Descriptions
Λ	Rate of immigration into the class of susceptible population
β	Contact rate of susceptibles with infected individuals
η	Transmission rate of susceptible to the infected class due to their
	interaction with bacteria present in the environment
L	Half saturation constant
δ	Rate of transfer of recovered individuals to
	susceptible class due to immunity loss
d	Natural death rate of human population
ν	Recovery rate of human population
α	Disease-induced death rate of human population
s	Self-growth rate of bacteria
$s_0$	Natural death rate of bacteria
$s_1$	Growth rate coefficient of bacteria due to the discharged by infective
	in the environment
$s_2$	Efficacy of chemicals to reduce the density of
	bacteria discharged by infected individuals
k	Per-capita use of chemicals to reduce the density of bacteria
	at source discharged by infected individuals
M	Half saturation constant
${\pi }_1$	Decay coefficient in the density of bacteria present
	in the environment due to the uptake of chemicals
${\varphi }_1$	Per-capita reduction of chemicals due to uptake by bacteria
ϕ	Growth rate of chemicals due to increase in the density of bacteria
${\varphi }_0$	Natural depletion of chemical

At first, we show that the model system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) is well-possed. In this regard, we present the following theorem [43].

Theorem 2.1

Let ${\mathbb{R}}_{+}^n=[0,\mathrm{\infty }{)}^n$ be the cone of non-negative vectors in ${\mathbb{R}}^n$. Let $F:{\mathbb{R}}_{+}^{n+1}\to {\mathbb{R}}^n$ defined by $F(t,X)=\left(F_1\right(t,X),F_2(t,X),\dots ,F_n(t,X\left)\right),$ where $X=(X_1,X_2,\dots ,X_n)$ be locally Lipschitz and satisfy $F_j(t,X){|}_{X_j=0}\ge 0$ whenever $t\ge 0$ and $X\in {\mathbb{R}}_{+}^n$. Then, for every $X^0\in {\mathbb{R}}_{+}^n$, there exists a unique solution of $X^{\mathrm{\prime }}=F(t,X)$, $X\left(0\right)=X^0$ with values in ${\mathbb{R}}_{+}^n$, which is defined on some interval $[0,b)$, b>0. If $b<\mathrm{\infty }$, then $\underset{t\to b}{lim}sup\sum _{j=1}^n{X}_j\left(t\right)=\mathrm{\infty }.$

Proof.

To apply this theorem for system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ), we have to show that $F(t,X)$ is locally Lipschitz and $F_j(t,X){|}_{X_j=0}\ge 0$. If $s_1>s_2$, then the condition $F_j(t,X){|}_{X_j=0}\ge 0$ is trivially satisfied. As the right hand of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) is continuously differentiable on ${\mathbb{R}}_{+}^5$, so the locally Lipschitz condition trivially holds. Let $\Vert X\Vert \le r$, i.e. $|S|+|I|+|R|+|B|+|C_h|\le r$. Taking $F_1:=\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS$ as an example, we have $\begin{array}{rl}\left|\frac{\mathrm{\partial }{F}_1}{\mathrm{\partial }S}\right|& \le |\text{β}I|+|\eta |+|d|=\text{β}r+\eta +d=r_1,\\ \left|\frac{\mathrm{\partial }{F}_1}{\mathrm{\partial }I}\right|& \le |\text{β}S|=\text{β}r\le r_1,\\ \left|\frac{\mathrm{\partial }{F}_1}{\mathrm{\partial }B}\right|& \le \left|\frac{\eta L}{(L+B{)}^2}S\right|\le \frac{\eta }{L}r=r_2,\\ \left|\frac{\mathrm{\partial }{F}_1}{\mathrm{\partial }R}\right|& \le |\delta |=\delta .\end{array}$ Notably the partial derivative of F exists and is bounded in $\Vert X\Vert \le r$, i.e. F is locally Lipschitz. This implies that there exists a unique solution which is defined on $[0,b)$, b>0. Assume that $b<\mathrm{\infty }$, then $\underset{t\to b}{lim sup}\sum _{j=1}^n{X}_j\left(t\right)=\mathrm{\infty }.$ This contradicts the boundedness of $\Vert X\Vert$. Thus, $b=\mathrm{\infty }$, and hence there exists a unique, non-negative solution on $[0,\mathrm{\infty })$.

Now, using the fact that N = S + I + R, we reduce system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) to the following equivalent form: (2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) Now onwards, we study the dynamics of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) in detail.

3. Mathematical analysis of system (2)

3.1. Positivity and boundedness of solutions

Theorem 3.1

If $s_1>s_2$, system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) has a unique non-negative solution. The region of attraction for all solutions of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) initiating in the positive orthant is given by the following set: $\begin{array}{rl}\mathrm{\Omega }& =\left\{\right(I,R,N,B,C_h)\in {\mathbb{R}}_{+}^5:0\le I+R\le N\le \frac{\mathrm{\Lambda }}{d},0\le B\le \frac{s_1\mathrm{\Lambda }}{d{s}_3},\\ & 0\le C_h\le \frac{s_1\varphi \mathrm{\Lambda }}{d{s}_3({\varphi }_0+k)}\}.\end{array}$ With respect to model system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), the region Ω is compact and positively invariant. Further, it is closed and bounded in the positive cone of the five-dimensional space.

Proof.

We can rewrite system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) in the following form: $\frac{\mathrm{d}X}{\mathrm{d}t}=CX+D,$ where $X=[S,I,R,B,C_h{]}^T$ and $\begin{array}{r}C=\left(\begin{array}{ccccc}-C_{11}& 0& \delta & 0& 0\\ C_{21}& -C_{22}& 0& 0& 0\\ 0& \nu & -(\delta +d)& 0& 0\\ 0& C_{42}& 0& -C_{44}& 0\\ 0& 0& 0& \varphi & -C_{55}\end{array}\right)\end{array}$ with $\begin{array}{rl}& C_{11}=\text{β}I+\eta \frac{B}{L+B}+d,C_{21}=\text{β}I+\eta \frac{B}{L+B},C_{22}=\nu +\alpha +d,\\ & C_{42}=s_1-s_2\frac{k{C}_h}{M+k{C}_h},C_{44}=s_3+{\pi }_1{\varphi }_1{C}_h,C_{55}={\varphi }_0+k+{\varphi }_{1}B.\end{array}$ The vector $D=[\mathrm{\Lambda },0,0,0,0{]}^T$ is positive. Notably, all the off-diagonal entries of the matrix $C\left(X\right)$ are non-negative. Therefore, the matrix $C\left(X\right)$ is Metzler for every $X\in {\mathbb{R}}_{+}^5$. Thus, system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) is positively invariant in $X\in {\mathbb{R}}_{+}^5$ [1]. Hence, all the solution trajectories of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) that originate from initial states in $X\in {\mathbb{R}}_{+}^5$, confine therein forever.

From the third equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), we have $\frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{\Lambda }-dN-\alpha I\le \mathrm{\Lambda }-\mathrm{d}N.$ Using the standard comparison theorem [18], we have $0\le N\left(t\right)\le \frac{\mathrm{\Lambda }}{d}+\left(N\right(0)-\frac{\mathrm{\Lambda }}{d})e^{-dt}.$ If $N\left(0\right)<\frac{\mathrm{\Lambda }}{d},$ then the above inequality gives $\underset{t\to \mathrm{\infty }}{lim sup}N\left(t\right)\le \frac{\mathrm{\Lambda }}{d}.$ Hence, for any larger values of t, $0\le N\left(t\right)\le \frac{\mathrm{\Lambda }}{d}.$ Now, we have $\frac{\mathrm{d}(N-I-R)}{\mathrm{d}t}=\mathrm{\Lambda }-\text{β}(N-I-R)I-\eta \frac{B}{L+B}(N-I-R)-d(N-I-R)+\delta R.$ At N−I−R = 0, we have $\frac{\mathrm{d}(N-I-R)}{\mathrm{d}t}=\mathrm{\Lambda }+\delta R>0$. This implies that $N-I-R\ge 0$. Therefore, $0\le I+R\le N\le \frac{\mathrm{\Lambda }}{d}$ for any large time. Further, from the fourth equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), as $I\le \frac{\mathrm{\Lambda }}{d}$, we get $\frac{\mathrm{d}B}{\mathrm{d}t}\le s_1\frac{\mathrm{\Lambda }}{d}-s_{3}B$ for any large value of t. From the theory of differential inequality [18], we obtain $\underset{t\to \mathrm{\infty }}{lim sup}B\left(t\right)\le \frac{s_1\mathrm{\Lambda }}{d{s}_3}.$ This gives, $0\le B\left(t\right)\le \frac{s_1\mathrm{\Lambda }}{d{s}_3}$ for large value of t.

Furthermore, from the fifth equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), as $B\left(t\right)\le \frac{s_1\mathrm{\Lambda }}{d{s}_3}$ for large value of t, we get $\frac{\mathrm{d}{C}_h}{\mathrm{d}t}\le \varphi \frac{s_1\mathrm{\Lambda }}{d{s}_3}-({\varphi }_0+k)C_h.$ From the theory of differential inequality [18], we obtain $\underset{t\to \mathrm{\infty }}{lim sup}{C}_h\left(t\right)\le \frac{s_1\varphi \mathrm{\Lambda }}{d{s}_3({\varphi }_0+k)}.$ This yields that $0\le C_h\left(t\right)\le \frac{s_1\varphi \mathrm{\Lambda }}{d{s}_3({\varphi }_0+k)}$ for large time. Thus, we find that all solutions of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) enter into the region Ω, indicating that this region is an attracting set.

3.2. Disease-free equilibrium and its stability

The disease-free equilibrium corresponding to model system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) can be obtained by equating the right-hand side of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) to zero, when none of people in the considered domain is affected by the disease. For system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), the disease-free equilibrium is $E_0=(0,0,\frac{\mathrm{\Lambda }}{d},0,0)$. Apparently, the equilibrium $E_0$ always exists in the system. Now, we will establish the local stability of the equilibrium $E_0$ in terms of the basic reproduction number (${\mathcal{R}}_0$), a potential measure that determines whether a specific disease can invade a population or not.

3.2.1. Basic reproduction number

In order to obtain an expression for the basic reproduction number (${\mathcal{R}}_0$) associated with system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), we use the approach of the next generation matrix [47]. To this, we find the transmission (new infection terms) and transition terms in the infected subsystem of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) as follows: $\begin{array}{rl}& \mathcal{F}=\left(\begin{array}{c}{\displaystyle \text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)}\\ 0\end{array}\right),\\ & \mathcal{V}=\left(\begin{array}{c}(\nu +\alpha +d)I\\ {\displaystyle -(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I+s_{3}B+{\pi }_1{\varphi }_1{C}_{h}B}\end{array}\right).\end{array}$ Thus, the transmission (new infection) and transition matrices corresponding to the equilibrium $E_0$ are respectively obtained as $\begin{array}{r}F=\left(\begin{array}{cc}{\displaystyle \frac{\text{β}\mathrm{\Lambda }}{d}}& {\displaystyle \frac{\eta \mathrm{\Lambda }}{dL}}\\ 0& 0\end{array}\right),V=\left(\begin{array}{cc}\nu +\alpha +d& 0\\ -s_1& s_3\end{array}\right).\end{array}$ It has been well established that the basic reproduction number of an epidemic model is the spectral radius of the next generation matrix $F{V}^{-1}$, i.e. ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)$. Thus, for model system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), we get (3) ${\mathcal{R}}_0=\frac{\text{β}\mathrm{\Lambda }}{d(\nu +\alpha +d)}+\frac{\eta s_1\mathrm{\Lambda }}{s_{3}dL(\nu +\alpha +d)}.$(3)

Remark 3.1

Note that the expression in (3(3) ${\mathcal{R}}_0=\frac{\text{β}\mathrm{\Lambda }}{d(\nu +\alpha +d)}+\frac{\eta s_1\mathrm{\Lambda }}{s_{3}dL(\nu +\alpha +d)}.$(3) ) stands for the basic reproduction number associated with the model system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) and represents the average number of secondary infections produced by a single infected person in a group of completely susceptible people. From the expression of ${\mathcal{R}}_0$, the impacts of disease transmission through the direct contacts of susceptible with the infected individuals and with bacteria present in the environment are clearly visible. It expresses the facts that the disease spreads in the community through the direct as well as indirect routes of transmission. It is also apparent from the expression of ${\mathcal{R}}_0$ that reducing the net growth rate of bacteria in the environment can help in controlling the prevalence of the disease.

Following [47], the local stability of the disease-free equilibrium $E_0$ of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) is stated in the following theorem.

Theorem 3.2

The disease-free equilibrium $E_0$ is always feasible. It is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof.

It can be easily proved by evaluating the Jacobian matrix of the system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) at the disease-free equilibrium $E_0$ and finding the corresponding eigenvalues.

Regarding global stability of the disease-free equilibrium $E_0$, we have the following theorem.

Theorem 3.3

The disease-free equilibrium $E_0$ of model system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) is globally asymptotically stable in the region ${\mathbb{R}}_{+}^2$ of I−B plane if ${\mathcal{R}}_0<1$.

Proof.

We obtain the following subsystem from the system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ): $\begin{array}{rl}& \frac{\mathrm{d}I}{\mathrm{d}t}=\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I:=f_1,\\ & \frac{\mathrm{d}B}{\mathrm{d}t}=(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B:=f_2.\end{array}$ Consider, $h(I,B)=\frac{1}{IB}$ and ${\mathrm{△}}_{E_0}(I,B)=\frac{\mathrm{\partial }}{\mathrm{\partial }I}\left(h{f}_1\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }B}\left(h{f}_2\right)$. Since $h(I,B)>0$ for all I, B>0, therefore ${\mathrm{△}}_{E_0}(I,B)=\frac{\mathrm{\partial }}{\mathrm{\partial }I}\left(h{f}_1\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }B}\left(h{f}_2\right)=-[\frac{\text{β}}{B}+\frac{\eta (N-R)}{(L+B)I^2}+\frac{1}{B^2}(s_1-s_2\frac{k{C}_h}{M+k{C}_h})]<0.$ Clearly, ${\mathrm{△}}_{E_0}(I,B)$ does not change its sign and also it is not identically zero in the positive quadrant of the I−B plane. Thus, by the Bendixson–Dulac criterion [50], system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) has no limit cycle in the positive quadrant of the I−B plane. As disease-free equilibrium $E_0$ is locally asymptotically stable whenever ${\mathcal{R}}_0<1$, so it will be globally asymptotically stable in the region ${\mathbb{R}}_{+}^2$ of the I−B plane for ${\mathcal{R}}_0<1$. A similar procedure can be used to prove the global stability of the disease-free equilibrium $E_0$ in the other planes.

3.3. Endemic equilibrium and its stability

Let $E^{\ast }=(I^{\ast },R^{\ast },N^{\ast },B^{\ast },C_h^{\ast })$ be an endemic equilibrium of the epidemic model (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), whose components are positive solutions of the system's equilibrium equations. From the second equilibrium equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), we have (4) $R=\frac{\nu I}{\delta +d}.$(4) From the third equilibrium equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), we get (5) $N=\frac{\mathrm{\Lambda }-\alpha I}{d}.$(5) Next, from the fifth equilibrium equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), we have (6) $C_h=\frac{\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}.$(6) Using Equation (6(6) $C_h=\frac{\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}.$(6) ) in the fourth equilibrium equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), we have (7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) Note that the left-hand side (LHS) of the above equation is a decreasing whereas the right-hand side (RHS) is an increasing function of B. Thus, if there is a positive solution of Equation (7(7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) ), then it is unique in the interval $(0,\mathrm{\infty })$. Let $B\to 0$ from above, then the RHS of Equation (7(7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) ) becomes zero whereas the LHS will be positive and finite for $I\ne 0.$ This means, LHS>RHS. If $B\to \mathrm{\infty }$ from below, then the RHS of Equation (7(7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) ) becomes $\mathrm{\infty }$ whereas the LHS is positive and finite. In this case, LHS<RHS. Therefore, there exists a unique positive solution of Equation  (7(7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) ), say $B=f_B\left(I\right)$ in the interval $(0,\mathrm{\infty })$. We may rewrite Equation (7(7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) ) as (8) $B=f\left(B\right)I,$(8) where $\begin{array}{rl}& f\left(B\right)=\frac{f_1\left(B\right)}{f_2\left(B\right)},f_1\left(B\right)=\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B},\\ & f_2\left(B\right)=s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}.\end{array}$ On differentiating Equation (8(8) $B=f\left(B\right)I,$(8) ) with respect to I, we get $\frac{\mathrm{d}B}{\mathrm{d}I}=\frac{f\left(B\right)}{1-f^{\prime }\left(B\right)I}.$ Notably, $\frac{dB}{dI}>0$ provided $f^{\prime }\left(B\right)<0$. We can conclude from the above that $f_1^{\prime }\left(B\right)=\frac{-s_{2}k\varphi M({\varphi }_0+k)}{\left[M\right({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B{]}^2}<0,f_2^{\prime }\left(B\right)=\frac{({\varphi }_0+k){\pi }_1{\varphi }_1\varphi }{({\varphi }_0+k+{\varphi }_{1}B{)}^2}>0.$ This yields that $f^{\prime }\left(B\right)=\frac{f_2\left(B\right)f_1^{\prime }\left(B\right)-f_1\left(B\right)f_2^{\prime }\left(B\right)}{\left[f_2\right(B){]}^2}<0.$ Thus, we find that $\frac{dB}{dI}>0$ and $f^{\prime }\left(B\right)<0$. Notably, $S=N-I-R=\frac{\mathrm{\Lambda }-\alpha I}{d}-I-\frac{\nu I}{\delta +d}.$ Thus, for a positive value of S, we must have $I\le I_0=\frac{\mathrm{\Lambda }(\delta +d)}{(\alpha +d)(\delta +d)+\nu d}.$ Further, by using Equations (4(4) $R=\frac{\nu I}{\delta +d}.$(4) ), (5(5) $N=\frac{\mathrm{\Lambda }-\alpha I}{d}.$(5) ) and (8(8) $B=f\left(B\right)I,$(8) ) in the first equilibrium equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), and noting that $B=f_B\left(I\right)$ is an appropriate positive solution of Equation (7(7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) ), we get (9) $F\left(I\right)=(\text{β}+\frac{\eta f\left(B\right)}{L+f\left(B\right)I})(\frac{\mathrm{\Lambda }-(\alpha +d)I}{d}-\frac{\nu I}{\delta +d})-(\nu +\alpha +d).$(9) The following properties of Equation (9(9) $F\left(I\right)=(\text{β}+\frac{\eta f\left(B\right)}{L+f\left(B\right)I})(\frac{\mathrm{\Lambda }-(\alpha +d)I}{d}-\frac{\nu I}{\delta +d})-(\nu +\alpha +d).$(9) ) can be easily observed.

1. $F\left(0\right)=(\nu +\alpha +d)({\mathcal{R}}_0-1)$, which is positive provided ${\mathcal{R}}_0>1$.

2. $F\left(I_0\right)=-(\nu +\alpha +d)<0$.

3. $F^{\prime }\left(I\right)<0$ ∀$I\in (0,I_0)$.

Thus, there is a unique positive root (say $I=I^{\ast }$) of the equation $F\left(I\right)=0$ in the interval $(0,I_0)$ whenever ${\mathcal{R}}_0>1$. Putting $I=I^{\ast }$ in Equations (4(4) $R=\frac{\nu I}{\delta +d}.$(4) )–(7(7) $\frac{s_{1}M({\varphi }_0+k)+[{\varphi }_1{s}_{1}M+(s_1-s_2\left)k\varphi \right]B}{M({\varphi }_0+k)+(M{\varphi }_1+k\varphi )B}I=[s_3+\frac{{\pi }_1{\varphi }_1\varphi B}{{\varphi }_0+k+{\varphi }_{1}B}]B.$(7) ), we can obtain the positive values of $R^{\ast }$, $N^{\ast }$, $C_h^{\ast }$ and $B^{\ast }$. Thus, the equilibrium $E^{\ast }$ exists in the system provided ${\mathcal{R}}_0>1$ and $I\in (0,I_0)$.

Regarding local stability of the endemic equilibrium $E^{\ast }$, we have the following theorem.

Theorem 3.4

The endemic equilibrium $E^{\ast }$, if feasible, is locally asymptotically stable.

Proof.

The Jacobian matrix J at the endemic equilibrium $E^{\ast }$ becomes, $\begin{array}{r}{J}_{E^{\ast }}=\left(\begin{array}{ccccc}-b_{11}& -b_{12}& b_{13}& b_{14}& 0\\ \nu & -(\delta +d)& 0& 0& 0\\ -\alpha & 0& -d& 0& 0\\ b_{41}& 0& 0& -b_{44}& -b_{45}\\ 0& 0& 0& b_{54}& -b_{55}\end{array}\right),\end{array}$ where $\begin{array}{rl}& b_{11}=\text{β}{I}^{\ast }+\frac{\eta B^{\ast }}{L+B^{\ast }}\frac{(N^{\ast }-R^{\ast })}{I^{\ast }},b_{12}=b_{13}=\text{β}{I}^{\ast }+\frac{\eta B^{\ast }}{L+B^{\ast }},\\ & b_{14}=\frac{\eta L(N^{\ast }-I^{\ast }-R^{\ast })}{(L+B^{\ast }{)}^2},\\ & b_{41}=s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }},b_{44}=s_3+{\pi }_1{\varphi }_1{C}_h^{\ast },b_{45}=\frac{s_{2}kM{I}^{\ast }}{(M+k{C}_h^{\ast }{)}^2}+{\pi }_1{\varphi }_1{B}^{\ast },\\ & b_{54}=\varphi -{\varphi }_1{C}_h^{\ast },b_{55}={\varphi }_0+k+{\varphi }_1{B}^{\ast }.\end{array}$ For the matrix $J_{E^{\ast }}$, we get the following characteristic equation: (10) $(\psi +b_{11}+\frac{b_{12}\nu }{\psi +\delta +d}+\frac{b_{13}\alpha }{\psi +d})(\psi +b_{44}+\frac{b_{45}{b}_{54}}{\psi +b_{55}})=b_{14}{b}_{41}.$(10) If $\psi =\chi +i\rho$, $\chi \ge 0$ is a root of Equation (10(10) $(\psi +b_{11}+\frac{b_{12}\nu }{\psi +\delta +d}+\frac{b_{13}\alpha }{\psi +d})(\psi +b_{44}+\frac{b_{45}{b}_{54}}{\psi +b_{55}})=b_{14}{b}_{41}.$(10) ), then $\begin{array}{r}|\chi +i\rho +b_{11}+\frac{b_{12}\nu (\chi +\delta +d-i\rho )}{(\chi +\delta +d{)}^2+{\rho }^2}+\frac{b_{13}\alpha (\chi +d-i\rho )}{(\chi +d{)}^2+{\rho }^2}|\\ \times |\chi +i\rho +b_{44}+\frac{b_{45}{b}_{54}(\chi +b_{55}-i\rho )}{(\chi +b_{55}{)}^2+{\rho }^2}|=b_{14}{b}_{41}.\end{array}$ Note that $\begin{array}{rl}& |\chi +i\rho +b_{11}+\frac{b_{12}\nu (\chi +\delta +d-i\rho )}{(\chi +\delta +d{)}^2+{\rho }^2}+\frac{b_{13}\alpha (\chi +d-i\rho )}{(\chi +d{)}^2+{\rho }^2}|\\ & =\{{(\chi +b_{11}+\frac{b_{12}\nu (\chi +\delta +d)}{(\chi +\delta +d{)}^2+{\rho }^2}+\frac{b_{13}\alpha (\chi +d)}{(\chi +d{)}^2+{\rho }^2})}^2\\ & +{{\rho }^2{(1-\frac{b_{12}\nu }{(\chi +\delta +d{)}^2+{\rho }^2}-\frac{b_{13}\alpha }{(\chi +d{)}^2+{\rho }^2})}^2\}}^{1/2}\\ & \ge \chi +b_{11}+\frac{b_{12}\nu (\chi +\delta +d)}{(\chi +\delta +d{)}^2+{\rho }^2}+\frac{b_{13}\alpha (\chi +d)}{(\chi +d{)}^2+{\rho }^2}\ge \chi +b_{11}\ge b_{11}.\end{array}$ Similarly, we can show that $|\chi +i\rho +b_{44}+\frac{b_{45}{b}_{54}(\chi +b_{55}-i\rho )}{(\chi +b_{55}{)}^2+{\rho }^2}|\ge b_{44}.$ Thus, the $|LHS|$ of Equation  (10(10) $(\psi +b_{11}+\frac{b_{12}\nu }{\psi +\delta +d}+\frac{b_{13}\alpha }{\psi +d})(\psi +b_{44}+\frac{b_{45}{b}_{54}}{\psi +b_{55}})=b_{14}{b}_{41}.$(10) ) $>b_{11}{b}_{44}.$ This shows that for $b_{11}{b}_{44}>b_{14}{b}_{41}$, Equation (10(10) $(\psi +b_{11}+\frac{b_{12}\nu }{\psi +\delta +d}+\frac{b_{13}\alpha }{\psi +d})(\psi +b_{44}+\frac{b_{45}{b}_{54}}{\psi +b_{55}})=b_{14}{b}_{41}.$(10) ) does not possess any root with positive real part. Note that the inequality $b_{11}{b}_{44}>b_{14}{b}_{41}$ holds provided $(\text{β}{I}^{\ast }+\frac{\eta B^{\ast }}{L+B^{\ast }}\frac{(N^{\ast }-R^{\ast })}{I^{\ast }})(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })>\frac{\eta L(N^{\ast }-I^{\ast }-R^{\ast })}{(L+B^{\ast }{)}^2}(s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }}).$ Obviously, $1>\frac{L}{L+B^{\ast }}$ and $(N^{\ast }-R^{\ast })>(N^{\ast }-I^{\ast }-R^{\ast })$. Thus, $\begin{array}{rl}& (N^{\ast }-R^{\ast })>\frac{L(N^{\ast }-I^{\ast }-R^{\ast })}{L+B^{\ast }}\\ & \Rightarrow \frac{\eta (N^{\ast }-R^{\ast })}{L+B^{\ast }}>\frac{\eta L(N^{\ast }-I^{\ast }-R^{\ast })}{(L+B^{\ast }{)}^2}\\ & \Rightarrow \text{β}\frac{(I^{\ast }{)}^2}{B^{\ast }}+\frac{\eta (N^{\ast }-R^{\ast })}{L+B^{\ast }}>\frac{\eta L}{(L+B^{\ast }{)}^2}(N^{\ast }-I^{\ast }-R^{\ast })\\ & \Rightarrow (\text{β}{I}^{\ast }+\frac{\eta (N^{\ast }-R^{\ast })}{L+B^{\ast }}\frac{B^{\ast }}{I^{\ast }})\frac{I^{\ast }}{B^{\ast }}>\frac{\eta L(N^{\ast }-I^{\ast }-R^{\ast })}{(L+B^{\ast }{)}^2}.\end{array}$ The fourth equilibrium equation of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) yields that $\frac{I^{\ast }}{B^{\ast }}=\frac{s_3+{\pi }_1{\varphi }_1{C}_h^{\ast }}{s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }}}$. This together with the above equation implies $b_{11}{b}_{44}>b_{14}{b}_{41}$. Thus, the $|LHS|$ of Equation (10(10) $(\psi +b_{11}+\frac{b_{12}\nu }{\psi +\delta +d}+\frac{b_{13}\alpha }{\psi +d})(\psi +b_{44}+\frac{b_{45}{b}_{54}}{\psi +b_{55}})=b_{14}{b}_{41}.$(10) ) $>b_{11}{b}_{44}>b_{14}{b}_{41}.$ The above discussions imply that Equation (10(10) $(\psi +b_{11}+\frac{b_{12}\nu }{\psi +\delta +d}+\frac{b_{13}\alpha }{\psi +d})(\psi +b_{44}+\frac{b_{45}{b}_{54}}{\psi +b_{55}})=b_{14}{b}_{41}.$(10) ) does not have any solution with positive real part. Therefore, it can be concluded that the endemic equilibrium $E^{\ast }$, if exists, is locally asymptotically stable.

3.4. Direction of bifurcation at ${\mathcal{R}}_0=1$

Theorem 3.2 indicates that there is an exchange of stability between equilibria $E_0$ and $E^{\ast }$ at the critical value ${\mathcal{R}}_0=1$. Thus, there is a possibility of transcritical bifurcation at ${\mathcal{R}}_0=1$. To determine the nature of bifurcation, involving equilibrium $E_0$ at ${\mathcal{R}}_0=1$, we employ central manifold theory [7].

Without loss of generality, we introduce the new variables $x_1=I$, $x_2=R$, $x_3=N$, $x_4=B$ and $x_5=C_h$. With these, system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) gets the following form: (11) $\begin{array}{rl}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}& =\text{β}(x_3-x_1-x_2)x_1+\eta \frac{x_4}{L+x_4}(x_3-x_1-x_2)-(\nu +\alpha +d)x_1,\\ \frac{\mathrm{d}{x}_2}{\mathrm{d}t}& =\nu x_1-(\delta +d)x_2,\\ \frac{\mathrm{d}{x}_3}{\mathrm{d}t}& =\mathrm{\Lambda }-d{x}_3-\alpha x_1,\\ \frac{\mathrm{d}{x}_4}{\mathrm{d}t}& =(s_1-s_2\frac{k{x}_5}{M+k{x}_5})x_1-s_3{x}_4-{\pi }_1{\varphi }_1{x}_4{x}_5,\\ \frac{\mathrm{d}{x}_5}{\mathrm{d}t}& =\varphi x_4-({\varphi }_0+k)x_5-{\varphi }_1{x}_4{x}_5.\end{array}$(11) Note that at ${\mathcal{R}}_0=1$, $\text{β}={\text{β}}^{\ast }=\frac{d(\nu +\alpha +d)}{\mathrm{\Lambda }}-\frac{\eta s_1}{L{s}_3}$. The disease-free equilibrium associated with system (11(11) $\begin{array}{rl}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}& =\text{β}(x_3-x_1-x_2)x_1+\eta \frac{x_4}{L+x_4}(x_3-x_1-x_2)-(\nu +\alpha +d)x_1,\\ \frac{\mathrm{d}{x}_2}{\mathrm{d}t}& =\nu x_1-(\delta +d)x_2,\\ \frac{\mathrm{d}{x}_3}{\mathrm{d}t}& =\mathrm{\Lambda }-d{x}_3-\alpha x_1,\\ \frac{\mathrm{d}{x}_4}{\mathrm{d}t}& =(s_1-s_2\frac{k{x}_5}{M+k{x}_5})x_1-s_3{x}_4-{\pi }_1{\varphi }_1{x}_4{x}_5,\\ \frac{\mathrm{d}{x}_5}{\mathrm{d}t}& =\varphi x_4-({\varphi }_0+k)x_5-{\varphi }_1{x}_4{x}_5.\end{array}$(11) ) is $x_0=(x_{01},x_{02},x_{03},x_{04},x_{05})=(0,0,\mathrm{\Lambda }/d,0,0).$ The linearized matrix corresponding to system (11(11) $\begin{array}{rl}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}& =\text{β}(x_3-x_1-x_2)x_1+\eta \frac{x_4}{L+x_4}(x_3-x_1-x_2)-(\nu +\alpha +d)x_1,\\ \frac{\mathrm{d}{x}_2}{\mathrm{d}t}& =\nu x_1-(\delta +d)x_2,\\ \frac{\mathrm{d}{x}_3}{\mathrm{d}t}& =\mathrm{\Lambda }-d{x}_3-\alpha x_1,\\ \frac{\mathrm{d}{x}_4}{\mathrm{d}t}& =(s_1-s_2\frac{k{x}_5}{M+k{x}_5})x_1-s_3{x}_4-{\pi }_1{\varphi }_1{x}_4{x}_5,\\ \frac{\mathrm{d}{x}_5}{\mathrm{d}t}& =\varphi x_4-({\varphi }_0+k)x_5-{\varphi }_1{x}_4{x}_5.\end{array}$(11) ) around the equilibrium $x_0$ at $\text{β}={\text{β}}^{\ast }$ is obtained as $\begin{array}{r}J(x_0,{\text{β}}^{\ast })=\left[\begin{array}{ccccc}{\displaystyle -\frac{\eta s_1\mathrm{\Lambda }}{dL{s}_3}}& 0& 0& {\displaystyle \frac{\eta \mathrm{\Lambda }}{dL}}& 0\\ \nu & -(\delta +d)& 0& 0& 0\\ -\alpha & 0& -d& 0& 0\\ s_1& 0& 0& -s_3& 0\\ 0& 0& 0& \varphi & -({\varphi }_0+k)\end{array}\right].\end{array}$ The above matrix has a simple eigenvalue zero and all other eigenvalues are negative at ${\mathcal{R}}_0=1$.

If $w=(w_1,w_2,w_3,w_4,w_5{)}^T$ is the right eigenvector corresponding to the zero eigenvalue at $(x_0,{\text{β}}^{\ast })$, then $w_1=1,w_2=\frac{\nu }{\delta +d},w_3=-\frac{\alpha }{d},w_4=\frac{s_1}{s_3},w_5=\frac{\varphi s_1}{s_3({\varphi }_0+k)}.$ Similarly, from the matrix $J(x_0,{\text{β}}^{\ast })$, we obtain a left eigenvector $v=(v_1,v_2,v_3,v_4,v_5{)}^T$ (associated with the zero eigenvalue) with components as $v_1=s_3,v_2=v_3=0,v_4=\frac{\eta \mathrm{\Lambda }}{dL},v_5=0.$ Now, we compute the following quantities: (12) $a=\sum _{k,i,j=1}^5{v}_k{w}_i{w}_j\frac{{\mathrm{\partial }}^2{f}_k}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}(x_0,{\text{β}}^{\ast }),b=\sum _{k,i=1}^5{v}_k{w}_i\frac{{\mathrm{\partial }}^2{f}_k}{\mathrm{\partial }{x}_i\mathrm{\partial }\text{β}}(x_0,{\text{β}}^{\ast })$(12) to get, (13) $\begin{array}{rl}a& =-2s_3[\frac{d(\nu +\alpha +d)}{\mathrm{\Lambda }}(1+\frac{\alpha }{d}+\frac{\nu }{\delta +d})+\frac{\eta \mathrm{\Lambda }}{d}{\left(\frac{s_1}{L{s}_3}\right)}^2]\\ & -\frac{2\eta \mathrm{\Lambda }}{dL}[\frac{\varphi s_1}{({\varphi }_0+k)s_3}\frac{s_{2}k}{M}+\frac{\varphi {\pi }_1{\varphi }_1}{{\varphi }_0+k}{\left(\frac{s_1}{s_3}\right)}^2],\\ b& =\frac{s_3\mathrm{\Lambda }}{d}.\end{array}$(13) Now, we present the following theorem.

Theorem 3.5

Consider model system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), and let a and b be as in (13(13) $\begin{array}{rl}a& =-2s_3[\frac{d(\nu +\alpha +d)}{\mathrm{\Lambda }}(1+\frac{\alpha }{d}+\frac{\nu }{\delta +d})+\frac{\eta \mathrm{\Lambda }}{d}{\left(\frac{s_1}{L{s}_3}\right)}^2]\\ & -\frac{2\eta \mathrm{\Lambda }}{dL}[\frac{\varphi s_1}{({\varphi }_0+k)s_3}\frac{s_{2}k}{M}+\frac{\varphi {\pi }_1{\varphi }_1}{{\varphi }_0+k}{\left(\frac{s_1}{s_3}\right)}^2],\\ b& =\frac{s_3\mathrm{\Lambda }}{d}.\end{array}$(13) ), where a<0 and b>0. The local dynamics of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) around the equilibrium $E_0$ can be stated as: when ${\mathcal{R}}_0<1$ with ${\mathcal{R}}_0\approx 1$, the equilibrium $E_0$ is locally asymptotically stable, and there exists a negative unstable equilibrium $E^{\ast }$. On the other hand, for ${\mathcal{R}}_0>1$ with ${\mathcal{R}}_0\approx 1$, the equilibrium $E_0$ is unstable and there exists a positive locally asymptotically stable equilibrium $E^{\ast }$. This signifies that at ${\mathcal{R}}_0=1$, system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) showcases a supercritical (forward) bifurcation.

Proof.

The proof follows from Theorem 4.1 p. 373 and Remark 1 p. 375 in [7].

The next theorem is for the global asymptotic stability of the endemic equilibrium $E^{\ast }$.

Theorem 3.6

The endemic equilibrium $E^{\ast }$, if feasible, is globally asymptotically stable in the region of attraction Ω provided, (14) $\begin{array}{rl}\left(\frac{\eta \mathrm{\Lambda }}{d(L+B^{\ast })I^{\ast }}\right)(s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }})& <\frac{1}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}),\end{array}$(14) (15) $\begin{array}{rl}{(\frac{s_{2}k\mathrm{\Lambda }}{d(M+k{C}_h^{\ast })}+\frac{{\pi }_1{\varphi }_1{s}_1\mathrm{\Lambda }}{d{s}_3})}^2& <\frac{1}{4}\frac{(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })({\varphi }_0+k)}{(\varphi -{\varphi }_1{C}_h^{\ast }{)}^2}.\end{array}$(15)

Proof.

To establish the global stability of the endemic equilibrium $E^{\ast }$, we consider a positive definite function as follows: (16) $\begin{array}{rl}G& =(I-I^{\ast }-I^{\ast }\ln \frac{I}{I^{\ast }})+\frac{m_1}{2}(R-R^{\ast }{)}^2+\frac{m_2}{2}(N-N^{\ast }{)}^2+\frac{m_3}{2}(B-B^{\ast }{)}^2\\ & +\frac{m_4}{2}(C_h-C_h^{\ast }{)}^2.\end{array}$(16) Here, $m_i$'s $(i=1-4)$ are some positive constants that will be determined later.

The time derivative of G along the solution of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ), for $m_1=\frac{1}{\nu }(\text{β}+\frac{\eta B^{\ast }}{(L+B^{\ast })I^{\ast }})$ and $m_2=\frac{1}{\alpha }(\text{β}+\frac{\eta B^{\ast }}{(L+B^{\ast })I^{\ast }}),$ is obtained as $\begin{array}{rl}\frac{\mathrm{d}G}{\mathrm{d}t}& =-(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}+\frac{\eta B(N-I-R)}{(L+B)I{I}^{\ast }})(I-I^{\ast }{)}^2\\ & -\frac{(\delta +d)}{\nu }(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }})(R-R^{\ast }{)}^2-\frac{d}{\alpha }(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }})(N-N^{\ast }{)}^2\\ & -m_3(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(B-B^{\ast }{)}^2-m_4({\varphi }_0+k+{\varphi }_{1}B\left)\right(C_h-C_h^{\ast }{)}^2\\ & +\frac{\eta L(N-I-R)}{(L+B)(L+B^{\ast })I^{\ast }}(I-I^{\ast })(B-B^{\ast })+m_3(s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }})(I-I^{\ast })(B-B^{\ast })\\ & -m_3(\frac{s_{2}kMI}{(M+k{C}_h)(M+k{C}_h^{\ast })}+{\pi }_1{\varphi }_{1}B)(B-B^{\ast })(C_h-C_h^{\ast })\\ & +m_4(\varphi -{\varphi }_1{C}_h^{\ast })(B-B^{\ast })(C_h-C_h^{\ast }).\end{array}$ Clearly, following inequalities ensure the negative definite of the derivative $\frac{\mathrm{d}G}{\mathrm{d}t}$ inside the region Ω: (17) $\begin{array}{rl}& {\left(\frac{\eta \mathrm{\Lambda }}{d(L+B^{\ast })I^{\ast }}\right)}^2<\frac{m_3}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}),\end{array}$(17) (18) $\begin{array}{rl}& m_3{(s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }})}^2<\frac{1}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}),\end{array}$(18) (19) $\begin{array}{rl}& m_3{(\frac{s_{2}k\mathrm{\Lambda }}{d(M+k{C}_h^{\ast })}+\frac{{\pi }_1{\varphi }_1{s}_1\mathrm{\Lambda }}{d{s}_3})}^2<\frac{m_4}{2}({\varphi }_0+k)(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast }),\end{array}$(19) (20) $\begin{array}{rl}& m_4(\varphi -{\varphi }_1{C}_h^{\ast }{)}^2<\frac{m_3}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast }\left)\right({\varphi }_0+k).\end{array}$(20) From inequalities (17(17) $\begin{array}{rl}& {\left(\frac{\eta \mathrm{\Lambda }}{d(L+B^{\ast })I^{\ast }}\right)}^2<\frac{m_3}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}),\end{array}$(17) ) and (18(18) $\begin{array}{rl}& m_3{(s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }})}^2<\frac{1}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}),\end{array}$(18) ), a positive value of $m_3$ can be chosen if the inequality (14(14) $\begin{array}{rl}\left(\frac{\eta \mathrm{\Lambda }}{d(L+B^{\ast })I^{\ast }}\right)(s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }})& <\frac{1}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}),\end{array}$(14) ) holds. Also, from inequalities (19(19) $\begin{array}{rl}& m_3{(\frac{s_{2}k\mathrm{\Lambda }}{d(M+k{C}_h^{\ast })}+\frac{{\pi }_1{\varphi }_1{s}_1\mathrm{\Lambda }}{d{s}_3})}^2<\frac{m_4}{2}({\varphi }_0+k)(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast }),\end{array}$(19) ) and (20(20) $\begin{array}{rl}& m_4(\varphi -{\varphi }_1{C}_h^{\ast }{)}^2<\frac{m_3}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast }\left)\right({\varphi }_0+k).\end{array}$(20) ), we obtain the positive value of $m_4$ in view of condition (15(15) $\begin{array}{rl}{(\frac{s_{2}k\mathrm{\Lambda }}{d(M+k{C}_h^{\ast })}+\frac{{\pi }_1{\varphi }_1{s}_1\mathrm{\Lambda }}{d{s}_3})}^2& <\frac{1}{4}\frac{(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })({\varphi }_0+k)}{(\varphi -{\varphi }_1{C}_h^{\ast }{)}^2}.\end{array}$(15) ). Thus, the derivative $\frac{\mathrm{d}G}{\mathrm{d}t}$ will be negative definite inside the region of attraction Ω if both the conditions (14(14) $\begin{array}{rl}\left(\frac{\eta \mathrm{\Lambda }}{d(L+B^{\ast })I^{\ast }}\right)(s_1-s_2\frac{k{C}_h^{\ast }}{M+k{C}_h^{\ast }})& <\frac{1}{2}(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })(\text{β}+\eta \frac{B^{\ast }}{(L+B^{\ast })I^{\ast }}),\end{array}$(14) ) and (15(15) $\begin{array}{rl}{(\frac{s_{2}k\mathrm{\Lambda }}{d(M+k{C}_h^{\ast })}+\frac{{\pi }_1{\varphi }_1{s}_1\mathrm{\Lambda }}{d{s}_3})}^2& <\frac{1}{4}\frac{(s_3+{\pi }_1{\varphi }_1{C}_h^{\ast })({\varphi }_0+k)}{(\varphi -{\varphi }_1{C}_h^{\ast }{)}^2}.\end{array}$(15) ) hold good.

4. Numerical simulations

Now, to support the analytically obtained results in the previous section and to get a better understanding of the behavior of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ), we perform some numerical simulations. Unless it is mentioned in the text, the parameter values utilized for the numerical observations of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) are the same as provided in Table 3.

Table 3. Values of parameters used for numerical simulations of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ).

Parameters	Units	Values	References
Λ	persons day${}^{-1}$	5	[29]
β	person${}^{-1}$day${}^{-1}$	0.000004	[29, 45]
η	day${}^{-1}$	0.0005	[36]
L	cells mm${}^{-3}$	1000	[36]
δ	day${}^{-1}$	0.002	Assumed
ν	day${}^{-1}$	0.2	[29, 45]
α	day${}^{-1}$	0.00001	[29, 45]
d	day${}^{-1}$	0.00004	[29, 45]
s	day${}^{-1}$	0.02	[36]
$s_0$	day${}^{-1}$	0.08	[36]
$s_1$	cells mm${}^{-3}$person${}^{-1}$day${}^{-1}$	0.05	[46]
$s_2$	cells mm${}^{-3}$person${}^{-1}$day${}^{-1}$	0.01	Assumed
k	day${}^{-1}$	0.1	Assumed
M	ml day${}^{-1}$	10	Assumed
${\pi }_1$	cells mm${}^{-3}$ ml${}^{-1}$	2	Assumed
${\varphi }_1$	mm${}^3$ cell${}^{-1}$ day${}^{-1}$	0.01	Assumed
ϕ	ml mm${}^3$ cell${}^{-1}$ day${}^{-1}$	0.05	Assumed
${\varphi }_0$	day${}^{-1}$	0.01	Assumed

4.1. Sensitivity analysis

With respect to some epidemiologically important model parameters, we calculate the basic reproduction number's normalized forward sensitivity indices following the method described in [25]. The ratio of the relative change in a variable to the relative change in a parameter is given by the normalized forward sensitivity index of that variable to the parameter. For instance, the normalized forward sensitivity index of a variable γ that relies differentiably on a parameter δ is defined by $X_{\delta }^{\gamma }=\frac{\mathrm{\partial }\gamma }{\mathrm{\partial }\delta }\times \frac{\delta }{\gamma }$. Figure 2 depicts the normalized forward sensitivity indices of the basic reproduction number ${\mathcal{R}}_0$ with respect to some of the model parameters. The figure demonstrates that the value of ${\mathcal{R}}_0$ grows with an increment in the value of either of the parameters Λ, β, η, $s_1$ and s as each of them has positive index with ${\mathcal{R}}_0$. Apparently, the figure depicts that the model parameters L, ν and $s_0$ have negative indices with ${\mathcal{R}}_0$. So, increments in these model parameters can lead to a decline in the value of ${\mathcal{R}}_0$. The figure also shows that the sensitivity indices of ${\mathcal{R}}_0$ with respect to the parameters Λ and ν are $+1$ and $-1$, respectively. This represents the fact that 1% rise in the value of Λ (ν) will result in 1% increment (decrement) in the value of ${\mathcal{R}}_0$. Importantly, a lesser value of ${\mathcal{R}}_0$ is preferable for control of the disease. This suggests to prevent increment in the values of parameters that have positive indices with ${\mathcal{R}}_0$ whereas parameters having negative indices with ${\mathcal{R}}_0$ should be encouraged. Hence, the policymakers must focus on the prevention strategies that reduce (enhance) the values of parameters having positive (negative) sensitivity indices with the basic reproduction number.

Figure 2. The normalized forward sensitivity indices of the basic reproduction number (${\mathcal{R}}_0$) with respect to the model parameters Λ, β, η, L, ν, $s_1$, $s_0$ and s. Here, the values of parameters are chosen from Table 3.

It is worthy to note that the parameter values used for simulations of the mathematical models may have some errors as they are determined from experiments. To tame this uncertainty in the selection of model parameters, we perform global sensitivity analysis in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) by implementing Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCCs) [6, 24]. Note that LHS permits to differ several model parameters simultaneously in an efficient way, whereas PRCC correlates the model input parameters with the response function by giving values between $-1$ and 1. Here, we have taken the response function as the infected population (I). Notably, the signs of PRCCs indicate the type of correlations between the input parameters and the output variable while the values represent the strength of their influences. We run 1000 simulations per LHS by letting uniform distributions for each of the input parameter. The baseline values of parameters are taken from Table 3 with $\pm 25\mathrm{\%}$ deviations, and the PRCC values are plotted as bar graphs in Figure 3. From the figure, we see that the model parameters Λ, β, η, $s_1$ and s have positive correlations with the infected population in the epidemic region. The figure also depicts that the infected population is negatively correlated with the parameters L, ν and $s_0$. All the parameters with significant correlations with the infected population are marked by ${}^{\ast }$ on the respective bars. The sensitivity results indicate that boosting up the parameters with negative PRCCs with the infected population can help in reducing the epidemic burden. Besides, controlling the parameters with positive PRCC values can play an important role in the disease eradication.

Figure 3. The uncertainty of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) on infected individuals (I). Baseline values of parameters are the same as in Table 3 except $\eta =0.005$. Significant parameters are marked by $\ast$ (p-value $<0.05$).

4.2. Impacts of some key parameters on disease control

We further explore the effects of the parameters $s_0$, s and ν on the basic reproduction number (${\mathcal{R}}_0$). By varying $s_0-s$ and ν simultaneously, we plot ${\mathcal{R}}_0$ as contour line in Figure 4. It is apparent from the figure that if the difference between $s_0$ and s increases, the basic reproduction number decreases. We can see from the figure that for the lower values of $s_0-s$, the basic reproduction number (${\mathcal{R}}_0$) is greater than unity. However, if the difference $s_0-s$ is above a certain value, ${\mathcal{R}}_0$ becomes less than unity, where the endemic equilibrium loses its feasibility and the disease-free equilibrium becomes feasible. It is also clear from the figure that for smaller values of ν, the value of ${\mathcal{R}}_0$ is higher and greater than one. But, ${\mathcal{R}}_0$ becomes less than unity after a certain range of ν. Moreover, we see from the figure that simultaneous increments in the values of ν and the difference $s_0-s$ decrease the value of ${\mathcal{R}}_0$, and push it back to a value smaller than one after certain ranges of these two. Notably, the disease will persist among the population for lower values of ν and $s_0-s$, but if their values increase, ${\mathcal{R}}_0$ becomes less than unity, and the endemicity will be lost in the population. Thus, it can be concluded from the figure that the prevalence of disease can be reduced by increasing the recovery rate of infected population and the decay rate of bacteria in the epidemic region. Importantly, attention must be paid on the reduction of the self growth rate of bacteria.

Figure 4. Contour plots of the basic reproduction number (${\mathcal{R}}_0$) as a functions of $s_0-s$ and ν. Rest of the parameters are at the same values as in Table 3. In the figure, dashed green line stands for ${\mathcal{R}}_0=1$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Now, we show the presence of forward bifurcation in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ). In view of expression (3(3) ${\mathcal{R}}_0=\frac{\text{β}\mathrm{\Lambda }}{d(\nu +\alpha +d)}+\frac{\eta s_1\mathrm{\Lambda }}{s_{3}dL(\nu +\alpha +d)}.$(3) ), we first write the model parameter β in term of the basic reproduction number (${\mathcal{R}}_0$). Then, we vary ${\mathcal{R}}_0$ in the interval [0.01,2.5], and plot the equilibrium number of infected population in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ), see Figure 5. In the figure, the blue line denotes the stable disease-free equilibrium ($E_0$), the red line represents unstable disease-free equilibrium ($E_0$) and the magenta line stands for the stable endemic equilibrium ($E^{\ast }$). The figure demonstrates that for ${\mathcal{R}}_0<1$, system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) settles at the stable disease-free steady state. But, as ${\mathcal{R}}_0$ exceeds the unit value, the endemic equilibrium $E^{\ast }$ bifurcates into the positive quadrant and showcases stable configuration whereas the disease-free equilibrium loses its stability. This indicates that the disease-free and endemic equilibria of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) are interlinked via a transcritical bifurcation. This manifests that the epidemic will persist (or die out) in the population whenever ${\mathcal{R}}_0>1$ (or ${\mathcal{R}}_0<1$). Thus, Figure 5 numerically illustrates the analytical findings of Theorem 3.2.

Figure 5. Bifurcation diagram of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) with respect to ${\mathcal{R}}_0$. Parameters are at the same values as in Table 3 except $\alpha =0.001$, d = 0.004, $s_2=0.04$ and ${\pi }_1=10$. Here, blue, red and magenta colours, respectively, denote the stable disease-free equilibrium ($E_0$), unstable disease-free equilibrium ($E_0$) and the stable endemic equilibrium ($E^{\ast }$). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In Figures 6 and  7, we simultaneously vary two parameters in model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ), and see their impacts on the infected population and bacterial density in the epidemic region. Here, the pairs of parameters are taken as $(\text{β},s_2)$, $(\eta ,{\pi }_1)$, $(\varphi ,{\varphi }_0)$ and $(s_1,k)$. Indeed, these model parameters play a significant role in either spread or control, of the disease. In the figures, the two surfaces represent the maximum and minimum values of the variables on the simultaneous variations of these parameters. If two surfaces are visible in the pictures, then fluctuations occur in the equilibrium values of the populations while when the two surfaces collide, a stable situation is attained. The impacts of direct and indirect ways of transmissions of disease on the population levels are apparent from Figure 6. We can see in Figure 6(a) that for lower values of β, the equilibrium values of both infected and bacterial populations become zero. But, their equilibrium values increase as the value of β increases. Additionally, we see that after a certain value of β, the infected population attains a plateau. The parameter η possesses similar effect on the persistence of disease in the population (see Figure 6(b)). The figure depicts that up to a certain range of η, the disease wipes out from the population. However, the disease persists if the value of η exceeds a fixed value. The number of infected persons increases up to a certain range of η, and it attains a plateau as η surpasses a certain value. Figure 7(a) shows that by adding more disinfectants for the control of bacteria, the number of infected population can be reduced. Further, it shows that the number of infected persons rapidly decreases with an increase in the growth rate of chemicals due to increased bacterial density in the environment. Indeed, if the chemicals are being used at a high rate, then the bacteria population dies out from the environment and the disease will be under control. An increased growth rate of bacteria is dangerous for the human population because it causes a rapid incline in the number of infectives. This behaviour can be seen in Figure 7(b). Overall, Figures 6 and 7 demonstrate that the increased disease transmission rates due to direct and indirect routes lead to an accelerated growth in the epidemic size. But, in that case, controlling the bacterial population through the use of chemicals can play a pivotal role in the disease eradication.

Figure 6. Surface plots of the infective population (first column) and bacterial density (second column) with respect to (a) β and $s_2$, and (b) η and ${\pi }_1$. Other parameters are at the same values as in Table 3.

Figure 7. Surface plots of the infective population (first column) and bacterial density (second column) with respect to (a) ϕ and ${\varphi }_0$, and (b) $s_1$ and k. Other parameters are at the same values as in Table 3.

Now, we see the role of disinfectants in the control of disease in the population. To this, we plot the infective population (I) and the density of bacteria (B) for four different combinations of $s_2$ and ${\pi }_1$, Figure 8(a). It can be seen from the figure that the infective population and bacteria population attain high equilibrium values when $s_2={\pi }_1=0$. But, for $s_2=0$ and ${\pi }_1=1$, the epidemic curve gets reduced. We see further decrements in the infective population and the bacterial density for $s_2=0.01$ and ${\pi }_1=4$. Moreover, we observe that for $s_2=0.04$ and ${\pi }_1=10$, the infection wipes out from the community and the bacteria disappear from the environment. The figure showcases that the efficacy of chemicals to reduce the density of bacteria discharged by infected individuals plays a major role in the control of bacterial density and elimination of the disease. Finally, we see the impacts of chemicals in reducing the density of bacteria at the source discharged by the infected individuals. A variation plot for the infected individuals (I) and the bacterial density (B) in Figure 8(b) show that the disease persists in the population for $k={\varphi }_1=0$. However, the number of infectives and the bacterial density fall on increasing the values of k and ${\varphi }_1$ to 0.5 and 0.0001, respectively. A comparatively more decrease in the infected population can be seen from the figure on increasing the value of ${\varphi }_1$ to 0.01. Importantly, the figure shows that both the infective and bacteria populations can be completely controlled for $k={\varphi }_1=0.8$. Thus, one can say that the disease can be totally controlled by enhancing the death of bacteria through the use of disinfectants.

Figure 8. Variations in the infective population (first column) and bacterial density (second column) with respect to time for different combinations of (a) $s_2$ and ${\pi }_1$, and (b) k and ${\varphi }_1$. Rest of the parameters are at the same values as in Table 3 except $\text{β}=0.000003$.

Figure 9 illustrates the global stability of the disease-free equilibrium $E_0$ in two different spaces. For this figure, we have chosen $\alpha =0.001$, d = 0.004, $s_2=0.04$ and ${\pi }_1=14$ while the values of other parameters are taken from Table 3. For this parametric set-up, we obtain the value of basic reproduction number as ${\mathcal{R}}_0=0.027<1$. In Figure 9(a), we observe that the solution trajectories of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) starting from four distinct values converge to the components of the disease-free equilibrium $E_0$ in I−R−B space. Likewise, Figure 9(b) demonstrates the global stability of the disease-free equilibrium $E_0$ in $I-R-C_h$ space. The figures suggest that for such parametric values, the disease can be wiped out from the population for a longer period of time.

Figure 9. Figure illustrating the global asymptotic stability of disease-free equilibrium $E_0$ in (a) I−R−B and (b) $I-R-C_h$ spaces. Parameters are at the same values as in Table 3 except $\alpha =0.001$, d = 0.004, $s_2=0.04$ and ${\pi }_1=14$. Rest of the parameters are at the same values as in Table 3.

For the endemic equilibrium $E^{\ast }$, the global stability is illustrated in Figure 10. For this figure, we take the values of parameters from Table 3 for which ${\mathcal{R}}_0=2.76$. Figure 10(a) clearly demonstrates that solution trajectories of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) starting within the region of attraction Ω converge to the components of the equilibrium $E^{\ast }$ in $I-R-C_h$ space. Similarly, we observe in Figure 10(b) that the solution trajectories of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\text{β}SI-\eta \frac{B}{L+B}S+\delta R-dS,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}SI+\eta \frac{B}{L+B}S-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =sB-s_{0}B+(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(1) ) originating in the region Ω ultimately converge to the components of the endemic equilibrium $E^{\ast }$ in I−R−B space. The global stability of the equilibrium point $E^{\ast }$ can also be shown in other spaces using phase diagrams. Thus, the analytical result of Theorem 3.6 is justified in Figure 10. Thus, one can say that for such values of model parameters, the disease persists among the population.

Figure 10. Global stability of the endemic equilibrium $E^{\ast }$ in (a) $I-R-C_h$ (b) I−R−B spaces.

5. Conclusion

Here, we have explored the impact of chemicals in controlling the diseases that spread in the population due to direct contacts between infective and the susceptible individuals, and also via bacteria present in the environment. The model equations described the time evolutions of susceptible human population, infected human population, recovered human population, bacteria present in the environment and chemical concentration of disinfectants/sanitizers. Mathematically, we showed that the proposed model exhibits disease-free and endemic equilibria which are connected via a forward bifurcation. We observed that for ${\mathcal{R}}_0<1$, system (2(2) $\begin{array}{rl}\frac{\mathrm{d}I}{\mathrm{d}t}& =\text{β}(N-I-R)I+\eta \frac{B}{L+B}(N-I-R)-(\nu +\alpha +d)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\nu I-(\delta +d)R,\\ \frac{\mathrm{d}N}{\mathrm{d}t}& =\mathrm{\Lambda }-dN-\alpha I,\\ \frac{\mathrm{d}B}{\mathrm{d}t}& =(s_1-s_2\frac{k{C}_h}{M+k{C}_h})I-s_{3}B-{\pi }_1{\varphi }_1{C}_{h}B,\\ \frac{\mathrm{d}{C}_h}{\mathrm{d}t}& =\varphi B-{\varphi }_0{C}_h-k{C}_h-{\varphi }_1{C}_{h}B.\end{array}$(2) ) exhibits a globally stable disease-free equilibrium. On the other hand, a unique endemic equilibrium exists in the system for ${\mathcal{R}}_0>1$. Moreover, sufficient conditions are derived under which the endemic equilibrium is globally asymptotically stable. These results advise to focus on some control strategies that lower the epidemic threshold below unity. This can be achieved by controlling the rate of immigration into the epidemic region and transmissions of disease due to direct/indirect contacts of the susceptible individuals with the infected ones and the bacteria. Importantly, accelerating the recovery rate of infected individuals and the death rate of bacteria present in the environment can be effective ways for controlling the disease.

We have implemented some sophisticated simulation techniques to discover the influences of some epidemiologically important model parameters on the disease dynamics. Our sensitivity results showed that the parameters Λ, β, η, $s_1$ and s have positive correlations with the infected individuals. Thus, increments in the values of these parameters can boost up the disease prevalence. On the other hand, the parameters L, ν and $s_0$ have negative influences on the infective cases. These results indicate that controlling the transmission rates of disease due to the contacts between susceptibles and infectives, and via bacteria present in the environment can help in reducing the disease prevalence. At the same time, fostering the rate of recovery and the natural death rate of bacteria can play significant roles in disease eradication. Additionally, we found that controlling the bacterial population at the sources discharged by the infected persons significantly reduces the epidemic size. Our simulation results showed the importance of disinfectants in the control of disease among the human population. The results evoked that the disinfectants of high efficacy can help in controlling the bacterial density and thereby eradicate the disease from the region of epidemic outbreak. Overall, the findings of this study suggest that using disinfectants of high quality can play a crucial role in the control of bacterial diseases in the human population.

Acknowledgments

All the authors express their gratitude to the associate editor and the reviewers whose comments and suggestions have helped in improving this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Rabindra Kumar Gupta is thankful to UGC Nepal for partial financial support in the form of ‘PhD Fellowship and Research Support’ (No. PhD/76-77 S&T-17).