Sensitivity analysis and optimal control of COVID-19 dynamics based on SEIQR model.

It is of great curiosity to observe the effects of prevention methods and the magnitudes of the outbreak including epidemic prediction, at the onset of an epidemic. To deal with COVID-19 Pandemic, an SEIQR model has been designed. Analytical study of the model consists of the calculation of the basic reproduction number and the constant level of disease absent and disease present equilibrium. The model also explores number of cases and the predicted outcomes are in line with the cases registered. By parameters calibration, new cases in Pakistan are also predicted. The number of patients at the current level and the permanent level of COVID-19 cases are also calculated analytically and through simulations. The future situation has also been discussed, which could happen if precautionary restrictions are adopted.

Introduction

The World Health Organization (WHO) was informed by the Wuhan Municipal Health Commission Hubei Province of China about a class of twenty-seven cases of unknown aetiological pneumonia, that were generally reported to the fish market in Wuhan city in December 2019. Also, it was observed that seven out of twenty-seven patients had critical condition. Symptoms in first case appeared on December 2019. Chinese authorities, On 7th January, 2020, recognized a new one strain of family virus which caused the outbreak that belonged to Corona = virus family. It was called SARS-CoV-2. The medical observation related to the virus had given name COVID-19. By 11th of March, the World Health Organization (WHO) announced the pandemic globally [1]. Currently, the globally declared COVID-19 cases are million with almost million mortalities.

Corona-virus is the viruses family which causes infection in some animals and humans. Corona-viruses are zoonosis/vector-born diseases, that is, from animals to humans they will pass on [2]. Human affected Corona-virus(HCoV) may produce clinical symptoms such as those caused by the Middle East Respiratory Syndrome (MERS-CoV) and extreme Severe Acute Respiratory Syndrome (SARS) from the common cold to serious illness [3]. SARS-COV-2 transfers from animal to human and human to human. The first one is still not yet known, but some researchers say that it may transfer through digestive system material and respiratory secretions. The second is close to infection with the other Corona-viruses with the excretions of infected persons, primarily by direct exposure with respiratory droplets and hands or object contaminated with these excretions, accompanied by exposure with the mucus of the eyes, nose and mouth [4], [6].

Officially, the novel Corona-virus (COVID-19) pandemic was confirmed to have reached in Pakistan on 26th of February, 2020. The very first patient was reported in Province Sindh and second was reported in the capital territory Islamabad [5]. In less than a week after the appearance of the first two cases, this virus started to grow in different regions of the country. On April, 29, 2020, confirmed number of cases in Pakistan were 15,827, with 4052 ( of the cumulative cases) recovery and 361 ( of the cumulative cases) deaths and Punjab was the province having the higher number of increasing cases more than 6000 [7]. This pandemic continuously increased till the start of July. July 2, 2020 was the peak day of active cases. There were 108,466 COVID-19 patients in Pakistan on July 2, 2020. After that day decline started and the infected number of patients from COVID-19 decelerate to 5525 till September 15, 2020. After the mid of September, the second wave started in Pakistan and the number of active cases again accelerate to 55,354 till December 7, 2020 [8].

Modeling is a science of artistic abilities linked to comprehensive learning in the form of mathematical models in a range of techniques to explain physical phenomena [9], [10], [11], [12]. In the current situation, the disease control authorities maintain and report all COVID-19 data on a regular basis. This data set includes the number of individuals who have positive Corona-virus test, total number of deaths, total number of recovered individuals, and total active number of cases, and whole combined data from every part of the world. Hence, at this point an acceptable model with greater accuracy is needed. Lower dimensional mathematical models, with lower numbers of parameters and population compartments, that could be better reliably determined by the actual data, are preferable for studying and forecasting the ongoing pandemic [13]. Any model with higher number of dimensions need a large number of parameters to be expressed it, but this large number of parameters are difficult to find with adequate accuracy [14], [15]. Compartmental epidemic mathematical models representing the overall behavior and fluctuation of the system may be an initial point because of a lack of data. Also the simplest models have several variables, that are difficult to quantify from the current data [16], [17]. The least SIR model illustrates the compartment of the susceptible , the infected , and the removed (recovered or deceased) populations [18], [19].

Yousaf et al. [20] performed statistical analysis for the month of May, for Pakistan, by using vector Auto Regressive Integrated Moving Average Model . According to their model prediction, the confirmed number of cases may increase up to 33,079 (Fig. 1 a). Chaudhry et al. [21] used the Social Sciences statistical package (SPSS) version 23 (IBM Corp., Armonk, NY), and forecast was given by utilizing a simple moving average method in time series modeler. In Pakistan, they predicted that active cases of COVID-19 may go over 35,000 at the end of May (Fig. 1b). Muhammad Aslam [22] used Kalman filter with ARIMA for COVID-19 pandemic in Pakistan. He gave prediction of five days from May 1, 2020 to May 5, 2020 (Fig. 1c). Dil et al. [23] predicted COVID-19 number of patients in the Eastern Mediterranean Region. They estimated that by the June 2020, cases will shoot up to 510,000 (Fig. 1d). Vector Auto-regressive model has been used by Khan et al. [24] to predict COVID-19 cases in Pakistan. Their 10 days forecasting (June 29, 2020–July 8, 2020), is shown in Fig. 1e. Ahmad and Asad [25] used very efficient technique, artificial neural network, to predict 7 days (July 11, 2020–July 17, 2020) active cases of COVID-19 in Pakistan (Fig. 1f). Mathematical analysis and prediction of COVID-19 cases has been discussed in Pakistan by Ali et al. [26] and Ozair et al. [27]. Their prediction show the decline of COVID-19 in December, 2020. Many other efficient studies can be seen in the literature for the better understanding of COVID-19 pandemic and to design effective control strategies. For example Boccaletti et al. [28] discussed three scientific communities that may play vital role in controlling this pandemic. Atangana [29] proposed a mathematical model based on fractal-fractional differential and integral operators to discuss the indication of spread profile through statistical figures of Italy. Khan and Atangana [30] discussed the dynamics of COVID-19 appeared in China by construction the mathematical model having Atangana-Baleanu fractional derivative. Although face mask is the powerful tool to protect against COVID-19. However, hypercapnia has been recognized which is related to symptoms of fatigue, muscular weakness, discomfort, headaches as well as drowsiness. Ernestine Atangana and Abdon Atangana [31] provided the literature review in order to determine which type of face masks should be used to avoid re-inhaling rejected . Furthermore, mathematical models representing the transport of COVID-19 spread through wind has also been designed. The comprehensive statistical analysis of COVID-19 cases in Turkey and South Africa has been given by Abdon Atangana and Seda İğret Araz [32]. Recently, Atangana and Araz [33] utilized stochastic and deterministic approaches for forecasting of COVID-19 in Africa and Europe. The interested readers may also go through [34], [35], [36], [37]. To evaluate the dynamics of second wave, we modeled the transmission through a deterministic model SEIQR for COVID-19 infection in Pakistan.

Fig. 1

Prediction of COVID-19 patients in Pakistan using various models and methods.

The Model structure

The whole population is distributed into five groups in the SEIQR style model, the susceptible , exposed , the infectious , quarantine and the recovered . We suppose that susceptible humans are recruited at the rate and leave the compartment through natural mortality rate as well as after getting infection. Humans belonging to all classes have natural mortality rate . The class having COVID-19 infectious patients (I) has additional mortality rate . Corona Virus is transmitted in susceptible humans, at the rate , as a result of exposure of infectious humans with susceptible ones. Exposed individuals switch with the rate to the infectious group. Transition of infectious individuals occurs at the rate . They move to either quarantine or recovered class. The portion of infectious individuals enter in the quarantine class (at the rate ) and the remaining portion , after recovery, join the recovered class. The detailed description of the parameters and their dimension is given in Table 1 . The assumption that people are homogeneously mixed leads to the following mathematical form of the model:with initial conditions and .

Table 1

Definition, dimension of the model and parameter values (1).

Parameter	Description	Dimension
β	The effectual contact rate among susceptible and infectious individuals	day −1
μ	Death rate of individuals (Natural)	day −1
α_1	The rate of exposed individuals becoming infectious	day −1
γ	The recovery rate of infected ones	day −1
1-δ	The portion of infected individuals
	that are quarantined or recovers (at the rate γ)	dimensionless
δ	The resting part of infected individuals who are quarantined	dimensionless
Λ	The recruitment rate of susceptible individuals	day −1
ϕ	Quarantine rate of infections	day −1
η	Disease related death rate	day −1

Analysis of the model

Positively invariant set

We can easily verify that all classes of individuals () are non-negative for all with non-negative initial conditions. Furthermore, the closed setis positively invariant.

Equilibrium points and basic reproduction number

The model has two types of equilibrium points, first one is disease free (DF), , and second is endemic equilibrium (EE), . We can write very easily our , by simply putting the disease classes and derivatives equal to zero in (1), it would have the form . The complete dynamics of the ailment is contingent on a threshold parameter, known as basic reproduction number, which is defined as the mean number of the secondary infections or secondary infections that are produced by an infected person during his total infectious period. By using this parameter, one can know that whether the disease persists in the community or not. With the help of next generation method [38], we derive the expression for basic reproduction number as

For , let = be endemic equilibrium, where,

Estimation of parameters

We estimate the values of parameters from available COVID-19 data and literature published. In given Fig. (2 ), only active cases with recovered ones and deaths have been plotted from of Sep., 2020 to December 15, 2020. Pakistan currently has a total population of nearly 220 millions [39] and life expectancy of 67 years [40]. Different studies show that the incubation time period for COVID-19 ranges from 2 to 14 days and almost people show disease symptoms within days of infection [41], [42], [43]. We estimate the remaining parameters by fitting the model with active number of COVID-19 patients. The second wave of Corona-virus appears in Pakistan on September 15, 2020. From that day, the COVID-19 patients have been continuously boosting up till December 6, 2020. It means that the peak period of second wave of the pandemic was December. Now the number of active cases are decreasing day by day. So we divide the analysis in two phases. In the first phase, we estimate the number of cases by fitting the model with the increasing trend and in second phase, we calibrate our model results by the deceleration of number of active cases in Pakistan.

Fig. 2

Number of COVID-19 patients, Recoveries and Deaths in Pakistan for second wave.

Phase-I

Here we study the increasing pattern of COVID-19 patients which starts from September 15, 2020 lasts up to December 7, 2020. Fig. (3) shows the number of actual cases that happened during the said period. We simulate the model and estimate the patients by finding suitable values of the parameters. The comparison between actual and fitted data is expressed in Fig. (4) . The estimated parameter values are shown in Table 2 .

Fig. 3

Actual COVID-19 patients in Pakistan till December 7, 2020.

Fig. 4

Comparison of actual and estimated cases of COVID-19 till December 7, 2020.

Table 2

Definition, dimension of the model and parameter values (1).

Parameter	Description	Dimension	Value
β	The effectual contact rate among susceptible and infectious individuals	Day −1	2.43269 ×10^-9
μ	Death rate of individuals (Natural)	day −1	0.0000408914
α_1	The rate of exposed individuals becoming infectious	day −1	0.0869565
γ	The recovery rate of infected ones	day −1	0.082
1-δ	The portion of infected individuals that are quarantined or recovers (at the rate γ)	Dimensionless	0.131
δ	The resting part of infected individuals who are quarantined	day −1	0.896
Λ	The recruitment rate of susceptible individuals	day −1	8996.11
ϕ	Quarantine rate of infections	day −1	0.906
η	Disease related death rate	day −1	0.262

Phase-II

The decline in number of COVID-19 patients in Pakistan occurs after December 7, 2020. Small variation appears between December 8, 2020 and December 15, 2020. After this date, the pandemic is continuously disappearing. We again simulate the model by consider active cases of December 16 as the initial value and estimate the parameter values. Fig. (5) show the actual deceleration in number of patients of COVID-19 and Fig. (6) shows the fitting data with the actual cases. The values of estimated parameters is shown in Table 3 .

Fig. 5

Actual Number of Patients from December 15 to 23.

Fig. 6

Comparison of Estimated and Real Cases from December 15 to 23.

Table 3

Definition, dimension of the mode and parameter values (1).

Parameter	Description	Dimension	Value
β	The effectual contact rate among susceptible and infectious individuals	day −1	00.000000000709267
μ	Death rate of individuals (Natural)	day −1	0.0000408914
α_1	The rate of exposed individuals becoming infectious	day −1	00.0869565
γ	The recovery rate of infected individuals	day −1	00.0916795
1-δ	The portion of infected individuals are quarantine or recovers (at the rate γ)	Dimensionless	00.12
δ	The resting part of infected individuals who are quarantined	day −1	00.88
Λ	The recruitment rate of susceptible individuals	day −1	8996.11
ϕ	Quarantine rate of infections	day −1	0.172
η	Disease related death rate	day −1	00.0000000000139611

Sensitivity analysis

To contemplate capacity of pandemic to get in a populace is our primary concern. Sensitivity analysis is carried out to analyse those factors due to which this malady disseminates and persists in the community. We target that parameters due to which larger deviation is produced in the value of basic reproduction number.

Sensitivity indices of

The corresponding variance in the state variable, by the variation of a parameter, can be calculated through the sensitivity indices. Calculation of these indices has been made through the definition given in [44]. The sensitivity index is defined, in the form of partial derivatives, as given below:

By definition given above, the analytical formulas for the sensitivity indices of can be calculated. We can see from Table 4 that the parameter which shows most sensitivity is the adequate contact rate among susceptible and infectious individual followed by the recovery rate of infectious individuals. These indices indicate that if we able to enhance the recovery rate and reduce the adequate contact rate by almost 10%, the value of basic reproduction number may reduce up to 10%.

Table 4

Sensitivity indices of corresponding to all parameters

Parameter	Sensitivity index	Parameter	Sensitivity index
Λ	1	μ	1.00092
α_1	0.00047003	β	1
γ	-0.999554	η	1.522136984456319×10^-10

Reduction in basic reproduction number for different parameters

By using the parameter values given in Table (3), we plotted the reproduction number as a function of and and and . Fig. (7a) shows the effect of effectual contact rate and incubation period of exposed patients on the reproduction number, we can see that these two parameters may decrease to less than unity if the value of remain up to and the incubation period remains between 2 and 14 days. The combined effect of transmission and recovery rate on is shown in Figs. (7b). As one parameter () is directly proportional and is inversely related to so we can see that the basic reproduction number value may be reduced below unity even if the transmission rate approaches to and the infectious period increases to 5 days. Fig. (7c) shows the effect of transmission rate and . We can see that these two parameters have considerable impact in reducing or increasing the value of basic reproduction number .

Fig. 7

Contour Plots of for different parameters.

Optimal control

Sensitivity analysis discussed in the previous section guides us to design effective control strategy in order to overcome this pandemic. As we applied the sensitivity analysis on the estimated parameters that have been obtained from the comparison of actual data and the outcomes of model. So, we can see almost clear picture of the future behaviour of this pandemic. By observing the most influential parameters, we reshape the model (1) to survey the influence of control measures regarding the future scenario. After inclusion of two controls, the model (1) gets the following form:

The control represents the precautionary measures such as social distancing, abundant use of face mask etc. to reduce the transmission of infection into healthy people. Control represents the quick action to isolate the people having COVID-19 virus. To discuss the ideal attempts to restraint the disease, we formulate the objective functional J. Its aim is to curtail the infectious individuals as well as the minimization of cost of applied controls and . One haswhere represent the positive weights. Our main objective is to curtail the number of infectious whilst reducing the cost of above mentioned controls with the help of described objective functional. Controls and can be found in such a way thatin which, U, which is called control set, and is given by

By the use of Pontryagin’s Maximum Principle [45], we can have the essential conditions and solution of (2).

Existence and analysis for optimal control

The existence and analysis of an optimal control could be verified with the help of powerful classical procedure. According to [46], it is necessary to check that the following hypotheses hold:

Controls set and corresponding described variables is nonempty.

U must fulfills the two properties the convexity and closeness.

Right hand side of the described system is bounded by a linear function in the control state.

Objective functional J has convex integrand on U and bounded below by , where and .

The result in Lukes ([47], Th 9.2.1, p 182) is utilized for existence of solutions of system (2). In this way, we verify the above hypothesis. is satisfied due to boundedness of coefficients and solutions, boundedness shows that U fulfils the second hypothesis. As the system (2) is bilinear in and due to the boundedness of solutions, right hand side of (2) fulfils the third hypothesis . Since the control functions and are convex and number of infectious hosts and vectors are also bounded at any time t by the positivity of solutions, so we can choose constants all are positive and such that the inequalityis satisfied. Thus the last condition is also verified. So, we have the theorem given below:

For objective functional where subject to Eqs. (2) with initial conditions, there exists an optimal control so that .

If we can have the Lagrangian and Hamiltonian of the system (2), then one can get the optimal solution. Lagrangian, for (2), is given by

One has to build up the least possible value given above. To get objective, we have to define the Hamiltonian H in this way:

Let us consider and , so we have:

The optimality system

Maximum Principle explored by Pontryagin as reported in [48] is used for finding essential optimal control conditions. This is calculated as follows:

For optimized solution of optimal control problem (2), a vector function , which is nontrivial, prevails that fulfils the ensuing conditions. The stated relation is

Now, mandatory conditions are incorporated for H.

For and solutions of the system (2) , variables occurs having, with transversality conditions . Additionally, are expressed as

By making use of H, we can get the adjoint equations and the transversality conditions. Set and differentiate the Hamiltonian H with respect to variables . Hence we havewith transversality conditions . Optimality conditions and the property set U of controls give the result

Outcome of the optimal controls on the transmission of COVID-19 can also be observed numerically. Horizontal axes represent time in days and vertical axes indicates individuals in the specific compartment. we start the control from December 16, 2020. The actual number of COVID-19 patients was . We take this number as an initial value and simulate the modified model for next 30 days. In Fig. 8a , we observe that the number of susceptible humans increases by applying the controls. Fig. 8b shows that the persons having Corona-virus decreases by applying the controls. Fig. 8c also represents the significant reduction in the number of infectious by adopting the control measures. There will be no considerable difference in the number of quarantined individuals. However, Fig. 8d shows that they decrease after 27 days. Fig. 8e shows that recovered people will decrease by applying the controls. It will happen due to significant reduction in the number of exposed and infectious individuals in the coming days. Control profile is shown in Fig. 8f. It shows that, to get the huge reduction in the pandemic, we have to apply both controls with almost efficacy.

Fig. 8

Use of Controls and .

Conclusions

In this study, we have interpreted the deterministic COVID-19 disease model carefully to determine the feasibility of emergence of COVID-19 patients throughout Pakistan. We compared the result of our model with the actual confirmed cases, it has been noted that there is appropriate correspondence between estimated values and actual confirmed numbers. If the current scenario is continuing then, according to our estimation, the active number of COVID-19 cases will decrease gradually. However, the scenario will change by the changing of parameters.

The model has been analyzed qualitatively. Explicit formula for the reproduction number has been calculated. Two types of equilibria, disease free and endemic, has been calculated. The model has been applied on the second wave of COVID-19 in Pakistan. We fitted the model both for increasing and decreasing pattern of COVID-19 data. The parameters have been estimated through the calibration of real data. Active patients in Pakistan are decreasing gradually for past few days. So we focus our analysis on the recent trend.

The parameters having vital impacts on the basic reproduction number have been identified. This sensitivity analysis has been performed by using the numerical values of the parameters obtained through fitting of the model with actual number of active cases showing decreasing pattern. First we calculate the sensitivity indices of basic reproduction number with respect to all the parameters. It has been observed that the most sensitive parameter (ignoring the sensitivity indices of demographic parameters) is the adequate contact rate among susceptible and infectious humans. Secondly, we observed the combined effect of every two parameters on the basic reproduction number.

On the basis of measuring sensitivity, we modified the formulated mathematical model by introducing the control functions. First is the attempt to maintain the social distancing or keeping the practice of wearing face mask and second is the quick action of making the quarantine of infected people. We numerically solve the optimal control problem. It has been calculated that, by the application of control measures, the number of exposed and infectious individuals will be and , respectively. However, these numbers may increase up to and , respectively, till the end of January, 2021, if proper interventions have not been applied. In future, we shall try to prove the global stability of disease free equilibrium as well as endemic equilibrium and, through the parameter calibration, we shall estimate the constant maximum number of COVID-19 patients for both cases that is for disease persistence and disease extinction. Mathematical model can also be modified for example by the addition of hospitalized class.

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.