Complex dynamics and control analysis of an epidemic model with non-monotone incidence and saturated treatment.

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 R 0 . 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 R 0 = 1 . 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.

Introduction

Epidemiology is a part of mathematical biology which mainly associates with the study of infectious and vector-borne 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 the middle of twentieth century, mathematical epidemiology strode widely around the world [4-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-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 and standard incidence rate in the form 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 . The authors established that if < 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-21]. To avoid unboundedness of infection, Capasso and Serio [22], Anderson and May [23] introduced saturated types of infection rate in the form where denotes force of infection and 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 to include psychological effect of population. The function has maximum value at , it is increasing first upto , 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 into epidemic model wherewhere 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, bywhere 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 is not remain as a threshold quantity. The authors in [14, 27, 29, 30] considered treatment function asfor discussing the limited number of available medical facilities. Above function h(y) has maximum available limit 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 . In Sect. 4, we examine local as well as global stability of the disease-free equilibrium point and also endemic equilibrium point . 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 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.

Model formulation

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 with vaccination and saturated treatment

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 non-monotone 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 with nonnegative model parameters are enlisted in Table 1, and flow diagram of the proposed model is given in Fig. 1.

Table 1

Model parameters and their descriptions

Parameter	Description
A	Birth rate of the population
β	Transmission rate of the disease
α	Parameter measuring inhibitory factors
d	Natural death rate of population
δ	Disease-induced death rate
μ	Natural recovery rate of infectives
u_1	Vaccinated control parameter
u_2	Treatment control parameter
b	Delayed parameter of treatment
a	Cure rate

Fig. 1

Flow diagram of the proposed model

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 .

Theorem 1

All solutions (x(t), y(t)) of model (2) with initial conditions are nonnegative for all .

Proof

From first equation of model (2), we getUsing the integrating factor, we obtain the solution of the above differential equation in the form:i.e., for all . Similarly, we can show .

This completes the proof.

In the next theorem, we shall establish the boundedness of the solutions for all .

Theorem 2

Let, M= , then the region is a positively invariant and attracting set for the model (2).

First, we assume that .

Adding both equations in (2), we getwhere . Integrating and taking limit as , we get . Thus, we have for sufficiently large t. So, all solutions of model (2) are bounded and enter the region .

In the next section, we shall describe the conditions for existence of different equilibrium points.

Existence of equilibria and basic reproduction number

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 equilibriumAs the disease-free equilibrium exist, the basic reproduction number of the model (2) exists. In the following theorem, we compute basic reproduction number of the proposed model using next-generation matrix method [37].

Theorem 3

Basic reproduction number of the proposed model system (2) is given by

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 and with containing newly infected terms and the remaining terms, respectively. Therefore, and . Differentiating and with respect to infected component y at disease-free equilibrium point , we get F = and V= . Now, next-generation matrix is given by Basic reproduction number is the spectral radius of the next-generation matrix . Here,

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

The endemic equilibria of the considered model system is whereand is the positive root of the cubic equationwhere .

Discriminant [38] of cubic equation (7) isLetting and , the expression of () can be written in the form:whereHere, coefficients and are always positive and sign of depends on in (7). So, we have the following cases for positive roots of Q(y) = 0 :

Case-1 > 1, then 0; in this case, unique positive solution of exists, regardless of the sign of = .

Case-2 < 1, in this case, 0. Summarizing the above discussions, we get following theorem.

If = 0, the equation has no positive solution.

If = < 0, the following cases arise,

, then has two positive solutions.

, there is unique positive solution of the equation .

, we get no positive solution.

Theorem 4

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

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

When > 1, model system (2) has a unique endemic equilibrium.

When < 1 and

, model system (2) has no endemic equilibrium.

and also , there exist two endemic equilibria and .

For and also , these two endemic equilibria coalesce into the same endemic equilibrium .

Stability analysis of equilibria

For discussing stability of the equilibrium point, we need to compute Jacobian matrix of model system (2) at any equilibrium point E (, ). The Jacobian matrix of model (2 ) is given byNext, we shall prove local stability of disease-free equilibrium point for in Theorem 5. Then, we shall show global stability of disease-free equilibrium for using Lyapunov stability theorem in Theorem 6. Further, we shall discuss local stability of endemic equilibrium for in Theorem 7, and next we shall prove global stability of endemic equilibrium using Dulac function criterion for in Theorem 8.

First, we discuss local as well as global stability at disease-free equilibrium point .

Theorem 5

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

Jacobian matrix of the proposed model (2) at disease-free equilibrium point is given byTherefore, characteristic roots of are and . Therefore, for < 1, both eigenvalues are negative, and for > 1, one of them is positive. So, we can conclude that is locally asymptotically stable for < 1 and unstable for > 1.

For , one eigenvalue of the Jacobian matrix is zero; therefore, in this case, disease-free equilibrium point is non-hyperbolic type, and to examine dynamical behavior of the disease-free equilibrium point , we have to use center manifold theorem [39].

First, we use a perturbation at disease-free equilibrium point using the transformation , , thus model system (2) becomes:Next, we use the transformation , ; then, system (12) becomes:From (13), it is clear that if , then disease-free equilibrium point is a saddle-node of co-dimension 1 [39]. Now if , then we have to use center manifold theorem to find the nature of the disease-free equilibrium point .

By definition of local center manifold of the systemsatisfying where is sufficiently small, Dh represents the derivative of h with respect to u. Using the conditions, we get the value of as

Therefore, local manifold of the system is given by , and thus system (13) becomesFrom (15), it is clear that if , then disease-free equilibrium point is a semi-hyperbolic attracting node of co-dimension 2 for

Thus, proof of the theorem is complete.

Theorem 6

If < , then disease-free equilibrium point is globally stable where is defined in the text.

To show the global asymptotic stability of the disease-free equilibrium point , we are using here Lyapunov’s stability theorem [40, 41].

Let is a Lyapunov function as it is positive definite. Then,where .

Therefore, < 0 if < .

So, by Lyapunov’s stability theorem, it follows that is globally stable if < .

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 for .

Theorem 7

If > 1, then model system (2) has a unique endemic equilibrium which is locally asymptotically stable for

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 byTherefore,Since, from second equation of (5), we haveNow, we consider a function f(y) defined asSo, for endemic equilibria, y component is the solution of the equation with .

It follows from the expression of f(0) thatDifferentiating f(y) in (19) w.r.t. y, we getFrom (17) and (21), we can easily draw a relationship between F(0) and , which givesThus, from (22), we see the following statements hold:Here, we take this implies (by (20)). We have already seen that if , then there is only one endemic equilibrium . Since, hence f(y) should decrease in some neighborhood of . So, in this case, we have and so .

Now, eigenvalues of the Jacobian matrix are solutions of the equation , i.e.,whereIt is clear from the expression of that if .

Therefore, all eigenvalues of Jacobian matrix have negative real part. So, for > 1, we see that is locally asymptotically stable if

Theorem 8

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

To prove global stability of the endemic equilibrium point , here we use Dulac criterion [28, 39]. First, we write our model system (2) in the following form:Now, we construct Dulac function as then we haveTherefore by Dulac criterion for , endemic equilibrium point is globally asymptotically stable if .

In view of epidemiology, above result indicates disease persists in the population. Therefore for , disease cannot be eradicated from the population when .

Bifurcation analysis

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.

Transcritical bifurcation

First of all, we discuss the transcritical bifurcation for model (2). Since for = = , i.e., one eigenvalue of 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 () crosses the critical value = .

Theorem 9

The model system (2) undergoes a transcritical bifurcation at disease-free equilibrium point () as the control parameter passes through the bifurcation value .

Proof

Here, we use Sotomayor’s theorem [39, 43] to verify the transversality conditions for transcritical bifurcation; we denote whereThen,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 and , respectively, where V and W are given by and

Furthermore, differentiating partially (26) and (27) w.r.t. , we getTherefore, and Thus, we have ,if .

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

Backward bifurcation

According to epidemiology, for , 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 . In this situation, eradication of disease not only depends on the value of 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 For backward bifurcation, there exists bistability of equilibrium points (one stable disease-free equilibrium and stable endemic equilibrium) for . Therefore, if , eradication of disease depends on initial infected population.

Theorem 10

The model system (2) experiences a backward bifurcation at = 1 when .

In Sect. 3, we have already seen that infected component y satisfies Eq. (7). Since from the expression of , we haveso the coefficients of the above equation can be expressed as a function of putting value of . So, the solution of the equation is a function of . Now, to derive a necessary and sufficient condition on the parameters such that the system experiences backward bifurcation with respect to , we have to compute value of .

Differentiating partially Eq. (4) w.r.to. and putting , we getThus, the endemic equilibrium curve bifurcates backward direction if the slope of the solution curve has negative value at which is equivalent to . Thus, the necessary and sufficient condition for backward bifurcation is .

From above, it is clear that the model system (2) has backward bifurcation due to saturated treatment. So, there is a range < < 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:

Theorem 11

If and < hold, then model system (2) has two endemic equilibria. The one with smaller number of infected () is unstable, while the other with higher number of infected () is locally asymptotically stable if

Since, the equation (where f(y) is given by (19)) has two distinct real roots , with < , then f(y) is increasing in neighborhood of that means and consequently f(y) is decreasing in neighborhood of that means . So, we have the following two cases.

Case-1 For the equilibrium , we have > 0, So F(0) < 0 (by (23)) . Since is a quadratic equation, we have as . Hence, there must be a real positive such that F() = 0. So Jacobian matrix has at least one positive eigenvalue. Therefore, is unstable.

Case-2 For the equilibrium , we have < 0, So F(0) > 0 (by Eq. (23)) . In this case, we proceed same as Theorem 7 and obtain that is locally asymptotically stable if

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.

Saddle-node bifurcation

Suppose at , Eq. (7) has coincident real positive root, and the corresponding endemic equilibrium point is .

Theorem 12

The model system (2) experiences saddle-node bifurcation at the coincident endemic equilibrium point when the system parameter crosses the critical value and

It can be easily shown that for = , we have = 0 and at the coincident endemic equilibrium point . Hence, one eigenvalue of the J(E) will be zero at . If V and W are eigenvectors corresponding to eigenvalue zero for the matrices J(E) and , respectively, thenFurthermore, we haveand = [A ( + ) + ) ] . Thus, we haveandThus, the transversality condition for the occurrence of saddle-node bifurcation at coincident equilibrium for = 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 = .

Hopf bifurcation

In this section, we consider the case when has a pair of purely imaginary eigenvalues. Hopf bifurcation occurs when this pair of complex conjugate eigenvalues around non-trivial equilibrium cross-imaginary axis [44]. In this case, endemic equilibria loses stability through Hopf bifurcation under some parametric restriction. Considering A as the bifurcation parameter, suppose at , we have = 0 and 0. The existence of a Hopf bifurcation around endemic equilibrium is established in the following theorem:

Theorem 13

The model system (2) experiences a Hopf bifurcation at any endemic equilibrium point when the model parameter A crosses the critical value . Moreover,

If , there is a family of stable periodic orbits of model system (2) as A decreases from ,

If , there is a family of unstable periodic orbits of model system (2) as A increases from , where will be defined in the proof of the theorem.

From above discussions, it is clear that and at when endemic equilibria exists. That means the corresponding Jacobian matrix J(E) has a pair of purely imaginary eigenvalues if .

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 isOur next target is to find the direction of Hopf bifurcation in the neighborhood of the endemic equilibrium point for the critical value . First, we translate the endemic equilibrium to origin using transformation , , then model system (2) becomeswhere the coefficients is given in “Appendix-I.”

Let, denote the origin of the X–Y plane, then = > 0 and let = . Then, under the transformation and , the model system (2) reduces to the normal form:whereandThe discriminating quantity () for finding the direction of Hopf bifurcation is defined below:where denotes andBy [11, 45], proof of the theorem is complete.

Bogdanov–Takens bifurcation

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 and , the Jacobian matrix of the coincident endemic equilibrium has a double zero eigenvalue. Suppose these two conditions occur for the parameter values and 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.

Condition for existence of Bogdanov–Takens bifurcation

To ensure the existence of BT bifurcation for the critical values of the parameters ( i.e., and ), we are using the transformation , and the condition , , the model system (2) reduces towhere 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 , then system (30) transforms towhere is a smooth function in at least of order three andNext, letting and neglecting higher terms of order greater than 2 the system (31) becomes,where is a smooth function in (x, s) at least of order three. One more time, we set the variables system (32) becomeswhere 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 transformation , then system (33) transforms towhere , is a smooth function in (u, v) at least of order three. From above, it is clear that 0 and 0 which implies model system (2) admits Bogdanov–Takens bifurcation.

Thus, we obtain the following theorem:

Theorem 14

Suppose that is an endemic equilibrium of model system (2) such that = 0 and = 0 with 0, 0 . Then, is a Bogdanov–Takens singularity of co-dimension 2 .

Bogdanov–Takens bifurcation 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-51]. First, we give a small perturbation to bifurcation parameters A and by and where around their critical values then the model system (2) becomes:Now translating endemic equilibrium point to origin using the transformation , Y = , then system (35) becomeswhere is a smooth function in (X, Y) at least of order 3 and the coefficients are given in “Appendix-II” (using maple software [46] ).

Letting ,we getwhere is a smooth function in (u, v) at least of order 3 and the coefficients ’s are given in “Appendix-III.”

Next, we introduce a new time variable by , and rewriting as t, we get from (37)Let , then system (38) reformats as,where is a smooth function in (X, Y) at least of order 3 and = , = , = . Now, we set , then the system (44) takes the form :where is a smooth function in (x, y) at least of order 3 and = , Finally, we change variables using X = Y = , = and denoting by t, system (40) becomeswhere is a smooth function in (X , Y ) at least of order 3 and = ,

= . If the transversality conditionholds, then the system (35) undergoes through Bogdanov–Takens bifurcation when is in the small neighborhood of (0, 0) . The local representation of the bifurcation curves are given below:

The saddle-node bifurcation curve is given by : .

The Hopf bifurcation curve is given by : .

The Homoclinic bifurcation curve is given by :  = .

Numerical simulation

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. to justify the Backward bifurcation (see Fig. 2).

Fig. 2

Backward bifurcation diagram with respect to and other parameters are given in Table 2 with , , the blue line corresponds to stable branch and red line is unstable

From Fig. 2, it is clear that disease-free equilibrium point is stable if , and it is unstable for . Further, for , there exists a range where two endemic equilibrium points exist. For , the endemic equilibrium point with lower infected density is unstable and the endemic equilibrium point with higher infected density is stable. Therefore for , there exists bistability of equilibrium points (one stable disease-free equilibrium point and one stable endemic equilibrium point). At , 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 .

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 . 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.

Fig. 3

Schematic bifurcation diagram in A– plane fixing other variables as given in Table 2

Table 2

Model parameters and their respective values

4Parameter	β	d	δ	μ	u_1	u_2	b	a
Value	0.35	0.1	0.2	0.6	0.25	0.37	5.0	6.5

In region , there is one disease-free equilibrium 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.

Fig. 4

Phase portrait for the parameter values a In region (, ) b On Saddle-node line(, ) c In region (, ) d In region (, ) and other parameters are given in Table 2

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 from region . In the region , there exists two endemic equilibrium points and the stable disease-free equilibrium point (). Among endemic equilibrium points, one is stable spiral and another is saddle (see the phase portrait in Fig. 4c). Therefore, in 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 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 by crossing Hopf bifurcation line (black line) from region . Here, we see that the number of equilibrium points is same as in region . 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.

Fig. 5

Phase portrait for the parameter value a On Homoclinic line (, ) b In region (, ) c In region (, ) d In region (, ) and other parameters are given in Table 2

Next, we enter in region by crossing homoclinic bifurcation line from region . In region , 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 , solution of the system (2) converges to disease-free equilibrium, i.e., disease-free equilibrium is globally stable here (see the phase 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 from region . Here saddle endemic equilibrium disappears, other endemic equilibrium point stay 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 from region 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 but the density of the infected population will oscillate in this case.

Sensitivity analysis

By sensitivity analysis, we can identify the influential model parameters which have more impact on the basic reproduction number 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 with respect to different model parameters. For calculating sensitivity index, we use normalized forward sensitivity index method [53-55]. According to normalized forward sensitivity index theory, sensitivity index of the basic reproduction number with respect to model parameter is defined by 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 , vaccinated control parameter , and cure rate (a). We enlist sensitivity indexes of all model parameters in Table 3.

Table 3

Sensitivity indexes of the model parameters

Parameter	Sensitivity index
A	0.9999999999
β	1.0000000000
d	- 0.3159714719
u_1	- 0.7142857144
μ	- 0.1815431165
δ	- 0.06051437216
a	- 0.7276853253
u_2	- 0.7276853253

To perform sensitivity analysis numerically, we have used partial rank correlation coefficient (PRCC) method [56]. This method correlates between output and input. Here, basic production number () 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 to , which indicates that corresponding model parameters induce positive or negative impact on output . The PRCC values of model parameters imply how much model parameters are sensitive to and the sign assigns a relation between model parameters (input) and (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 and model parameters have negative influence on .

Fig. 6

Diagram of sensitivity index of each model parameter on using PRCC

Application of optimal control technique

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 and treatment . 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.

Formulation of 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) = , where L(w, u) is Lagrangian of optimal control problem. Therefore, optimal control problem is given by:whereIn L(t, w(t), u(t)), are loss due to susceptible and infected populations, respectively, whereas are weight costs for control functions 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,

Existence and uniqueness of optimal control problem

Theorem 15

For given an optimal control problem (43) subject to (1) admits an optimal pair and such that

where U is given by (44).

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

Adding three equations of model (1), we haveTherefore, , i.e., as , i.e., all state variables are bounded. Also integrand of objective functional is convex w.r.t. and .

Now, we rewrite system (1) as:

where

. Therefore,i.e., .

Now, if , then . Thus, Lipschitz condition is satisfied for all state variables. Hence, there exists optimal control pair and such that

This completes proof of the theorem.

Characterization of optimal control problem

Theorem 16

Optimal control pair which minimizes over U is given bywhere is the nonnegative root of .

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

Solving (45), we observe adjoint variables satisfywithAt , we have . Thus, we get , , where is the nonnegative root of .

Now, , and . at .

Hence, optimality condition is satisfied at .

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 , and , , and with initial conditions , , and . In Fig. 7a–c, solution curves of susceptible, infected and recovered populations, respectively, are given with and without controls and time series of control , 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.

Fig. 7

Time series of the population with control (blue line), without control (red line) and control variables; a Susceptible population b Infected population c Recovered population d e . (Color figure online)

Efficiency analysis

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 and treatment for infectives . 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) (ii) (iii) . We call this cases as Strategy 1, Strategy 2 and Strategy 3, respectively. Calculating efficiency index, we compare the applied strategies. The strategy with higher index is the best strategy [59]. Efficiency index (E. I.) is where and denote cumulative density of infected population without and with controls, respectively. We calculate value of and using Simpson’s rule. The value of is determined by . We get and values of and E.I. are given in Table 4.

Table 4

Values of efficiency indexes

Strategy	Applied controls	A^c	E.I.
Strategy 1	u_1≠0,u_2≠0	146.6684	0.7037
Strategy 2	u_1≠0,u_2=0	295.0079	0.4040
Strategy 3	u_1=0,u_2≠0	168.7416	0.6591

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.

Conclusion

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 , 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 , 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 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 In order to identify most significant model parameters, we estimate variation of the basic reproduction number 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.