Qualitative analysis and optimal control of an SIR model with logistic growth, non-monotonic incidence and saturated treatment

This paper describes an SIR model with logistic growth rate of susceptible population, non-monotonic incidence rate and saturated treatment rate. The existence and stability analysis of equilibria have been investigated. It has been shown that the disease free equilibrium point (DFE) is globally asymptotically stable if the basic reproduction number is less than unity and the transmission rate of infection less than some threshold. The system exhibits the transcritical bifurcation at DFE with respect to the cure rate. We have also found the condition for occurring the backward bifurcation, which implies the value of basic reproduction number less than unity is not enough to eradicate the disease. Stability or instability of different endemic equilibria has been shown analytically. The system also experiences the saddle-node and Hopf bifurcation. The existence of Bogdanov-Takens bifurcation (BT) of co-dimension 2 has been investigated which has also been shown through numerical simulations. Here we have used two control functions, one is vaccination control and other is treatment control. We have solved the optimal control problem both analytically and numerically. Finally, the efﬁciency analysis has been used to determine the best control strategy among vaccination and treatment.


Introduction
Mathematical modelling has been increasingly recognized as an important tool for understanding the transmission processes of different infectious diseases. A good mathematical model is able to evaluate the effective control strategies for the infectious diseases. It suggests effective control and preventive measures, and provides an estimate for the severity of the epidemic [1,2,3].
To formulate a compartmental model in mathematical epidemiology the researchers consider the constant growth rate, exponential growth rate, logistic growth rate of the population. The logistic growth rate is used for a relatively long-lasting disease or a disease with high death rate [17]. It is also considered due to limited space capacities or sources [21]. Again, the incidence rate plays a very important role in epidemiological modelling studies [22,25,26,27,28,29]. Authors use the bilinear incidence rate βSI (where the parameter β is transmission rate of infection and the variables S, I are, respectively, the number of susceptible and infected population) where the force of infection f (I) = βI is an increasing function of I and f (I) −→ ∞ as I −→ ∞. This is not realistic for many infectious diseases in large population because it does not include crowding effect of infected population and awareness factors of susceptible population [6,8,15]. Capasso and Serio in [4] proposed a saturated incidence rate βSI 1+αI , where α is defined as inhibitory or awareness factors. Here the force of infection f (I) = βI 1+αI is an increasing function of I and f (0) = 0. Also, f (I) −→ β α as I −→ ∞. Again, the force of infection decreases as the parameter α increases. Here the effect of α stems from epidemic control (taking awareness or protection measures) [22]. Xiao and Ruan [19] proposed a non-monotonic incidence rate with psychological effect of the form βSI 1+αI 2 , which includes the effects of psychological factors, protection measures and intervention policies when a serious disease emerges.
Here the infection force f (I) = βI α . This non-monotonic incidence rate is very effective for a new disease in some countries or regions, where initially the contact rate and the infection probability increase as people have poor knowledge about the disease but when I becomes large, the disease becomes more serious [28,29]. Then people are aware about the disease, so they take appropriate preventive measures and awareness and so infection force will then decrease.
To recovery from infection, treatment of infected population is the most important method. In classical epidemic models, the authors use the treatment function as T (I) = rI, I ≥ 0 (where r is a positive constant). If the number of infected population is very large then it is not always possible to provide such type of treatment which is proportional to the number of infected population. To avoid this Wang and Ruan [5]  , which is non-constant and bounded above by the upper bound kI 0 . Recently, Zhang and Liu [6] introduced a continuously differentiable treatment function of the form T (I) = aI 1+bI (where a, b are positive constants), which is clearly an increasing function of I and is bounded above by the least upper bound a b . This saturated treatment function is a better alternative for outbreak disease such as SARS, Dengue, etc in a new region or area [18], because at the beginning of the outbreak there is small of effective treatment due to negligence or lack of knowledge about the disease or lack of awareness. But then the treatment rate is increased with the improving of hospital's treatment conditions including skilful treatment techniques and effective medicines [22,26]. Finally, the treatment rate is reached to its maximum due to boundedness of medical resources of any country or community. Here the parameter b measures the extent of the effect of the infected being delayed for treatment [6].
On the other hand, optimal control theory is applied to the epidemiological models to determine the control strategies [7,12,22,26]. Since economical resources of any country or community are limited, hence the main objective for applying optimal control in epidemiological models is to minimize total loss occurs due to the presence of infection as well as the total cost due to implementing of the control(s). Generally, researchers use treatment or vaccination or both controls in SIR,SIRS,SEIR models [22,26]. Optimal control theory helps to address the question of how to optimally combine the control strategies for minimizing the infection in a community with limited resources.
In this article, we have considered an SIR model with logistic growth rate of susceptible population, a non-monotonic incidence rate of the form βSI 1+αI 2 and a saturated treatment function of the form T (I) = au 2 I 1+bu 2 I . We have also considered two controls, one is vaccination control u 1 and other is treatment control u 2 . These two controls have been used to address the question of how to optimally combine the vaccination and treatment strategies to minimize the susceptible and infected population as well as the cost of implementation of these two interventions. In our work, we have studied the stability and bifurcation analysis of co-dimension one and two in the neighbourhood of equilibrium points. It is important to note that we shall deal with not only qualitative analysis of the model but also the optimal control of the disease. Efficiency analysis shall be performed to determine the best control strategy. It should be mentioned that the stability or instability of endemic equilibrium points shall be analysed by applying different technique described in [9].
Let, S(t), I(t) and R(t) be, respectively, the abundance of susceptible, infected and recovered population at time t.
Incorporating all the assumptions, our proposed SIR model can be formulated as with initial conditions S(0) ≥ 0, I(0) ≥ 0, R(0) ≥ 0. Here, the intrinsic growth rate r should be greater than u 1 , otherwise there is no biological sense. All the model parameters used in the system (1) are non-negative and interpreted in Table 1. b : Delayed parameter of treatment.
The control variable, be the percentage of susceptible individuals being vaccinated per unit of time.
The treatment control parameter.
The organization of this paper is as follows. Boundedness of solutions is given in Section 2. Section 3 is devoted to the equilibria, the basic reproduction number, stability and bifurcation analysis for the model. The optimal control problem and the efficiency analysis are discussed in Section 4. Final section gives the conclusions.

Boundedness of Solutions
Before the investigation of boundedness of solutions of the system (1) first we study the positivity of the solutions of the considered system. The positivity is important for feasibility of the solutions and boundedness implies the finiteness of the the solutions.
From the first equation of (1), we have Integrating the above inequality and using initial conditions, we obtain Similarly, from other equations of (1) we can obtain the following It is clear from the above expressions that all the solutions of the system (1) are feasible. Now, we analyse the boundedness of solutions of the system (1). To show the boundedness of solutions of the system (1), first we add all three equations of the system (1) and we have , which implies that for any real number λ is a positive number. Integrating and taking limsup as t → ∞, we get lim sup t→∞ N (t) ≤ M λ . Thus, we can summarize the details in the following: Since R does not involve in first two equations of (1), hence the system (1) is equivalent to investigate the system (2). Since the exact solution of the non-linear autonomous system (2) is impossible to find, hence we shall analyse the qualitative behaviour of the solutions in the neighbourhood of the equilibrium points of the system (2).

Equilibria, Basic Reproduction Number, Stability and Bifurcation Analysis
In this section, we shall investigate the existence of equilibria and basic reproduction number of the model (2). The stability analysis of equilibrium points will be discussed here. Different bifurcation analysis, namely transcritical bifurcation, backward bifurcation, saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation will be analysed. Throughout this section we assume that two controls u 1 and u 2 are fixed constants.

Equilibria and Basic Reproduction Number
The system (2) has always the trivial equilibrium point A 0 (0, 0). Since r > u 1 , hence the disease free equilibrium r always exists. Since the system (2) has a DFE, hence the system has a threshold parameter R 0 , known as the basic reproduction number.
Proof. : The system (2) has only one infected compartment, that is, the variable I and A 1 (S 1 , 0) is the disease free equilibrium point. The basic reproduction number R 0 is defined as the spectral radius of the next generation matrix F V −1 with large domain [10,24], where F = . Thus, the basic reproduction number R 0 of the model (2) is βk(r−u 1 ) r(d+δ+γ+au 2 ) . Hence the lemma is proved.

Stability Analysis of Trivial and Disease Free Equilibrium Point
To investigate the stability analysis of equilibria, we use variational matrix of the system (2), which is Theorem 1. The trivial equilibrium point A 0 is a saddle point.
Proof. : The eigenvalues of the variational matrix J at the equilibrium point A 0 are (r −u 1 ) and −(d+δ +γ +au 2 ).
Since r > u 1 , hence the trivial equilibrium point A 0 is a saddle point. Hence the theorem is proved.
Theorem 2. If R 0 < 1 then DFE A 1 is asymptotically stable and if R 0 > 1 then it is unstable.
Proof. : The eigenvalues of the variational matrix J at the equilibrium point A 1 are −(r − u 1 ) and (d + δ + γ + Hence the theorem is proved. Proof. : For R 0 = 1, the eigen values of the variational matrix (4) at the equilibrium point A 1 (S 1 , 0) are 0, −(r − u 1 ). So, A 1 is a non-hyperbolic equilibrium point and Centre Manifold Theory [11] will be applied to determine its stability.
Putting S = S − S 1 , I = I in the system (2) and using Taylor's expansion we get (omitting the 'dash' sign) (Neglecting the terms of order ≥ 3) . By using the transforma- (5) can be transformed into the form (omitting the 'dash' sign) where By the Centre Manifold Theory, there exists a centre manifold of the system (6) which can be expressed by Theorem we have the flow on the centre manifold W c (0) defined by the differential equation Since βh 1 = − βr k(r−u 1 ) < 0, hence A 1 is asymptotically stable. Hence the theorem is proved.
) then DFE A 1 is globally asymptotically stable.
Proof. : Here we consider the Dulac function B(S, I) = 1+bu 2 I SI and we get Again, Hence the theorem is proved.
This result is biologically very significant, because DFE is globally asymptotically stable i.e. the disease will be eliminated if the transmission rate of infection does not exceed a threshold.

Transcritical Bifurcation
Theorem 5. If βkr > βku 1 + (d + δ + γ)r, then the system (2) experiences a transcritical bifurcation at DFE A 1 as the cure rate a varies through the bifurcation value a 0 = βkr−{βku 1 +(d+δ+γ)r} has a simple eigenvalue λ = 0. Thus, we shall use Sotomayor theorem [11] to establish the existence of transcritical bifurcation. Now,  and a eigenvector of Therefore, all the conditions for transcritical bifurcation in Sotomayor theorem are satisfied. Hence, the system (2) experiences a transcritical bifurcation at the equilibrium point A 1 as the parameter a varies through the bifurcation value a = a 0 . Hence the theorem is proved.
Thus, the unstable DFE becomes stable as the cure rate crosses the bifurcation value. So, cure rate plays a crucial role in the changing of dynamics of infectious disease.

Backward Bifurcation and Stability Analysis of Endemic Equilibria
In this section, we shall establish the existence of backward bifurcation and the conditions for stability of endemic equilibria of the system (2).
Theorem 6. The system (2) experiences a backward bifurcation at R 0 = 1 if and only if bru Proof. : We have already established that the infected component I of any endemic equilibrium point (S, I) of the system (2) satisfies the Eq. (3). Differentiating Eq. (3) implicitly with respect to R 0 we get The system (2) has a backward bifurcation at R 0 = 1 if and only if the slope of the curve I = I(R 0 ) at the point (1, 0) in R 0 − I plane must be negative. Thus, we obtain a necessary and sufficient condition that the backward bifurcation occurs in the form bru 2 (d + δ + γ) + β 2 k < βkbu 2 (r − u 1 ).
Hence the theorem is proved.
Here it is important to note that if b = 0 there can not be a backward bifurcation i.e. the phenomenon of the backward bifurcation occurs due to the saturated treatment function. In case of backward bifurcation, there exists a positive number R * 0 < 1 such that a stable endemic equilibrium point arises with the stable DFE for R * 0 < R 0 < 1 i.e. the bi-stability arises for R * 0 < R 0 < 1, which has been shown by the backward bifurcation diagram (see Figure  1). In Figure 1, red and blue lines denote the lines of unstable and stable equilibrium points, respectively. Now, the stability or instability of endemic equilibria shall be analytically established by the following theorem.
Proof. The characteristic equation of the variational matrix J at any endemic equilibrium point (S, I) of the system (2) is given by So, Again, the component I of any endemic equilibrium point (S, I) satisfies the equation Thus, R 0 > 1 if and only if H(0) > 0 and R 0 < 1 if and only if H(0) < 0.
The numerical verification of the above theorem has been presented in Figure 1.

Saddle-node and Hopf Bifurcation
The number and nature of the equilibrium points of the system (2) depend on the model parameters. In next two theorems we shall show that the system (2) experiences the saddle-node as well as the Hopf bifurcation with respect to some model parameter.
Let f α denote the vector of partial derivatives of the components of f with respect to α.
Therefore, by Sotomayor theorem we can say that the system (2) experiences a saddle-node bifurcation at the endemic equilibrium point A 2 (S * , I * ) when the parameter α passes through α = α [SN ] .
Hence the theorem is proved.
Theorem 9. The system (2) exhibits a Hopf bifurcation leading to a family of periodic solutions that bifurcates from the endemic equilibrium point A 2 (S * , I * ) for suitable values of intrinsic growth rate r in a neighbourhood of r [HB] .
In addition, if Γ < 0 then the system is said to be supercritical and if Γ > 0 then the system is said to be subcritical, where Γ is the Liapunov number.
Therefore, we can easily compute the Liapunov number Γ, which is as follows:

By using Hopf bifurcation theorem, we obtain that Hopf bifurcation is supercritical if Γ < 0 and subcritical if Γ > 0.
Hence the theorem is proved.

Bogdanov-Takens Bifurcation
We have already analysed the co-dimension 1 bifurcations considering α (for saddle-node bifurcation) and r (for Hopf bifurcation) as bifurcation parameter for the system (2). Here, we shall analyse a co-dimension 2 bifurcation considering α and r as the bifurcation parameter for the system (2).
Since the matrix A 0 has two zero eigenvalues, hence a 11 + a 22 = 0 and a 11 a 22 = a 12 a 21 .
where SN , H, and HL are the saddle-node, the Hopf and the homoclinic bifurcation curve, respectively.
Hence the theorem is proved.
To interpret different bifurcations numerically, we have considered r and α as the bifurcation parameters. In Figure 2, we have presented the schematic bifurcation diagram in r −α plane of system (2) with different bifurcation curves. In this figure, the points BT and GH, respectively, denote the Bogdanov-Takens bifurcation point and global Hopf bifurcation point. By different bifurcation curves, the feasible region of r − α plane is divided into six subregions in the basis of existence and character of equilibrium points and limit cycle. First we consider the region R 1 , where the trivial equilibrium point is unstable and DFE is stable and no endemic equilibrium point exists (see Figure 3(a)). This region is biologically significant because disease will not persist for the parametric values in the region R 1 . Thus, if the intrinsic growth rate of the population is low then disease will be eliminated easily. Now, we consider the parametric values on the saddle-node line. Then only one endemic equilibrium point arises with unstable mode along with a unstable trivial equilibrium point and a stable DFE (see Figure 3(b)). Now crossing the saddle-node line, we enter the region R 2 where two endemic equilibrium points arise along with unstable trivial equilibrium point and stable DFE (see Figure 3(c)). Among two endemic equilibria, one is saddle point and other is stable focus. Thus, bi-stability arises. Now crossing the Hopf bifurcation line, we enter the region R 3 where the number of equilibrium points is same as R 2 . But, the stable focus becomes unstable focus as we enter R 3 from R 2     and a stable limit cycle arises. Thus, here two attractors exist, one is stable DFE and other is stable limit cycle (see Figure 3(d)). The stable manifold of the saddle endemic equilibrium point separates the basins for two attractors, so it is called the basin boundary (green line). Two trajectories that form the stable manifold of the saddle point are also called separatrices. If the initial point is outside the separatrices then the disease will be eradicated and if the initial point is inside the separatrices then the number of infected population oscillates and so the disease is difficult to control. Thus, for the parametric values in region R 3 , the initial values of susceptible and infected population play a crucial role to control the disease. Now, we consider the set of values of r and α on homoclinic bifurcation line then the stable limit cycle moves closer and closer to the saddle point and it touches the saddle point and a homoclinic loop arises (see Figure 4(a)). Now crossing the homoclinic bifurcation line, we enter the region R 4 then the number and nature of equilibrium points is same as R 3 but the stable limit cycle disappears through homoclinic bifurcation and in this case the only attractor is DFE and so it is globally asymptotically stable (see Figure 4(b)). Now we move from R 2 to R 6 through the transcritical bifurcation line, then the saddle endemic equilibrium point disappears and the DFE becomes unstable (see Figure 4(d)). Thus for these values of parameters in R 6 , only a stable endemic equilibrium point exists and other points are unstable. So the stable endemic equilibrium point is globally asymptotically stable i.e. disease always persists. Now we enter from R 6 to R 5 by crossing the Hopf bifurcation line, then the endemic equilibrium point becomes unstable focus and a stable limit cycle arises (see Figure 4(c)). So, in this case there is only one attractor which is stable limit cycle. In Figure 2, there exists values of the parameters in the vicinity of GH point for which two limit cycles arise around the unstable endemic equilibrium point. Among these two limit cycles the smaller one is stable and larger one is unstable. The corresponding phase portrait is given in Figure 4(e).

Optimal Control Problem and Efficiency Analysis
In this section, we shall investigate the optimal control of the model (1) on the assumption that the vaccination control u 1 and the treatment control u 2 are the functions of time t. Here, we shall also perform an efficiency analysis to determine the best control strategy between vaccination of susceptible population and treatment of infected population.

Optimal Control Problem
The system (1) is reformulated as an optimal control problem, where both vaccination and treatment controls are time dependent as they are applied according to the necessity. Here, we construct the objective functional as follows where the constants A 1 and A 2 are, respectively, per capita loss due to presence of susceptible and infected individual. The constants B 1 and B 2 represent the costs associated with vaccination of susceptible population and treatment of infected population, respectively. Our goal is to minimize the total number of susceptible and infected individuals as well as the costs associated with the implement of vaccination and treatment controls on the time interval [0, T ]. The optimal control problem (1) is to find optimal functions (u * 1 (t), u * 2 (t)) in such a manner that is Lebesgue measurable on [0, 1] and 0 ≤ u 1 (t), u 2 (t) ≤ 1 for all t ∈ [0, T ]}. Now, we show the existence of an optimal control for the system (1).
Proof. : The integrand of the objective functional J(u 1 , u 2 ) is a convex function of u 1 and u 2 , because the constants B 1 and B 2 are positive. Moreover, the control space U is also closed and convex region. Hence the optimal control is bounded and therefore there exists an optimal pair (u * 1 , u * 2 ) which minimizes J for t ∈ [0, T ] with the help of the system of differential equation (1) [25,26]. Hence the theorem is proved. Now applying Pontryagin's Maximum Principle [13], we convert (1) and the objective cost functional J(u 1 , u 2 ) into a problem of minimizing a Hamiltonian H, which is given by where λ 1 , λ 2 and λ 3 satisfy the adjoint equations dλ 1 (t) = − ∂H ∂R with the transversality conditions λ i (T ) = 0, i = 1, 2, 3 i.e. λ i (i = 1, 2, 3) satisfy the system of equations with the transversality conditions Now, we differentiate the Hamiltonian H partially w.r.t. u 1 , u 2 and we obtain the optimality conditions that follows and 2B 2 u 2 (1 + bu 2 I) 2 = (λ 2 − λ 3 )aI. From these two, we obtain the optimal pair (u * 1 , u * 2 ) as stated below γ 0.7 [22] u * 2 = max 0, min u 2 , 1 , where u 2 is the non-negative root of the equation Here, S * , I * , R * are, respectively, the optimum values of S, I, R and (λ * 1 , λ * 2 , λ * 3 ) is the solution of the system (19) with the condition (20). Thus, we summarize the details in the following: Theorem 12. The optimal pair (u * 1 , u * 2 ) that minimizes J over the region U are given by u * 1 = max 0, min and u * 2 = max 0, min u 2 , 1 , where u 2 is the non-negative root of the equation 2B 2 u 2 (1 + bu 2 I * ) 2 = aI * (λ * 2 − λ * 3 ).
To justify the theoretical findings of the optimal control problem (1), we have solved it numerically by applying forward-backward sweep method that combines the forward application of a fourth order Runge-Kutta method for the state system (1) with the backward application of a fourth order Runge-Kutta method for the adjoint system (19) with the transversality conditions (20). Here, we assume T = 20 and so vaccination and treatment are stopped after 20 units of time. We also assume that the parametric values and values of weight constants A i , B i (i = 1, 2) for simulation of optimal control problem are given in Table 2 with the initial conditions S(0) = 50, I(0) = 4 and R(0) = 0.01.
In Figure 5(a)-(c) we can compare the susceptible (S), infected (I) and recovered (R) population at any time 20] for no control and with controls. Figure 6(a) and 6(b) represent the optimal control functions u * 1 and u * 2 . Thus, vaccination and treatment are highly effective for reducing both susceptible and infected population and application of these two controls gives more number of recovered population than the no control system.

Efficiency Analysis
When two or more controls are used in a optimal control problem, then the efficiency analysis is applied to compare the efficiency of different control strategies for reducing the infection of any disease. The efficiency index (E.I.) [14] of any control strategy is defined as E.I. = (1 − A c A o ) × 100, where A c is the cumulated number of infected individuals when the control strategy is applied and A o is the cumulated number of infected individuals without use of the strategy. The best strategy will be the one whom efficiency index will be bigger. In this paper, two control functions, namely vaccination control u 1 and treatment control u 2 , are considered. Here, we distinguish two control strategies Strategy1 and Strategy2 where Strategy1 is the strategy where only vaccination is used (i.e.u 1 = 0 , u 2 = 0) and Strategy2 is the strategy where only treatment is used (i.e.u 1 = 0 , u 2 = 0). To determine the best control strategy among these two, we have to calculate E.I. for each strategy. It is noted that the cumulated number of infected individuals during the time interval [0, 20] is defined as A = 20 0 I(t)dt and Simpson's 1 3 rule has been applied to evaluate the value of integration. Then, we have A 0 = 98.8610. The values of A c and efficiency index (E.I.) for Strategy1 and Strategy2 have been given in Table 3. From Table 3, we can conclude that Strategy1 is better than Strategy2. Epidemiologically, for the parametric values given in Table   2 the vaccination control is more effective than treatment to reduce infection.

Conclusions
This paper deals with an SIR model with logistic growth rate of susceptible population for a long-lasting disease or a disease with high death rate and a non-monotonic saturated incidence rate has been considered which includes psychological effects, protection measures and intervention policies when a serious disease emerges. A saturated treatment rate has also been used to represent the boundedness of medical resources of any country or community.
The DFE is locally asymptotically stable if the basic reproduction number R 0 is less than 1 and is unstable if R 0 > 1. If R 0 = 1, the DFE is non-hyperbolic equilibrium point and it is shown that it is stable by using Centre Manifold Theory. It is also shown that DFE is globally asymptotically stable if R 0 < 1 and the transmission rate of infection less than some quantity, which signifies that the disease will die out for low transmission rate along with basic reproduction number R 0 < 1. The system undergoes a transcritical bifurcation at DFE w.r.t. cure rate, which indicates the cure rate plays an important role to eliminate the infection of a disease. We have obtained a necessary and sufficient condition that the backward bifurcation occurs and it has been proved that the backward bifurcation occurs due to saturated treatment rate. The phenomenon of backward bifurcation is biologically important because making the basic reproduction number less than unity is not sufficient to eradicate the disease if the condition for backward bifurcation holds. Here, it is important to mention that we have analytically established in Theorem 7 the stability or instability of different endemic equilibria as the basic reproduction number R 0 varies. The system experiences saddle-node and Hopf bifurcation of co-dimension 1 at endemic equilibrium point w.r.t. the inhibitory factor α and intrinsic growth rate r, respectively. The system also undergoes a BT bifurcation of co-dimension 2 at endemic equilibrium point w.r.t. α and r.
We have also made the SIR model to an optimal control problem by considering vaccination control u 1 and treatment control u 2 . The optimal control to minimize the susceptible, infected individuals and costs for implementation of these two controls has been obtained and numerical results show the positive impacts for implementing vaccination to susceptible individuals and treatment for infected individuals. We have also done an efficiency analysis, which clears that the vaccination control is more effective than the treatment control to control the disease. This work is theoretical modelling and it can be further justified by using experimental results.

Conflict of Interest
Conflict of Interest: The authors declare that they have no conflict of interest.