key: cord-0908042-xb1a8n33
authors: Brozak, Samantha J.; Pant, Binod; Safdar, Salman; Gumel, Abba B.
title: Dynamics of COVID-19 pandemic in India and Pakistan: A metapopulation modelling approach
date: 2021-10-15
journal: Infect Dis Model
DOI: 10.1016/j.idm.2021.10.001
sha: 92e4e9052c63621d370042462e13cd117a62e3cc
doc_id: 908042
cord_uid: xb1a8n33

India has been the latest global epicenter for COVID-19, a novel coronavirus disease that emerged in China in late 2019. We present a base mathematical model for the transmission dynamics of COVID-19 in India and its neighbor, Pakistan. The base model was rigorously analyzed and parameterized using cumulative COVID-19 mortality data from each of the two countries. The model was used to assess the population-level impact of the control and mitigation strategies implemented in the two countries (notably non-pharmaceutical interventions). Numerical simulations of the basic model indicate that, based on the current baseline levels of the control and mitigation strategies implemented, the pandemic trajectory in India is on a downward trend. This downward trend will be reversed, and India will be recording mild outbreaks, if the control and mitigation strategies are relaxed from their current levels. By early September 2021, our simulations suggest that India could record up to 460,000 cumulative deaths under baseline levels of the implemented control strategies, while Pakistan (where the pandemic is comparatively milder) could see over 24,000 cumulative deaths at current mitigation levels. The basic model was extended to assess the impact of back-and-forth mobility between the two countries. Simulations of the resulting metapopulation model show that the burden of the COVID-19 pandemic in Pakistan increases with increasing values of the average time residents of India spend in Pakistan, with daily mortality in Pakistan peaking in mid-August to mid-September of 2021. Under the respective baseline control scenarios, our simulations show that the back-and-forth mobility between India and Pakistan could delay the time-to-elimination of the COVID-19 pandemic in India and Pakistan to November 2022 and July 2022, respectively.

yet infectious E(t)), pre-symptomatic infectious individuals (P (t)), symptomatically-infectious individuals (I(t)), asymptomatically-infectious individuals (A(t)), hospitalized individuals (H(t)), and recovered indi-128 viduals (R(t)), so that N (t) = S(t) + E(t) + P (t) + I(t) + A(t) + H(t) + R(t). The basic model is given by 129 the following deterministic system of nonlinear differential equations (where a dot represents differentiation 130 with respect to time t): individuals progress to the symptomatically-infectious class (at a rate rσ 2 ) and the remaining proportion, 140 1 − r, progress to the asymptomatically-infectious class (at a rate (1 − r)σ 2 ). Therefore, the intrinsic 141 incubation period of the disease is 1/σ 1 + 1/σ 2 . The parameter γ I (γ A )(γ H ) represents the recovery rate for Table 1 and Table 2 , respectively. Some of the main assumptions made in (c) the timescale of the novel pandemic is much shorter than the demographic timescale (so that births 156 and natural death processes can be ignored);

(d) natural recovery induces permanent immunity against future infection. Although there is currently 158 limited evidence regarding the duration of infection-acquired immunity, Dan et al.

[48] recently 159 estimated that natural immunity due to recovery from COVID-19 could last up to eight months. 160 For housekeeping purposes, it is convenient to let D(t) be the total number of individuals who died from 161 COVID-19 in the population. It follows from the basic model (2.1) that the equation for the rate of change 162 of D(t) isḊ = δ I I + δ H H. The basic qualitative properties of the single-patch model (2.1) will now be explored. Progression rate from exposed to pre-symptomatic class (latent period) r

Proportion of individuals who show clinical symptoms of the disease rσ 2

Progression rate from the pre-symptomatic to the symptomatic infectious class (1 − r)σ 2 Progression rate from the pre-symptomatic to the asymptomatic infectious class φ

Hospitalization rate for the symptomatically-infectious individuals γ I (γ A )(γ H )

Recovery rate for individuals in the I(A)(H) class δ I (δ H ) Disease-induced mortality rate for individuals in the I(H) compartment where N (0) is the initial total population. For the basic model (2.1) to be mathematically-and biologically-166 meaningful, it is necessary that the solutions of the basic model (2.1) remain non-negative for all non-167 negative initial conditions. That is, solutions that start in Ω remain in Ω for all time t > 0 (i.e., Ω is where N (0) is the initial total population, 0 < S * ≤ N (0), 0 ≤ R * < N (0), and 0 < S * + R * ≤ N (0). The asymptotic stability property of ξ 0 can be explored using the next generation operator method [49]. Using the notation in [49] , it follows that the associated non-negative matrix of new infection terms (F ) and the M-matrix of all linear transition terms (V ) are given, respectively, by

It is convenient to define the quantity (where ρ is the spectral radius):

where,

and N * = S * + R * (with S * and R * as defined above). The quantity R 0 is the basic reproduction number of the basic model (2.1). It measures the average number of new cases generated by a typical infected individual if introduced into a completely susceptible population (i.e., a population where no one has 179 immunity due to previous exposure to the disease and no public health intervention and mitigation measures 180 are implemented in the community). The quantities R P , R I , R A and R H are the constituent basic 181 reproduction numbers for the pre-symptomatic, symptomatic, asymptomatic, and hospitalized individuals, 182 respectively. The result below follows from Theorem 2 of [49] . 183 Theorem 2.2. The continuum of disease-free equilibria (ξ 0 ) of the basic model (2.1) is locally-asymptotically 184 stable (LAS) if R 0 < 1.

Theorem 2.2 implies that a small influx of COVID-19 cases will not cause a significant outbreak in the 186 community if R 0 < 1. In the case of epidemic models, such as the basic model (2.1), the epidemiological 187 requirement of having R 0 < 1, while sufficient, is not necessary for the elimination of the disease from 188 the community. This is owing to the fact that Kermack-McKendrick-type epidemic models (such as (2.1)) 189 do not include vital dynamics (i.e., no birth and natural death processes) and no immigration (since Remark. In the case that R 0 > 1, the epidemic will grow to a peak and then eventually decline to zero The quantity R P in the expression for R 0 is the product of the infection rate by pre-symptomatically 205 infectious individuals near the continuum of disease-free equilibria, given by β P ( S * N * ), and the average 206 duration in the pre-symptomatically infectious class ( 1 σ 2 ). Similarly, R I is the product of the infection rate 207 by symptomatically-infectious individuals near the disease-free equilibria, given by β I ( S * N * ), the probability 208 of surviving the pre-symptomatic class and moving to the symptomatic compartment ( rσ 2 σ 2 = r) and the 209 average duration in the symptomatically-infectious class ( 1 φ+γ I +δ I ). The quantity R A is the product of the 210 infection rate by asymptomatically-infectious individuals near the disease-free equilibria, given by β A ( S * N * ), The sum of R P , R I , R A and R H gives R 0 . When referring to the reproduction number R 0 obtained via the basic model (2.1) in India and Pakistan, we use the notation R 0i and R 0p , respectively.

In order to ensure that elimination of the pandemic when R 0 is less than unity is independent of the 220 initial size of the sub-populations of the basic model (2.1), it is necessary to show that the continuum of 221 disease-free equilibria (ξ 0 ) is globally-asymptotically stable. We claim the following result: Pakistan, respectively. Using the parameter values listed in Table 3 and . Subscript i denotes a parameter estimated using the cumulative mortality data for India, and subscript p denotes a parameter estimated using the cumulative mortality data for Pakistan. lockdown relaxation and strengthening scenarios discussed above, are further summarized in Table 5a . Table 6 . Time-to-elimination is defined as the date when new cases of COVID-19 are asymptotically Table 3 . Values of the estimated parameters prior to June 6, 2021 are listed in Table 4b . Projected cumulative mortality are rounded to the nearest hundred. To achieve the above objective, the basic model (2.1) will be extended to account for the back-and-forth mobility between the two countries. We assume that the mobility pattern between the two countries is largely Lagrangian in nature (i.e., it is for a limited period of time, and not permanent migration) [46] . In order to extend the basic model (2.1) to account for the Lagrangian mobility, we consider each of the two countries to be a single patch. Let the subscripts i and p represent the residents in the patch for India and Pakistan, respectively. Furthermore, for each patch k, we stratify the total population at time t, denoted by N k (t), into the mutually-exclusive compartments of susceptible (S k ), exposed/latent (E k ), pre-symptomatic (P k ), symptomatic (I k ), asymptomatic (A k ), hospitalized (H k ), and recovered (R k ) individuals, so that:

The lowercase indexing notation, k = {i, p}, is not to be confused with the uppercase indexing notation 417 P and I, which refer to the populations of individuals in the pre-symptomatic (P ) and symptomatic (I) 418 infectious compartments, respectively.

Let β P,k , β I,k , β A,k , and β H,k represent, respectively, the effective contact rates for pre-symptomatic, symp- to choose to travel to a patch with high daily COVID-19 mortality, and that those living in a patch with 428 higher daily COVID-19 mortality may opt to go to a patch with lower daily mortality. To account for these assumptions, the following inequalities must hold [47]:

Differentiating q kk + q kl = 1 (for k = l), with respect to D i and D p , reduces the inequalities in (3.1) to:

Based on the above assumptions and analyses, we choose the following functional forms for q kl (k, l = {i, p}), 432 which satisfy the two inequalities in (3.2): 

To keep track of the total number of residents in each of the two patches (countries) at time t, we let 440 M kl represent the total number of residents of patch k who are now in patch l (see Table 7 for the four 

is the total number of residents of India who can travel (to Pakistan) at any

is the total number of residents of India who can travel (but choose to remain in India), while I i + H i is the total number of symptomatic and hospitalized residents of India who 

Let λ kk be the rate at which residents of patch k acquire infection in patch k. Hence,

Similarly,

Furthermore, let λ kl represent the rate at which residents of patch k acquire infection while in a different patch l. Thus, Tables 1 and 2, respectively) .

(3.8)

For bookkeeping purposes, we let D k (t) represent the number of COVID-19 mortality in of residents of patch k at time t (with k = i, p). Hence, the rate of change of the population D k is given by: Populations of susceptible (S), exposed (E), pre-symptomatic (P ), asymptomatic (A), and recovered (R) compartments can travel between the two countries (hence, they are placed along the hypothetical "border" of the two nations). Symptomaticallyinfectious (I) and hospitalized individuals (H) cannot travel due to their ill health (hence, they remain in their respective patches).

The derivation of the (overall) reproduction number of the metapopulation model (denoted by R 0m ), using the next generation operator method, is presented in Appendix C. These analyses reveal that R 0m is expressed as the maximum of the constituent reproduction numbers for the residents of India (R 0ri ) and Pakistan (R 0rp ). That is, R 0m = max{R 0ri , R 0rp }. Table 4 , as well as the fixed parameters in Table 3 , are used). It should be 489 noted that we again assume S * N * ≈ 1 for each country, and estimate values for S * and N * using the population Table 4 , together with the fixed parameters tabulated in Table 3 ). (b) Lockdown measures are further strengthened to the extent that the baseline values of the community contact rate parameters (β P , β I , β A , β H ) are decreased by 20%. (c) Lockdown measures are mildly relaxed to the extent that the community contact rate parameters are increased by 20%. The reproduction number is computed using the next generation method in Appendix C. Values for fixed parameters for both countries are listed in Table 3 . Values for parameters estimated during lockdown (used as the baseline) are listed in Table 4a (for India) and 

The metapopulation model (3.8) will now be simulated to assess the impact of the back-and-forth mobility Table 8 . For simulation purposes, we consider the following four 529 mobility/dispersal scenarios (tabulated in Table 8 ), depending on the preference to stay in one patch, in 

Here, daily mortality in both countries declines steadily (Figures 7(b) and (d), brown curves).

On the other hand, if the lockdown measures are relaxed, corresponding to a 20% increase in the baseline 584 values of the community contact rate parameters, our results depicted in (Figures 7(a) and (c), gold curves) 585 show that while India will record a projected additional 46, 000 deaths (a 10% increase from the projected overall results obtained for these scenarios are qualitatively very similar to those obtained for the semi-626 symmetric 1 case (see Table 9 ). Specifically, we show that relaxing the control and mitigation measures, 

In the uni-directional 2 scenario, a brief outbreak occurs almost immediately after mobility begins at the 648 beginning of June 2021, and a small increase in daily mortality is observed later that same month with 649 continuous declines in cases and mortality thereafter (Figure 9 , purple curves).

The simulation results described in this section are tabulated in Tables 9 and 10 the likelihood for significant outbreaks (and higher mortality) in Pakistan. However, after an initial surge of infections due to increasing residence time residence of one patch spend in the other, the disease will eventually die out in both countries regardless of the mobility scenario (since the reproduction number, 665 R 0m , is less than unity in each mobility scenario). This result is consistent with the results reported in 666 Section 2.2 for the single-patch model (2.1).

667 Table 11 summarizes the projected time-to-elimination of the COVID-19 pandemic in India (Table 11a ) 668 and Pakistan (Table 11b) 2) and three effectiveness levels of lockdown measures (baseline, 20% reduction in β and 20% increase in 686 β), the reproduction number of the model is always less than unity (see Table 11 ). Furthermore, increased 687 accessibility and uptake of vaccines may shorten the time-to-elimination. Strengthening or relaxing of lockdown measures is measured in terms of a percent reduction or increase in the baseline values of the community contact rate parameters (β P , β I , β A and β H ) from their baseline values during lockdown (given in Table 4b ). Projected cumulative mortality are rounded to the nearest thousand. Table 8 . The reproduction number is computed using the next generation operator matrices given in Appendix C. Values for fixed parameters are listed in Table 3 . Values for parameters estimated during lockdown are listed in Table 4a (for India) and Table 8 .

This wave, which thoroughly overwhelmed India's healthcare system (in addition to making India the Proof. Consider the second equation of the basic model (2.1),Ė = λS − σ 1 E, which can be re-written in terms of the inequalityĖ ≥ −σ 1 E, with solution E(t) ≥ E(0) exp(−σ 1 t) ≥ 0. Similarly, the following inequalities can be established for the solutions of the next three state variables of the basic model:

It follows from the above inequalities that, P (t), I(t), A(t), and H(t) are non-negative for all time t > 0; thus it can be seen that the following inequality holds for S(t) as well:

Noting that the number of recovered individuals (R(t)) is increasing (sinceṘ depends on I, A, and H, all 764 of which are shown to be non-negative), then R is also non-negative for all t > 0). Thus, we have shown L = E + a 1 P + a 2 I + a 3 A + a 4 H, where a 1 = R 0 , a 2 = 1 φ+γ I +δ I φβ H γ H +δ H + β I , a 3 = β A γ A , and a 4 = β H γ H +δ H . The Lyapunov derivative is given by:L =Ė + R 0Ṗ + a 2İ + a 3Ȧ + a 4Ḣ , which can be simplified tȯ L = β P S N − R 0 σ 2 + a 2 rσ 2 + a 3 (1 − r)σ 2 P + β A S N − a 3 γ A A + β I S N − a 2 (φ + γ I + δ I ) + a 4 φ I in Ω whenever R 0 ≤ 1.

Here, too, the next-generation operator method will be used to compute the reproduction number for the 779 metapopulation model (3.8), denoted by R 0m . It can be seen that the associated matrix of new infection 780 terms of the metapopulation model (denoted by F ) is given by: .

Similarly, the associated matrix of linear transition terms in the infected compartments of the metapopulation model (denoted by V ) is given by: 

Publisher: World Health Organization

World Health Organization. Novel Coronavirus (2019-nCoV) Situation Report -1

World Health Organization

COVID-19: a comparison to 794 the 1918 influenza and how we can defeat it

Publisher: The Fellowship of Postgraduate Medicine preprint

World Health Organization. WHO Coronavirus (COVID-19) Dashboard, 2021. URL https: 799 //covid19.who.int

Worldometer. India population

First confirmed case of COVID-19 infection in India: A case report

CSSEGISandData/COVID-809 19

Library Catalog: www.who.int

World Health Organization. Coronavirus disease 2019 (COVID-19) Situation Report -46

Attributes and predictors of 821 long COVID

Long COVID and myalgic encephalomyelitis/chronic fa-825 tigue syndrome (ME/CFS)-A systemic review and comparison of clinical presentation and symp-826 tomatology

Number: 5 Publisher: Multidisciplinary Digital Publishing diseases/novel-coronavirus-2019/covid-19-vaccines

COVID-19) update: FDA authorizes additional monoclonal 876 antibody for treatment of COVID-19

Impact of non-pharmaceutical interventions (NPIs) to reduce 885 COVID19 mortality and healthcare demand

Micro-smart lockdowns imposed in different areas of Karachi's District Central

The News International

UP government imposes Sunday lockdown in entire state -The Economic Times

Maharashtra govt may impose 15-day lockdown starting Wednesday night

Lockdown news: Unlocking begins in some states, COVID curbs to remain in many others | India 899

Delhi extends lockdown by another week, restrictions to remain in force till

Unlocking: India states start reopening amid dip in COVID cases

ANI. Pak imposes 10-day countrywide lockdown from May 8 to control COVID. Business 909 Standard India

The Hindu. Pakistan to impose lockdown in Lahore, other cities as coronavirus cases witness 913 surge. The Hindu

Pakistan increases COVID restrictions amid a third wave

Umar Farooq. Pakistan deploys army in 16 cities to enforce COVID-19 precautions

922 [36] PTI. Several provinces in Pakistan reopen schools after drop in COVID-19 cases

Doctor in India: Emergency room is so crowded, 'it's nearly impossible to walk

The system has collapsed': India's descent into Covid hell

SARS-CoV-2 variant classifications and definitions

Substantial 938 undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

Publisher: American Association for 941 the Advancement of Science Section

Estimates of the severity of coronavirus 948 disease 2019: a model-based analysis. The Lancet Infectious Diseases

The India-Pakistan border is so closely guarded that it can be seen from space

Worldometer

Pakistan cuts last remaining transport link to India over Kash-1011 mir dispute. Reuters

Which country sends the most tourists to India?

About 0.9m tourists visited Pakistan in 10 months

Emirates extends flight ban to UAE from India, Pakistan, two other countries 1020 till

Estimated SARS-CoV-2 seroprevalence 1027 in the U.S. as of

Farzad Mostashari, Don Ol-1031 son

Estimation of excess deaths associated with the COVID-19 pandemic in the United 1033

Pakistan bans travellers from India for 2 weeks amid spread of new coronavirus vari-1037 ant

Travellers entering Pakistan as per guidelines by NCOC

Wing Yin Venus Lau, transmission of COVID-19 prior to symptom onset. eLife

Presymptomatic transmission of SARS-CoV-2 -Singapore

CDC Cruise Ship Response Team, California De-1059 partment of Public Health COVID-19 Team, Solano County COVID-19 Team

California Department of Public Health COVID-19 Team

Solano County COVID-19 Team

Public health responses to COVID-1072 19 outbreaks on cruise ships -Worldwide

Coronavirus 1077 pandemic (COVID-19) vaccinations

Appendix A Proof of Theorem 2.1

None.