key: cord-0698743-i9sb02mu authors: Saha, Pritam; Ghosh, Uttam title: Global dynamics and control strategies of an epidemic model having logistic growth, non-monotone incidence with the impact of limited hospital beds date: 2021-06-21 journal: Nonlinear Dyn DOI: 10.1007/s11071-021-06607-9 sha: e6c817bdde0e0afeaf06af0e4599b0490e819d2b doc_id: 698743 cord_uid: i9sb02mu 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. 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] [4] [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] [9] [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 (βx y) 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 βx y N 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 β (1−cy) x y to describe measles outbreak and also Levin and Iwasa [16] considered in the form λx p y q . 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 βx y 1 + αy where βy is the force of infection and 1 1 + αy 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 β α x 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 βx y 1 + αy 2 . 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 y = 1 √ α then starts to decrease and tends to zero as t → ∞. Further, Ruan and Wang [21] , Lu and Ruan [8] studied different types of nonlinear incidence rates, namely kx y 2 1 + αy 2 , kx y 2 1 + βy + αy 2 , 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 βx y p 1 + αy q where β, p, q > 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 T (y) = r, y > 0 0, y = 0 maximum medical resources is used here to cure the disease. Further, Wang [25] proposed a piecewise linear recovery function in the form T (y) = ky, 0 y y 0 ky 0 , y > y 0 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) = r y 1 + ky 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 r k 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 T (I ) = αy ω + y 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 form is the minimum per capita recovery rate and μ 1 (> 0) 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 μ (b, y) 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-19infected 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 β (y) = βy 1 + αy 2 and treatment function 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. 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 βx y 1 + αy 2 , and the treatment rate is of the form . Based on the above discussion, we have proposed the deterministic epidemic model in the following form: 0 for all t 0. The model parameters of the system (1) are enlisted in Table 1 . 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 x(0) 0, y(0) 0 for all t 0. 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). First, we shall show that all solutions of model system (2) satisfying the initial conditions x(0) 0, y(0) 0 are non-negative for all t 0, and then, we shall prove that the solutions are bounded for all t 0. 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 t (> 0). 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) Integrating and simplifying, we get Thus, we have established x(t) > 0 ∀ t > 0, i.e., the susceptible populations are always positive. Similarly, we can show that y(t) ≥ 0 ∀ t > 0. This completes proof of the theorem. Next we show solutions of model system (2) are uniformly bounded for all t 0. We establish our result in the following theorem : is a positively invariant attracting region for the epidemic model given by the system (2) in R 2 + . Proof Adding both equations of model system (2), we get Integrating and taking limit as t −→ ∞, we get 0, y(0) 0. 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 . In this section, we shall discuss the existence of different equilibrium points and basic reproduction number (R 0 ) of model system (2) . Basic reproduction number (R 0 ) 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 E 00 (0, 0) and a disease-free equi- Hence, R 0 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. . Let F = dF dy and V = dV dy , then the next-generation matrix is K = . Therefore, basic reproduction number of the model system (2) is given by The components of endemic equilibrium point and Since p(y) = 0 is a fifth-order equation, it is very difficult to find the expressions of positive roots of the equation explicitly in terms of R 0 . 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 R 0 and the thresholds 1 = r (μ 1 − μ 0 ) K bβ 2 Considering b = 0.57, r = 0.58, μ 1 = 1.2 and keeping all other parameters fixed as shown in Table 3 , we are varying the disease transmission rate β (consequently R 0 will vary) and have drawn the solution curve for infected population (i.e., for existence of endemic equilibria) with respect to R 0 (see Fig. 1 ). The number of equilibrium points for different values of R 0 is summarized in the following lemma: 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 E 00 (0, 0) is stable if r < (d + u 1 ) and saddle point if r > (d + u 1 ). Proof The Jacobian matrix (4) at the trivial equilibrium point E 00 (0, 0) is: Eigenvalues of the above Jacobian matrix are given by r − (d + u 1 ) and −(d + μ + δ + μ 1 ). Clearly, one eigenvalue is negative and other will be negative if r < (d + u 1 ), and hence, under this restriction E 00 (0, 0) is stable. If r > (d + u 1 ) holds, then one eigenvalue is positive and other value is negative; consequently, E 00 (0, 0) is a saddle point. In biological point of view, r < (d +u 1 ) 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., u 1 is nearer to r such that (u 1 + d) > r ) 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, r > (d + u 1 ) means birth rate of the population is higher than sum of vaccination and death rate of population which is realistic. Theorem 5 If R 0 < 1, then the disease-free equilibrium point E 0 is locally asymptotically stable, and for and a semi-hyperbolic attracting Proof Jacobian matrix of the system (2) at the diseasefree equilibrium point E 0 is . If R 0 < 1, then both the eigenvalues are negative, and hence, E 0 is stable and is unstable(saddle) for R 0 > 1. For R 0 = 1, one of the above eigenvalues vanishes, and hence, the disease-free equilibrium point E 0 is a non-hyperbolic equilibrium point. Use center manifold theorem to determine the nature of equilibrium point E 0 . First we translate disease-free equilibrium point , 0 to origin using the trans- ; then, the model system (2) takes the following form: Introducing the transformation S = − β K r x + y, I = x, the System (5) becomes : It is clear from (6) that if , then use center manifold theorem [36] for investigating nature of the disease-free equilibrium point E 0 . Let the local center manifold of the system is given by 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 y = h(x) = αβ K r x 3 , and then, the system (6) reduces to From system (7), it is clear that disease-free equilibrium point E 0 is a semi-hyperbolic attracting node of co- . Hence, the theorem is proved. Next we shall establish the local stability of the endemic equilibrium E(x * , y * ) in the following theorem: Proof Jacobian matrix of the model system (2) at the endemic equilibrium E(x * , y * ) is The corresponding characteristic equation is Roots of the above equation will have negative real Hence, endemic equilibrium E(x * , y * ) is locally 3 and is unstable otherwise. Thus, we arrive at our desired result. 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 R 0 . Eradication and persistence of the disease can be determined by stability of disease-free equilibrium point , 0 and the endemic equilibrium E(x * , y * ). First we shall show the disease-free equilibrium point E 0 is globally stable when R 0 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 E 0 . Proof Using Lyapunov's stability theorem [37, 38] , we shall establish the global stability of disease-free equilibrium point. Let, V : i.e., So, by Lyapunov's Stability Theorem [37, 38] The above result is epidemiologically most important because if the value of R 0 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 R 0 is less than that threshold value. Now we shall prove global stability of endemic equilibrium point E(x * , y * ). 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. satisfied. Proof We write our model system (2) as given below: To show global stability of endemic equilibrium E(x * , y * ), we use Dulac function criterion [25, 36] of stability analysis. Considering Dulac function as . Combining the above two conditions, we get Thus, by Dulac criterion we can conclude that system (2) has no closed orbit that means endemic equilibrium point E(x * , y * ) is globally stable for R 0 > 1 under the conditions stated in the theorem. we have shown in the above theorem that for R 0 > 1 there exist a certain region of the model parameter β where endemic equilibrium E(x * , y * ) 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. 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 (μ 1 ), and it gives the upper limit of μ 1 above which disease will eliminate from the system and below which disease will persist. Again Backward bifurcation is studied with respect to R 0 . 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 R 0 but also on initial number of infection. If the value of intrinsic growth rate r passes through its critical value r * , 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 E 1 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. Theorem 9 The model system (2) experiences Transcritical bifurcation at the disease-free equilibrium Proof Since at R 0 = 1, i.e., μ 1 = μ * 1 , one of the eigenvalues of J E 0 is zero and other is negative, to investigate the nature of the solution at E 0 , we have to use Sotomayor's theorem [36, 39] . For this purpose, we construct a function with two components given by Let corresponding to zero eigenvalue J E 0 (μ * 1 ) and Therefore, Hence, the model system (2) satisfies all conditions of Sotomayor's theorem for Transcritical bifurcation with respect to parameter μ 1 . Thus, model system (2) experiences Transcritical bifurcation in the neighborhood of the critical value of the maximum value of the recovery parameter μ 1 = μ * 1 . Here, numerically we shall interpret significance of the maximum value of the per capita recovery rate ( μ 1 ) as the Transcritical bifurcation parameter for lower and upper value of the rate of infection(β). We consider the lower value as β = 0.09 and the upper value as β = 0.9. For β = 0.09 and other parameter values are described as shown in Table-3 with b = 1.24, r = 0.98, we get the critical value of Transcritical bifurcation parameter as μ * 1 = 1.689795918. (i) For μ 1 = 1.63 < μ * 1 , one saddle disease-free equilibrium point E 0 and one stable endemic equilibrium point E 1 exist (see the corresponding phase portrait in Fig. 2a) . (ii) for μ * 1 < μ 1 = 1.75, there exist one stable disease-free equilibrium point E 0 and two endemic equilibrium points E 1 and E 2 where E 1 is saddle point and E 2 is stable spiral (see the phase portrait in Fig. 2b) . Thus, when μ 1 passes from left side to right side through its critical value μ * 1 , then the disease-free equilibrium point E 0 and the endemic equilibrium point E 1 exchange their stability through Transcritical bifurcation with creation of one stable endemic equilibrium point. (iii) If we take value of μ 1 far from the transcritical value in right side of μ * 1 say μ 1 = 5, then only a stable disease-free equilibrium point exists (see the phase portrait in Fig. 2c ). It is clear from the above discussion that there exists a critical value of the maximum value of treatment rate (μ 1 ) below which the disease will persist in the system for any initial population density. But for the values of μ 1 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. For β = 0.9 and other parameter values are given in Table 3 with b = 1.24, r = 0.98, we get the critical value of the maximum cure rate as μ * 1 = 24.99795918. Now, for μ 1 = 28 > μ * 1 , 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 μ 1 = 20 < μ * 1 , there exist a unstable (saddle) diseasefree 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 μ 1 more in the left side of μ * 1 , say μ 1 (= 4) < μ * 1 , 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 μ * 1 . From the above, we have seen that for higher treatment rate disease can be eradicated from the population. In just left side of μ * 1 , i.e., for μ 1 (= 20) < μ * 1 , disease will persist in the system with oscillatory population density, and for very lower value of μ 1 , i.e., for μ 1 (= 4) < μ * 1 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. Theorem 10 The model system (2) experiences Back- Proof 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 and left eigen vector, respectively, of J E 0 corresponding to the eigenvalue λ = 0. In order to apply Castillo-Chavez and Song's bifurcation theorem [40] , we need to compute the value of φ and ψ which are given by Thus, the model system (2) experiences Backward In Fig. 4 , we have plotted the bifurcation diagram with respect to R 0 varying the disease transmission rate β. It is clear from the figure that there is a critical value of R 0 , say R * 0 , such that in 0 < R 0 < R * 0 there exists only one stable disease-free equilibrium point, for R * 0 < R 0 < 1, there exist two endemic equilibrium points along with the stable disease-free equilibrium point, and for R 0 > 1 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 Table 3 with r = 0.58, b = 0.57, μ 1 = 1.2, the blue line corresponds to stable branch and red line is for unstable branch means for R * 0 < R 0 < 1, 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 R 0 as well as on initial infection density also. At the same time for 0 < R 0 < R * 0 , only stable disease-free equilibrium point exists which indicates disease eradicates from the system. For R 0 > 1, 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. In this subsection, we shall examine the Saddle-node bifurcation with respect to the intrinsic growth rate (r ). Let for the critical value r = r * the equation (3) has a coincident root. The model system (2) experiences a Saddle-node bifurcation at the coincident endemic equilibrium point E(x * , y * ) for r = r * if the conditions stated in the text are satisfied. Proof It can be shown that for r = r * one of the eigenvalues of J E vanishes and other one is negative real, these occur if det (J E )| r =r * = 0 and tr(J E )| r =r * > 0. 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 J E and J T E , respectively, then V = Hence, all the conditions for Saddle-node bifurcation are satisfied for r = r * at the coincident endemic equilibrium point E(x * , y * ) 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 E(x * , y * ) 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 r * , and for higher value of disease transmission rate (β) there exist two Saddle-node bifurcation points, say r * 1 and r * 2 . 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 b = 0.57, μ 1 = 1.2. 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 β = 0.091 keeping other parameters fixed for Saddle-node bifurcation as given in Table 3 with b = 0.57, μ 1 = 1.2. For β = 0.091, we get only one critical value of r which is r * = 0.5835215509 (see Fig. 5a ). To interpret biological significance of the results, we consider following cases. (i) First, consider r (= 0.50) < r * ; 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 ). (ii) For r = r * , one stable disease-free equilibrium point and one saddle endemic equilibrium point exist (see the phase portrait in Fig. 6b) . (iii) Next, we consider r * < r (= 0.62); 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) . Thus, we observe that the number of endemic equilibrium point changes from zero to two through Saddlenode bifurcation as bifurcation parameter r passes through its critical point r * 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 (r ) 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. Table 3 Secondly, we take β = 0.91 and other parameter values are described in Table 3 with b = 0.57, μ 1 = 1.2. For β = 0.91, we get two critical values of population birth rate; those are r * 1 = 1.164851284 and r * 2 = 1.283695472 (from Fig. 5b , we see that two critical values of r exist for higher value of β). Now we consider three cases (i) r < r * 1 (ii) r * 1 < r < r * 2 and (iii) r * 2 < r . In each of the above three cases, we find the number of equilibrium points and their behavior. (i) For r (= 1.10) < r * 1 , the system contains two equilibrium points: one is disease-free equilibrium point (E 0 ) and other one is endemic equilibrium point (E 1 ). Disease-free equilibrium point is saddle, whereas endemic equilibrium point is stable spiral (see the phase portrait in Fig. 7a ). For r = r * 1 , saddle disease-free equilibrium point E 0 and two endemic equilibrium points exist among which one is stable and other one is saddle (see the phase portrait in Fig. 7b) . (ii) r * 1 < r (= 1.22) < r * 2 saddle disease-free equilibrium point and three endemic equilibrium points E 1 , E 2 , E 3 exist where E 1 , E 3 are stable spiral and E 2 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 r = r * 1 from left side to right side (see the phase portrait in Fig. 8a ). For r = r * 2 , saddle disease-free equilibrium, one stable endemic equilibrium and one saddle endemic equilibrium point exist (see the phase portrait in Fig. 8b ). (iii) Now for r * 2 < r (= 1.35) saddle disease-free equilibrium point and one stable endemic equilibrium point E 1 exist. Here two endemic equilibrium points diminish through the Saddle-node bifurcation when r crosses r = r * 2 from left side to right side (see the phase portrait in Fig. 8c ). From the above three cases, we observe that when r passes from left side to right side through critical value r * 1 , the number of endemic equilibrium points increases, and when r crosses left side to right side through critical value r * 2 , the number of endemic equilibrium points decreases through Saddle-node bifurcation. Moreover, we see that for higher values of dis- . Now putting values of λ in equation (8) and separating the real and imaginary parts we get Differentiating (11) with respect to b and putting b = b * , we get Eliminating dG(b * ) db from the above two equations, we obtain Summarizing 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 E(x * , y * ) is that there where λ is complex root of Eq. (8). 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] [42] [43] and direction of Hopf bifurcation. Proof First, we translate the origin to the endemic equilibrium point E(x * , y * ) of the system (2) using the transformation By Taylor's series expansion about (x * , y * ), from (13) we get where A = a 10 a 01 b 10 b 01 , φ 1 (X, Y ) = a 20 X 2 + a 11 XY + a 02 Y 2 + a 12 XY 2 + a 03 Y 3 and φ 2 (X, The characteristic roots of A are in the form F ± i G in the neighborhood of b = b * . The eigen vector corresponding to the eigenvalue Using the transformation Under the above transformation the system (14) reduces to where f (u, v) = (a 20 + N a 11 + N 2 a 02 )u 2 + (Pa 11 + 2N Pa 02 )uv + P 2 a 02 v 2 + (N 2 a 12 + N 3 a 03 )u 3 + (2N Pa 12 +3N 2 Pa 03 )u 2 v+(3N P 2 a 03 + P 2 a 12 )uv 2 + P 3 a 03 v 3 System (15) can be written in the polar form using the methodology developed in [41] [42] [43] which is the normal form of Hopf bifurcation. Expanding by Taylor's series expansion about b = b * , from (16) we get where Thus, for existence of the periodic solution we must have [43] Again if c 1 (b * ) = 0, then the system follows a surface of periodic solution in the center manifold [41, 43] , the periodic solution is stable c 1 (b * ) < 0 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 P 1 (for b = b * 1 = 0.3608114828) through supercritical Hopf bifurcation about the endemic equilibrium point. Again at P 2 (for b = b * 2 = 1.185564499), the system becomes stable through subcritical Hopf bifurcation. For the parameter values given in Table 3 with β = 0.2, r = 0.5, μ 1 = 1.8, we get two critical values of the bifurcation parameter b namely b * 1 = 0.3608114828 and b * 2 = 1.185564499 for which tr(J E ) = 0 with det (J E ) > 0. 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: In the first case, for b = 0.3, the system contains one saddle diseasefree equilibrium point and a stable endemic equilibrium Table 3 with β = 0.2, r = 0.5, μ 1 = 1.8 point (see the phase portrait in Fig. 10a) . For the second case, we consider b = 0.75; 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 b = 1.25, the endemic equilibrium point again becomes stable and the disease-free equilibrium point remains same (see the phase portrait in Fig. 10c ). For b < b * 1 , disease-free equilibrium point is saddle and endemic equilibrium point is stable which indicates disease persists in this case. For b * 1 < b < b * 2 , diseasefree 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 b * 2 < b, disease-free equilibrium point is saddle and endemic equilibrium point is stable which means disease also persists in this case. 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. 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 (b). We have already seen that the value of per capita cure rate lies between μ 0 and μ 1 , 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 u 2 in has minimum and maximum value μ 0 and μ 1 corresponding to u 2 = 0 and u 2 = 1 with large number of infected individuals, respectively. Thus, the control variable u 2 is bounded in the function (u 2 ). To formulate optimal control problem, we are considering the control variables u 1 (t) and u 2 (t) as function of t. The problem is to minimize objective functional given by satisfying x(0) 0, y(0) 0, z(0) 0 and the controls u 1 (t) and u 2 (t) are Lebesgue measurable functions such that control constraint U is given by In the optimal objective functional, constants A 1 (> 0) and A 2 (> 0) are the weighted cost for minimizing total number of susceptible and infected individuals, whereas A 3 (> 0) and A 4 (> 0) 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 (u 1 (t), u 2 (t)) ∈ U . Thus, we seek an optimal control pair u * 1 , u * 2 such that where U is defined in (20) . Theorem 14 For optimal control problem (18) subject to the condition (19) , there exist an optimal control pair (u * 1 , u * 2 ) and corresponding optimal state variables (x * , y * , z * ) which minimize objective functional M(u 1 , u 2 ) over set of optimal constraint U = {(u 1 (t), u 2 (t)) | 0 Proof 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 u 1 , u 2 ∈ U 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 M(u 1 , u 2 ) is also convex with respect to control variables u 1 , u 2 . Finally, we can obtain w 1 , w 2 > 0 and ρ such that integrand of objection functional M(u 1 , u 2 ) satisfies Therefore, 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. Theorem 15 Given optimal control pair (u * 1 , u * 2 ) and corresponding optimal state variables (x * , y * , z * ) of state system (18) and (19) which minimize M(u 1 , u 2 ) over U , there exist adjoint variables λ 1 , λ 2 and λ 3 satisfying with transversality condition λ i (T ) = 0, i = 1, 2, 3 and the optimal control pair u * (t) = (u * 1 (t), u * 2 (t)) is given by where u 2 is the non-negative root of the equation Proof To characterize the optimal control, we use Pontryagin's Maximum Principle [45] . First we define Lagrangian of optimal control problem (18) as and Hamiltonian H as satisfies adjoint equations with transversality condition λ i (T ) = 0, i = 1, 2, 3 Solving equations in (26), we see that adjoint variables λ 1 , λ 2 , λ 3 satisfy the following equations : with λ 1 (T ) = λ 2 (T ) = λ 3 (T ) = 0. where u 2 is non-negative root of the equation A 4 u 2 (b+ u 2 y * ) 2 = b 2 y * (λ 2 − λ 3 )(μ 1 − μ 0 ). It is clear that . Therefore, optimal control problem is minimum at u * (t) = (u * 1 (t), u * 2 (t)). 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 diseaseinduced 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 A i , i = 1, 2, 3, 4 as given in Table 4 with initial population x(0) = 100000, y(0) = 10000, z(0) = 0. 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 β = 0.1 (for Influenza) with d = 0.0147, δ = 0.07 and β = 0.9 (for COVID-19) with d = 0.0147, δ = 0.038, 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., β = 0.1). 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. In our study, we have used two types of control parameter u 1 and u 2 where u 1 is applied for susceptible individuals, whereas u 2 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.= 1 − A c A 0 × 100 where A c and A 0 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 u 1 and u 2 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., u 1 = 0 but u 2 = 0, and (ii) when we use only treatment for infected individuals, but we do not use any vaccination for susceptible, i.e., u 1 = 0 but u 2 = 0. (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 u 1 = 0 but u 2 = 0, strategy 2 is u 1 = 0, but u 2 = 0 and strategy 3 is u 1 = 0 but Table 5 . For β = 0.9, A 0 = 96.5638 and values of A c 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. 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 R 0 < 1, the disease-free equilibrium point is locally stable, and for R 0 > 1 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 R 0 = 1. 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 R 0 < 1, smaller one is unstable, whereas bigger one is stable which shows that R 0 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 R 0 , 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. A stochastic SIRS epidemic model with infectious force under intervention strategies Complex dynamics in a unified SIR and HIV disease model: a bifurcation theory approach The Mathematical Theory of Infectious Diseases and Its applications Mathematical Models in Population Biology and Epidemiology Infectious Diseases of Humans: Dynamics and control An attempt at a new analysis of the mortality caused by smallpox and of the advantages of inoculation to prevent it, reprint A contribution to mathematical theory of epidemics Bifurcation analysis of an SIRS epidemic model with a generalized non-monotone and saturated incidence rate Hopf bifurcation of a delayed SIQR epidemic model with constant input and nonlinear incidence rate Stability analysis of an SIR epidemic model and effect of control strategies with constant recruitment Stability and bifurcation analysis of an SIR epidemic model with logistic growth and saturated treatment Analysis of an SIR model with bilinear incidence rate Bifurcation analysis of a discrete SIRS epidemic model with standard incidence rate An Introduction to Mathematical Epidemiology Recurrent outbreaks of measles, chickenpox and mumps Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models A generalization of the Kermack-Mckendric deterministic epidemic model Complex dynamics of an SIR epidemic model with saturated incidence rate and treatment Complex dynamics of an SEIR epidemic model with saturated incidence rate and treatment Global analysis of an epidemic model with non-monotone incidence rate Dynamical behavior of an epidemical model with a nonlinear incidence rate Dynamical behavior of epidemiological models with nonlinear incidence rates Some epidemiological models with nonlinear incidence Bifurcation in an epidemic model with constant removal rate of the infectives Backward bifurcation of an epidemic model with treatment Backward bifurcation of an epidemic model with saturated treatment function Dynamics of an SIR epidemic model with limited medical resources revisited Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds Modeling the spread and control of dengue with limited public health resources Optimal control applied to biological model A theoretical study on mathematical modelling of an infectious disease with application of optimal control Stability analysis and optimal vaccination of an SIR epidemic model Optimal treatment of an SIR epidemic model with time delay Optimal control applied to vaccination and treatment strategies for various epidemiological models Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease trans-mission Introduction to Applied Nonlinear Dynamical System and Chaos Stability of motion, with a contribution by V. A. Pliss and an introduction by V Global stability of infectious disease models using Lyapunov functions Differential Equations and Dynamical Systems Dynamical model of tuberculosis and their applications The Hopf Bifurcation and Its Applications Applied Nonlinear Dynamics: Analytical Melnikov Functions and Bifurcations of Limit Cycles Optimal control in epidemiology The Mathematical Theory of Optimal Processes A Filippov-type lemma for functions involving delays and its application to time delayed optimal control problems Modelling the epidemic spread of an H1N1 influenza outbreak in a rural university town Parameter estimation, sensitivity and control strategies analysis in the spread of influenza in Mexico Modeling influenza transmission dynamics with media coverage data of the 2009 H1N1 outbreak in Korea COVID-19 pandemic in India: a mathematical model study Forecasting the daily and cumulative number of cases for the COVID-19 pandemic in India A simulation of a COVID-19 epidemic based on a deterministic SEIR model Modelling the dynamics of dengue real epidemics Mathematical modeling of dengue epidemic: control methods and vaccination strategies Data availability Data of this study will be made available from the corresponding author on reasonable request. Conflict of interest The authors of this paper declare no conflicts of interest.