A model based study on the dynamics of COVID-19: Prediction and control.

As there is no vaccination and proper medicine for treatment, the recent pandemic caused by COVID-19 has drawn attention to the strategies of quarantine and other governmental measures, like lockdown, media coverage on social isolation, and improvement of public hygiene, etc to control the disease. The mathematical model can help when these intervention measures are the best strategies for disease control as well as how they might affect the disease dynamics. Motivated by this, in this article, we have formulated a mathematical model introducing a quarantine class and governmental intervention measures to mitigate disease transmission. We study a thorough dynamical behavior of the model in terms of the basic reproduction number. Further, we perform the sensitivity analysis of the essential reproduction number and found that reducing the contact of exposed and susceptible humans is the most critical factor in achieving disease control. To lessen the infected individuals as well as to minimize the cost of implementing government control measures, we formulate an optimal control problem, and optimal control is determined. Finally, we forecast a short-term trend of COVID-19 for the three highly affected states, Maharashtra, Delhi, and Tamil Nadu, in India, and it suggests that the first two states need further monitoring of control measures to reduce the contact of exposed and susceptible humans.

Introduction

The recent pandemic, which is commonly known as COVID-19, is an infectious disease caused by the virus, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). In late December 2019, the disease COVID-19 was first identified in Wuhan, the capital Hubei province, China, and causing the first pandemic of this century. Since the start of the year 2020, the infectious disease COVID-19 has started to spread globally and resulting in almost three million positive cases till today (till 27th April, 2020 there with about one million closed cases having a 20% death rate). As we understand, the disease COVID-19 has become a pandemic due to several reasons. Some of them are (i) unavailability of both suitable vaccine and exact medicine, (ii) high disease transmission rate, and (iii) the precise nature of the virus SARS-CoV-2 being still unknown. As many as 210 countries and territories around the globe have been suffering from COVID-19 today [see [16]].

The primary step taken by most of the governmental agencies to control COVID-19 is the implementation of lockdown to maintain social distance. This procedure is an excellent measure to control the spreading of the disease. Still, from an economic point of view, the complete lockdown may be the cause of a significant financial crisis for the near future. In particular, lockdown in high dense countries may reduce the disease transmission rate, although complete control may not be achievable. Hence to alive the economic status of a country, a full lockdown for an indefinite period is not desirable at all in any circumstances. Therefore there should be a suitable balance between the two different characteristics of governmental policies complete lockdown and healthy free conditions.

Therefore, in the present scenario, a qualitative analysis of COVID-19 is more significant compare to quantitative analysis. Hence a suitable mathematical model would not only able to represent the whole disease system but also the study of the model would undoubtedly derive the precise nature of the disease. It may forecast the behavioral aspect of the disease shortly. Although the primitive mathematical models on theoretical epidemiology (see Bernoulli [1], Hamer [13], Ross [40]. Kermack and McKendric [24], etc) look quite simple, still, from a new perspective, those works are the major works on mathematical epidemiology. With the advancement of computational tools and software, we can develop a complex mathematical model and analyze it thoroughly in a scientific manner. In the history of literature, we see that many model-based studies have successfully achieved the global dynamics of the corresponding infectious disease (see the references Keeling and Rohani [23], Wang and Zhao [41], Buonomo et al. [3], Zhou et al. [43], Jana et al. [18], [20], Li et al. [29], Zegarra and Hernandez [42] etc).

According to the information received, it may take around one week to two weeks for the exposure of symptoms of COVID-19 of an infected person, although during this period, that person able to infect other susceptible persons. However, there may be some infected persons whose infection is so mild that the person would recover due to innate immunity even before the hospitalization. Thus in this article, by the term ’infected’ person, we will mean those persons who are hospitalized. Further, we assume that the medical personals assisting COVID-19 positive hospitalized individuals have taken necessary protective items. Thus to keep simplicity, we believe that only exposed persons and asymptomatic infected persons can spread the disease. Still, the infected persons who are hospitalized in due course are not spreading the virus. Although the effect of COVID-19 started not more than six months ago, some significant works have already been done or are undergoing to expose COVID-19 with the help of suitable mathematical models. Kucharski et al. [26] have proposed a mathematical model-based analysis of COVID-19, where the authors have considered all the positive cases of Wuhan, China, till 5th March 2020. In another recent work, Ndairou et al. [35] have presented an autonomous system of mathematical model to study the spread of COVID-19 in the Wuhan city, China. Using mathematical modeling, Prem et al. [38] have analyzed the controlling status of COVID-19 of Wuhan city. Hellewell et al. [15] have described the effective procedure of COVID-19 disease using isolation. In another work, Mizumoto and Chowell [33] have studied the transmission dynamics of COVID-19 in the international conveyance, Diamond Princess Cruises Ship. The basic reproduction number is an important parameter to analyze the nature of an infectious disease. In their work, Liu et al. [30] have investigated the possible numerical value of the basic reproduction number of COVID-19. In a recent article, Fanelli and Piazza [7] have analyzed and predicted the nature of COVID-19 in three most affected countries till March 2020 with the help of mathematical modelling. Ribeiro et al. [39] have used some stochastic based regression models to forecast the phenomena in as many as ten most affected states of Brazil. In other recent work, Chakraborty and Ghosh [4] have considered a hybrid ARIMA-WBF model to forecast various COVID-19 affected countries throughout the globe.

Since the nature and destruction of COVID-19 depends on various parameters (namely personal immunity, history of visiting into a COVID-19 pandemic country, maintaining the required hygiene, etc.) of the affected system, using a single model we cannot describe the whole disease system throughout the globe. Addition of detail and complexity can make models more accurate but this also complicates their mathematics. However, even this kind of simple model is very helpful when formal vaccination or proper treatment control is not available. Motivated by this, in the present paper, we propose a mathematical model introducing a quarantine class and governmental intervention measures like lockdown, media coverage on social distancing, and improvement of public hygiene, etc to mitigate disease transmission. The main purpose of this work is to explore the role of quarantine and the governmental intervention strategies on COVID-19 control and elimination. We study a thorough dynamical behavior of the model in terms of the basic reproduction number. Further, we perform the sensitivity analysis of the essential reproduction number and found that reducing the contact of exposed and susceptible humans is the most critical factor in achieving disease control. To lessen the infected individuals as well as to minimize the cost of implementing government control measures, we formulate an optimal control problem, and optimal control is determined. Finally, we predict a short-term trend of COVID-19 for the three most affected states, Maharashtra, Delhi, and Tamil Nadu in India, and it suggests that the first two states need further monitoring of control measures to reduce the contact of exposed and susceptible humans.

The rest portion of this paper has been organized in the following way: In Section 2, we describe the theoretical analysis of the article, which includes the basic reproduction number, existence, and asymptotic stability criteria of the two equilibria, namely disease-free and endemic equilibrium point. In Section 3, we formulate an optimal control problem and solve the problem analytically. In Section 4, we provide some numerical examples to illustrate the theoretical counterparts. Next, in Section 5, we consider real-world examples on COVID-19 of three most COVID-19 affected states in India and derive some short term forecasting based on the proposed model. In the final section, we provide a thorough discussion of our article.

Model formulation

In this section, we formulate a mathematical model on COVID-19 based on some realistic assumptions. At any time instant t, the human populations are subdivided into five time-dependent classes, namely Susceptible S(t), Exposed E(t), Hospitalized infected I(t), Quarantine Q(t) and Recovered or Removed R(t). Based on those five state variables, we aim to form an autonomous system using first-order differential equations. Let A be the constant recruitment rate to the susceptible population, and β be the disease transmission rate. However the disease transmission from vulnerable to infected persons (here the class E) depends on several parameters, namely, precautions (use of face mask, social distancing, not rubbing face and nose using hand, etc.) and hygienic environment (use of soap and sanitizer, hand washing, cleaning, etc.) taken by both susceptible as well as infected persons. Since here, we have assumed that the virus of COVID-19 is spreading when a vulnerable person comes into contact with an exposed person; therefore we think that ρ 1(0 < ρ 1 < 1) portion of susceptible human would maintain proper precaution measure and ρ 2(0 < ρ 2 < 1) portion of the exposed class would take proper precaution measure for disease transmission (i.e., use of face mask, social distancing and implementing hygiene). Therefore the disease can only be transmitted to the portion of susceptible individuals due to the contact of portion of exposed individuals with a bi-linear disease transmission rate β. We know that a person is whether infected by the SARS-CoV-2 virus or not can be clinically detected using RT-PCR examination and a person with negative results in the RT-PCR test may still be COVID-19 positive as it may take some days (from 7 to 21 days) to express infection. Therefore, the portion with positive COVID-19 of the class of population E is considered as infected, and they are hospitalized. Let α and b 2 be the portions of the exposed class goes to the infected class and quarantine class, respectively. It should be noted that since it would take quite a long time to get the output of the RT-PCR test, and sometimes it requires more than one RT-PCR analysis for a single person for confirmation of COVID-19. Let among the quarantine classes of populations, cQ portion of communities move to infected level, and the b 1 Q part would become susceptible to the disease after the quarantine period. Let η and σ be respectively recovery rate of the hospitalized infected populations I and exposed class E. Let d be the natural death rate, which is common to all classes of communities and δ be the COVID-19 induced death rate. Also, it is statistically observed a person once recovered from the disease COVID-19 has very little chance to become infected again for the same disease. Hence, we assume that no portion of the recovered population moves to the susceptible class back. Also, in the formulation of a mathematical model, we understand that the people of the hospitalized infected class (i.e., I(t)) would not be spreading the disease or spreading a negligible amount of disease since they are kept completely isolated from the susceptible individuals. However, to control the pandemic COVID-19, suitable governmental measures (like complete or semi lockdown, rationing system, continue media coverage on social isolation and improvement of public hygiene, home deliveries of essential commodities, to make an alternate source of income for job-losers during the lockdown, etc.) have been implemented by the various governmental and non-governmental agencies. Thus this policy may be considered as one of the effective control tools, and mainly the susceptible population of COVID-19 cases would be benefited due to this policy. Mathematically we represent this policy by the parameter M and let this policy be implemented at a rate p. Therefore due to the implementation of this policy, the portion, pSM moves from susceptible to the recovered or removed class (R(t)). It should be noted that usually, the parameter M should be time-dependent as for the optimum result, the policy should be applied according to the necessity of the situation and keeping a proper balance between two different states, namely lockdown and complete free state. However, in the qualitative analysis of this model, we assume the parameter M as time-independent, but in formulating an optimal control problem to keep the infected individuals at the minimum level, M is assumed as a time-dependent control M(t), where indeed 0 ≤ M(t) ≤ 1. Based on the above assumptions, we present a schematic diagram of the proposed COVID-19 system in Fig. 1 . Further, based on the above-stated assumptions and the flow diagram of COVID-19, we formulate an autonomous a dynamical system consisting of five first-order differential equations shown as below:

Fig. 1

Flow diagram of disease transmission.

On the basis of above stated assumptions and the flow diagram of COVID-19, we formulated an autonomous dynamical system consisting of five first order differential equations shown as below:

The initial condition of the system (1) are S(0) ≥ 0, E(0) ≥ 0, Q(0) ≥ 0, I(0) ≥ 0, R(0) ≥ 0.

Analysis of the system for fixed control

In this section we assume the fixed value for control parameter M. Here we mainly study the uniformly boundedness of the solutions and subsequently the basic reproduction number, different equilibria and their stability criteria, sensitivity etc.

Boundedness of the system

Here we examine the boundedness property of the system (1).

The solutions of the system (1) are uniformly bounded.

We assume that

Therefore Integrating the above inequality and by applying the theorem of differential equation due to Birkhoff and Rota [2], we get

Now for t → ∞,

Hence all the solutions of (1) that are initiating in are confined in the regionfor any ϵ > 0 and for t → ∞. Hence the theorem. □

Basic reproduction number

In any epidemic model, the basic reproduction number is the most important epidemiological parameter for determining the nature of a disease. Generally, it is denoted by R 0 and is defined as the “number of secondary infected individuals caused by single infected individuals in the whole time interval” (see, van den Driessche and Watmough [6]). Therefore the dimensionless quantity R 0 refers as the expectation of the spreading disease.

In literature, several techniques are available to evaluate R 0 for an epidemic spread. In our present research article we use the next generation matrix approach [5], [9], [22]. Now the classes which are directly involved for spread of disease is only E, Q, I. Therefore from system (1) we have

The above system can be written as where

It is clear from the system (1) that is a disease free equilibrium. Now the Jacobian matrix of Φ and Ψ at the disease free equilibrium are respectively given by, and

The basic reproduction number (R 0) is the spectral radius of the of the matrix and for the present model it is given by

Equilibria

The system has two possible equilibria. One is disease free equilibrium point, where infection vanishes from the system. It is given by E 0(S 0, 0, 0, 0, R 0), where . The other equilibrium point is E 1(S*, E*, Q*, I*, R*), where infection is always present in the system is called endemic equilibria, where

.

Note It is observed from the expression of the above two equilibrium point is that the disease free equilibrium E 0 is always feasible but the endemic equilibrium E 1 is feasible if R 0 > 1.

Stability analysis

In the present section we investigate the local asymptotic stability criteria of the different equilibria.

The disease free equilibrium E 0 is locally asymptotically stable if R 0 < 1 and it is unstable if R 0 > 1.

The Jacobian matrix at the disease free equilibrium of the system (1) is given by

Now the characteristic equation of the system (1) at its disease free equilibrium is given by

Clearly all the eigen value of the Jacobian matrix are negative if and only if R 0 < 1. Hence the system is locally asymptotically stable if R 0 < 1 and it is unstable if R 0 > 1. Hence the theorem. □

Note Here we see that the disease free equilibrium E 0 losses its stability when the R 0 increases to its value greater than 1. So, we may conclude that at R 0 the system (1) passes through a bifurcation around its disease free equilibrium which are discussed in the next theorem..

The system (1) passes through a transcritical bifurcation around its disease free equilibrium when

From the above analysis, it has been observed that when R 0 < 1 between the two equilibria, only the disease free equilibrium exists and locally asymptotically stable where as R 0 > 1 is the threshold condition for both existence and asymptotic stability criteria of the endemic equilibrium point although the disease free equilibrium reduces to unstable nature at the threshold R 0 > 1. Hence we may conclude that there is change of feasibility as well as stability occurs at Following the articles Gukenheimer and Holmes [12], Kar and Jana [21], Jana et al. [19] etc, it is concluded that the system undergoes a transcritical bifurcation at and in Fig. 4, we represent graphically the phenomenon of transcritical bifurcation. □

Fig. 4

The transcritical bifurcation diagram depicts the exchange of stability at

Now we study the local asymptotic stability of the endemic equilibrium E 1.

The Endemic equilibrium E 1 is locally asymptotically stable for R

The jacobian matrix for the system (1) is given by where, The characteristic equation of the system (1) around its endemic equilibrium E 2 is where

It is clear from the equation (5) that first two roots are negative real numbers and remaining roots are the roots of the cubic polynomial. It is also observe that here C 1, C 2, C 3 and all are positive for any parametric value. Hence following the Routh-Hurwitz criterion we may conclude that the system (1) is locally asymptotically stable around its endemic equilibrium E 1. □

Optimal control problem

Here we focus on the time varying control M(t), which represents the awareness due to media coverage. Now it is very important to find out a strategy which minimizes the number of infected persons as well as the associated cost. In this regard, the optimal control theory is a very powerful tool to figure out such policy. Therefore we consider the optimal control problem to minimize the objective functional. Following [27], [31], we construct the objective functional as follows: subject to the proposed model (1). The parameters c 1 and c 2 corresponds as the weight constraints for the infected population and the control respectively. Here the objective functional is linear in the control with bounded states. Therefore it can be be showed by using standard results that an optimal control and corresponding optimal states exist [8]. Now we need to find out the value of the optimal control M*(t) such thatwhere .

Here we use the Pontryagin’s Maximum Principle [8], [28], [37] to derive the necessary conditions for our optimal control and corresponding states. The Lagrangian is given by

The Hamiltonian is defined as follows

For the optimal control M*(t), there exist adjoint variables, corresponding to the states S, E, Q, I ans R such that: where the adjoint variables satisfy the transversality conditions

We minimize the Hamiltonian with respect to the control variable M*(t). Now, since the Hamiltonian is linear in the control parameter, so we consider if the optimal control is bang-bang or singular. First we find the switching function as

Now the singular control occurs when the switching function vanishes on non-trivial interval of time. Also the optimal control would take its upper bound or its lower bound according as or  > 0.

Now to investigate the singular case, we assume that on some non-trivial interval. In this case we calculate .

After some simplifications of the time derivative of we obtain

Using the equations of the system (2) and (5), we obtain

We observe that the control parameter M does not explicitly occur in the above expression, so next we calculate the second derivative with respect to time. where .

Using the state and co-state equations of systems (1) and (9), we simplify the Eq. (14) and finally obtain

The above equation can be written in the form and then we can solve the singular control as provided Φ1(t) ≠ 0 and where and

Moreover in order to satisfy the Generalized Legendre-Clebsch Condition for the singular control to be optimal, we require to be negative [25]. Therefore we summarize the control profile on a nontrivial interval in the following way:

Hence the control is optimal provided Φ1(t) < 0 and .

Numerical simulation

We study numerical results in two different cases, first for fixed control and second when the control has been applied optimally. First, we consider the values of parameters in Table 1 , for numerical simulations. Since δ is the disease induced mortality rate and d is the natural death rate, hence δ > d. Using these parameters and the initial conditions as

Table 1

Description and values/ranges of parameters.

Parameters	Description	Value/Ranges
A	Total recruitment	50
β	Disease transmission rate	[0.5, 2.3]
ρ1	Portion of S contact with E	(0, 1)
ρ2	Portion of E contact with S	(0, 1)
d	Natural death rate	0.2
b1	The rate that Q becomes S	0.25
b2	The rate that E becomes quarantine	0.8
α	The rate that E becomes I	0.3
η	The rate that I becomes R naturally	0.25
σ	The rate that E becomes R naturally	0.2
c	The rate that Q becomes I	0.12
δ	The mortality rate for I	0.25
M	Policy parameter	0.8
p	Implementation rate of policy	0.78

Fig. 2 verifies the numerical result when R 0 < 1 and in this case the solutions of the model (1) converge to the DFE, and Fig. 3 verifies when R 0 > 1 and the corresponding solutions of the model (1) converge to the EE, . These results are discussed in the Theorem 3.2 and Theorem 3.4. Further, in the Theorem 3.3, we have shown that the system (1) may undergo through a transcritical bifurcation at the threshold parametric condition Therefore in Fig. 4, we demonstrate the scenarios when and it has been observed that the model (1) possesses a transcritical bifurcation there.

Fig. 2

Dynamical behavior around DFE.

Fig. 3

Dynamical behavior around EE.

We now perform the sensitivity analysis to determine the higher and lower impact of some parameters on the basic reproduction number, R 0 and hence impact of such parameters on the dynamics of the proposed model (1). To perform the sensitivity analysis [22], we calculate the normalized forward sensitivity index of R 0 and these indices measures the relative change in a R 0 with respect to the relative change in its parameters (Table 1).

The normalized forward sensitivity index of a function, F(x 1, x 2, ⋅⋅⋅, x), for x (1 ≤ i ≤ n), is defined as: .

To find the sensitivity of R 0, we consider the parameters A, β, ρ 1, ρ 2, α, d, p, M, b 2, σ as R 0 is the functions of these parameters. The sensitivity index of R 0 with respect to the parameter β is given by . Similarly, we can find the sensitivity indices of R 0 with respect to the other parameters.

Positive index indicates that R 0 is an increasing function of the corresponding parameter and negative index implies that R 0 is a decreasing function of that parameter. For example, as , it shows that if β is increased by 10% then the R 0 is also increased by 10%. Again, as implies that 10% increment in d will decrease R 0 by 4.09%. From Table 2 , we see that the most sensitive parameter which has positive impact on R 0 are the disease transmission rate β and recruitment rate A, and the parameters ρ 1, ρ 2 produce more negative impact on R 0, as compared to the parameters d, p, M, b 2, σ. We have shown this behavior in Fig. 5 .

Table 2

Sensitivity index of R0 against parameters.

Parameters	A	β	ρ1	ρ2	α	d	p	M	b2	σ
Sensitivity index	1	1	−3.545	−7.50	−0.138	−0.386	−0.757	−0.757	−0.575	−0.144

Fig. 5

Sensitivity index of R0 against some parameters.

Table 3

Number of active, recovered, death and confirmed classes of COVID 19 cases in the states of Maharashtra, Delhi and Tamil Nadu.

Date	States of India
	Maharashtra				Delhi				Tamil Nadu
March/April	Active	Recov	Death	Conf	Active	Recov	Death	Conf	Active	Recov	Death	Conf
02	0	0	0	0	2	0	0	2	0	0	0	0
03	0	0	0	0	2	0	0	2	0	0	0	0
04	0	0	0	0	3	0	0	3	0	0	0	0
05	0	0	0	0	3	0	0	3	0	0	0	0
06	0	0	0	0	4	0	0	4	0	0	0	0
07	0	0	0	0	4	0	0	4	1	0	0	1
08	0	0	0	0	4	0	0	4	1	0	0	1
09	2	0	0	2	4	0	0	4	1	0	0	1
10	5	0	0	5	5	0	0	5	1	0	0	1
11	11	0	0	11	6	0	0	6	1	0	0	1
12	11	0	0	11	8	0	0	8	1	0	0	1
13	19	0	0	19	8	0	0	8	1	0	0	1
14	31	0	0	31	9	0	0	9	1	0	0	1
15	33	0	0	33	9	0	0	9	1	0	0	1
16	39	0	0	39	9	0	0	9	1	0	0	1
17	41	0	0	41	10	0	0	10	1	0	0	1
18	45	0	0	45	12	0	0	12	2	0	0	2
19	48	0	0	48	14	0	0	14	3	0	0	3
20	52	0	0	52	14	0	0	14	3	0	0	3
21	64	0	0	64	18	0	0	18	6	0	0	6
22	74	0	0	74	26	0	0	26	7	0	0	7
23	95	0	0	95	29	0	0	29	8	0	0	8
24	104	0	0	104	30	0	0	30	14	0	0	14
25	124	1	3	128	35	0	0	35	21	0	0	21
26	120	6	3	129	38	0	0	38	27	0	0	27
27	111	15	4	130	38	0	0	38	36	0	0	36
28	150	25	5	180	38	0	0	38	39	2	0	41
29	155	25	6	186	47	0	0	47	47	3	0	50
30	165	25	8	198	95	0	0	95	62	4	0	66
31	168	39	9	216	95	0	0	95	117	6	0	123
01	254	39	9	302	144	6	2	152	227	6	1	234
02	280	42	13	335	207	8	4	219	227	6	1	234
03	277	42	16	335	207	8	4	219	302	6	1	309
04	424	42	24	490	424	15	6	445	403	6	2	411
05	424	42	24	490	478	18	7	503	476	6	3	485
06	647	56	45	748	497	19	7	523	558	8	5	571
07	764	56	48	868	548	21	7	576	608	8	5	621
08	875	79	64	1018	546	21	9	576	664	19	7	690
09	946	117	72	1135	639	21	9	669	709	21	8	738
10	1142	125	97	1364	660	25	13	698	805	21	8	834
11	1276	188	110	1574	684	25	14	723	859	44	8	911
12	1426	208	127	1761	1025	25	19	1069	915	44	10	969
13	1619	217	149	1985	1103	27	24	1154	1014	50	11	1075
14	1948	229	160	2337	1452	30	28	1510	1104	58	11	1173
15	2250	259	178	2687	1501	30	30	1561	1111	81	12	1204
16	2437	295	187	2919	1504	42	32	1578	1110	118	14	1242
17	2711	300	194	3205	1551	51	38	1640	1072	180	15	1267
18	2791	331	201	3323	1593	72	42	1707	1025	283	15	1323
19	3075	365	211	3651	1779	72	42	1893	992	365	15	1372
20	3473	507	223	4203	1668	290	45	2003	1051	411	15	1477
21	3865	572	232	4669	1603	431	47	2081	1046	457	17	1520
22	4248	722	251	5221	1498	611	47	2156	943	635	18	1596
23	4594	789	269	5652	1476	724	48	2248	949	642	18	1629
24	5307	840	283	6430	1518	808	50	2376	911	752	20	1683
25	5559	957	301	6817	1604	857	53	2514	867	866	22	1755
26	6229	1076	323	8068	1702	869	54	2625	838	960	23	1821
27	6538	1188	342	8068	1987	877	54	2918	841	1020	24	1885

Here, in table 3 the phrases ‘Recov’ and ‘Conf’ represents recovered and confirmed infected class respectively.

Figs. 6 and 7 represent the contour plot of the basic reproduction number R 0 with respect to the parameters β vs ρ 1 and ρ 1 vs ρ 2 respectively. It can be noticed from Fig. 6 that for increasing value of β and decreasing value of ρ 1, R 0 increases. Again, for decreasing values for ρ 1 and ρ 2, R 0 increases, which has been demonstrated in Fig. 7.

Fig. 6

Contour plot of R0 as a function of β and ρ1.

Fig. 7

Contour plot of R0 as a function of ρ1 and ρ2.

We know that the numerical value of the basic reproduction number R 0 determines the exact nature of secondary infections due to the disease. From the Table 1 and Fig. 5, Fig. 6, Fig. 7; we claim that the parameters with negative indices can be increased then the numerical value of R 0 would be decreased significantly. Also the parameters namely A and β has positive indices and hence reduction of these two parameters would ultimately control the disease. Further, from the sensitivity indecies, we see that the two parameters related to isolation namely ρ 1 and ρ 2 are the most sensitive parameters and one can control easily these two parameters. Thus we may conclude from here that the best procedure to control the disease is to increase the value of isolation.

Next we present the numerical solution of the formulated optimal control problem. To find the numerical solution of the objective functional (6) subject to the system of differential Eq. (1), we use fourth order Runge-kutta method with the help of MATLAB 16 software. First, we find the solution of the model (1) by forward fourth order Runge-Kutta method with an initial size of population. Then solve the representing adjoint Eq. (9) by backward fourth order Runge-Kutta method by applying the solutions of the system (1) and the transversality conditions (10). Here, we choose the following parameters set:

And the corresponding initial guess of population size is . Fig. 8 depicts the solutions of all state variables in the presence and absence of control parameter. It is noticed that in presence of control, exposed, quarantine and infected population increases slower than the absence of the control. Similarly, in presence of control, susceptible population decreases faster than the absence of the control and recovered population increases faster than the absence of the control. These real phenomena justify our proposed model. In Fig. 9 , we plot the evolution of the variations of adjoint or costate variables and the control parameter are depicted in the Fig. 10 . Clearly the optimal control follows a combination of singular and bang-bang control and this type of representation of an optimal control is a very rare phenomenon in theoretical epidemiology.

Fig. 8

Variation of population in presence and absence of control strategy.

Fig. 9

Variation of the adjoint variables when control applied optimally.

Fig. 10

Variation of the control strategy of control parameter M.

Application of model to forecast COVID-19 at some states in India

In this portion, we present an example of a real case where we can apply the proposed model. The World Health Organization (WHO) first confirmed an incident of COVID-19 on 12th January 2020 that was the induce respiratory illness in a cluster of people in Wuhan, Hubei Province, China, and this case was actually reported to the WHO on 31st December 2019. India, the second most populated country, has also been affected severely due to the global pandemic COVID-19. In this paper, we analyze the effect of COVID-19 cases in three states in three different parts of India. These states are the Western state Maharashtra (where the commercial capital of India, Mumbai situated), the northern state Delhi (where the capital of India, New Delhi is located), and the southern state Tamil Nadu (where another highly populated mega-city, Chennai is situated) and the some of total populations of these three states crossing two hundred millions. We apply the model system (1) to study the pandemic situation due to COVID-19 in these three states of India from 2nd March, 2020 to 27th April, 2020 and predict the future behavior of the disease in a short term basis. Although the first case of COVID-19 in India was confirmed in the state Kerala on 30th January 2020, the first cases of COVID-19 were reported in the state of Delhi, Tamil Nadu, and Maharashtra respectively on 2nd, 7th, and 9th March 2020. The active, recovered, death and confirmed COVID-19 cases of these three states are collected from the official websites of the official updates of coronavirus, COVID-19 in India, Government of India [36], Government of Delhi [11], Government of Maharashtra [10] and Government of Tamil Nadu[14] and presented during this period in the Table 3.

We fit the proposed model (1) to the daily active infected and confirmed (cumulative) infected COVID-19 cases in those three states of India using the set of parameters as given in Table 4 and the initial size of the population from the Table 5 . To fit these real data, we use the software Mathematica and then predict the behavior of COVID-19 for those three states on a short term basis. In Figs. 11–Fig. 13, we respectively present the active COVID-19 cases in Maharashtra, Delhi, and Tamil Nadu for 91 days starting from 2nd March, 2020, till the 31st May 2020. Also, in Figs. 14, 15 and in Fig. 16, we present the cumulative confirmed (i.e., the sum of active cases, recovered and death) COVID-19 cases of Maharashtra, Delhi, and Tamil Nadu, respectively, for the same period.

Table 4

Values of parameters in the states of Maharashtra, Delhi and Tamil Nadu.

Parameter Value	Maharashtra	Delhi	Tamil Nadu	References
A	28112	2144	14205	[34]
d	0.0058	0.0036	0.0071	[34]
δ	0.0685	0.0232	0.0122	[10], [11], [14], [36]
α	1/5.2	1/4	1/5.4	[10], [11], [14], [36]
η	1/14	1/14	1/14	[32], [36][17]
ρ1	0.64	0.72	0.62	Estimated
ρ2	0.78	0.82	0.75	Estimated
b1	0.07122	0.045185	0.062856	Estimated
b2	1.11013	0.78529	0.157832	Estimated
σ	0.119732	0.14029	0.36186	Estimated
p	0.96657	0.923769	0.90076	Estimated
M	0.916499	0.925161	0.91462	Estimated

Table 5

Initial values of population in the states of Maharashtra, Delhi and Tamil Nadu.

Initial Values	Maharashtra	Delhi	Tamil Nadu	References
S(0)	1374333	134555	947030	[10], [11], [14], [36]
E(0)	312	214	346	Estimated
Q(0)	945	845	710	Estimated
I(0)	0	2	0	[10], [11], [14], [36]
R(0)	0	0	0	[10], [11], [14], [36]

Fig. 11

Active COVID-19 cases in Maharashtra.

Fig. 13

Active COVID-19 cases in Tamil Nadu.

Fig. 14

Confirmed COVID-19 cases in Maharashtra.

Fig. 15

Confirmed COVID-19 cases in Delhi.

Fig. 16

Confirmed COVID-19 cases in Tamil Nadu.

Discussion and conclusion

After world war II, the globe has never undergone like the present scenarios that arrived due to the pandemic COVID-19. Starting from Wuhan city in December 2019, the virus SARS-CoV-2 has spread throughout as many as 210 countries and territories and continues to increase its pandemic nature. In the absence of proper vaccination and treatment of COVID-19, one would rely on qualitative control of the disease rather than complete eradication. Since these types of situations are unknown to present society, in these circumstances, various governments have taken some policies to control the destructive nature of the disease as much as possible. To study the dynamics of the disease, in this article, we have proposed and analyzed a classical SEIR type mathematical model to incorporate the COVID-19 scenarios in the system (1). A detailed analysis shows that the proposed system posses two equilibria, namely one disease-free and one endemic whose existence and asymptotic stability criteria depend on the numerical value of basic reproduction number R 0. We have also established that the proposed system (1) undergoes a transcritical bifurcation at the threshold To study the impact of the government control measures to prevent the extensive transmission of COVID-19, we have introduced a control parameter M. Further, to reduce the infected individuals as well as to minimize the cost of implementing government control measures, an objective functional has been formed and solved using Pontryagin’s maximum principle. The optimal control follows a combination of bang-bang and singular control during its application. We also have discussed the sensitivity index of the threshold parameter R 0 and found the most sensitive parameter, which has a positive impact on R 0 is the disease transmission rate.

Next, we consider three cases of populations, namely (i) active cases of COVID-19, (ii) confirmed cases of COVID-19, and (iii) recovered cases of COVID-19 in three states of India, namely, western state Maharashtra, northern state Delhi, and southern state Tamil Nadu where there are more than two hundred people residing. In Table 3, we provide the data of COVID-19 confirmed, active, recovered, and death cases of these three states starting from 2nd March 2020, till 27th April 2020. In Tables 4 and 5, we estimate the parametric values associated with the model system (1). Finally, using the software Mathematica, we try to fit our model (1) to estimate the nature of COVID-19 in those three states of India. In Fig. 11, Fig. 12, Fig. 13 , we present the active COVID-19 cases in Maharashtra, Delhi, and Tamil Nadu, respectively, for 91 days starting from 2nd March 2020, till the 31st of May 2020. In particular, we use the existing data for the first 57 days (till 27th April 2020) collecting from the websites of the Government of India and WHO to forecast the probable active COVID-19 cases in those three states for the rest of 34 days. It should be noted that here we can forecast the nature of the COVID-19 for the short term only as the Governmental policy would change in time to time resulting in the corresponding changes in the associated parameters of the proposed model system. From those three figures, we can claim that on 31st May 2020, the active cases of the three states Maharashtra, Delhi, and Tamil Nadu, will be around 46183, 3344, and 411 respectively. Similarly, in Fig. 14, Fig. 15 and in Fig. 16, we respectively present the cumulative confirmed (i.e., the sum of active, recovered and death) COVID-19 cases of Maharashtra, Delhi, and Tamil Nadu starting from 2nd March 2020, till the 31st May 2020. Further, using the existing data for the first 57 days (till 27th April 2020), the cases for the next 34 days has been predicted. It has been forecasting from those three figures that till 31st May 2020, the cumulative confirmed instances of those three states Maharashtra, Delhi, and Tamil Nadu will be around 83958, 19667, and 5133 respectively. Moreover, from above figures (Fig. 11 to 16), we can estimate the number of active cases and number of cumulative confirmed cases of those three states at any day of May, 2020.

Fig. 12

Active COVID-19 cases in Delhi

The primary finding of this article is that we have derived a mathematical model that can be used to study the qualitative dynamics of COVID-19. The basic reproduction number and its sensitivity analysis would determine the controlling procedure of the disease. Also, we have incorporated the governmental policy in our mathematical model and proposed a linear objective functional considering that governmental policy as the time-dependent control parameter. On the other hand, it should be noted that the estimated values of COVID-19 in the three different states of India in three different parts are forecasted using existing parametric space only. In the case of the first two states, Maharashtra and Delhi we see that the graph shows an increasing trend whereas for the state Tamil Nadu, the predicted data shows that the disease can be controlled with the existing parametric space. The main two burdens behind the control of COVID-19 are (i) the unconsciousness of the people about the disease and (ii) the high population density exposed the infection when the common people are out for their essential commodities. Moreover, the weakness of forecasting cumulative COVID-19 cases in the three states of India is that here the forecast is done on totally some existing parametric conditions. But human behaviour is the most uncertain phenomena. Hence, if the corresponding parametric space has been altered then there may be some changes in the graphs of COVID-19 cases. Therefore, in this article, we are targeting to forecast the COVID-19 cases on a short term basis and henceforth there would be very little chance to change in corresponding parametric space. But the framework of the current model gives some important insights into the dynamics of COVID-19 spread and control. Further, our simulation works suggest that both the quarantine and governmental intervention strategies like lockdown, media coverage on social distancing and public hygiene can play an important role in diminishing COVID-19 transmission. However, the successes of these strategies mainly rely on the proper implementation of the process. In our future work, we need to extend our model by incorporating some essential biological as well as epidemiological factors. To be given some examples, we may mention the acquired immunity in humans, differential susceptibility and infectivity of humans to COVID-19 infection, and spatial heterogeneity.

CRediT authorship contribution statement

Manotosh Mandal: Conceptualization, Methodology, Software, Formal analysis, Investigation, Validation, Writing - original draft. Soovoojeet Jana: Conceptualization, Validation, Methodology, Writing - original draft, Writing - review & editing, Visualization. Swapan Kumar Nandi: Validation, Formal analysis, Investigation. Anupam Khatua: Investigation, Resources, Data curation. Sayani Adak: Resources, Data curation. T.K. Kar: Conceptualization, Writing - review & editing, Visualization, Supervision.

Declaration of Competing Interest

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