Mathematical modeling of COVID-19 transmission: the roles of intervention strategies and lockdown

1.   Introduction

Many cases of common cold (flu) are due to different coronaviruses. However, in the past two decades, two of these viruses have left their impact at large scale: (ⅰ) Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV) in 2002 (ⅱ) Middle East Respiratory Syndrome Coronavirus (MERS-CoV) in 2012. Infection with these viruses have become serious because they are known to be evolved with animals then jump to humans via intermediate hosts ^[1]. For example, in SARS-CoV, the host animal was Palm Civet while intermediate animal was Raccon Dog. On 12th December 2019, 27 new cases of viral pneumonia including seven of them being critically ill, were reported by Wuhan Municipal Health Commission (WMHC). Most of them had a recent history of exposure to wildlife animals at the Huanan Seafood Wholesale Market in Wuhan, China, where snake, poultry, bats, and farm animals were also sold ^[2]. This was recognized to be a cause of a new type of coronavirus, officially named COVID-19 by the WHO ^[3]. In case of current global pandemic COVID-19, coronavirus is believed to jump from Bat to Pangolin to Human. The coronavisues are important pathogens of mammals (including humans) and birds and belong to the Coronaviridae family of enveloped, capped, positive sense (single stranded) RNA viruses ^[4]. Due to positive sense, the Coronavirus can easily replicate in host cells causing fast spread of severe COVID-19 cases. A typical genetic structure of Coronavirus is made of several parts containing different types of proteins: nucleocapsid, membrane and spike proteins ^[5,6,7,8]. Among all these proteins, spike proteins are very important in case of COVID-19, because these spike proteins are responsible for the interactions with human cells and it also help into entry process to the host cells.

The COVID-19 outbreak has resulted 7,082,263 confirmed cases and 405,081 deaths in 213 countries and territories all over the world by June 7, 2020. According to WHO, most of people who are infected with COVID-19 virus experience mild to moderate respiratory illness and recover without needing special treatment. Older people or who have underlying medical problems, for example, diabetes, cardiovascular disease, cancer, and chronic respiratory disease are more severe to develop the illness. Common symptoms of COVID-19 are dry cough, fever, tiredness, sore throat, aches, and shortness of breath. In general, many times it is possible to have infection without any symptoms like parasitic infection by Cryptosporidia. Also, in case of current pandemic COVID-19, according to New York times ^[9], some individuals who are infected with the coronavirus can spread it even though they have no symptoms. It is also called incubation period (time from exposure to the development of symptoms), reported between 2–14 days ^[10]. However, incubation period found as long as 27 days ^[11]. But in the case of COVID-19 spread, it is more reasonable to call it as asymptomatic stage. According to Center for disease control and Prevention (CDC) ^[12], after the incubation period, illness causes mild symptoms to severe conditions and death. People who are healthy or have mild symptoms should keep themselves in self-quarantine and contact COVID-19 information line for guidance on testing and referral. People with cough, fever or difficulty in breathing should seek medical treatment.

Coronavirus is a single-stranded RNA virus and many RNA viruses have already been adopted by our body like HIV. Therefore, sometimes drugs used for the diseases with RNA virus are also being used in case of COVID-19. However, the treatment strategies like, antibiotics and other medicines are just supporting systems not any specific treatment for COVID-19. The good point is that throughout the world, different research and development section's researchers/doctors/consultants are working and performing different kinds of experiments to investigate specific treatment/vaccine particularly targeted to COVID-19. However, the development of vaccine/medicine must take sufficient time due to different steps involved in the drug/vaccine development procedure. Therefore, at present, the key question is how to prevent/control the spread of COVID-19 until there is no treatment/vaccines available. In general, there are several ways for minimizing the transmission rate of a particular epidemic like, vaccination, use of antibiotics, maintaining proper hygiene, avoiding crowded places, washing hands, wearing protective masks, awareness programs etc. In case of COVID-19, social distancing has emerged as one of the most broadly adopted intervention strategies (e.g., social distancing, quarantine of exposed individuals, isolation of infected individuals, promoting social consensus on self-protection like wearing a face mask in the public area, washing hands regularly etc.) to reduce the infection risk/transmission rate and control the spread of COVID-19 through the reduction of social contacts. More importantly, at this stage of the outbreak, it is important to understand the human to human transmission and deployment of different control strategies such as quarantine of susceptible individuals, isolation of seriously ill individuals, quarantine of asymptomatic individuals, etc. Human-to-human transmission of COVID-19 virus was confirmed by Chan et al. ^[13], who had reported a case of five patients in a family cluster.

Isolation and quarantine are two important measures by which asymptomatic or infected individuals could be detached from the population to stop further spread of the disease. Quarantine is generally used for seemingly healthy but possibly infected individuals, while isolation applies to already infected individuals. Quarantine was also applied as one of the effective intervention strategies during the SARS epidemic of 2002–2003 ^[14]. {More than this, to control the outbreak of COVID-19, different governments are actively restricting the movement of people by imposing lockdown, which may be known as one of the largest quarantine in history. The central government of India implemented a 14-hours public curfew on March 22, 2020, after that, the Prime Minister of India ordered a complete lockdown (Phase 1) of 21 days on March 24, 2020. At the end of first lockdown period, the Indian government recommended the extension of lockdown (Phase 2) until May 3, 2020, with some conditional relaxation after April 20, 2020. Furthermore, on May 1, 2020 the lockdown (Phase 3) was again extended for two weeks until May 17, 2020. In between, the government also divided the nation into three zones- red, orange and green, with some relaxation accordingly. Finally, on May 17, 2020, lockdown (Phase 4) has further been extended upto May 31, 2020. After that government has announced lockdown (Phase 5) only for containment zones from June 1, 2020 to June 30, 2020 (it is termed as Unlock 1). Except lockdown, Indian goverment along with different countries/territories is also adopting various steps and imposing different types of intervention strategies, for instance, social distancing, washing hands for at least 30 seconds, wearing masks on public places, tracing close contacts. Therefore, in respect of transmission control, investigation of role of different intervention strategies remain an important problem.

The transmission potential is often measured in terms of the basic reproduction numbers, the outbreak peak value, time and duration under current and evolving intervention measures ^[15,16]. The basic reproduction number ($R_0$) is described as the expected number of secondary infections appearing from a single infectious individual throughout his/her whole infectious period, in the entire susceptible population ^[15,17,18]. In the study of epidemiology, the fundamental concept of reproduction number ($R_0$) is one of the most valuable ideas that the mathematical thinking has conveyed to epidemic theory ^[15]. Most significantly, $R_0$ frequently used as a threshold value that forecasts whether the disease will die out or spread. Estimation of $R_0$ by mathematical modeling can be effective for determining the potential and severity of an outbreak and providing crucial information for recognizing the disease intensity and interventions. From the explanation of $R_0$, it is obvious that if $R_0 < 1,$ each infected individual makes, on average, less than one new infected individual, and we can predict that the infection will die out from the population. If $R_0 > 1,$ the pathogen is capable to attack the susceptible population and human-to-human transmission with continued transmission chains will appear. In an epidemic disease, it could be determined that which control measures (intervention strategies) would be most helpful for suppressing $R_0$ below one and which may also provide important advice for public health initiatives. More importantly, the $R_0$ is also called a controlled reproduction number when it depends on the control strategies, it is computed for mathematical models including control strategies ^[16]. In case of COVID-19, the estimation of $R_0$ is different from different research teams and has continuously been updated as new information arises time to time. Using early evidence, WHO estimated $R_0$ to be between 1.4–2.5 ^[20]. In another preliminary study done by Zhao et al. ^[21], the $R_0$ has been estimated between 3.6–4.0 and 2.24–3.58. Further, $R_0$ has also been estimated between 1.5–3.5 in ^[22,23,24].

In 1927, Kermack and Mckendrick ^[25] introduced a prominent compartmental model to analyze the plague disease in Mumbai and succeeded in revealing its epidemiology. After that mathematical modeling have been playing a significant role in analyzing the spread and control of different infectious diseases ^[26,27]. A number of compartmental models have already been proposed and analyzed for the COVID-19 outbreak in different countries ^{[28,29,30,31,32,33,34,35,36]}. In particular, Yang et al. ^[29] proposed a mathematical model for COVID-19 incorporating multiple transmission pathways, including both human-to-human and environment-to-human transmission routes. The authors employed a bilinear incidence rate based on the law of mass action and fitted the model with the data of Wuhan city of China and estimated the reproduction number. They also found that the contribution of the environmental reservoir is significant in shaping the overall disease risk. Their results also indicate that the COVID-19 infection remain endemic, which necessitates intervention programs and long-term disease prevention policies. Li et al. ^[30] proposed an SEIQR difference-equation model for COVID-19 outbreak of Shanxi Province of China and made predictions, risk analysis and performed assessment. The authors also revealed the effects of city lockdown date on the final scale of cases. More importantly, it was also found that a little earlier lockdown in Wuhan city may result the fewer infectious cases in Shanxi and other nearby provinces. Ngonghala et al. ^[31] developed a mathematical model of COVID-19 pandemic in US (particularly, in New York) for assessing the population-level impact of the mitigation strategies. The authors performed the rigorous analysis of the model and the impacts of non-pharmaceutical intervention strategies, social distancing, quarantine, contact-tracing, isolation, face mask, etc. Eikenberry et al. ^[32] proposed a compartmental model of COVID-19 to estimate the community-wide effect of mask used by the general, asymptomatic public, a portion of which may be asymptotically infectious. The authors used the data of US states of New York and Washington for model simulation. The findings indicate that use of mask decreases the effective transmission rate, reduce the community transmission, decreases deaths and peak hospitalizations. Garba et al. ^[34] proposed a compartmental model to analyze the dynamics of COVID-19 in South Africa. The model system in ^[34] was used to estimate the effect of mitigation strategies and various control. The results of this particular study was twofold: (ⅰ) the disease may die out if control measures are implemented early and for a sustainable period of time (ⅱ) effectiveness of self-isolation reduces the number of cases.

Several compartmental models of COVID-19 outbreak in India, have also been studied ^{[37,38,39,41]}. Khajanchi et al. ^[37] proposed a compartmental model with quarantine for the transmission dynamics of COVID-19 and calibrated the model with daily and cumulative cases for the four provinces of India. The authors have performed a detailed theoretical analysis in terms of the basic reproduction number and predicted the cumulative cases. Moreover, the study suggests that quarantine, unreported and reported individuals as well as intervention policies like social distancing, lockdown, and media effect can play an important role in controlling the transmission of COVID-19. Sarkar et al. ^[38] proposed a mathematical model that predicts the dynamics of COVID-19 in India along with its 17 provinces. Their findings revealed the fact that the contact rate between susceptible and infected individuals could be reduced by a strict isolation imposed for susceptible individuals. Moreover the numerical evaluations of the model system ^[38] suggested the complete elimination of COVID-19 via suitable combination of contact tracing and restrictive social distancing. Further the authors also indicated that the accurate course of epidemic largely depends on how and when precautionary measures, isolation, and quarantine are enforced. In this direction, Sardar et al. ^[39] also considered a mathematical model on COVID-19 to analyze the impact of social distancing and lockdown. The authors have done a detailed analysis and validated the model with the data of India and its five different states. Their results suggested that the lockdown will be effective in those locations where a higher percentage of symptomatic infections exist.

As we have discussed that the government of India implemented different levels of lockdown for different time periods, e.g., full lockdown followed by lockdown with some relaxation. Therefore in Indian perspective, the analysis of different lockdown strategies on COVID-19 transmission dynamics in the presence of different intervention schemes becomes significant. In particular, the partial vs. full lockdown with other intervention strategies, like, quarantine of asymptomatic individuals and hospitalization of symptomatic individuals}. In the present work, taking care the significant role of intervention strategies, we propose a new epidemic model system with different intervention strategies of COVID-19 in a homogeneous host population to control the spread of COVID-19. The general aim of the present work is to infer significant epidemiological characteristics by investigating the role of social distancing, lockdown, quarantine, and isolation for the proposed epidemic system. Here we will mainly investigate the impact of social distancing, quarantine, and isolation via parameter estimation through the best empirical data fitting. The objectives of the present study are:

1. To analyze the transmission dynamics of COVID-19 among humans via mathematical modeling.

2. To investigate the impact of control strategies such as lockdown, quarantine and isolation to control the spread of the global pandemic (COVID-19).

3. The impact of lockdown/quarantine of susceptible individuals. What would be the outcomes if there is partial/full lockdown?

The remaining part of the paper is organized as follows: Section 2 describes the proposed system; Section 3 represents the non-negativity and boundedness; Section 4 discusses the dynamics of the model system including the reproduction number, local and global stability of disease-free and endemic equilibrium, bifurcation analysis; Section 5 presents the data fitting and parameter estimation; Section 6 demonstrates numerical simulation and sensitivity analysis; and Section 7 concludes the paper.

2.   A dynamical model

In this section, we propose a compartmental mathematical model “$SQ_1 A Q_2 IMR$” on epidemiology of COVID-19 virus infection. Here, the population is divided into seven compartments: $S$ is the compartment of susceptible individuals; $A$ is the compartment of asymptomatic individuals but having infection, these individuals also contribute to the distribution of the disease; $Q_1$ is compartment of quarantined susceptible individuals due to lockdown; $Q_2$ is the compartment of self-quarantined individuals who are asymptomatic; $I$ is the compartment of seriously ill individuals; $M$ is the compartment of individuals who are isolated via medical care; $R$ is the compartment of recovered individuals. Let $N(t)$ be the total population at time $t$ and $N(t) = S(t)+A(t)+Q_1(t)+Q_2(t)+I(t)+M(t)+R(t)$. The proposed system is given by the following system of nonlinear ordinary differential equations:

with nonnegative initial conditions $S(0) > 0, Q_1(0) \geq 0, A(0) \geq0, Q_2(0) \geq0, I(0) \geq0, M(0) \geq0, R(0) \geq0$ and $N(0) > 0.$ It is assumed that all the parameters are nonnegative. The biological interpretations of parameters and schematic diagram of the system (2.1) are given in Table 1 and Figure 1, respectively.

In system (2.1), we assume that susceptible individuals under quarantine (those in compartment $Q_1$) due to lockdown make contact with infected individuals but at a lower rate than susceptible (those in compartment $S$). This particular idea has been modeled via multiplying the effective contact rate $\beta$ by a scaling factor $\sigma$ ($0\leq \sigma \leq 1,$ where $1-\sigma$ describes lockdown effectiveness). Here, $\sigma = 0$ represents complete lockdown scenario which means susceptible individuals are under quarantine, while $\sigma = 1$ tells no lockdown situation.

To simplify our notations in system (2.1), the occupation variables associated with different compartments ($S, Q_1, A, Q_2, \; I, M, R$) are defined as the respective fractions of population ($N$) that belong to each of the corresponding compartments:

However, to avoid complications in presentation, we still write the occupation variable of compartments as $S, Q_1, A, Q_2, I, M, R$. Therefore, the system (2.1) can be rewritten in the following form:

Here, now onward, we consider the system (2.2) for further analysis.

3.   Non-negativity and boundedness of solutions

It is important to show that all the population variables are nonnegative for all $t \geq 0,$ which implies that any trajectory which starts with positive initial condition will remain positive for $t \geq 0$. It is an important feature of an epidemiological model. From the first equation of system (2.2), we have

integrating the above inequality and using initial condition, we obtain

Thus $S(t) > 0.$ Similarly, one can show that all the variables are non-negative for all $t \geq 0$.

Theorem 3.1. The closed region

is a positively invariant set for the system (2.2).

Proof. Consider $N(t) = S(t)+A(t)+Q_1(t)+Q_2(t)+I(t)+M(t)+R(t).$ Here we obtain

which yields

This implies that $\frac{dN}{dt}$ is bounded above by $\Lambda - \mu$ and bounded below by $\Lambda - (\mu + \delta).$ Now integrating the above inequality and using initial conditions, we obtain

Considering $t \rightarrow \infty$, we have

Hence the feasible region for the system (2.2) is

Hence the region $\Omega$ is positively invariant so that no solution path moves beyond the boundary of $\Omega$. Thus above theorem ensures that the proposed model is feasible both epidemiologically and mathematically.

4.   Dynamics of system (2.2)

4.1.   Disease free equilibrium and basic reproduction number

System (2.2) always has a disease-free equilibrium

Using next generation method ^[45], the basic reproduction number $R_0$ can be calculated from the relation $R_0 = \rho(V^{-1}F)$. Let $x = (A, Q_2, I, M)^{T}$, then the system (2.2) can be rewritten as

where

The Jacobian of $F(x)$ and $V(x)$ at $E^0$ are

and $F = DF(E^0)$ and $V = DV(E^0)$. Therefore $V^{-1}F$ is the next generation matrix of system (2.2) and the spectral radius of matrix $V^{-1}F$ is

According to the Theorem 2 of ^[19], the basic reproduction number $R_0$ is

Remark 4.1. For $\sigma = 0,$ system (2.2) reduces to “$SAQ_2 IMR$” system with complete lockdown. The basic reproduction number $R_1$ for the model “$SAQ_2 IMR$” is given by

From the Remark 4.1, $R_0$ can be rewritten as

4.2.   Local and global stability of $E^0$

The following theorems discuss the local and global stability of $E^0.$

Theorem 4.1. The disease free equilibrium $E^0$ is locally asymptotically stable if $R_0 < 1$ and unstable if $R_0 > 1$.

Proof. The Jacobian matrix of system (2.2) at $E^0$ is given by

It is obvious that $\lambda_1 = -\mu$ and $\lambda_2 = -\eta-\mu$ are two negative eigenvalue of $J_{E^0}$. Remaining eigenvalues are given by the following two block matrices

The eigenvalues of block matrix $J_{1E^0}$ are $\lambda_3 = -\zeta - \kappa - \mu$ and $\lambda_4 = -\mu$ and the characteristic polynomial of block matrix $J_{2E^0}$ is given by

with

It is obvious that $a_3 > 0.$ Since $R_0 < 1,$ then $a_0, a_1, a_2 > 0.$ Using the Routh-Hurwitz stability criterion ^[16], it can be shown that all the eigenvalues of matrix $J_{2E^0}$ have negative real parts i.e. $\Re(\lambda_4), \Re(\lambda_5), \Re(\lambda_6) < 0.$ If $R_0 > 1,$ then $a_0 < 0,$ thus matrix $J_{2E^0}$ has at least one eigenvalue with positive real part. Hence, disease free equilibrium ($E^0$) of the system (2.2) is locally asymptotically stable if $R_0 < 1$ and unstable if $R_0 > 1.$

Theorem 4.2. If $R_0 \leq 1, \; E^0$ is globally stable in its feasible region.

Proof. Consider the following Lyapunov function:

The derivative of $L$ along the solution of system (2.2) gives

It is easy to verify that the maximal compact invariant set in $\left\lbrace (S, Q_1, A, Q_2, I, M, R) \in \Omega : \frac{dL}{dt} = 0 \right\rbrace$ is $\left\lbrace E^0 \right\rbrace$ when $R_0 \leq 1$. Hence from the LaSalle invarience principle ^[43], $E^0$ is globally stable if $R_0 \leq 1$.

4.3.   Existence of endemic equilibrium

In this subsection, we investigate the existence of endemic equilibrium $E^* = (S^*, Q_1^*, A^*, Q_2^*, I^*, M^*, R^*)$ of system (2.2). For this, we need to solve the following system of equations:

By solving the above system of equations, we obtain

From the second and third equation, we have $\beta (q_1q_3+(\mu +q_3+q_4) (\gamma +\delta +\mu +q_2))(S + \sigma Q_1) = (\gamma +\delta +\mu) (\mu +q_1+q_2) (\mu +q_3+q_4)$. Since $S + \sigma Q_1 < \frac{\Lambda}{\mu},$ then there is no endemic equilibrium if $\beta \Lambda (q_1q_3+(\mu +q_3+q_4) (\gamma +\delta +\mu +q_2)) \leq \mu(\gamma +\delta +\mu) (\mu +q_1+q_2) (\mu +q_3+q_4)$ i.e., $R_1 \leq 1$. For $R_1 > 1,$ the endemic equilibria can be obtained by solving the following equation:

where

Here Eq (4.3) can have zero/one/two positive roots, depending on the parameter values. For the case $0 < \sigma \leq 1$, $P_3 < 0$ if $R_0 > 1$, $P_3 = 0$ if $R_0 = 1$, $P_3 > 0$ if $R_0 < 1$. Since Eq (4.3) is quadratic, therefore if $R_0 > 1$, then Eq (4.3) has a unique positive root and there is a unique endemic equilibrium. If $R_0 = 1,$ then $P_3 = 0$ and there is unique non-zero solution of (4.3) given by $I = -\frac{P_2}{P_1}$, which is positive if and only if $P_2 < 0.$ {If $R_0 = 1, P_3 = 0,$ then

The condition $P_2 < 0$ gives

with $R_1$ is determined by Eq (4.4). Further, we have

which gives

which is not possible. Hence, if $R_0 \leq 1$, there is no endemic equilibrium.

Remark 4.2. If $\sigma = 0$ (complete lockdown), then from Eq (4.3), we can easily observe that $P_1\; = \; 0, P_2 > 0$ and $R_0 = R_1$. Hence Eq (4.3) has a unique positive root if $R_1 > 1$ and no positive root if $R_1 \leq 1$. Thus for the case of $\sigma = 0$, the system (2.2) has a unique endemic equilibrium if $R_1 > 1$ and no endemic equilibrium if $R_1 \leq 1$.

4.4.   Transcritical bifurcation

From the above discussion, we observe that the system (2.2) may undergo a transcritical bifurcation at $E^0$ when $R_0 = 1$. In this subsection, we establish the conditions on the parameters using Theorem 4.1 from Castillo-Chavez and Song ^[26] and center manifold theory ^[44]. Here we omit the variable {$R$ as $R$ does not play} any role in the remaining six equations in system (2.2). We choose $\beta$ as a bifurcation parameter. By solving $R_0 = 1,$ we obtain

It can easily be obtained that the Jacobian $J_{(E^0, \beta^*)}$ evaluated at $E^0$ and $\beta = \beta^*$ has a simple zero eigenvalue and other eigenvalues have negative sign. Hence $E^0$ is a non-hyperbolic equilibrium, when $\beta = \beta^*$. Now, we calculate a right eigenvector $W = (w_1, w_2, w_3, w_4, w_5, w_6)$ and a left eigenvector $V = (v_1, v_2, v_3, v_4, v_5, v_6)$ associated to the zero eigenvalues. Here

Now from Theorem 4.1 of ^[26], we need to calculate the bifurcation constants $a$ and $b$. Choosing $f_3$ as third equation in system (2.2) and calculating partial derivatives of $f_3$ (evaluated at $E^0$, $x_1\; = \; S, x_2\; = \; Q_1, x_3 = A, x_4 = Q_2, x_5 = I, x_6 = M$), we obtain

Here $w_1$, $w_2$ are negative and $w_3$, $w_5$ are positive so that $a$ is negative and $b$ is positive. Therefore, from Theorem 4.1 of ^[26], $E^0$ changes its stability from stable to unstable at $\beta = \beta^*$ and there exists a positive equilibrium as $\beta$ crosses its critical value. Hence the system (2.2) undergoes a transcritical bifurcation at $\beta = \beta^*$. Thus the transmission rate $\beta$ plays an important role in the disease spread. If $\beta$ is less than the critical value then its easy to control the disease but if the transmission rate $\beta$ is above the critical value then the society will experience endemic disease spreading.

4.5.   Global stability of endemic equilibrium

In this subsection, we investigate the global stability ^[40] of endemic equilibrium by constructing a suitable Lyapunov function when $R_0 > 1.$ Define the Lyapunov function:

Clearly, $L(E^*) = 0$ and $L(E)\neq 0$ when $E\neq E^*$. The upper right derivative of $L(S, Q_1, A, Q_2, I, M, R)$ is given by

In the above equation, there are 14 types of situations depending on the size of $S$ and $S^*$, $Q_1$ and $Q_1^*$, $A$ and $A^*$, $Q_2$ and $Q_2^*$, $I$ and $I^*$, $M$ and $M^*$, $R$ and $R^*$. It is sufficient to analyze for $S > S^*, Q_1 > Q_1^*, A > A^*, Q_2 > Q_2^*, I > I^*, M > M^*, R > R^*$, similarly one can do for the other cases. Here

Integrating from $t_0$ to $t$ both sides, we obtain

The boundedness of $S, Q_1, A, Q_2, I, M, R$ implies that all $S, Q_1, A, Q_2, I, M, R$ have bounded derivatives for $[t_0, +\infty)$. Therefore, $L$ is uniformly continuous on $[t_0, +\infty)$. From Barbalat's Lemma ^[46], $\lim\limits_{t \rightarrow +\infty}L(t) = 0$. Hence $D^{+}L < -\mu L < 0,$ which implies that $E^*$ is globally stable.

5.   Data fitting, cumulative and daily new cases of India

In this section, we fit data by least-square approach by taking the daily new cases and cumulative cases of COVID-19 in India. We have collected data for the period April 1, 2020 to June 7, 2020, i.e. data for 68 days from Coronavirus-Worldometer ^[42]. Using system (2.2), we simulate cumulative number of cases for the period April 1, 2020 to June 7, 2020 in India. In our paper, $Z(t)$ represents cumulative number of cases, where

We assume and calculate some parameters and initial conditions except $\beta$ and $q_{3}$ by least-square approach, which are shown in Tables 2 and 4. In order to estimate the values of parameters $\beta$ and $q_{3}$, we use extensive Markov-chain Monte-Carlo (MCMC) simulations based on the adaptive combination Delayed Rejection and Adaptive Metropolis (DRAM) algorithm ^[30,48,49] for system (2.2). Using 100000 sample realizations, we can acquire the parameter values for $\beta$ and $q_{3}$ with MCMC chain histogram and the time evolution of both infection cases and comparison with the confirmed cases of COVID-19 and cases with three different stages in Figure 2. We further compute the mean values, the standard deviation, Geweke values of $\beta$ and $q_{3}$ and the mean value and standard deviation of the reproduction number for these three different stages, which are shown in Table 3.

On the other hand, the fitness of the model system (2.2) can be verified by computing the residual. The residuals are defined as

where $Y_j$ is the daily cumulative infection data and $I(t_j)$ is the model predictive cumulative data of the same day. If the residuals are small and randomly distributed, then we can say that the fitness is reasonably good.

The system (2.2) is fitted with respect to the cumulative number of cases in India in Figure 3(a). Figure 3(b) represents the residuals of fit in Figure 3(a). In Figure 4(a)–(b), we have performed the current scenario of infected individuals and we make a short time prediction for cumulative number of cases and new daily cases for 140 days. It is clear from Figure 4(a) that the cumulative number of cases are increasing for 140 days. From Figure 4(b), it is clear that the burden of disease will continue in the month of June and July. Thus we can conclude that if all the conditions remain same then in the month of July, India will observe the uncontrollable number of infections.

6.   Numerical simulation and sensitivity analysis

In this section, we perform global sensitivity analysis and numerical simulation to support our analytical results. The values of parameters are given in Table 2 and it is to be noted that the unit of parameters (rate constants) is per day. From the viewpoints of biological significances, $R_0$ plays a significant role in determining the severity (burden of disease), outcome and process of the infection. Firstly, we perform sensitivity analysis using the methodology of Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCCs) ^[47] to investigate the dependence of $R_0$ on the parameters. We also examine that the effect of sensitivity of the parameters on the population size $Q_1$ in the presence of complete lockdown $(\sigma = 0)$, partial lockdown $(\sigma = 0.25)$ and no lockdown $(\sigma = 1).$ From Figure 5(a–c), we observe that recruitment rate $\Lambda$ and transmission rate $\beta$ are the most sensitive parameters in every situation e.g., complete/partial and no lockdown. Also from Figure 5(d), we observe that recruitment rate $\Lambda$ and transmission rate $\beta$ are the most sensitive parameters for $R_0.$ To generate the LHS matrices, we assume that all the model parameters are uniformly distributed. Then using the baseline values from Table 2, a total of 1000 simulations of the system (2.2) per LHS run were carried out.

Further we show the variation of $R_0$ with respect to transmission rate $\beta$ in Figure 6. It can easily be observed that $R_0$ increases with the value $\beta$ and after a certain value of $\beta$, $R_0$ becomes greater than 1. It implies that up to a certain value of $\beta$, disease-free equilibrium is stable and beyond that value of $\beta$, disease-free equilibrium becomes unstable. In Figures 7, 8, 9 and 10, we have demonstrated the impact of different parameters on $R_0.$ Form Figure 7, it is clear that the transmission rate $\beta$ has more impact on $R_0$ as compared to $\kappa$ which is the rate that susceptible individuals move into quarantine class $Q_1$. From Figure 8, we can see that the parameter $\beta$ has more impact on $R_0$ than $\zeta$ which is the rate that quarantine individuals move into susceptible class. From Figure 9, it is easy to ensure that parameter $q_1$ which is the rate that asymptomatic individuals become self-quarantined have more influence on $R_0$ rather than the rate at which self-quarantined individuals become infected ($q_3$). Again, from Figure 10, we obtain that the parameters $q_2$ and $q_3$ have almost equal influence on $R_0$.

From Figure 11, it is clear that when $R_0 < 1$, the system (2.2) has no endemic equilibrium and the disease-free equilibrium is stable. When $R_0 > 1$, a stable endemic equilibrium appears and the disease-free equilibrium becomes unstable, i.e. exchange of stability of the equilibriums (transcritical bifurcation) arises.

To perform the stability analysis of system (2.2), we consider the parameter values from Table 2. For $E^0$, the reproduction number $R_0 = 0.183$ and $E^0(928.5, 4.25\times10^4, 0, 0, 0, 0, 0)$. Figure (12)(a) depicts the stability of $E^0.$ For $E^*$, the reproduction number $R_0 = 2.149$ and $E^*(417.2, 5.749\times10^5, 4.6\times 10^5, 8.025\times 10^4, 9.238\times 10^4, 3.746 \times 10^4, 6.557\times 10^4)$. Figure (12)(b) depicts the stability of $E^*.$ Figure (12)(c) demonstrates the variation of infected individuals with different values of $\sigma$. We observe that if the lockdown will open then endemic level will increase (i.e., the disease burden increases). In particular, different values of $\sigma$ may be interpreting the lockdown situations with certain kind of relaxation. Since $\beta$ is most sensitive parameter for our system, therefore, we investigate the dynamics of infected individuals with different values of $\beta$. Figure (12)(d) shows that an increase in the value of transmission rate $\beta$ could lead a high endemic level.

7.   Conclusion

In this work, a compartmental epidemic model “$SQ_1 AQ_2 IMR$” for the transmission dynamics of the pandemic COVID-19 is proposed and analyzed. Since the transmission of COVID-19 virus is between human to human and the daily confirmed cases are rising day by day, therefore, the prediction of infected individuals is of utmost importance for health care arrangements and control the spread of COVID-19 virus. Taking care this important issue, we incorporate different intervention strategies such as lockdown, quarantine, isolation in our proposed system (2.2) and investigate the different dynamics with respect to these strategies. The basic reproduction number $R_0$ has been calculated for the proposed model. It is proved that the dynamics of the system totally depends on $R_0$. Based on the reproduction number $R_0$, it is revealed that whenever $R_0 < 1,$ the system (2.2) has only disease free equilibrium $E^0$ which is locally stable. By constructing a suitable Lyapunov function it has also been shown that $E^0$ is globally stable when $R_0 < 1.$ When $R_0 > 1$, the system (2.2) has a unique endemic equilibrium $E^*$ and $E^0$ becomes unstable. It has also been shown that when $R_0 = 1,$ the system (2.2) undergoes the transcritical bifurcation at $E^0$ and $E^*$ is globally stable when $R_0 > 1$.

Furthermore, we fitted data for the period of 68 days to our system by least-squares approach and predicted the cumulative cases and daily new cases in India. Model parameters have also been estimated and the residuals have been plotted for the associated fitting of data. The randomness of the residuals depicts that the fitness is good. We have also made a short time prediction for 140 days and observe that the peak of infection will not reach after 73 days from June 7, 2020 that means till the middle of August if the all restriction remain same. From the short time prediction we observe that India will experience to 6,000,000 infected individuals within 140 days (see Figure 4).

Furthermore from sensitivity analysis of $R_0$, we observe that recruitment rate $\Lambda$ and transmission rate $\beta$ are most sensitive parameter. PRCCs reveal that transmission rate $\beta$ and rates $\zeta$ and $\sigma$ are positively correlated and rates $\kappa$ and $q_1$ are negatively correlated to $R_0$. This indicates that increase of quarantine and lockdown and decrease in transmission rate will reduce $R_0$ and subsequently will reduce the disease load. We have also performed sensitivity analysis for the population $Q_1$ for different cases for lockdown such as complete lockdown $(\sigma = 0)$, partial lockdown $(\sigma = 0.25)$, and no lockdown $(\sigma = 1).$ We observe that recruitment rate $\Lambda$ and transmission rate $\beta$ are most sensitive parameters to our system (2.2). After that we have demonstrated the impact of different parameters on $R_0$. Numerical results support the fact that decrease in the transmission rate $\beta$ causes the decrease in the value of $R_0$ and after a certain level of $\beta, \; R_0$ become less than one. Also decrease in the parameter $\zeta$ (the rate at which the quarantined individuals due to lockdown move into susceptible compartment) also shows positive impact on $R_0$ as it causes decrease in the value of $R_0$. Finally, an increase in the parameter $q_1$ which is the rate that asymptomatic individuals become self-quarantined has a positive influence to decrease the value of $R_0$. Therefore, COVID-19 is controllable by reducing the contacts, increasing the efficacy of lockdown and quarantine of asymptomatic individuals.

Therefore, by analyzing all the results of the proposed system (2.2), we predict that India may face a crucial phase in near future due to pandemic COVID-19. To prevent/control the pandemic, Indian Government and population must consider some strategies other than usual lockdown/quarantine/isolation. Moreover to prevent/control the pandemic, we need to decrease the transmission and recruitment rates and to achieve this, the people should make rather less contacts with infected individuals. To ensure lesser contacts, we may also spread awareness among population about the COVID-19. One of the other strategies may be faster testing (to identify the infection quickly, more peoples to be tested) to ensure that infected individuals do not spread the disease further.

In the last few months, several mathematical models have been proposed for COVID-19 outbreak. Therefore, it is interesting and necessary to compare our study with some of relevant previous studies ^{[37,38,39,41]}. In ^[37], the authors only considered quarantine of asymptomatic individuals as intervention strategies. Sarkar et al. ^[38] considered quarantine of susceptible individuals and isolation of infected individuals. The authors have not considered quarantine of asymptomatic individuals. The authors in ^[39] have considered lockdown and hospitalization of symptomatic individuals. Further, in ^[41], the authors considered only quarantine of symptomatic individuals. However, in our work, we have considered lockdown, quarantine of asymptomatic individuals, and isolation of symptomatic individuals. We have investigated the impact of these intervention strategies on $R_0$ and have also compared the impact of no lockdown, partial lockdown and full lockdown.

Mathematical modeling with spatial effects play a significant role in characterizing and understanding the spread of a particular infectious disease. Understanding the spatial spread of COVID-19 is essential for clarifying mechanisms of transmission and targeting control interventions ^{[50,51,52,53]}. More precisely, the spatial heterogeneity affects the transmission of dynamics and spatially explicit models are more effective in evaluating control strategies. To prevent an exponential spread over India, it is important to detect the spatial spread of COVID-19. In this direction, to understand the precise impact of spatial heterogeneity on the dynamics of COVID-19, we need to build a mathematical model following multi-group or multi-patch approach with different intervention strategies (for instance, lockdown, partial lockdown or no lockdown) which will be our future work. For this, one idea may be improving the existing model systems with different patches because due to partial lockdown, people usually travel among different regions/countries.

Acknowledgments

The research work of first author (Sarita Bugalia) is supported by the Council of Scientific & Industrial Research (CSIR), India [File No. 09/1131(0025)/2018-EMR-I]. The research work of second author (Vijay Pal Bajiya) is supported by the Council of Scientific & Industrial Research (CSIR), India [File No. 09/1131(0006)/2017-EMR-I]. The research work of third author (Jai Prakash Tripathi) is supported by Science and Engineering Research Board (SERB), India [File No. ECR/2017/002786] and UGC-BSR Research Start-Up-Grant, India [No. F.30-356/2017(BSR)]. This work is also supported by National Natural Science Foundation of China under Grants (11671241, 11801398), Outstanding Young Talents Support Plan of Shanxi province, and Selective Support for Scientific and Technological Activities of Overseas Scholars of Shanxi province. The authors are grateful to the handling editor and reviewers for their helpful comments and suggestions that have improved the quality of the manuscript.

Conflict of interest

All authors declare no conflicts of interest in this paper.