Mathematical modeling and optimal intervention strategies of the COVID-19 outbreak.

34,354,966 active cases and 460,787 deaths because of COVID-19 pandemic were recorded on November 06, 2021, in India. To end this ongoing global COVID-19 pandemic, there is an urgent need to implement multiple population-wide policies like social distancing, testing more people and contact tracing. To predict the course of the pandemic and come up with a strategy to control it effectively, a compartmental model has been established. The following six stages of infection are taken into consideration: susceptible (S), asymptomatic infected (A), clinically ill or symptomatic infected (I), quarantine (Q), isolation (J) and recovered (R), collectively termed as SAIQJR. The qualitative behavior of the model and the stability of biologically realistic equilibrium points are investigated in terms of the basic reproduction number. We performed sensitivity analysis with respect to the basic reproduction number and obtained that the disease transmission rate has an impact in mitigating the spread of diseases. Moreover, considering the non-pharmaceutical and pharmaceutical intervention strategies as control functions, an optimal control problem is implemented to mitigate the disease fatality. To reduce the infected individuals and to minimize the cost of the controls, an objective functional has been constructed and solved with the aid of Pontryagin's maximum principle. 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. Extensive numerical simulations show that the implementation of intervention strategy has an impact in controlling the transmission dynamics of COVID-19 epidemic. Further, our numerical solutions exhibit that the combination of three controls are more influential when compared with the combination of two controls as well as single control. Therefore, the implementation of all the three control strategies may help to mitigate novel coronavirus disease transmission at this present epidemic scenario. Supplementary Information: The online version supplementary material available at 10.1007/s11071-022-07235-7.

Introduction

As of November 06, 2021, the COVID-19 pandemic has already crossed 250,227,541 infected cases and more than 5,059,368 deaths throughout the world, becoming the considerable public health problem facing after the Second World War. It has a great impact to forecast the trend of SARS-CoV-2 pandemic circumstances. At this highly advanced stage, it is of utmost priority to reply some questions immediately: when the COVID-19 infection can be entirely removed by proper vaccination or immunization, and what is the interaction among the threshold of vaccination and disease reproduction number [1]. The SARS-CoV-2 has an important research as it was identified as pandemic by the World Health Organization (WHO) on March 11, 2020, and obtaining a lot of recognition of researchers/scientists in different areas of people knowledge.

Dynamical viewpoint is a fascinating outlook to deal with the bio-medical systems. Several literature presents different outlooks regarding dynamical models for an epidemic [2-4], generally utilized to delineate various infectious diseases. A fascinating and convenient outlook is relied on the population dynamics, where various individuals are employed to indicate the diseases, representing their evolution and interplays. Kermack and McKendrick [5] were one of the pioneers in developing three individuals in epidemic dynamic models, namely susceptible, infectious and recovered (SIR). Their mathematical model describes the course of epidemic exhibiting that it never be terminated necessarily by all populations becoming infectious. Anderson and May [6] modified the well-established SIR model developed by Kermack and McKendrick [5], incorporating an additional compartment: exposed class (E). Presently, most of the epidemic model is relied on the susceptible–exposed–infected–removed (SEIR) structure, widely used to delineate various infectious diseases as Hepatitis, HIV, HTLV, Dengue, Typhoid, among others.

The first case of COVID-19 was identified at Thrissur district, Kerala in India on January 30, 2020, when a student came from Wuhan, China. Due to the India COVID-19 tracker [7], there are 34,285,612 total confirmed cases, 33,661,339 recovered people and 458,470 total deaths are recorded in the Republic of India on November 06, 2021. The Govt. of India announced to maintain social distance or to adopt self-quarantine policy in order to avoid community transmission of the novel coronavirus diseases among the population. On March 22, 2020, the Indian Govt. has implemented a “Janata Curfew” during 14 h, while the total number of infected cases became more than 500 individuals. After that the Indian Govt. announced a nationwide lockdown for 21 days starting from March 25, 2020, in order to mitigate the spread of the novel coronavirus transmission. Later on, the Indian Govt. increased the nationwide lockdown to slow down the transmission of novel coronavirus pandemic. Due to the limited vaccine for COVID-19 diseases, the people must have to maintain social distance, wearing mask, avoid mass gathering or applying self-isolation to prevent the transmission of novel coronavirus [8]. The nationwide lockdown incorporates prohibition on individuals from stepping out of their houses, stop cultural functions, mass gathering, shutdown of all stores/shops excluding medicine shops, school/colleges, banks, hospitals, nursing home, etc., postponement of all the academic college/university/institution (allowed work from home only), postponement of all the private and public cars and also the closure of all cultural activities, political gathering, sports, entertainments, etc. This epidemic has significantly affected the life of all the people from economically as well as health. According to the report from the World Bank and the Reserve Bank of India (RBI), after the Second World War, this is the first time when the economic growth rate of India and the World as a whole will be diminished by 1.5–2.8% due to the effect of COVID-19 outbreak. Indeed, the expansion of this novel coronavirus outbreak has disrupted the daily life of the people, economic and the health of the people. Thus, this is one of the major issues for all the people how long this situation will continue and when the COVID-19 disease will be under maintained. In order to control the novel coronavirus diseases, the scientists around the world has been trying to improve the vaccine quality such that the disease can be eliminated.

Due to the fast expansion of novel coronavirus, which lead to necessary beds in medical colleges saturating in many countries [9], virtually many countries executed some form of non-pharmaceutical interventions (NPIs) [8, 10]. The aim of these intervention strategies was mainly to “flatten the curve” indicating that a mitigation of the coronavirus outbreak could be attained through a change of behavior in the human [11]. Flaxman et al. [10] quantifies the role of different non-pharmaceutical interventions on the transmission number, , utilizing real data from the 11 different countries of Europe that were dealing with the novel coronavirus disease during the first wave of infections. An important interpretation is that the close progression of the non-pharmaceutical intervention strategies make it hard to obtain a clear answer as to which non-pharmaceutical intervention strategy has the most impact. Also, non-pharmaceutical intervention strategy take sometime to be effective. Thus, since the real alteration of the expansion happens much later than the intervention dates, the cumulative alteration tends to become aggregated in latest possible time window. The outbreak of COVID-19 has presented unprecedented challenges for public health and disease prevention. The transmission of COVID-19 is very efficient, in which one patient could reportedly transmit the disease to 2–3.5 other patients at the early stage [12, 13].

Nowadays, we control the epidemic associated with the ability to generate meaningful data, which in turn will allow us to use same data for mathematical modeling purposes that are the best framework to deal with forthcoming scenarios [14]. Many attempts have been taken in this field to understand the dynamics of epidemiological cycle of the epidemic [12, 15–17]. The adjustment of the system parameters in a dynamic way, through the imposition of limits on the system in order to optimize a given function, can be implemented through the optimal control theory [18].

Different types of mathematical models have been developed and observed to study the transmission dynamics of the novel coronavirus epidemic with their dynamical behaviors [19-29]. Data-driven model for the COVID-19 transmission dynamics with the incorporation of distributed time delay has been investigated by Liu et al. [30]. Khyar et al. [31] proposed a multi-strain SEIR compartmental model to investigate the interactive dynamics of the novel coronavirus outbreak and showed that their model produces global dynamics. The theory of optimal control strategies has been utilized to investigate the COVID-19 disease transmission, and the authors showed that optimal doses are needed to curtail the spread of the novel coronavirus diseases [32-34]. Different types of qualitative behavior has been investigated in different mathematical models [35-37]. Data-driven modeling approach is a powerful tool to investigate the dynamics of the infectious diseases [38]. To understand the arithmetic optimization algorithm in case of COVID-19 pandemic model, reader can see the articles [39, 40]. Also, different kinds of mathematical modeling has been studied to better understand the complicated biological phenomena [41, 42].

The use of optimal control theory in infectious disease is very common: while mathematical modeling of the infectious diseases has expressed that amalgamations of hospitalization, isolation and vaccination are usually essential in order to obliterate infectious diseases, the theory of optimal control tells us how they should be managed, by giving the proper times for intervention with proper dosages [43]. This optimization technique has also been utilized in some research works within the span of novel coronavirus. The optimal control theory of a modified susceptible–exposed–infected–recovered (SEIR) mathematical model has been investigated with the goal to aim to study the effectiveness of two inflexible lockdown techniques on the United Kingdom population [44]. Another novel coronavirus case studies incorporate the utilization of the theory optimal control in the USA [45]. An optimal use of an hypothetical vaccine for novel coronavirus has been also studied [46]. Non-pharmaceutical intervention strategies can be used as a key parameter to maintain the active number of infected populations.

Here, we are interested in applying the theory of optimal control which has a tool to understand ways to curtail the progression/development of the novel coronavirus in a case study of India by designing optimal disease intervention techniques. Furthermore, we consider some significant questions that are not completely studied in the existing literature. Our model demonstrates the implementation of optimal control technique, to minimize the symptomatic infected class and hospitalized individuals where the response applicability of public health facilities is perpetuated. As the COVID-19 epidemic has exhibited, the public health concern is not only a medical problem, but also affects the country as a whole [34]. The qualitative behavior of the interactive dynamics of the novel coronavirus epidemic has also been investigated in this manuscript. We employed three control strategies to better understand the complicated dynamics of COVID-19 transmission and to mitigate transmission of COVID-19.

The organization of the article is as follows: In the next section, we describe the proposed SAIQJR model. In the same section, we describe the qualitative behavior of the proposed model incorporating positivity, boundedness, computation of basic reproduction number, normalized forward sensitivity analysis for , existence and stability of the biologically feasible equilibrium points. We also performed the global asymptotic stability of the disease-free equilibrium points. Moreover, we implement the optimal control strategy with three different controls and solve the control problem analytically in Sect. 3. In Sect. 4, we performed extensive numerical simulations without optimal control and with the implementation of optimal control to validate our theoretical analysis. The paper ends with a through discussion from the view point of novel coronavirus pandemic.

The mathematical model

We investigate an extended SIR (susceptible–infected–recovered)-type epidemiological model that incorporates the interventions such as quarantine and hospitalization. The epidemiological model stratify the dynamics of six individuals, namely susceptible (S(t)), asymptomatic infected (A(t)), clinically ill and symptomatic infected (I(t)), quarantined (Q(t)), isolation or hospitalization (J(t)) and the recovered (R(t)). Thus, the total population N(t) is given by . In our model, quarantine represents the separation of SARS-CoV-2-infected individuals from the susceptible population before the progression of the clinical symptoms, while hospitalization delineates the separation of COVID-19 infected populations with such clinical symptoms. A schematic diagram of the proposed novel coronavirus model based on the present situation of India is shown in Fig. 1. The transmission dynamics of SARS-CoV-2 is modeled by the following system of coupled nonlinear ordinary differential equations:subject to nonnegative initial conditions:In our model, we include the demographic effects by assuming a proportional natural death rate in each of the six individuals. Furthermore, we include the net inflow of susceptible population at the rate per unit time. This parameter indicates the new births and immigration of the individuals. The susceptible individuals become asymptomatic infected individuals, due to effective contact with exposed class at a rate . The transmission coefficients for the asymptomatic infected, symptomatic infected and hospitalized or isolated infected populations are , , and respectively. Among the quarantined individuals, portions move to the isolation or hospitalization due to seriousness of the diseases, and the part become susceptible to the novel coronavirus diseases after the quarantined period. Here, represents the rate at which asymptomatic infected individuals become clinically ill or infectious, is the rate at which asymptomatic infected individuals become symptomatic infected individuals by a proportion of , and another portion move to the quarantined individuals. The term delineates the natural decay rate of all the individuals, whereas , and denotes the symptomatically infected individuals move to the isolated class, and quarantined individuals also move to the isolated class, respectively. The terms , , and describe the rates at which asymptomatic infected individuals, clinically ill or symptomatic infected populations, and isolated infected individuals move to the recovered stages, respectively. Our epidemiological model is a system of coupled nonlinear ordinary differential equations. In exercise, some of the parameters may alter in time as control measures are executed or altered. The analytical calculations are studied for the model with constant parameters, and all the parameters are assumed to be positive.

Fig. 1

Schematic diagram of the novel coronavirus model (1). The epidemiological model includes six sub-populations: susceptible S, asymptomatic infected A, symptomatic infected I, quarantined Q, isolated or hospitalized J and recovered R individuals

Positivity and boundedness

The system of equations (1) is related to novel coronavirus, thus we must have to check whether the state variables are nonnegative for all time t with given initial values (2). Next, we shall investigate the solutions of the system (1) are uniformly bounded. In this purpose, we study the following theorem.

Theorem 1

All the solutions of the system (1) with initial conditions (2) remain nonnegative and are uniformly bounded in the positively invariant region for all time .

Proof

Let us consider be a solution of (1) for , where . Now, we assume that there exists some time such that and for the first time , . But this assumption ends with a contradiction, since for the model system (1). Therefore, , for all .

From the first equation of (1), we obtainwhere . After integration, we getTherefore, for all , we obtain that . The second equation of (1) can be written asthis leads toThird equation of (1) giveswhich implies thatIn the similar manner, we can write the fourth and fifth equations of (1) can be expressed asThen after integration, we haveIn the same way, the last equation of (1) can be expressed aswhich givesTherefore, the solutions of (1) satisfying the initial conditions for all , remain nonnegative in the region .

Now, we prove the boundedness of the solutions (S, A, I, Q, J, R) of (1). The positivity of the solutions implies thatFrom the above equation, we can write thatand consequently,Here, we consider the total population N as . Taking derivative of (1), we havethis leads toThus, we obtain and consequently

Therefore, all the solution trajectories satisfying the initial conditions are uniformly bounded in and in the regionTherefore, all the solutions of (1) initiating in will enter into the region in a finite time frame.

Basic reproduction number

The system of equations (1) has a unique disease-free equilibrium (DFE) point where and , that is, is defined by where and . Now, we determine the basic reproduction number of our compartmental model (1) using the next-generation matrix and the approach developed by Diekmann et al. [47], van den Driessche and Watmough [48], Khajanchi et al. [49]. To do this, we consider the matrix and a non-singular M matrix , represents the creation of new-infection and transition part, respectively, for the equations (1). We linearized infectious subsystem at the DFE as , where are described by are given byThe Jacobian matrix of the model (1) can be calculated at the disease-free steady state , and we getThus, we can derive the next-generation matrix , and the basic reproduction number is the spectral radius of the next-generation matrix, and we haveThe basic reproduction number is an epidemiologic measurement, and it helps us to convey predictions about transmissibility of emergent infectious disease like COVID-19.

Sensitivity analysis of

In this subsection, we study the normalized forward sensitivity analysis of the system parameters associated with the basic reproduction number to analyze their robustness corresponding to the SAIQJR model system prediction as the initial transmission of any epidemic is associated with the reproduction number . Forward sensitivity analysis helps to recognize the most effective parameters with respect to . Since a little perturbation to the value of the most sensitive parameters can impose a great effect in any epidemiological system. Thus, to convey predictions about the prevalence of the diseases, persistence and reduction in the transmission of COVID-19 infection, we performed the normalized forward sensitivity analysis.

Definition 1

The normalized forward sensitivity index of any model parameter, m related with differentially is defined as:The normalized forward sensitivity index may be dependent on the model parameters or may be independent (constant value).

Now, we calculate the normalized forward sensitivity indices of the system parameters associated with , and the indices are listed in Table 1.

Table 1

Sensitivity indices of the SAIQJR model parameters associated with

Parameter	Λ_s	β_s	α_a	α_j	α_q	ρ_i	β_i
Sensitivity
Index	1.00000	1.00000	0.67680	0.32315	0.31271	0.00482	0.00029
Parameter	α_i	γ_i	β_q	σ_a	γ_a	γ_j	δ
Sensitivity
Index	0.00005	-0.00003	-0.01188	-0.22556	- 0.25519	- 0.32154	- 0.92141

From Table 1, it can be noticed that some of the sensitivity indices are positive and the rest are of negative. The positive or negative sign represents that the perturbation to these parameters can increase or decrease the value of , respectively. As for example, it can be observed that the normalized forward sensitivity index of is . This sensitivity index indicates that if the value of would be increased by 10%, as a result the value of would be increased by 6.768%. Again, it is obtained that the sensitivity index of is given by . The index is indicating the fact that 10% increase to the value of will draw 9.2141% decrease to the value of . From Table 1, it is found that the normalized forward sensitivity indices of the parameters , , , , , , and have positive influence on (that is, increase the value of ) and the rest of the indices (the parameters , , , , and ) create negative influence on (that is, decrease the value of ). In this scenario, we notice that the most effective parameters of the system (1) are the disease transmission rate () and the recruitment rate of susceptible individuals () which are inducing positive impact on . Besides,the natural mortality rate of all the individuals () is also very influencing parameter imposing negative impact on . We have to pay more attention to these three highly sensitive parameters to predict the COVID-19 disease transmission. We have shown this behavior of the parameters with reference to in Fig. 5.

Fig. 5

Tornado plot for normalized forward sensitivity indices of the model parameters related with basic reproduction number are listed in Table 1 indicating that the most sensitive parameters are , and

Local stability in terms of the basic reproduction number

In this subsection, we studied the local asymptotic stability of the system (1) with respect to the state variables are A, I, Q, J, that is, the variational matrix related to these variables of the system (1) is as follows:From the above Jacobian matrix , we assume that be the eigenvalue of the matrix. Thus, the eigenvalues of the Jacobian matrix are the roots of the following characteristic equationwherewith . Due to the well-known Lienard–Chipard test [50, 51], all the roots are negative or have negative real part, provided that the following conditions hold: To verify these expressions of the Lienard–Chipard test, we express the coefficients , , and of the above characteristic equation (5) in terms of and is written as:Furthermore, we also calculate the expression in terms of as follows:From the above expression, we can conclude that if , then the Lienard–Chipard test holds successfully and thus the disease-free steady state 0, 0, 0, 0) is stable. For , we obtain due to Descartes’ rule of sign, we decide that at least one of eigenvalues is nonnegative. Hence, the model equations (1) is unstable. The result can be summarized in the following theorem.

, ;

.

Theorem 2

The disease-free steady state of (1) is locally asymptotically stable for and unstable for .

Existence of endemic steady state

In this subsection, we determine the endemic steady state with their feasibility conditions, whereThe system (1) has a unique biologically feasible interior equilibrium point , if the following conditions hold:

(i) ,   (ii) ,    (iii) .

Global stability analysis of the disease-free steady state

Theorem 3

The disease-free steady state of (1) is globally asymptotic stable if and unstable if .

From the system of equations (1), it can be observed that the disease-free compartments of (1) are S, R and the infected compartments are A, I, Q and J. We can arrange the system equations (1) as follows:whereFollowing the method developed by Castillo-Chavez [52], we shall try to investigate the global stability of the disease-free steady state and to compute this, the system equations (1) must have to satisfy the following two conditions: Therefore, the system (1) can be rewritten asTherefore, it can be observed that for , and also we can find that is globally asymptotic stable steady state of the limiting system, . Accordingly, the two conditions are satisfied and in this way the disease-free steady state is globally asymptotic stable, when .

, is globally asymptotically stable.

, where is the Metzler Matrix and , the biologically feasible and attracting region where all the solutions trajectories will enter into the interior of this region and will never leave it.

Local stability of endemic steady state

Theorem 4

The unique interior steady state is locally asymptotically stable if .

The standard linearization technique of the SAIQJR model around the unique interior equilibrium point is very laborious mathematically. Thus, we use the center manifold theory [52] to investigate the local asymptotic stability of the interior infection steady state . To apply this method, we use the following simplification and change of variables are made at first. Let , , , , , , using this following vector notation, the system (1) can be written in this following form as:The variational matrix of (1) around the disease-free fixed point is given byHere, we used , , , and that are already mentioned as earlier before and . Considering as the bifurcation parameter, we calculate by using . Now, we compute as and is given byIt is worthy in mentioning that one of the eigenvalues of the above Jacobian is 0 (simple eigenvalue) and the other four eigenvalues are negative or have negative real parts. Therefore, apply the center manifold theorem [52] near and investigate the dynamics of the system (1).

We obtain the right eigenvector and the left eigenvector corresponding to the Jacobian . The components of the two eigenvectors w and v are given byandThus, we can compute the constants a and b are as follows:Considering the nonzero second order partial derivatives measured by replacing the value at the disease-free equilibrium point at , we haveTherefore, and at bifurcation parameter a transcritical bifurcation happens at . Therefore, the unique interior fixed point is locally asymptotically stable if when is near to 1.

The optimal control problem

In this section, we employ the theory of an optimal control, which is a powerful tool for theoretical analysis corresponding to the infectious diseases. By using this mathematical tool, we have produced conditions which will point out the appropriate steps to mitigate the span of novel coronavirus infection in the human population. Our main aim is to minimize the infected individuals and the hospitalized individuals by using 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 SARS-CoV-2 virus and to minimize the disease burden and their implementations. The definition of the parameters and assumptions lead to the model (1), which implies a coupled system of nonlinear differential equations with six state variables, that is, S(t), A(t), I(t), Q(t), J(t) and R(t). In this control problem, we have introduced three control variables that externally control the number of isolated or hospitalized cases and clinically ill or symptomatic infected cases over a specified time frame.

Enhancing the response of susceptible individuals by modifying the isolation rate

Susceptible population begins to become aware of the diseases and their control strategies when they are given sufficient statistics and as a result the behavior of the population is changed. Due to the systematic news, the population has begun to wear masks, maintaining social distances with proper hygiene, remaining at quarantine and even following home isolation strategy. Thus, we have modified the isolation rate by introducing the control factor that quantifies the external attempt keep in preventing or controlling the of interplays between isolated (J) and susceptible (S) individuals. Therefore, in the term , the parameter represents the maximum transmission rate per isolated and per susceptible populations with representing the depletion in created by the outside isolation or hospitalized control; quantifies the usefulness of the control (where ). The most expected framework is , that is, the condition when interplays between isolated and hospitalized populations are almost “perfectly” avoided (prevented via successful isolation) decreasing the total transmission rate successfully to 0. Here, 0 describes no response and 1 describes fully response due to the given treatment. Therefore, this response intensity is explicitly associated with population’s behavioral response. We have taken into account this response intensity as a control variable. To increase the response among population via isolation or hospitalization, so that they alter their behavior and the cost will be associated as a nonlinear function of . We wish to obtain the optimal response for susceptible population via isolation or hospitalization.

Better treatment strategy for infected individuals

To provide the better treatment policy to the symptomatic infected individuals not only decreases the disease prevalence but also influences its development. At the present scenario in novel coronavirus, we assume that the treatment is accessible and is given to the clinically ill population. As the accessibility of resources associated to medical treatments cannot be unbounded, thus in the model system (1), we introduce the treatment in the form . Here, denotes the treatment rate with intensity rate . There are different costs associated such as medicines, kits, diagnosis, health protocols and other associated costs when medication is given. Therefore, we assume that the treatment intensity as a control variable lies between 0 and 1. The controlling attempt changes the fraction of asymptotic infected cases obtaining antivirals per unit of time.

Provided treatment to isolated individuals

In order to give the better treatment strategy against novel coronavirus to the isolated or hospitalized individuals at the very beginning of infection can lower the disease mortality. It also influences the development of the diseases. It can also 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 equipments are restricted. To keep this mind, we modified the model system (1) by introducing the treatment in the form . Here, denotes the treatment rate with intensity rate . There are different cost related to the medicines, health problem, isolation, vaccination, medical kits, etc., at the time of medication period that need to be taken into account. Thus, we consider as a control variable that lying between 0 and 1, where 0 represents no response and 1 represents full response after the treatment. The controlling effort changes the fraction of isolated population () obtaining the antiviral treatment per unit of time. Let us consider that the controls attain maximum per capita effectiveness vaccination rates when the controls become 1.

The main aim of this section is to obtain the optimal treatment strategy and optimal controls through clinically ill infected population and isolated individuals with minimum cost with the aid of vaccination. Thus, the admissible set for three control variables , and is defined as follows:The following optimal control problem is investigated with the above discussed control strategies along with the objective functional to be minimized:subject to minimize the objective functionalwith the initial conditions are stated in Eq. (2). Here, the objective functional represents the total incurred cost, that is, the sum of of the costs stated in the integrand:represents the current value of the cost at any time t. All the parameters , , , , and are nonnegative. They are used as a weight constants and balancing the units of the integrand. Also they represents the measure of the relative costs of the control functions over [0, T]. The control functions and exists in the admissible control set that mainly minimized the objective functional .

Existence of optimal control

In this subsection, the existence of an optimal control of (9) is analyzed. Then, we employ the Pontryagin’s maximum principle [53] to characterize the optimal control. The Pontryagin’s maximum principle is being utilized to get necessary conditions for the three controls. To prove the existence of an optimal control, we utilize the following theorem which assured that the existence of the controls and that minimizes the objective functional

Theorem 5

There exist optimal controls in that minimizes the objective functional corresponding to the control system (9)–(10).

In Sect. 2, we have already proved that all the solutions of (1) are nonnegative and uniformly bounded for the given initial conditions (2). Nonnegativity of the objective functional is obvious and the nonnegativity implies the boundedness of the objective functional. Hence, minimizing the sequence of controls exists in such a way thatTherefore, the controls belonging to are uniformly bounded in and consequently they are uniform bounded in . Since the space is reflexive [54], then there exists a subsequence in such a way thatThen, it is obvious that the state sequence is uniformly bounded corresponding to . It can be seen that the uniformly boundedness of the right-hand side of the system (1) implies the uniformly boundedness of the derivatives for and equi-continuity of the corresponding state sequence . Due to the Arzel–Ascoli theorem, there exist such that on a subsequence,Next, we choose a proper subsequence passing through the limit of the system (1) and corresponding to the controls , we can achieve the state solution . Then, the lower semi-continuity of the —norm with respect to weak convergence implies thatTherefore, is an optimal control.

Characterization of optimal control

Since the existence of the controls are established for minimizing the objective functional (10) subject to the state equations (9), then we employed Pontryagin’s maximum principle [53] to obtain necessary conditions for the optimal control problem [55, 56].

Theorem 6

There exists optimal controls and the corresponding state solutions of the system (9) that minimizes the objective functional , then there exist adjoint variables , , , , , and satisfying the following canonical equations:with the transversality conditions:Moreover, the corresponding optimal controls are represented by

We investigate necessary conditions for the control functions with the aid of the Pontryagin’s maximum principle for the model system (9). In order to do this, we construct the Hamiltonian which is defined for all as follows:where , , , , , and are the co-states or adjoint variables. With the aid of Pontryagin’s maximum principle, we obtain the minimized Hamiltonian which minimizes the objective functional . Pontryagin’s maximum principle plays an important role in adjoining the objective functional with the state variables including co-states or adjoint variables.

Due to Pontryagin’s maximum principle, there exist co-states , , , , and that satisfies the following canonical expressions:with the transversality conditions: and . Now, we get the characterization of optimal controls are given byFrom , , and , we getBy considering the upper and lower bounds for , , and , we have following the characterization of the optimal controls:It can also be noted thatwhich assured that our optimal controls minimize the Hamiltonian.

Now, we derive the corresponding optimality system utilizing the optimal controls that has already been characterized in the above. The optimality system that minimized Hamiltonian at , , is given as follows:with initial values: , , , , , and the corresponding co-states arewith the transversality conditions are , , , , , , and the three controls , and are the same as (12)–(14). Hence, the proof of the theorem.

Model calibration and data source

In this section, we validate our theoretical analysis based on the numerical solutions with real data from India. In order to do this, we estimate the most sensitive parameters for our SAIQJR model.

Sensitivity analysis

We have calibrated our novel coronavirus model (1) by using the daily confirmed new SARS-CoV-2 viruses and cumulative confirmed SARS-CoV-2 viruses for India. We have taken data for novel coronavirus diseases from World Health Organization (WHO) registered cases for the time period February 01, 2021, to May 21, 2021 [58]. Our proposed model (1) has a system of six coupled ordinary differential equations with 14 parameters, among them six key parameters have been estimated by using semi-relative sensitivity analysis [59], namely , , , , and for the novel coronavirus model (1). These six parameters value and the initial individuals size plays a critical role in the model analysis.

Sensitivity analysis depicts the significance of parameter alters in the model output. Due to the lines of Martin Fink [60], we draw the sensitivity graph by utilizing the code myAD, established by Fink [60] (see Fig. 2). To quantify the sensitivities of the parameters from Fig. 2, we calculate the sensitivity coefficient by calculating -norm that can be defined as follows:From the computation of -norm, it is clear that the most effective parameters are , , , , and in descending order (see Fig. 3) by comparing and ranking the sensitivity function that represents the other model parameters are less sensitive with respect to the clinically ill or symptomatic infected population, thus they will affect the dynamics of novel coronavirus a little.

Fig. 2

Semi-relative sensitivities of the system parameters by using automatic differentiation. From the figures, it can be observed that the parameters , , , , , , are most effective with respect to the symptomatic infected individuals. The biologically realistic time window is [0, 400]

Fig. 3

Quantification of semi-relative sensitivity analysis by computing sensitivity coefficient using norm

A parameter is known as practically identifiable if an individual estimate can be determined from various initial conditions utilizing the available data. We consider Fisher’s information matrix , where S represents the normalized sensitivity matrix calculated by automatic differentiation [60]. It can be noticed that m parameters are said to be locally identifiable iff the column rank of the matrix S = m, or identically . After that, we execute QR factorization method with column pivoting in MATLAB routine qr, . This procedure find out matrix P in such a way that (QR is the factorization of FP). This indices in first k columns of P recognize k parameters are most estimable. In our study, , , , , and are the model parameters that are most estimable from the novel coronavirus data. We have estimated the most sensitive parameters by using least square method, after identifying them. The values of these parameters are provided in Table 2. We have estimated the data for 110 days (from February 01, 2021 to May 21, 2021). The numerical procedures and data used for parameter estimation are given in the supplementary material. Model fitting is shown in Fig. 4.

Table 2

Description of system parameter values used for novel coronavirus model system (1)

Parameter	Description	Values (Unit)	Source
Λ_s	Inflow rate of susceptible individual	1.6 × 10^6 personday^-1	[29]
β_s	Disease transmission rate	10^-7  day^-1	Estimated
α_a	Modification factor for asymptomatic population	0.001	Estimated
α_i	Modification factor for symptomatic population	0.0001	[57]
α_j	Adjustment factor for isolated individuals	0.036	[57]
β_q	The rate that quarantined becomes susceptible	0.0135  day^-1	Fit to data
δ	Natural mortality rate of all the populations	0.005  day^-1	Estimated
γ_a	Asymptomatic individuals recovered from the disease	0.0066  day^-1	Estimated
σ_a	Progression rate of infectious	0.014  day^-1	Estimated
ρ_i	Fraction of asymptotic individuals showing clinical symptoms	0.031  day^-1	Estimated
γ_i	Infected population progress to the recovered class	0.002  day^-1	[26]
β_i	Infected population progress to the isolated class	0.653  day^-1	[25]
α_q	Rate at which quarantined becomes isolated	0.337  day^-1	Assumed
γ_j	Isolated population progress to the recovered class	0.997  day^-1	Fit to data

Fig. 4

Parameter estimation for the novel coronavirus model system (1) based on data from India. The SAIQJR model fitted with the real data for the daily new confirmed COVID-19 cases and the cumulative confirmed COVID-19 cases of India for the time period February 01, 2021 to May 21, 2021 (110 days). Real data are shown in blue circles, whereas the red curve represents the best fitting curve of the system (1). Parameters are used for numerical solution and specified in Table 2. The initial size of the individuals used for numerical simulation are given in the text

Model calibration

We calibrated the proposed SAIQJR model by fitting the daily COVID-19 infections and cumulative COVID-19 cases in India from February 01, 2021, to May 21, 2021 (110 days). We used the nonlinear least squares method to minimize the root-mean-squared error (RMSE) between model-simulated and reported data. Due to the varied orders of magnitude of data, we standardized these data for the same orders of magnitude by dividing by the max value of the data sets, which could avoid the fitting process automatically skewing larger orders of magnitude data and better match each type of data. We estimated six most effective parameters (identified through sensitivity analysis) by data fitting and obtained the other model parameters from the published literatures (see Table 2).

To measure the accuracy of our model fitting, we have computed mean absolute error and root-mean-square error . The and are defined as follows:where C(i) denotes the observed value, S(i) denotes the model simulation and n is the sample size of the observed data. The value of for India is given in Fig. 4. Our model performs well in the case study of India.

Model simulation without control strategy

We performed the numerical simulations in twofolds, namely without administration of optimal control and with the administration of optimal control theory. For the validation of our theoretical analysis, we conduct some numerical simulations for the model parameters that are specified in Table 2. Due to the present situation of the COVID-19 pandemic in India, the newly infected population per unit time is represented by . The total population in India in August 2021 is approximately 138 million, and according to the data from World Health Organization (WHO) [58], the daily confirmed number of infected cases on February 01, 2021, is 11,427. Thus, we have used the following set of the initial size of the populations:

Figures (Fig. 6a–d) represent the contour plots that demonstrate the dependence of the basic reproduction number with respect to (disease transmission rate) and (progression rate of infected individuals); (adjustment factor for clinically ill population) and ; (fraction of isolated or hospitalized individuals) and ; and (infected individuals progress to isolated class) for India. From Fig. 6a, we noticed that for an increasing value of disease transmission rate and , the value of the basic reproduction number increases. Similar phenomena has been noticed from Fig. 6b–d. Thus, from Fig. 6a, we can conclude that in order to mitigate the transmission of SARS-CoV-2, the basic reproduction number need to be maintained and that can be possible by maintaining the social distancing as well as contact tracing.

Fig. 6

Contour plots of the basic reproduction number for the case of India. a Contour plot of as a function of and ; b Contour plot of as a function of and ; c Contour plot of as a function of and ; d Contour plot of as a function of and . All parameter values other than (a) and ; (b) and ; (c) and ; (d) and are specified in Table 2

As we are familiar with that the numerical value of describes the behavior of the infectious diseases. From Table 1 and Figs. 5, 6a–d, we can see that parameters with non-positive indices can be improved, and then the magnitude of would be diminished substantially. The model parameters, namely , and , has nonnegative indices, and thus mitigation of these three model parameters can control the outbreak of the COVID-19 diseases. Moreover, from the sensitivity indices, we can observe that three parameters associated with susceptible population, namely , and , are the most effective parameters, and we can restrict positively these three parameters. Therefore, we may deduce that the most favorable method to maintain the COVID-19 diseases is to decrease the value of the parameters.

To better understand the threshold parameter , we have drawn a surface plot for with respect to (disease transmission rate) and (rate by which asymptomatic populations become recovered), and the other parameters are defined in Table 2. The expression for basic reproduction number is defined in (4), and we note that and . If we increase , the basic reproduction number will increase and cross the value , thus leading to the substantial outbreak of the SARS-CoV-2 virus. Moreover, when the recovery rate of asymptomatic individuals increases, the basic reproduction number will decrease, and if becomes lower than 1, the infected populations will no longer persist. The surface plot shown in Fig. 7 represents the behavior of in the range of and . Therefore, our model illustration proposed that the recovery rate of asymptomatic individuals aids to eliminate the novel coronavirus by decreasing the basic reproduction number . Also, our simulations suggest that the reproduction number can be maintained by decreasing the value of the disease transmission rate and that can be obtained by maintaining social distancing and contact tracing.

Fig. 7

Surface plot of as a function of (disease transmission rate) and (the rate at which asymptomatic individuals become recovered). All parameter values other than and are given in Table 2

Model simulation with control strategy

Now, we introduce additional recovery through the medication by using the theory of optimal control problem. In the model system (9), we have introduced three different control policies to mitigate the transmission of COVID-19 burden and to minimize the total cost by obtaining the optimal control paths. We performed numerical simulations to investigate these outcomes and also to observe the influence of different parameters and intervention strategies on the transmission dynamics of COVID-19. To investigate the numerical simulation of the objective functional (10) subject to the coupled system of ordinary differential equations (9), we utilize fourth-order Runge–Kutta method with the aid of MATLAB. We solve the co-states (11) by using backward Runge–Kutta method by employing the solutions of (9) with the transversality conditions defined in Theorem 6. After that, these controls are upgraded for the optimality utilizing the Hamiltonian of the optimal system. To do this, we used the steepest descent method of optimization that has been studied in [55, 56]. The time frame of investigation and implementation of the optimal control strategies is considered to be 400 days. To solve the controlled system and the corresponding adjoint system, we use the parameters set defined in Table 2.

To analyze the impact of optimal treatment policies for the transmission dynamics of novel coronavirus, we employ different control policies for the disease compartments, namely asymptomatic population, clinically ill population, quarantine populations and the isolated individuals. In general, high basic reproduction number () create epidemics with high epidemic peak. Figure 8a–d represents the impact of each treatment strategies on the disease state variables in the presence of optimal control as well as without optimal controls when = 1.8469. These plots demonstrates the asymptomatic cases, daily number of infected cases, quarantine individuals and the isolated or hospitalized cases under no control , with the combination of three different controls and for single controls where . The solid red curves for epidemic (without administration of intervention strategies) are exhibited to highlight the difference from those created via implementation of optimal treatment polices. The solid blue epidemic curves (with the implementation of three controls) are shown the significant reduction in the asymptomatic cases, clinically ill cases, quarantine cases and isolated cases when compared with the no intervention scenario. The solid magenta curves, solid green curves and solid black curves represents the single controls, namely and , respectively, shows the reduction in the asymptomatic cases, clinically ill cases, quarantine cases and isolated cases when compared with the no intervention strategies as well the implementation of three controls. Implementation of single control is better in comparison with the without control policy. Also, in the cases for asymptomatic cases (only ), clinically ill cases (only ), quarantine cases (only ) and isolated or hospitalized cases (only ) show the significant reduction in the individuals when compared with the no intervention scenario and compared with the implementation of three controls. Thus, from Fig. 8a–d, we can conclude that the implementation of intervention strategy has an impact in controlling the transmission dynamics of novel coronavirus epidemic.

Fig. 8

Comparisons of the corresponding to the (a) asymptomatic infected cases, (b) clinically ill or infected individuals, (c) quarantine individuals, (d) hospitalized or isolated individuals without intervention strategies, with the implementation of intervention strategies (combination of three controls and single controls). Optimal treatment strategy (solid blue line) demonstrates substantial reduction in all the state solutions when compared with the no controls and the application of single controls. The parameter values are given in Table 2 when

Now, we employ different treatment strategies for the disease compartments, namely asymptomatic individuals, symptomatic or clinically ill population, quarantined cases and the isolated individuals. The time series Fig. 9a–d 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 for . These plots demonstrate the asymptomatic cases, daily number of infected cases, quarantine individuals and the isolated or hospitalized cases under no control , with the combination of three different controls and the combination of two controls, namely and ; and ; and . The solid red curves for epidemic (without implementation of intervention strategies) are shown to highlight the difference from those created via administration of optimal treatment polices. The solid blue epidemic curves (with the administration of the combined of three controls) are exhibited the significant reduction in the asymptomatic cases, clinically ill cases, quarantine cases and isolated cases when compared with the no intervention strategies. The solid magenta curves, solid green curves and solid black curves denotes the combination of two controls, namely and ; and ; and , respectively, shows the reduction in the asymptomatic cases, clinically ill cases, quarantine cases and isolated cases when compared with the no intervention strategies as well the implementation of three controls. Implementation of the combination of two controls is better when compared with the no control policy. Also, in the cases for asymptomatic cases (when and ), clinically ill cases (when and ), quarantine cases (when and ) and isolated or hospitalized cases (when and ) shows the significant reduction in the individuals when compared with the combination of three controls. Thus, from the numerical simulations (Fig. 9a–d), we can infer that the implementation of intervention strategy has an impact in controlling the transmission dynamics of novel coronavirus epidemic. Also, we can conclude that the combination of two controls are more effective when compared with the implementation of the single control. It can also be noted that combination of three controls are more effective when compared with the combination two controls as well as single controls.

Fig. 9

Comparisons of the corresponding to the (a) asymptomatic infected cases, (b) clinically ill or infected individuals, (c) quarantine individuals, (d) hospitalized or isolated individuals without intervention strategies, with the implementation of intervention strategies (combination of three controls and combination of two controls). Optimal treatment strategy (solid blue line) demonstrates substantial reduction in the cases in all the state solutions when compared with the no controls. The parameter values are given in Table 2 when

The control functions ) and ) represents the antiviral treatment effectiveness aimed at clinically ill or infected cases and hospitalized or isolated populations, respectively. Another control function ) maintains the mitigation of the creation of secondary infected cases resulting from the isolated individuals. The numerical results for the controls are shown in Fig. 10a–d. The numerical solutions of the implementation of three control effectiveness are simultaneously shown in Fig. 10a. Corresponding to the implementation of the combinations of two controls, namely controls ) and ); ) and ); ) and ) are shown in Fig. 10b–d, respectively. The plot of the three graphs (see Fig. 11) represents the optimal control functions where single control strategies, namely , and have been implemented. The optimal control strategies based on single controls need the administration of in-depth effectiveness at the starting of an outbreak (Fig. 11a–c) of the novel coronavirus followed by quick reductions (most likely the outcome of high reduction in the levels of susceptible individuals).

Fig. 10

Optimal control functions as a function of time corresponding to the (a) three controls (, and ); (b) two controls ( and ); (c) two controls ( and ); (d) two controls ( and ). Parameters are given in Table 2

Fig. 11

Optimal control functions as a function of time corresponding to the single control (a) , (b) and (c) . Parameters are given in Table 2

Discussion

Mitigating the impact and spread of the ongoing novel coronavirus pandemic was on the minds of every newscaster, Govt. official, public health departments, aspiring politicians and literally millions of populations throughout the world within a couple of months ago. According to the report by WHO, about 22 crore populations are infected with the novel coronavirus throughout the world [58]. Though the Govt. of different countries and clinicians of each and every country are trying to give safeguards to the individuals, still the rate of infection is high due to absence of proper vaccination for the COVID-19 virus. After the Second World War, the world has never gone through like the current situation that arrived due to novel coronavirus pandemic. Beginning from the Wuhan city of China in December 2019, the SARS-CoV-2 has outspread throughout the world. Due to the absence of proper antiviral or effective treatment of novel coronavirus, one can rely on the qualitative control of novel coronavirus disease in spite of the elimination. In the absence of proper vaccination and with an incomplete understanding of epidemiological cycle, mathematical models may help in understanding of both novel coronavirus transmission and its control.

To keep this pandemic scenario in our mind, we proposed and studied a compartmental SAIQJR model for novel coronavirus. Theoretical as well as numerical analysis of model system has been investigated using stability theory. Positivity and boundedness of the state equations assured that the proposed model (1) is well-defined. Our model system (1) possess two equilibrium points, namely a unique disease-free fixed point and a unique endemic equilibrium point whose feasibility and stability criteria depends on the numerical value of . Using the theory of center manifold, we showed that the unique endemic steady state is locally asymptotically stable when basic reproduction number exceeds 1. We also performed the global stability criteria for the disease-free fixed point [61]. From the parameters in Table 2, we obtained the basic reproduction number that demonstrates the significant outbreak of the SARS-CoV-2 virus.

We calibrated our proposed model based on the daily new novel coronavirus cases and cumulative number of novel coronavirus cases of India. Our proposed model has 14 parameters, and we estimated only six parameters depends on the sensitivity analysis. As the sensitivity analysis exhibited that only six parameters are most sensitive with reference to the clinically ill or infected individuals. We also performed the normalized forward sensitivity analysis for our proposed model which reveals that the disease transmission rate , recruitment rate of susceptible individuals, modification factor for asymptomatic individuals and the mortality rate of all the individuals have an impact effect on the clinically ill individuals and . Our model illustration shows that to control the disease transmission must have to maintain the social distancing, contact tracing, avoid mass gathering and use proper sanitization.

Moreover, we have modified our SAIQJR model by introducing optimal control problem. Implementation of intervention policies aids to mitigate the COVID-19 disease burden. To investigate the influence of non-pharmaceutical and pharmaceutical intervention policies and to stop the community transmission of the SARS-CoV-2 virus, we implemented three different control measures for . We established the existence of the optimal control functions. Moreover, to mitigate the infected populations and to minimize the cost of the controls, an objective functional has been formed and interpreted with the aid of Pontryagin’s maximum principle.

The behavioral changes in the susceptible individuals alters with time, and to prevent or limit the rate of interaction between susceptible (S) and isolated (J) individuals, we put the external effort 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 controlling effort that changes the fraction of clinically ill cases. Again, we employ the controlling effort that taking into account the cost regarding proper diagnosis, use of test kits and life-saving medicines that given to the hospitalized or isolated individuals. Thus, we considered all these three control strategies in our proposed model.

Finally, through extensive numerical simulations, we conclude certain observations. Numerical simulations shows in Fig. 11a that the intervention strategy presenting behavioral response efforts with high magnitude immediately after introduction but slowly it reduces with time. Another way, the control strategy representing clinically ill or infected cases efforts with its high magnitude for quite a long time and then it reduces (see Fig. 11b). Later on, the controls indicating better treatment strategy for isolated individuals efforts with high magnitude nearly about 3 months following a mitigation at a later stage, though this magnitude is higher than the magnitude of (see Fig. 11c). From the extensive numerical simulations (see Fig. 9a–d), we can infer that the implementation of intervention strategy has an impact in controlling the transmission dynamics of novel coronavirus epidemic. Also, we can conclude that the combination of two controls are more effective when compared with the implementation of the single control. It can also be noted that combination of three controls is more influential when compared with the combination two controls as well as single controls. Thus, implementing all the three control strategies may help to mitigate novel coronavirus disease transmission at this present epidemic scenario. From medical perspective, we presume that our discussion sheds light on the dynamic model of COVID-19 transmission and its control.

It is important to point out that further monitoring is essential since other factors such seasonality of respiratory diseases might play an important role on COVID-19 transmission and control. More classified data are available, specifically the age-structured data, then more complex mathematical models would be helpful for in-depth investigation of the COVID-19 pandemic. We plan to address the new progression of the epidemic data in another piece of work in near future. We also plan to expand our modeling efforts by introducing vaccination and the effect of environment contamination. Most probably, this COVID-19 epidemic will remain in the communities for many years, so the changes we see in public health services and in human habits have to go on over time. Actions as simple as washing hand, space hygiene, social distance, contact tracing and use of face masks in closed spaces should be incorporated into education for health. The containment measures, which should be essential when epidemics arise, must be rigorously investigated and worked with families.

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 516 KB)