key: cord-0079122-h36zjrv1
authors: Saha, Pritam; Ghosh, Uttam
title: Complex dynamics and control analysis of an epidemic model with non-monotone incidence and saturated treatment
date: 2022-05-26
journal: Int J Dyn Control
DOI: 10.1007/s40435-022-00969-7
sha: 2c1844f2fec31bd596e868ab82032747cb28f210
doc_id: 79122
cord_uid: h36zjrv1

In this manuscript, we consider an epidemic model having constant recruitment of susceptible individuals with non-monotone disease transmission rate and saturated-type treatment rate. Two types of disease control strategies are taken here, namely vaccination for susceptible individuals and treatment for infected individuals to minimize the impact of the disease. We study local as well as global stability analysis of the disease-free equilibrium point and also endemic equilibrium point based on the values of basic reproduction number [Formula: see text] . Therefore, disease eradicates from the population if basic reproduction number less than unity and disease persists in the population if basic reproduction number greater than unity. We use center manifold theorem to study the dynamical behavior of the disease-free equilibrium point for [Formula: see text] . We investigate different bifurcations such as transcritical bifurcation, backward bifurcation, saddle-node bifurcation, Hopf bifurcation and Bogdanov–Takens bifurcation of co-dimension 2. The biological significance of all types of bifurcations are described. Some numerical simulations are performed to check the reliability of our theoretical approach. Sensitivity analysis is performed to identify the influential model parameters which have most impact on the basic reproduction number of the proposed model. To control or eradicate the influence of the emerging disease, we need to control the most sensitive model parameters using necessary preventive measures. We study optimal control problem using Pontryagin’s maximum principle. Finally using efficiency analysis, we determine most effective control strategy among applied controls.

Epidemiology is a part of mathematical biology which mainly associates with the study of infectious and vectorborne disease and its control strategies in a human population. To study the influence and effect of the communicable disease, researchers often used mathematical models which lead a significant role in theory of epidemiology [1, 2] . Mathematical models are used to provide understanding the dynamical behavior of the disease and to improve plans for eradicating the infection. In epidemiology, first known mathematical model is invented and worked by Daniel Bernoulli [3] . In B Uttam Ghosh uttam_math@yahoo.co.in Pritam Saha pritamsaha1219@gmail.com 1 Department of Applied Mathematics, University of Calcutta, Kolkata 700009, India the middle of twentieth century, mathematical epidemiology strode widely around the world [4] [5] [6] . The father of modern mathematical epidemiology, Ronald Ross, derived a threshold quantity, called as basic reproduction number, which helps the planners to conclude when the infection will eliminate. In 1927, Kermack and McKendrick wielded mathematical epidemiology to a new level by his contribution for infectious disease [7] . Later on, researchers formulate various mathematical models for developing strategies to control disease outbreak and making a trade-offs in choosing a best possible treatment for seizing the infection [8] [9] [10] .

In mathematical modeling, recruitment rate of susceptible population, incidence rate of disease transmission and recovery rate have been playing a crucial role in deterministic approach of mathematical epidemiology. Mainly two types of recruitment rates are considered to formulate mathematical models. Constant birth rate has been taken into consideration for discussing compartmental model by the authors in [11, 12] , which is not realistic in many certain cases. Under this circumstances, logistic growth rate of susceptible has been used to make mathematical models more realistic in case of limited resources and space capacities, disease with long lasting disease [13, 14] .

Incidence rate has a vital role for discussing nature of a disease in theory of epidemiology. In early, bilinear incidence rate in the form βx y and standard incidence rate in the form βx y N are used in epidemic model where β is the rate of disease transmission per contact, x and y denote number of susceptible population and number of infected population respectively, and N is the total number of population [15, 16] . In such cases, epidemic model has relatively simple dynamical behavior that means model does not have any periodic oscillations and condition for persistence and eradication of the infection can be determined simply by calculating basic reproduction number (R 0 ). The authors established that if R 0 < 1, then disease will eliminate and otherwise disease will permanent [17] . But these types of incidence rates are not realistic as they do not include social awareness of susceptible population and crowding effect of infected population. To get rich dynamical phenomenon of the infectious disease, many authors studied nonlinear types of incidence rate [18] [19] [20] [21] . To avoid unboundedness of infection, Capasso and Serio [22] , Anderson and May [23] introduced saturated types of infection rate in the form βx y 1 + α y where β y denotes force of infection and 1 1 + α y measures inhibition factors due to crowding effect of infected class and also for social awareness of the population, α is inhibition coefficient. Later on, to get rich dynamics of infectious disease that means to get periodic oscillations, multiple peaks of infection, Ruan and Xiao [24] introduced a non-monotone type of infection rate in the form βx y 1 + α y 2 to include psychological effect of population.

The function y 1 + α y 2 has maximum value at y =

is increasing first upto y = 1 √ α , then decreases to zero as y goes to infinity. Here, initially infection rate increases due to lack of knowledge about the disease, but taking appropriate awareness and protection measures, incidence rate gradually diminishes by avoiding unnecessary contacts [10, 25] . Thereafter, many researchers used general types of nonlinear incidence rate for describing critical situation of disease spreading [21, 26] . Spreading of the disease and influence of the infection can be controlled by taking proper treatment and giving suitable vaccine to infected class. Usually, in classical epidemic model, treatment function of infection is taken as proportional to number of infectives, i.e., every infected should take up treatment for a curable disease. But there is a maximum capacity of medical resources for treating against the disease. So, for largely populated countries, this type of recovered rate is unusable for treatment of the disease. Therefore, saturated types of treatment functions are more realistic as they tend to a finite limit when infected populations increase. In order to include limited number of medical resources such as hospital facilities, medicines, doctors etc., Wang and Ruan [27] considered constant treatment function T (y) into epidemic model where

where a constant capacity of medical resources is used to cure from the infection. They studied that model experiences Hopf bifurcation, saddle-node bifurcation, Bogdanov-Takens bifurcation varying two parameters of the system. Later on, Wang [28] introduced a piecewise treatment function which is proportional to the number of infected before reaching its maximum capacity, by

where ky 0 is maximum number of medical resources available for curing the disease. In this model, authors established that model undergoes with backward bifurcation and also observed that basic reproduction number (R 0 ) is not remain as a threshold quantity. The authors in [14, 27, 29, 30] considered treatment function as h(y) = r y 1 + α y for discussing the limited number of available medical facilities. Above function h(y) has maximum available limit r α as medical treatment. For a large number of population, treatment may be delayed to receive infected class. This function includes the parameter α as a extent of effect of infective as delaying the treatment [31] .

Our objective is to reduce spread of infection and eradicate disease from the community. Therefore, we have to find an optimal control strategy which is suitable for minimizing the invade of infection, i.e., we should find an optimal value of applied controls to reduce the infection. The aim of applying optimal control is to reduce the density of infected population and minimize control cost [32] . Sometimes, researchers used only vaccination to reduce the infection [33, 34] , whereas researchers considered only treatment function to control the invade of infection [35] . Further both treatment and vaccination used to minimize the spread of infection [13, 36] . In this paper, we consider both vaccination and treatment to control the spread of the disease.

The motivation behind this paper is the extension of the work [19, 24] . In [19] , authors proposed an epidemic model having constant birth rate of susceptible population with saturated incidence rate and saturated treatment rate. They studied biological properties of the considered model and also proved local as well as global stability of different equilibrium points. They showed model system exhibited through backward bifurcation; moreover, they studied optimal control of the proposed model. In [24] , authors considered an epidemic model having constant birth rate of susceptible population with non-monotone incidence rate. They showed stability analysis of equilibrium points. In both papers, authors gave less attention to study critical dynamical behavior like Hopf bifurcation, Bogdanov-Takens bifurcation, etc. In this paper, we have considered an epidemic model having constant birth rate of susceptible population with non-monotone incidence rate and saturated treatment rate. After model formulation, we are interested to examine more critical dynamical behavior of the model such as transcritical bifurcation, saddle-node bifurcation, Hopf bifurcation, Bogdanov-Takens bifurcation, etc. Moreover, we try to find an optimal paths of applied controls to reduce the spread of infection by studying optimal control theory.

Organization of this paper: In Sect. 2, we formulate model and discuss basic properties of the model system such as positivity and boundedness of the solution of the model system. In Sect. 3, we examine different types of equilibrium points and also find basic reproduction number (R 0 ). In Sect. 4, we examine local as well as global stability of the disease-free equilibrium point (R 0 ) and also endemic equilibrium point (E * ). In Sect. 5, bifurcation of co-dimension 1, namely transcritical bifurcation, backward bifurcation, saddle-node bifurcation, Hopf bifurcation and bifurcation of co-dimension 2, namely Bogdanov-Takens bifurcation of co-dimension 2 are investigated. In Sect. 6, numerical simulations are performed to check the reliability of theoretical findings. Sensitivity analysis is done to identify the most influential model parameters which have most affect on basic reproduction number R 0 in Sect. 7. Optimal control problem and efficiency analysis is performed in Sect. 8 and 9. Finally, critical dynamical behavior of the model system, and its biological significance are discussed in Sect. 10.

Let the total populations are classified into three compartments, namely susceptible individuals (x(t)), infected individuals (y(t)) and recovered individuals (z(t)) at any time t. In this study, we consider an epidemic model having constant birth rate of susceptible class with non-monotone incidence (affected by inhibitory factors, namely social awareness, psychological factors, etc.) in the form βx y 1 + α y 2 with vaccination and saturated treatment in the form au 2 y 1 + bu 2 y .

In reality, at first disease spreads rapidly among the population due to their unconsciousness. But for psychological effect as well as awareness of susceptible population and also for crowding effect of infected population, transmission rate of the disease starts to decrease after someday. Therefore, considering non-monotone incidence rate is appropriate for describing the spread of any disease. Also there is a maximum capacity for treatment facility. When the number of infected peoples are very small, then every people gets treatment as saturated treatment function shows linear type behavior for this case. But for huge number of infected peoples, the saturated treatment tends to its maximum capacity. Therefore, in this case, every people does not get treatment, which is realistic. For these reasons, we have considered here nonmonotone incidence rate and saturated treatment function. Two types of death are considered, one for normal death and another for disease induced death. Some infected peoples are recovered for natural immunity from infected class at a rate μ. Considering all of above assumptions, the differential equation of the proposed model can be placed in the following form:

satisfying x(0) ≥ 0, y(0) ≥ 0, z(0) ≥ 0 with nonnegative model parameters are enlisted in Table 1 , and flow diagram of the proposed model is given in Fig. 1 . We notice that recover class (z) does not appear in the first two equations of the model (1), thus we can omit third equation for analysis purpose of the system and focus only on the following subsystem of (1) which is proposed below:

First, we show that all the solutions of model (2) with nonnegative initial condition remain nonnegative for t ≥ 0. 

Using the integrating factor, we obtain the solution of the above differential equation in the form:

i.e., x(t) > 0 for all t ≥ 0. Similarly, we can show y(t) ≥ 0 ∀t ≥ 0. This completes the proof.

In the next theorem, we shall establish the boundedness of the solutions for all t ≥ 0. Proof First, we assume that N (t) = x(t) + y(t).

Adding both equations in (2), we get In the next section, we shall describe the conditions for existence of different equilibrium points.

The equilibria of the model (2) are solutions of the following equations:

It is clear from system (5) that the model (2) always has a disease-free equilibrium

As the disease-free equilibrium E 0 exist, the basic reproduction number (R 0 ) of the model (2) exists. In the following theorem, we compute basic reproduction number of the proposed model using next-generation matrix method [37] .

.

Proof Basic reproduction number plays a crucial role in the disease spreading dynamics. The number of secondary infections produced by a single infected individual in its entire life span as an infectious host is defined as basic reproduction number. Based on the values of basic reproduction number, epidemiologist can predict when the disease persists or eradicates from the population. Now, we find basic reproduction number of the proposed model using next-generation matrix approach [37] . In next-generation matrix approach, we consider two matrices F and V with containing newly infected terms and the remaining terms, respectively. Therefore, F = βx y 1 + α y 2 and V = (d + μ + δ)y + au 2 y 1 + bu 2 y . Differentiating F and V with respect to infected component y

. Now, nextgeneration matrix K is given by K = F V −1 . Basic reproduction number is the spectral radius of the nextgeneration matrix K . Here,

Therefore, basic reproduction number of the proposed model system (2) is given by

The endemic equilibria of the considered model system is

and y * is the positive root of the cubic equation

where

Discriminant [38] of cubic equation (7) is

the expression of ( ) can be written in the form:

Here, coefficients A 2 and A 3 are always positive and sign of A 0 depends on R 0 in (7) . So, we have the following cases for positive roots of Q(y) = 0 : Summarizing the above discussions, we get following theorem.

Theorem 4 (1) Model system (2) always has a disease-free equilibrium point E 0 . (2) When R 0 > 1, model system (2) has a unique endemic equilibrium.

. (c) For A 1 < 0 and also 0 = 0, these two endemic equilibria coalesce into the same endemic equilibrium E * (x * , y * ).

Otherwise, model system (2) has no endemic equilibrium.

For discussing stability of the equilibrium point, we need to compute Jacobian matrix J (E(x * , y * )) of model system (2) at any equilibrium point E (x * , y * ). The Jacobian matrix of model (2 ) is given by

Next, we shall prove local stability of disease-free equilibrium point for R 0 < 1 in Theorem 5. Then, we shall show global stability of disease-free equilibrium for R 0 < 1 using Lyapunov stability theorem in Theorem 6. Further, we shall discuss local stability of endemic equilibrium for R 0 > 1 in Theorem 7, and next we shall prove global stability of endemic equilibrium using Dulac function criterion for R 0 > 1 in Theorem 8.

First, we discuss local as well as global stability at diseasefree equilibrium point E 0 .

Theorem 5 If R 0 < 1, then disease-free equilibrium point E 0 is locally asymptotically stable and unstable for R 0 > 1.

Proof Jacobian matrix of the proposed model (2) at diseasefree equilibrium point E 0 is given by

. Therefore, for R 0 < 1, both eigenvalues are negative, and for R 0 > 1, one of them is positive. So, we can conclude that E 0 is locally asymptotically stable for R 0 < 1 and unstable for R 0 > 1. For R 0 = 1, one eigenvalue of the Jacobian matrix J E 0 is zero; therefore, in this case, disease-free equilibrium point E 0 is non-hyperbolic type, and to examine dynamical behavior of the disease-free equilibrium point E 0 , we have to use center manifold theorem [39] .

First, we use a perturbation at disease-free equilibrium

, Y = y, thus model system (2) becomes:

Next, we use the transformation

From (13), it is clear that if 2 2 , then we have to use center manifold theorem to find the nature of the disease-free equilibrium point E 0 .

By definition of local center manifold of the system

satisfying h(0) = 0, Dh(0) = 0 where δ is sufficiently small, Dh represents the derivative of h with respect to u. Using the conditions, we get the value of a 1 as

Therefore, local manifold of the system is given by v = h(u) = abu 2 2 d + u 1 u 2 , and thus system (13) becomes

From (15), it is clear that if 2 2 , then diseasefree equilibrium point E 0 is a semi-hyperbolic attracting node of co-dimension 2 for α A ≥ abu 2 2 . Thus, proof of the theorem is complete.

Proof To show the global asymptotic stability of the diseasefree equilibrium point E 0 , we are using here Lyapunov's stability theorem [40, 41] . Let L(x, y) = y is a Lyapunov function as it is positive definite. Then,

This result is epidemiologically more significant since if basic reproduction number is less than a threshold quantity, then disease eradicates from the population. Therefore to control or eradicate the emerging disease, our aim will be to maintain the value of the basic reproduction number less than that threshold value.

Next, we study local as well as global stability of endemic equilibrium point E(x * , y * ) for R 0 > 1.

If R 0 > 1, then model system (2) has a unique endemic equilibrium which is locally asymptotically stable for β ≥ max abu 2 2 ,

Proof To prove the above-mentioned theorem, here we use Jacobian approach of matrix [8, 42] . The characteristic equation of the Jacobian matrix (11) at any endemic equilibrium point E(x, y) is given by

Therefore,

Since, from second equation of (5), we have

Now, we consider a function f (y) defined as

So, for endemic equilibria, y component is the solution of the

It follows from the expression of f (0) that

From (17) and (21), we can easily draw a relationship between F(0) and f (y), which gives

Thus, from (22), we see the following statements hold:

Here, we take R 0 > 1 this implies f (0) > 0 (by (20)).

We have already seen that if R 0 > 1, then there is only one

should decrease in some neighborhood of y * . So, in this case, we have f (y) < 0 and so F(0) > 0. Now, eigenvalues of the Jacobian matrix are solutions of the equation F(λ) = 0, i.e.,

where

It is clear from the expression of

. Therefore, all eigenvalues of Jacobian matrix have negative real part. So,

Theorem 8 For R 0 > 1, endemic equilibrium point is globally asymptotically stable if u 1 > abu 2 2 holds.

Proof To prove global stability of the endemic equilibrium point E(x * , y * ), here we use Dulac criterion [28, 39] . First, we write our model system (2) in the following form:

Now, we construct Dulac function as D(x, y) = 1 y , then we

Therefore by Dulac criterion for

In view of epidemiology, above result indicates disease persists in the population. Therefore for R 0 > 1, disease cannot be eradicated from the population when u 1 > abu 2 2 .

Bifurcation means qualitative change of the solution of a model system. In this section, we shall discuss different kinds of bifurcations of the solutions in the neighborhood of different equilibrium points. Among them, some are co-dimension 1 such as transcritical bifurcation, backward bifurcation, saddle-node bifurcation, Hopf bifurcation and the co-dimension two bifurcation likely the Bogdanov-Takens bifurcation of co-dimension 2.

First of all, we discuss the transcritical bifurcation for model (2) . Since for u 2 

is zero, and consequently, the usual eigenanalysis method fails here. Using Sotomayor theorem [39, 43] , we shall now establish the system experiences transcritical bifurcation when the control parameter (u 2 ) crosses the critical value u 2 = u 0 2 . Theorem 9 The model system (2) undergoes a transcritical bifurcation at disease-free equilibrium point (E 0 ) as the control parameter u 2 passes through the bifurcation value

Proof Here, we use Sotomayor's theorem [39, 43] to verify the transversality conditions for transcritical bifurcation; we

Then,

and therefore,

has one negative eigenvalue and other eigenvalue is zero. Let, V and W be the two eigenvectors corresponding to the zero eigenvalue of the matrix

Furthermore, differentiating partially (26) and (27) w.r.t. u 2 , we get

and

Thus, we have W T f u 2 (E 0 ; u 0 2 ) = 0,

Thus, all conditions of Sotomayor's theorem for transcritical bifurcation are satisfied; therefore, the model system (2) experiences transcritical bifurcation at E 0 when control parameter passes through the critical value u 2 = u 0 2 . Here, the stable disease-free equilibrium point become unstable and one stable endemic equilibrium point generates through transcritical bifurcation.

According to epidemiology, for R 0 < 1, disease can be eradicated from the system. But due to saturated treatment, the model system (2) may exhibit backward bifurcation. Epidemiologically the backward bifurcation is most important [28, 30] because in this case the system generates another stable endemic equilibrium point for R 0 < 1. In this situation, eradication of disease not only depends on the value of R 0 but it also depends on the initial number of infected populations. If system undergoes through backward bifurcation, we cannot conclude about eradication of disease for R 0 < 1. For backward bifurcation, there exists bistability of equilibrium points (one stable disease-free equilibrium and stable endemic equilibrium) for R 0 < 1. Therefore, if R 0 < 1, eradication of disease depends on initial infected population. Theorem 10 The model system (2) experiences a backward bifurcation at R 0 = 1 when (d + μ + δ + au 2 ) < Aabu 2 2 .

Proof In Sect. 3, we have already seen that infected component y satisfies Eq. (7). Since from the expression of R 0 , we have

so the coefficients of the above equation can be expressed as a function of R 0 putting value of β. So, the solution of the equation is a function of R 0 . Now, to derive a necessary and sufficient condition on the parameters such that the system experiences backward bifurcation with respect to R 0 , we have

Differentiating partially Eq. (4) w.r.to. R 0 and putting R 0 = 1, y = 0, we get

Thus, the endemic equilibrium curve bifurcates backward direction if the slope of the solution curve has negative value

the necessary and sufficient condition for backward bifurcation is (d + μ + δ + au 2 ) < Aabu 2 2 .

From above, it is clear that the model system (2) has backward bifurcation due to saturated treatment. So, there is a range R * 0 < R 0 < 1, for which model system (2) has two endemic equilibria. Now, we concentrate on discussing the stability of these two endemic equilibria when backward bifurcation condition is satisfied. So, we have the following theorem:

2 hold, then model system (2) has two endemic equilibria. The one with smaller number of infected (E 1 ) is unstable, while the other with higher number of infected (E 2 ) is locally asymptotically stable if β ≥ max abu 2 2 ,

Complex dynamics and control analysis of an epidemic model with non-monotone incidence and… Proof Since, the equation f (y) = 0 (where f (y) is given by (19) ) has two distinct real roots y * 1 , y * 2 with y * 1 < y * 2 , then f (y) is increasing in neighborhood of y * 1 that means f (y * 1 ) > 0 and consequently f (y) is decreasing in neighborhood of y * 2 that means f (y * 1 ) < 0. So, we have the following two cases.

Case-1 For the equilibrium E 1 (x * 1 , y * 1 ), we have f (y * 1 ) > 0, So F(0) < 0 (by (23)) . Since F(λ) is a quadratic equation, we have F(λ) → ∞ as λ → ∞ . Hence, there must be a real positive λ * such that F(λ * ) = 0. So Jacobian matrix has at least one positive eigenvalue. Therefore, (23)) . In this case, we proceed same as Theorem 7 and obtain that E 2 (x * 2 , y * 2 ) is locally asymptotically stable if β ≥ max abu 2 2 ,

Thus, we prove the above stated theorem.

In the next theorem, we shall establish the condition of saddle-node bifurcation. This type of bifurcation occurs when two fixed points create and annihilate depending on the model parameter.

Suppose at α = α [SN] , Eq. (7) has coincident real positive root, and the corresponding endemic equilibrium point is E(x * , y * ). Theorem 12 The model system (2) experiences saddle-node bifurcation at the coincident endemic equilibrium point E(x * , y * ) when the system parameter α crosses the critical value α = α [SN] and

Proof It can be easily shown that for α = α [SN ] , we have det(J (E))| α=α [SN] = 0 and tr(J (E))| α=α [SN] = 0 at the coincident endemic equilibrium point E(x * , y * ). Hence, one eigenvalue of the J (E) will be zero at α = α [SN] . If V and W are eigenvectors corresponding to eigenvalue zero for the matrices J (E) and J T (E), respectively, then

Furthermore, we have

Thus, the transversality condition for the occurrence of saddle-node bifurcation at coincident equilibrium E(x * , y * ) for α = α [SN] is satisfied for system (2) . Hence, by Sotomayor's theorem [39, 43] , we can conclude that system (2) has saddle-node bifurcation when model parameter α passes through its critical value α = α [SN] .

In this section, we consider the case when J (E(x * , y * )) has a pair of purely imaginary eigenvalues. Hopf bifurcation occurs when this pair of complex conjugate eigenvalues around non-trivial equilibrium E(x * , y * ) cross-imaginary axis [44] . In this case, endemic equilibria E(x * , y * ) loses stability through Hopf bifurcation under some parametric restriction. Considering 

, where ζ will be defined in the proof of the theorem.

when endemic equilibria E(x * , y * ) exists. That means the corresponding Jacobian matrix J (E) has a pair of purely imaginary eigenvalues if A = A [HB] . Now, in order to ensure the changes of stability through Hopf bifurcation, we have to check the transversality condition for Hopf bifurcation . The transversality condition is

Our next target is to find the direction of Hopf bifurcation in the neighborhood of the endemic equilibrium point E(x * , y * ) for the critical value A = A [HB] . First, we translate the endemic equilibrium E(x * , y * ) to origin using transformation X = x − x * , Y = y − y * , then model system (2) becomes

where the coefficients a i j is given in "Appendix-I." Let, denote the origin of the X -Y plane, then det(J ( )) = a 11 a 22 − a 12 a 21 > 0 and let Ω = √ det(J ( )). Then, under the transformation u = −X and v = a 11 X Ω + a 12 Y Ω , the model system (2) reduces to the normal form: The discriminating quantity (Γ ) for finding the direction of Hopf bifurcation is defined below:

where f uv denotes ∂ 2 f ∂u∂v (0, 0) and By [11, 45] , proof of the theorem is complete.

As so far, we have studied the bifurcations of co-dimension 1, i.e., considering only one parameter as the bifurcation parameter. The model system (2) experiences the bifurcations of co-dimension 1 such as transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation. We are now interested to investigate co-dimension 2 bifurcation such as Bogdanov-Takens bifurcation of co-dimension 2. Under the condition tr(J (E)) = 0 and det(J (E)) = 0, the Jacobian matrix of the coincident endemic equilibrium E(x * , y * ) has a double zero eigenvalue. Suppose these two conditions occur for the parameter values A = A 0 and α = α 0 when other parameters are keeping as constant. So the model system (2) may admit a Bogdanov-Takens (BT) bifurcation in vicinity of the coincident endemic equilibrium point. Now, first we investigate the conditions for which BT bifurcation will occur.

To ensure the existence of BT bifurcation for the critical values of the parameters ( i.e., A = A 0 and α = α 0 ), we are using the transformation X = x − x * , Y = y − y * and the condition tr(J (E)) = 0, det(J (E)) = 0, the model system (2) reduces to

where Q i (X , Y ), i = 1, 2 are smooth functions of (X , Y ) at least of order three and the coefficients are given in "Appendix-I" (using maple software [46] ).

We set x = X , y = a 11 X + a 12 Y , then system (30) transforms to

where R i (x, y), i = 1, 2 is a smooth function in (x, y) at least of order three and Next, letting x = x, s = y + c 3 y 2 and neglecting higher terms of order greater than 2 the system (31) becomes,

where R i (x, s), i = 3, 4 is a smooth function in (x, s) at least of order three. One more time, we set the variables x = x, z = s − c 6 xs, system (32) becomes

where R i (x, z), i = 5, 6 is a smooth function in (x, z) at least of order three. Finally, to obtain the normal form of Bogdanov-Takens bifurcation [47, 48] , we perform the trans-

where P i (u, v), i = 1, 2 is a smooth function in (u, v) at least of order three. From above, it is clear that c 4 = 0 and c 5 + 2 c 1 = 0 which implies model system (2) admits Bogdanov-Takens bifurcation. Thus, we obtain the following theorem:

Suppose that E(x * , y * ) is an endemic equilibrium of model system (2) such that tr(J E ) = 0 and det(J E ) = 0 with c 4 = 0, c 5 +2 c 1 = 0 . Then, E(x * , y * ) is a Bogdanov-Takens singularity of co-dimension 2 .

Now, we find analytical expressions for saddle-node bifurcation curve, Hopf bifurcation curve and homoclinic bifurcation curve in a small neighborhood of BT point, by determining the versal unfolding of the model system (2) as discussed in [49] [50] [51] . First, we give a small perturbation to bifurcation parameters A and α by A = (A 0 + λ 1 ) and α = (α 0 + λ 2 ) where |λ i | << 1, i = 1, 2 around their critical values then the model system (2) becomes:

Now translating endemic equilibrium point E(x * , y * ) to origin using the transformation X = x − x * , Y = y − y * , then system (35) becomes

where φ i (X , Y ), i = 1, 2 is a smooth function in (X , Y ) at least of order 3 and the coefficients b i are given in "Appendix-II" (using maple software [46] ).

Letting

where φ 3 (u, v, λ 1 , λ 2 ) is a smooth function in (u, v) at least of order 3 and the coefficients d j 's are given in "Appendix-III." Next, we introduce a new time variable τ by dt = (1 − d 5 u) dτ , and rewriting τ as t, we get from (37)

where φ 4 (X , Y , λ 1 , λ 2 ) is a smooth function in (X , Y ) at least of order 3 and e 1 = (

then the system (44) takes the form :

where φ 4 (x, y, λ 1 , λ 2 ) is a smooth function in (x, y) at least of order 3 and

Finally, we change variables using X = 

where φ 5 (X , Y , λ 1 , λ 2 ) is a smooth function in (X , Y ) at least of order 3 and ψ 1 (λ 1 ,

. If the transversality condition

holds, then the system (35) undergoes through Bogdanov-Takens bifurcation when (λ 1 , λ 2 ) is in the small neighborhood of (0, 0). The local representation of the bifurcation curves are given below:

(1) The saddle-node bifurcation curve is given by SN = {(λ 1 , λ 2 ): ψ 1 (λ 1 , λ 2 ) = 0, ψ 2 (λ 1 , λ 2 ) = 0}. (2) The Hopf bifurcation curve is given by H = {(λ 1 , λ 2 ):

= {(λ 1 , λ 2 ) : ψ 1 (λ 1 , λ 2 ) = − 49 25 ψ 2 2 (λ 1 , λ 2 ) + o(ψ 5 2 2 ), ψ 2 (λ 1 , λ 2 ) > 0}.

In this section, we shall justify the theoretical findings of our proposed model using some numerical approaches. For this purpose, first, we draw one parameter bifurcation diagram w.r.t. R 0 to justify the Backward bifurcation (see Fig. 2 ).

From Fig. 2 , it is clear that disease-free equilibrium point is stable if R 0 < 1, and it is unstable for R 0 > 1. Further, for R 0 < 1, there exists a range R * 0 < R 0 < 1 where two endemic equilibrium points exist. For R 0 < 1, the endemic equilibrium point with lower infected density is unstable and 

Backward bifurcation diagram with respect to R 0 and other parameters are given in Table 2 with A = 3.24, α = 0.35, the blue line corresponds to stable branch and red line is unstable the endemic equilibrium point with higher infected density is stable. Therefore for R 0 < 1, there exists bistability of equilibrium points (one stable disease-free equilibrium point and one stable endemic equilibrium point). At R 0 = 1, disease-free equilibrium changes its stability from stable to unstable with creation of one stable endemic equilibrium point; therefore, model system (2) undergoes through transcritical bifurcation at disease-free equilibrium for R 0 = 1. Next, we draw schematic bifurcation diagram for Bogdanov-Takens bifurcation (using [52] ) as given in Fig. 3 for model system (2) which divides first quadrant of A-α plane (considering other parameters values in Table 2 ) in six sub-regions, namely R 1 , R 2 , R 3 , R 4 , R 5 , R 6 . In Fig. 3 , the red line is saddle-node bifurcation (SN) line, black line represents Hopf bifurcation (H) line, cyan line indicates homoclinic bifurcation (HL) line, and magenta line denotes transcritical bifurcation line. Now, we draw the phase portraits for each of the sub-region of A-α plane considering different parameter values for interpreting number of equilibrium points and their behavior.

In region R 1 , there is one disease-free equilibrium (E 0 ) point which is stable node (see the Fig. 4a) . From epidemiological point of view, this region is most important because disease is not persisting here.

Next, we consider the parameter values on saddle-node line (the red line). For parametric values on the SN-line, along with the stable disease-free equilibrium point, one coincident endemic equilibrium point arises here which is unstable (saddle point) in nature (see Fig. 4b ). Since in this case, disease-free equilibrium is stable, disease eradicates from the population, and also endemic equilibrium is unstable; therefore, disease does not persist in the community. Now, crossing SN-line, we enter in region R 2 from region R 1 . In the region R 2 , there exists two endemic equilibrium points and the stable disease-free equilibrium point (E 0 ). Among endemic equilibrium points, one is stable spiral and another is saddle (see the phase portrait in Fig. 4c) . Therefore, in R 2 region, there exist bi-stability behavior of equilibrium points, one is disease-free equilibrium and other stable point is endemic equilibrium point. Stable manifold of saddle endemic equilibrium point (the green line in Fig. 4c ) divides the whole R 2 region into two sub-region. We see, if we take initial point inside the basin of attractor then solution of curve approaches to stable endemic equilibrium point, and if we take the initial point outside, the basin of attractor then solutions are converging to disease-free equilibrium point. Thus, we can conclude that for lower or much higher number of infected population at an initial stage in a system, disease will be eliminated and for a moderate number of infected population at an initial stage disease will persist in the system.

Next, we enter region R 3 by crossing Hopf bifurcation line (black line) from region R 2 . Here, we see that the number of equilibrium points is same as in region R 2 . But stable endemic equilibrium point becomes unstable and there arises a stable limit cycle. Here, also bistability arises: one is the stable disease free equilibrium and other is stable limit cycle around the endemic equilibrium point (see the phase portrait in Fig. 4d ). Stable manifold divides basin of attractor into two region, outside basin of attractor disease will eradicate and inside basin of attractor number of infected oscillates which means it is difficult to control the disease. That means for elimination of disease our aim is to minimize the density of infected individuals such that at any time it falls outside the basin of attraction. Now, we consider parameter values on the homoclinic line (cyan line). For these values of the parameters, the stable limit cycle connect the saddle endemic point, therefore a homoclinic loop appears (See the phase portrait in Fig. 5a) . Nature of other equilibrium points remain same. Thus, if initial density of populations are inside homoclinic loop, disease will persist in the system with oscillatory behavior and for lying outside homoclinic loop disease will eliminate from the system.

Next, we enter in region R 4 by crossing homoclinic bifurcation line from region R 3 . In region R 4 , stable limit cycle vanishes through homoclinic bifurcation, and it becomes unstable spiral. Only disease-free equilibrium is stable here. If we take initial point anywhere in region R 4 , solution of the system (2) converges to disease-free equilibrium, i.e., disease-free equilibrium is globally stable here (see the phase P. Saha, U. Ghosh Table 2 portrait in Fig. 5b ). That means disease will eliminate in the system if we choose the parameters in this region. Now, crossing transcritical bifurcation line (magenta line), we enter in region R 5 from region R 2 . Here saddle endemic equilibrium disappears, other endemic equilibrium point stay Table 2 stable and the disease-free equilibrium becomes unstable (saddle point). As there is no any other stable equilibrium, the endemic equilibrium point is globally stable here (see the phase portrait in Fig. 5c ). From biological point of view, this region is dangerous because disease permanently stay in the system.

At last, we enter in region R 6 from region R 5 by crossing Hopf bifurcation line (black line). Here, disease-free equilibrium is unstable (saddle point) and endemic equilibrium is unstable spiral. A stable limit cycle arises around endemic equilibrium (see the phase portrait in Fig. 5d ). Therefore, disease will persist here like the region R 5 but the density of the infected population will oscillate in this case.

By sensitivity analysis, we can identify the influential model parameters which have more impact on the basic reproduction number R 0 of the proposed model. The model parameters having higher sensitive index have more significant role in disease spreading dynamics among the population. By controlling such more sensitive parameters, the impact of the infection can be diminished or controlled. In order to identify such sensitive model parameters, we need to estimate the variation of basic reproduction number R 0 with respect to different model parameters. For calculating sensitivity index, we use normalized forward sensitivity index method [53] [54] [55] . According to normalized forward sensitivity index theory, sensitivity index of the basic reproduction number R 0 with

Positive sign of the sensitive index of the model parameters indicates basic reproduction number increases with model parameter increases and negative sign indicates basic reproduction number decreases with model parameter increases. In our proposed model most sensitive model parameters are susceptible recruitment rate(A), disease transmission rate (β), treatment control parameter (u 2 ), vaccinated control parameter (u 1 ), and cure rate (a). We enlist sensitivity indexes of all model parameters in Table 3 .

To perform sensitivity analysis numerically, we have used partial rank correlation coefficient (PRCC) method [56] . This 

Diagram of sensitivity index of each model parameter on R 0 using PRCC method correlates between output and input. Here, basic production number (R 0 ) is treated as output, and all model parameters are treated as inputs. PRCC measures robustness of sensitivity to determine the monotonic and nonlinear relationship between output and input; therefore, it plays a crucial role in sensitivity analysis. To carry out PRCC, we have taken sample size 100 and the corresponding PRCC diagram is placed in Fig. 6 . From Fig. 6 , it is clear that PRCC values lie from −1 to +1, which indicates that corresponding model parameters induce positive or negative impact on output R 0 . The PRCC values of model parameters imply how much model parameters are sensitive to R 0 and the sign assigns a relation between model parameters (input) and R 0 (output). Positive sign implies relation between input and output is proportional and negative sign implies relation is disproportional. It is obvious from Fig. 6 that model parameters A, β have positive influence on R 0 and model parameters d, u 1 , μ, δ, a, u 2 have negative influence on R 0 .

Our intention is to prevent spread of infection and eradicate disease from the population. For this purpose, we need to adopt suitable treatment policy and perfect vaccination campaign. In this paper, we consider two controls namely vaccination (u 1 ) and treatment (u 2 ). Now, our target is to identify best control among vaccination and treatment, i.e., we try to find out an optimal values of applied controls. Therefore, we have to use optimal control theory to get optimal system. To verify the existence of optimal control problem, we use Pontryagin's maximum principle. The main motive for applying optimal control is to decrease infected population, at the same time reduce applied control cost [32] .

For analyzing purpose, we consider controls as constant but implemented cost may be high for constant controls; therefore, we have to use control functions as function of time t.

In theory of optimal control, we consider time-dependent control functions. Now, we formulate an optimal control problem.

In this section, we formulate an optimal control problem of the proposed model. First, we consider an cost functional which will be minimum in the considered interval. Here, we construct integral of the cost functional as L(w, u) = A 1 x + u) is Lagrangian of optimal control problem. Therefore, optimal control problem is given by:

where

In L(t, w(t), u(t)), A 1 (> 0), A 2 (> 0) are loss due to susceptible and infected populations, respectively, whereas B 1 , B 2 are weight costs for control functions u 1 (t), u 2 (t), respectively. Here, we consider quadratic term in control to present severity of side effect of vaccination and treatment [57] . The control functions are Lebesgue measurable which are given by,

Theorem 15 For given an optimal control problem (43) subject to (1) admits an optimal pair u * 1 (t) and u * 2 (t) such that (44) .

Proof The state and control variables are non-empty and also nonnegative. Control constraint U is also convex.

Adding three equations of model (1), we have

, all state variables are bounded. Also integrand of objective functional J (u 1 (t), u 2 (t)) is convex w.r.t. u 1 and u 2 . Now, we rewrite system (1) as:

Thus, Lipschitz condition is satisfied for all state variables. Hence, there exists optimal control pair u * 1 (t) and u * 2 (t) such that

This completes proof of the theorem.

Theorem 16 Optimal control pair u * 1 (t), u * 2 (t) which minimizes J (u 1 (t), u 2 (t)) over U is given by

Proof For characterization of optimal control problem, we apply Pontryagin's maximum principle [13, 58] . First, we define Hamiltonian H as

with ρ i (T ) = 0, i = 1, 2, 3. Solving (45), we observe adjoint variables ρ 1 , ρ 2 , ρ 3 satisfy

with

For numerical simulation of optimal control problem, here we use forward-backward sweep method where in this method forward application of fourth-order Runge-kutta method of system (1) is combined with backward application of fourth-order Runge-kutta method of system (46) with transversality condition (47) . We take time interval as [0, 20], i.e., after 20 units time all controls will be stopped. To draw numerically, we use parameter values as in Table 2 with A = 100, α = 0.5 and A 1 = 0.05, A 2 = 0.05, B 1 = 0.6 and B 2 = 0.8 with initial conditions x(0) = 100, y(0) = 10, and z(t) = 0. In Fig. 7a -c, solution curves of susceptible, infected and recovered populations, respectively, are given with and without controls and time series of control u 1 , u 2 are shown in Fig. 7d and e. Our proposed model can be used for describing disease dynamics such as SARS, MERS, Dengue disease spreading. Since these types of disease spread rapidly in community at first, but after someday, the disease transmission starts to decrease; therefore, the consideration of such model is appropriate for such diseases.

Efficiency analysis is performed to find out best effective strategy when more than one control is applied. In this paper, we consider two controls namely vaccination for susceptibles (u 1 ) and treatment for infectives (u 2 ). Now, we have to know the effective strategy among applied controls to reduce the density of infectives and minimize the implemented control cost. So now, we perform efficiency analysis. Since two controls are used, therefore three cases may arise: (i) u 1 = 0, u 2 = 0 (ii) u 1 = 0, u 2 = 0 (iii) u 1 = 0, u 2 = 0. We call this cases as Strategy 1, Strategy 2 and Strategy 3, Table 4 .

From our analysis, Strategy 1 is the best strategy among three strategies that if we give both vaccine and treatment to population then disease will easily die out. But between vaccination and treatment, treatment is the more effective strategy, i.e., biologically if we use treatment to infected population, then density of infectives is more reduced.

In this paper, we have considered a deterministic epidemic model with constant birth rate and a non-monotone incidence rate affected by inhibitory factors, crowding effects, etc. Saturated type of treatment rate is taken into consideration for limited hospital resources. To discuss the complex dynamical behavior of the model system, we compute basis reproduction number R 0 , which has crucial role for eradicating the disease infection. We have established that when basic reproduction number is less than unity, then the disease will eradicate from the system for any initial population size, i.e., disease-free equilibrium is globally asymptotically stable under certain parametric restriction. Further we have seen that for R 0 > 1, endemic equilibrium point is globally asymptotically stable under some restrictions, i.e., disease will persist in the system. We use center manifold theorem to study the dynamical behavior of the disease-free equilibrium point at R 0 = 1. In our study, we observe that disease-free equilibrium point may be a saddle-node of co-dimension 1 or may be a semi-hyperbolic attracting node of co-dimension 2.

For discussing complex dynamical phenomenon, we examine the transcritical bifurcation, backward bifurcation, saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation of co-dimension 2. When backward bifurcation occurs, then the system shows bistability, one is stable disease-free equilibrium point and other is the stable endemic equilibrium point. By Hopf bifurcation, a stable limit cycle arises around the endemic equilibrium point. There exists a basin of attractor inside which all solution curves trend to that limit cycle and outside which all solution curves trend to disease-free equilibrium point which indicates disease oscillates inside the basin of attractor and outside of the basin of attractor disease eradicates from the community. Depending on the parameter value, the system experiences the homoclinic bifurcation, i.e., the stable limit cycle connecting the saddle point disappear through homoclinic bifurcation. Finally, numerical simulations are used to justify the theoretical findings. The numerical simulations show that there is a threshold range of the initial density of the infected population outside of which the disease will eradicate from the population.

Finally sensitivity analysis is performed to find out the influential model parameters which have most impact on the basic reproduction number R 0 . In order to identify most significant model parameters, we estimate variation of the basic reproduction number R 0 with respect to different model parameters. To control or eradicate the influence of the emerging disease, we need to take proper preventive measures to control the most sensitive model parameters. In our study, we observe that disease can be controlled or eradicated in some regions. Also there is some region where disease persists in the population. Optimal control is studied to reduce the density of infected population and minimize control cost. Using efficiency analysis, we identify most effective strategy. In our study, the most effective strategy is treatment strategy. Therefore, if we use treatment to infected population, then density of infected population is more reduced. In our view, study of controls may give some light in the layout.

Bifurcation analysis of an SIRS epidemic model with a generalized non-monotone and saturated incidence rate

Complex dynamics in a unified SIR and HIV disease model: a bifurcation theory approach

An attempt at a new analysis of the mortality caused by smallpox and of the advantages of inoculation to prevent it

The mathematical theory of infectious diseases and its applications

Infectious diseases of humans: dynamics and control

Mathematical models in population biology and epidemiology

A contribution to mathematical theory of epidemics

Dynamical modeling and analysis of epidemics

Compartmental models in epidemiology

A stochastic SIRS epidemic model with infectious force under intervention strategies

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

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

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

On the computation of R 0 and its role on global stability

An SIRS model with a non-linear incidence rate

Complex dynamics of an SIR epidemic model with saturated incidence rate and treatment

Dynamics of a SIR model with nonlinear incidence rate and treatment rate

Bifurcations of an SIRS model with generalized non-monotone incidence rate

A generalization of the Kermack-Mckendric deterministic epidemic model

Regulation and stability of hostparasite interactions: I. Regulatory processes

Global analysis of an epidemic model with non-monotone incidence rate

The effect of backward bifurcation in controlling measles transmission by vaccination

Bistability in a SIRS model with general non-monotone and saturated incidence rate

Bifurcation in an epidemic model with constant removal rate of the infectives

Backward bifurcation of an epidemic model with treatment

Dynamical behaviors of a modified SIR model in epidemic disease using non-linear incidence and recovery rates

Backward bifurcation of an epidemic model with saturated treatment function

Nilam: Dynamical model of epidemic along with time delay; Holling type II incidence rate and Monod-Haldane type treatment rate

Optimal control applied to biological model. Mathematical and compulational biology series

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 transmission

Cubic formula

Introduction to applied nonlinear dynamical system and chaos

Pliss and an introduction by V

Global stability of infectious disease models using Lyapunov functions

An introduction to mathematical epidemiology

Differential equations and dynamical systems

Hopf and backward bifurcations induced by immune effectors in a cancer oncolytic virotherapy dynamics

Hopf bifurcation of an epidemic model with a non-linear birth in population and vertical transmission

Practical computation of normal forms of the Bogdanov-Takens bifurcation

Practical computation of normal forms on center manifolds at degenerate Bogdanov-Takens bifurcations

Bifurcations of the limit circle of a family of plane vector fields

Versal deformations of a singular point on the plane in the case of zero eigen-values

Forced oscillation and bifurcation applications of global analysis

Matcont: A Matlab package for numerical bifurcation analysis of ODEs

Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model

Dynamics of SEIR model: a case study of COVID-19 in Italy

COVID-19 pandemic in India: a mathematical model study

A methodology for performing global uncertainty and sensitivity analysis in systems biology

Optimal control in epidemiology

Optimal vaccination strategy of a constrained time-varying SEIR epidemic model

Mathematical modeling of dengue epidemic: control methods and vaccination strategies

Authors have contributed equally in simulation, preparation of the article.Funding UGC, Govt. of India fellowship through UGCJRF scheme.Availability of data and material No specific data or unique material is used for this work.code availability Not applicable.

The authors of this paper declare no conflict of interest.