Global dynamics and control strategies of an epidemic model having logistic growth, non-monotone incidence with the impact of limited hospital beds.

In this paper, we have considered a deterministic epidemic model with logistic growth rate of the susceptible population, non-monotone incidence rate, nonlinear treatment function with impact of limited hospital beds and performed control strategies. The existence and stability of equilibria as well as persistence and extinction of the infection have been studied here. We have investigated different types of bifurcations, namely Transcritical bifurcation, Backward bifurcation, Saddle-node bifurcation and Hopf bifurcation, at different equilibrium points under some parametric restrictions. Numerical simulation for each of the above-defined bifurcations shows the complex dynamical phenomenon of the infectious disease. Furthermore, optimal control strategies are performed using Pontryagin's maximum principle and strategies of controls are studied for two infectious diseases. Lastly using efficiency analysis we have found the effective control strategies for both cases.

Introduction

The main objective of mathematical epidemiology is to study of the mechanism that obsesses the disease transmission and controls its impact in the population [1, 2]. Mathematical modeling becomes an important tool to investigate the dynamical evolution of infectious disease and effective measures as needed to reduce the influence of the disease [3-5]. Using mathematical models, researchers can identify trends of the disease, analyze epidemiological surveys and make general forecasts about the disease. First classical model for smallpox spreading was designed and analyzed by Bernoulli in 1760 [6]. Thereafter, many researchers have used various types of mathematical models to express the dynamical changes in the population due to influence of the disease and to prevent its impact. Among them, in 1927 Kermack and Mckendrick’s compartmental SIR model is a milestone for deterministic process in theory of epidemiology [7].

Population birth rate and disease transmission rate play important role in behavioral change of disease dynamics. Two types of birth rates, namely constant birth rate [8-10] and logistic growth rate [11], are used to formulate compartmental model. Constant birth rate is not more realistic for large populated countries. To prevent huge number of spreading infection, logistic type of birth rate is considered in mathematical modeling for limited medical resources, long-lasting disease, disease with high death rate and large number of the population [3, 5]. Various types of disease incidence rates are considered in mathematical modeling to describe the disease spreading. Bilinear disease transmission rate () is used based on mass action property where is the rate of disease transmission per contact and x and y represent density of susceptible and infected individual, respectively [12]. For largely populated countries, consideration of bilinear incidence rate is not realistic as rate of disease transmission increases with the increase in susceptible population. So for largely populated countries standard form of incidence is used in the form where N is total number of population [13, 14]. For standard incidence rate, disease transmission rate does not tend to infinity with number of susceptible increases. In case of bilinear and standard incidence rate, we get simple dynamical behavior of the model system. To get rich dynamics, nonlinear incidence rates are used in mathematical modeling. Yorke and London [15] used the incidence rate in the form to describe measles outbreak and also Levin and Iwasa [16] considered in the form . For concerning inhibition effect of susceptible population and crowding effect of infected population, saturated type of incidence rate was first used by Capasso and Serio [17] in the form where is the force of infection and is inhibition effect due to change in behavior of susceptible population and also crowding effect of the infected individuals. At first, disease transmission increases with the increase in infected population, and for huge number of infected individuals it ultimately tends to which indicates for social awareness of susceptible population, disease transmission gradually decreases [18, 19]. Later on for describing disease dynamics of some short-term but highly infectious diseases such as SARS, influenza Xiao and Ruan [20] proposed non-monotone type of incidence rate in the form . Here, at the beginning due to lack of knowledge about the infection, the disease spreads among the people in short time. After some days, when the disease is spreading highly among the people, then due to social awareness and other psychological factors people take necessary actions to control the infection, and then, the rate of infection will start to decrease. The above-defined incidence rate increases first with y and reaches the maximum value when then starts to decrease and tends to zero as . Further, Ruan and Wang [21], Lu and Ruan [8] studied different types of nonlinear incidence rates, namely , , respectively, to get more complicated dynamical phenomenon of the disease such as periodic oscillations, multiple peaks of the infection, different types of bifurcations. Generalized type of nonlinear transmission rate where 0 and is suppression effect co-efficient is investigated by Liu [22], Hethcote and Driessche [23] to observe the effect of behavioral changes.

To get recovery from infectious disease, treatment is an important method. Treatment has effective role to prevent and control the spreading of the disease. Several researchers studied different types of treatment functions such as vaccination, quarantine, isolation and medicine. Proper and timely treatment can reduce the effect of the spreading disease. In classical epidemic model, rate of treatment is assumed as proportional to the number of infected individuals. But this type of treatment function is not realistic for large populated countries since treatment is generally based on medical resources such as medicines, doctors, hospital beds, vaccines and isolation places. Every country has a limited capacity in medical resources, and due to this limitation there can be delayed in getting the treatment, so it is very important to adopt a suitable recovery function. Wang and Ruan [24] introduced a constant recovery function T(y) into epidemic model which is in the form = maximum medical resources is used here to cure the disease. Further, Wang [25] proposed a piecewise linear recovery function in the form = where recovery rate is proportional to number of infected population before reaching its maximum capacity. Later, Zhang and Liu [26] considered a saturated type recovery function in the form T(y) = where delay in receiving treatment and limited medical resources are included. Consideration of this type of saturated treatment function has an advantage that T(y) has linear type characteristic when the number of infected population is very low, whereas it tends to a constant limit for higher value of infected population. They established that the system experiences Backward bifurcation in the presence of saturated type treatment. Additionally, Zhou and Fan [27] modified this type of treatment function by introducing in the form and studied existence, stability of equilibria, Backward bifurcation and also Hopf bifurcation. Later on, obeying the report of World Health Organization (WHO) Statistical Information System, health planners used hospital bed population ratio (HBPR), number of hospital beds available per 10,000 population as a method for estimating medical resource ability to the public. Considering WHO’s rationale, Shan and Zhu [28] introduced removal rate as a function of hospital beds (b) and number of infected individuals y in the formwhere is the minimum per capita recovery rate and is the maximum per capita recovery rate. They have established that the system with standard incidence and limited bed recovery rate shows rich and interesting dynamical phenomenon such as Saddle-node bifurcation, Backward bifurcation, Hopf bifurcation and Bogdanov–Takens bifurcation of co-dimension 2 and 3. In [29], Abdelrazec studied recovery rate to investigate available medical resources as spreading and controlling dengue fever which is helpful for health planners to allocate resources for controlling dengue transmission. Recently, the limitation of hospital beds is highly applicable for providing treatment to the COVID-19-infected population.

On the other hand, use of optimal control strategy in epidemic model is an important mathematical tool. Using optimal control in epidemic model the health planners can find out the best strategy which is more effective to reduce the disease spreading among the population with minimum cost [30]. Sometimes, only vaccination is applied to control the disease [31, 32], whereas only treatment is also used to reduce the influence of the infection [33]. Both vaccination and treatment are used in [34] to minimize the impact of the infection as well as reduce the cost of implementation due to applying treatment and vaccination to the individuals. Our intention will be fulfilled if a number of infective individuals are reduced by applying optimal control strategies taking the best effective preventive measures for controlling the disease.

The proposed model in this paper is extension of the work of Xiao and Ruan [20], Shan and Zhu [28]. We have extended the model considering the nonlinear incidence rate = and treatment function as = + in the presence of logistic type birth rate for susceptible population. We have studied different types of bifurcations such as Transcritical, Backward, Saddle-node and Hopf bifurcations and finally investigate the optimal control policy and efficiency analysis of the applied controls.

The rest of the paper is organized as follows : In Sect. 2, we have formulated the proposed model and discussed the basic properties of the solutions. Then, we have analyzed different cases for existence of equilibria and obtained basic reproduction number in Sect. 3. Next in Sect. 4, we have examined the necessary and sufficient conditions for local stability of disease-free equilibrium and endemic equilibrium point. In Sect. 5, we have studied persistence and extinction of the infectious disease under some suitable conditions. We have investigated bifurcations such as Transcritical bifurcation, Backward bifurcation, Saddle-node bifurcation and Hopf bifurcation in Sect. 6. In Sect. 7, we have formulated optimal control of the problem and used Pontryagin’s maximum principle for obtaining necessary conditions for optimality. Finally, in Sect. 8 we have summarized the results and given epidemiological significance for number of limited hospital beds and interpreted strategies for controlling the emerging disease.

Model formulation and boundedness

We assume that the total population is divided into three classes, namely susceptible, infected and recovered classes. Let, x(t), y(t) and z(t) represent density of the susceptible, infected and recovered population at any time t, respectively. Let the susceptible population have logistic growth rate with two types of death rates, namely normal death and disease-induced death. The rate of incidence is non-monotonic type in the form , and the treatment rate is of the form . Based on the above discussion, we have proposed the deterministic epidemic model in the following form:satisfying initial conditions , , for all . The model parameters of the system (1) are enlisted in Table 1.

Table 1

Model parameters and their descriptions

Parameter	Description
r	Intrinsic growth rate of susceptible individuals
K	Carrying capacity of the population
β	Disease transmission rate
α	Parameter measuring inhibitory factors
d	Natural death rate of population
δ	Disease-induced death rate
μ	Natural recovery rate of infected individuals
u_1	Control parameter denoting vaccination
b	Number of available hospital beds
μ_0	Minimum number of per capita recovery
μ_1	Maximum number of per capita recovery
R_0	Basic reproduction number

In the proposed epidemic model, we assume that the total population remains constant or approaches asymptotically a constant which is reasonable for quickly spreading disease or rarely death-caused disease [28]. Since the third equation of (1) is independent from the other two, for stability analysis purpose we omit the last equation, but to investigate the optimal control problem we will consider all the three equations.with initial conditions , for all . In the next sub-section, we shall investigate the basic properties such as positivity and boundedness of the solutions of the reduced model system (2).

Basic properties of the model

First, we shall show that all solutions of model system (2) satisfying the initial conditions , are non-negative for all , and then, we shall prove that the solutions are bounded for all .

Theorem 1

All solutions (x(t), y(t)) of model system (2) satisfying the non-negative initial conditions are non-negative for any value of .

Proof

First, we show x(t) is non-negative for all values of t. To prove this, we proceed with first equation of model system (2)where M = Integrating and simplifying, we getThus, we have established , i.e., the susceptible populations are always positive. Similarly, we can show that .

This completes proof of the theorem.

Next we show solutions of model system (2) are uniformly bounded for all . We establish our result in the following theorem :

Theorem 2

The set is a positively invariant attracting region for the epidemic model given by the system (2) in .

Adding both equations of model system (2), we geti.e.,Integrating and taking limit as , we getThus, with . Using this result and the result of Theorem 1, we can state that the solution of the model (2) is positively bounded in the region .

Existence of equilibria and basic reproduction number

In this section, we shall discuss the existence of different equilibrium points and basic reproduction number () of model system (2). Basic reproduction number is the average number of secondary infections produced by an infective during an entire period of infection. The model system (2) has always a trivial equilibrium point and a disease-free equilibrium point = will exist when . Hence, can be calculated using next-generation matrix approach [35]. In the following theorem, we compute basic reproduction number using Driessche and Watmough’s [35] technique.

Theorem 3

Basic reproduction number of the model system (2) is = .

Proof

The model system (2) has only one infected component, namely y. Let and be the new infection terms and remaining terms, respectively, then and Let and , then the next-generation matrix is = = . Therefore, basic reproduction number of the model system (2) is given by

The components of endemic equilibrium point of the system (2) are = and is positive root of p(y) = 0and .

Since is a fifth-order equation, it is very difficult to find the expressions of positive roots of the equation explicitly in terms of . Using the Descartes’ rules of signs, we can get idea about the number of positive roots of Eq. (3) which is summarized in Table 2 in terms of and the thresholds and .

Table 2

Number of positive roots for different values of

Considering and keeping all other parameters fixed as shown in Table 3, we are varying the disease transmission rate (consequently will vary) and have drawn the solution curve for infected population (i.e., for existence of endemic equilibria) with respect to (see Fig. 1). The number of equilibrium points for different values of is summarized in the following lemma:

Table 3

Model parameters and their respective values

Parameter	u_1	d	δ	μ	K	α	μ_0
Value	0.175	0.1	0.2	0.6	40	0.246	0.1

Fig. 1

Solution curve with respect to for a no root for , b two roots for , c one root for , d three roots for

Lemma 1

Model system (2) always has a trivial equilibrium point and a disease-free equilibrium point .

For , model system (2) may have no, one coincident or two endemic equilibrium points.

For , model system (2) may have one, two with one coincident or three endemic equilibrium points.

Local stability analysis

To investigate the nature of the equilibrium points, we first need to compute Jacobian matrix of model system (2) at any equilibrium point E(x, y) which is given by

Theorem 4

Trivial equilibrium point is stable if and saddle point if .

The Jacobian matrix (4) at the trivial equilibrium point is:Eigenvalues of the above Jacobian matrix are given by and . Clearly, one eigenvalue is negative and other will be negative if , and hence, under this restriction is stable. If holds, then one eigenvalue is positive and other value is negative; consequently, is a saddle point.

In biological point of view, means birth rate of susceptible is less than sum of vaccination rate and death rate of the susceptible population. Two cases may arise: (i) high vaccination (i.e., is nearer to r such that ) and (ii) high death rate. The first case implies most of the people are vaccinated to get recovery from a particular disease. As for example for preventing polio, tuberculosis disease, etc., vaccination is given step by step from the childhood. As a result, these types of diseases are eliminated now from the population. The second case is totally impossible because in this case all susceptible population will die out in the community. On the other hand, means birth rate of the population is higher than sum of vaccination and death rate of population which is realistic.

Theorem 5

If , then the disease-free equilibrium point is locally asymptotically stable, and for it is unstable. For , disease-free equilibrium point is a saddle-node of co-dimension 1 if and a semi-hyperbolic attracting node of co-dimension 2 if .

Jacobian matrix of the system (2) at the disease-free equilibrium point isEigenvalues of are and . If , then both the eigenvalues are negative, and hence, is stable and is unstable(saddle) for .

For , one of the above eigenvalues vanishes, and hence, the disease-free equilibrium point is a non-hyperbolic equilibrium point. Use center manifold theorem to determine the nature of equilibrium point .

First we translate disease-free equilibrium point to origin using the transformation , ; then, the model system (2) takes the following form: Introducing the transformation , the System (5) becomes : It is clear from (6) that if , then disease-free equilibrium point is a saddle-node of co-dimension 1 [36]. Now, if , then use center manifold theorem [36] for investigating nature of the disease-free equilibrium point .

Let the local center manifold of the system is given by satisfying where is sufficiently small and Dh is the derivative of h with respect to x. Using the conditions, we get , .

Thus, local manifold of the system is given by , and then, the system (6) reduces toFrom system (7), it is clear that disease-free equilibrium point is a semi-hyperbolic attracting node of co-dimension 2 if . Hence, the theorem is proved.

Next we shall establish the local stability of the endemic equilibrium in the following theorem:

Theorem 6

The characteristic endemic equilibrium is locally asymptotically stable if and is unstable otherwise.

Jacobian matrix of the model system (2) at the endemic equilibrium isThe corresponding characteristic equation iswhere and .

Roots of the above equation will have negative real part if and .

Now,   if and =   if   .

That means and if .

Hence, endemic equilibrium is locally asymptotically stable if and is unstable otherwise. Thus, we arrive at our desired result.

Persistence and extinction of the infection

Basic reproduction number is playing a crucial role in dynamical change of infection. In this section, we shall establish the eradication and persistence of disease is dependent on the value of . Eradication and persistence of the disease can be determined by stability of disease-free equilibrium point and the endemic equilibrium . First we shall show the disease-free equilibrium point is globally stable when less than a threshold value then infection will be eradicated from the system. In the following theorem, we shall prove the global stability of disease-free equilibrium point .

Theorem 7

Disease-free equilibrium point is globally stable if where .

Using Lyapunov’s stability theorem [37, 38], we shall establish the global stability of disease-free equilibrium point.

Let, V : R where defined by is a Lyapunov function as it is positive definite.

Then, .

i.e.,

i.e.,

i.e.,

i.e.,

Therefore, if where .

So, by Lyapunov’s Stability Theorem [37, 38] one can conclude that is globally stable if .

The above result is epidemiologically most important because if the value of is less than a threshold quantity, then the disease will be eliminated from the population. Therefore, to control the influence of the infection, our target is to adopt proper preventive policies such that value of is less than that threshold value.

Now we shall prove global stability of endemic equilibrium point . For discussing global stability of endemic equilibrium point, use Dulac criterion [25, 36]. In the following theorem, we shall establish the condition of global stability of the endemic equilibrium point.

Theorem 8

For , endemic equilibrium point is globally asymptotically stable if satisfied.

We write our model system (2) as given below:To show global stability of endemic equilibrium , we use Dulac function criterion [25, 36] of stability analysis. Considering Dulac function as , we haveAgain implies .

Combining the above two conditions, we get if   .

Thus, by Dulac criterion we can conclude that system (2) has no closed orbit that means endemic equilibrium point is globally stable for under the conditions stated in the theorem.

we have shown in the above theorem that for there exist a certain region of the model parameter where endemic equilibrium is globally asymptotically stable that means infection persists in the system. If new infection produced by single susceptible during its entire infectious period is greater than one, then there exist a certain number of infected individuals. Therefore, the number of infected individuals gradually diminishes if new infection produced by single susceptible during lifespan less than one. Thus, the necessary condition to eliminate or persist of infection is to maintain the value of basic reproduction number less or greater than one, respectively.

Bifurcation analysis and epidemiological interpretation

In this section, we shall discuss various types of bifurcations about different equilibrium points considering different model parameters as the bifurcation parameters and interpret the biological consequences of the results at the corresponding bifurcation values of the model parameters. We have studied Transcritical bifurcation with respect to maximum number of per capita recovery rate , and it gives the upper limit of above which disease will eliminate from the system and below which disease will persist. Again Backward bifurcation is studied with respect to . We have found the threshold number of available hospital beds (b) below which the system shows Backward bifurcation, i.e., disease eradication not only depends on values of but also on initial number of infection. If the value of intrinsic growth rate r passes through its critical value , then new endemic equilibrium points are created and destroyed through Saddle-node bifurcation. Thus, we get different stability criteria for controlling the spreading of the disease. The endemic equilibrium point may lose its stability through Hopf bifurcation when number of available hospital beds (b) passes through its critical value. Thus, we have shown the condition for persistence of the disease may change through Hopf bifurcation.

Transcritical bifurcation

Theorem 9

The model system (2) experiences Transcritical bifurcation at the disease-free equilibrium as the recovery parameter passes through the critical value if where .

Since at , i.e., , one of the eigenvalues of is zero and other is negative, to investigate the nature of the solution at , we have to use Sotomayor’s theorem [36, 39]. For this purpose, we construct a function with two components given by whereLet corresponding to zero eigenvalue and have eigenvectors V and W, respectively, where and From (10) we have , and =

Therefore,Hence, the model system (2) satisfies all conditions of Sotomayor’s theorem for Transcritical bifurcation with respect to parameter . Thus, model system (2) experiences Transcritical bifurcation in the neighborhood of the critical value of the maximum value of the recovery parameter .

Significance of as a Transcritical bifurcation parameter

Here, numerically we shall interpret significance of the maximum value of the per capita recovery rate ( ) as the Transcritical bifurcation parameter for lower and upper value of the rate of infection(). We consider the lower value as and the upper value as .

(a) For lower value of disease transmission rate

For and other parameter values are described as shown in Table-3 with , we get the critical value of Transcritical bifurcation parameter as . For , one saddle disease-free equilibrium point and one stable endemic equilibrium point  exist (see the corresponding phase portrait in Fig. 2a). for , there exist one stable disease-free equilibrium point and two endemic equilibrium points and where is saddle point and is stable spiral (see the phase portrait in Fig. 2b). Thus, when passes from left side to right side through its critical value , then the disease-free equilibrium point and the endemic equilibrium point exchange their stability through Transcritical bifurcation with creation of one stable endemic equilibrium point. If we take value of far from the transcritical value in right side of say , then only a stable disease-free equilibrium point exists (see the phase portrait in Fig. 2c).

Fig. 2

Phase portrait at different values of in the neighborhood of the Transcritical bifurcation value for lower

It is clear from the above discussion that there exists a critical value of the maximum value of treatment rate below which the disease will persist in the system for any initial population density. But for the values of greater than the critical value persistence of the disease will depend on the population density of the initial spreading. In this situation, two stable equilibrium points exist: one is disease-free equilibrium and other is with disease. In this case, initial population density plays a crucial role for persistence or elimination of the disease. Again for higher value of the maximum value of cure rate, disease is eradicated from the system which indicates applying high level of treatment disease can be eradicated from the population. Therefore, population density and treatment both have vital role in disease controlling.

(b) For higher value of disease transmission rate

For and other parameter values are given in Table 3 with , we get the critical value of the maximum cure rate as . Now, for , only one stable disease-free equilibrium point exists, i.e., disease will not persist in the system (see the phase portrait in Fig. 3a) and for , there exist a unstable (saddle) disease-free equilibrium point and one unstable endemic equilibrium point with a stable limit cycle (see the phase portrait in Fig. 3b, the red curve is stable limit cycle). Taking the value of more in the left side of , say , there exist one unstable disease-free equilibrium point and one endemic equilibrium point which is stable spiral (see the phase portrait in Fig. 3c). Hence, for higher values of the disease transmission rate, the transcritical value of maximum cure rate is also high; then, the disease will be eradicated from the system for higher amount of treatment. For investigating disease dynamics, we observe the behavioral change of equilibrium point in left side as well as right side of .

Fig. 3

Phase portrait at different values of in the neighborhood of the Transcritical bifurcation value for higher

From the above, we have seen that for higher treatment rate disease can be eradicated from the population. In just left side of , i.e., for , disease will persist in the system with oscillatory population density, and for very lower value of , i.e., for disease will persist in the system with steady state. Thus, the treatment has important role in controlling disease transmission.

Therefore, it is clear from the above analysis that for lower rate of the infection if we take proper action, i.e., apply moderate or high medical facility, then disease can be easily eradicated from the system, but if the rate of the infection is higher, then to eliminate the disease, we need to apply high medical facility to control the disease. If we provide lower medical facility, then disease will persist in the system either it oscillates or in steady state.

Backward bifurcation

Theorem 10

The model system (2) experiences Backward bifurcation if when crosses the critical value and Forward bifurcation if

To determine the condition for Backward bifurcation of the model system (2), we use the methodology as developed in [40]. Let us consider as the bifurcation parameter and implies , and for this value of one of the eigenvalues of vanishes. Let and be the right eigen vector and left eigen vector, respectively, of corresponding to the eigenvalue . In order to apply Castillo–Chavez and Song’s bifurcation theorem [40], we need to compute the value of and which are given byThus, the model system (2) experiences Backward bifurcation at if , i.e., , and for bifurcation is forward.

In Fig. 4, we have plotted the bifurcation diagram with respect to varying the disease transmission rate . It is clear from the figure that there is a critical value of , say , such that in there exists only one stable disease-free equilibrium point, for , there exist two endemic equilibrium points along with the stable disease-free equilibrium point, and for one stable endemic equilibrium point with unstable disease-free equilibrium point exists. Among the two endemic equilibrium points, the equilibrium with lower infected level is saddle, but point with higher infected level is stable, which means for , one stable disease-free equilibrium point as well as one stable endemic equilibrium point exist. Epidemiologically this result is highly significant because eradication of disease from the system depends on values of as well as on initial infection density also. At the same time for , only stable disease-free equilibrium point exists which indicates disease eradicates from the system. For , we see that disease-free equilibrium point is unstable and endemic equilibrium point is stable (see Fig. 4). It is clear from theoretical analysis that there is a critical value of the number of hospital beds (b) below which this type of phenomenon occurs. Hence, to eliminate the disease from any system we have to increase the number of hospital beds, i.e., the treatment facility.

Fig. 4

Backward bifurcation diagram with respect to varying the disease transmission rate , and other parameters are given in Table 3 with , the blue line corresponds to stable branch and red line is for unstable branch

Saddle-node bifurcation

In this subsection, we shall examine the Saddle-node bifurcation with respect to the intrinsic growth rate (r). Let for the critical value the equation (3) has a coincident root.

Theorem 11

The model system (2) experiences a Saddle-node bifurcation at the coincident endemic equilibrium point for if the conditions stated in the text are satisfied.

It can be shown that for one of the eigenvalues of vanishes and other one is negative real, these occur if and . We have to check the transversality conditions for Saddle-node bifurcation using Sotomayor’s theorem [36, 39]. Let V and W be the eigen vectors corresponding to the zero eigenvalue for the matrices and , respectively, then and .

Now, and Therefore, = and if

Hence, all the conditions for Saddle-node bifurcation are satisfied for at the coincident endemic equilibrium point of the model system (2). So by Sotomayor’s theorem for Saddle-node bifurcation, we can conclude that model system (2) experiences a Saddle-node bifurcation at the coincident endemic equilibrium point when the intrinsic growth rate crosses its critical value.

In Fig. 5, we have presented the one parameter bifurcation diagram of infected individuals considering r as the bifurcation parameter for the proposed model. It is clear from the figures that for lower value of one Saddle-node bifurcation point exists, say , and for higher value of disease transmission rate there exist two Saddle-node bifurcation points, say and .

Fig. 5

Bifurcation diagram with respect to r considering different values of and other parameters are given in Table 3 with

Effect of population birth rate and disease transmission rate on model dynamics

In this subsection, numerically we shall show that population birth rate (r) and rate of infection are playing crucial role in disease spreading dynamics where other parameters are fixed as shown in Table 3 with . First, we investigate the dynamics considering the lower value of disease transmission rate then for higher value of and finally compare the results.

First, we consider keeping other parameters fixed for Saddle-node bifurcation as given in Table 3 with . For , we get only one critical value of r which is (see Fig. 5a). To interpret biological significance of the results, we consider following cases. First, consider ; in this case, there exists only one stable disease-free equilibrium point and no endemic equilibrium point exists (see the phase portrait in Fig. 6a). For , one stable disease-free equilibrium point and one saddle endemic equilibrium point exist (see the phase portrait in Fig. 6b). Next, we consider then, one stable disease-free equilibrium point exists along with two endemic equilibrium points. Among them, one endemic equilibrium point is saddle and another one is stable (see the phase portrait in Fig. 6c).

Fig. 6

Phase portrait at different values of r in the neighborhood of the Saddle-node bifurcation value for lower

Thus, we observe that the number of endemic equilibrium point changes from zero to two through Saddle-node bifurcation as bifurcation parameter r passes through its critical point from left side to right side. Here always one stable disease-free equilibrium point exists. That means in this case disease can be eradicated from the population for lower value of r, whereas for higher value of r disease persists in the system. So, we need to take necessary actions such that values of intrinsic birth rate become lower to control the infection in the population. Here we observe that for lower value of disease transmission rate the influence of the disease spreading can be decreased; moreover, disease can be eliminated from the population.

Secondly, we take and other parameter values are described in Table 3 with . For , we get two critical values of population birth rate; those are and (from Fig. 5b, we see that two critical values of r exist for higher value of ). Now we consider three cases and . In each of the above three cases, we find the number of equilibrium points and their behavior. For , the system contains two equilibrium points: one is disease-free equilibrium point and other one is endemic equilibrium point . Disease-free equilibrium point is saddle, whereas endemic equilibrium point is stable spiral (see the phase portrait in Fig. 7a). For , saddle disease-free equilibrium point and two endemic equilibrium points exist among which one is stable and other one is saddle (see the phase portrait in Fig. 7b). saddle disease-free equilibrium point and three endemic equilibrium points , , exist where , are stable spiral and is saddle. In this case, two more endemic equilibrium points arise through Saddle-node bifurcation among them one is saddle and another one stable endemic equilibrium point when r crosses from left side to right side (see the phase portrait in Fig. 8a). For , saddle disease-free equilibrium, one stable endemic equilibrium and one saddle endemic equilibrium point exist (see the phase portrait in Fig. 8b). Now for saddle disease-free equilibrium point and one stable endemic equilibrium point exist. Here two endemic equilibrium points diminish through the Saddle-node bifurcation when r crosses from left side to right side (see the phase portrait in Fig. 8c).

Fig. 7

Phase portrait at different values of r in the neighborhood of the Saddle-node bifurcation value for higher

Fig. 8

Phase portrait at different values of r in the neighborhood of the Saddle-node bifurcation value for higher

From the above three cases, we observe that when r passes from left side to right side through critical value , the number of endemic equilibrium points increases, and when r crosses left side to right side through critical value , the number of endemic equilibrium points decreases through Saddle-node bifurcation. Moreover, we see that for higher values of disease transmission rate , there always exists a stable endemic equilibrium point. It is also clear from the analysis that if values of disease transmission rate increases, then the threshold value of population birth rate for Saddle-node bifurcation also increases. That means for higher value of r one stable endemic equilibrium point exists which indicates disease always persists in this case.

Hopf bifurcation

To investigate Hopf bifurcation of the model system (2) about the endemic equilibrium point , we choose b as the bifurcation parameter fixing other model parameters. Let us suppose that roots of the equation (8) become purely imaginary at . Suppose there exists a neighborhood () of in which roots of the equation (8) are of the form with where . Now putting values of in equation (8) and separating the real and imaginary parts we getDifferentiating (11) with respect to b and putting , we getEliminating from the above two equations, we obtainSummarizing the above results, we have the following theorem:

Theorem 12

A necessary and sufficient condition for model system (2) will experience Hopf bifurcation at the endemic equilibrium point is that there exists such that

(i) and (ii)

where is complex root of Eq. (8).

Direction of Hopf bifurcation

In the above, we have derived the necessary and sufficient condition for which model system (2) experiences Hopf bifurcation. Now we deduce the normal form of Hopf bifurcation to investigate the sufficient conditions for stability of Hopf bifurcated periodic solutions [41-43] and direction of Hopf bifurcation.

Theorem 13

If , then the model system (2) undergoes through Hopf bifurcation about the endemic equilibrium point when the number of bed (b) crosses the critical value . Furthermore, (i) If , then periodic solution is stable; therefore, Hopf bifurcation is supercritical, and

(ii) , then periodic solution is unstable; therefore, Hopf bifurcation is subcritical where is defined in the text.

First, we translate the origin to the endemic equilibrium point of the system (2) using the transformation , which givesBy Taylor’s series expansion about , from (13) we getwhere , and with , , = , = ,   = , = , = = , = , = , = .

The characteristic roots of A are in the form in the neighborhood of . The eigen vector corresponding to the eigenvalue is where and P = - .

Using the transformation where which gives .

Under the above transformation the system (14) reduces towhere System (15) can be written in the polar form using the methodology developed in [41-43]which is the normal form of Hopf bifurcation.

Expanding by Taylor’s series expansion about from (16) we getwherewith

Thus, for existence of the periodic solution we must have [43]Again if , then the system follows a surface of periodic solution in the center manifold [41, 43], the periodic solution is stable and unstable for .

This completes proof of the theorem.

In Fig. 9, we have presented the bifurcation diagram with respect to the number of available hospital beds (b). It is clear from the figures that switching of instability occurs here. We see from the figures that the system losses stability at (for ) through supercritical Hopf bifurcation about the endemic equilibrium point. Again at (for , the system becomes stable through subcritical Hopf bifurcation.

Fig. 9

Bifurcation diagram with respect to parameter b and other parameter values given in Table 3 with

Significance of b as a Hopf bifurcation parameter

For the parameter values given in Table 3 with , we get two critical values of the bifurcation parameter b namely and for which with . The considered system experiences Hopf bifurcation at these two critical points. To investigate the dynamics for different values of b, we consider the following three cases: (i) (ii) (iii) In the first case, for , the system contains one saddle disease-free equilibrium point and a stable endemic equilibrium point (see the phase portrait in Fig. 10a). For the second case, we consider ; then, the endemic equilibrium point becomes unstable and a stable limit cycle arises around the endemic equilibrium point, whereas the disease-free equilibrium point remains same (see the phase portrait in Fig. 10b, the red curve is the stable limit cycle). In the third case for , the endemic equilibrium point again becomes stable and the disease-free equilibrium point remains same (see the phase portrait in Fig. 10c).

Fig. 10

Phase portrait at different values of b in the neighborhood of the Hopf bifurcation value

For , disease-free equilibrium point is saddle and endemic equilibrium point is stable which indicates disease persists in this case. For , disease-free equilibrium point is saddle and endemic equilibrium point is unstable which gives stable limit cycle that means disease oscillates within and on the limit cycle, so disease cannot be eradicated here. Further, for , disease-free equilibrium point is saddle and endemic equilibrium point is stable which means disease also persists in this case.

Optimal control technique and its application

In this section, our aim is to find the best policy among vaccination and treatment for controlling the disease as well as minimizing the cost using the technique of optimal control theory. To minimize the disease transmission, generally vaccination and treatment are performed to susceptible and infected individuals, respectively. But if constant rate of vaccination and treatment are taken into consideration for eradication of the disease, then cost for implementation of vaccine and treatment may be very high. So for eradication of the disease with minimum cost in a finite time, we need to consider control variable as a function of time t. The main purpose for applying optimal control strategy is to reduce the number of susceptible and infected individuals with minimum cost for vaccination and treatment [30]. In the proposed model, we have considered vaccination to the susceptible individuals and treatment to the infected individuals as disease controlling variable. Now we need to formulate an optimal control problem for minimizing the implemented cost for vaccine and treatment and also at same time to reduce susceptible and infected individuals.

Formulation of optimal control problem

In model system (2), we consider the treatment rate as a function of limited hospital resources that means in terms of available number of hospital beds . We have already seen that the value of per capita cure rate lies between and , whereas the value of b is finite. In order to consider the effect of treatment control, we are reconstructing treatment rate introducing new control variable in the form which has minimum and maximum value and corresponding to = 0 and = 1 with large number of infected individuals, respectively. Thus, the control variable is bounded in the function . To formulate optimal control problem, we are considering the control variables and as function of t.

The problem is to minimize objective functional given bysubject to the system of equationssatisfying , , and the controls and are Lebesgue measurable functions such that control constraint U is given byIn the optimal objective functional, constants and are the weighted cost for minimizing total number of susceptible and infected individuals, whereas and are weighted cost associated with vaccination and treatment function. We use quadratic term of the control variable in objective functional to reflect severity of side effects of treatment and vaccination for giving too much amount to an individual [44]. Our intention is to reduce susceptible and infected individuals, at the same time to maximize recovered individuals with minimize cost due to implementation of vaccine and treatment on a finite interval [0, T]. In optimal constraint U, vaccinated and treated percentage of susceptible and infected individuals per unit time is represented by control variable . Thus, we seek an optimal control pair such thatwhere U is defined in (20).

Existence and uniqueness of optimal control

Theorem 14

For optimal control problem (18) subject to the condition (19), there exist an optimal control pair and corresponding optimal state variables which minimize objective functional over set of optimal constraint such that

The existence of optimal control problem can be verified using the result of Filippove–Cesari theorem [46]. The control and corresponding state variables are non-negative and non-empty. Control variables are closed and bounded. Compactness of control variables is needed for existence of optimal control. The control constraint set U is also convex set and the state variables (x, y, z) are also bounded. Furthermore, integrand of the objective functional is also convex with respect to control variables . Finally, we can obtain and such that integrand of objection functional satisfiesTherefore, all conditions for existence of optimal control problem are satisfied and uniqueness property are guaranteed by Lipschitz property of the state variables x, y, z. Hence, the theorem is proved.

Characterization of optimal control problem

Theorem 15

Given optimal control pair and corresponding optimal state variables of state system (18) and (19) which minimize over U, there exist adjoint variables and satisfyingwith transversality conditionand the optimal control pair = is given bywhere is the non-negative root of the equation

To characterize the optimal control, we use Pontryagin’s Maximum Principle [45]. First we define Lagrangian of optimal control problem (18) asand Hamiltonian H assatisfies adjoint equationswith transversality condition

Solving equations in (26), we see that adjoint variables satisfy the following equations :with .

Now using optimality conditions and at , , we getwhere is non-negative root of the equation . It is clear that , and . > at .

Therefore, optimal control problem is minimum at = .

We have solved the optimal control problem numerically using forward–backward sweep method which combines forward application of fourth-order Runge–Kutta method of the system (19) with backward application of fourth-order Runge–Kutta method of the system (22) satisfying the condition (23). For model simulation purpose, we take the vaccination and treatment time interval as [0, 20] (where unit of time may be day or week depending on nature of the disease), i.e., vaccination and treatment continue up to 20 units of time, and after this time, both policy will automatically stop. To study the effect of control strategies here we are considering two infectious disease like influenza (low rate of infection) and COVID-19 (high rate of infection). We observe from the literature that the rate of infection (), natural death rate of human (d) and disease-induced death rate of human () of populations are in the following range: For model simulation purpose, we consider value of other model parameters and weight functions as given in Table 4 with initial population , , .

Table 4

Values of the parameters used for optimal control problem

Parameter	r	μ	K	α	μ_0	μ_1	b	A_1	A_2	A_3	A_4
Value	1.0	0.7	200	10.5	0.1	1.8	0.15	0.1	0.2	0.2	1.5

For influenza: , and [47–49, 53].

For COVID-19: , and [50-53].

In Figs. 11 and 12, we have presented the time series of susceptible, infected, recovered individuals and control strategies in presence and absence of controls for (for Influenza) with and (for COVID-19) with , , respectively. It is clear from first three figures of both cases that in the presence of control, density of infected population is lower compared to the absence of controls. In Figs. 11d, e and 12d, e we have presented the time series of control variables. It is clear from Fig. 11d, e that full strength of vaccination will continue up to 18.98 unit and full strength of treatment policy is not required for lower rate of infection (i.e., ). The controls will reduce following the path as shown in Fig. 11d, e and it will stop after 20 unit time. But for high rate of infection (see Fig. 12d, e) the full strength vaccination and treatment policy will continue up to 18.48 unit and from 1.03 unit to 4.43 unit, respectively, and it will stop after 20 unit time. If same control parameters are used beyond t = 20 unit time, then also the density of infected individuals will reduce in both cases. Thus, it is clear from the above analysis that the vaccination policy is more or less same for both lower and higher rate of infection, but higher rate of treatment is necessary for higher rate of infection.

Fig. 11

Time series of the population and control variables with control (blue line), without control (red line) for ; a susceptible individuals b infected individuals c recovered individuals d vaccination control() e Treatment control ()

Fig. 12

Time series of the population and control variables with control (blue line), without control (red line) for ; a susceptible individuals b infected individuals c recovered individuals d vaccination control () e treatment control ()

Efficiency analysis

In our study, we have used two types of control parameter and where is applied for susceptible individuals, whereas is used for infected individuals. If we want to apply only one type of control among vaccination and treatment to reduce cost as well as time, one question arises ’which one will better among them to reduce the infection?’ When two or more than two control are used in mathematical models, to know best efficient policy among them to reduce the infection, then we perform efficiency analysis. In efficiency analysis, efficiency index (E.I.) is given by the formula E.I.= where and are cumulative number of infective with and without control parameter, respectively. Higher value of efficiency index (E.I.) number is the best strategy [54].

Two types of controls and are used in our model system. Now three cases may arise. (i) When we use only vaccination for susceptible, but we do not use treatment for infected individuals, i.e., but , and (ii) when we use only treatment for infected individuals, but we do not use any vaccination for susceptible, i.e., but . (iii) When we use both vaccination for susceptible and treatment for infected individuals. We rename these three cases as strategy 1, strategy 2 and strategy 3, respectively. So strategy 1 is but , strategy 2 is , but and strategy 3 is but . Now to evaluate the value of efficiency index number (E.I.), we first calculate the values of and . We compute the value of , cumulative number of infected individuals without control during time t, by using Simpson’s rule. For , and values of and efficiency index (E.I.) are given in Table 5. For , and values of and efficiency index (E.I.) are given in Table 6. From Tables 5 and 6, it is clear that between strategy 1 and strategy 2, strategy 1 is better than strategy 2 for both cases. In view of epidemiology, vaccination for susceptible is more efficient than using treatment to infected individuals to reduce the influence of infection. But strategy 3 is more effective among three strategies for both cases. So in biological sense of view, if both vaccination and treatment are applied in the community, then the number of infected individuals is more reduced.

Table 5

Values of efficiency index for (for Influenza)

Strategy	Applied controls	A^c	E.I.
Strategy 1	u_1≠0,u_2=0	17.7283	0.5361
Strategy 2	u_1=0,u_2≠0	34.3304	0.1017
Strategy 3	u_1≠0,u_2≠0	14.5984	0.6179

Table 6

Values of efficiency index for (for COVID-19)

Strategy	Applied controls	A^c	E.I.
Strategy 1	u_1≠0,u_2=0	34.5129	0.6425
Strategy 2	u_1=0,u_2≠0	88.7237	0.0812
Strategy 3	u_1≠0,u_2≠0	28.6617	0.7032

Conclusion

In this paper, we have proposed an epidemic model with logistic growth rate of susceptible population, non-monotone incidence rate and nonlinear treatment rate with the effect of limited hospital beds. We have carried out stability and bifurcation analysis of the system in the neighborhood of different equilibrium points with respect to different model parameters. For , the disease-free equilibrium point is locally stable, and for the endemic equilibrium point is always locally stable. Using center manifold theorem, we have established that the disease-free equilibrium point shows two types of behaviors: one is saddle-node point of co-dimension one and other is an attracting node of co-dimension 2 for the threshold value . We have studied different types of bifurcations, namely Transcritical bifurcation, Backward bifurcation, Saddle-node bifurcation and Hopf bifurcation for different parameter values. Numerically we have shown bifurcation diagrams and the corresponding phase portraits for each type of bifurcation. We have shown that for lower values of both the rate of infection and birth rate of susceptible population, disease will be eradicated from the system, but in this situation if the birth rate becomes higher, then disease will persist in the system. On the hand for higher value of the rate of infection, disease persists in the system with more infected population density. We have also shown that for higher influence of disease transmission, disease cannot be eradicated from the system even higher level of treatment provided. According to our study, we see that there exists two endemic equilibrium points for , smaller one is unstable, whereas bigger one is stable which shows that is not enough to conclude the dynamical behavior of the model system. To diminish the impact of the disease spreading, we need to control the density of the population. For our proposed model, multi-stable state, stable limit cycle, stable spiral arise which indicates to control the disease, we should take necessary actions such that initial infected population is at lower level. We observed that rich dynamical behavior not only depends on , it also depends on the number of hospital beds, vaccination, inhibitory factors of susceptible, crowding effect of infected population etc. Therefore, to control or eradicate the disease, we should combined government intervention policies as well as hospitalized condition. Optimal control strategies and efficiency analysis are performed to find out the best policy to minimize the impact of the disease. Improving medical techniques and investing large number of medical resources such as medicines, doctors, isolation places and quarantine, disease can be controlled in the community.