Stability analysis and optimal control of covid-19 with convex incidence rate in Khyber Pakhtunkhawa (Pakistan).

The dynamic of covid-19 epidemic model with a convex incidence rate is studied in this article. First, we formulate the model without control and study all the basic properties and results including local and global stability. We show the global stability of disease free equilibrium using the method of Lyapunov function theory while for disease endemic, we use the method of geometrical approach. Furthermore, we develop a model with suitable optimal control strategies. Our aim is to minimize the infection in the host population. In order to do this, we use two control variables. Moreover, sensitivity analysis complemented by simulations are performed to determine how changes in parameters affect the dynamical behavior of the system. Taking into account the central manifold theory the bifurcation analysis is also incorporated. The numerical simulations are performed in order to show the feasibility of the control strategy and effectiveness of the theoretical results.

Introduction

In December 2019, the capital of Hubei Province i.e., Wuhan city gains attention worldwide, at the end of 2019 as an unknown virus attack people and start killing. Later on, the causative agent was identified as a novel coronavirus and now named as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The Chinese government did its best to minimize the spread but unluckily the virus spreads to other countries. The latency period of covid-19 is from 11 to 14 days whereas myalgia, diarrhea, cough, fever, and shortness of breath are common symptoms. Aged people are highly exposed to covid-19 as well as patients having comorbid conditions. The covid-19 transmission rate is very high and the basic reproduction number lies between 2.2 and 3.58. That’s why it spread throughout the world and affected 213 countries. Hence, the world health organization (WHO) declared covid-19 as one of the global pandemics on January 30, 2020. The two main routes of transmission is direct physical contact with the infected person and respiratory droplet coming out from covid-19 patient [1].

At first, Iran and Italy were the highly exposed countries to covid-19 after China. Iran is the neighboring country of Pakistan and every year thousands of people came from Pakistan to Iran for religious ceremonies. While coming back from Iran some of these people were infected and they become the cause for spreading covid-19 in Pakistan. Although the government closed the Pakistan-Iran border but yet the first case was confirmed officially on February 26, 2020, from Karachi. Avoiding the spread and causalities government decided to quarantine most of the people at home [1].

Mathematical modeling is considered an effective tool for describing the dynamical behavior of infections [2], [3], [4], [5], [6]. Mathematicians frequently using the tools of mathematical modeling from the last century. More recently, the field of mathematical modeling got considerable attention and a number of authors put their valuable contributions in this area. For realizing and controlling the outbreak of transmissible diseases in a group, many researchers have formulated models. The second leading source of death around the globe is infectious diseases. The application of mathematical modeling has been in vogue for the study of transmissible infectious diseases. Many scholars have been discussing over the last few years infectious diseases and their dynamics using various approaches. “To study the dynamics of various infectious diseases, mathematical modeling is considered as one of the best techniques to formulate the phenomenon in the system of equations. Several researchers have worked on different infectious diseases. They have developed different mathematical models for these epidemic diseases, and then studied the stability analysis and optimal control of these epidemic models (see e.g., [7], [8], [9], [10], [11], [12]), which not only helps in the control/spread of infectious diseases but also helps in prevention of these diseases in daily life. Many researchers have worked on epidemic models to analyze and control different diseases, for example, hepatitis B, avian influenza, leishmaniasis, tuberculosis, etc. Modeling of epidemic models is helpful to academia as well as to daily life. Mathematical modeling has a rich literature on the transmission dynamics of infectious diseases. Mathematical modeling is one of the powerful tools to describe the dynamical behavior of different diseases [22], [23], [24], [25], [26], [27].

The incidence rate plays a major role in epidemic models. The incidence rate of was introduced by Yorke and London [13] whereas with was investigated by Liu and coworkers [14]. While this incidence rate is more difficult to deal with than the bilinear one, it seems more reasonable. One of the interesting examples of nonlinear incidence rate is convex incidence rate i.e , where and are positive constants. For a short time period it corresponds to an increased rate of infection which causes due to two exposure. The contact rate results from a single contact, while the term corresponds to a new infective which comes into being due to the double exposure. Until now researchers established stability analysis of epidemic models by geometric approach with non-convex (respect to I) incidence rate. For convex incidence of the approach seems to give very restrictive conditions [33].

One of the reasons why we chose the convex incidence rate is because this function could model the increase in infection, which occurs due to a double exposure situation. This is in parallel with the uncertain mechanisms of infection for a covid-19 pandemic: it has been reported that, until this research work is conducted, there is no evidence whether people have heard or long-term immunity against this covid-19 attack [16], which makes an additional exposure to this novel coronavirus could result in the possibility of infection.

Optimal control is one of the most powerful and effective mathematical tools through which one can formulate the control program for controlling various contagious diseases. We apply the optimal control methods to formulate control mechanisms for the spread of the covid-19 SEIQR model. This technique has been applied over the last few years to prevent various contagious diseases from the community. There are different methods for determining the optimal control for a given mathematical model. For instance, Pontryagin’s maximum principle [15] allows the optimal control calculation for an ordinary equation model system with a given constraint. The aim of this work is to set up an optimal control problem relative to the covid-19 SEIQR model.

To the best of the author’s knowledge, very limited attention has been given to study the convex incidence rate of an epidemic model with control strategies and sensitivity analysis. Therefore, motivated by the above studies, we reformulated the covid-19 SEIQR model by incorporating the convex incidence rate and develop the control strategy for the proposed model [17]. In Section 2, we present the formulation of the model and its basic reproduction number. The local asymptotic stability (LAS) and the global asymptotic stability (GAS) of disease-free equilibrium (DFE) and endemic equilibrium (EE) are discussed in Section 3. In Section 4, we discuss optimal control theory using two control variables. Numerical simulation evaluated in Section 5. The bifurcation analysis sensitivity analysis is performed in Sections 6, 7 respectively while concluding remarks are discussed in Section 8.

Model formulation

Due to a large number of reported (covid-19) cases in the KPK province of Pakistan. The authors are interested to formulate a mathematical model for the covid-19 transmission mechanism. The authors divided the total population into five subclasses namely; susceptible individuals which denotes individuals vulnerable to the infection; exposed class which denotes individuals which take the covid-19 infection but not yet infected i.e., latent class; infectious which represent the individuals that are fully infected and ready for spreading the infection to the others; quarantined-treated and refers to recovered as well as those individuals whose immune system is strong enough to be safe from covid-19 infection at time t.

Further, the incidence rate plays a major role in epidemic models. According to reports of WHO the covid-19 spread very past in a short time worldwide, that’s why the authors considered that the convex incidence rate is suitable for our model. Because for a short time period it corresponds to an increased rate of infection in which cases due to the two exposures. The contact rate results from a single contact, while the term correspond to a new infective which comes into being due to the double exposure.subject to the conditions

The parameter b is the constant birth rate, and are positive constants and are the disease transmission coefficients, is the natural death rate, is the death rate for the disease of infectious individuals and are the state transition rates. represents those individuals whose immune system is strong enough to prevent them from covid-19 infection (see Fig. 1 ).

Fig. 1

The flowchart represents the covid-19 epidemic model.

Letting

Then

Now as , we obtain . We study in the following closed set

Equilibrium of the model

The DFE of the model is denoted by , where

This equilibrium is used to compute the threshold number , which is actually the basic reproduction number. This number is then conveniently used to study the endemic steady state as well. To calculate the number for our model, it is preferred to use van den Driessche and Watmough method [18].

Hence

Endemic equilibrium point

System (1) is rearranged to get and in terms of . Thus

Using the above values in the following equation,gives uswhere

By the Descarte’s rule of sign if then (4) has one positive one negative and one zero roots, and implies that . Hence, for , a unique positive equilibrium exists for the model.

Local stability

We establish the local stability of the system (1) in this section at disease free point as well as at endemic equilibrium point .

The disease free equilibrium (DFE) point of the system (1) is locally asymptotically stable if .

Jacobian matrix of the system (1) at isThe characteristic equation of takes the following form:whereIt is clear from the characteristic equation that the first three eigenvalues are negative whereas for the remaining factors we need to find condition under which the real parts of the two eigenvalues are negative.

We see that , and the value of C can either be positive or negative. Therefore, considering the value of C we have the following cases.

When , then the nature of the roots is real with one root positive and the other is negative. Thus DFE is unstable.

When , then the root are real, negative and zero, DFE is stable (not asymptotically), some trajectories will not be heading to the DFE for .

When then the roots are real, negative and repeated, CFE asymptotically stable.

When we obtain two real and negative roots, hence the DFE is asymptotically stable.

When then the roots are complex conjugate with negative real part, hence DFE is asymptotically stable.

Again from the five cases we can decide that when then DFE is stable, when , the DFE is stable (not asymptotically), and if (cases ) then the DFE is asymptotically stable. All trajectories are moving towards the DFE. So the stability of the DFE point depends only on the value of C.

Furthermore, the value of if This shows that for , the disease free equilibrium point is asymptotically stable. □

At endemic equilibrium point

Let M be a real matrix. If and are all negative, then all eigenvalues of M have negative real part.

If , then the endemic equilibrium of the system (1) is locally asymptotically stable.

The Jacobian matrix of the system (1) at isHere and for the remaining three eigenvalues we have the following reduced matrix:From the Jacobian matrix we havealso Further the second additive compound matrix iswhereHenceTherefore, by Lemma 1, the endemic equilibrium point is locally asymptotically stable. □

Global asymptotic stability

At disease free equilibrium point

The Castillo Chavez approach [19] is summarized as, let the proposed model (1) is reduced into the following two subsystems given bywhere and are the number of uninfected and infected individuals respectively. Thus and . Let us denote the disease-free equilibrium by and define as . For disease-free equilibrium the existence of global stability depends onwhere is an M-matrix having the positive off-diagonal entries and represents the feasible region. Thus the following statement holds.

If is globally asymptotically stable.

, where for ,

The equilibrium point of the system (1) is globally asymptotically stable, if the above conditions are satisfied and .

Now for proving the global stability of (1) at disease-free equilibrium, we apply the above technique.

At covid-19 free equilibrium , the model (1) is globally asymptotically stable if and unstable otherwise.

Consider , and define , where

From the model (1), we have

For , and , we get

From Eq. (13) as , . Thus is globally asymptotically stable.

Now

As , hence . Clearly, B is M-matrix and hence both the conditions are proved, so by Lemma 2, the disease free equilibrium point is stable globally asymptotically.

Endemic equilibrium (global stability)

For the global stability of (1) at endemic equilibrium , we use the geometrical approach [20]. To investigate the sufficient condition through which the is globally asymptotically stable. Therefore, consider the differential equationwhere the open set is simply connected and is a function such that . Assuming that is the solution of Eq. (15) and for , the following are true.

3. There exist a compact absorbing set .

4. System (15) has a unique equilibrium.

The solution is said to be globally asymptotically stable in U, if it is locally asymptotically stable and all trajectories in U converges to the equilibrium . For , a condition satisfied for f, which precludes the existence of non-constant periodic solution of Eq. (15) known is Bendixson criteria. The classical Bendixson criteria for is robust under [20]. Furthermore a point is wandering for Eq. (15), if there exist a neighborhood N of and , such that is empty for all . Thus the following global stability principle established for autonomous system in any finite dimension.

If the conditions and Bendixson criterion are satisfied for Eq. (15) (i.e. robust under local perturbation of f at all non equilibrium, non wandering point for Eq. (15) ), then is globally asymptotically stable (GAS) in U provided it is stable.

Define a matrix valued function P on U by

Eq. (16) is a matrix valued function on U. Further assume that exist and is continuous for . Now define a quantity define, such thatwhere is the second additive compound matrix of J i.e. and . Let be the Lozinski measure of the matrix B with respect to the norm in [20] defined by

Hence if , which shows that the presence of any orbit that give rise to a simple closed rectifiable curve, such as periodic orbits and heterocyclic cycles.

Let U is simply connected and the condition are satisfied, then the unique equilibrium of Eq. (15) is globally asymptotically stable (GAS) in U, if .

If , then the model (1) is globally asymptotically stable (GAS) at endemic equilibrium and unstable otherwise.

Now examine the sub-system of (1), that is

The Jacobian matrix becomes

The compound second additive matrix iswhere

The function , then

So, it follows

The matrix , can be expressed in the formwhere

Suppose be a vector in , with the norm defined by

Let be the Lozinski measure with respect to this norm, we choosewhere

are matrix norms with respect to the vector norm, and refers the Lozinski measure with respect to this norm. Then and

Hence, we have

From , we get

Thus, we have

Moreover, we getthenwhich implies , based on Theorem of [21]. We prove that positive equilibrium is GAS.

Examine the sub system of system the limit system is given by

The solution of last equation yields,implies thatwe get that is GAS.

Optimal control strategies

To develop the control strategy, we use optimal control theory [28], [29]. To assess the conditions necessary to indicate the possible optimal control in the proposed epidemic model, we use Pontryagin’s maximum principal [15].

It should be noticed that the social distancing of the vulnerable population and treatment of the infected people makes sense. To incorporate control strategy into the study, we will update the system by including the controls in the form of treatment and, social distancing. The first control function refers to show the fraction of the population practicing social distancing. The second control function represents treatment of infected individual. For minimization of the objective functional, we havesubject to the state system,with the initial condition for the above system asIn the objective functional (24), and and are constants and nonnegative which represent the relative weights of the infected individuals. Similarly, in the objective functional (24), and are positive constants and it measures the associated cost on social distancing and treatment, respectively. It is to be noted that (24) translates a clear objective: to decrease the number of susceptible, exposed and infected persons while keeping in view the costs on controls. Our goal is to find the control function, such that,subject to the system (25), where the control set is defined as,

The conditions that an optimal solution must satisfy is obtained by using Pontryagin’s maximum principal [15]. This principal translates Equations (24), (25) into a problem characterizes with minimizing the following Hamiltonian H, with regard to control variableswhere and made up of the adjoint variables. The system solution is determined by taking the partial derivatives of the Hamiltonian (28) with respect to the associated state variable.

In view of the optimal controls be the solution of the above control system (24), (25) , then we can find the adjoint variables , for , satisfying where and with transversality conditions

Further more, the optimal controls variables , are defined by

where, .

We use the condition (29), and differentiated the function H with respect to the control variables, to get the corresponding adjoint system. We therefore have the following adjoint system

Moreover, we apply the condition to get the optimal controls given in (31), for .

Numerical simulation

This section is devoted to the numerical solutions of the optimal system and the associated results of different optimal controls scenarios and . The numerical solutions of the optimal system are determined with the help of MATLAB. The adjoint system, was determined to achieve the optimal control solution of the proposed system. A fourth-order Runge–Kutta iterative scheme is used to achieve the solution of the optimal system. To be more specific, we will code state system (25) using forward RK4 method and backward scheme for the adjoint Eq. (32). The scheme will also include the initial data related to state variables and the transversality conditions (30) for the adjoint system. The desired time interval for simulation is taken . We set some parameter values for the system as mentioned in Table 1 . Moreover, the weight constant are assumed to be and . Consequently, we get the following findings from Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6 . Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6 represents the graphs of susceptible, exposed, infected, quarantine-treated and recovered-treated populations. We consider the most highly sensitive parameter i.e and change its value from to . The red color depicts “without control” graphs while the black color depicts “with control” graphs. The dotted graphs are the results when , while the line graphs are the results when .

Table 1

Parameters value.

Notation	value	Reference
b	0.04day^-1	[36]
δ	0.2day^-1	assumed
μ	0.1day^-1	[36]
d_1	0.2day^-1	[36]
d_2	0.3day^-1	[36]
d_3	0.4day^-1	[36]
λ	0.3day^-1	[36]
β	0.5day^-1	assumed
γ	0.3day^-1	[36]
η	0.6day^-1	[36]
τ	0.1day^-1	assumed
∊	0.3day^-1	[36]

Fig. 2

The plot demonstrates the dynamic of susceptible population.

Fig. 3

The plot demonstrates the dynamic of exposed population.

Fig. 4

The plot demonstrates the dynamic of infected population.

Fig. 5

The plot demonstrates the dynamic of quarantine population.

Fig. 6

The plot demonstrates the dynamic of recovered population.

The simulations clearly illustrate the effect as well as the desired goal of the controls, i.e., to reduce the number of exposed, infected, individuals, and to increase . There is obviously a significance difference between the controlled and without control cases. This positive influence suggests that over the intervention period, the control technique is useful in controlling the disease.

Case Study: Khyber Pakhtunkhawa (Pakistan)

As other provinces of Pakistan, the Khyber Pakhtunkhawa province is also effected by covid-19 virus. We choose the Khyber Pakhtunkhawa province as a case study to fit our model to the real data of covid-19 cases from this province from Feb to May 2020. We use Matlab minimization technique and consider the following initial value in which and are estimated while the remaining values are taken from [35]

We set some parameter values as mentioned in Table 2 . Fig. 7 depicts the total cases of covid-19 in the Khyber Pakhtunkhawa province from Feb to May 2020, which is a period of two and a half months. In Fig. 8 we fit real data with the infected class of our covid-19 model which clearly shows the appropriateness of behavior of the infected class. Fig. 9 shows the long time behavior of the covid-19 cases vs time (months). We can see that the data fits accurately to the model curve, and further, one can observe that the number of cases grows exponentially with time. This case could be alarming that the incidence may increases further in the coming months if the government does not apply the proper optimal control strategies.

Table 2

Parameters value.

Notation	Value	References
b	0.028	[35]
α	0.2	Estimated
μ	0.011	[35]
d_1	0.2	Estimated
d_2	0.08	Estimated
d_3	0.06	[35]
λ	0.3	Estimated
β	0.5	Estimated
γ	0.3	[35]
η	0.6	Estimated
τ	0.1	[35]
∊	0.3	[35]

Fig. 7

The incidence data of covid-19 from Khyber Pakhtunkhwa, Pakistan.

Fig. 8

Covid-19 model comparison with real data of Khyber Pakhtunkhawa, Pakistan.

Fig. 9

The incidence data of covid-19 from Khyber Pakhtunkhwa, Pakistan.

Bifurcation analysis

Bifurcation analysis plays a pivotal role in controlling and decimating epidemic diseases. In modeling epidemic diseases, to eradicate the disease, a classical condition is imposed on basic reproduction number i.e. must be less than unity [30], [31]. But on the other hand, if a bifurcation takes place, the existence of an endemic equilibria for would imply that the disease will persist over time. This goes to show the importance of bifurcation analysis in epidemiology and preventive medicines.To this end, we follow [32], [33] to study the bifurcation phenomenon in our model.

Let , where and assume a general system of ODE’s with a parameter

Without loss of generality, we assume that is an equilibrium for system (1).

If represents the linearized matrix of system (1) around the equilibrium point with evaluated at 0. Zero is a simple eigenvalue and all other eigenvalues having negative real parts.

Matrix A has a non-negative right eigenvector w and a left eigenvector v corresponding to zero eigenvalues.

Let denotes the component of g, and

Then, the local dynamics of system (1) around are determined by the coefficients a and b.

. When with , is locally asymptotically stable and there exists a positive unstable equilibrium, when is unstable and there exist a negative and locally asymptotically equilibrium;

. When with is unstable; when is locally asymptotically stable and there exist unstable positive equilibrium;

. When with is unstable and there exist a locally asymptotically stable negative equilibrium, while in case of is stable and there exist a positive unstable equilibrium;

. When varies from negative to positive, changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable.

The above Theorem prescribes the role of the coefficients a and b of the normal form representing the system dynamics on the central manifold, in deciding the direction of the transcritical bifurcation occurring at .

More precisely, if and , then the bifurcation is forward; if and then the bifurcation is backward.

Taking into account the above theorem we show that our proposed model (1) exhibits bifurcation by using central manifold theory [34]. We put to find as

We get the Jacobian at disease free equilibrium as

The corresponding eigenvalues of are and .

Since one of the eigenvalue is zero, hence it confirms the presence of bifurcation. Following [34], we find the left eigenvectors of as

Also, we find the right eigenvectors as

Setting , and , in (1), we get

Taking into account the left and right eigenvectors we calculate the coefficients a and b asimplies

Alsoimplies

As and , hence it confirms that our proposed model exhibits forward bifurcation (1). The forward bifurcation is highlighted in Fig. 10 .

Fig. 10

“The solid line (–) denotes stability; the dashed line (- -) denotes instability”.

Sensitivity analysis

Determining the parameters which are helpful in decreasing the spread of infectious disease is carried out by sensitivity analysis. Forward sensitivity analysis is considered as vital component of disease modeling although its computation become tedious for complex biological model. Sensitivity analysis of have received much attention from the ecologist and epidemiologist.

The basic reproduction number is the normalized forward sensitivity index that depends differentiability on a parameter is defined as

Three methods are normally used to calculate the sensitivity indices, (i) by direct differentiation, (ii) by a Latin hypercube sampling method (iii) by linearizing system (1) and then solving the obtain set of linear algebraic equations. We will apply the direct differentiation method as it gives analytical expressions for the indices. The indices not only shows us the influence of various aspects associated with the spreading of infectious disease but also gives us important information regarding the comparative change between and different parameter. Consequently, it helps in developing the control strategies.

Table 3 shows that the parameters , and have a positive influence on the reproduction number , which describe that the growth or decay of these parameters say by 10 percent will increase or decrease the reproduction number by 10 percent, 1.0 percent and 9.0 percent, respectively. But on the other hand, the index for parameters and illustrates that increasing their values by 10 percent will decrease the values of reproduction number by 1.4 percent, 7.1 percent, 8.5 percent, 2.9 percent, 8.1 percent, 0.1 percent and 0.5 percent, respectively.

Table 3

Sensitivity indices of the reproduction number against mentioned parameters.

Parameters	S.Index	Value	Parameters	S.Index	Value
β	S_β	1.00000000	d_2	S_d_2	−0.8593749998
b	S_b	1.00000000	μ	S_μ	−0.2965237119
γ	S_γ	−0.1428571429	d_3	S_d_3	−0.8196721309
d_1	S_d_1	−0.7142857142	η	S_η	−0.01639344262
λ	S_λ	0.9062499999	∊	S_∊	−0.05714285714

Fig. 11 shows sensitivity of against and . From Table 3, it is clear that both have negative impact on , which is depicted in Fig. 11 as it varies very slowly from to its maximum value . It shows that the death rate of infectious individuals and those individuals whose immune system is strong enough to prevent them from covid-19 infection has no impact in variation of basic reproduction number.

Fig. 11

The plot demonstrates the variation of against and .

Fig. 12 shows sensitivity of against and . From Table 3, it is clear that has negative while has positive impact on , which is depicted in Fig. 12 as no variation is seen when is increasing while very fast increasing in is felt when is increasing. It shows that the death rate of infectious individuals have no impact on while the transmission rate of covid-19 infection from susceptible people to infected have a very high impact on basic reproduction number.

Fig. 12

The plot demonstrates the variation of against and .

Fig. 13 shows sensitivity of against and . From Table 3, it is clear that and both have positive impact on , which is depicted in Fig. 13 as very big variation of is seen when and is increasing. It shows that the transmission rate of covid-19 infection from susceptible and exposed people to infected people have a very high impact on basic reproduction number.

Fig. 13

The plot demonstrates the variation of against and .

Fig. 14 shows sensitivity of against and . From Table 3, it is clear that has negative while has positive impact on , which is depicted in Fig. 14 as no variation is seen when is increasing while very fast increasing in is felt when is increasing. It shows that the transmission rate of covid-19 infection from infected people to recovered people have no impact on while the transmission rate of covid-19 infection from exposed people to infected have a very high impact on basic reproduction number.

Fig. 14

The plot demonstrates the variation of against and .

Conclusion

In this article, we performed the stability and bifurcation analysis of covid-19 SEIQR epidemic model with convex incidence rate. First, we investigated the basic reproductive number for the system (1) which completely characterized the stability of the DFE and EE. We obtained sufficient conditions for LAS and GAS. The approach we used for the global stability of DFE was the Lyapunov function theory, whereas, the global stability of EE was carried out using the geometrical approach. Furthermore, we have developed an optimal control strategy with two control variables to minimize the infection in the host population. The controls under consideration are: social distancing and treatment. The sensitivity analysis complemented by simulations are also performed to assess how changes in the parameters affect the dynamical behavior of system (1). It is concluded that with the implementation of such controls strategy for the long run, it is possible to eradicate epidemic diseases from community.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.