Impact of cross border reverse migration in Delhi–UP region of India during COVID-19 lockdown

Shubhangi Dwivedi; Saravana Keerthana Perumal; Sumit Kumar; Samit Bhattacharyya; Nitu Kumari

doi:10.1515/cmb-2022-0151

Article Open Access

Impact of cross border reverse migration in Delhi–UP region of India during COVID-19 lockdown

Shubhangi Dwivedi , Saravana Keerthana Perumal , Sumit Kumar , Samit Bhattacharyya and Nitu Kumari

Published/Copyright: July 19, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 11 Issue 1

Abstract

The declaration of a nationwide lockdown in India led to millions of migrant workers, particularly from Uttar Pradesh (UP) and Bihar, returning to their home states without proper transportation and social distancing from cities such as Delhi, Mumbai, and Hyderabad. This unforeseen migration and social mixing accelerated the transmission of diseases across the country. To analyze the impact of reverse migration on disease progression, we have developed a disease transmission model for the neighboring Indian states of Delhi and UP. The model’s essential mathematical properties, including positivity, boundedness, equilibrium points (EPs), and their linear stability, as well as computation of the basic reproduction number ( R 0 ) , are studied. The mathematical analysis reveals that the model with active reverse migration cannot reach a disease-free equilibrium, indicating that the failure of restrictive mobility intervention caused by reverse migration kept the disease propagation alive. Further, PRCC analysis highlights the need for effective home isolation, better disease detection techniques, and medical interventions to curb the spread. The study estimates a significantly shorter doubling time for exponential growth of the disease in both regions. In addition, the occurrence of synchronous patterns between epidemic trajectories of the Delhi and UP regions accentuates the severe implications of migrant plight on UP’s already fragile rural health infrastructure. By using COVID-19 incidence data, we quantify key epidemiological parameters, and our scenario analyses demonstrate how different lockdown plans might have impacted disease prevalence. Based on our observations, the transmission rate has the most significant impact on COVID-19 cases. This case study exemplifies the importance of carefully considering these issues before implementing lockdowns and social isolation throughout the country to combat future outbreaks.

Keywords: SARS-CoV-2; NPIs failure; migrated labours; doubling time; phase synchronization; parameter estimation; scenario analysis; PRCC analysis

MSC 2010: 34D06; 34D20; 37M05

1 Introduction

In the past, humanity has faced various epidemics that have posed threats to our existence. The most dire situation in this context is a pandemic. In March 2020, COVID-19, caused by the SARS-CoV-2 virus, which originated in Wuhan, China, in late 2019, had spread to over 203 countries. The outbreak was characterized by symptoms such as cough, fever, fatigue, and shortness of breath, leading to severe pneumonia and fatalities [29]. Despite this, it took the World Health Organization (WHO) 4 months to officially declare it as a pandemic on March 11, 2020 [17]. As a result of this delayed declaration, planning and implementation of strict precautionary measures against the spread of the virus in India were also delayed. In addition, China’s refusal to provide accurate data on 174 positive cases it had identified made it difficult for WHO to extract meaningful information. With limited information on disease symptoms and medications, WHO published a manual on setting up treatment centers and infection screening facilities to optimize patient care in affected regions.

India, with its high population density, has a high rate of human-to-human social contact, making it challenging to control the pandemic in its early stages [4]. To prevent the spread over a larger area, the Government of India (GoI) initiated screening of travelers at the seven busiest airports in the country [9]. However, with the limited clinical information about the disease, thermal screening was the only option for detecting infected individuals. The first case of COVID-19 in India was reported on January 30, 2020, in Kerala’s Thrissur district, by a student who had returned home from Wuhan University in China. In February, the spread of the disease among international immigrants reached three cases [7]. While infection rates were relatively low in February [25], positive cases spiked significantly in mid-March, reaching 100. The primary cause for the surge in India was international travelers, many of whom were asymptomatic and undetectable through available detection techniques. Consequently, the disease spread across the country. National media reports pointed out that due to WHO’s delayed characterization of COVID-19 as a pandemic, India may have reacted too late in its decision to cut down on flights from other countries.

Due to India’s high population density and inadequate health infrastructure (with less than 1 doctor per 1,000 people), the authorities were unable to effectively respond to the unexpected situation due to a lack of testing facilities, isolation wards, and healthcare professionals [24]. As a result, strict non-pharmaceutical interventions (NPIs) such as quarantine, mobility controls, and social distancing were implemented. The Government of India (GoI) declared the first nationwide lockdown with only 4 hours’ notice, starting from midnight on March 25, 2020, and ending on midnight on April 15, 2020. However, the timeframe (March 11 to March 25, 2020) was insufficient for the authorities to educate socially and economically disadvantaged populations on precautionary measures. The lack of disease awareness and information led to a large-scale migration of workers returning to their native places (referred to as “reverse migration”), which significantly increased the number of positive cases in India. Panic buying of essentials, the Tablighi Jamaat event, and liquor sales during the lockdown also contributed to the rise in cases. Nevertheless, reverse migration posed a major challenge for India, leading to a surge of patients across the country.

This study aims to analyze the lessons that India can learn from this experience. We assess the impact of unidirectional reverse migration on disease propagation in the capital city of Delhi and the most populous state of Uttar Pradesh (UP). To date, no specific mathematical model has been studied to describe the impact of Delhi–UP migration during the first lockdown. Based on a report [13], the total population influx into Delhi from UP was reported as 45.16% from 1991–2001, which was a crucial factor behind Delhi’s economic growth. However, the disease-induced lockdown resulted in unemployment and uncertainty about livelihood, compelling migrants to cross the Delhi border to reach their native places in order to survive the ordeal.

After all this, GoI increased testing, quarantine facilities, treatment facilities, and contact tracing and extended the first lockdown three times until May 31, 2020, to mitigate the virus spread. However, these efforts fell short in suppressing the pandemic. Despite the 68-day lockdown, statistics confirmed a significant spike in these two neighbouring states and across India. Using COVID-19 data [19], a demonstration of the spiked positive cases across Delhi, UP, and India is given in Figure 1.

Figure 1

Spike in confirmed COVID-19 cases.

The statistics presented in Table 1 confirm that Uttar Pradesh (UP) had a higher case load than Delhi during the second phase of the lockdown (LD 2.0). During LD 2.0, Delhi had approximately 1.7 times the case load compared to LD 1.0, whereas UP had 3.3 times the case load compared to LD 1.0. These statistics also suggest that the spread of COVID-19, which was initially concentrated in urban areas of India until the first lockdown, likely reached rural areas of UP as migrant workers returned to their native places [22]. The rural division of UP faced significant challenges, including shortages of medical equipment and health facilities, limited medical staff, fear of mistreatment, social stigma, and quarantine concerns, while dealing with COVID-19. In this study, we take the opportunity to model the impact of migration and anticipate the potential implications of lockdown policies in India. Since the end of 2019, researchers worldwide have formulated various epidemiological models to forecast the quantitative growth of COVID-19 in different regions. In contrast to previous studies, this work mathematically analyzes the post-lockdown and post-pandemic situation in the Delhi–UP region, providing insights into disease dynamics through synchronization, which can elucidate how viruses evolve beyond their basic biological characteristics.

Table 1

Spike in number of cases in Delhi–UP region and across India [19]

Phase	Start date	End date	Days	Surge in India	Surge in Delhi	Surge in UP
LD 1.0	25 March,20	14 April, 20	21	9,714	1,480	525
LD 2.0	15 April, 20	3 May, 20	19	28,541	2,561	1,752
LD 3.0	4 May, 20	17 May, 20	14	48,394	4,435	1,492
LD 4.0	18 May, 20	31 May, 20	14	85,974	8,495	2,840

2 Model

The mathematical model used to describe the return of laborers from Delhi to UP is based on the SIR model. The model equations take into account the population of Delhi and UP, which differ in population density, economic condition, and healthcare facilities, but have similar living environments and disease transmission probabilities. The model assumes that the population is homogeneous and has homogeneous mixing. Other significant assumptions of the model are as follows:

The population is assumed to grow at a constant rate, b 0 , and decline at the natural death rate, d 0 .
The total population of each state N i ( i = 1 , 2 ) is distributed into seven compartments with
N i = S i + E i + Q i + I A i + I S i + I I i + R i .
The description of each compartment is mentioned in Table 2.
The standard incidence rate is used as β S I N , where β measures the combined effect of the infectiousness of the disease and the contact transmission rates. The incidence function is assumed to be independent of media coverage, as the role of media was negligible in communicating preventive measures to the underprivileged population in the early stage of the pandemic.
A parameter a is used to represent the fraction of susceptible population who are likely to roam around and have effective contact with infectious individuals. Susceptible individuals acquire infection at a rate of β S ( I + η I I ) a N , where 0 < η < 1 is a modification parameter that accounts for the assumed reduction in disease transmission by home-quarantined individuals. η measures the efficacy of isolation or treatment given to the home-quarantined individuals, with perfect isolation if η = 0 , leaky isolation if 0 < η < 1 , and completely ineffective isolation if η = 1 [23].
The exposed individuals are captured by the term β S ( I + η I I ) a N . These individuals are removed from the susceptible compartment and added to the exposed compartment.
The individuals in the exposed compartment remain there for the incubation period before becoming infectious. The latent period is considered as 1 ε , where ε is the rate of removal from the exposed compartment to the infectious compartment [6]. The value of ε is determined by the incubation period, which is 5 days for COVID-19.
The population in the infected and quarantine compartments increases over time but then reduces due to recovery and death. The recovered population always increases by susceptible and migrant individuals after they leave the quarantine class. However, a low disease-induced death rate δ is considered in the model [21].
Migration between states occurs with a finite speed, and there is an associated delay of Ω 1 , 2 in the migration rate m , which is proportional to the geographic distance between the states.
The recovery rate χ is calculated as 1 (recovery period) , with a recovery period of 14 days after contact with a confirmed case according to the literature [15].

Table 2

Description of compartments

Compartment	Description
Susceptible (S)	Entire population is assumed to be susceptible.
Exposed (E)	Fraction who came in contact with COVID-19-infected person
Quarantined (Q)	Fraction who were tested positive and facilitated the service from
	government quarantine centers for 14 days
Asymptomatic infected ( I A )	Fraction who had the disease but did not show symptoms
Symptomatic infected ( I S )	Fraction with the disease symptoms
Home-quarantined infected ( I I )	Fraction who isolated themselves at their home for incubation period due to limited number of government quarantine centres
Recovered ( R )	Removed population either by death or recovery

A schematic diagram for the COVID-19 mathematical model is shown in Figure 2. With the aforementioned assumptions, the model equations for reverse migration from Delhi to UP are presented as follows:

(1) S ˙ 1 ( t ) = b 1 − β 1 S 1 ( I A 1 + I S 1 + η I I 1 ) a 1 N 1 − d 1 S 1 − m S 1 ( t − Ω ) E ˙ 1 ( t ) = β 1 S 1 ( I A 1 + I S 1 + η I I 1 ) a 1 N 1 − ( ε + k 1 ) E 1 − m E 1 ( t − Ω ) Q ˙ 1 ( t ) = k 1 E 1 + α 1 I A 1 − χ Q 1 I ˙ A 1 ( t ) = k 1 ε E 1 − ( α 1 + χ ) I A 1 − m I A 1 ( t − Ω ) I ˙ S 1 ( t ) = ( 1 − k 1 ) ε E 1 − ϕ 1 I S 1 − m I S 1 ( t − Ω ) I I 1 ˙ ( t ) = ϕ 1 I S 1 − ( χ + δ 1 ) I I 1 R ˙ 1 ( t ) = χ ( I A 1 + I I 1 + Q 1 ) − d 1 R 1 − m R 1 ( t − Ω ) S ˙ 2 ( t ) = b 2 − β 2 S 2 ( I A 2 + I S 2 + η I I 2 ) a 2 N 2 − d 2 S 2 + m S 1 ( t − Ω ) E ˙ 2 ( t ) = β 2 S 2 ( I A 2 + I S 2 + η I I 2 ) a 2 N 2 − ( ε + k 2 ) E 2 + m E 1 ( t − Ω ) Q ˙ 2 ( t ) = k 2 E 2 + α 2 I A 2 − χ Q 2 I ˙ A 2 ( t ) = k 2 ε E 2 − ( α 2 + χ ) I A 2 + m I A 1 ( t − Ω ) I ˙ S 2 ( t ) = ( 1 − k 2 ) ε E 2 − ϕ 2 I S 2 + m I S 1 ( t − Ω ) I I 2 ˙ ( t ) = ϕ 2 I S 2 − ( χ + δ 2 ) I I 2 R ˙ 2 ( t ) = χ ( I A 2 + I I 2 + Q 2 ) − d 2 R 2 + m R 1 ( t − Ω ) .

For description of other parameters used in the proposed model, we provide Table 3, which lists the explanation and numerical values of the chosen parameters.

Figure 2

A schematic diagram of COVID-19 model (1). All arrows indicate transition from one compartment to another.

Table 3

Details of parameters used in the model simulation

Parameter	Description of parameters	Numerical value (sources)
b 1 , 2	Birth rate	0.000024246, 0.0000717
d 1 , 2	Death rate	0.000024246, 0.0000717
β 1 , 2	Contact rate	0.191, 0.158
a 1 , 2	Fraction of population which is susceptible	0.2, 0.3 (Calibrated)
k 1 , 2	Rate of quarantine rates of exposed individuals	0.059, 0.09
ε	Incubation period	1 5
α 1 , 2	Rate of detection through contact tracing	0.069, 0.086
χ	Rate of recovery	0.09871 [3]
ϕ 1 , 2	Rate of quarantine of symptomatic individuals	0.122, 0.111
δ 1 , 2	Disease-induced death rate	0.0026, 0.0025
Ω	Delay in traveling	4.5
m	Migration rate	0.0005

3 Well-posedness

3.1 Positive invariance

Proposition 1

The non-negative space R ≥ 0 14 is positively invariant by the system.

Proof

System (1) can be written in the following form:

(2) S ˙ ( t ) = b − S I β ( a N ) − 1 − d S − m S ( t − Ω ) E ˙ ( t ) = S I β ( a N ) − 1 − ( ε I i d + k ) E − m E ( t − Ω ) Q ˙ ( t ) = k E + α I A − χ I i d Q I ˙ A ( t ) = k ε I i d E − ( α + χ I i d ) I A − m I A ( t − Ω ) I ˙ S ( t ) = ( 1 − k ) ε I i d E − ϕ I S − m I S ( t − Ω ) I I ˙ ( t ) = ϕ I S − ( χ I i d + δ ) I I R ˙ ( t ) = χ I d ( I − η I I + Q ) − d R − m R 1 ( t − Ω ) ,

where S = [ S 1 , S 2 ] T , E = [ E 1 , E 2 ] T , Q = [ Q 1 , Q 2 ] T , I A = [ I A 1 , I A 2 ] T , I S = [ I S 1 , I S 2 ] T , I I = [ I I 1 , I I 2 ] T , and R = [ R 1 , R 2 ] T . Also a = [ a 1 , a 2 ] T , I = [ I A 1 + I S 1 + η I I 1 , I A 2 + I S 2 + η I I 2 ] T , m = ( m , − m ) T , and I i d is the identity matrix. Moreover, b = [ b 1 , b 2 ] T , N = [ N 1 , N 2 ] T , k = [ k 1 , k 2 ] T , α = [ α 1 , α 2 ] T , ϕ = [ ϕ 1 , ϕ 2 ] T , and δ = [ δ 1 , δ 2 ] T .

Again rewriting the system (2) in a more compact form:

(3) x ˙ ( t ) = b − ∑ i = 1 2 ⟨ β i ˜ ∣ y ⟩ ⟨ e i ∣ x ⟩ diag − 1 ( a N ) + A x x + B x x ( t − Ω ) y ˙ ( t ) = ∑ i = 1 2 ⟨ β i ˜ ∣ y ⟩ ⟨ e i ∣ x ⟩ diag − 1 ( a N ) + A y y + B y y ( t − Ω ) ,

where

x = ( x 1 , x 2 ) T = ( S 1 , S 2 ) T ∈ R ≥ 0 2 , x ( t − Ω ) = ( S 1 ( t − Ω ) , S 2 ( t − Ω ) ) T ,

y = ( ( y j ) j = 1 12 ) T = ( E 1 , Q 1 , I A 1 , I S 1 , I I 1 , R 1 , E 2 , Q 2 , I A 2 , I S 2 , I I 2 , R 2 ) T ∈ R ≥ 0 12 ,

e 1 = ( 1 , 0 ) , e 2 = ( 0 , 1 ) ,

β 1 ˜ = ( 0 , 0 , β 1 , β 1 , η β 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) , β 2 ˜ = ( 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , β 1 , β 1 , η β 1 , 0 ) .

A x = − d 1 0 0 − d 2 , B x = − m 0 0 m , A y = A y 1 O 6 × 6 O 6 × 6 A y 2 ,

where

A y 1 = − ( ε + k 1 ) 0 0 0 0 0 k 1 − χ α 1 0 0 0 k 1 ε 0 − ( α 1 + χ ) 0 0 0 ( 1 − k 1 ) ε 0 0 − ϕ 1 0 0 0 0 0 ϕ 1 − ( χ + δ 1 ) 0 0 0 χ χ χ − d 1 ,

A y 2 = − ( ε + k 2 ) 0 0 0 0 0 k 2 − χ α 2 0 0 0 k 2 ε 0 − ( α 2 + χ ) 0 0 0 ( 1 − k 2 ) ε 0 0 − ϕ 2 0 0 0 0 0 ϕ 2 − ( χ + δ 2 ) 0 0 0 χ χ χ − d 2 ,

and

B y = diag ( − m , 0 , − m , − m , 0 , − m , m , 0 , m , m , 0 , m )

with zero delay, system (3) can be rewritten as follows:

(4) x ˙ ( t ) = b − ∑ i = 1 2 ⟨ β i ˜ ∣ y ⟩ ⟨ e i ∣ x ⟩ diag − 1 ( a N ) + ( A x + B x ) x ( t ) y ˙ ( t ) = ∑ i = 1 2 < β i ˜ T x i diag − 1 ( a N ) + A y + B y y ,

where ∑ i = 1 2 < β i ˜ T x i diag − 1 ( a N ) + A y + B y and A x + B x are Metzler matrices since x ≥ 0 . A Metzler matrix is a matrix with off-diagonal entries non-negative [12,26]. Then, the plane R ≥ 0 2 is positively invariant by

x ˙ ( t ) = b − ∑ i = 1 2 ⟨ β i ˜ ∣ y ⟩ ⟨ e i ∣ x ⟩ diag − 1 ( a N ) + ( A x + B x ) x ( t ) ,

since it is well known that a linear Metzler system let invariant the non-negative space. This proves the positive invariance of the non-negative space R ≥ 0 12 by system 1, and this achieves the proof.□

3.2 Boundedness of trajectories

In order for the proposed system to hold epidemiological significance, it is crucial to establish that all of its state variables are non-negative and remain within a bounded range throughout the entirety of the system’s lifespan. To substantiate this claim, we propose the following proposition to demonstrate the boundedness of trajectories.

Proposition 2

For ψ ≥ 0 , the simplex:

N ψ = ( S 1 , 2 , E 1 , 2 , Q 1 , 2 , I A 1 , 2 , I S 1 , 2 , I I 1 , 2 , R 1 , 2 ) ∈ R ≥ 0 14 , N ( t ) ≤ ( b 1 + b 2 ) min ( d 1 , d 2 ) + ψ

is a compact forward invariant set for system (1).

Proof

Let the initial data be

S 1 , 2 ( 0 ) > 0 , E 1 , 2 ( 0 ) > 0 , Q 1 , 2 ( 0 ) > 0 , I A 1 , 2 ( 0 ) > 0 , I S 1 , 2 ( 0 ) > 0 , I I 1 , 2 ( 0 ) > 0 , R 1 , 2 ( 0 ) > 0 .

Then, it can be shown that the solutions

( S 1 , 2 , E 1 , 2 , Q 1 , 2 , I A 1 , 2 , I S 1 , 2 , I I 1 , 2 , R 1 , 2 )

of system (1) are positive for all t > 0 . Now, adding all equations in the differential system (1) gives

N ˙ ( t ) = b 1 − d 1 S 1 − δ 1 I I 1 − d 1 R 1 + b 2 − d 2 S 2 − δ 2 I I 2 − d 2 R 2 N ˙ ( t ) = ( b 1 + b 2 ) − d 1 S 1 − 0 d 1 Q 1 + 0 d 1 I A 1 + δ 1 d 1 I I 1 + R 1 − d 2 S 2 − 0 d 2 Q 2 + 0 d 2 I A 2 + δ 2 d 2 I I 2 + R 2 .

Since d 1 > 0 and d 2 > 0 , therefore, it can be deduced that

N ˙ ( t ) ≤ b − d N ( t ) ,

where b = b 1 + b 2 and d = min ( d 1 , d 2 ) . It then follows that lim t → ∞ N ( t ) = b d , which implies that the trajectories of system (1) are bounded.

On the other hand, solving the differential inequality N ˙ ( t ) ≤ b − d N ( t ) gives

N ( t ) ≤ N ( 0 ) exp ( − d t ) + b d ( 1 − exp ( − d t ) ) .

In particular, N ( t ) ≤ b d if N ( 0 ) ≤ b d . Thus, the defined simplex set is absorbing for ψ > 0 .□

4 Existence of disease-free and endemic equilibria

The occurrence of a disease-free equilibrium or an endemic equilibrium in a population is dependent on the value of the basic reproduction number ( R 0 ). However, it is possible to demonstrate this independently through the use of a biologically meaningful initial condition.

(5) ( S i ( 0 ) , E i ( 0 ) , Q i ( 0 ) , I A i ( 0 ) , I S i ( 0 ) , I I i ( 0 ) , R i ( 0 ) ) ∈ { ( S i , E i , Q i , I A i , I S i , I I i , R i ) ∈ [ 0 , N i ] 14 , S i > 0 , E i > 0 , Q i > 0 , I A i > 0 , I S i > 0 , I I i > 0 , R i > 0 , S i + E i + Q i + I A i + I S i + I I i + R i = N i for i ∈ { 1 , 2 } } .

(6) R 0 ≤ 1 ⇒ lim t → ∞ ( S 1 , 2 ( t ) , E 1 , 2 ( t ) , Q 1 , 2 ( t ) , I A 1 , 2 ( t ) , I S 1 , 2 ( t ) , I I 1 , 2 ( t ) , R 1 , 2 ( t ) ) = DFE , R 0 > 1 , I A 1 , 2 ( 0 ) , I S 1 , 2 ( 0 ) , I I 1 , 2 ( 0 ) > 0 ⇒ lim t → ∞ ( S 1 , 2 ( t ) , E 1 , 2 ( t ) , Q 1 , 2 ( t ) , I A 1 , 2 ( t ) , I S 1 , 2 ( t ) , I I 1 , 2 ( t ) , R 1 , 2 ( t ) ) = EE .

In this section, the governing equations of the model, which involve 14 compartments, are considered without incorporating any delay term for the purpose of simplified calculations and enhanced comprehension. Specifically, Ω is taken as 0, indicating the absence of delays. It is known that for the steady state equilibrium, we follow

d P d t = 0 ∀ P ∈ { S i , E i , Q i , I A i , I S i , I I i , R i } , for i ∈ { 1 , 2 } .

The analytic solution of equilibrium points is obtained by equating the model equations (1) to zero. There exists an unique non-trivial equilibrium point

E E E = ( S 1 ∗ , E 1 ∗ , Q 1 ∗ , I A 1 ∗ , I S 1 ∗ , I I 1 ∗ , R 1 ∗ , S 2 ∗ , E 2 ∗ , Q 2 ∗ , I A 2 ∗ , I S 2 ∗ , I I 2 ∗ , R 2 ∗ ) ,

where each component of E E E is non-zero. However, the model (1) does not possess any disease-free equilibrium for active reverse migration.

Suppose that the influx of migrants in UP from Delhi was nil, i.e., m = 0 . In this case, the model (1) possess its four non-trivial equilibrium points E j , for j ∈ { 1 , 2 , 3 , 4 } in explicit form, which are as follows:

(7) EE in Delhi , DFE in UP ( E 1 ) = ( S 1 E 1 , E 1 E 1 , Q 1 E 1 , I A 1 E 1 , I S 1 E 1 , I I 1 E 1 , R 1 E 1 , b 2 d 2 , 0 , 0 , 0 , 0 , 0 , 0 ) , DFE in both regions ( E 2 ) = b 1 d 1 , 0 , 0 , 0 , 0 , 0 , 0 , b 2 d 2 , 0 , 0 , 0 , 0 , 0 , 0 , EE in both regions ( E 3 ) = ( S 1 E 3 , E 1 E 3 , Q 1 E 3 , I A 1 E 3 , I S 1 E 3 , I I 1 E 3 , R 1 E 3 , S 2 E 3 , E 2 E 3 , Q 2 E 3 , I A 2 E 3 , I S 2 E 3 , I I 2 E 3 , R 2 E 3 ) , DFE in Delhi , EE in UP ( E 4 ) = b 1 d 1 , 0 , 0 , 0 , 0 , 0 , 0 , S 2 E 4 , E 2 E 4 , Q 2 E 4 , I A 2 E 4 , I S 2 E 4 , I I 2 E 4 , R 2 E 4 .

where components other than zero entries are non-negative. Two of the non-trivial fixed points E 2 and E 4 represent disease-free equilibrium (DFE) and endemic equilibrium (EE), respectively. There may also be other implicit equilibrium points in the model. However, their computation is complex due to the large dimension of the model.

Remark

Disease-free equilibrium for model (1) with m ≠ 0 offers initial evidence that a pandemic like COVID-19 could be suppressed by strict lockdown only.

To gain a better understanding of the dynamics and stability of the system at various equilibrium points, we plot the nullclines of the system’s trajectories with m = 0 . As m = 0 divides the model (1) into two independent patches (which are identical but with distinct parameter values), it is adequate to plot the nullclines of just one isolated patch to gain insight into how the intersection of state variables leads to an equilibrium point. Furthermore, to plot the nullclines, we draw the susceptible and exposed compartments and consider the other compartments as quasi-steady states. Figure 3 displays the nullclines for the S 1 and E 1 compartments, along with an equilibrium point resulting from the intersection of these two nullclines.

Figure 3

Nullcline plot for susceptible and exposed compartments.

4.1 Disease-free equilibrium of model with no migration

On the basis of the analysis conducted in Section 4, we consider the model equations with m = 0 to assess the stability of the disease-free equilibrium. Therefore, we analyze model (1) without the delay term, and compute the Jacobian matrix to investigate its stability around the disease-free equilibrium. The Jacobian matrix is obtained as follows:

(8) Jacobian 14 × 14 = J 7 × 7 Delhi O 7 × 7 J 1 7 × 7 UP J 7 × 7 UP ,

(9) J 7 × 7 Delhi = − ( d 1 + m ) − β 1 ( I A 1 + I S 1 + η I I 1 ) N 1 a 1 0 0 − S 1 β 1 N 1 a 1 − S 1 β 1 N 1 a 1 − S 1 β 1 η N 1 a 1 β 1 ( I A 1 + I S 1 + I I 1 η ) N 1 a 1 − ε − k 1 − m 0 S 1 β 1 N 1 a 1 S 1 β 1 N 1 a 1 S 1 β 1 η N 1 a 1 0 0 k 1 − χ α 1 0 0 0 0 ε k 1 0 − α 1 − χ − m 0 0 0 0 − ε ( k 1 − 1 ) 0 0 − ( m + ϕ 1 ) 0 0 0 0 0 ϕ 1 − ( χ + δ 1 ) 0 0 0 χ χ 0 χ − ( d 1 + m ) ,

(10) J 1 7 × 7 UP = m 0 0 0 0 0 0 0 m 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m 0 0 0 0 0 0 0 m 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m ,

and

(11) J 7 × 7 UP = − d 2 − β 2 ( I A 2 + I S 2 + η I I 2 ) N 2 a 2 0 0 − S 2 β 2 N 2 a 2 − S 2 β 2 N 2 a 2 − S 2 β 2 η N 2 a 2 β 2 ( I A 2 + I S 2 + I I 2 η ) N 2 a 2 − ( ε + k 2 ) 0 S 2 β 2 N 2 a 2 S 2 β 2 N 2 a 2 S 2 β 2 η N 2 a 2 0 0 k 2 − χ α 2 0 0 0 0 ε k 2 0 − ( α 2 + χ ) 0 0 0 0 ε ( 1 − k 2 ) 0 0 − ϕ 2 0 0 0 0 0 0 ϕ 2 − ( χ + δ 2 ) 0 0 0 χ χ 0 χ − d 2 .

It is worth noting that the matrix (8) is a lower triangular block matrix, composed of (9) and (11). Consequently, the eigenvalues of the Jacobian matrix (8) will be a combination of the eigenvalues of these individual submatrices. The stability of the equilibrium points can be verified by calculating the Eigenvalues of the Jacobian matrix (8) of equations (1) at these points.

With strict lockdown (even a single individual cannot slip in UP from Delhi), the uncoupled Jacobian matrices (9) and (11) (see Supplementary information) at the disease-free equilibrium ( b 1 d 1 , 0 , 0 , 0 , 0 , 0 , 0 , b 2 d 2 , 0 , 0 , 0 , 0 , 0 , 0 ) are as follows:

(12) J 7 × 7 Delhi , UP ∣ DFE = − d 1 , 2 0 0 0 − c 1 , 2 − c 1 , 2 − c 1 , 2 0 − ε − k 1 , 2 0 c 1 , 2 c 1 , 2 c 1 , 2 0 0 k 1 , 2 − χ α 1 , 2 0 0 0 0 ε k 1 , 2 0 − α 1 , 2 − χ 0 0 0 0 ε ( 1 − k 1 , 2 ) 0 0 − ϕ 1 , 2 0 0 0 0 0 0 ϕ 1 , 2 − ( χ + δ 1 , 2 ) 0 0 0 χ χ 0 χ − d 1 , 2 ,

where c 1 , 2 = − b 1 , 2 β 1 , 2 N 1 , 2 a 1 , 2 d 1 , 2 .

Next, model (1) with m = 0 will be stable at disease-free equilibrium only if all the e-values of J Delhi ∣ DFE and J UP ∣ DFE will have negative real part. Therefore, since uncoupled patches Delhi and UP have identical dynamics, we compute the e-values of J Delhi ∣ DFE as the same computation will work for the other patch. The characteristic polynomial of J Delhi ∣ DFE is expressed as follows:

(13) Det ( J Delhi ∣ DFE − λ I ) = − ( λ + d 1 ) 2 ( λ + χ ) [ p 4 λ 4 + p 3 λ 3 + p 2 λ 2 + p 1 λ + p 0 ] ,

where

(14) p 4 = 1 , p 3 = ( ε + a 1 + k 1 + δ 1 + 2 χ − ϕ 1 ) , p 2 = ( − ϕ 1 ( χ + δ 1 ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 − ϕ 1 ) + ( ε + k 1 ) ( a 1 + χ ) − ε c 1 ) , p 1 = ( − ϕ 1 ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( χ + δ 1 − ϕ 1 ) ( ( ε + k 1 ) ( a 1 + χ ) − ε c 1 k 1 ) − ε c 1 ( 1 − k 1 ) ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) , p 0 = − ϕ 1 ( χ + δ 1 ) ( ( ε + k 1 ) ( a 1 + χ ) − ε c 1 k 1 ) − ε c 1 ( 1 − k 1 ) ( a 1 + χ ) ( δ 1 + χ + ϕ 1 ) .

Obviously, equation (13) has three negative roots λ = − d 1 , − d 1 , − χ . Further, we consider the expression in the square brackets

(15) p 4 λ 4 + p 3 λ 3 + p 2 λ 2 + p 1 λ + p 0 .

To ensure that all roots of equation (15) have negative real parts, the Routh-Hurwitz stability criterion [8] requires

(16) p 3 > 0 , p 1 > 0 , p 0 > 0 , p 1 ( p 2 p 3 − p 1 ) > p 0 p 3 2 .

The inequalities in equation (16) hold when R 0 ≤ 1 . On the basis of subsection, we show the inequalities in equation (16) at the disease-free equilibrium.

First, we rewrite p 3 as follows:

(17) p 3 = ( ε + a 1 + k 1 + δ 1 + 2 χ ) − ϕ 1 ,

Thus, p 3 > 0 , when

(18) [ C 1 ] : ( ε + a 1 + k 1 + δ 1 + 2 χ ) > ϕ 1 .

Further, let us re-write p 1 as follows:

(19) p 1 = − ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) + ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) .

Thus, p 1 > 0 , when

[ C 2 ] : ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) > ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) .

Next, for p 0 ,

(20) p 0 = ε c 1 k 1 ( ϕ 1 + ( a 1 + χ ) ( δ 1 + χ + ϕ 1 ) ) − ( a 1 + χ ) ( ϕ 1 ( χ + δ 1 ) ( ε + k 1 ) + ε c 1 ( δ 1 + χ + ϕ 1 ) ) .

Hence, p 0 > 0 if

[ C 3 ] : ε c 1 k 1 ( ϕ 1 + ( a 1 + χ ) ( δ 1 + χ + ϕ 1 ) ) > ( a 1 + χ ) ( ϕ 1 ( χ + δ 1 ) ( ε + k 1 ) + ε c 1 ( δ 1 + χ + ϕ 1 ) ) .

As a first step toward proving the last inequality in equation (16), i.e., p 1 p 2 p 3 > p 0 p 3 2 , it is sufficient to establish the following two inequalities:

(21) p 1 p 2 p 3 > 2 p 1 2 , p 1 p 2 p 3 > 2 p 0 p 3 2 .

with condition [ C 2 ] , we check the condition of positivity for the first expression p 2 p 3 − 2 p 1 of equation (21)

(22) p 2 p 3 − 2 p 1 = 2 ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) + ( ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) ( ε + a 1 + k 1 + δ 1 + 2 χ ) + ϕ 1 ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) ) , − ( ϕ 1 ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) + ( ε + a 1 + k 1 + δ 1 + 2 χ ) × ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) ) 2 ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) , = κ 1 − κ 2

where

(23) κ 1 = 2 ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) + ( ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) ( ε + a 1 + k 1 + δ 1 + 2 χ ) + ϕ 1 ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) ) , κ 2 = ( ϕ 1 ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) + ( ε + a 1 + k 1 + δ 1 + 2 χ ) ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) ) 2 ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) .

If [ C 4 ] : κ 1 > κ 2 , then the expression p 1 p 2 p 3 > 2 p 1 2 will be positive with conditions [ C 2 ] and [ C 4 ] .

Finally, to show the second condition of equation (21), we write p 1 p 2 p 3 > 2 p 0 p 3 2 as follows:

(24) p 3 ( p 1 p 2 − 2 p 0 p 3 ) > 0 .

With condition [ C 1 ] if the expression ( p 1 p 2 − 2 p 0 p 3 ) > 0 holds, then all the condition of equation (16) will hold. Thus,

(25) p 1 p 2 − 2 p 0 p 3 = ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) × ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) + ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) + 2 ( a 1 + χ ) ( ϕ 1 ( χ + δ 1 ) ( ε + k 1 ) + ε c 1 ( δ 1 + χ + ϕ 1 ) ) × ( ε + a 1 + k 1 + δ 1 + 2 χ ) + 2 ϕ 1 ε c 1 k 1 ϕ 1 + ( a 1 + χ ) ( δ 1 + χ + ϕ 1 ) − 2 ( ε c 1 k 1 ( ϕ 1 + ( a 1 + χ ) ( δ 1 + χ + ϕ 1 ) ( ε + a 1 + k 1 + δ 1 + 2 χ ) + 2 ϕ 1 ( a 1 + χ ) ( ϕ 1 ( χ + δ 1 ) ( ε + k 1 ) + ε c 1 ( δ 1 + χ + ϕ 1 ) ) ) + ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) + ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) ) , = κ 3 − κ 4 ,

where

(26) κ 3 = ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) × ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) + ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) + 2 ( a 1 + χ ) ( ϕ 1 ( χ + δ 1 ) ( ε + k 1 ) + ε c 1 ( δ 1 + χ + ϕ 1 ) ) × ( ε + a 1 + k 1 + δ 1 + 2 χ ) + 2 ϕ 1 ( ε c 1 k 1 ( ϕ 1 + ( a 1 + χ ) ( δ 1 + χ + ϕ 1 ) , κ 4 = ( 2 ( ( ε c 1 k 1 ( ϕ 1 + ( a 1 + χ ) ( δ 1 + χ + ϕ 1 ) ( ε + a 1 + k 1 + δ 1 + 2 χ ) + 2 ϕ 1 ( a 1 + χ ) ( ϕ 1 ( χ + δ 1 ) ( ε + k 1 ) + ε c 1 ( δ 1 + χ + ϕ 1 ) ) ) + ( ϕ 1 ( ( χ + δ 1 ) ( ε + k 1 + a 1 + χ ) + ( ε + k 1 ) ( a 1 + χ ) ) + ε c 1 ( k 1 ( χ + δ 1 ) + ( a 1 + 2 χ + ϕ 1 + δ 1 ) ) ) ( ( ε + k 1 ) ( a 1 + χ ) + ( ε + k 1 + a 1 + χ ) ( χ + δ 1 ) ) + ( ϕ 1 ( 2 χ + δ 1 + ε + k 1 + a 1 ) + ε c 1 ) ( ε c 1 k 1 ( a 1 + 2 χ + 2 ϕ 1 + δ 1 ) + ( χ + δ 1 ) ( ε + k 1 ) ( a 1 + χ ) ) ) .

If [ C 5 ] : κ 3 > κ 4 , then the expression p 1 p 2 p 3 > 2 p 0 p 3 2 will be positive with conditions [ C 1 ] and [ C 5 ] . It is apparent that the isolated patch (Delhi) of the model with m = 0 is stable around disease-free equilibrium point under conditions [ C j ] ∀ j ∈ { 0 , 1 , 2 , 3 , 4 , 5 } .Thus, the disease would have eradicated if all the conditions derived below were satisfied.

Note. Using a similar procedure, the stability of the endemic equilibrium of the model with m ≠ 0 can also be verified by checking the signs of the eigenvalues at the endemic equilibrium point. However, due to the complexity and lengthiness of the calculations involved in analyzing the stability of the endemic equilibrium, the analytical calculations have been omitted in this study.

4.2 Basic reproduction number

The first half of the model represents the disease transmission for Delhi, while the second half, combined with the first, represents transmission for UP. To find R 0 Delhi , we consider the initial seven equations of the proposed model system. In view of the aggressive movement of migrants during the first lock-down period, we consider the Delhi patch with non-zero migration parameter. The disease-free equilibrium ( E 1 ) of model (1) for Delhi patch can be considered in a partitioned form as E 1 Delhi = ( b d + m , 0 , 0 , 0 , 0 , 0 , 0 ) . In the proposed model, progression from E to I A or I S is not considered a new infection but rather a progression through various compartments. Following the next-gen matrix technique [28], the transmission and transition matrices for Delhi can be given by

ℱ = 0 β 1 S 1 ( I A 1 + I S 1 + η I I 1 ) a 1 N 1 0 0 0 0 0 V = − b 1 + β 1 S 1 ( I A 1 + I S 1 + η I I 1 ) a 1 N 1 + d 1 S 1 − m S 1 ( ε + k 1 ) E 1 + m E 1 − k 1 E 1 − α 1 I A 1 + χ Q 1 − k 1 ε E 1 + ( α 1 + χ ) I A 1 + m I A 1 − ( 1 − k 1 ) ε E 1 + φ 1 I S 1 + m I S 1 − φ 1 I S 1 + ( χ + δ 1 ) I I 1 − χ 1 ( I A 1 + I I 1 + Q 1 ) + d 1 R 1 + m R 1 .

The infected compartments are E 1 , Q 1 , I A 1 , I S 1 , and I I 1 , which gives us the order m = 5 . Without loss of generality, it can be assumed that at the disease-free state S 0 = N . This implies

F = 0 0 β 1 a 1 β 1 a 1 β 1 η a 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V = ε + k 1 + m 0 0 0 0 − k 1 χ − α 1 0 0 − k 1 ε 0 α 1 + χ + m 0 − ( 1 − k 1 ) ε 0 0 φ 1 + m 0 0 0 0 − φ 1 χ + δ 1 .

This gives

V − 1 = 1 k 1 + m + ε 0 0 0 0 k 1 ε α 1 + k 1 ( m + α 1 + χ ) ( k 1 + m + ε ) χ ( m + α 1 + χ ) 1 χ α 1 χ ( m + α 1 + χ ) 0 0 k 1 ε ( k 1 + m + ε ) ( m + α 1 + χ ) 0 1 m + α 1 + χ 0 0 − ( − 1 + k 1 ) ε ( k 1 + m + ε ) ( m + φ 1 ) 0 0 1 m + φ 1 0 − ( − 1 + k 1 ) ε φ 1 ( k 1 + m + ε ) ( m + φ 1 ) ( δ 1 + χ ) 0 0 φ 1 ( m + φ 1 ) ( δ 1 + χ ) 1 δ 1 + χ .

Now, we know that R 0 = ρ ( F V − 1 ) , where ρ ( M ) represents the spectral radius of a matrix A.

Thus,

(27) R 0 Delhi = β 1 ε a 1 ( k 1 + m + ε ) k 1 m + α 1 + χ + ( k 1 − 1 ) ( m + φ 1 ) ( η φ 1 − δ 1 − χ ) ( δ 1 + χ ) .

Using the parameter values provided in Table 3, R 0 Delhi is estimated to be 2.5078.

Next, we consider the whole compartmental model to derive the expression of R 0 UP . New infections appear in precisely three compartments: exposed ( E 2 ), infected asymptomatic ( I A 2 ), and infected symptomatic ( I S 2 ). In a disease-free state, the susceptible population of Uttar Pradesh is considered equal to the state’s total population. Hence, by using the next-generation matrix technique, we have

ℱ = 0 0 0 β 2 S 2 ( I A 2 + I S 2 + η I I 2 ) a 2 N 2 + m E 1 0 0 0 m I A 1 0 m I S 1 0 0 0 0 V = − b 1 + d 1 S 1 + m S 1 + β 1 S 1 ( I A 1 + I S 1 + η I I 1 ) a 1 N 1 − b 2 + d 2 S 2 + m S 2 + β 2 S 2 ( I A 2 + I S 2 + η I I 2 ) a 2 N 2 − β 1 S 1 ( I A 1 + I S 1 + η I I 1 ) a 1 N 1 + ( ε + k 1 ) E 1 + m E 1 ( ε + k 2 ) E 2 − k 1 E 1 − α 1 I A 1 + χ Q 1 − k 2 E 2 − α 2 I A 2 + χ Q 2 − k 1 ε E 1 + ( α 1 + χ ) I A 1 + m I A 1 − k 2 ε E 2 + ( α 2 + χ ) I A 2 − ( 1 − k 1 ) ε E 1 + φ 1 I S 1 + m I S 1 − ( 1 − k 2 ) ε E 2 + φ 2 E 2 − φ 1 I S 1 + ( χ + δ 1 ) I I 1 − φ 2 I S 2 + ( χ + δ 2 ) I I 2 − χ ( I A 1 + I I 1 + Q 1 ) + d 1 R 1 + m R 1 − χ ( I A 2 + I I 2 + Q 2 ) + d 2 R 1 − m R 1 .

Further at the disease-free state, the transmission part F and transition part V can be given as follows:

F = 0 0 0 0 0 0 0 0 0 0 m 0 0 0 0 β 2 a 2 0 β 2 a 2 0 η β 2 a 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .

V = k 1 + m + ε 0 0 0 − β 1 a 1 0 − β 1 a 1 0 − η β 1 a 1 0 0 ε + k 2 0 0 0 0 0 0 0 0 − k 1 0 χ 0 − α 1 0 0 0 0 0 0 − k 2 0 χ 0 − α 2 0 0 0 0 − ε k 1 0 0 0 m + χ + α 1 0 0 0 0 0 0 − ε k 2 0 0 0 χ + α 2 0 0 0 0 ( − 1 + k 1 ) ε 0 0 0 0 0 m + φ 1 0 0 0 0 ε ( − 1 + k 2 ) 0 0 0 0 0 φ 2 0 0 0 0 0 0 0 0 − φ 1 0 χ + δ 1 0 0 0 0 0 0 0 0 − φ 2 0 χ + δ 2 .

To avoid the complications in the evaluation, we have substituted Γ = m + φ 1 , ϒ = m + χ + α 1 , Φ = k + m + ε , Ω = χ + δ 1 , and k 1 0 = k 1 − 1 . Then V − 1 can be partitioned as follows:

V 1 = ϒ Ω Γ ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 0 0 β 1 Ω Γ ε k 1 0 ϒ η β 1 φ 1 + Ω ( ε k 1 0 ϒ β 1 + ( ϕ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 1 ε + k 2 0 0 0 k 1 ( m + χ + ( 1 + ε ) α 1 ) Ω Γ χ ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 1 χ 0 ε k 1 0 α 1 η β 1 φ 1 − Ω ( − ε k 1 0 α 1 β 1 − ( Φ a 1 α 1 + k 1 β 1 ) Γ ) ε χ k 1 0 ϒ η β 1 φ 1 + χ Ω ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 k 2 ( χ + ( 1 + ε ) α 2 ) χ ( ε + k 2 ) ( χ + α 2 ) 0 1 χ 0 ε k 1 Ω Γ ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 0 0 ε k 1 0 η β 1 φ 1 a 1 + Ω ε k 1 0 β 1 a 1 + Φ Γ ε k 1 0 ϒ η β 1 φ 1 a 1 + ( χ + δ 1 ) ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 ε k 2 ( ε + k 2 ) ( χ + α 2 ) 0 0 0 − ε k 1 0 ϒ Ω ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 0 0 − ε k 1 0 β 1 Ω ε k 1 0 ( m + χ η β 1 φ 1 + Ω ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 − ε ( − 1 + k 2 ) ( ε + k 2 ) φ 2 0 0 0 − ε k 1 0 ϒ φ 1 ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 0 0 − ε k 1 0 β 1 φ 1 ε k 1 0 ϒ η β 1 φ 1 + Ω ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 − ε ( − 1 + k 2 ) ( ε + k 2 ) ( χ + δ 2 ) 0 0 0

V 2 = 0 ϒ ( β 1 Ω + η β 1 φ 1 ) ε k 1 0 ϒ η β 1 φ 1 + Ω ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 ϒ η β 1 Γ ε k 1 0 ϒ η β 1 φ 1 + Ω ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 0 0 0 0 0 0 k 1 ( m + χ + ( 1 + ε ) α 1 ) ( β 1 ( χ + δ 1 ) + η β 1 φ 1 ) ε χ k 1 0 ϒ η β 1 φ 1 + χ ( χ + δ 1 ) ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 k 1 ( m + χ + ( 1 + ε ) α 1 ) η β 1 Γ ε χ k 1 0 ϒ η β 1 φ 1 + χ ( χ + δ 1 ) ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 α 2 χ 2 + χ α 2 0 0 0 0 0 ε k 1 ( β 1 ( χ + δ 1 ) + η β 1 φ 1 ) ε k 1 0 ϒ η β 1 φ 1 + ( χ + δ 1 ) ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 ε k 1 η β 1 Γ ε k 1 0 ϒ η β 1 φ 1 + Ω ( ε k 1 0 υ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 1 χ + α 2 0 0 0 0 0 Φ ϒ − ε k 1 β 1 a 1 Ω ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 − ε k 1 0 ϒ η β 1 ε k 1 0 ϒ η β 1 φ 1 + Ω ( ε k 1 0 ϒ β 1 + ( Φ a 1 ϒ − ε k 1 β 1 ) Γ ) 0 0 0 1 φ 2 0 0 0 Φ ϒ − ε k 1 β 1 a 1 φ 1 ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ 0 ε k 1 0 ϒ β 1 a 1 + Φ ϒ − ε k 1 β 1 a 1 Γ ε k 1 0 ϒ η β 1 φ 1 a 1 + Ω ε k 1 0 ϒ β 1 a 1 + Φ ϒ ϒ − ε k 1 β 1 a 1 Γ 0 0 0 1 χ + δ 2 0 1 χ + δ 2

Now,

R 0 = ρ ( F V − 1 ) ,

here, ρ is the slandered notation for spectral radius. Further calculation leads us to

(28) R 0 UP = ε β 2 ( ε + k 2 ) a 2 k 2 χ + α 2 + ( k 2 − 1 ) ( η φ 2 − χ − δ 2 ) φ 2 ( χ + δ 2 ) ,

By using the estimated parameter values provided in Table 1 and the aforementioned expression, R 0 UP is estimated to be 1.5218.

Remark

The higher R 0 for Delhi indicates a more severe outbreak in Delhi at the beginning. The reason might be the influx of international migrants in Delhi from COVID-19-infected countries before March 2020. However, the disease prevalence in UP was relatively low in the beginning due to the significantly low influx of international migrants.

5 Doubling time during initial phases of lockdown

As per the previous section, both regions have a reproduction number greater than one, ensuring exponential growth of diseases [19]. Consequently, it is evident that the number of cases will double in both regions after a certain period of time. Therefore, we estimate the doubling time for the exponential growth of disease in the Delhi–UP region using the data of positive cases for the first two phases of lockdown [19]. Following the methodology explained in [1], the SAS DATA step estimates the doubling time by using the slope estimates at the end of the first lockdown (21 days) and the second lockdown (19 days), assuming that the number of cases continues to grow at the estimated rate on the last day. The doubling time and epidemic strength are inversely related, so if an epidemic declines, the doubling time increases.

In Figure 4, the count for positive cases in Delhi is predicted to double rapidly as the value of M is higher for the region, whereas the count for UP is predicted to double slowly. The table corresponding to Figure 4 indicates the doubling time for Delhi at 3.980 ≈ 4 days and for UP at 4.647 ≈ 5 days. According to the State Census 2011 [2], Delhi is the nation’s most densely populated state with a density of 11,320 per square kilometer, which might have significantly contributed to the COVID-19 surge with a shorter doubling time.

Figure 4

Cases during first phase of lockdown in Delhi–UP region.

We next examine the effectiveness of the second phase of complete lockdown by estimating the doubling time during April 15, 2020 May 3, 2020, in the Delhi–UP region. The table in Figure 5 indicates that Delhi’s doubling time increased significantly due to the 19-day second lockdown, during which mobility restrictions on the urban population had slowed down the spread ( M = 0.037 ). It is noteworthy that the doubling time in UP and Delhi were equal, irrespective of Delhi’s high population density.

Figure 5

Cases during second phase of lockdown in Delhi–UP region.

Even after two phases of lockdown, the cases in both regions continued to increase due to reverse migration. This situation highlights that despite being the most populous state with a large rural division and facing challenges such as lack of resources and reverse migration of workers, the UP government effectively managed the case load. However, there seems to be a discrepancy between the estimated doubling time for UP and the actual prevalence of the disease, as mentioned in Section 5.1. Several factors, as outlined in the mentioned Section 5.1, could be responsible for this mismatch.

5.1 Factors for mismatch between real scenario and collected data

Considering that Delhi is the capital city and has a smaller area, the available data can be considered relatively more reliable. On the other hand, UP, with a population of 19.96 crores, faced challenges in accurately consolidating data on active cases, particularly in its rural division.
As a union territory, Delhi is characterized by a high concentration of urban areas and a high literacy rate. This has contributed to a higher awareness among the urban population, which in turn facilitated the consolidation of data on COVID-19 cases. However, in UP, there has been a lack of awareness about COVID-19 among rural populations, which, coupled with reluctance to undergo testing and treatment, has hindered the proper collection of positive cases data.
The Arogya Setu application, which is used for contact tracing, was found to be effective in Delhi, largely due to its usage among the privileged urban class. However, in UP, the majority of migrants were daily wage workers, with limited access to smartphones. As a result, the Arogya Setu application was not as effective in tracking down contacts among the rural and underprivileged populations of UP.

6 Migration-induced unintended synchrony

Previous studies, such as those by Borrego-Salcido et al. [5], Earn et al. [10], Qian et al. [20], and Zhang et al. [30], suggest that infected individuals can spread epidemics from one community to another, resulting in long-term synchrony. In the case of individuals passing from Delhi to UP, they create interlinks between the two states. To investigate the synchrony of epidemic trajectories between the Delhi–UP region, we used the methodology proposed by Pikovsky et al. [18]. We explored the change in mobility level during late March 2020 and considered both symptomatic and asymptomatic populations from both regions. We plotted the phase synchronization coefficient (PSC) as the absolute value of the mean exp ( j Δ θ ) , where a higher value of PSC indicates anti-synchronization, and a lower value quantifies phase synchronization. By using Hilbert transformation, we plotted PSC as a function of the migration rate for two sets of epidemic compartments: ( I A 1 , I A 2 ) , and ( I S 1 , I S 2 ) . After calculating the phases of the time series, we characterize the synchronization that contributes to more complex forms, such as anti-phase synchrony between the coupled regions.

In Figure 6(a), it can be observed that the phase synchronization coefficient (PSC) for trajectories of asymptomatic infected individuals is relatively high even for weak levels of migration rate, but it starts declining at high levels. This suggests that as the migration rate increased, asymptomatic individuals who slipped through screening and reached UP contributed to phase synchrony between the regions. Since Delhi had a higher R 0 -value with a significantly shorter doubling time, the Government of India (GoI) should have responded more effectively to stop the migration and control the synchrony between I A 1 and I A 2 .

$Figure 6 Phase synchronization coefficient versus migration rates at various levels { 0 , 0.002 , 0.004 , 0.006 , 0.008 } \left\{0,0.002,0.004,0.006,0.008\right\} , { 0.01 , 0.02 , 0.04 , 0.06 , 0.08 , 0.1 , 0.2 , 0.4 , 0.6 , 0.8 } \left\{0.01,0.02,0.04,0.06,0.08,0.1,0.2,0.4,0.6,0.8\right\} (a) ( I A 1 , I A 2 ) \left({I}_{{A}_{1}},{I}_{{A}_{2}}) and (b) ( I S 1 , I S 2 ) \left({I}_{{S}_{1}},{I}_{{S}_{2}}) .$

Figure 6

Phase synchronization coefficient versus migration rates at various levels { 0 , 0.002 , 0.004 , 0.006 , 0.008 } , { 0.01 , 0.02 , 0.04 , 0.06 , 0.08 , 0.1 , 0.2 , 0.4 , 0.6 , 0.8 } (a) ( I A 1 , I A 2 ) and (b) ( I S 1 , I S 2 ) .

In Figure 6(b), a critical point is displayed approximately at m = 0.2 , below which simulations characterize phase synchronization ( P S C < 0.525 ), and above which they are mainly anti-phase due to a higher PSC. Due to a higher level of migration rate ( m ), symptomatic individuals held a significantly higher fraction at one location than at the other. However, the predicted scenario highlights the burden of COVID-19 in UP due to the influx of migrants who escaped from Delhi and introduced the disease to rural parts of UP. This resulted in community transmission due to a lack of awareness about the disease, subsequently increasing the case load in UP compared to Delhi.

6.1 Effective plan of action

The lack of coordination between the governments of Delhi and Uttar Pradesh in arranging transportation for migrant workers resulted in a significant failure, leading to mass gatherings and chaos at the Delhi bus stand. This setback undermined the effectiveness of non-pharmaceutical interventions (NPIs) in controlling the spread of COVID-19. The consequences could have been mitigated if the central government had intervened with effective strategies and precautions to restrict movement across Delhi. Such strategies could have included raising awareness about the disease, improving healthcare infrastructure, enhancing testing capacity, ensuring adequate food supply, and setting up large-capacity centralized quarantine centers in Delhi. Instead of hastily arranging bus transportation, migrants could have been directed to stay at designated quarantine centers in Delhi for 14 days during the incubation period. Through testing, tracking, and tracing, asymptomatic migrants could have been released from the quarantine centers in Delhi and transported back to Uttar Pradesh via properly arranged means of transportation.

7 Simulation and parameter estimation

7.1 Migration and epidemiological factors shaping disease dynamics

In Section 4.2, we derived the expression for the basic reproduction number and calculated its exact value for both states using estimated parameter values. Here, we will investigate the combined impact of migration and other parameters on the basic reproduction number through numerical simulations.

Based on the analysis of Figure 7(a), it is evident that an increase in the migration coefficient, coupled with the contact rate, directly amplifies the spread of the disease by increasing the basic reproduction number in Delhi. However, as shown in Figure 7(b), when the rate of quarantine for exposed individuals is increased along with the migration rate, the spread of the disease is mitigated compared to the previous case. Further, Figure 8(a) displays the simultaneous impact of rate of detection through contact tracing in Delhi patch and contact rate on basic reproduction number whereas Figure 8(b) displays the simultaneous impact of rate of quarantine rates of exposed individuals in Delhi patch and contact rate on basic reproduction number. Figure 9 shows the simultaneous impact of the contact rates of Delhi and UP on the overall basic reproduction number. It can be observed that as the values of contact rates β 1 and β 2 increase, the value of R 0 also increases. Our simulation results provide theoretical evidence that to curb the COVID-19 pandemic, the most feasible approach is to control the mixing of susceptible populations by reducing contact as much as possible.

$Figure 7 Simultaneous impact of (a) migration rate m m and contact rate β 1 {\beta }_{1} on the basic reproduction number and (b) migration rate m m and k 1 {k}_{1} on the basic reproduction number.$

Figure 7

Simultaneous impact of (a) migration rate m and contact rate β 1 on the basic reproduction number and (b) migration rate m and k 1 on the basic reproduction number.

$Figure 8 Simultaneous impact of (a) α 2 {\alpha }_{2} and contact rate β 2 {\beta }_{2} on the basic reproduction number and (b) β 2 {\beta }_{2} and k 2 {k}_{2} on the basic reproduction number.$

Figure 8

Simultaneous impact of (a) α 2 and contact rate β 2 on the basic reproduction number and (b) β 2 and k 2 on the basic reproduction number.

$Figure 9 Simultaneous impact of contact rates β 1 {\beta }_{1} and β 2 {\beta }_{2} on the basic reproduction number.$

Figure 9

Simultaneous impact of contact rates β 1 and β 2 on the basic reproduction number.

7.2 Global sensitivity analysis of basic reproduction numbers

Section 7.1 highlights the dependency of R 0 on two specific parameters, although the model (1) actually depends on a total of 21 parameters. To identify the crucial parameters affecting R 0 , we conducted a global sensitivity analysis using Latin hypercube sampling (LHS) - partial rank correlation coefficient (PRCC) approach, which allows us to explore the entire parameter space of the model. LHS evenly samples the parameter space, and PRCC identifies the most significant and sensitive parameters in terms of their impact on R 0 Delhi and R 0 UP . The PRCC values, ranging from − 1 to 1, quantify the variability in R 0 , with positive values indicating that as a parameter increases, so does R 0 , and negative values indicating the opposite. PRCC values of large magnitude ( 0.5 ≤ ∣ PRCC ∣ ≤ 1 ) characterize the most significant parameters [16]. Since the response functions are given by equations (27) and (28), we assumed a uniform distribution for each parameter in model (1). We conducted 300 simulations per LHS, using a 10 % deviation from the baseline values of the parameters mentioned in (3).

Figure 10 identifies five parameters that influence R 0 Delhi and R 0 UP the most, on which feasible control could have curbed the virus spread.

χ : rate of recovery,
η : efficacy of isolation to home-quarantined individuals,
ϕ i : rate of quarantine of symptomatic individuals, ∀ i ∈ { 1 , 2 } ,
β 2 : effective contact rate for UP.

The five parameters mentioned earlier, namely, χ = 1 recovery period , η , and ϕ 1 , 2 , exhibit a common positive correlation with the values of R 0 for both UP and Delhi, as shown in Figure 10. The dependency of R 0 on χ highlights that the longer an infected person remains sick, the higher the risk of transmission, resulting in larger values of R 0 for UP and Delhi. The significance of η is evident in indicating that effective home quarantine measures could nullify the exposure of infected individuals to susceptible. Given the high population density in these states, home quarantine was a sustainable option that balanced individual freedom and public protection while conserving resources. However, initially, home quarantine was less effective than centralized quarantine due to low awareness about the disease among the population [14]. The parameters ϕ 1 , 2 show a negative correlation coefficient with R 0 for both states, suggesting that a higher quarantine rate for symptomatic individuals could have effectively curbed the infection. Finally, the positive PRCC value of β 2 indicates a correlation with disease prevalence in UP, implying that a higher fraction of capable contacts leading to infection would result in an increased estimated value of R 0 . Overall, this analysis emphasizes the importance of improved treatment, better quarantine facilities, and enhanced detection methods in curbing the spread of the disease. Therefore, it is essential to carefully consider these five parameters while estimating R 0 from data [14].

Figure 10

Sensitivity of the basic reproduction numbers for Delhi and UP to changes in the model parameters using PRCC index.

7.3 Parameter estimation using case reports

We estimate the parameters by fitting the model outputs to the observed counts of COVID-19 cases in the Delhi and UP states. We used the case data from March 1 to May 31, 2020, which includes a countrywide lockdown, and there was a migration of labours from Delhi to the state of UP. The lockdown occurred from March 25 to May 31, carried out in four phases: March 25 to April 14, April 15 to May 3, May 4 to May 17, and finally, May 18 to May 31, when the national lockdown was lifted. No further dates are considered because both states significantly differ in their lockdown policies post-May 31. Moreover, since migration only began after lockdown, we consider no migration ( m = 0 ) before March 25 to fit the model to the data.

The following are the parameters we have estimated:

The reporting probability, p , was different for the two states. We calibrate this using results from Unnikrishnan et al. [27], which reflects the underreporting of COVID-19 cases in Indian states. The fraction of cases reported is estimated from different case fatality ratios (CFRs), assuming that the number of disease-induced deaths is correctly reported. We take the case where the CFR is 0.66%. For this assumption, the percentage of cases reported for Delhi and UP are 31 and 39%, respectively. As depicted in the model fit Figure 11, the model nicely captures the exponential growth pattern shown in the case data. We measure the goodness-of-fit R 2 between the actual data and model output. The R 2 values for Delhi output and UP output are 0.857 and 0.881, respectively, which means that the model can explain 85.7 and 88.1% of the data. These values indicate that the model fits the actual data reasonably well. The estimated values of the parameters are given in Table 4.

Figure 11

The model fits the daily confirmed COVID-19 cases in Delhi and Uttar Pradesh. The solid lines represent the model simulated predictions, while the dots represent the actual number of daily confirmed cases. The initial conditions used for the model are S 1 ( 0 ) = ( a 1 − 37 ) , E 1 ( 0 ) = 0 , Q 1 ( 0 ) = 0 , I A 1 ( 0 ) = 35 , I S 1 ( 0 ) = 2 , I I 1 ( 0 ) = 0 , R 1 ( 0 ) = 0 , S 2 ( 0 ) = ( a 2 − 42 ) , E 2 ( 0 ) = 0 , Q 2 ( 0 ) = 0 , I A 2 ( 0 ) = 35 , I S 2 ( 0 ) = 7 , and I I 2 = 0 , R 2 = 0 . The initial values of the number of symptomatic infected were taken from actual data while the numbers of asymptomatic infected were estimated since there are no accurate data for asymptomatic cases, and given their low detectability, the actual number of cases would be much higher.

Table 4

Estimated parameter values

Parameter	Description	Estimated values
β 1 , 2	Contact rate	0.191, 0.158
k 1 , 2	Rate of quarantine of exposed individuals	0.059, 0.09
α 1 , 2	Rate of detection through contact tracing	0.069, 0.086
ϕ 1 , 2	Rate of isolation of symptomatic individuals	0.122, 0.111
δ 1 , 2	Disease induced death rate	0.0026, 0.0025
Ω	Delay in traveling	4.5
m	Migration rate	0.0005

7.4 Scenario analysis

In this subsection, we conduct scenario analysis based on two different scenarios picked from the estimated parameters to judge migration’s impact and lockdown effectiveness for each state. Scenario (I) is the initiation of lockdown on different dates. Four different scenarios were taken to observe how these impact the number of cases: 20 days before official lockdown; 10 days before; 10 days after; and 20 days after, i.e., on March 5, March 15, April 1, and April 11, respectively. Our output shows that an earlier lockdown would have resulted in a lower COVID-19 caseload, and there is a higher caseload when the lockdown implementation is later than the actual official lockdown. Also, there is a more significant effect of changing lockdown dates in Delhi than in UP (Figure 12).

Figure 12

Scenario I analysis by taking different lockdown dates (a): for Delhi, (b): for UP.

Scenario (II) is the absence of migration, and it was applied to each of the different lockdown scenarios. For the dates on and after the official lockdown, the number of cases when there is no migration is lower for both the states, though we can see that there is a slightly more significant effect on UP when there is no influx of migrants. On the day of the official lockdown, the UP shows an 8.04% decrease without migration, and Delhi shows a 3.98% decrease. Similarly, the percentage decrease in cases for 10 days after lockdown and 20 days after lockdown is 7.01 and 5.51% for UP and 5.15 and 4.68% for Delhi. Ten days before the lockdown, the number of cases in Delhi without migration was slightly higher. The retention of migrants at this particular date might have caused the increase in cases by 5.99%, while UP decreased by 6.42%. The best results can be seen in the 20 days before lockdown. Delhi and UP case loads go down by 18.43 and 18.37%, respectively, which is the most significant decrease in infections among all the cases considered. Therefore, having an earlier lockdown has the most advantageous effect in curbing the virus (Figure 13).

Figure 13

Scenario analysis by taking different lockdown date (a) 20 days before lockdown, (b) 10 days before lockdown, (c) 10 days after lockdown, (d) 20 days after lockdown.

8 Discussion and conclusion

The emergence of the COVID-19 pandemic posed numerous challenges for governments worldwide due to the incomplete understanding of the clinical presentation of the disease. In India, despite implementing a stringent lockdown soon after the outbreak was declared a pandemic, there were failures in effectively slowing down and curbing the spread of the virus, particularly due to the unknown symptoms and the unplanned return of migrant workers from Delhi to their hometowns in Uttar Pradesh (UP) due to their employment and economic situations.

To analyze the impact of reverse migration on the spread of COVID-19, we developed a deterministic model based on the migration between Delhi and UP during the first lockdown. We mathematically demonstrated the well-posedness and existence of equilibrium points in the model. The analysis of COVID-19 data revealed short doubling times for exponential growth in Delhi, despite the complete closure of 21 days during the lockdown. The synchrony effects confirmed the potential transmission risk to the UP region due to population outflow from Delhi. Immediate control of migrant flow from Delhi was found to be necessary to curb the spread within the community, as highlighted by the synchrony between infected individuals from both regions. Remarkably, the analytical results for the transmission potential of the disease in both states were in complete agreement with the actual scenario. In addition, the sensitivity analysis using partial rank correlation coefficient (PRCC) identified that the effectiveness of home quarantine and recovery rate could have reduced the R 0 value.

For parameter estimation, we focused on key parameters, such as the contact rate ( β ) and migration rate ( m ) , which were used to model the lockdown and movement of migrants, respectively, emulating the real-life situation. Maximum likelihood estimation (MLE) was used assuming a Poisson distribution for the number of daily cases, considering the underreporting of cases using a reporting probability. Furthermore, a stochastic element was introduced in the model by modeling it as a Poisson process. Our numerical simulations aligned well with the analytical results, supporting the conclusion that reducing contact among the population could have controlled the spread of the virus.

To conduct scenario analysis, we considered two cases: scenario (I) involved setting different dates for the lockdown and scenario (II) involved the absence of migration at each of these different dates. These scenarios were chosen to assess the effectiveness of an earlier lockdown in slowing down the number of cases and to evaluate the role of migration in accumulating COVID-19 cases. Different values of β were used to represent social behaviors such as gathering, mask-wearing, and social distancing, which could vary before and after the lockdown. The findings indicated that implementing an earlier lockdown resulted in fewer cases, suggesting that timely intervention could have been more effective in slowing the spread of the virus. Furthermore, the cases were lower in the absence of migration.

During the course of this study, we recognized the limitations of data collection, particularly in matching theoretical studies with real-time situations. With an estimated 83.33% rural population in UP, data consolidation of cases was limited to the urban division by late March, which posed challenges in designing short-term and long-term welfare measures to address future outbreaks [11].

Funding information: Research of Nitu Kumari is funded by Council of Scientific & Industrial Research (CSIR) under the grant 25[22832]/2022. Research of Sumit Kumar is funded by University Grants Commission (UGC) with file 1167/(CSIRNETJUNE2019).
Author contributions: All authors designed, analyzed, and prepared the presented work.
Conflict of interest: The authors state no conflict of interest.
Ethical approval: This research did not required ethical approval.

References

[1] blogs.sas.com. https://blogs.sas.com/content/iml/2020/04/01/estimate-doubling-time-exponential-growth.html. Accessed: April 2020. Search in Google Scholar

[2] census2011.co.in. https://www.census2011.co.in/states.php. Accessed: 2011. Search in Google Scholar

[3] Alqarni, M. S., Alghamdi, M., Muhammad, T., Alshomrani, A. S., & Khan, M. A. (2022). Mathematical modeling for novel coronavirus (covid-19) and control. Numerical Methods for Partial Differential Equations, 38(4), 760–776. 10.1002/num.22695Search in Google Scholar PubMed PubMed Central

[4] Biswas, S. K., Ghosh, J. K., Sarkar, S., & Ghosh, U. (2020). Covid-19 pandemic in india: A mathematical model study. Nonlinear Dynamics, 102(1), 537–553. 10.1007/s11071-020-05958-zSearch in Google Scholar PubMed PubMed Central

[5] Borrego-Salcido, C., Juárez-del Toro, R., & Fonseca-Zendejas, A. (2023). The waves and cycles of covid-19 pandemic: A phase synchronization approach. AJS, 52, 25–38. 10.17713/ajs.v52i3.1450Search in Google Scholar

[6] Chatterjee, K., Chatterjee, K., Kumar, A., & Shankar, S. (2020). Healthcare impact of covid-19 epidemic in India: A stochastic mathematical model. Medical Journal Armed Forces India, 76(2), 147–155. 10.1016/j.mjafi.2020.03.022Search in Google Scholar PubMed PubMed Central

[7] Kerala Defeats Coronavirus. India’s three covid-19 patients successfully recover. The Weather Channel. Archived from the original on, 18, 2020. Search in Google Scholar

[8] DeJesus, E. X., & Kaufman, C. (1987). Routh-hurwitz criterion in the examination of eigenvalues of a system of nonlinear ordinary differential equations. Physical Review A, 35(12), 5288. 10.1103/PhysRevA.35.5288Search in Google Scholar

[9] Duddu, P. (2020). Coronavirus outbreak: Safety measures at major international airports. Search in Google Scholar

[10] Earn, D. J., Rohani, P., & Grenfell, B. T. (1998). Persistence, chaos and synchrony in ecology and epidemiology. Proceedings of the Royal Society of London. Series B: Biological Sciences, 265(1390), 7–10. 10.1098/rspb.1998.0256Search in Google Scholar PubMed PubMed Central

[11] Kumar, A., Nayar, K. R., & Koya, S. F. (2020). Covid-19: Challenges and its consequences for rural health care in India. Public Health in Practice, 1, 100009. 10.1016/j.puhip.2020.100009Search in Google Scholar PubMed PubMed Central

[12] Laohombé, A., NgningoneEya, I., Tewa, J. J., Bah, A., Bowong, S., & Oukouomi Noutchie, S. C. (2014). Mathematical analysis of a general two-patch model of tuberculosis disease with lost sight individuals. Abstract and Applied Analysis, 2014, 1–14.10.1155/2014/263780Search in Google Scholar

[13] Lusome, R., & Bhagat, R. Internal Migration in India 1971–2011. Search in Google Scholar

[14] Ngonghala, C. N., Goel, P., Kutor, D., & Bhattacharyya, S. (2021). Human choice to self-isolate in the face of the covid-19 pandemic: A game dynamic modelling approach. Journal of Theoretical Biology, 521, 110692. 10.1016/j.jtbi.2021.110692Search in Google Scholar PubMed PubMed Central

[15] Government of India Ministry of Health and Family Welfare Directorate General of Health Services. Guidelines for home quarantine. https://www.mohfw.gov.in/pdfGuidelinesforhomequarantine.pdf. Search in Google Scholar

[16] Dr. Denise Kirschner University of Michigan. Uncertainty and sensitivity analysis. http://malthus.micro.med.umich.edu/lab/usadata/. Search in Google Scholar

[17] World Health Organization, et al. (2020). Coronavirus Disease 2019 (COVID-19): Situation Report, 88. Search in Google Scholar

[18] Pikovsky, A., Rosenblum, M., & Kurths, J. (2000). Phase synchronization in regular and chaotic systems. International Journal of Bifurcation and Chaos, 10(10), 2291–2305. 10.1142/S0218127400001481Search in Google Scholar

[19] prsindia. Details on cases. Search in Google Scholar

[20] Qian, X., & Ukkusuri, S. V. (2021). Connecting urban transportation systems with the spread of infectious diseases: A trans-seir modeling approach. Transportation Research Part B: Methodological, 145, 185–211. 10.1016/j.trb.2021.01.008Search in Google Scholar

[21] Rauta, A. K., Rao, Y. S., & Behera, J. (2021). Spread of covid-19 in Odisha (India) due to influx of migrants and stability analysis using mathematical modeling. Artificial Intelligence for COVID-19 (pp. 295–309). 10.1007/978-3-030-69744-0_17Search in Google Scholar

[22] Sahoo, P. K., Biswal, S., Kumar, H., & Powell, M. (2022). Urban to rural covid-19 progression in india: The role of massive migration and the challenge to india’s traditional labour force policies. The International Journal of Health Planning and Management, 37(1), 528–535. 10.1002/hpm.3327Search in Google Scholar PubMed PubMed Central

[23] Sahu, G. P., & Dhar, J. (2015). Dynamics of an seqihrs epidemic model with media coverage, quarantine and isolation in a community with pre-existing immunity. Journal of Mathematical Analysis and Applications, 421(2), 1651–1672. 10.1016/j.jmaa.2014.08.019Search in Google Scholar PubMed PubMed Central

[24] Sharma, G. D., & Mahendru, M. (2020). Lives or livelihood: Insights from locked-down india due to covid19. Social Sciences & Humanities Open, 2(1), 100036. 10.1016/j.ssaho.2020.100036Search in Google Scholar PubMed PubMed Central

[25] Siddiqui, A. F., Wiederkehr, M., Rozanova, L., & Flahault, A. (2020). Situation of india in the covid-19 pandemic: India’s initial pandemic experience. International Journal of Environmental Research and Public Health, 17(23), 8994. 10.3390/ijerph17238994Search in Google Scholar PubMed PubMed Central

[26] Tewa, J. J., Bowong, S., & Mewoli, B. (2012). Mathematical analysis of two-patch model for the dynamical transmission of tuberculosis. Applied Mathematical Modelling, 36(6), 2466–2485. 10.1016/j.apm.2011.09.004Search in Google Scholar

[27] Unnikrishnan, J., Mangalathu, S., & Kutty, R. V. (2021). Estimating under-reporting of covid-19 cases in indian states: An approach using a delay-adjusted case fatality ratio. BMJ Open, 11, e042584. 10.1136/bmjopen-2020-042584. Search in Google Scholar PubMed PubMed Central

[28] Van den Driessche, P., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1–2), 29–48. 10.1016/S0025-5564(02)00108-6Search in Google Scholar

[29] Wiersinga, W. J., Rhodes, A., Cheng, A. C., Peacock, S. J., & Prescott, H. C. (2020). Pathophysiology, transmission, diagnosis, and treatment of coronavirus disease 2019 (covid-19): A review. Jama, 324(8), 782–793. 10.1001/jama.2020.12839Search in Google Scholar PubMed

[30] Zhang, B., Liang, S., Wang, G., Zhang, C., Chen, C., Zou, M., ..., Du, X. (2021). Synchronized nonpharmaceutical interventions for the control of covid-19. Nonlinear Dynamics, 106(2), 1477–1489. 10.1007/s11071-021-06505-0Search in Google Scholar PubMed PubMed Central

Received: 2023-03-20

Revised: 2023-05-07

Accepted: 2023-05-09

Published Online: 2023-07-19

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2022-0151

Keywords for this article

SARS-CoV-2; NPIs failure; migrated labours; doubling time; phase synchronization; parameter estimation; scenario analysis; PRCC analysis

Creative Commons

BY 4.0