Optimal strategic pandemic control: human mobility and travel restriction

1.   Introduction

The ongoing deterioration of 2019 novel coronavirus (COVID-19) pandemic highlights the challenges and trade-offs that many countries face in mitigating the epidemic. This paper studies an extended SIR model to find optimal mitigation policy when considering human mobility between the severely-affected and mildly-affected regions. Our work is motivated by one missing aspect in the related economic literature: the analysis of the economic impact of cross-region human mobility on pandemic mitigation. Human mobility across regions is an important factor behind the transmission of infectious diseases and has been modelled in some epidemiology literature (Sattenspiel and Dietz [1], Ferguson et al. [2], Hollingsworth et al. [3], Robertson [4], Lai et al. [5]). Cross-region human mobility, the movement of population from a severely-affected area to a mildly-affected one, increases the spread of disease during the outbreak of epidemics; for example, Sattenspiel and Dietz [1], Rvachev and Longini Jr [6], Wang and Zhao [7], Seno [8], Biswas and Mandal [9]. Consequently, policymakers have imposed travel restrictions as a response to pandemic threats (Bajardi et al. [10], Wang and Taylor [11], Charu et al. [12], Fang et al. [13], Kraemer et al. [14], Yabe et al. [15], Xiong et al. [16]). However, such restrictions on human mobility are not only controversial, but also cannot last long because of the negative economic impacts. Therefore, it is vital to quantify and understand the impact of human mobility, so far neglected in the economic literature. The result is critical for improving epidemic control and tailoring cooperative policy among regions.

In the COVID-19 crisis, due to the extremely fast trend after the initial phase of the epidemic (Remuzzi and Remuzzi [17], Yue et al. [18], Tang et al., [19]), many countries face much more serious constraints in the availability of resources and accessibility to medical equipment than during normal times (Duke et al. [20], Siow et al. [21]). This paper addressed the harsh reality faced by many policymakers—that they need to find the most effective way to jointly implement feasible mitigation policies with a limited budget of resources. This problem is significant not only because it is almost impossible to solely rely on one particular policy, such as imposing a complete lockdown to quarantine everyone, or complete testing to immediately find all infections, but also because the combination of different policies could be more efficient. Normally, two types of mitigation policies are jointly implemented in practice: non-pharmaceutical interventions (NPIs) like social distancing, quarantine, and lockdown; and pharmaceutical interventions (PHIs) such as testing and hospitalisation (Adhikari et al. [22], Charpentier et al. [23], Hou et al. [24], Jones and Adida [25], Tognotti [26], Nicolaides et al. [27], Piguillem and Shi, [28]). The NPIs are mainly implemented via isolation while the PHIs rely on the healthcare system. The COVID-19 pandemic has aroused economists' attentions on studying mitigation policies to reduce the economic costs as much as possible. For instance, Berger et al. [29] and Roques et al. [30] investigated the role of testing policy during the COVID-19 crisis. Atkeson [31] studied a lockdown planning problem to analyse the intensity and duration of the optimal lockdown policy. Bonaccorsi et al. [32], Hadjidemetriou et al. [33], Galeazzi et al. [34] discussed the economic and social consequences of the policy of travel restrictions in several countries. Jones et al. [35] studied the economic impact of the pandemic when the planner optimally imposed the policy of social distancing and encouraged working from home. However, most works mainly focused on the effectiveness of or optimal strategies regarding one particular type of mitigation policy, instead of the optimal combination under some constraints like limited resources. Few studies like Charpentier et al. [23], Piguillem and Shi [28], and Huang et al. [36] investigated the optimal combination of mitigating interventions. Thus, we contribute to literature by considering the cross-region human mobility between the severely- and mildly-impacted regions, which makes optimal mitigation policy more complicated.

We extend the SIR epidemiology model by considering several features: individuals' infection status can either be known or unknown, resources available for testing and the infrastructure for lockdown are limited, and there is human mobility across regions. Specifically, we followed a similar but simplified approach of Sattenspiel and Dietz [1] that fit cross-region human mobility into the SIR model. The human mobility is captured as a net population outflow from the severely-affected region to the mildly-affected one—a non-negligible factor in optimal policymaking. We consider how the policymaker of each region optimally allocates its limited resources into testing and travel restriction policies to contain the pandemic under the objective of minimising economic cost. Furthermore, we provide a comprehensive numerical analysis on the optimal resource allocation policy and mitigation outcome. Our analysis is based on parametrising the model using estimates from the COVID-19 pandemic. The results of benchmark scenario and sensitivity analysis show that the efficient use of resources depends on many factors such as pandemic progression, duration of mitigation policies, and cross-region human mobility.

Our study makes some potential contributions to the economic literature on mitigating pandemics like the most recent COVID-19. First, we extend the literature that mainly focuses on particular types of epidemic intervention policies, by formulating a pandemic control problem with limited resources to be allocated to different mitigation strategies. Second, although cross-region human mobility has been investigated in epidemiology literature, it lacks a counterpart in the economic context. This paper attempted to bridge human mobility and optimal mitigation through an economic objective. Lastly, our study is effective in analysing the economic impact of human mobility and ongoing travel restrictions during the COVID-19 pandemic.

This work is first linked to the ongoing discussion on the optimal policy response to the COVID-19 pandemic, such as Shen et al. [37], Tang et al. [38], Tang et al. [39], Barro et al. [40], Dewatripont et al. [41], Eichenbaum et al. [42], Hall et al. [43], and Baldwin and di Mauro [44]. Our paper is especially closely related to Berger et al. [29], Atkeson [31], and Jones et al. [35] that aim to find the optimal way of implementing mitigation policies to control infectious diseases. In particular, we share a common scope with Charpentier et al. [23] and Piguillem and Shi [28], who investigated the optimal balance between testing and quarantine policy interventions. Charpentier et al. [23] considered a modified SIR model and used optimization algorithms to compute the 'best' level of lockdown and testing. They showed that, when massive amounts of resources are introduced to detect infected individuals, the pressure on social distancing can be eased. Piguillem and Shi [28] attempted to discover the optimal COVID-19 quarantine and testing policies from the perspective of a planner who tries to stop the diffusion of the disease. However, the distinction lies in that we consider the constraint of limited resources for implementing a mitigation approach, while Charpentier et al. [23] included the limited intensive care unit (ICU) capacity and Piguillem and Shi [28] optimized the strategy to maximise social welfare.

Second, we also join the economic literature that extend epidemiology models like the SIR model to analyse pandemic transmission with specific economic objectives. A detailed introduction to compartment models can be found in Hethcote et al. [45], Anderson and May [46], Brauer et al. [47], Lefevre [48], Tang et al. [49], Xue et al. [50]. Most studies supplemented economic considerations into their epidemiological models to study mitigation policies and economic impacts; for example, Tognotti [26], Atkeson [31], Feng and Garrido [51], Chen et al. [52], Ferguson et al. [53].

The remainder of this paper is organised as follows. Section 2 describes our model setup based on the SIR model and Section 3 presents the pandemic mitigation setting. Section 4 studies the optimal mitigation policy in the benchmark scenario and sensitivity analysis. Section 5 further discusses the economic impact of travel restrictions. Section 6 contains the main conclusions of the manuscript.

2.   A two-region SIR model with cross-region mobility

2.1.   Model setup

We build our model based on the SIR model that characterises the interaction among three compartments. Consider the number of population inhabiting a region at the time $t$, $N(t)$. The SIR model separates the whole population into three groups: susceptible people $S(t)$, infected people $I(t)$, and recovered people $R(t)$. $S(t)$ denotes the group of susceptible individuals who have no immunity to the studied disease; $I(t)$ represents infected individuals who are able to transmit the disease when susceptible individuals come into contact with them; and individuals removed from the group $I(t)$ owing to recovery are in class $R(t)$. The recovered group $R(t)$ is immune after recovery and does not include those who died from the pandemic. Here, we follow the common practice of literature to assume that the recovered group $R(t)$ is immune after recovery, see for example, Sun et al. [54], Gatto et al. [55].

We extend the baseline SIR model by considering some major extensions as follows:

Infection status: unknown infected and known infected. The fast transmission of contagious disease makes it impossible to distinguish all infected individuals, especially the asymptomatic ones, from the susceptible ones. This is mainly because some infected are asymptomatic or weakly symptomatic in the early stages of the disease and some others may exhibit no or indistinguishable symptoms throughout their whole infection status. In our model, we allow the infection status of the infected to be unknown at first and further assume that the authentic infection status is then revealed after the individual gets tested or becomes symptomatic. Thus, the infected group $I(t)$ is then further divided into two subgroups:

● Unknown infected (UI): Infected individuals that are asymptomatic and not yet tested. The UI group is indistinguishable from the susceptible group, but is able to transmit the disease.

● Known infected (KI): Infected individuals that are tested or deteriorate to show symptoms. We assume that the KI group is well quarantined and under treatment in hospital. They are not infectious as they cannot physically come in contact with the susceptible people.

There is a special case for KI group, which is the situation of home quarantine. In many countries, people who are diagnosed as asymptomatic or infected with very mild symptoms are usually required to be isolated at home. Compared with hospital quarantine which is the main way to well isolate susceptible and infected people, the situation of home quarantine will cause a high exposure risk due to cohabitation and activities to obtain basic daily necessities. However, due to the implementation of epidemic controlling measures, the home quarantined infected cannot travel cross regions and contact susceptible people as frequently as UI group. Therefore, for simplicity, once an infected is transferred from UI to QI, whether due to the onset of symptoms or the positive testing result, we assume that the infected person can be completely isolated in the hospital. This assumption is reasonable. On the one hand, the symptomatic infections can receive the timely appreciate treatment in normal hospitals. On the other hand, people can take measures such as building mobile cabin hospitals or temporary field hospitals to provide safe isolation places and necessary medical interventions for asymptomatic and very mildly symptomatic infections.

Region difference: severely-affected and mildly-affected regions. We consider there are two regions: a severely-affected region (region A) and a mildly-affected region (region B). The concept of region in our model may correspond to cities, states, or even countries in the real world, whose geographic segmentation only allows the residents to make physical contact with those in the same region. At time $0$, the observed epidemic situation in region A is worse than that in region B. Namely, the KI cases in region A, $I_{A}^{K}(t)$, is more than that in region B, $I_{B}^{K}(t)$. This assumption is realistic because the pandemic situation in the region that first experiences the outbreak remains worse than in other regions that experience the outbreak for a short time horizon like a few months. Therefore, at any time $t$, the total population in a region is

Cross-region human mobility. Human mobility across different regions, the fact that people travel from one region to another, contributes to the transmission of infectious diseases. The outbreak of epidemic in one region is often accompanied by panicked residents fleeing from the more severely-impacted region (Ferguson et al. [53], Shi [56]). Indeed, many policymakers have imposed restrictions on human mobility flows as one of their mitigation response (Bajardi et al. [10], Wang and Taylor [11], Charu et al. [12], Kraemer et al. [14], Yabe et al. [15], Xiong et al. [16]).

Following an approach similar to that of Sattenspiel and Dietz [1], which fits mobility process among discrete geographic regions into SIR models, we include a process to describe the net human mobility flow of S and UI groups from region A to region B. This is because it is reasonable to assume that during a pandemic outbreak, more people temporarily move out from the more severely-impacted region than the other way around.

We assume that the net mobility flow increases with the contrast of QI cases between regions A and B. Moreover, we assumed that the initial population of regions A and B are equivalent, and thus the comparison of QI cases between regions A and B is consistent with the contrast of QI population density. Moreover, in order to focus on the economic impact of human mobility and limited resources, when performing simulation analysis, we set several parameters such as the initial population, infection fatality rate, and natural newborn population in the two regions to the same value (as shown in the Table 2). The similar settings can also be found in Hethcote [57], Arino and Van den Driessche [58], Sanders et al. [59]. This assumption is reasonable because it is echoed by the realistic observation (Kraemer et al. [14], Read [60], Tian et al. [61], Xu et al. [62]). The similar model setting can be found in Sattenspiel and Dietz [1], Rvachev and Longini Jr [6], Seno [8], Lee et al. [63], Wang and Zhao [7,64], Nakata and Rost [65]. In this article, we measure the intensity of mobility rate of S and UI groups by $\lambda_{t}$. The amount of net cross-region flows of S and UI population per unit of time at time $t$ are $\lambda_{t}S(t)$ and $\lambda_{t}I_{A}^{U}(t)$, respectively. Specifically, $\lambda_{t}$ is driven by the gap of infectious cases between A and region B. We assumed the human mobility rates for S and UI groups are

$\overline{\lambda^{U}}$ and $\overline{\lambda^{S}}$ are the upper limits of UI and S groups' mobility rate between two regions. $f(x)$ is a monotone increasing function with $f(x)\in\lbrack-1, 1]$ and $f(0) = 0$. $c$ is a sufficiently small positive constant. $c$ is set to ensure that we can use the next generation matrix technique (Van Den Driessche and Watmough [66]) to calculate the basic reproduction number $R_0$ in Section 3.3. The calculation of $R_0$ is based on the disease free equilibrium (DFE): $(I_{A}^{U}, I_{A}^{K}, I_{B}^{U}, I_{B}^{K}, S_{A}, R_{A}, S_{B}, R_{B}) = (0, 0, 0, 0, S_{A}^{0}, 0, S_{B}^{0}, 0)$ where $S_{A}^{0} = \frac{A_{S}}{d+\underline{\lambda}}$, $S_{B}^{0} = \frac{B_{S}+\frac{\underline{\lambda}A_{S}}{d+\underline{\lambda}}}{(d+\underline{\lambda})d}$. Therefore, to ensure $\lambda_{t}$ exists when $I_{A}^{K} = I_{B}^{K} = 0$, $c$ is set to be a sufficiently small positive constant. More details can be found in Appendix A. For simplicity, our baseline model first considered that the mobility rates of UI and S groups are equal, $\overline{\lambda^{U}} = \overline{\lambda^{S}} = \overline{\lambda}$, and thus

Equation (2.1) shows $\lambda_{t}$ increases with the gap of KI cases between regions A and B, and a negative $\lambda_{t}$ represents the case when there is a net flow of population from region B to A. The cross-region human mobility process has the following features: 1) the cross-region mobility happens between different regions; 2) the net flow is from the more severely-impacted region to the less impacted one; 3) the intensity of mobility flow is dynamic, and it varies with the pandemic processes in both regions; and 4) the maximal mobility rate is restricted by the transportation capability. Now, we show that these features of our human mobility process are rooted in reality and have foundations in literature.

First, many models have attempted to incorporate the contact and human mobility of populations among discrete geographic regions. There are mainly two types of model structures found in the literature: 1) the structures with cross-regional contact rate (Gatto et al. [55], Takeuchi et al. [67], Keeling et al. [68], Cui et al. [69]); 2) the structures with cross-regional mobility rate rather than contact rate (Rvachev and Longini Jr [6], Sattenspiel and Dietz [1], Seno [8], Wang and Zhao [7,64], Nakata and Rost [65]). In addition to these, some works adopt a form of joint structure that combines the two methods, see for example, Sattenspiel and Herring [70] and Sattenspiel and Herring [71]. Generally, the first method does not consider changes in the regional population size, while the second structure, which our paper adopts, allows such changes.

Second, one major commonality of the epidemiological literature with cross-region mobility is that the spread of pandemic increases with the deterioration of situations in some severely-impacted regions, given the constant (or dynamic) contact rates and mobility rates among different regions they adopt. This is echoed by the realistic observation that the outbreak of pandemic always accompanies the accelerated net mobility outflow from the severely-impacted region to mildly-impacted region (Kraemer et al. [14], Read [60], Tian et al. [61], Xu et al. [62]). For example, Tian et al. [61] noted that, before the lockdown policy of Wuhan, China, the number of recorded movements from Wuhan to other regions of China increased with the number of COVID-19 cases. Moreover, some literature describes the form of mobility outflow in their mathematical epidemic models (see for example, Sun et al. [54] and Hethcote [57]). In particular, in the Hethcote [57] model, the mobility rate is linearly related to the relative infection size between two regions.

Third, the cross-region mobility process has also long been captured as a dynamic process in the literature (Sattenspiel and Dietz [1], Rvachev and Longini Jr [6], Seno [8], Sun et al. [54], Hethcote [57], Lee et al. [63], Wang and Zhao [7,64], Nakata and Rost [65]). These studies model the mobility rates to be related to time, infections size, transportation restrictions, and so on (see for example, Lee et al. [63]). More complicated models even combine the ordinary differential equations or numerical network dynamics models to capture the mobility or migration process (O'Sullivan et al. [72]).

Last, the cross-region human mobility flows between regions are largely set by the intensity of transportation links between regions (see e.g., O'Sullivan et al. [72]). The maximal mobility rate is measured by the ratio of population in a region that can be transferred. During the COVID-19 crisis, the intensity of cross-region transportation like flights and trains are largely reduced. To study the effect of mobility rate, our model considered a range of possible rates instead of a particular rate.

2.2.   Transmission dynamic

To begin with, we introduce the fatality and recovery pathways of our model in Figure 1. In our model, the recovery pathway of the three groups of population introduced above is

Any infection is first asymptomatic and unknown, and then possibly turns into the KI state. However, both UI and KI groups may recover from infection. Moreover, regarding the fatality path, we are particularly interested in death caused by infection. While all groups face a common natural mortality rate $d$, the UI and KI groups face an additional fatality rate due to infection, $d^{U}$ and $d^{K}$. Thus, the fatality pathway due to the infection is

Next, we introduce the transmission dynamic of different groups considering human mobility between two regions. For simplicity of notation, we leave out the notations of time $t$. However, all the group population parameters are dynamically varying with time. First, we introduced the dynamic motion of susceptible population within region A. First, the new infections are caused by the effective contact between the susceptible people and infected people. Consistent with the SIR model, we measure the number of susceptible that become infected per unit of time as $\beta_{A}S_{A}I_{A}^{U}$, because the KI group are not infectious. Here, $\beta_{A}$ is the infection rate when there is no lockdown and people make contacts at the normal rate. The motion of the population of susceptible agents per unit of time is

where $A_{S}$ is the increment of susceptible individuals owing to natural birth in region A, and $d$ is the natural death rate.

Second, we consider the dynamic motion of UI and KI populations. The susceptible people that get infected first acquire UI status, who later change to KI status with an intensity of $\varepsilon$ fraction per unit of time owing to testing or deterioration of symptoms. The $\varepsilon$, called the testing rate hereafter, could be reinforced with policymakers' increased testing. Moreover, both the UI and KI groups face an additional fatality rate, $d^{U}$ and $d^{K}$ on top of the natural one, and their recovery rates are $v^{U}$ and $v^{K}$, respectively. After considering the cross-region mobility, that is, outflows for regions A and B, we utilised the following:

where $\lambda_{t}$ is the mobility intensity rate of S and UI groups, respectively. Here, we assume that the fatality rate of QI infected is lower than that of UI infected. This assumption is reasonable because, as we mentioned above, both asymptomatic infections and symptomatic infections can be treated promptly in the hospital once they are diagnosed.

Lastly, the population of the recovery group increases with the improvement of the infected and decreases owing to a natural death,

Then, the transmission dynamic of region B can be easily implied by that of region A. The global pandemic transmission process can be captured by the following nonlinear dynamic system:

We summarise the model transmission dynamic in Figure 2 and notation of all parameters in Table 1. Note that Figure 2 illustrates the case when region A is more severely-impacted and omits the process of natural newborn and death.

3.   Pandemic mitigation strategy with limited resources

We consider that both regions are planning to mitigate the pandemic during $[0, T]$ with the endowed resources $\Omega_{i}$. $T$ is the duration of the mitigation period, measured in days. The amount of resources $\Omega_{i}$ is divided evenly for each day and then the daily resources of region $i$ is $E_{i}$ with $\Omega_{i} = E_{i}T$. The planner's resource allocation is to implement an efficient use of resources by combining the available mitigation approaches during $[0, T]$.

As seen during past pandemic outbreaks, both PHIs and NPIs are effective mitigation methods (Adhikari et al. [22], Hou et al. [24]). PHIs include the testing, medical treatment, and quarantined hospitalisation to detect and cure the infections, as well as vaccination. NPIs are mainly via reduction of public exposure to possible infection chances, such as wearing masks, social distancing, and government forced self-quarantine or lockdown. Though the PHIs such vaccines and medicines are more effective, their availability is unlikely to be adequate during the early stages of a pandemic (Martinez and Das [73]). By contrast, NPIs are among the best ways of controlling pandemic outbreak when PHIs are not yet massively available.

The implementation of both PHIs and NPIs approaches are costly and require the input of resources. For instance, according to the experience of many countries during the COVID-19 crisis, the cost of PHIs such as purchasing medical equipment, building hospitals, and hiring personnel, is heavily covered by the government. NPIs like social distancing, quarantine, and school and workplace closure also incurs the costs of distributing necessary food, preventive materials, and making financial compensations to ordinary households and the unemployed.

Therefore, with a limited budget of resources, each region considers implementing a policy combining both PHIs and NPIs, to optimally mitigate the pandemic. Mathematically put, we assume $\alpha_{A}$ ($\alpha_{B}$) ratio of $E_{A}$ ($E_{B}$) is allocated into implementing the PHIs in region A (B), and the remaining resources is then used for implementing NPIs (see Figure 3).

The resources $\Omega_{i}$ is restricted and thus the maximal duration $T$ of implementing PHIs and NPIs is also limited. The assumption of $T$ is realistic because most countries cannot afford an inordinately long period, such as six months, during which they implement some intensive mitigation policies like a lockdown. Specifically, we considered the value of $T$ to be around $150$ days in the later benchmark scenario, based on the experience of the policy horizon in many countries during the COVID-19 crisis.

3.1.   Mitigation approaches

3.1.1.   PHIs: testing & hospitalisation

In this study, PHIs refer to the testing & hospitalisation (hereafter abbreviated as "testing" for brevity) policy only. We do not include vaccination because of the inability of the current technology to develop and deploy a vaccine in a short amount of time (Martinez and Das [73]). Here, we specially focus on the policymaking in the relative early stage of the epidemic where the vaccination is unavailable. This is because the epidemic control measures in the initial phase of the infection outbreak are critical and would significantly influence the subsequent evolution of the pandemic. According to the experience of many countries, the initial phase of the infection outbreak was followed by an extremely fast trend (Remuzzi and Remuzzi [17], Yue et al. [18], Tang et al. [19]). The rapid deterioration of the epidemic is considered as an unprecedented and tough challenge to the health system's capacity of most countries and will push the national health system to full capacity. Facing the fast developing epidemic, the severe shortage of medical resources is arising in more and more countries, and it has seriously hindered the well quarantine and timely hospitalization of the explosively increasing infected. Testing policy can distinguish the UI group from the S population. If tested positive, the UI individual turns to KI state and is also quarantined in the hospital to receive medical treatment. We assume that the fatality rate of the KI group is lower than that of the UI group. Thus, testing policy helps identifying the infected individuals and eliminating their transmission capacity.

During the mitigation period $[0, T]$, the testing rate $\varepsilon_{i}$ in region $i$ increases with the amount of resources allocated to it, $\alpha_{i}E_{i}$. Specifically, we further assumed the effect of resources on the rate $\varepsilon_{i}$ is linear, and set the coefficient to be

where $\varepsilon_{0}$ is the bottom rate at which infected individuals naturally transit from the UI to the KI group, due to an incubation period for the infected individuals to transit from asymptomatic to symptomatic cases (Wu et al. [74]). We do not consider the congestion in the healthcare system because we want to illustrate our main points without complicating other settings. The result of our benchmark scenario in Section 3 shows that the peak of infections, around 1 case per thousand, is within the capacity in many countries during the COVID-19 crisis. Here, we need to distinguish the difference between the parameter $\varepsilon_i$ in our model and the transmission rate from Exposed (E) to Infected (I) in SEIR model. As we mentioned above, by implementing policies of increasing testing, $\varepsilon_i$ can be dramatically increased. Therefore, $\varepsilon_i$ is not determined by the nature of the infected individuals solely. Through continuous, large-scale testing, even those infected who would be asymptomatic in the whole infectious period can be found. Thus, this parameter $\varepsilon_i$ is essentially different from the transmission rate in SEIR model. $k_{\varepsilon}$ is the sensitivity coefficient describing the marginal effect of daily resource input on increasing the testing rate. We consider the maximal $\varepsilon$ is always lower than $1$ even though all resources are allocated into it.

3.1.2.   NPIs: lockdown & quarantine

NPIs, also known as lockdown & quarantine (hereafter abbreviated as "lockdown" for brevity) policy in our study, restrict people's movement and prevent their potential exposures. Consistent with Alvarez et al. [75], we assumed that the lockdown policy in region $i$ puts $l_{i}\in\lbrack0, L]$ fraction of the S and UI population, $(S_{i}+I_{i}^{U})$, quarantined at home during the mitigation period $[0, T]$. Note that $L\leq1$ because the government cannot lock down the whole population as the economic activities, such as energy and basic food supply must continue. We follow the parameter that $L = 0.7$ in Alvarez et al. [75]. We assume that those recovered can be identified and are not under the lockdown policy. As the $l_{i}(S_{i}+I_{i}^{U})$ population is quarantined, the number of new infections per unit of time becomes

Thus, as those in lockdown cannot get infected nor infect others, the effective infection rate $\beta_{i}$ in this region decreases with the lockdown ratio $l_{i}$ in a quadratic function,

where $\beta_{0}$ is the infection rate without implementing the lockdown policy, namely $l = 0$.

Furthermore, we consider the effect that resources put on the lockdown policy during $[0, T]$. For ease of illustration, we assume that the lockdown ratio $l_{i}$ linearly increases with its allocated resources $(1-\alpha_{i})E_{i}$,

where $k_{l}$ is the sensitivity coefficient of lockdown ratio $l_{i}$.

3.2.   Basic reproduction number

The concept of the basic reproduction number in epidemiology, denoted by $R_{0}$, measures the contagion rate of an infectious disease. $R_{0}$ indicates the expected number of infected cases directly generated by one existing infected case (Diekmann et al. [76]). We emphasise that $R_{0}$ is important for policy intervention. It is expected that the virus would gradually disappear when $R_{0} < 1$; otherwise, it never dies.

For preparation, we denote the net human mobility when $I_{A}^{K} = I_{B}^{K}$ by $\underline{\lambda}$, the net outflow rate when the KI cases are equal in both regions. According to Eq (2.1), we have

which indicates people's "confidence" in region A. For instance, $\underline{\lambda} > 0$ is the case where people have less faith in the region A and still flee to region B though they observe the same KI cases in both regions. Now, we calculate the basic reproduction number in both regions based on our disease transmission dynamic.

Theorem 1. The basic reproduction numbers of each region, $R_{0}^{A}$ and $R_{0}^{B}$, are

Some parameters of $R_{0}^{A}$ and $R_{0}^{B}$ are determined by the characteristics of the disease or by nature, and others can be changed by disease mitigation policies. The former is difficult to change in a short period of time by people's efforts like the natural and infection fatality rate; while the latter, such as the infection and testing rates can be significantly reduced via mitigation measures. As expected, both $R_{0}^{A}$ and $R_{0}^{B}$ are positively proportional to $\beta$, while they are negatively proportional to $\varepsilon$. This suggests that both the PHIs and NPIs are effective in slowing down the pandemic transmission.

Corollary 1. Consider the effect of human mobility rate $\underline{\lambda}$, we have

Regarding the impact of the human mobility rate, Corollary 1 reveals an interesting result that in general $\underline{\lambda}$ relieves the pressure of region A, and causes challenge for region B. In the particular case when $\underline{\lambda} = 0$, $R_{0}^{A}$ and $R_{0}^{B}$ can be expressed in a unified form. In this case, $R_{0}^{A}$ and $R_{0}^{B}$ are equivalent to the cases without human mobility between two regions. Thus, from the perspective of $R_{0}$, it is the minimal human mobility rate $\underline{\lambda} = 0$ instead of the maximal rate $\overline{\lambda} = 0$, that plays an effective role in mitigating the pandemic transmission.

Proposition 1. The basic reproduction number $R_{0}$ of the whole two-region system is

The form of $R_{0}$ is consistent with epidemiological literature with similar discrete population segmentation, see for example, Van den Driessche and Watmough [66]. This discrete feature of $R_{0}$ may be hard to understand at first glance, however, it has a reasonable explanation. The concept of $R_{0}$ indicates the infection ability of disease in the whole two-region population, and thus is more determined by its infectiousness in the worse region. For instance, in the extreme case when region B fully controls the spread of the pandemic with $R_{0}^{B} = 0$, the infection ability of the disease is then $R_{0} = R_{0}^{A}$ because it only exists in region A.

3.3.   Optimization under limited resources

To begin with, we formulate policymakers' optimization problem from an economic perspective. At the current time $0$, the planner of region $A$ or $B$ needs to determine the optimal allocation of resources for implementing the PHIs and NPIs. They solve for the optimal mitigation strategy by minimising an objective of economic cost,

where $\alpha_{i}\in[0, 1]$ is the resources allocated into PHIs in region $i\text{} = \text{ }A, B$. For the planner of region $i$, the choice $\alpha_{i}$ is the control variable. Regarding the objective function, the economic cost of pandemic mitigation policy includes the lockdown cost and death cost :

Lockdown cost refers to the economic cost of production suspension owing to lockdown; and the death cost is the value of statistical life (v.s.l.) caused by fatalities. The planner faces trade-offs when trying to minimise the overall economic cost. On the one hand, the planner needs to optimally balance the resource put into the lockdown policy and testing policy. On the other hand, there is also a trade-off in the time horizon. For instance, the lockdown policy increases economic cost today as quarantined people cannot produce. However, it lowers the death cost in the future as it slows the transmission of the pandemic.

In the optimization problem, the planners face the constraints of budget and maximal duration of mitigation policy. The goal of optimization problem is to maximise the effect of limited resources under the planner's objective function. In the unrealistic case of unlimited available resources to be allocated to mitigate the pandemic, the optimal solution would become a simple but extreme strategy; that is, complete testing or full quarantine. However, we observe that many countries are facing budget constraints and costs to take measures for slowing the spread of the pandemic. For instance, many countries lacked medical equipment and fiscal funding for mitigating COVID-19. It is also a common assumption in the literature on pandemic interventions that the social planer's available resources are limited. Moreover, our assumption of a limited time horizon comes from the observation that it is impossible for planners to implement mitigation policies like lockdown for a long period, such as more than 12 months.

3.4.   Economic cost

We adopt the economic setting of Alvarez et al. [75], in which each individual produces $w$ units of output per day when there is no lockdown. First, the economic cost of implementing a lockdown policy in region $i$ that puts $l_{i}$ fraction of S and UI groups in quarantine equals $wl_{i}(S_{i}+I_{i}^{U})$. This is because those in quarantine, that is, $l_{i}(S_{i}+I_{i}^{U})$, are not able to produce. Second, the economic cost of fatalities caused by infection is measured by the v.s.l., denoted by $\eta$. It is the discounted value of fatalities' future economic output and often measured in units of annual GDP per capita. For instance, Alvarez et al. [75] adopt $20$ times annual GDP per capita for the fatalities. Thus, the planner of $i$ region faces the following economic cost at time $t$

where $\dot{D}_{i}$ is the number of deaths caused by infection per unit of time in region $i$.

The planner aims to minimise the economic cost within a time horizon $[0, T]$, where the future costs are discounted at the rate $r > 0$. For region $i$, the economic cost that the planner aims to minimise is

$\hat{T} = \inf \{ t | I_{A}(t)+I_{B}(t) = I_{min} \}$, where $I_{min}$ is a constant between 0 and 1. $\hat{T}$ can be regarded as the end of the pandemic period ^*. Intuitively, the time span $\widetilde{T}$ represents the periods taken into consideration by the planner. Here, the reason for not considering an infinite time horizon is because we assume that the pandemic will eventually disappear after the deployment of a vaccine at some point. In other words, $\widetilde{T}$ is the time point when a vaccine is produced. Our time horizon setting is consistent with Djidjou-Demasse et al. [78] who study the optimal COVID-19 pandemic control until vaccine deployment. We optimistically considered $\widetilde{T} = 365$ days in our simulation in Section 3.3. This estimation is based on the fact that it took nine months to develop and produce a vaccine during the 2009 H1N1 outbreak (Halder et al. [79], Aunins et al. [80]). However, instead of implementing mitigation policy till $\widetilde{T}$, we considered the mitigation period for policymakers to be within $[0, T]$, where $T \in(0, \widetilde{T}]$. This is mainly because the stock of available resources are always limited and cannot last till the development of a vaccine. It is not realistic that the mitigation policy like lockdown could last a whole year. Obviously, the choice of $T$ affects the economic cost to a large extent. In Section 4, we conduct some sensitivity analysis to investigate the effects of different $T$.

^*Following an approach similar to that of Hansen and Day [77].

We explain the motive of the consideration that the social planners adopt time-invariant optimal mitigation policies during the time horizon of mitigation. The major reason for this setting lies in the fact that it is almost impossible for the planner to adopt a highly frequent dynamic mitigation policy in reality. For instance, though the planning under optimal control technique can determine time-varying mitigation strategy, the government can hardly accurately vary the daily number of workers in lockdown. That is, though technically we can have a time-varying allocation of resources, the constant mitigation strategy is more practical, especially for a short mitigation horizon.

4.   Strategic optimal pandemic mitigation

4.1.   Strategic optimization settings

In this study, we consider the simultaneous mitigation policy-making setting that the planners of two regions make their decisions concurrently without knowing in advance the actions chosen by the other. Mathematically put, we aim to find the optimal solution pair $(\alpha_{A}^{\ast}, \alpha_{B}^{\ast})$ in $[0, 1]\times[0, 1]$ satisfying

We explain the reason for our optimization setting. The decision-making setting that one social planner per region is responsible for making mitigation policy in the region can be supported by some observations of different geographical regions (e.g. cities, states, or countries) during the COVID-19 pandemic. Very often, in reality, these geographically-separated regions have different fiscal budgets and are responsible for making their own pandemic contagion policies. We observed that many states, provinces, and countries responded differently to the ongoing COVID-19 crisis. Thus, our model considered the optimization setting according to the fact that the planners are responsible for mitigating policies in their own regions when facing the same outbreak. It is noteworthy that some other model settings, like the existence of a central planner for coordinating the policies of two regions, may also correspond to some realistic observations. However, owing to the limited scope and length of our paper, we are not able to investigate all patterns and thus only focus on the most common institutional setting that each region makes their own decisions. In our concluding remarks, we further note how future studies can extend our setting by considering other optimization settings.

The problem (4.1) is a canonical simultaneous move game. In the search for the optimal ratio of resources to be allocated into PHIs, the planners of region $A$ and $B$ face constraints from the two-region model and limited resources budget. Specifically, given the settings, region $i$'s objective function in model (3.4) becomes

The solution is the Nash equilibrium in which no planner has a profitable deviation given the actions of the other planner. Thus, we use the classical best response function technique to find the equilibrium in which $\alpha_{A}^{\ast}$ is a best response to $\alpha_{B}^{\ast}$ and $\alpha_{B}^{\ast}$ is a best response to $\alpha_{A}^{\ast}$.

The process for solving the pair $(\alpha_{A}^{\ast}, \alpha_{B}^{\ast})$ taken within the set $[0, 1]\times[0, 1]$ is given in the Appendix (see e.g., Cave [81], Osborne et al. [82], Morris [83]). For our model setting, the numerical result shows that the equilibrium exists and is unique.

4.2.   Parameter calibration of the benchmark scenario

We parameterise the benchmark scenario using a range of estimates on the COVID-19 pandemic from some existing related literature. Except the initial pandemic progression conditions, we considered the following benchmark parameters to be the same in region A and region B.

Initial conditions. Suppose the initial populations in both regions are equivalent, $N_{A}(0) = N_{B}(0) = 8\times10^{6}$. The 8 million population size is roughly that of Wuhan in February 2020. Moreover, we considered that the KI cases at time $t = 0$ in regions A and B are $7500$ and $2500$, respectively. The UI groups of two regions are assumed to be $4$ times of KI groups based on the practical experience of pandemic outbreak. Thus, the region A is relatively more severely impacted than region B at time $t = 0$.

Disease parameters. To calibrate the infection rate $\beta_{0}$ without lockdown, we use the value of $6.25\times10^{-8}$ based on the estimation of Tang et al. [19] using data from Wuhan, as our population size $N_{i}(0)$ is also based on that of Wuhan.^† Regarding $v^{K},$ $v^{U},$ $d^{K}$ and $d^{U}$, we refer to Zhou et al. [84] and Gatto et al. [55] for COVID-19 pandemic parameters. Zhou et al. [84] finds the median duration of recovery in hospitalisation is 8 days, and Gatto et al. [55] reports that without hospitalisation is 14 days. Thus, we adopt the estimations that $v^{U} = 1/14$ and $v^{K} = 1/8$. In addition, according to the COVID-19 data from Wuhan, the average fatality rate of those with and without hospitalisation are 2% and 20%, respectively. As the median duration from infection to death is 11 days (Zhou et al. [84]), we set $d^{K} = 0.02/11$ and $d^{U} = 0.2/11$. Last, we adopt the annual natural death rate (0.7%) and birth rate (1.2%) in China to calibrate the daily death rate $d$, newborn $A_{S}$ and $B_{S}$ in both regions.

^†Note that $\beta_{0}$ in our model is equivalent to the $\beta c_{0}$ in the model of Tang et al. [19].

Human mobility. For simplicity of illustration, we consider a simple linear form of $f(x) = x$ in Eq (2.1), thus

The adoption of the linear form results from both the literature foundations and its simplicity. As introduced earlier, the infection risk in a mildly-impacted region is increasing with the number of infection cases in other regions (see e.g., Hethcote [57], Rvachev and Longini Jr [6], Sattenspiel and Dietz [1] Wang and Zhao [7,64], Nakata and Rost [65], Seno [8], Sun et al. [54]). Specifically, Hethcote [57] adopted a linear form to describe the human mobility process, which is determined by the relative scale of infections between the mildly-impacted regions and severely-impacted regions. In Hethcote [57], the human outflow from the mildly-impacted region $i$ to the severely-impacted region $j$ was in the following form:

where $\theta$ is a positive coefficient, $I_j(t)$ and $I_i(t)$ are the number of infections in each region, and $N(t)$ is the region $i$'s population size. Our linear function is inspired by the literature though we model the net human mobility flow. As robustness tests, we also examine some simple non-linear functions function of $f(x)$, including the convex function $f(x) = sign(x)\cdot x^2$ and concave function $f(x) = sign(x)\cdot |x|^{1/2}$, where $sign(x)$ indicates the positiveness or negativeness of $x$. The results show that our main conclusions remain robust which are not reported owing to the length of paper.

Moreover, we estimate the maximal human mobility rate $\overline{\lambda}$ based on the net population outflow of Wuhan within January 2020. As of January 23, 2020, when Wuhan announced its lockdown, an estimated 5 million of the total 14 million residents had already left (Ai et al. [85], Jia et al. [86]). This surprisingly high outflow rate, 1/3 of the total population in one to two months, might be exaggerated by the effect of Chinese Spring Festival and being the first epicentre of the COVID-19 pandemic. Thus, we set a much milder rate of $\overline{\lambda} = 40\%/365$. In this case, the maximal net outflow is at a rate of annually $40\%$ population of region A. In the sensitivity analysis in Section 4, we further study the outcomes under some different $\overline{\lambda}$.

Resource endowment. Regarding the available resources, we considered the endowments for regions A and B were the same in our simulation, $\Omega_{A} = \Omega_{B} = \Omega$. In addition, the daily resources in both regions is then equal $E_{A} = E_{B} = E$. We further normalise $E = 1$ unit and set the duration of mitigation period as $T = 150$ days. Thus, we have $\Omega = 150$ units.

Mitigation policy parameters. Without directly available data on the cost and effect of testing and lockdown policy, we made rough estimations on the mitigation policy parameters. On the premise that the daily resources have been normalised as $E = 1$ unit, we need to properly set the ratio of $k_{\varepsilon}$ and $k_{l}$. This is because the absolute values of $k_{\varepsilon}$ and $k_{l}$ depend on the unit of $E$. For instance, enlarging the unit of $E$ is equivalent to multiplying $k_{\varepsilon}$ and $k_{l}$ to the same extent. Now, we consider the marginal effect of putting $\Delta$ amount of resources into PHIs and NPIs.

● As we do not consider the congestion in the health system in our benchmark scenario, all $\Delta$ is allocated into enlarging testing ability. According to the public data that the cost of COVID-19 testing is around 120 USD per person in the U.S., ^‡ then the funding $\Delta$ can additionally test the following ratio

^‡See for instance https://www.ctvnews.ca/health/coronavirus/u-s-fda-approves-first-at-home-covid-19-test-1.4905616.

● The lockdown policy incurs the cost of subsidising those in quarantine. Some countries have planned paying quarantine subsidy benefits during the lockdown period and the amount varies across countries. For instance, the U.S. has provided checks to over 80 million households and more plans are on the way.^§ We take the widely held number of 66 USD per day (2000 USD per month) as the cost of putting an individual in quarantine, then the funding $\Delta$ could allow the following ratio in lockdown

^§See for instance https://www.businessinsider.com/countries-offering-direct-payments-or-basic-income-in-corona-crisis-2020-4.

Thus, we have

and set $k_{\varepsilon} = 0.3$ and $k_{l} = 0.6$ as the daily resources $E$ is normalised. This assumption is consistent with the observation that few governments can test the majority of residents in a short time, but most implement different degrees of lockdown. This observation is partly owing to the fact that it is too costly or impossible to rapidly expand the capacity of testing in a short time. In addition, we adopt $\varepsilon_{0} = 1/11$ according the public data from Wuhan. Thus, in our benchmark scenario, the testing rate $\varepsilon$ and lockdown ratio $l$ are always lower than $1$.

Value of statistical life. The daily output $w$ works as a nominal amount in measuring the economic cost, and thus we set $w = 1$. Our choice of v.s.l. is $20$ times annual per capita GDP, that is, $\eta = 20\times365$, which is in line with Hall et al. [43] and Alvarez at al. [75].

The adopted parameter values are given in Table 2. In the simulation, the time is measured in days, and thus the rates in our model are then in the unit of "per day".

4.3.   Simulated results

As shown in the Section 4.2, the parameters of our simulation is mainly based on the COVID-19 pandemic. In fact, although there are many epidemiological literature that have studied the cross-regional movement of people in the context of COVID-19, the results of analyzing the impact of human mobility and travel restrictions on the optimal mitigation strategies in the economic context are rare. For example, Kraemer et al. [14], Yabe et al. [15], and Xiong et al. [16] demonstrated the effect of human mobility and control measures on the COVID-19 epidemic in China, Japan and the United States, respectively. But these results rarely involves discussion of economic costs. Bonaccorsi ey al. [32], Hadjidemetriou et al. [33], and Galeazzi et al. [34] discussed the economic and social consequences of travel restrictions in Italy, UK, and France, respectively. However, they didn't study the optimal control strategy. Therefore, this paper bridges human mobility and optimal mitigation through an economic objective.

We first simulate the evolution of contagion and economic costs in both regions when there is no control policy. Figure 4 presents the result of this no-control scenario from day 0 to 365. As no measure interferes with the pandemic dynamics, the infection cases evolve exponentially quickly in the beginning, and peak at around day 40. The overall duration of the pandemic is short because the population of susceptible people vanishes fast. The evolution of no-control scenario is similar with that in Charpentier et al. [23]. Moreover, the pandemic evolution in region B keeps pace with that in region A owing to the human mobility. The vertical dotted lines represent the end of mitigation period ($T = 150$) in Figures 4, 5, E1–E4, and that in Figure F5 ($T = 100$), Figure F6 ($T = 175$).

Then, we investigate the mitigation outcome under the optimal mitigation strategy in the benchmark scenario. The optimal resource allocations for regions A and B are $\alpha_{A}^{\ast} = 0.80$ and $\alpha_{B}^{\ast} = 0.50$, respectively. This implies region A should relatively more rely on the testing policy than region B, mainly because it has accumulated more UI cases. In our simulation, the lockdown policy is lifted in a region when its KI cases is less than $1$.

Infection and death. Panels (a)–(c) of Figure 5 plot the daily current number of UI, KI cases ($I_{i}^{U}$ and $I_{i}^{K}$), cumulative number of infection cases ($CI_{i}$), and cumulative death ($CD_{i}$) in region $i = A, B$, measured in the percentage of population. Both regions successfully contain the pandemic during the mitigation period. As expected, region A has more infection and cases of death than region B. After implementing the mitigation strategies since day $0$, the $I_{A}^{U}$ and $I_{B}^{U}$ both soon peak and then quickly decrease. $CI_{A}$ and $CD_{A}$ are much lower than $CI_{B}$ and $CD_{B}$, indicating that region B remains relatively well contained. This result shows that the pandemic control policy has a more significant effect in the initial stages of the pandemic progression, namely when the fraction of infected is low. Moreover, compared with the no-control scenario, the infections and deaths in region B remain much lower than region A, indicating that the effect of human mobility between two regions is much eased as it has been considered in the control measure.

Economic cost. Panels (d)–(f) of Figure 5 report the cumulative economic loss in region $i = A, B$, including the cumulative lockdown economic cost ($CEC_{i}^{L}$), cumulative death economic cost ($CEC_{i}^{D}$), and cumulative total economic cost ($CEC_{i}^{T}$). The economic loss is measured in the percentage of annual GDP in that region. Both regions have considerable scale of $CEC_{i}^{T}$, around 7.5–15% of annual GDP. Please note that these figures of our benchmark scenario are the economic losses when the pandemic is controlled. The death cost $CEC_{A}^{D}$ is much larger than $CEC_{B}^{D}$, as region A has much more loss of lives. However, a larger proportion of $CEC_{i}^{T}$ results from loss of lockdown policy, $CEC_{i}^{L}$. It is notable that the lockdown cost in region B is higher than region A because it assigns more weight to lockdown. Consequently, the slope of $CEC_{B}^{L}$ is higher than $CEC_{A}^{L}$, and thus region B faces a higher lockdown cost though it lifted its lockdown policy earlier. This explain the observation that $CEC_{B}^{T}$ is higher than $CEC_{A}^{T}$.

In sum, human mobility is the main reason for explaining region B's higher economic cost. Facing the risk caused by the population inflow from region A, the optimal response of region B is to put more people in lockdown. Such a mitigation policy is effective in preventing more infections, but also causes more economic loss.

4.4.   Sensitivity analysis of human mobility

We conduct sensitivity analyses to further explore some variations of our benchmark scenario. In the sensitivity scenarios, we change the value of investigated parameters while the others remain unchanged.

4.4.1.   Mobility rate $\overline{\lambda}$

We investigate the effect of human mobility $\overline{\lambda}$ and consider a range from $0\%/365$ to $100\%/365$. We summarize the contagion and economic costs of both regions in Figure 6. Consistent with our intuition, in general, the mobility rate $\overline{\lambda}$ benefits the pandemic control in region A, but deteriorates that in region B. For instance, the larger $\overline{\lambda}$ increases the $CI_{B}(365)$ and $CD_{B}(365)$, but decreases those in region A. Consequently, we observe a similar pattern of $CEC_{i}^{D}(365)$ in both regions. Regarding the economic cost, the lockdown cost $CEC_{B}^{L}(365)$ first increases with $\overline{\lambda}$ in $[0\%/365, 60\%/365]$, and then fell with $\overline{\lambda} > 60\%/365$. The monotonicity of $CEC_{A}^{L}(365)$ is exactly the opposite. It is intuitive that the mobility rate $\overline{\lambda}$ is for (against) region A (B) in the first stage. When $\overline{\lambda}$ is sufficiently large, the net mobility flow can be in reverse, namely from region B to region A, because region B may experience a later outbreak than region A. Regarding the optimal allocation, Figure G7 illustrates the variation of $\alpha^{*}_{i}$ over $\overline{\lambda}$. In general, $\alpha^{*}_{i}$ is consistent with the pressure of infections and deaths in region $i$. Thus, we observe that generally, $\overline{\lambda}$ reduces $\alpha^{*}_{A}$ and increases $\alpha^{*}_{B}$. Moreover, panel (f) illustrates the overall economic costs of two regions, $OEC(365)$, which reach a maximum around $75\%/365$.

As it is not likely that the planners would choose to enlarge the mobility $\overline{\lambda}$ to be larger than the normal level ($40\%/365$) during the outbreak of pandemic, we are more interested in the left part of each panel in Figure 6. The results of figure (e) partly help to explain why the planners of mildly-affected regions consider travel restrictions and closing borders as a response to the pandemic outbreak in other regions.

4.4.2.   Mobility of UI group

So far, we have considered that the UI group is indistinguishable from the S group and thus the mobility rate of both $S_{i}$ and $I^{U}_{i}$ is $\lambda_{t}$. Now, we further explore the effect of mobility of the UI group. The motive of this investigation comes from the idea that the S and UI populations' accessibility to public transportation and the possibility to travel between regions could become different during the pandemic. On the one hand, the infection groups could have more motivation to move to regions with less pressure of the healthcare system. On the other hand, some travel restriction measures that access the level of the travellers' contagion risk may restrict the mobility of infections. For instance, during the COVID-19 crisis many international flights adopt special boarding regulations, like body temperature measurement and pre-departure COVID-19 testing. Upon arrival, many local regions also require international travellers to have an immediate quarantine, e.g., 14 days. Specifically, the mobility rate of UI individuals can be different from that of S and the upper limits of mobility for these two groups are $\overline{\lambda^{U}}$ and $\overline{\lambda^{S}}$, respectively. In the following, we study the outcomes with the mobility rate $\overline{\lambda^{U}}$ varying in the range between $0\%/365$ and $100\%/365$ and $\overline{\lambda^{S}} = 40\%/365$ unchanged. Figure 7 illustrates that the effects of $\overline{\lambda^{U}}$ on regions A and B are similar to $\overline{\lambda}$, but in a much milder way. Region B has an increase in infections, deaths, and economic costs as $\overline{\lambda^{U}}$ rises, though the growth trend is slight. Similarly, $CI_{A}(365)$ and $CD_{A}(365)$ undergo a mild decreasing trend with $\overline{\lambda^{U}}$. These patterns are also reflected in the economic costs of both regions, as illustrated in panels (c)–(e). Moreover, the panel (f) exhibits an increasing effect of $\overline{\lambda^{U}}$ on $OEC(365)$.

This result is helpful for analysing the policy of restricting the travel of UI individuals by some means.^¶ In general, the mobility of the UI group is harmful for region B without benefitting region A too much. Thus, from the perspective of region B, its planner would favour the case with $\overline{\lambda^{U}} = 0$.

^¶We do not discuss the particular ways to vaguely identify the UI group here. In practice, there are some means to partly restrict the mobility of the likely UI individuals such as body temperature measurement and mobile phone tracking.

5.   Discussion: travel restrictions

The COVID-19 pandemic has sparked an unprecedented shutdown of borders and airlines in many countries to contain the pandemic. So far, around 90% of commercial passenger flights have been affected and more than 130 countries have introduced some form of travel restrictions since the COVID-19 outbreak (Devi [87]). For instance, in February 2020, some countries such as the USA banned travel to and flights from China; in mid-March, European countries in the Schengen area introduced border controls to restrict people's non-essential travel. The restrictions on travel and transportation from high-risk areas helps delay and contain the COVID-19 pandemic (Chinazzi [88]), but hampers the economic communication of personnel in the meanwhile.

In this section, we further discuss the economic impact of travel restrictions. We assume that the policymakers can jointly lower the human mobility rate between two regions by restricting the transportation to $\tau\in\lbrack0, 1]\ $ratio of its normal level. The restricted human mobility rate is

where $\tau$ is the transportation level. Here, $\tau = 0$ is the case when the transportation between two regions is fully cut off, and $\tau = 1$ represents the case with no restrictions on the transportation.

We consider the economic cost of travel restrictions. Personnel communication is critical for the economy production in two regions, and thus we further complicate the economy output in region $i$ as

where the productivity rate $\varphi(\tau)$, common for both regions, increases with $\tau$. Thus, the economic loss due to travel restrictions is

In particular, we consider the transportation level is in the form as follows:

with $\delta > 0$. Thus, $\varphi(1) = 1$ indicates that the economic productivity is not harmed as there are no travel restrictions. Further, $\varphi(0) = e^{-\delta}$ is the economic productivity level when the cross-region transportation is fully cut off. Considering the economic impact of transportation, the economic cost or region $i$ is given as

The term $(1-\varphi(\tau))w(1-l_{i})(S_{i}+I_{i})$ measures the loss of economic output due to travel restrictions.

We parameterize the importance of personal communication on the economic output. Specifically, we estimate $e^{-\delta}$, the economic loss ratio due to full cross-region travel restrictions. It is clear that banning the cross-region personnel communication harms the economic output mainly via reducing the volume of exports and imports (E & Is). The economic output of E & Is are impacted but not terminated by travel restrictions. The economic loss of full travel restrictions is estimated as follows:

We refer to some observations. First, the ratio of E & Is in GDP of many countries fall into the range of 20–60%; e.g., 26.6% (U.S.), 37.8% (China), and 59.5% (Italy).^|| We adopt the middle-level ratio of 33%, which is slightly larger than U.S. but lower than China. Second, we estimated the impact of full travel restrictions on the E & Is by the evidence of the gross transaction volume decrease of China Import and Export Fair (CIEF) during the COVID-19 crisis. CIEF, established in 1957, is the largest foreign-trade expo in China, which was changed to be held online for the first time in June 2020 due to many countries' travel restriction policies. According to the statistics, the gross transaction volume has fallen by around 60%, ^** which we adopt as the maximal loss rate of E & Is caused by complete travel restrictions. Taken together, we estimate that the maximal impact of complete lockdown is around 20% (33%*60%) and thus adopt $e^{-\delta} = 0.8$ in our simulation. It is notable that this ratio corresponds to economic loss under a complete banning of cross-border mobility, which hardly happens in reality. Without particularly specifying the geographic scale and economic structure of the two regions in our model, the estimation represents a mid-level impact of full travel restrictions.

^||Refer to the statistics from https://ourworldindata.org/grapher/trade-as-share-of-gdp?tab=table&time=2010..latest.

^**The disclosed amount of transaction volume of 2020 Spring Fair fall to 135 million USD from 297 million USD in 2019, see the link (in Chinese): https://www.thepaper.cn/newsDetail\_forward\_7999469 and https://www.cantonfair.org.cn/about/deal-stats.

Figure 8 displays the outcomes of optimal mitigation policies under varying transportation level $\tau$ between two regions. In general, $CI_{B}(365)$, $CD_{B}(365)$ and $CEC_{B}^{D}(365)$ are less sensitive to $\tau$ than those in region A. Lower $\tau$ reduces the lockdown cost of region B, however, increases that of region A. Panel (e) shows the cumulative economic cost due to travel restrictions ($CEC_{i}^{TR}(365)$) in both regions, suggesting both regions experiences economic loss due to travel restriction. For region B, Panel (c) and (e) show that it faces a trade-off under the objective of minimising $CEC_{B}^{T}(365)$. On the one hand, the travel restrictions help ease the pandemic pressure in region B as they reduce the inflow of infections and $CEC_{B}^{L}(365)$. On the other hand, the travel restrictions are also harmful for the economic output of region B and pushing up $CEC_{B}^{TR}(365)$. Figure 8 shows that the downside of travel restriction is slightly more pronounced in this setting and that region B favours $\tau = 1$ for the sake of economic cost. From the perspective of region A, it is also optimal to maintain the normal transportation level $\tau = 1$. This is not only because the net human outflow relieves the stress of region A to some extent, but also because travel restrictions bring additional economic costs $CEC_{B}^{TR}(365)$. Therefore, we observe that $CEC_{A}^{T}(365)$ is decreasing with $\tau$. Thus, under our benchmark scenario, both regions agree on that the optimal transportation level should be $\tau = 1$.

6.   Concluding remarks

In this paper, we extend the SIR model to find optimal pandemic control policy considering limited resources and human mobility. We included the constraint of available resources and cross-region human mobility, that is, from a severely-affected region to a mildly-affected one. In our model, the planner of each region optimally allocates the resources into testing and lockdown policies to minimise economic cost. Our model reveals two key insights.

First, we demonstrate that the optimal way of jointly implementing lockdown and testing policies is complicated, depending on many factors. We show how efficient resource allocation depends on the mitigation duration, available resources, and human mobility. Our sensitivity analysis further suggested that the intensive mitigation strategies, with higher daily resource input and shorter duration, are more effective in under-resourced regions.

Second, this study is applicable for drawing policy implications on travel restrictions. Human mobility across regions plays an important role in determining the optimal response strategy as well as the mitigation outcomes. As we write this paper, many countries are in different stages of the COVID-19 pandemic, and some have closed their borders. We further elucidate the impact of human mobility on pandemic transmission and optimal mitigation policy. Furthermore, our model is effective for analysing the economic impact of travel restriction policies in regions with distinct pandemic progressions.

Future research can extend our model to examine other policy-relevant questions. For instance, our model did not consider the time-varying daily resources or time-varying allocation strategies within the mitigation period. At the cost of model complexity, future studies can investigate varying mitigation strategies for different pandemic progression stages. Since our work has only considered the simultaneous move game setting, future works can investigate the extension like the sequential move game setting or the existence of a central planner for coordinating the policies of two regions.

Acknowledgements

This work is financially supported by National Natural Science Foundation of China (Grant No. 72101256), MOE (Ministry of Education in China) Project of Humanities and Social Sciences (Project No.21YJC790016), National Key R & D Program of China (Grant No. 2018YFA0703900), National Natural Science Foundation of China (Grant Nos. 11871309 and 11371226), and National Social Science Fund of China (Grant No. 21AZD028). This work is also supported by Public Computing Cloud Platform, Renmin University of China.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.   Approximation of the human mobility rate

Both the calculation of $R_0$ and the necessary conditions of $(\alpha_{A}^{\ast}, \alpha_{B}^{\ast})$ require $\lambda^{U}_{t}(I_A^K, I_B^K)$ and $\lambda^{S}_{t}(I_A^K, I_B^K)$ to be differentiable with respect to $I_A^K$ and $I_B^K$. Obviously, the denominator item $max(I_{A}^{K}(t), I_{B}^{K}(t)$ is not differentiable when $I_A^K = I_B^K$. Thus, for preparation, we do an approximation of the human mobility rate.

We assumed the human mobility rates for S and UI groups are

where

$M$ is a sufficient large positive constant. If $I_{A}^{K}-I_{B}^{K} > 0$, then $g(I_{A}^{K}, I_{B}^{K}) \approx I_{A}^{K}$. If $I_{A}^{K}-I_{B}^{K} < 0$, then $g(I_{A}^{K}, I_{B}^{K}) \approx I_{B}^{K}$. If $I_{A}^{K}-I_{B}^{K} = 0$, then $g(I_{A}^{K}, I_{B}^{K}) = I_{A}^{K} = I_{B}^{K}$. Thus when $I_{A}^{K}-I_{B}^{K}\neq0$

When $I_{A}^{K}-I_{B}^{K} = 0$, $\lambda_{t} = 0$. Moreover, the second derivative of $\lambda_{t}$ with respect of $I_{A}^{K}, I_{B}^{K}$ is bounded in $[0, \frac{A_S+B_S}{d}] \times [0, \frac{A_S+B_S}{d}]$.

B.   Proof of Theorem 1 and Proposition 1

Proof. First, we verify that

with $\lim_{t\rightarrow \infty} N(t) = N_{0}$ if and only if $\lim _{t\rightarrow \infty} I(t) = 0$. Then the set

is the positive invariant set of Model (2.2).

By Model (2.2) and $N = I_{A}^{U} + I_{A}^{K} + I_{B}^{U} + I_{B}^{K} + S_{A} + R_{A} + S_{B} + R_{B}$, we have

For the homogeneous equations $\dot{N} = -dN$, we have $N = C e^{-dt}$. For the non-homogeneous equations $\dot{N} = -dN+A_{S}+B_{S}$, we have $N = \frac{A_{S}+B_{S} }{d}+\frac{C}{e^{dt}}$, and then $\lim_{t\rightarrow \infty}N = \frac{A_{S}+B_{S}}{d}$. Therefore $\lim_{t\rightarrow \infty}N(t)\le \frac{A_{S}+B_{S}}{d}$. It is easy to find that $\boldsymbol{\Omega}$ is the positive invariant set of Model (2.2).

Then it is easy to find that the Model (2.2) in $\boldsymbol{\Omega}$ has a unique disease free equilibrium (DFE).

where $S_{A}^{0} = \frac{A_{S}}{d+\underline{\lambda}}$, $S_{B}^{0} = \frac {B_{S}+ \frac{\underline{\lambda} A_{S}}{d+\underline{\lambda}} }{(d+\underline{\lambda})d}$.

Let $A_1 = d^U+d+\varepsilon_A+v^U$, $B_1 = d^K+d+v^K$, $A_2 = d^U+d+\varepsilon_B+v^U$ and $B_2 = d^K+d+v^K$. Denote $\lambda_t = \lambda_t(I_A^K, I_B^K)$, then when $I_A^K = I_B^K = 0$, $\lambda_t(I_A^K, I_B^K) = \underline{\lambda}$. We apply the next generation matrix technique as proposed in [66]. New infections come from the process from $S_{A}$ to $I_{A}^{U}$ and $S_{B}$ to $I_{B}^{U}$. $I_{A}^{U}, I_{A}^{K}, I_{B}^{U}, I_{B}^{K}$ are infected compartments. $S_{A}, R_{A}, S_{B}, R_{B}$ are uninfected compartments. Rewrite Model (2.2) as the following system of equations

then

Denote ${\rm{D}}f(x^{0})$ as the Jacobian matrix at the DFE, $x^{0}$. Then it is easy to verify that if $\mathcal{F}(x)$ is set to zero, then all eigenvalues of $\mathrm{D}f(x^{0})$ have negative real parts. According to Lemma 1 of [66],

For region A, we have

and then

For region B, we have

and then

If we consider regions A and B, denote

where

Denote

where

Then

and then

C.   Necessary conditions of $(\alpha_{A}^{\ast}, \alpha_{B}^{\ast})$

Denote $x^T = (I_{A}^{U}, I_{A}^{K}, I_{B}^{U}, I_{B}^{K}, S_{A}, R_{A}, S_{B}, R_{B})$. We can rewrite model (2.2) as $\dot{x} = f(t, x, \alpha_A, \alpha_B)$. Denote $\gamma^T = (1, 1, 1, 1, 0, 0, 0, 0)$, then $I_A+I_B = I_{min} \Leftrightarrow \gamma^T x = I_{min}$. $\hat{T}(x) = \inf \{ t | \gamma^T x = I_{min} \}$. Define

and

where $\eta_i(t)\in\mathbb{R}^8$, $i\in\{A, B\}$. Denote the optimal $(\alpha_{A}, \alpha_{B})$ by $(\alpha_{A}^{\ast}, \alpha_{B}^{\ast})$ and the corresponding $x$ by $x^{\ast}$, then we have the following Theorem.

Theorem 2. The optimal $(\alpha_{A}, \alpha_{B})$ and its corresponding equilibrium $x^{\ast}$ satisfies the following conditions:

Case (1). If $\hat{T}(x^{\ast}) > \widetilde{T}$, then $\forall t \in [0, \widetilde{T}]$,

Case (2). If $\hat{T}(x^{\ast}) < \widetilde{T}$, then $\forall t \in [0, \hat{T}(x^{\ast})]$,

Proof. If $\hat{T}(x^{\ast}) > \widetilde{T}$, given $\alpha_B^{\ast}$, let $\alpha_A^{\ast}$ be the optimal $\alpha_A$. Define

and the corresponding $x(t)$ is $x^{\varepsilon}(t)$. We set $\varepsilon$ as a sufficient small constant which satisfied $\hat{T}(x^{\varepsilon}) > \widetilde{T}$. By Pontryagin's maximum principle, $\forall t \in [0, \widetilde{T}]$,

Given $\alpha_A^{\ast}$, let $\alpha_B^{\ast}$ be the optimal $\alpha_B$. Repeat the above process, then Case (1) is proved. If $\hat{T}(x^{\ast}) < \widetilde{T}$, given $\alpha_B^{\ast}$, let $\alpha_A^{\ast}$ be the optimal $\alpha_A$. Define

and the corresponding $x(t)$ is $x^{\varepsilon}(t)$. We set $\varepsilon$ as a sufficient small constant which satisfied $\hat{T}(x^{\varepsilon}) < \widetilde{T}$. Because $L_A(t, x, \alpha_A, \alpha_B) \geqslant 0$, thus for optimal $\alpha_A^{\ast}$, the terminal condition $t = \inf \{ t | I_{A}(t)+I_{B}(t) = I_{min} \}$ is equivalent to the condition that $\gamma^T x(t) = I_{min}$. By Pontryagin's maximum principle, $\forall t \in [0, \hat{T}(x^{\ast})]$,

Repeat the above process, then Case (2) is proved.

D.   Equilibrium solving process

We followed a two-step procedure for solving the Nash equilibrium in which no player of a region has a profitable deviation given the actions of the other planner. First, we divide $[0, 1]$ evenly into $N = 100$ intervals, and $\alpha_{i}, i = \text{A or B},$ takes values on the boundary of the intervals. That is,

where $\alpha_{i}^p-\alpha_{i}^{p-1} = 1/N, k = 1 \cdots N$.

Next, we solve for the planner's best response function. For every $\alpha_{B}^p \in \{\alpha_{B}^p \}_{p = 0}^N$, the value of the best response function of region A's planner is

which is the region A planner's best response function $\alpha_{A}^{\ast}(\alpha_{B}^p), \alpha_{B}^p \in \{\alpha_{j}^p \}_{p = 0}^N$. Similarly, we can also get the region B planner's best response function

Finally, the Nash equilibrium $(\alpha_{A}^{\ast}, \alpha_{B}^{\ast})$ is the intersection of the two best response functions $\alpha_{A}^{\ast}(\alpha_{B}^p), \alpha_{B}^p \in \{\alpha_{B}^p \}_{p = 0}^N$ and $\alpha_{B}^{\ast}(\alpha_{A}^q), \alpha_{A}^q \in \{\alpha_{A}^q \}_{q = 0}^N$. Our numerical result shows that the equilibrium pair $(\alpha_{A}^{\ast}, \alpha_{B}^{\ast})$ exists and is also unique.

E.   Example illustration in sensitivities analysis of Ω

F.   Example illustration in sensitivities analysis of T

G.   Optimal allocation strategy under different human mobility λ