Modeling the dynamics of COVID-19 pandemic with implementation of intervention strategies.

The ongoing COVID-19 epidemic spread rapidly throughout India, with 34,587,822 confirmed cases and 468,980 deaths as of November 30, 2021. Major behavioral, clinical, and state interventions have implemented to mitigate the outbreak and prevent the persistence of the COVID-19 in human-to-human transmission in India and worldwide. Hence, the mathematical study of the disease transmission becomes essential to illuminate the real nature of the transmission behavior and control of the diseases. We proposed a compartmental model that stratify into nine stages of infection. The incidence data of the SRAS-CoV-2 outbreak in India was analyzed for the best fit to the epidemic curve and we estimated the parameters from the best fitted curve. Based on the estimated model parameters, we performed a short-term prediction of our model. We performed sensitivity analysis with respect to R 0 and obtained that the disease transmission rate has an impact in reducing the spread of diseases. Furthermore, considering the non-pharmaceutical and pharmaceutical intervention policies as control functions, an optimal control problem is implemented to reduce the disease fatality. To mitigate the infected individuals and to minimize the cost of the controls, an objective functional has been formulated and solved with the aid of Pontryagin's maximum principle. This study suggest that the implementation of optimal control strategy at the start of a pandemic tends to decrease the intensity of epidemic peaks, spreading the maximal impact of an epidemic over an extended time period. Our numerical simulations exhibit that the combination of two controls is more effective when compared with the combination of single control as well as no control.

Introduction

The novel coronavirus is bizarre due to several causes, making its transmission unpredictable and hard to control due to its peculiar epidemiological traits. In December 2019, the COVID-19 pandemic has started in Wuhan city of Hubei Province, the Republic of China, caused by the severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) and China became the epicenter, and has now spread throughout the world [1]. As of November 30, 2021, the total number of confirmed cases 26,14,35,768 total number of deaths 52,07,634 [2]. The worldwide pandemic resulting from SARS-CoV-2 was preceded by other two epidemics of human coronavirus, namely Middle East respiratory syndrome coronavirus (MERS-CoV) and severe acute respiratory syndrome coronavirus (SARS-CoV) infections [3], posing the tremendous menace to the economy and global public health after the 2nd World War [4]. Presently one of the major problems is that there is no clear picture or consensus on the future progression of the COVID-19 pandemic. The possibilities of the source of the COVID-19 transmission incorporate (but not limited to) animals, human to human and intermediate animal vectors [3].

Due to non-appearance of any pharmaceutical interventions, government of different countries is incorporating different policy to control this pandemic and the most common one is the implementation of lockdown to maintain the social distancing, and contact tracing [5]. This is an excellent measure to control the spreading of the COVID-19, but from an economic view point, the nationwide lockdown may be the cause of an important financial catastrophe in the near future. Specifically, the lockdown in the 2nd most populated country reduces the contact rate of the diseases, but complete control on disease transmission may not be obtained. Thus, to dynamic the economic situation of a country, a complete nationwide lockdown for an uncertain period is not acceptable at all in any circumstances. The index case for COVID-19 epidemic in India was reported on January 30, 2020, in Trissur district of Kerala when a student returned from Wuhan, China [6]. The Govt. of India implemented one-day ‘Janata Curfew,’ to control the SARS-CoV-2 epidemic and maintain the social distancing on March 22, 2020 [7], and the first lockdown throughout the country for 21 days was implemented on and from March 25, 2020.

Mathematical models of the infectious diseases through a system of differential equations are now ubiquitous. Mechanistic mathematical models based on the system of nonlinear ordinary differential equations (ODEs) may give an important information regarding the transmission dynamics of COVID-19 outbreak and its control as well as mitigation strategies. Kermack and McKendrick [8] constructed a fundamental epidemic SIR model to study the human-to-human disease transmission dynamics of individuals through three mutually exclusive phages of infection, namely susceptible (S), infected (I) and recovered class (R) in the year 1927. More complex mathematical models have been introduced to study the COVID-19 transmission dynamics. Kucharski and colleagues [1] investigated a mathematical model to study the transmission dynamics of COVID-19 outbreak for cases in Wuhan, China (including cases that originated there), using stochastic approach and use data-driven analysis and estimate the data based on likelihood of the outbreak taking place in other geographical locations. Tang and colleagues [9] constructed a compartmental mathematical model to investigate the transmission dynamics of COVID-19 and compute the basic reproduction number , which is very high for COVID-19 outbreak as well as for the infectious diseases. Khajanchi and Colleagues [10] constructed a compartmental model to control and predict the COVID-19 outbreak for the overall India and for the four states of India. The authors also introduce the effect of media to mitigate the spread of COVID-19 outbreak. Sarkar and Khajanchi [11] developed a compartmental mathematical model to investigate the transmission dynamics and forecast the COVID-19 outbreak in seventeen states of India and the overall India. Wu and colleagues [12] investigated a simple SIR model to delineate the transmission dynamics of SARS-CoV-2 virus and also estimate the clinical severity for the coronavirus diseases. Giordano and colleagues [13] developed a new compartmental mathematical model for SARS-CoV-2 epidemics to mitigate the transmission of COVID-19 among the human and suggest to restrictive social distancing including lockdown, closing educational institutions contact tracing, etc.

In an extended version of the classical SEIR (susceptible–exposed–infected–recovered) model to study the intervention strategies for infectious diseases incorporates the fact that asymptomatic and pre-symptomatic infected individuals play a key role in the transmission dynamics of COVID-19 pandemic [14-23]. One oddity is how simply individual can get infected by someone without any clinical symptoms. But there is a distinction between pre-symptomatic and asymptomatic spread. The transmission of virus by persons who do not have symptoms and will never get symptoms from their infection is known as asymptomatic spread of the virus. Thus, the asymptomatic individuals are the carrier of the virus and they could still get others very sick. The transmission of the virus by persons who don’t look or feel sick, but will eventually get symptoms later is known as pre-symptomatic spread of the virus. The pre-symptomatic individuals can also infect other person without knowing it.

The course of an epidemic can be delineated by a series of significant factors, but some of which are very difficult to understand the present SARS-CoV-2 disease. Also, the basic reproduction number is a key identifier for the transmissibility of the infectious diseases, as quantifies how contagious the disease is. The average number of secondary infections by a single infected person in a whole susceptible class is known as , or more precisely the area under the pandemic curve. Again for , the disease due to infection by the virus is expected to stop spreading, but for an infected person can infect on an average single person; that is, the spread of the disease is stable. Hence, the disease due to infection by the virus can spread and become epidemic if becomes larger than unity. can aid in understanding the effectiveness of the diseases, that is, under what condition the disease can spread or stop [24].

Curbing the worldwide spread of COVID-19 needs implementation of multiple population-wide policies; however, how the stringency and timing of such measures will influence ‘flattening the curve’ remains unknown. In order to answer some of the important issues, we have proposed a compartmental mathematical model that forecasts the evolution of outbreaks and aids to assess the effect of various policies to restrain the spread of the infection. We estimate the system parameters, and the data have been taken from the World Health Organization (WHO) [2]. We perform a short-term prediction of the pandemic based on the estimated model parameters, and the results show the increasing trend of the COVID-19 in India.

At the present situation, we maintain the outbreak related to the potentiality to obtain actual data that in turn will capable to use real data for modeling work purposes, which are the perfect skeleton in agreement with upcoming phenomena [25]. Several efforts have been taken into consideration to better understand the kinetics of epidemiological cycle of the coronavirus diseases [13, 26–30]. The rearrangement of the system parameters in a practicable way, through the imposition of limits on the system in getting the optimization of a given function that can be implemented using the theory of optimal control.

Various kinds of mathematical models have been created and established to investigate the interactive dynamics of the SARS-CoV-2 outbreak with their dynamical behaviors [10, 11, 31–39]. The study related to data-driven model for the COVID-19 transmission dynamics with the implementation of distributed time lag has been observed by Liu et al. [40]. Khyar et al. developed a multi-strain SEIR model to study the complicated dynamics of the SARS-CoV-2 epidemic and observed that their model showed global dynamics [41]. The theory of optimal control policies has been used to study the SARS-CoV-2 viral dynamics, and the authors optimal doses are required to curtail the SARS-CoV-2 viruses [42, 43]. Different kinds of qualitative behavior have been observed in various mathematical models [44-47].

The implementation of the control theory in infectious diseases is very familiar: While modeling of the infectious diseases has used that combinations of quarantine, isolation, clinically ill cases and vaccination are necessary to eliminate the infectious diseases, the theory of optimal control suggests us how they should be maintained, by providing actual times for intervention with actual dosages [48]. This optimization procedure has also been used in some advanced research the span of SARS-CoV-2 epidemic. The theory of optimal control of a modified susceptible–exposed–infected–recovered (SEIR) mathematical model has been studied with the aim to understand the impact of two inflexible lockdown policies on the UK [49].

Herein, we are interested in employing the control theory procedure to better understand the ways to maintain the progression of the SARS-CoV-2 outbreak in a case study of India by designing optimal disease intervention strategies. Moreover, we raised some important questions that are not entirely investigated in the existing literature. Our proposed model exhibits the incorporation of optimal control policy, to minimize the clinically ill individuals and isolated individuals where the response applicability of public health policies is perpetuated. The SARS-CoV-2 outbreak has shown us that the public health facility is not only a medical issues but also affects the country as a whole [43]. The theoretical analysis of the transmission dynamics of the SARS-CoV-2 pandemic has also been studied in the present manuscript.

The article is organized in the following way: Sect. 2 briefly describes the assumptions and formulation of the compartmental model. The mathematical analysis is shown in Sect. 2.1. In the next section, we performed rigorous model simulations to validate our theoretical analysis by considering the estimated model parameters for the Republic of India. The system parameters are estimated based on real world example on SARS-CoV-2 for India and performed a short-term prediction of the pandemic based on the estimated parameter values. A discussion in Sect. 6 concludes the manuscript.

The mathematical model

Herein, we proposed a mathematical model that represents the transmission dynamics of SARS-CoV-2 or COVID-19 epidemic. The total population is classified into nine compartments of the COVID-19 disease: S, susceptible (uninfected); E, exposed; A, asymptomatic; P, pre-symptomatic; U, symptomatic infected (undetected); D, symptomatic infected (detected) or severely symptomatic; I, isolation or hospitalization; R, recovered or healed; H, dead or extinct; collectively termed SEAPUDIRH. The SEAPUDIRH dynamical system consists of nine nonlinear ordinary differential equations, delineating the evolution of the individuals in each compartment over time t:We assume that the detected symptomatic or severely symptomatic infected individuals will no longer involve into the infections because they are isolated and move to the hospitalized class. Thus, the infectious individuals belonging to pre-symptomatic (P), asymptomatic (A) and symptomatic undetected individuals (U) spread the disease. We consider that the COVID-19 pandemic situation usually persists for a short time period, and thus, we neglect the demographic factors (that is, birth, death and immigration) in our model. The parameters , and represent the probability of disease transmission rate due to contact between a uninfected individuals subject to an asymptomatic, a pre-symptomatic and an undetected symptomatic infected individuals, respectively. These parameters can be recast by social distancing strategies, for example, lockdown, closing educational institutions and remote working. The exposed individuals develop to infectious class at the rate , implying the average time spent in the exposed class is days. We assume that a fraction of individual move from the exposed compartment (E) to the infectious compartments (A) and (P), respectively. Pre-symptomatic infected individuals develop to infectious class (undetected symptomatic and severely symptomatic) at a rates and , and thus, the average time spent in the pre-symptomatic compartment is days. Isolated or hospitalized infectious individuals either recover at a rate or dead at a rate , indicating the mean infectious period is days. The mortality rate (subject to the life-threatening symptoms) can be mitigated by means of improved therapeutics. The parameters , , and indicate the recovery rate for the four compartments of infected individuals, namely isolation, asymptomatic, symptomatic undetected and pre-symptomatic, respectively. They may differ substantially if an appropriate treatment for the COVID-19 disease is known and adopted for the diagnosed patients. The interactions among the nine stages of COVID-19 compartmental mathematical model are shown in Fig. 1.

Fig. 1

Graphical schematic diagram represents the interactions among different compartments of infection in the mathematical model SEAPUDIRH: S, susceptible or uninfected; E, exposed; A, asymptomatic; P, pre-symptomatic; U, symptomatic undetected; D, symptomatic detected or severely symptomatic; I, isolation or hospitalization; R, recovered or healed; H, dead or extinct

In our model formulation, we ignore the probability of becoming susceptible individual again, after having already recovered from the COVID-19 infection, as this appears to be ignorable based on early evidence [52].

Analysis of the mathematical model

Our SEAPUDIRH model system (1) is a bilinear system of nine ordinary differential equations. The model system is positive for all as all the nine state variables are initialized at time 0 with nonnegative values. It can be noted that R(t) and H(t) are the cumulative state variables, which depends only on the other state variables and their own initial values.

The SEAPUDIRH model is compartmental model and shows the mass conservation property, which can be checked by summing all the state variables as follows:and thus, the sum of the states (total population) is constant, say N, that is,including the dead class. For the given initial values S(0), E(0), A(0), P(0), U(0), D(0), I(0), R(0), H(0) summing to N, we can show that the state variables converge to an equilibrium statewith . Therefore, only the uninfected, the recovered and the death individuals are eventually present, indicating that the epidemic scenario is no more. Thus, the possible equilibria are given by , with .

To understand the dynamics of the SEAPUDIRH system (1), we divide the entire population into three subclasses: The first subclass incorporates just single variable S (representing to the uninfected population), the second subclass incorporates E, A, P, U, D and I (the infectious population) that are nonzero only during the transient, and the last subclass incorporates R and H (corresponding to recovered and dead class). We focus on the infected (second) subsystem that can be denoted as the EAPUDI subsystem. It is worthy to mention that when (and only when) the infectious individuals become zero, the remaining state variables S, R and H are at an equilibrium. The uninfected variable S is monotonically decreasing and converges to , whereas the state variables R and H are monotonically increasing and converge to their asymptotic values and , respectively, which happen if and only if E, A, P, U, D and I converge to zero.

The whole system can be rewritten in the feedback structure, in which the EAPUDI subsystem can be shown as a positive linear system subject to a feedback signal f as follows. Denoting , we can recast the EAPUDI subsystem aswherewhere and . The rest of the three variables satisfy the following differential equationsThe time-varying uninfected individual S(t) ultimately converges to a constant value , and we can proceed with a parametric investigation with respect to the asymptotic convergence of . An important characteristics is stated in the following theorem.

Theorem 1

The EAPUDI subsystem with uninfected individuals is locally asymptotically stable if and only if

Proof

Since at early stages and all the other individuals E, A, P, U, D, I, R, D N, one can conduct a linearization of the above system of equations. This informs us about the early stage growth of the COVID-19 outbreak, in particular the exponential growth rate. The variational matrix of the linearized system around the equilibrium point is given byThe above Jacobian matrix has three zero eigenvalues, and the six eigenvalues are the roots of the given polynomialwhereThe transition function from f to for the system of equations (2)-(6) is Since the SEAPUDIRH system (1) is positive, thus the norm of is equal to

Therefore, by the standard root locus (small gain argument) on the positive SEAPUDIRH system , we can conclude that the polynomial is Hurwitz (as all the roots lie in the left-hand plane) if and only if the expression (8) holds, where , which proves the theorem.

It can be observed that, therefore, we define the basic reproduction number and we are well justified to define the basic reproduction numberand stability of the steady state occurs for .

Sensitivity analysis of basic reproduction number

To determine how best to mitigate the human mortality and morbidity due to COVID-19, it is very important to understand the importance of multiple factors responsible for the transmission of coronavirus. At the beginning, transmission of the virus is directly associated to the basic reproduction number . The main reason for sensitivity analysis is to investigate the robustness of model predictions to the system parameter values, as there are errors in collected data and presumed system parameter values. We perform the sensitivity analysis to delineate the impact of parameters which are related to the basic reproduction number , and hence the impact of those parameters on the system dynamics. Also, sensitivity analysis gives an idea about the intervention strategies.

More precisely, when a system parameter alters, then the sensitivity indices allow us to quantify the relative change in a variable. In order to do the sensitivity analysis, we compute the normalized forward sensitivity index of and these indices quantify the relative change in with respect to the relative change in its parameter. Thus, the normalized forward sensitivity index of [51], with respect to the parameter , is given byIn similarly way, we can find the sensitivity indices of with respect to the other system parameters.From the above expression, we obtain that , if and if . We also obtain that . Notice that does not depend on the system parameters , so .

Since and are the probability of disease transmission rate, , which implies that the increment of and will cause to increase. Here, and , which implies that increment of rate of pre-symptomatic to symptomatic undetected and fraction of asymptomatic carriers will cause to increase. Again, , so the increment of recovery rates and will cause to decrease. Moreover, the increment of the rate at which pre-symptomatic infected individuals become symptomatic detected will cause to decrease as . These conclusions are biologically acceptable.

From the expressions of , and , it is easy to show that , . Here, , which implies that the perturbations in the parameters and produce equal and opposite change in . Also here, , which implies that the perturbations in the parameters and produce equal and opposite change in . Again, , which implies that the perturbation in produces a relatively low changes in comparing to the perturbation in . The sensitivity index becomes most positive if the following inequality holds:

Table 1

Description of the SEAPUDIRH model parameters used for numerical simulations

Parameters	Description	Values & Unit	References
β_p	Disease transmission rate due to contact between susceptible and pre-symptomatic	0.3649 day^-1	Estimated
β_a	Disease transmission rate due to contact between susceptible and asymptomatic	0.2788 day^-1	Estimated
β_u	Disease transmission rate due to contact between susceptible and undetected symptomatic	0.3690 day^-1	Estimated
γ_e	Transition rate from exposed to infectious	0.15 day^-1	[16]
ρ	Fraction of asymptomatic carriers	0.67 Dimensionless	Estimated
γ_a	Rate at which asymptomatic infected individuals become recovered	0.1428 day^-1	Estimated
γ_p	Rate at which pre-symptomatic infected individuals become symptomatic detected	0.93575 day^-1	Estimated
ξ_p	Rate at which pre-symptomatic infected individuals become symptomatic undetected	0.01425 day^-1	Assumed
η_p	Transition rate of pre-symptomatic to recovery	0.0112 day^-1	[16]
γ_u	Recovery rate of undetected symptomatic infected	0.428 day^-1	Assumed
γ_d	Transition rate of symptomatic detected to isolation or hospitalization	0.736 day^-1	[13]
γ_i	Transition rate of isolation to recovery	0.008 day^-1	[10]
ξ_i	Mortality rate of isolated infected subject to the life-threatening symptoms	0.000037 day^-1	Assumed

The sensitivity index becomes most negative if the following inequality holds:Hence, the absolute value of the normalized forward sensitivity index becomes largest for the parameters and . Positive sensitivity indices imply that is an increasing function of the corresponding system parameter and negative indices indicate that is a decreasing function of that system parameter. As for example, indicates that if is increased by 10%, then the is also increased by 9.09% Table 1. Again, for indicates that 10% increment in will decrease by 9.09%. From Table 2, it can be observed that the disease transmission rate and the fraction of asymptomatic carriers are most sensitive parameter which has positive impact on and the parameters and show more negative impacts on Fig. 2 shows sensitivity of the model parameters with respect to .

Table 2

Sensitivity indices of evaluated at the baseline parameter values listed in Table 1

Parameters	β_a	β_p	β_u	ρ	γ_a	γ_u	ξ_p	η_p	γ_p
Sensitivity index	0.9099	0.0872	0.0094	0.7271	-0.9099	-0.0089	0.0056	-0.0081	-0.0877

Fig. 2

The plot represents the normalized forward sensitivity indices of the basic reproduction number with respect to the baseline system parameter values specified in Table 1

Optimal control problem

The control measures play an important role in mitigating the transmission of the COVID-19 dynamics. Thus, in this section we formulate an optimal control problem using and as time-dependent control variables corresponding to the mathematical model (1). We are mainly interested to investigate the effect of these intervention strategies on the transmission of the COVID-19 virus. We minimize the infected population and the hospitalized population by using the theory of optimal intervention strategies [53]. We performed the theoretical analysis as well as numerical illustrations to show how these control strategies make an evident on the transmission of COVID-19 pandemic and to minimize the disease burden with their implementations. The definition of the parameters and supposition leads to the model (1) which implies a coupled system of nonlinear differential equations with nine state variables, that is, S(t), E(t), A(t), P(t), U(t), D(t), I(t), H(t) and R(t). In our control problem, we implemented two control variables that externally control the number of isolated or hospitalized cases and clinically ill or symptomatic infected cases over a specified time window. Now, we delineate these two control strategies and then investigate their impact on COVID-19 transmission dynamics.

To provide better treatment to the clinically ill or infected individuals

Better treatment strategy to the symptomatic infected population not only weakens the disease prevalence but also influences its development. It is very essential to the individuals to consult with a doctor if the mildly symptom is observed in this body. The individuals should not neglect by considering it as a mildly cases. Consulting to the doctors at the beginning can reduce the disease transmission, and thus, a saturated isolation rate is implemented in our original model (1) where represents the rate at which clinically ill cases move to the hospitals for medical treatment without neglecting the symptoms with magnitude and half saturation constant . There are different costs related to the medicines, kits, diagnosis, health protocols, etc., during the period when an individual detected with coronavirus symptoms is taken into account. Hence, we consider that the treatment intensity as a control variable lies between 0 and 1 where 0 represents case when an individual ignores his/her mild symptoms and 1 indicates the case when an individual consults with doctor neglecting the mild symptoms.

Better treatment strategy to the isolated individuals

The disease mortality can be lower by providing proper antivirals or vaccination against novel coronavirus to the isolated or hospitalized individuals at the very beginning of infection. It influences the development of the diseases too. It is to be noted that the vaccine is accessible and given to the population who are admitted in the nursing home or hospitals in a limited amount. As the accessibility of resources associated with the medical treatments, financial crisis, disease diagnosis, etc., all the equipment is restricted. To keep this mind, we implemented intervention strategy in our original model (1) by introducing a saturated treatment rate function where represents the treatment rate with magnitude along with half saturation constant . Different costs related to the medicines, health problem, isolation, vaccination, medical kits, etc., at the time of medication period are need to be taken into account. Thus, we consider as a control variable lying between 0 and 1, where 0 represents no response and 1 represents full response after the treatment.

Our main aim is to obtain the optimal treatment strategy for the clinically ill population and isolated individuals with minimum cost. To do this, we consider the admissible set for two control variables and defined as follows:After introducing the controls in the original system (1), we get the following system of equations:subject to the minimize the objective functionalHere, J is the total cost and at any time t, the integrand represents the current value of the cost. The positive parameters and represent the weight constants, which balance the units of the integrands [44, 46, 54].

Existence of optimal control

In this section, the existence of optimal solution of the system (10) with objective functional (11) is verified. We assure the existence of the control pair and inwhich minimizes the cost functional J.

Theorem 2

There exist optimal controls and in corresponding to the control problem (10) and (11) on a fixed interval [0, T].

We will use Lukes [53] results to proof the theorem. In order to ensure the existence of optimal controls, the conditions listed below must be satisfied: The integrand in the cost functional J, is convex on the control set . The exist constants and such tat where depends upon the upper bounds of D and I and .

The boundedness of the solution of the control system (10) confirms the existence of the solution of control system (10).

The optimal control set and the corresponding state variables are non-empty.

The solution of the control system (10) is bounded above by a linear function in the control as well as state variables.

Characterization of optimal control

Now, we need to estimate the value of optimal controls and such thatwhereThe Lagrangian is given byThe Hamiltonian is defined as follows:where the adjoint variables can be obtained as the solution of the following system of differential equationssatisfying , for , that is the transversality conditions. Now, using Pontryagin’s maximum principle [53, 54] we characterize the optimal control pair and as follows:

Theorem 3

Let and be optimal control variables and and be corresponding optimal state variables of the system (10)-(12). Then, there exists adjoint variable that satisfies the canonical equations (15) with transversality conditions , for . The corresponding optimal controls and are given as,and

Let and be the given optimal control variables and and be corresponding optimal state variables of the system (10) which minimize the cost functional (11). Hence, by the Pontryagin’s maximum principle, there exist adjoint variables and , which satisfy canonical equations (15) with transversality conditions , for . Now, using the optimality condition, we obtain

at and at .

Thus, we obtain and .

Hence, using the above results along with the characteristics of control set , we obtainandwhich can be equivalently written as (16) and (17).

We numerically solved the control system (10) subject to the objective functional with system parameter values which are listed in Table 1, and we assume the initial population as  A(0),  P(0),  U(0),  = . We use forward–backward sweep method to solve the optimality system numerically. The state equations are solved forward in time and after that the adjoint state system solved backward in time using Runge–Kutta fourth-order iterative scheme. The values of the controls are updated at each iteration. The process is repeated to reach a desired convergence criteria. We have compared our model with the implementation of optimal control strategy and without administration of optimal control policy, and the variation of the state variables in the presence of optimal controls and in the absence of optimal control is shown in Figs. 3, 4 and 5. It is observed that there is no significant effect of the optimal controls for the S, E, A, P and U populations.

Fig. 3

The figure represents the comparison of the corresponding susceptible (S), exposed (E), asymptomatic (A), pre-symptomatic (P), symptomatic undetected (U), clinically ill or symptomatic detected (D), isolated or hospitalized (I), recovered (R) and dead (H) classes without intervention strategies, with the implementation of intervention strategies (only control). Optimal treatment strategy (solid blue line) demonstrates reduction of the (D), (I) and (H) classes and increase in the (R) class when compared with the no controls (solid black curves). The parameter values are given in Table 1

Fig. 4

The figure represents the comparison of the corresponding susceptible (S), exposed (E), asymptomatic (A), pre-symptomatic (P), symptomatic undetected (U), clinically ill or symptomatic detected (D), isolated or hospitalized (I), recovered (R) and dead (H) classes without intervention strategies, with the implementation of intervention strategies (only control). Optimal treatment strategy (solid blue line) demonstrates substantial reduction of the (D), (I) and (H) classes and substantial increment of the (R) class when compared with the no controls (solid black curves). The parameter values are given in Table 1

Fig. 5

The figure represents the comparison of the corresponding susceptible (S), exposed (E), asymptomatic (A), pre-symptomatic (P), symptomatic undetected (U), clinically ill or symptomatic detected (D), isolated or hospitalized (I), recovered (R) and dead (H) classes without intervention strategies, with the implementation of intervention strategies (combination of the two controls and ). Optimal treatment strategy (solid blue line) demonstrates significant reduction of the (D), (I) and (H) classes and significant increment of the (R) class when compared with the no controls (solid black curves). The parameter values are given in Table 1

We employ two different optimal treatment strategies, namely and The solid black curves for epidemic (without administration of intervention strategies) are exhibited to highlight the difference from those generated via implementation of optimal treatment polices. The solid blue epidemic curves (with the implementation of single control only) are shown the reduction of the symptomatic detected cases, isolated cases and dead population when compared with the no intervention scenario. It can also be noted that the recovered population increases after the implementation of single optimal control when compared without administration of controls, that is, Thus, the implementation of single control is better in comparison with the without control policy. The impact of the time-dependent optimal control only [ here is shown in Fig. 3. This numerical simulation shows significant decrease in the density of I, H and D populations and significant increase in R population.

Again, the impact of the time-dependent optimal control only [ here is plotted in Fig. 4. The solid blue epidemic curves (with the implementation of single control only) show the reduction of the symptomatic detected cases, isolated cases and dead population when compared with the no intervention scenario (solid black curves). It can also be noted that the recovered population increases after the implementation of single optimal control when compared without administration of controls, that is, From the numerical simulation (see Fig. 4), we can infer that the effect of the optimal control is more effective than the optimal control .

The time series Fig. 5 represents the impact of each treatment strategies on the disease state variables in the presence of the optimal control strategy and in the absence of optimal treatment strategy. These plots demonstrate the cases under no control , with the combination of two different controls . The solid black curves for epidemic (without implementation of intervention strategies) are shown to highlight the difference from those generated through the administration of optimal treatment polices. The solid blue epidemic curves (with the administration of the combined of two controls) exhibit the significant reduction of the symptomatic detected cases, isolated or hospitalized cases and dead individuals when compared with the no intervention strategies. The reduction of the symptomatic infected (detected) class (D) becomes maximum when the optimal controls and are applied simultaneously. It is worthy in mentioning that the recovered individuals (R) increased substantially when compared with the no control cases as well as the implementation of single controls and Implementation of the combination of two controls is better when compared with the no control policy as well as the single control policy. Thus, from the numerical simulations (Fig. 5) we can infer that the employment of intervention strategy has an effect in controlling the spread of novel coronavirus epidemic. Also, we can conclude that the combination of two controls is more effective when compared with the implementation of the single control.

The time evolution of the control variables and is shown in Fig. 6 with respect to the disease transmission rate due to contact between susceptible and asymptomatic population. Figure 7 shows the time evolution of the control variable with respect to the parameter (left panel) and the time evolution of the control variable with respect to the parameter (right panel). The time evolution of the control variable and is shown in Fig. 8 with respect to the rate at which asymptomatic infected individuals become recovered.

Fig. 6

Time profile of the optimal controls (left panel) and (right panel) for different values of . Rest of the parameters are defined in Table 1

Fig. 7

Time profile of the optimal controls (left panel) for different values of and (right panel) for different values of

Fig. 8

Time profile of the optimal controls (left panel) and (right panel) for different values of

We designed an analysis with only optimal control policies, that is, implementing of optimal control strategy on isolated or hospitalized individuals. Numerical simulations of the proposed COVID-19 model, collectively SEAPUDIRH, have been shown with the application of only optimal control and without implementation optimal control in Fig. 3. Implementation of only control decreases the symptomatic detected individuals (D), isolated or hospitalized population (I), dead cases (H) and consequently the recovered class (R) has been increased. More people become recovered, and less number of people become infected after the successful implementation of only optimal control. Also, the number of death cases decreased after the introduction of only optimal control.

Next we implement the single optimal control for the saturated treatment rate with . Here, the value of lies between 0 and 1, where 0 indicates no response for COVID-19 medical treatment, while 1 indicates full response for COVID-19 medical treatment. Numerical simulation of the proposed COVID-19 model, collectively SEAPUDIRH, has been shown with only optimal control and without implementation of optimal control in Fig. 4. Implementation of only optimal control decreases the symptomatic detected (D) cases, isolated or hospitalized (I) cases, dead classes (H) and increases the number of recovered (R) population. More population become recovered, and less number of population become infected after the successful implementation of only optimal control. Here, also the number of death cases reduced after the implementation of only optimal control.

Now, we implement both the optimal controls and simultaneously in our proposed SEAPUDIRH model. Numerical solution of the proposed COVID-19 model, collectively SEAPUDIRH, is shown with both and optimal controls in Fig. 5. Significant decrease in the clinically ill or symptomatic detected (D) cases, isolated or hospitalized (I) individuals, dead classes (H) is observed and significantly increases the recovered (R) individuals. The number of clinically ill or symptomatic detected (D), isolated (I), dead (H) classes becomes less, while the optimal controls and are applied simultaneously when compared to the implementation of single optimal controls.

Numerical simulations

In this section, we performed some numerical illustrations for the COVID-19 pandemic in India based on the SEAPUDIRH model system (1). We estimated the SEAPUDIRH model parameters based on the daily observed COVID-19 cases, cumulative number of COVID-19 cases and the dead cases due to SARS-CoV-2 virus infection, for the time period March 01, 2020, to August 15, 2020. The data were collected from the WHO situation report [2]. We estimated six parameter values, namely , , , , and among 13 system parameters. The estimated value of those six system parameters is listed in Table 1. Based on the estimated model parameters and other parameters, we performed short-term prediction of the epidemic. We also performed some mathematical analysis to get an idea about the basic reproduction number . We also showed that our system is asymptotically stable for .

Preliminary insight from COVID-19 observed data

First, we perform a simple mathematical analysis that can be endeavored to obtain some insight into the COVID-19 pandemic and to generate iterative time-lag maps. The motivation is to build a relationship between some individuals at time (in days) and the same individuals at day n, corresponding to the time lag of k days. Thus, we have investigated the relation between the populations of and days for the daily new COVID-19 infection cases and daily new deaths due to COVID-19 infection with a time lag of t days. Now considering , we have graphically shown the days population along horizontal axis and days population along the vertical axis (see Fig. 9). Both the populations follow the power lawThe values of the parameters a and b are 1.192 and 0.9842, respectively, for the daily new COVID-19 cases (D) with coefficient of determination = 0.9928. Also, we consider the values of the parameters a and b are 1.821 and 0.8954, respectively, for the daily new deaths (H) with coefficient of determination = 0.8950.

Fig. 9

The figure represents the recurrence plots for the daily new COVID-19 cases and the daily new deaths due to COVID-19 infection. The best fit of the power law of the kind (18) is shown by solid blue curve. The coefficient of determination for the daily new COVID-19 cases is 0.9928, and the coefficient of determination for the daily new deaths is 0.8950. The best fit curves of the power law are increasing for both the populations

Model calibration

We have calibrated our proposed COVID-19 model system (1) with the daily observed COVID-19 infection cases, cumulative number of COVID-19 cases as well as the number of daily deaths due to COVID-19 infection. The data of daily new COVID-19 infection cases, cumulative number of COVID-19 infection cases, daily deaths and cumulative number of deaths due to COVID-19 infection in India were obtained from the situation report of World Health Organization (https://www.who.int) [2]. First COVID-19 cases have been observed in India on January 30, 2020, but March 01, 2020, onward, the daily new COVID-19 infection cases were reported continuously. Hence, we have collected the COVID-19 data from March 01, 2020, to August 15, 2020. Model simulation depends on both the system parameter values and the initial population sizes. As the population in India is and the 3 persons were infected by COVID-19 as on March 01, 2020, we have considered the initial population sizes as = and varied the system parameter values to fit with the observed data. The model simulations were fitted with the observed daily new COVID-19 cases, cumulative number of COVID-19 cases as well as daily deaths due to COVID-19 infection, using the least square method by minimizing the sum of errors. Six model parameters, namely , , , , and , have been estimated by this least square method out of the thirteen system parameters [50]. The value of the estimated parameters is listed in Table 1, and the best fitted curves are shown in Fig. 10. The daily new COVID-19 infection cases, cumulative number of COVID-19 cases and cumulative number of deaths have shown in first, second and third row, respectively. Filled red circles are the observed cases, and the blue curves are best fitted model simulation. Our SEAPUDIRH model system (1) nicely captured the increasing trend in all the populations. To compute the accuracy of model simulation, we have computed the root-mean-square error and mean absolute Error , which can be defined as follows:where R(i) represents the observed COVID-19 infection cases, M(i) represents the model simulation on the day, and is the sample size of the observed COVID-19 data. The values of and are 3266.6 and 2204.0, respectively.

Fig. 10

The figure shows the model simulations fitted with the observed daily new COVID-19 cases, cumulative confirmed COVID-19 cases and cumulative deaths due to COVID-19 infection in India. Observed cases are shown by the solid red circle, and the best fitting curve for the SEAPUDIRH Model system (1) is shown by blue curve. The first row represents the daily new COVID-19 cases, the second row represents the cumulative number of confirmed COVID-19 cases, and the third row represents the cumulative number of deaths due to COVID-19 infection. The initial values are used to solve the system of ordinary differential equations = . The estimated parameter values are listed in Table 1

Impact of basic reproduction number

An important identifier for the transmissibility of infectious diseases such as SARS-CoV-2 virus is the basic reproduction number , which can be defined as the average number of secondary infections created by a single infected person, in a population whose all members are uninfected or susceptible. The numeric value of the basic reproduction number, , is computed using the estimated system parameter values, and the system parameter values are listed in Table 1, and for our proposed SEAPUDIRH model, we obtained , which shows the substantial outbreak of COVID-19 in India.

Figure 11a represents the basic reproduction number of the SEAPUDIRH system (1) with respect to (contact rate between susceptible and asymptomatic) and (contact between susceptible and pre-symptomatic). The value of increases if any one of the parameters or increases. The disease-free equilibrium point of the model system (1) is asymptotically stable for , and for lower values of or , the reproduction number is less than unity. Thus, the coronavirus outbreak in India can be controlled by reducing contact rate between the individuals, which can be done by maintaining social distancing and contact tracing. Figure 11c shows that if both and increase, then the value of will also increase significantly.

Fig. 11

The figure shows the basic reproduction number of the system (1) in terms of a and , b and , c and d . In the figures (c) and (d), green shading region indicates and red shading region indicates . Parameters value are , , , , , , , , , , , and . The black dashed contours in the figures (a) and (b) are the contours for

Figure 11b represents the basic reproduction number of the SEAPUDIRH system (1) with respect to (the rate at which pre-symptomatic infected person becomes symptomatic detected) and (the rate at which asymptomatic infected individuals become recovered). For larger value of , the value of becomes less than 1 if increases; that is, if the rate of conversion of pre-symptomatic infected individuals (P) to symptomatic detected (D) is higher, then higher rate of recovery makes the disease-free steady-state stable. But for a lower value of , if increases, then the value of remains greater than unity, that is, for lower rate of conversion from pre-symptomatic infected individuals (P) to symptomatic detected (D) will not the disease-free equilibrium even if the rate of recovery is high enough. This suggests that, to make the system disease free, identification of pre-symptomatic infected individuals (P) is utmost important which can be done by increasing the number of COVID-19 test. Figure 11d represents that for higher value of the value of becomes less than one if increases and the disease-free steady state becomes stable. Thus, to mitigate the coronavirus disease from the system, the value of need to increase, which can be done by increasing the number of COVID-19 test.

The time evolution of the proposed mathematical model system (1) is shown in Fig. 12. The system parameters are chosen in such a way that for the first column, for the second column and for the third column of Fig. 12. The time series solution in the first column of Fig. 12 shows that the disease-free equilibrium point of the system (1) is stable when . The third column in Fig. 12 shows that the endemic equilibrium point of the system (1) is stable when . The second column shows the time evolution of the system (1) for , which indicates that there is no exponential growth of the number of active cases with time. However, it can still be growing at a slower than exponential rate, for example, linearly or logarithmically. Thus, we observed from our numerical simulation that the model system (1) is locally asymptotically stable around the disease-free equilibrium if and the endemic equilibrium point is stable if .

Fig. 12

The figure represents the numerical simulation of the system (1) for (first column), (second column) and (third column). Parameter values are specified in Table 1

Sensitivity analysis of basic reproduction number shows that the system parameters , and are most sensitive with reference to . The time evolution of the system (1) is shown in Fig. 13 for different values of the system parameters and . Other system parameters for the numerical simulation are same as in Table 1. Here, for , ; for , ; for , ; for , . Thus, the number of infected cases is more for (or ) but lower value of . Therefore, if the disease transmission rate () (contact between susceptible and asymptomatic individuals) remains constant but the rate () at which asymptomatic infected individuals become recovered is decreased, then the number of infected cases will increase with time.

Fig. 13

Time profile of the populations of the system (1) for different values of and . Other parameters are same as in Table 1

Again the time evolution of the model system (1) is shown in Fig. 14 for different values of the system parameters and . Other system parameters for the numerical simulation are same as in Table 1. Here, for , ; for , ; for , ; for , . Here, the number of infected cases become more for (or ) and for the higher value of . Thus, if the disease transmission rate () (contact between susceptible and asymptomatic individuals) remains constant but fraction of asymptomatic carriers () is increased, then the number of infected cases will also increase. In this scenario, we noticed that the most effective parameter of the system (1) is the disease transmission rate () and has an impact in controlling the basic reproduction number .

Fig. 14

Time profile of the populations of the system (1) for different values of and . Other parameters are same as in Table 1

Short-term prediction

As on August 15, 2020, the total number of COVID-19 infected cases is 25,26,192 and 49,036 people died due to COVID-19 infection. To control the transmission of the COVID-19 disease, Govt. of India has implemented nationwide four phases of lockdown and two phases of unlock. Due to the absence of licensed vaccine, antivirals and any therapeutics, Govt. is taking lockdown and unlock policies to control the community level transmission. Infection disease mathematical modeling provides an efficient prediction of the COVID-19 outbreak for short time duration. But, due to an inaccurate of data, the prediction based on the estimated data and model assumptions may not always be true. This predictions are essential to take necessary administrative strategies by the government and some other non-governmental agencies. The second phase of unlock starts from July 01, 2020, and the second phase unlock will continue up to July 31, 2020. In this study, we have considered COVID-19 data from March 01, 2020, to August 15, 2020, and our model provides short-term prediction for 21 days from August 16, 2020, to September 05, 2020 (see Fig. 15 for short-term prediction). The black dashed curve in Fig. 15 is the predicted cases, first row predicts the daily new COVID-19 cases, second row predicts the cumulative number of confirmed COVID-19 cases, and third row predicts the total deaths due to COVID-19 infection in India. Our model predicts 138,744 daily new cases, 4,676,105 cumulative cases and total 103,729 deaths on September 05, 2020. This prediction gives us an estimation of the COVID-19 outbreak in India. This prediction highly depends on the system parameter values and the initial population sizes. If the administrative strategies changed or effective vaccine or therapeutics introduced, the prediction may differ from the actual cases.

Fig. 15

The figure represents the short-term prediction (21 days from July 01, 2020, to July 23, 2020) of the model (1) for the daily new COVID-19 cases (first row), the cumulative number of COVID-19 cases (second row) and the cumulative number of deaths due to COVID-19 infection (third row) in India. Our model predicts increasing trends for all the classes

Conclusion

In this article, we have analyzed the COVID-19 pandemic data in India made available to the researchers by the WHO situation report [2] and referring to the time period March 01, 2020, to August 15, 2020. Everyday the total number of daily confirmed cases rises with more than 2,500 reported from several states and union territories in India. This is a worrying situation as with a second most populated country in the world, and within few days, India will enter in community-wide transmission of COVID-19 epidemic. Due to the absence of any effective treatment, antivirals or licensed vaccine and with an incomplete understanding of epidemiological traits, predictive compartmental mathematical models can aid in understanding of both the disease spread and its control [13, 17].

In our present study, we have formulated and analyzed a compartmental epidemiological model of SARS-CoV-2 viruses to control and forecast the pandemic. We have computed the basic reproduction number to get an idea under what condition the disease will be die out or spread into the population. We also showed that whenever , the disease-free steady state will be asymptotically stable. In our proposed model, we have considered the transmission variability between asymptomatic and pre-symptomatic infected persons to better understand the dynamics of the COVID-19 epidemic [18]. Also, we have investigated the sensitivity indices of the threshold quantity and obtained that the disease transmission rate is the most sensitive system parameter, which has a positive influenced on the basic reproduction number .

We have calibrated our proposed mathematical model to fit the daily confirmed new COVID-19 cases, cumulative number of confirmed COVID-19 cases and the cases of daily new deaths due to COVID-19 infection of India. Six model parameters are estimated from the best curve fitting of our SEAPUDIRH model. Using the estimated system parameter values, we then obtain the basic reproduction number , which shows the substantial outbreak of COVID-19 pandemic in India. The policymakers and health care agencies should focus on the successful implementation of control mechanisms to mitigate the burden of the COVID-19 diseases. From the contour plot (see Fig. 11), it can be concluded that to mitigate the disease transmission must maintain the social distancing and contact tracing by implementation of nationwide lockdown, closing of educational institutions, etc.

The calibrated model is then used for short-term prediction of India for the time period August 16, 2020, to September 05, 2020 (see Fig. 15). However, the increasing pattern of daily newly COVID-19 cases, cumulative number of COVID-19 cases and the deaths cases of COVID-19 is well captured by our model for India. From the short-term prediction, we can conclude that 138,744 daily new confirmed cases, 4,676,105 cumulative number of confirmed cases and the total 103,729 deaths on September 05, 2020. It is worthy to mention that here we can predict the nature of the SARS-CoV-2 for the short-term period as the Govt. strategy can be changed in time to time resulting in the interrelated changes in the associated system parameters of our SEAPUDIRH model.

Moreover, we have modified our SEAPUDIRH model system by implementation of the theory of optimal control. Implementation of intervention strategies helps to mitigate the novel coronavirus burden. To study the impact of non-pharmaceutical and pharmaceutical intervention strategies and to reduce/stop the community transmission of the novel coronavirus, we implemented two different control measures, namely and . Analytically, we investigated the existence of the optimal control functions. Furthermore, to mitigate the clinically ill populations and to give the better treatments for isolated or hospitalized individuals and to minimize the cost of the control functions, an objective functional has been formulated and interpreted with the aid of Pontryagin’s maximum principle.

The behavioral changes in the susceptible population alter with time, and to prevent or limit the rate of interaction between susceptible (S) and isolated (J) individuals, we put the external effort as an intervention strategy. Also, the clinically ill or infected individuals become aware about the coronavirus disease fatality and consult with medical officers or admit to the hospitals or nursing home if the symptom occurs and thus we use the control that alters the fraction of infected cases. Again, we apply the controlling effort that taking into account the cost regarding medical treatments, use of test kits and life-saving medicines is provided to the hospitalized or isolated individuals. Thus, we considered all these two control policies in our proposed model.

Finally, through an extensive numerical illustrations, we conclude certain observations. Numerical solutions are demonstrated in Fig. 3 that the control strategy represents clinically ill or infected cases efforts with its intensity for quite a long period of time and then it reduces. Later on, the control indicates better treatment strategy for isolated individuals efforts with high magnitude nearly about 300 days following a reduction at a later stage, though this magnitude is higher than the magnitude of without control strategy. From the extensive numerical simulations (see Fig. 5), we can conclude that the implementation of intervention strategy has an important effect in controlling the spread of novel coronavirus pandemic. Also, we can conclude that the combination of two controls is more effective when compared with the implementation of the single control as well as without control. Thus, introducing the combination of two control strategies may help to reduce the spread of novel coronavirus disease at this present pandemic scenario. We plan to expand our COVID-19 modeling by introducing vaccination and the effect of environmental contamination. Most probably, this COVID-19 epidemic will remain in the human population for many years, so the changes we see in public health services and in human habits have to go on over time.