Dynamics of a nonlinear epidemic transmission model incorporating a class of hospitalized individuals: a qualitative analysis and simulation

Whenever an epidemic outbreak occurs in society, it is of paramount importance for the government and healthcare workers to take swift action to control the spread of the disease and provide medical treatment to infected individuals. Mathematical modeling plays a vital role in understanding and predicting the spread of the disease, with numerous models such as SIS, SIR, SIRS, SEIR, SEIRS, SEIQ, SVIRS, ${\text{S}}{{\text{I}}_{\text{D}}}{\text{AIR}}$, and others that have been proposed in the epidemiology literature [1–15]. These models use different variables such as susceptible, vaccinated, infected, exposed, quarantined, infodemic, and recovered individuals to represent different stages of disease progression. The dynamics of these models are largely determined by the rates of incidence, contact, and recovery. Amongst these, the incidence rate, which represents the rate of new infections, is of critical importance. Generally, the incidence rate is proportionate to the number of susceptible and infectious people. In 1927, Kermack and McKendrick [16] proposed a SIR epidemic model with a bilinear incidence rate $\left( {\beta SI} \right)$; however, this model is ineffective in cases of overcrowded infectives. The bilinear incidence rate assumes that the contact rate and infection probability per contact remain constant over time, which does not reflect reality. The number of infective individuals greatly influences individual mobility as it represents the risk of infection. As an example, during the SARS outbreak in 2003, as the number of SARS-infected individuals increased, various protective measures and control policies such as closing schools, restaurants, and postponing conferences were implemented. These measures significantly reduced the number of contacts per unit of time, thus lowering the incidence rate [17]. Hence, it is essential to consider the infection force, which includes individuals' adaptation to infection risks, in the modeling of disease dynamics. Therefore, to address the limitations of traditional models and more accurately capture the complexities of disease transmission dynamics, Capasso and Serio [18] introduced a revolutionary concept in 1978: the saturated incidence rate $k\left( I \right)S$ into the disease transmission model. This innovative approach takes into account the real-world phenomenon of saturation, where the rate of infection tends to level off as the number of infected individuals increases. The saturated transmission infection force, represented by the equation $k\left( I \right) = \frac{{\alpha I}}{{1 + \beta I}}$, effectively accounts for the decrease in contact rate and infection probability per contact as the number of infective individuals increases. This provides a more accurate representation of the actual spread of the disease. The $\,\alpha I$ in $k\left( I \right)$ estimates the infection force, while $\frac{1}{{1 + \beta I}}$ measures the impact of behavioral changes in susceptible individuals due to increasing numbers or the density effect of infective individuals. The idea of incorporating the saturated incidence rate $k\left( I \right)S$ into the disease transmission model represents a significant advancement in the field of epidemiology. It is a more sensible approach than the traditional bilinear incidence rate $\beta SI$ [9] as it accounts for the crucial real-world phenomenon of behavioral change and the crowding effect of infective individuals. This approach also retains the unboundedness of the contact rate by choosing suitable parameters, providing a more accurate representation of the actual spread of the disease. Furthermore, the work by Capasso and Serio has been widely adopted by many researchers in the field, with numerous studies applying the saturated incidence rate to disease transmission dynamics in various epidemic models (see, for example [6, 19, 20]).

The recovery rate of infected individuals is a crucial factor in controlling the spread of disease and minimizing its impact on the population. This rate is primarily determined by the length of the recovery process, which is, in turn, related to the availability of good health infrastructure, such as the number of doctors, nurses, medicines, isolation areas, and available hospital beds per million populations. The availability of these resources is a key consideration for health planners in determining the level of treatment that can be provided to patients [21]. Unfortunately, limited health infrastructure or resources can have dire consequences, such as an increase in new infection cases, a rise in the death rate, and psychological effects on individuals' mental health. To characterize the saturation phenomenon of limited health infrastructure, a continually differentiable hospitalization rate function (Holling functional type II) is considered, represented by the equation $G\left( I \right) = \frac{{\gamma I}}{{1\, + \,\delta I}}$, where γ and δ are non-negative constants. The term $\gamma I$ measures the hospitalization rate of severely infected individuals, whereas the term $\frac{1}{{1 + \delta I}}$ gauges the limitation in the hospitalization of infected individuals due to the limited number of hospitals, isolation places, equipment, limitations of hospitals with required facilities, etc. The Holling type II functional form has been widely adopted in disease modeling. In 2008, Zhang and Liu [22] considered the treatment rate in saturated function form within their study. They proposed a SIR model using the following system of differential equations to describe the disease transmission dynamics:

$$\begin{equation*}\frac{{dS}}{{dt}} = A - dS - \frac{{\lambda SI}}{{1 + kI}},\end{equation*}$$

$$\begin{equation*}\frac{{dI}}{{dt}} = \frac{{\lambda SI}}{{1 + kI}} - \left( {d + \epsilon + \mu } \right)I - \frac{{rI}}{{1 + \alpha I}},\end{equation*}$$

$$\begin{equation*}\frac{{dR}}{{dt}} = \mu I + \frac{{rI}}{{1 + \alpha I}} - dI.\end{equation*}$$

In the above system of equations, $S,\,I,$ and $R$ stand for the susceptible, infected, and recovered populations, respectively whereas the parameters $A,\,d,\epsilon ,\,\mu ,\,\lambda ,\,k,\,r,$ and $\alpha$ represent the constant recruitment rate of the susceptible population, the natural mortality rate, the disease-induced mortality rate, the natural recovery rate of infected individuals, the transmission rate of infection, the psychological effects or inhibitory effects, the cure rate, and the extent of the effect of delayed treatment for the infected individuals, respectively. Moreover, $f\left( {S,I} \right) = \frac{{\lambda SI}}{{1 + kI}}$, and $T\left( I \right) = \frac{{rI}}{{1 + \alpha I}}$ represent the saturated incidence and treatment rates, respectively. Subsequently, in 2012, Zhou and Fan [23] modified the model of [22] by transforming the saturated treatment rate $T\left( I \right)$ into the form $\frac{{rI}}{{\alpha + I}}$ to better represent the realistic aspect of capturing the effects of limited medical resources for infected individuals and they studied the qualitative behavior of the model with this revised term. Despite these advancements, previous studies have not addressed the significance of the hospitalized population class and its impact on disease transmission dynamics. It is more reasonable to incorporate a dedicated class into the model, as the inclusion of this class enables a more precise representation of disease propagation and progression. This approach can facilitate tracking the population under hospital care at any time, the mortality count during hospital treatment, and the population that recovers after receiving medical treatment in hospitals. Therefore, this study aims to bridge this gap by introducing a novel model that incorporates a distinct class for hospitalized individuals within the SIR (susceptible-infected-recovered) framework, utilizing a saturated hospitalization rate. Furthermore, this representation acknowledges the potential discrepancy in growth rates between individuals requiring hospitalization and those becoming infected. Aligning with the evolving understanding of disease transmission and the role of hospitals in disease management enhances the effectiveness of this approach.

The goal of this study is to introduce and mathematically analyze an epidemic model that considers the hospitalized population class for severely infected individuals, limitations in the hospitalization of the infected population, and inhibitory effects by the population for disease transmission dynamics. In the present study, the model is proposed by dividing the total population into four subpopulation compartments: susceptible, infected, hospitalized, and recovered individuals, based on their disease status. By incorporating a saturated incidence rate and a saturated hospitalization rate, the standard SIR model is extended into the $SI{Q_H}R$ (susceptible-infected-hospitalized-recovered) model. We investigate the basic properties of the model solution. Moreover, the identification of equilibrium points and the computation of the basic reproduction number (BRN), denoted as ${R_0}$ [24, 25], are performed through Descartes's rule of signs and the Next-generation matrix method, respectively. The study examines the local and global stabilities of the equilibria using ${R_0}$, the Routh–Hurwitz criterion, and the Lyapunov function method. Additionally, we perform a bifurcation analysis for the model's stability behavior. Furthermore, numerical results are provided to corroborate the theoretical findings.

The paper is structured as follows: section 2 presents the proposed mathematical model and its properties. Section 3 examines the local stability of the model equilibria, while section 4 specifically analyzes the global stability of the endemic equilibrium (EE) due to the study's revelation of a backward bifurcation at ${R_0} = 1$. Section 5 supports the theoretical findings with numerical simulations, and section 6 explores the sensitivity of the BRN. Finally, section 7 concludes the paper with a summary of the findings and a discussion of their implications.

In this section, we present a nonlinear epidemic transmission model for disease transmission dynamics. To achieve this, the total population size $\left( P \right)$ is divided into four distinct subpopulations: susceptible, infected, hospitalized, and recovered. Each subpopulation represents a different stage in the progression of the disease. The model considers the complex dynamics of disease transmission, assuming that the disease can only be transmitted to susceptibles through close contact. Furthermore, individuals who have recovered and gained immunity no longer contribute to the transmission process. The susceptible subpopulation comprises individuals susceptible to contracting the infection through close contact with infected individuals. Conversely, infected individuals possess the capacity to transmit the disease under specific circumstances. The hospitalized subpopulation encompasses individuals who suffer severe illness due to the infection, necessitating specialized care for recovery. Recovered individuals, conversely, have successfully fought off the infection and regained their health, either through auto-recovery or medical treatment. To illustrate the process of disease transmission among these subpopulations, a block diagram is presented in figure 1.

Zoom In Zoom Out Reset image size

The transitions between subpopulations are mathematically expressed by the following system of ordinary differential equations, under the aforementioned assumptions:

$$\begin{align} \frac{{dS}}{{dt}} & = p - \mu S - \frac{{\alpha SI}}{{1 + \beta I}}, \nonumber \\ \frac{{dI}}{{dt}} & = \frac{{\alpha SI}}{{1 + \beta I}} - \frac{{\gamma I}}{{1 + \delta I}} - \left( {\mu + {d_1} + \theta } \right)I, \nonumber \\ \frac{{d{Q_H}}}{{dt}} & = \frac{{\gamma I}}{{1 + \delta I}} - a{Q_H} - \left( {\mu + {d_2}} \right){Q_H}, \nonumber \\ \frac{{dR}}{{dt}} & = a{Q_H} + \theta I - \mu R. \nonumber \\ \end{align} \tag{ 1 }$$

with initial conditions:

$$\begin{equation}S\left( 0 \right) > 0,\,I\left( 0 \right) \unicode{x2A7E} 0,\,{Q_H}\left( 0 \right) \unicode{x2A7E} 0,\,R\left( 0 \right) \unicode{x2A7E} 0.\end{equation} \tag{ 2 }$$

In model (1), the rate of change in each subpopulation is represented by $\frac{d}{{dt}}$, and the constant recruitment rate into the susceptible population is denoted by $p$. The transmission rate of infection is represented by $\alpha$, while the rate of inhibition by the infected population is denoted by $\beta$. The natural death rate in all subpopulations is symbolized by $\mu$, and the hospitalization rate of severely ill infected individuals is indicated by $\gamma$. The rate of limitation in the hospitalization of critically ill patients due to inadequate health infrastructure is denoted by $\delta$. The disease-induced death rate of infected individuals is denoted by ${d_1}$, and the disease-induced death rate of infected individuals in hospitals during treatment is denoted by ${d_2}$. The recovery rate of hospitalized individuals after medical treatment is represented by $a$, and the auto-recovery rate of infected individuals is represented by $\theta$. The term $K\left( {S,I} \right) = \frac{{\alpha SI}}{{1 + \beta I}}$ represents the saturated incidence rate at which susceptibles become infected when they come into close contact with infected individuals. In this term, $\alpha I$ represents the force of infection, and $\frac{1}{{1 + \beta I}}$ measures the inhibitory effects. If we consider no inhibitory effects (i.e. $\beta = 0$), then $K\left( {S,I} \right)$ simplifies to $\alpha SI$, which is a well-known bilinear incidence rate [16, 24]. The term $G\left( I \right) = \frac{{\gamma I}}{{1 + \delta I}}$ represents the saturated hospitalization rate of severely infected individuals, where the parameters $\gamma$ and $\delta$ are both non-negative. Below table 1 is given for a comprehensive description of the variables and parameters of the model (1).

Table 1. Notation and corresponding explanations.

Symbol	Description
$S$	Susceptible population.
$I$	Infected population.
${Q_H}$	Hospitalized population.
$R$	Recovered population.
$p$	Constant recruitment rate.
$\mu$	Natural death rate.
$\alpha$	Transmission rate of infection.
$\beta$	Inhibition rate or effect.
$\gamma$	Hospitalization rate of severely infected individuals.
$\delta$	Limitation rate for hospitalization of infected individuals.
${d_1}$	Mortality rate due to disease among infected people.
${d_2}$	Mortality rate due to disease among hospitalized people.
$\theta$	Auto-recovery rate of infected individuals.
$a$	Recovery rate of hospitalized individuals after medical treatment.

The variables and parameters in model (1) are assumed to be non-negative and positive, respectively, in line with biological reasons. Furthermore, we have identified a region of attraction for the solutions of model (1) as described in lemma 1.

Lemma 1. The set ${{\Gamma }} = \left\{ {\left( {S,\,I,\,{Q_H},\,R} \right) \in \mathbb{R}_ + ^4:0 < S + I + {Q_H} + R \unicode{x2A7D} \frac{p}{\mu }} \right\}$ is a positive invariant set of the system (1) under condition (2).

Proof. By adding all equations of the model (1), we obtain

$$\begin{align}& \frac{{dS}}{{dt}} + \frac{{dI}}{{dt}} + \frac{{d{Q_H}}}{{dt}} + \frac{{dR}}{{dt}} = p - \mu \left( {S + I + {Q_H} + R} \right) - {d_1}I - {d_2}{Q_H}\\ & \quad \Rightarrow \frac{d}{{dt}}\left( {S + I + {Q_H} + R} \right) = p - \mu \left( {S + I + {Q_H} + R} \right) - {d_1}I - {d_2}{Q_H}\\ & \quad \Rightarrow \frac{{dP}}{{dt}} = p - \mu P - {d_1}I - {d_2}{Q_H}, \quad\left( {\because P = S + I + {Q_H} + R} \right)\\ & \quad \Rightarrow \frac{{dP}}{{dt}} \unicode{x2A7D} p - \mu P.\end{align}$$

On integrating, we get

$$\begin{equation*}P \unicode{x2A7D} P\left( 0 \right){e^{ - \mu t}} + \frac{p}{\mu }\left( {1 - {e^{ - \mu t}}} \right).\end{equation*}$$

Therefore, $\mathop {\lim }\limits_{t \to \infty } Sup\,P \unicode{x2A7D} \frac{p}{\mu }$. Further, $\frac{{dP}}{{dt}} < 0$ if $P > \frac{p}{\mu }$. It shows that solutions of the model (1) approach ${{\Gamma }}$, ${{\Gamma }}$ is a positively invariant set and solutions of the model (1) are bounded, all provided $P > \frac{p}{\mu }$.

Excluding the fourth equation from the aforementioned model (1) for mathematical analysis does not have a significant impact on the overall system dynamics, as variable $R$ is absent in the first three equations of the model (1). Therefore, it is sufficient to focus on the following reduced system for further examination:

$$\begin{align} \frac{{dS}}{{dt}} & = p - \mu S - \frac{{\alpha SI}}{{1 + \beta I}}, \nonumber\\ \frac{{dI}}{{dt}} & = \frac{{\alpha SI}}{{1 + \beta I}} - \frac{{\gamma I}}{{1 + \delta I}} - \left( {\mu + {d_1} + \theta } \right)I, \nonumber\\ \frac{{d{Q_H}}}{{dt}} & = \frac{{\gamma I}}{{1 + \delta I}} - a{Q_H} - \left( {\mu + {d_2}} \right){Q_H}. \end{align} \tag{ 3 }$$

We now determine the equilibria of the model system (3) and analyze their stability behaviors.□

In this section, we examine the dynamics of the model by analyzing its equilibrium points. To do this, we set the right-hand side of all equations in the system (3) to zero, which yields the following two equilibria:

• (i)  
  Disease-free equilibrium (DFE) $\left( {{E_1} = \left( {\frac{p}{\mu },0,0} \right)} \right)$.
• (ii)  
  EE $\left( {{E_2} = \left( {{S^*},{I^*},Q_H^*} \right)} \right)$.

To investigate equilibria stability, we first compute the BRN [25] for our model using the next-generation matrix method [5, 25], as shown below:

3.1.  Computation of the BRN ${{R}_0}$

Consider the equation:

$$\begin{equation*}\dot G = F\left( x \right) - U\left( x \right).\end{equation*}$$

Here, $x = {\left( {I,{Q_H},S} \right)^T},\,F\left( x \right)$ represents the matrix of new infection terms, and $U\left( x \right)$ represents the matrix of transfer terms into and out of the compartments. The Jacobian matrices $\mathcal{F}$ and $\mho$ of $F\left( x \right)$ and $U\left( x \right)$, respectively, at ${E_1}\left( {\frac{p}{\mu },0,0} \right)$ are calculated as follows:

$$\begin{equation*}\mathcal{F} = \left( {\begin{array}{*{20}{c}} {\frac{{p\alpha }}{\mu }}&0&0 \\ 0&0&0 \\ 0&0&0 \end{array}} \right),\,\mho = \left( {\begin{array}{*{20}{c}} {{d_1} + \gamma + \theta + \mu }&0&0 \\ { - \gamma }&{a + {d_2} + \mu }&0 \\ {\frac{{p\alpha }}{\mu }}&0&\mu \end{array}} \right).\end{equation*}$$

Now,

$$\begin{align*}\mathcal{F}{\mho ^{ - 1}} & = \left( {\begin{array}{*{20}{c}} {\frac{{p\alpha }}{\mu }}&0&0 \\ 0&0&0 \\ 0&0&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\frac{{a\mu + {d_2}\mu + {\mu ^2}}}{{\mu \left( {a + {d_2} + \mu } \right)\left( {{d_1} + \gamma + \theta + \mu } \right)}}}&0&0 \\ {\frac{\gamma }{{\left( {a + {d_2} + \mu } \right)\left( {{d_1} + \gamma + \theta + \mu } \right)}}}&{\frac{1}{{a + {d_2} + \mu }}}&0 \\[4pt] { - \frac{{p\alpha }}{{{\mu ^2}\left( {{d_1} + \gamma + \theta + \mu } \right)}}}&0&{\frac{1}{\mu }} \end{array}} \right)\; \\ & = \left( {\begin{array}{*{20}{c}} {\frac{{p\alpha }}{{\mu \left( {{d_1} + \gamma + \theta + \mu } \right)}}}&0&0 \\ 0&0&0 \\ 0&0&0 \end{array}} \right).\end{align*}$$

The BRN, denoted as ${R_0}$, is defined as the spectral radius of the next-generation matrix $\mathcal{F}{\mho ^{ - 1}}$ [25]. Therefore, for our model, the BRN is:

$$\begin{equation*}{R_0} = \rho \left( {\mathcal{F}{\mho ^{ - 1}}} \right) = \frac{{p\alpha }}{{\mu \left( {{d_1} + \gamma + \theta + \mu } \right)}}.\end{equation*}$$

3.2.  Local stability behavior of $\,\,{{{E}}}_1$

In this subsection, the local stability behavior of DFE $\left( {{E_1}} \right)$ is investigated via the corresponding linearized matrix $\left( {{M_0}} \right)$ of the system (3) at ${E_1} = \left( {\frac{p}{\mu },0,0} \right)$. The linearized matrix $\left( {{M_0}} \right)$ is as given below:

$$\begin{equation*}{M_0} = \left( {\begin{array}{*{20}{c}} { - \mu }&{ - \frac{{p\alpha }}{\mu }}&0 \\ 0&{\left( {\mu + \theta + \gamma + {d_1}} \right)\left( {{R_0} - 1} \right)}&0 \\ 0&\gamma &{ - a - {d_2} - \mu } \end{array}} \right).\end{equation*}$$

The characteristic equation of ${M_0}$ is given by

$$\begin{equation}\left( {\lambda + \mu } \right)\left( {\lambda + a + {d_2} + \mu } \right)\left( {\lambda + \left( {1 - {R_0}} \right)\left( {\gamma + \theta + \mu + {d_1}} \right)} \right) = 0.\end{equation} \tag{ 4 }$$

It can be easily observed from (4) that two roots of equation (4) are negative, namely ${\lambda _1} = - \mu$ and ${\lambda _2} = - \left( {a + {d_2} + \mu } \right)$. The sign of the third root, ${\lambda _3}$, depends on the value of ${R_0}$. Consequently, if ${R_0} < 1$, then ${\lambda _3}$ is negative, and if ${R_0} > 1$, then ${\lambda _3}$ is positive. Based on this analysis and applying the Routh–Hurwitz criterion, we present the following theorem to describe the local stability behavior of ${E_1}$:

Theorem 1. When ${R_0} < 1$, the disease-free equilibrium $\left( {{E_1}} \right)$ of the system (3) is locally asymptotically stable. Conversely, if ${R_0} > 1$, the disease-free equilibrium is unstable.

Now, we discuss the stability behavior of ${E_1}$ when ${R_0} = 1$.

As the linearization technique cannot be used to determine the stability behavior of ${E_1}$ when ${R_0} = 1$, the center manifold theory approach outlined in [26, 27] is employed to investigate the stability behavior. To apply center manifold theory, it is supposed that $S = {x_1},\,I = {x_2}$ and ${Q_H} = {x_3}$, then, the system (3) is:

$$\begin{align} \frac{{d{x_1}}}{{dt}} & = p - \mu {x_1} - \frac{{\alpha {x_1}{x_2}}}{{1 + \beta {x_2}}} \equiv {\phi _1}, \nonumber \\ \frac{{d{x_2}}}{{dt}} & = \frac{{\alpha {x_1}{x_2}}}{{1 + \beta {x_2}}} - \frac{{\gamma {x_2}}}{{1 + \delta {x_2}}} - \left( {\mu + {d_1} + \theta } \right){x_2} \equiv {\phi _2}, \nonumber \\ \frac{{d{x_3}}}{{dt}} & = \frac{{\gamma {x_2}}}{{1 + \delta {x_2}}} - a{x_3} - \left( {\mu + {d_2}} \right){x_3} \equiv {\phi _3}. \nonumber \\ \end{align} \tag{ 5 }$$

The linearized matrix $M_0^*$ evaluated at ${R_0} = 1$ and bifurcation parameter $\alpha = {\alpha ^*} = \frac{{\mu \left( {{d_1} + \gamma + \theta + \mu } \right)}}{p}$ is obtained as

$$\begin{equation*}M_0^* = \left( {\begin{array}{*{20}{c}} { - \mu }&{ - \frac{{p{\alpha ^*}}}{\mu }}&0 \\ 0&{\,0}&0 \\ 0&\gamma &{ - a - {d_2} - \mu } \end{array}} \right).\end{equation*}$$

Let $u = {\left( {{u_1},{u_2},{u_3}} \right)^T}$ denotes the left eigenvector and presents the right eigenvector of $M_0^*$ corresponding to the zero eigenvalue. Then, we obtain

$$\begin{equation*}{u_1} = 0,\,{u_2} = 1,\,{u_3} = 0\;\;{\text{and}}\;\;{w_1} = - \frac{{p{\alpha ^*}}}{{{\mu ^2}}},\,{w_2} = 1,\,{w_3} = \frac{\gamma }{{a + {d_2} + \mu }}.\end{equation*}$$

At ${R_0} = 1$ and $\alpha = {\alpha ^*}$, the following non-zero partial derivatives, associated with ${\phi _2}$ of the system (5) are calculated as:

$$\begin{equation*}{\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial {x_2}\partial {x_1}}}} \right)_{{E_1}}} = {\alpha ^*} = {\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial {x_1}\partial {x_2}}}} \right)_{{E_1}}}\,,\,{\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial x_2^2}}} \right)_{{E_1}}} = 2\left( {\gamma \delta - \frac{{p{\alpha ^*}\beta }}{\mu }} \right),\,{\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial {x_2}\partial {\alpha ^*}}}} \right)_{{E_1}}} = \frac{p}{\mu }.\end{equation*}$$

The bifurcation coefficients ${C_*}$ and ${C_{**}}$ defined in theorem 4.1 of [26] are given by:

$$\begin{equation*}{c_*} = \mathop \sum \limits_{k,i,j = 1}^3 {u_k}{w_i}{w_j}{\left( {\frac{{{\partial ^2}{\phi _k}}}{{\partial {x_i}\partial {x_j}}}} \right)_{{E_1}}}\end{equation*}$$

$$\begin{align*}{c_*} & = {u_2}\left( {{w_1}{w_2}{{\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial {x_1}\partial {x_2}}}} \right)}_{{E_1}}} + {w_2}{w_1}{{\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial {x_2}\partial {x_1}}}} \right)}_{{E_1}}} + w_2^2{{\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial x_2^2}}} \right)}_{{E_1}}}} \right)\\ & = \left( {2\,{\alpha ^*}{w_1}{w_2} + 2\,w_2^2\left( {\gamma \delta - \frac{{p{\alpha ^*}\beta }}{\mu }} \right)} \right) \nonumber\\ & = - 2\left( {\frac{{p{\alpha ^*}^2}}{{{\mu ^2}}} + \frac{{p{\alpha ^*}\beta }}{\mu } - \gamma \delta } \right)\end{align*}$$

$$\begin{equation*}{c_{**}} = \mathop \sum \limits_{k,i = 1}^3 {u_k}{w_i}{\left( {\frac{{{\partial ^2}{\phi _k}}}{{\partial {x_i}\partial {\alpha ^*}}}} \right)_{{E_1}}}\end{equation*}$$

$$\begin{align*}{c_{**}} & = {w_2}{\left( {\frac{{{\partial ^2}{\phi _2}}}{{\partial {x_2}\partial {\alpha ^*}}}} \right)_{{E_1}}} \nonumber\\ & = \frac{p}{\mu } > 0. \end{align*}$$

As shown above, the bifurcation coefficient ${c_{**}}$ is positive. Accordingly, based on [26], the behavior of local stability at ${R_0} = 1$ depends on the sign of the second bifurcation constant ${c_*}$. If ${c_*}$ is negative, the system (3) demonstrates a forward transcritical bifurcation. Conversely, if ${c_*}$ is positive, the system (3) exhibits a backward transcritical bifurcation at ${R_{0}} = 1$.

Understanding forward and backward bifurcation is crucial in epidemiology. These concepts help us understand disease transmission dynamics and enhance disease control strategies. In epidemiology, these bifurcation types show distinct behaviors that profoundly affect disease management. Forward bifurcation becomes pivotal when the BRN $\left( {{R_0}} \right)$ is below unity. In such cases, eradicating the disease is possible if ${R_0}$ remains under unity. In contrast, with backward bifurcation, reducing ${R_0}$ below unity does not necessarily eliminate the disease. This difference significantly impacts disease control planning. The emergence of a backward bifurcation at ${R_0} = 1$ considerably complicates disease control. In the case of forward bifurcation at ${R_0} = 1$, if ${R_0}$ is below unity, the disease cannot establish itself within the population. Even introducing a few infected individuals leads to the number of infected reverting to zero. However, when ${R_0}$ is slightly above unity, the situation changes. The DFE $\left( {{E_1}} \right)$ shifts from stable to unstable, indicating controlled yet persistent disease within the population. On the contrary, backward bifurcation behaves differently. When ${R_0}$ is below unity, a small, unstable, positive EE arises, while the DFE and a larger EE remain locally stable. There are instances where both a stable DFE and a larger, stable positive equilibrium coexist. However, if ${R_0}$ exceeds 1, an unstable DFE and a stable positive EE coexist [28]. Hence, determining the range for bifurcation occurrence is crucial.

To investigate the region for the occurrence of the forward and backward transcritical bifurcations; the bifurcation coefficient ${c_*}$ is evaluated at ${\alpha ^*} = \frac{{\mu \left( {{d_1} + \gamma + \theta + \mu } \right)}}{p}$, which gives:

$$\begin{equation*}{c_*} = - \frac{2}{p}\left( {{\gamma ^2} + {c_1}\gamma + {c_0}} \right),\end{equation*}$$

where,

$$\begin{equation*}{C_0} = \left( {{{{d}}_1} + \theta + \mu } \right)\left( {{{{d}}_1} + p\beta + \theta + \mu } \right),\,{C_1} = \left( {2{{{d}}_1} + p\beta + 2\theta + 2\mu - p\delta } \right),\,\,{C_0} > 0.\end{equation*}$$

Let

$$\begin{equation}H\left( \gamma \right) = {\gamma ^2} + {c_1}\gamma + {c_0}.\end{equation} \tag{ 6 }$$

$\gamma _1^*$ and $\gamma _2^*$ be the two roots of $H\left( \gamma \right) = 0,$ and $D$ be the discriminant. Then, we obtain

$$\begin{align*}D = C_1^2 - 4{C_0} = p\left( {p{{\left( {\beta - \delta } \right)}^2} - 4\delta \left( {{{\text{d}}_1} + \theta + \mu } \right)} \right);\,\gamma _1^* = \frac{{ - {C_1} - \sqrt D }}{2}\,,\gamma _2^* = \frac{{ - {C_1} + \sqrt D }}{2}.\end{align*}$$

Hence, from [26], the following results are stated in the form of a theorem:

Theorem 2. The system (3) at ${R_0} = 1$ and bifurcation parameter $\alpha = {\alpha ^*} > 0$ has the following transcritical bifurcations:

• (i)  
  A backward transcritical bifurcation if $\left\{ {\begin{array}{*{20}{c}} {{C_1} < 0,} \\ {D > 0,} \\ {\gamma _1^* < \gamma < \gamma _2^*.} \end{array}} \right.$
• (ii)  
  A forward transcritical bifurcation if ${C_1} > 0\,$ or $\left\{ {\begin{array}{*{20}{c}} {{C_1} < 0,} \\ {D > 0,} \\ {\gamma < \gamma _1^*\,} {\text{or}}\,{\gamma > \gamma _2^*.\,} \end{array}} \right.$

To numerically validate the results of theorem 2, simulations are conducted using the experimental parameter values listed in table 2:

Table 2. The parameter's numerical value.

Parameter	$p$	$\alpha$	$\beta$	$\delta$	${d_1}$	$\mu$	$a$	${d_2}$	$\theta$
Value	$9$	$0.003$	$0.04$	$0.2$	$0.001$	$0.02$	$0.004$	$0.002$	$0.003$

With the above data, the value of the discriminant $\left( D \right)$ and coefficient $\,({C_1})$ of equation (6) are calculated as $1.9008\,$ and $- 1.392$, respectively, which give $\gamma _1^* = 0.00665248$ and $\gamma _2^* = 1.385335$. Then, the forward and backward transcritical bifurcations are shown by figures 2 and 3 respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, the conditions for the existence of an EE are derived and its stability behavior is discussed.

3.3. Investigation of the presence and local stability of the EE

This subsection investigates and discusses the local stability behavior of an equilibrium point of the system (3) for which the disease is endemic. By setting all equations of the system (3) to rest, i.e.,

$$\begin{align} p - \mu {S^*} - \frac{{\alpha {S^*}{I^*}}}{{1 + \beta {I^*}}} = 0, \nonumber \\ \frac{{\alpha {S^*}{I^*}}}{{1 + \beta {I^*}}} - \frac{{\gamma {I^*}}}{{1 + \delta {I^*}}} - \left( {\mu + {d_1} + \theta } \right){I^*} = 0, \nonumber \\ \frac{{\gamma {I^*}}}{{1 + \delta {I^*}}} - aQ_H^* - \left( {\mu + {d_2}} \right)Q_H^* = 0. \nonumber \\ \end{align} \tag{ 7 }$$

We get,

$$\begin{align}{S^*} = \frac{{p\left( {1 + {I^*}\beta } \right)}}{{{I^*}\alpha + \mu + {I^*}\beta \mu }},\,Q_H^* = \frac{{{I^*}\gamma }}{{\left( {1 + {I^*}\delta } \right)\left( {a + {d_2} + \mu } \right)}}, \text{ and } {I^*},\end{align}$$

where ${I^*}$ is a solution to the following equation:

$$\begin{equation}{A_2}{I^*}^2 + {A_1}{I^*} + {A_0} = 0\end{equation} \tag{ 8 }$$

and,

$$\begin{equation*}{A_0} = \mu \left( {{d_1} + \gamma + \theta + \mu } \right) - p\alpha = \mu \left( {{d_1} + \gamma + \theta + \mu } \right)\left( {1 - {R_0}} \right),\end{equation*}$$

$$\begin{equation*}{A_1} = \alpha \left( {\gamma - p\delta + \theta + \mu } \right) + {d_1}\left( {\alpha + \left( {\beta + \delta } \right)\mu } \right) + \mu \left( {\delta \left( {\theta + \mu } \right) + \beta \left( {\gamma + \theta + \mu } \right)} \right),\end{equation*}$$

$$\begin{equation*}{A_2} = \delta \left( {{d_1} + \theta + \mu } \right)\left( {\alpha + \beta \mu } \right).\end{equation*}$$

Now, for ${R_0} > 1$, ${A_0}$ is always negative and ${A_2}$ is positive. Thus, by Descartes's rule of signs [29], equation (8) has a unique positive value of ${I^*}$ for any ${A_1}$. Hence, the following result is presented as a theorem for the existence of the unique positive EE:

Theorem 3. The system (3) always has a unique positive equilibrium $\left( {{E_2} = \left( {{S^*},{I^*},Q_H^*} \right)} \right)$ whenever ${R_0} > 1$.

Now, we discuss the local stability behavior of ${E_2}$.

For the local stability behavior of ${E_2}$, first, we obtain the corresponding linearized matrix $\left( {{M_*}} \right)$ of the system (3) at ${E_2} = \left( {{S^*},{I^*},Q_H^*} \right)$ as given below:

$$\begin{equation*}{M_ * } = \left( \begin{array}{*{20}{ccc}} { - \frac{{{I^ * }\alpha }}{{1 + {I^ * }\beta }} - \mu }&{\frac{{{S^ * }{I^ * }\alpha \beta }}{{{{\left( {1 + {I^*}\beta } \right)}^2}}} - \frac{{{S^ * }\alpha }}{{1 + {I^*}\beta }}}&0 \\ {\frac{{{I^ * }\alpha }}{{1 + {I^ * }\beta }}}&{ - {d_1} - \frac{{{S^ * }{I^ * }\alpha \beta }}{{{{\left( {1 + {I^*}\beta } \right)}^2}}} + \frac{{{S^ * }\alpha }}{{1 + {I^*}\beta }} + \frac{{{I^ * }\gamma \delta }}{{{{\left( {1 + {I^ * }\delta } \right)}^2}}}\frac{\gamma }{{1 + {I^ * }\delta }} - \theta - \mu }&0 \\ 0&{\frac{{{I^ * }\gamma \delta }}{{{{\left( {1 + {I^ * }\delta } \right)}^2}}} + \frac{\gamma }{{1 + {I^ * }\delta }}}&{ - a - {d_2} - \mu } \end{array} \right).\end{equation*}$$

The characteristic equation of the matrix ${M_*}$ is

$$\begin{equation}\left( {\lambda + a + {d_2} + \mu } \right)\left( {{B_2}{\lambda ^2} + {B_1}\lambda + {B_0}} \right) = 0\end{equation} \tag{ 9 }$$

where

$$\begin{equation*}{B_2} = {\left( {1 + {I^*}\beta } \right)^2}{\left( {1 + {I^*}\delta } \right)^2},\,\,{B_1} = {b_1} - {b_2},\,\,{B_0} = {b_3} - {b_4},\end{equation*}$$

$$\begin{align*}{b_1} & = {d_1}{\left( {1 + {I^*}\beta } \right)^2}{\left( {1 + {I^*}\delta } \right)^2} + \left( {1 + {I^*}\beta } \right)\left( \gamma + \theta + 2\mu + {I^*}\left( \alpha {{\left( {1 + {I^*}\delta } \right)}^2} \right. \right. \nonumber\\ & \quad \left. \left. + \delta \left( {2 + {I^*}\delta } \right)\left( {\theta + 2\mu } \right) + \beta \left( {\gamma + {{\left( {1 + {I^*}\delta } \right)}^2}\left( {\theta + 2\mu } \right)} \right) \right) \right),\end{align*}$$

$$\begin{equation*}{b_2} = {S^*}\alpha {\left( {1 + {I^*}\delta } \right)^2},\end{equation*}$$

$$\begin{align*}{b_3} & = {I^*}\alpha \left( {1 + {I^*}\beta } \right)\left( {\gamma + {{\left( {1 + {I^*}\delta } \right)}^2}\theta } \right) + \left( {1 + {I^*}\beta } \right)\left( \gamma + \theta + {I^*}\left( \alpha {{\left( {1 + {I^*}\delta } \right)}^2} \right. \right. \nonumber\\ & \left. \left.\,\,\,\,\,\,\, + \left( {1 + {I^*}\beta } \right) + \delta \left( {2 + {I^*}\delta } \right)\theta + \beta \left( {\gamma + \theta } \right) \right) \right)\mu + {\left( {1 + {I^*}\beta } \right)^2}{\left( {1 + {I^*}\delta } \right)^2}{\mu ^2} \nonumber\\ & \quad + {d_1}\left( {1 + {I^*}\beta } \right){\left( {1 + {I^*}\delta } \right)^2}\left( {\mu + {I^*}\left( {\alpha + \beta \mu } \right)} \right),\end{align*}$$

$$\begin{equation*}{b_4} = {S^*}\alpha {\left( {1 + {I^*}\delta } \right)^2}\mu .\end{equation*}$$

From equation (9), it is clear that one root is negative, i.e. ${\lambda _1} = - \left( {a + {d_2} + \mu } \right)$ and the sign of other roots can be observed from the following quadratic equation:

$$\begin{equation}{B_2}{\lambda ^2} + {B_1}\lambda + {B_0} = 0.\end{equation} \tag{ 10 }$$

In equation (10), ${B_2}$ is always positive. Now, using the Routh–Hurwitz criterion, the following results are stated regarding the local stability behavior of ${E_2}$ in the form of two theorems.

Theorem 4. The endemic equilibrium ${E_2}$ of the system (3) exhibits local asymptotical stability under the simultaneous conditions:

• (i)  
  ${b_1} - {b_2} > 0$,
• (ii)  
  ${b_3} - {b_4} > 0.$

Theorem 5. The endemic equilibrium ${E_2}$ of the system (3) becomes unstable if any of the following conditions occur:

• (i)  
  ${b_1} - {b_2} < 0$;
• (ii)  
  ${b_3} - {b_4} < 0$;
• (iii)  
  ${b_1} - {b_2} < 0$ and ${b_3} - {b_4} < 0.$

In this section, we focus on comprehending the global stability behavior of the EE. To achieve this goal, we assume that:

$$\begin{equation*}F\left( I \right) = \frac{{\alpha I}}{{1 + \beta I}},\,G\left( I \right) = \frac{{\gamma I}}{{1 + \delta I}}.\end{equation*}$$

We now proceed to state and prove the following theorem regarding the global stability behavior of the EE:

Theorem 6. The endemic equilibrium $\left( {{E_2}\, = \,\left( {{S^*},\,{I^*},\,Q_H^*} \right)} \right)$ is globally asymptotically stable if ${R_0} > 1$.

Proof. At ${E_2}\, = \,\left( {{S^*},\,{I^*},\,Q_H^*} \right)$, the following conditions are derived from the system (3):

$$\begin{equation*}p = \mu {S^*} + {S^*}F\left( {{I^*}} \right);\,\left( {\mu + {d_1} + \theta } \right) = \frac{{{S^*}F\left( {{I^*}} \right) - G\left( {{I^*}} \right)}}{{{I^*}}};\,\left( {a + \mu + {d_2}} \right)\, = \frac{{G\left( {{I^*}} \right)}}{{Q_H^*}}.\end{equation*}$$

Now, to establish the global stability of $\,{E_2}$, we define the Lyapunov function:

$$\begin{equation*}{X_*}\left( t \right) = \left( {S - {S^*} - {S^*}\ln \frac{S}{{{S^*}}}} \right) + \left( {I - {I^*} - {I^*}\ln \frac{I}{{{I^*}}}} \right) + \left( {{Q_H} - Q_H^* - Q_H^*\ln \frac{{{Q_H}}}{{Q_H^*}}} \right).\end{equation*}$$

Differentiation of ${X_*}\left( t \right)\,$ along the solution of the system (3) is

$$\begin{align*}\frac{{d{X_*}\left( t \right)}}{{dt}} & = \left( {1 - \frac{{{S^*}}}{S}} \right)\frac{{dS}}{{dt}} + \left( {1 - \frac{{{I^*}}}{I}} \right)\frac{{dI}}{{dt}} + \left( {1 - \frac{{Q_H^*}}{{{Q_H}}}} \right)\frac{{d{Q_H}}}{{dt}} \nonumber\\ & = \left( {1 - \frac{{{S^*}}}{S}} \right)\left( {\mu {S^*} + {S^*}F\left( {{I^*}} \right) - \mu S - SF\left( I \right)} \right) + \left( {1 - \frac{{{I^*}}}{I}} \right)\left( SF\left( I \right) - G\left( I \right) \right. \nonumber\\ & \quad \left. - \left( {\frac{{{S^*}F\left( {{I^*}} \right) - G\left( {{I^*}} \right)}}{{{I^*}}}} \right)I \right) + \left( {1 - \frac{{Q_H^*}}{{{Q_H}}}} \right)\left( {G\left( I \right) - \frac{{G\left( {{I^*}} \right){Q_H}}}{{Q_H^*}}} \right) \nonumber\\ & = \left( {1 - \frac{{{S^*}}}{S}} \right)\left( {\mu {S^*}\left( {1 - \frac{S}{{{S^*}}}} \right) + {S^*}F\left( {{I^*}} \right)\left( {1 - \frac{{SF\left( I \right)}}{{{S^*}F\left( {{I^*}} \right)}}} \right)} \right) \nonumber\\ & \quad + \left( {1 - \frac{{{I^*}}}{I}} \right)\left( {{S^*}F\left( {{I^*}} \right)\left( {\frac{{SF\left( I \right)}}{{{S^*}F\left( {{I^*}} \right)}} - \frac{I}{{{I^*}}}} \right) + G\left( {{I^*}} \right)\left( {\frac{I}{{{I^*}}} - \frac{{G\left( I \right)}}{{G\left( {{I^*}} \right)}}} \right)} \right) \nonumber\\ & \quad + \left( {1 - \frac{{Q_H^*}}{{{Q_H}}}} \right)\left( {\frac{{G\left( I \right)}}{{G\left( {{I^*}} \right)}} - \frac{{{Q_H}}}{{Q_H^*}}} \right)G\left( {{I^*}} \right).\end{align*}$$

After some arrangements, we have

$$\begin{align*}\frac{{d{X_*}}}{{dt}} & = \mu {S^*}\left( {2 - \frac{S}{{{S^*}}} - \frac{{{S^*}}}{S}} \right) + {S^*}F\left( {{I^*}} \right)\left( {2 - \frac{I}{{{I^*}}} - \frac{{{S^*}}}{S} - \frac{{F\left( I \right)}}{{F\left( {{I^*}} \right)}}\left( {\frac{{{I^*}S}}{{{S^*}I}} - 1} \right)} \right) \nonumber\\ & \quad + G\left( {{I^*}} \right)\left( {\frac{I}{{{I^*}}} - \frac{{{Q_H}}}{{Q_H^*}} + \frac{{G\left( I \right){I^*}}}{{G\left( {{I^*}} \right)I}} - \frac{{G\left( I \right)Q_H^*}}{{G\left( {{I^*}} \right){Q_H}}}} \right).\end{align*}$$

Now, from the property of arithmetic mean, we have

$$\begin{equation*}\left( {2 - \frac{S}{{{S^*}}} - \frac{{{S^*}}}{S}} \right) \unicode{x2A7D} 0,\end{equation*}$$

and if

$$\begin{equation*}\left( {2 - \frac{I}{{{I^*}}} - \frac{{{S^*}}}{S} - \frac{{F\left( I \right)}}{{F\left( {{I^*}} \right)}}\left( {\frac{{{I^*}S}}{{{S^*}I}} - 1} \right)} \right) \unicode{x2A7D} 0,\end{equation*}$$

$$\begin{equation*}\left( {\frac{I}{{{I^*}}} - \frac{{{Q_H}}}{{Q_H^*}} + \frac{{G\left( I \right){I^*}}}{{G\left( {{I^*}} \right)I}} - \frac{{G\left( I \right)Q_H^*}}{{G\left( {{I^*}} \right){Q_H}}}} \right) \unicode{x2A7D} 0,\end{equation*}$$

then, $\frac{{d{X_*}}}{{dt}} \unicode{x2A7D} 0$. Hence, Using LaSalle's Invariance Principle [30], it can be deduced that if ${R_0} > 1$, the singleton set $\left\{ {{E_2}} \right\}$ is the largest invariant set in $\{ \left( {S,I,{Q_H}} \right)|\frac{{d{X_*}}}{{dt}}\, = \,0\}$ and is globally asymptotically stable.

Sensitivity analysis of ${R_0}$ helps to understand how uncertain the calculated value of ${R_0}$ is, and how much the outcome of the disease spread can be affected by parameter variations. To check the ${R_0}$ sensitivity, we calculate its derivatives as follows:

$$\begin{equation*}\frac{{\partial {R_0}}}{{\partial p}} = \frac{\alpha }{{\mu \left( {{d_1} + \gamma + \theta + \mu } \right)}}.\end{equation*}$$

$$\begin{equation*}\frac{{\partial {R_0}}}{{\partial \alpha }} = \frac{p}{{\mu \left( {{d_1} + \gamma + \theta + \mu } \right)}}.\end{equation*}$$

$$\begin{equation*}\frac{{\partial {R_0}}}{{\partial \mu }} = - \frac{{p\alpha \left( {{d_1} + \gamma + \theta + 2\mu } \right)}}{{{\mu ^2}{{\left( {{d_1} + \gamma + \theta + \mu } \right)}^2}}}.\end{equation*}$$

$$\begin{equation*}\frac{{\partial {R_0}}}{{\partial {d_1}}} = - \frac{{p\alpha }}{{\mu {{\left( {{d_1} + \gamma + \theta + \mu } \right)}^2}}}.\end{equation*}$$

$$\begin{equation*}\frac{{\partial {R_0}}}{{\partial \gamma }} = - \frac{{p\alpha }}{{\mu {{\left( {{d_1} + \gamma + \theta + \mu } \right)}^2}}}.\end{equation*}$$

$$\begin{equation*}\frac{{\partial {R_0}}}{{\partial \theta }} = - \frac{{p\alpha }}{{\mu {{\left( {{d_1} + \gamma + \theta + \mu } \right)}^2}}}.\end{equation*}$$

Because all parameters are positive, therefore, $\frac{{\partial {R_0}}}{{\partial p}} > 0,\,\frac{{\partial {R_0}}}{{\partial \alpha }} > 0,\,$ and $\frac{{\partial {R_0}}}{{\partial \mu }} < 0,\,$ $\frac{{\partial {R_0}}}{{\partial {d_1}}} < 0,\,\frac{{\partial {R_0}}}{{\partial \gamma }} < 0,\,\frac{{\partial {R_0}}}{{\partial \theta }} < 0$. Thus, ${R_0}$ is increasing with p, $\,\alpha$ and decreasing with $\mu ,\,{d_1},\,\gamma ,\theta$.

In this section, system (3) is numerically simulated using Mathematica 11. The simulations are conducted using the numerically experimental parameter values listed in table 3.

Table 3. The parameter's numerical value.

Parameter	$p$	$\alpha$	$\beta$	$\delta$	${d_1}$	$\mu$	$a$	${d_2}$	$\theta$	$\gamma$
Value	$4$	$0.003$	$0.05$	$0.0002$	$0.01$	$0.02$	$0.002$	$0.02$	$0.003$	$0.03$

The initial subpopulations are:

$$\begin{equation*}S\left( 0 \right) = 190,\,\,I\left( 0 \right) = 9,\,\,{Q_H}\left( 0 \right) = 1.\end{equation*}$$

Using the dataset given in table 3, the EE ${E_2}\left( {S,I,{Q_H}} \right)$ and the BRN $\,{R_0}$ are calculated as $\left( {65.6863,42.8123,30.3206} \right)$ and $9.52$, respectively. Figure 4 illustrates the dynamics of the susceptible $\left( S \right)$, infected $\left( I \right)$, and hospitalized $\left( {{Q_H}} \right)$ subpopulations when ${R_0} > 1$. It is evident from the figure that all subpopulations reach their endemic steady state, ${E_2}\,\left( {65.6863,42.8123,30.3206} \right)$, over time. This observation confirms the local asymptotic stability result of the EE ${E_2}$.

Zoom In Zoom Out Reset image size

Figure 5 illustrates the number of infected individuals at various initial values of infected individuals $\left( {I\left( 0 \right)} \right)$. The figure demonstrates that, regardless of the initial number of infected individuals, the infected population eventually reaches the same steady state. This observation confirms the global stability result of our study.

Zoom In Zoom Out Reset image size

The impact of the increased transmission rate $\left( \alpha \right)$ and inhibition rate $\left( \beta \right)$ on the infected population is illustrated in figures 6 and 7, respectively. As shown in figure 6, when the transmission rate rises, the number of infected individuals also increases. Conversely, figure 7 illustrates that as the inhibition rate increases, the number of infected individuals decreases. These findings suggest that controlling the spread of disease in society can be achieved by focusing on reducing the rate of infection transmission within the community and encouraging infected individuals to take necessary precautions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The impact of an increased hospitalization rate $\left( \gamma \right)$ and a limitation rate in hospitalization $\left( \delta \right)$ on the infected population is illustrated in figures 8 and 9, respectively. As depicted in figure 8, when the hospitalization rate increases, the number of infected individuals decreases. Conversely, figure 9 shows that as the limitation rate in hospitalization increases, the number of infected individuals also increases. These findings suggest that by enhancing the capacity to treat patients through measures such as adding hospital beds and increasing healthcare workers, doctors, and medical equipment, more infected individuals can receive treatment. This, in turn, reduces the risk of infection and lowers the number of fatalities in society.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 10 and 11 demonstrate the impact of an increased hospitalization rate $\left( \gamma \right)\,$ and the limitation rate in hospitalization $\left( \delta \right)$ on the number of hospitalized individuals, respectively. As depicted in figure 10, as the value of $\gamma$ increases, the number of infected individuals who are severely ill also rises. Conversely, figure 11 illustrates that by minimizing $\delta$, more severely ill infected individuals can receive treatment in hospitals, which consequently leads to a notable improvement in the recovery of the infected population.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 illustrates the importance of incorporating the hospitalized population class into an epidemic model, especially concerning the number of individuals who have recovered from infection. The figure demonstrates that when the recovery of infected individuals through hospital treatment is considered alongside natural recovery, simultaneously, it results in a significantly greater increase in the total recovered population compared to considering a single recovery process such as relying solely on medical treatment or the auto-immune responses.

Zoom In Zoom Out Reset image size

The goal of this study is to propose an epidemic model that includes a hospitalized class and analyze the stability behavior of the model. The incidence rate and hospitalization rate are modeled using the Holling type II functional form. Additionally, the BRN $\left( {{R_0}} \right)\,$ is calculated using the next-generation matrix method. The paper concludes that the DFE is locally asymptotically stable if ${R_0} < 1$. This implies that the disease can be eradicated by maintaining the threshold parameter, ${R_0}$, below unity. The center manifold theory is employed to analyze the model's behavior at ${R_0} = 1$. The analysis reveals that a forward (or backward) transcritical bifurcation occurs under specific conditions. Specifically, a forward transcritical bifurcation occurs when ${R_0}$ crosses unity from below, resulting in the emergence of a small, positive, asymptotically stable equilibrium. Conversely, a backward transcritical bifurcation occurs when ${R_0} < 1$. In this case, a small, positive unstable equilibrium appears alongside the locally asymptotically stable DFE and a stable larger positive equilibrium. From an epidemiological standpoint, a backward bifurcation implies that reducing the BRN to less than unity is insufficient to eliminate the disease. The paper also derives the conditions for the existence of the EE and demonstrates that the EE is locally asymptotically stable under certain conditions as stated in theorem (4), but unstable under others as given in theorem (5). Furthermore, the paper investigates the global stability of the EE and determines its global asymptotic stability under specific conditions when ${R_0} > 1$.

The theoretical findings of this study are strongly validated through numerical simulations. Figures 2 and 3 effectively illustrate the range and occurrence of forward and backward transcritical bifurcations in the model, utilizing the data from table 2. These simulations provide robust evidence, numerically validating both theorems 1 and 2. Figure 4 confirms the local asymptotical stability of the EE, as all subpopulations reach their respective steady states over time. Figure 5 is included to confirm the global stability result of the study. The evidence from figure 5 is clear: irrespective of the initial infected population, the system converges to the same equilibrium point. Moreover, the simulations demonstrate that adopting effective measures, such as wearing masks, regular hand sanitization, glove usage, hygiene maintenance, etc, by the population can significantly control the spread of infection. The importance of considering inhibitory effects is highlighted in figure 7. Furthermore, advocating for improvements in medical infrastructure allows severely infected individuals to access specialized treatment, resulting in more people receiving proper care, as depicted in figures 8 and 9. The significance of incorporating the 'Hospitalized' individuals class and considering a saturated hospitalization rate is demonstrated in figures 10 and 11, respectively. The importance of incorporation of hospitalized individual's class in the context of the increasing recovered population is also depicted in figure 12.

Beyond stability, the conclusion is that governments and health professionals should swiftly develop programs such as contact tracing and quarantine, while also ensuring an adequate number of spaces for quarantine. Additionally, ensuring the availability of medical infrastructure and adequate treatment facilities, including more hospitals, ventilators, and ICU beds, can help reduce the spread and accelerate the recovery of the infected population. These methods work to limit virus transmission and provide necessary care for the infected. Moreover, such programs should be sufficiently long-lasting to ensure the complete cessation of infection emergence. By increasing health education among individuals, developing a robust health infrastructure, and improving hospital facilities, we can effectively reduce the spread of infection and consequently decrease the loss of human lives. Health education plays a crucial role in increasing awareness and understanding of the virus and its transmission, leading to better compliance with measures such as contact tracing, quarantine, etc. This, in turn, helps alleviate the strain on healthcare systems during any epidemic outbreak.

The authors sincerely appreciate and thank the editors and anonymous referees for their invaluable comments and suggestions on the manuscript. The authors highly value their insightful feedback and helpful input, which significantly improved the manuscript's overall quality.

All data that support the findings of this study are included within the article (and any supplementary files).

This work was not funded by any specific source.

The authors confirm no conflicts of interest.

All authors contributed equally to conceptualization, formal analysis, writing the original draft, and reviewing/editing of the manuscript.