Spatially explicit models for exploring COVID-19 lockdown strategies.

This article describes two spatially explicit models created to allow experimentation with different societal responses to the COVID-19 pandemic. We outline the work to date on modeling spatially explicit infective diseases and show that there are gaps that remain important to fill. We demonstrate how geographical regions, rather than a single, national approach, are likely to lead to better outcomes for the population. We provide a full account of how our models function, and how they can be used to explore many different aspects of contagion, including: experimenting with different lockdown measures, with connectivity between places, with the tracing of disease clusters, and the use of improved contact tracing and isolation. We provide comprehensive results showing the use of these models in given scenarios, and conclude that explicitly regionalized models for mitigation provide significant advantages over a "one-size-fits-all" approach. We have made our models, and their data, publicly available for others to use in their own locales, with the hope of providing the tools needed for geographers to have a voice during this difficult time.

INTRODUCTION

The COVID‐19 pandemic has brought into sharp relief the complexities of managing a coordinated strategy to minimize human health impacts whilst at the same time minimizing disruption to economic and other social systems. Most responses to date have been applied uniformly, without consideration of the variance in risk or in case numbers that occurs regionally. Our hypothesis here is that geography not only matters, but using it can give a distinct advantage of allowing a tailored response that is more effective at minimizing harmful side‐effects.

In this article, we aim to provide geographers with the tools needed to contribute to the challenge of managing a pandemic through various stages of lockdown, or societal isolation. As our own country (New Zealand) began to consider how to move progressively out of lockdown, it seemed important to us to find a way to bring some geography into epidemiological models, so we responded by experimenting with spatially explicit (regional) models that include information on the connectedness of places, as well as the usual disease parameters such as the reinfection rate, asymptomatic rate, and incubation period. We describe these models below.

Rather than dwelling on the theory of explicitly spatial disease models, we choose instead here to give a more practical account of the software we developed to experiment with: lockdown strategy, regionalization, interaction between regions, disease parameters, test and isolation strategies, and exploring disease clusters. We believe that such models may be useful in other places at the time of writing, and being spatially explicit, they allow other aspects of geography to be included in the planning, thus giving geography (and geographers) a voice at this critical time.

The models and data we have used here have been placed in the public domain so that they can be reused in other locales, and customized to meet the needs experienced by other countries. They can be accessed from the following URLs :

https://southosullivan.com/misc/distributed‐seir‐RC‐web.html

https://southosullivan.com/misc/distributed‐branching‐process‐RC‐web.html

Further details about the modeling environment, and how to customize the models with different geography, are given in the Appendix.

In response to the COVID‐19 pandemic, most countries have employed some kind of lockdown strategy, in an effort to reduce social interactions to only those necessary for essential services such as food delivery and acute health services to continue. These strategies have met with varying degrees of success, such that most countries are faced with the difficult task of needing to reduce lockdown to help revive their economies and restore “normality,” whilst still battling the disease (Kupferschmidt, 2020). It has been remarked many times that it is easier to enter lockdown than it is to leave it. An excellent tool for studying the effectiveness of repression on a per‐country basis is the CovidTrends tracker (https://aatishb.com/covidtrends/), shown in Figure 1. Using this application shows the large number of countries that are now faced with pressure to ease lockdown whilst still struggling to repress the reproduction rate (R 0) below 1.0, hence the need for some very specific modeling of this aspect of the pandemic.

FIGURE 1

Image captured from the CovidTrends website showing the state of contagion in a selection of nations. The plot graphs the number of confirmed cases against the total number of cases over the past week. All the countries shown have had some success with slowing down the pandemic, but only Italy displays a strategy that is headed toward eradication

The modeling of infectious disease is a central focus of mathematical epidemiology, with a variety of mathematical approaches to epidemic modeling explored in some detail over the years (Bailey, 1975). The dominant approach is compartment modeling of the susceptible–infected–recovered (SIR) type. In this article we present one model of this type, and a second based on a stochastic branching process (Allen, 2017; Jacob, 2010), which could be readily extended to an individual‐based model. The former type of model is more suited to established epidemics and modeling their overall evolution, while the latter approach is better suited to the initial stages of an epidemic when the numbers of cases are small, and, we surmise, during the latter stages when an epidemic is more or less under control. Our SIR and branching process models are based on two models developed in New Zealand in the context of a rapidly changing situation (see James, Hendy, Plank, & Steyn, 2020; Plank, Binny, Hendy, Lustig, & James, 2020, respectively).

Our primary focus throughout has been to extend these models to include a spatial dimension. A literature search reveals limited consensus on best practice for spatially explicit SIR models (although see Bartlett, 1957), and “in contrast to the widespread use of other types of models [i.e., SIR], models based on branching processes are seldom used in practice for infectious disease control” (Farrington, Kanaan, & Gay, 2003), so there is scant work on spatially explicit versions of such models. Our branching process model is continuous‐time‐based and readily extensible to a more thoroughly individual‐based approach by the addition of individual attributes (such as demography) to disease cases. True individual‐based models of epidemics remain a rarity (see Willem, Verelst, Bilcke, Hens, & Beutels, 2017). In any case it is not productive to get lost in the semantics of model types. Regardless of the epidemic modeling approach adopted, we know that many regionalized variables will affect both the spread of the virus and our attempts to minimize it—for example, population density, overcrowding of public transport, distance to essential services, age profile, and many more (see Beam‐Dowd et al., 2020). It appears evident to us that smaller spatial units, with limited connections between them, might lead to more success in coming out of lockdown, and managing disease outbreaks that accompany the easing of restrictions on personal movement. We demonstrate later in Section 4 how our results back up these intuitions, while suggesting that much more work is required on this aspect of epidemic control.

CORE EPIDEMIC MODELS AND ADDING IN GEOGRAPHY

In this section we briefly describe two approaches to epidemic modeling and provide some details of the specific instances of these adopted in our work. We also describe previous work on making epidemic models geographical, which is our primary focus in this work. We describe platform‐specific aspects of model implementation, and details of the incorporation of spatial aspects in Section 3.

SIR models

First introduced by Kermack and McKendrick (1927), SIR (also called “removed”) models are compartment models in which the population transitions between different states over time. Initially, the entire population at risk (N) is considered susceptible (S), and the rate of progression toward the infectious (I) then recovered (R) states over time is defined by differential equations governing the transition rates:

The parameters β and γ are respectively rates of infection and recovery per unit of time (in our case per day). These parameters control the dynamics of an epidemic (Kermack & McKendrick, 1927, 1932, 1933). The ratio of the two gives us the basic reproduction number:

Additional insight is gained here if we consider that the expected number of time steps spent in the infectious state is given by 1/γ and since β is the rate of infections per time step, R 0 may be interpreted as the expected number of secondary infections due to each infectious case. Given Equation (4), we can rewrite Equation (2) as:

Thus if

then

meaning that there will be an epidemic until such time as S and/or R are significantly reduced.

The SIR model is readily extended to include additional compartments that might relate to stages of the disease progression or perhaps to different population subgroups. The SIR‐style model we have implemented is based on a five‐compartment model developed for the New Zealand context (see Gatto et al., 2020; James, Hendy, et al., 2020), and influential on the government's decision to go into an early lockdown. This model introduces an exposed (E) compartment and a pre‐symptomatic (P) compartment, as well as the traditional SIR compartments. So technically we could call this an SEPIR model, but since it is a variant of a well‐established class of models, we will stick with the SIR abbreviation here. Future refinements of the compartment model to better fit the observed dynamics of the COVID‐19 pandemic could be substituted into this model as needed at a later date.

In the exposed state cases are not infectious, while in the pre‐symptomatic state they are not yet symptomatic, but are infectious at a lower rate than when they become fully infectious. Using these five compartments allows us to separate the pre‐symptomatic population, which is useful because those who have symptoms are likely to seek medical advice and may well be tested, but those who are pre‐symptomatic will likely not be tested, unless as part of some random background testing. The equations governing this model are extended accordingly, with some minor elaboration to include deaths and hospitalization, which we omit for clarity:

The additional parameters α and β denote, respectively, the per‐time‐step rates of transition from exposed to pre‐symptomatic, and from pre‐symptomatic to infectious, while ε is infectiousness in the pre‐symptomatic relative to the infectious state.

For this model, the relationship between the rate of infection β and R 0 is given by:

which can again be interpreted as spreading the expected number of secondary infections from each case evenly across the expected period of infectiousness.

It is important to note that in many SIR implementations, the model is deterministic and compartment values are real‐valued (i.e., fractional values are allowed). In our implementation we make the model stochastic by treating compartment totals as integers, and the per‐time‐step rates of change as probabilities for drawing random binomial variates from the compartment totals. Thus, for example, the number of cases progressing from exposed (E) to pre‐symptomatic (P) is given in our model, not by αE but by Bin(E, α). This is a departure from mathematical epidemiological theory, where stochastic versions of SIR models are considered Markov processes with a time step sufficiently short that no more than one infection or recovery event may occur in a single time step (see Allen, 2017). Our approach is a good approximation to the Markov model and is a pragmatic choice that allows our model to operate on daily time steps, and reflects a focus on understanding potential variation in outcomes over time and space rather than longer‐term evolution of equilibrium model states.

Branching process and individual‐based models

SIR models are inherently aggregate based as they simulate overall transition rates in a population among disease states. By contrast, branching process and individual‐based models focus attention more closely on individual exposure or infection events. The terminology here is confusing, as many models referred to as branching process models in epidemiology are principally concerned still with the trajectory through time of aggregate counts (S, I, and R), but considered one event at a time (see Allen, 2017; Jacob, 2010).

Zooming in closer still, we can focus on individual cases and their time of occurrence, tracking those individual events over time and aggregating back up to total cases in various stages of disease progression, to provide similar outcome variables as SIR models. This is the approach adopted in our second model. Here, an individual case i is initiated either exogenously or by exposure to a previous case, at exposure time t. A statistical model of the evolution of disease infectiousness—and hence the probability of secondary infections over time from exposure—deploys a probability density function W(t) to generate new exposure times subsequent to an existing case. The number of secondary infections that each case gives rise to is directly parameterized in the model by setting R 0 as a model parameter. In our model, following Plank et al. (2020), the infectiousness progression is represented by a Weibull distribution with peak infectiousness around 5 days after initial exposure, and infectiousness decreasing rapidly after 10 days. In aggregate, this model is effectively the summation over time of many probability density functions each with different initial exposure times t 0.

This approach allows us to model individual case progressions and potentially to incorporate effects such as individual levels of infectivity. This might extend to considering aspects such as an individual case's age, ethnicity, or occupation, or simply stochastic variability between cases. In the present version of the model (reflecting Plank et al., 2020), the significant variation in infectivity is that cases are considered either clinical—in which case they have maximum R c, or subclinical—when infectivity R of individual i is suppressed by some factor k sub. Based on each case's individual R value, in our model, when a case is initialized, it generates Poisson(R) additional exposures set to occur at times t = t + t, where t is drawn from the infectiousness probability density function. This is the maximum number of secondary exposures that case i may give rise to in the absence of any controls. Lockdown and other controls are discussed in Section 3 and the Appendix, where we describe the implementation of our models in more detail.

It is instructive to consider the relationship between the branching process and SIR representations of infectiousness (see Figure 2). The smooth curve shows the distribution Weibull (2.83, 5.67) used in the branching process model, with peak infectivity at around 5 days after exposure. The stepped line shows the transmission coefficients εβ and β associated with cases in the P and I compartments, and the time periods that cases are expected to spend in these states, as a result of our chosen settings of α (= 1), δ (= 0.66667), and γ (= 0.2). The effective infectiousness profile of our SIR model is shown by the bar chart based on simulating 1,000 individual cases progressing through the states in our model. Note that our SIR model is parameterized differently than the model by James, Plank, et al. (2020) on which it is based, to produce a reasonable match between these profiles to assist in model comparisons. The relevant parameters in the SIR model are, however, user configurable to potentially produce very different outcomes, though the versions made available for download have our parameter choices as default. During the SIR model development and validation, we found that parameter combinations that extend the period of infectiousness of cases (i.e., their time spent in the I compartment—easily accomplished by reducing γ) make simulated epidemics much harder to control. However, having high infectiousness persist over prolonged periods does not reflect the situation prevailing in New Zealand during model development, where case identification and isolation has been highly effective, hence our decision to make adjustments in the SIR model parameters relative to the model on which it is based. A more complex SIR model with additional compartments to represent hospitalized and isolated states would allow for a more robust solution of this problem (see e.g., Legrand, Grais, Boelle, Valleron, & Flahault, 2007).

FIGURE 2

Comparison of the infectiousness profiles over time in the two models. The compartment SIR model is shown by the stepped blue line, the branching model by the continuous red line

We recognize that there are many other approaches to epidemic modeling than the two broad classes of model presented here. Of particular interest to the GIS community are agent‐ or individual‐based models (see Hunter, MacNamee, & Kelleher, 2018; Willem et al., 2017), which may incorporate richly detailed representations of economic, social, and other activities in addition to centrally representing the progress of an epidemic through the population. Our focus on epidemic control has led to our emphasis here on more “epidemic‐centered” approaches. The specifics of the development of our models also owe a great deal to the evolving crisis response in New Zealand in March–April 2020, when the two models we have highlighted here (James, Hendy et al., 2020; Plank et al., 2020) were deployed in advising the government on pandemic decision‐making.

Adding geography to epidemic models

If the population and the disease were distributed evenly through a region, there would be no value in a geographical approach. However, since the distribution of both vary a great deal with place, we can leverage this to change the lockdown strategy so that areas free from infection can also be free of restrictions as soon as it is deemed safe to do so. Lockdown strategies harm the economy and reduce access to healthcare services for other medical conditions, both of which are also causes of suffering and death, so minimizing this harm where it can safely be done is vital.

To facilitate a geographical approach to lockdown, we need to decompose the affected area into regions, and then limit non‐essential travel between them. Then we can tailor the management strategy to a more specific risk: a local mean for the infective population will be much more accurate than a global one, and thus will allow a more targeted approach.

Given that the mobility of the susceptible population and their interaction with infected people is vital for any infectious disease to spread, it is perhaps surprising that fundamental geographic processes have largely been excluded from SIR models and their variants (Allen, 1994; James, Plank et al., 2020; Kermack & McKendrick, 1927, 1932, 1933; Plank et al., 2020). By convention, these compartment models assume that the population is homogenous, and that the risk of a susceptible individual becoming infectious is a function of the probability of being a contact and the effective R at a given time. Elaborations of SIR models to handle “meta‐populations” categorized by age, ethnicity, or other characteristics commonly become mathematically complex rather quickly.

While SIR models are ubiquitous in the epidemic modeling literature, there appears to be no consistency in how best to develop spatially explicit epidemic models. For example, cellular automata (CA) models have been used to extend some of the grid‐based epidemic analyses proposed by Bartlett (1957) and Hagerstrand's (1953/1967) diffusion of innovation theory to infectious diseases. Munshi, Roy, and Balasubramanian (2020) used CA models to model the spread of COVID‐19 in New York City and to demonstrate how various control measures can potentially mitigate the epidemic. They developed a spatially explicit susceptible, latent, quarantine, infectious, isolated, recovered (CA‐SLQIIR) model that factored in measures such as social distancing and population mobility. Results suggested that prolonged (e.g., 180‐day) periods of lockdown, localized rather than regionalized mobility, and ensuring 95% of the population were tested for COVID‐19 would reduce the number of cases and fatalities substantially.

One potential limitation of grid‐based models such as Bartlett's or Munshi et al.’s is that modeling hierarchical connectivities, such as high‐volume transportation routes between cities, or arterial roads within cities, is challenging. However, network‐based approaches can overcome this limitation. Cliff, Haggert, and Versey (1981) proposed network‐based models of disease diffusion in an attempt to add spatial structure into disease modeling, focusing particularly on measles in island settings. Gould (1993), Gog et al. (2014), and Kissler et al. (2019) describe approaches to pandemic modeling based on gravitationally weighted diffusion over a transportation network, but with no local epidemic mechanism per se. Saito et al. (2013) used two sets of SIR models, coupled with a parameter representing the migration between two regions, to model the spread of the H1/N1 pandemic in Japan during 2009. Bertuzzo, Casagrandi, Gatto, Rodriguez‐Iturbe, and Rinaldo (2010) tested different network configurations representing hydrological connections between communities and implemented SIR models in which the spread of cholera was modeled through a biased random‐walk process on the network.

Stochastic spatial SIR models are another form of spatially explicit model that enable individuals to be “tracked” through the four stages of disease over discrete time and across discrete space. Brown, Porter, Oleson, and Hinman (2018) used an approximate Bayesian computing (ABC) approach to stochastic spatial susceptible, exposed, infectious, recovered (SEIR) modeling to simulate the spread of Chikungunya, a virus transmitted by mosquitoes with a typical incubation/latent period of 3–7 days. In their spatial SIR models, the authors use two approaches to calculate the transition probabilities representing disease latency (E to I) and infectious period (I to R), respectively. First, they use an exponentially distributed approach, which is very common in Bayesian SIR models, and then, as an alternative, a path‐specific SEIR (PS‐SEIR) approach that models individual paths through exposed and infectious compartments. The strength of the PS‐SEIR method is that the transitions need not be exponential.

Balcan et al. (2010) developed the global epidemic and mobility (GLEaM) computational model to simulate global pandemics. GLEaM is one of the few examples that leverage population mobility and sociodemographic data at the sub‐national level, combined with disease modeling parameters to simulate the spatial diffusion of infectious disease on a global scale. For population mobility they used air travel networks with origin–destination matrices of commuting patterns to ensure comparability between and within countries for their models. Their spatial disease model uses a susceptible–latent–infectious–recovered (SLIR) compartment approach, which is estimated for each geographic census area, and can be modified to account for age structure or seasonal variation, for example. The progression between states in the SLIR model is governed by transition probabilities similar to those for conventional SIR models, with additional probabilities included to account for the inclusion of national or international air travel. Thus, for each geographic census area, the disease is transmitted between adjacent areas when people commute locally, whereas the air travel component dictates the rate of spread nationally and internationally. The GLEaM models are also available as an online tool—GLEAMviz (http://www.gleamviz.org/model/), which allows users to re‐configure the models and simulate disease dispersion.

Klepac, Kissler, and Gog (2018) describe the implementation of an outbreak model using data obtained from a mobile phone app over a 1‐month period in the UK, as part of a BBC documentary to mark the 100th anniversary of the 1918 influenza epidemic. They captured user characteristics (age, gender, max‐distance‐traveled) and location coordinates for each hour during a 24‐h period. For most analytical purposes they restricted the sample to distances traveled below 100 km, however they used larger commuting distances (>200 km, with population density> 100,000/km2), with a small trickle effect to spread their contagion further afield. Age‐specific transmission rates were calculated and the authors combined geographic boundary data from Scotland, England, Wales, and Northern Ireland to create their geographical base—known as “patches.” They modeled within‐patch variations as a discrete‐time SIR‐style model with a realistic infectious profile, based on the age‐structured transmission rates. In addition, they used a gravity‐style model, with a stochastic implementation using the real movement data to model between‐patch patterns.

Initially, Klepac et al.'s model simulates a classic contagion pattern, in which the spread of an influenza‐like illness (ILI) encompasses the country in concentric waves over 14+ weeks, although the impact of London's long‐distance connectivity becomes obvious reasonably early in the epidemic. A subsequent model emulated the effect an increased rate of hand washing had on the dispersion of disease. Whereas the initial model took around 7 weeks to cover the length of Britain, washing hands reduced the rate of spread by about 4–5 weeks. While the control measures did not eliminate the pandemic, there were an estimated 10 million fewer cases.

The current COVID‐19 pandemic has highlighted a need to extend the SIR models to incorporate different levels of disease severity and disease management/recovery. To explore the impact of control measures on the spread of COVID‐19 through communities in Italy, Gatto et al. (2020) therefore proposed a spatially explicit model that comprised susceptible (S), exposed (E), pre‐symptomatic (P), infected with heavy symptoms (I), asymptomatic/mildly symptomatic (A), hospitalized (H), quarantined at home (Q), recovered (R), and deceased (D) individuals. Geographic relationships between communities were represented as a network of 107 nodes (i.e., major cities/towns) using a bi‐directional connectivity matrix, facilitating local transmission among adjacent nodes. A particular strength of this study is its use of publicly available information relating to human mobility, confirmed COVID‐19 diffusion through Italy between February 21 and March 25, 2020, and adjustments in the models to reflect four levels of government‐imposed control to mitigate transmission. Overall, while their models tended to be underestimates of the actual hospitalizations or mortality, their results were able to demonstrate that the control measures that the sequential restrictions imposed on population mobility reduced transmission of COVID‐19 by 45%.

Despite their wide utility, there appears to be no agreement in how SIR models should be implemented in a spatially explicit context. The brief review of models here demonstrates a lack of consistency in configurations of space and spatial structure, the compartmentalization of the epidemic phases, and measuring population mobility and rates of disease transmission between compartments. This concurs with criticisms by Getz, Salter, and Mgbara (2019), who critique the adequacy of SIR models for their (in)ability to model diffusion of diseases with spatial structure. Most of the currently available applications only model specific locations, and none allow explicit and detailed experimentation with lockdown strategies, hence our efforts here to fill this gap.

MODELING REGIONAL EPIDEMIC CONTROL

The model interface

We implemented two epidemic models using the NetLogo platform (Tisue & Wilensky, 2004; Wilensky, 1999). More details on the platform and the models are provided in the Appendix, but we provide a brief overview of the model interface here to help readers interpret the figures in the sections that follow. Figure 3 shows the user dashboard from the SIR version of the model; the branching model has a similar, albeit simpler, user interface. We have annotated the screenshot with boxes and labels to show the position and functions of different groups of controls. Within these groups:

FIGURE 3

The user dashboard for the web version of our geographically explicit SIR model. The dashboard has been annotated with outlines and labels to show the position and functions of different groups of controls. The branching model interface is largely similar. See text for more details

Green boxes indicate parameters that the user can set or change.

Yellow boxes show various working values and outcomes.

Blue buttons perform actions.

The map at center shows the chosen geography; the size of the population in a region is denoted by a proportional symbol (a circle). A square box surrounds each population circle, the border of which indicates the lockdown level, from most to least severe: red, orange, yellow, and green (these will be explained in more detail in Section 3.4). The color inside the circle indicates the number of recovered cases, and changes through pink and dark red to black as the numbers increase.

Setting up the model

Details of how to configure the epidemic properties, the starting conditions, the regional geography, and associated population and mobility attributes are provided in the Appendix (Sections A.2, A.3, A.4). These are important details to understand if you intend to use the model, but here we concentrate on demonstrating how a regional approach to managing an epidemic can be supported.

Regionalizing epidemics: Locales and connectivity

Central to our approach is that we implement a collection of local core epidemic models (whether SIR or branching process), across a collection of “locales,” with each locale represented as a separate unit in the model. Thus, in the SIR variant, locales hold a collection of attributes representing the total population in each compartment and manage the updating of these counts over time, with respect to a localized value of R 0, which in turn depends on a local alert‐level at each locale. These local SIR models operate exactly as a global whole‐system model might, with the exception of the geographical connectivity discussed below.

In the branching process variant, each locale maintains a collection of clinical and subclinical cases (each themselves an agent). The cases each maintain a list of initial and secondary infection exposure times used to determine whether a case has progressed to isolation, hospitalization, and so on. These in turn may affect the probability of an exposure event in that case’s list of exposures, actually resulting in a new case, when the time comes. When a case is initialized, its list of exposure times is populated with the maximum possible number of secondary infections it might cause. When an exposure time is drawn from the list, depending on the current alert‐level of the locale, and whether or not the case has been isolated in its disease progression, that infection event might not occur (for example, if the locale in which the case is located is in a high state of lockdown, the exposure will occur with relatively low probability). The possibility that a new exposure will occur not in the current locale of the case, but in another one elsewhere in the system is dependent on the connectivity among locales, to which we now turn.

A collection of entirely independent local epidemic models would not be particularly interesting. In our system, locales are connected into a network, by directed edges with associated edge weights w. Weights can be pre‐calculated and read in, or are calculated according to the common gravity formulation, such that:

where p and p are the populations of the two locales, and d is the distance between them. These weights are rescaled so that the total weight of the outbound links from each locale sums to 1. The w values are used in combination with a list of = {f} values to control the degree of interaction between locales in each time step. In the SIR variant of the model, outbound flows from each locale decrease and inbound flows increase the effective population of the pre‐symptomatic compartment in the locale‐specific model, according to:

where denotes the lower (i.e., more restricted) of the alert‐level flow levels between the two locales, and Pois(∙) is a Poisson‐distributed random deviate. In the branching process variant, the local alert‐level control f is the probability that a new infection arising from a case in this locale will give rise to a new case non‐locally (i.e., in another locale). When a non‐local infection occurs, the target locale is chosen randomly weighted by the relative outbound weights of the connections to other locales w. These approaches are not identical, but produce qualitatively similar behavior across the two model variants, and do not require populating with detailed empirical data about actual rates of movement between locations. In future work this aspect of the model can be refined to allow for initializing connections with average number of daily trips between locales, and the mechanics of how this is implemented in each model variant will be adjusted accordingly.

In New Zealand, health services are provided by 20 district health boards (DHBs), which are overseen by the national Ministry of Health. Since each DHB is tasked with managing the resources for their district, from an administrative perspective, one potential partitioning of New Zealand for regional management of COVID‐19 is by DHB. For the sub‐national government, New Zealand is also divided into 67 territorial authorities (TAs). These TAs serve as an alternative way of organizing regional management of COVID‐19, and are a particularly attractive solution because they correspond closely to the areas managed by local civil defense groups, which are responsible for emergency management. In most cases the TAs cover smaller areas than the DHBs, but not always. For example, the Auckland supercity region comprises a single TA, though it is served by three different DHBs: Waitemata, Auckland, and Counties Manukau. GIS data representing the geographic boundaries for the DHBs and TAs is supplied by Statistics NZ Tatauranga Aotearoa (https://www.stats.govt.nz/topics/population) (see Figure 4).

FIGURE 4

New Zealand DHBs (left) and TAs (right), used to evaluate our lockdown strategies. The populations in each region are represented by the size of the white circles. The three overlapping large circles in the left image represent the three DHBs serving Auckland. The blue arrows show connectivity, calculated as described in the text. Both depictions have the control set to 300 km

Simulating different lockdown strategies

On March 21, 2020 the New Zealand government introduced a four‐level alert system to communicate different degrees of lockdown on the mobility of people in the country (https://web.archive.org/web/20200505213339/https://covid19.govt.nz/alert‐system/covid‐19‐alert‐system/). Level 1 indicates that the disease is contained in New Zealand but uncontrolled overseas. At level 2 it is still contained but risk of community transmission is present. There is “high risk” that the disease is not contained at level 3 (community transmission might be occurring), and it is “likely” that it is not contained at level 4 (community transmission is occurring). Progressive measures of lockdown are implemented at each level, starting with social distancing measures in public places and leading to the closing of most businesses and schools, and strong restrictions on movement outside of home bubbles and local areas. Some countries (such as France) have implemented a three‐level alert system, but the principles are common—reduce infection by progressively removing or redefining the social interactions between people (Chinazzi et al., 2020; Kraemer et al., 2020).

To model these alert‐levels, we must make assumptions about how these progressive restrictions of social interactions work to modify R 0. Since this is the key epidemic parameter common to both models, this is done by setting up vectors of R modifiers, c = {c 1, c 2, c 3, c 4} for four different scenarios for each lockdown:

[1 0.8 0.6 0.36]

[1 0.72 0.52 0.32]

[1 0.64 0.44 0.28]

[1 0.7 0.4 0.16]

The first three scenarios are estimates taken from the work of Plank et al. (2020). The “other” scenario represents our ongoing “model experiential” estimate of what best reflects the observed performance of the lockdown levels in New Zealand, which were particularly restrictive. The values are ordered from level 1 to level 4 lockdown (which gives the biggest reduction in R). Thus, if R 0 is set at 2.5, the modified R value for level 4 in an optimistic lockdown strategy is c 4 R 0 (i.e., 0.28 × 2.5 = 0.7), and this value displays in the (ective) box, at model initialization (subsequently this value may change as local lockdown levels change to reflect a population‐weighted average of local R values). Different scenarios are needed because at this stage of the pandemic we are still unsure what effect different lockdown rules will ultimately have on transmission: it depends how much social interaction still persists in delivery of priority services, and how diligently the population adheres to the rules.

R values in these scenarios should be modified as relevant data become available, for example the Google COVID‐19 Community Mobility Reports (https://www.google.com/covid19/mobility/), or by carefully observing how R behaves under lockdown. The user can choose a scenario via the drop‐down.

To simulate changing lockdown levels, we need rules that activate different alert‐levels when various conditions are met. We have chosen to use threshold values for positive test results set in the control (one for each lockdown level). The testing component of the model is straightforward (and probably optimistic in many settings). The number of tests conducted is calculated from a baseline overall rate of testing (in model runs reported here set to 0.001) and a higher probability of testing of infectious cases. In the SIR variant, the probability that cases in the infectious stage will be tested is set by the parameter, while in the branching process model it is controlled by the parameter. After infection, a case determines a random time to isolation from the sum of two random number draws: t 1 = Gamma(5.5, 0.95) and t 2 = Exp(t) (again based on Plank et al., 2020). The first of these, t 1, is the time to onset of symptoms from exposure in clinical cases, while the second is the further time taken for detection (i.e., testing and isolation) of the case. In the day during which the time since exposure reaches t 1 + t 2, a test is assumed to have occurred (results reported in Section 4.3 show how important this setting can be and reflect the significance attached to effective testing and contact tracing). With some probability in both models a false negative test is possible, and occurs with probability . Calculations of the positive test rate are determined for every locale and can be aggregated up to overall global values.

Depending on the , the positive test rate triggers changes to the alert‐level based on the trigger levels. A variety of alert policies have been implemented:

“static” results in all locales have their alert‐level set to that specified by the setting;

“global‐max” means that the highest locally calculated positive rate triggers the corresponding alert‐level for all locales. This is a conservative approach where a single outbreak in only one locale will lead to increasing the lockdown level across all locales;

“global‐mean” sets the alert‐level in all locales based on the population‐weighted mean positive test rate across the whole system. This is a lax policy that may reduce alert‐levels in locales with current outbreaks, because there are few outbreaks elsewhere;

the “local” policy causes the alert‐level setting in each locale to be changed based on the positive test rate specific to the locale;

a “scripted” policy causes alert‐levels globally to be changed at the specified number of days after model initialization listed in the input textbox. This enables approximate matching of actual policies enacted in particular places to be simulated.

Figure 5 contrasts the results of a local and global lockdown strategy. Note how the epidemic outcomes (total infections, recovered, and so on) are very similar, but the lockdown levels are very different. Section 4.1 provides more detail on these differences and particularly on the lighter restrictions resulting from a regional approach.

FIGURE 5

Two approaches to lockdown. Both screenshots show progress after 170 days from the start of a 2,000‐case outbreak, using population data from New Zealand's TAs. The top image shows a static, global lockdown strategy using the highest level of lockdown. The bottom image shows an adaptive, local strategy where lockdown changes in response to the local infection count. In terms of containment, both approaches perform similarly, and the epidemic outcomes are similar. But notice how much more of the country has moved to lower lockdown levels in the bottom image, shown by the color of the square outline around the population centers (red = level 4, orange = level 3, yellow = level 2, green = level 1)

The frequency and onset time of changes between alert‐levels are controlled by the and settings. The former dictates the first day after model initialization when changing alert‐levels is even considered, while the latter controls both the frequency with which alert‐levels may be changed, and the time period over which test results are aggregated to trigger potential changes. It is assumed that it would be bad policy to change alert‐levels on a daily basis, and that weekly or less frequent updating is likely. When alert‐levels change they are “sticky” downwards, meaning that they only move down one level at a time, even if testing suggests a lower level might be viable. Thus, a locale at level 4 will move to level 3 even if testing suggests level 2 would be okay, and another cycle of duration given by will be required before the level might move down to level 2. Moving up levels is “responsive”: if testing dictates level 4 is appropriate, the alert‐level moves immediately to that setting. Setting alert‐levels is one thing, having the population respond is another, and this is controlled by the setting. The actual control setting of locales is based on a moving average window over the previous number of days specified by this setting. This means that when the alert‐level is changed, it will take this length of time for the lockdown change to be reflected in the local control levels. When we look at the data from implemented lockdowns, there seems to be a smoothing that occurs during transitions, hence the potential need for this control. One possible reason when moving to higher levels is that within isolation bubbles, new infections can still happen if one or more members carry the virus, so R is reduced more slowly. We have seen with horror how disastrous this can be in retirement communities and nursing homes. When moving down lockdown levels, there may similarly be a reticence to exploit all the new freedoms immediately. When set to zero, changes in R happen instantaneously.

Measuring disruption caused by lockdown

Like it or not, countries need to balance their epidemic response with other needs relating to the economy and the social interactions between people. So, as well as the usual outcome measures relating to infection, hospitalization, and death, we can also add in measures that describe the amount of time spent in a given lockdown level, and also the number of transitions that occur between different lockdown levels, since these are both indicators of the likely social and economic disruption. For example, we might decide to trade additional freedom of movement for a willingness to move back up to harsher lockdowns if the disease flares up again. The and displays show the number of people at a given lockdown level and the number of regions at a given lockdown level. The effective mean ( ) score also provides a sense of the relative freedom that the population experiences. This score is a population‐weighted sum of the R values for every step in the simulation. For the top image, this value remains fixed at 0.4 throughout the simulation, because the R‐modifier under this scenario for level 4 is 0.16 (0.16 * 2.5 = 0.4). For the bottom image, it changes as different regions move to lower lockdown levels as the outbreak becomes locally contained. The value of 1.228 shown in the bottom image reflects the fact that time is being spent at lower lockdown levels, with more activities restored.

Of course, should a flare‐up occur because of inter‐region movement, the lockdown level may need to respond by moving higher again, and this happens automatically in the simulation. This of course offers an increased risk of additional new cases, but appears to be more than offset by the more targeted use of local rates, since a global mean may be quite out of step with local conditions. The score gives a convenient measure of how many times (on average) a region changes level. A local adaptive strategy will usually lead to more changes, but also to more time at lower levels. In the adaptive strategy used in Figure 5, the average number of alert‐level changes per region is 2.576 as opposed to 0 for the static strategy. This disruption, and the need to manage communicating and enforcing different alert‐levels in different regions, are the downside of a locally adaptive strategy.

A large population center is problematic since it cannot easily be divided into regions that can effectively be isolated. Auckland, with its population of 1.5 million, is often the last TA region to clear in a local simulation. A different approach may be needed here (such as isolating particularly susceptible sub‐communities such as retirement villages inside an internal “bubble”), but is beyond our current scope.

New Zealand was fortunate to have begun lockdown early, and to have used a very stringent lockdown policy (described above). If we had waited until we had 5,000 cases in progress (that equates to about 1,800 known cases, assuming perfect testing), then the same lockdown strategy would have led to a resurgence as soon as we moved to level 2, as shown in Figure 6.

FIGURE 6

A simulation of the current New Zealand lockdown strategy of 5 weeks at level 4, followed by 3 weeks at level 3, and then down to level 2. In this case the strategy was applied later in the progression of the epidemic, assuming 5,000 infections in progress before beginning, then played out for an additional 8 weeks. The resulting surge is stronger than the initial outbreak.

More detailed results and findings pertaining to lockdown strategies and their effects are presented in Section 4.

DETAILED EXPERIMENTAL RESULTS

Lifting lockdown and controlling epidemic spread

A primary focus of this model is to explore the feasibility of localized lockdowns as a way to maintain control over epidemic spread while reducing the impact of disruptive lockdowns on the population.

Figure 7 shows a time series of the outworking of different lockdown levels under different alert‐policy settings for the total population. These results are from SIR variant model runs in the New Zealand setting, initialized with 2,000 cases, simulated 30 times, with regionalization by DHBs (n = 20). Under a strict “static” policy, the whole population remains in level 4 lockdown until eradication, and this takes almost an entire year. The “global‐max” policy yields a similar result, so that even though the option to reduce lockdowns is in place in these model runs, positive testing rates did not indicate such changes were desirable. The “global‐mean” policy not only takes longer (sometimes much longer) to eradicate the epidemic, it also sees many changes in alert‐levels. By contrast, under the “local” policy, many people are able to move into lower lockdown levels relatively early and there is general progress toward widespread lowering of lockdown restrictions while retaining similar or better levels of control over the epidemic. The results in this figure and the next one are based on “realistic” alert‐level control settings [1, 0.72, 0.52, 0.32]. These settings are probably pessimistic relative to apparent outcomes in New Zealand, but may be more realistic in other contexts.

FIGURE 7

Population counts in different lockdown levels over time, as a 2,000‐case outbreak is brought under control under different alert‐policy settings. The graphs show 30 different realizations, to show the degree of variation experienced from one run to the next with no parameter changes, due to stochastic variation

Figure 8 shows the temporal changes to the daily number of infections for the same set of runs shown in Figure 7, confirming that the “local” policy is similar in effectiveness to both the “static” and “global‐max” policies. Again, the results show 30 different realizations of the same settings, to avoid false results caused by minor stochastic variations. It is important here to note that none of the policies differ greatly in terms of control of the epidemic. The long‐duration tail on the “global‐mean” plot sees the number of new cases per day persist at levels as low as one or two cases, with no new major outbreak ever establishing across the 30 simulated runs.

FIGURE 8

Time series of daily reported new infections for the same set of model runs shown in Figure 7. The mean result is shown by the thin black line, the variance by the shaded gray area

A possibly more useful summary of these results is seen in Figure 9, which compares times to eradication across the same 30 realizations under different alert‐policy settings and different regionalizations, either DHBs (n = 20) or TAs (n = 67). Based on these results, it is not possible to say anything with certainty about the relative effectiveness of managing lockdown at different regional scales. However, it is again clear that a policy based on regional control (“local”) can perform on a par with the more conservative strategies (“static” and “global‐max”), albeit perhaps with more variance, but also (as seen in Figure 7) while allowing more freedom of movement.

FIGURE 9

Boxplots comparing times to eradicate a 2,000‐case COVID‐19 outbreak in the New Zealand context under a range of lockdown (alert‐level) policies across 30 different stochastic realizations of the SEIR model, and for two different regionalizations, by DHBs (n = 20) and TAs (n = 67)

It is important to be clear that the results of various experiments run to investigate trade‐offs between lockdown and epidemic control are highly complex and demand considerably more investigation. We consider the more optimistic settings, specifically the setting [1, 0.7, 0.4, 0.16] to probably be more representative of New Zealand’s lockdown regime. Using such settings, the differences in the success of different alert‐level policies are less marked, with even the “global‐mean” policy working relatively well. Of course, under those settings both levels 3 and 4 imply stringent controls on social interaction, and there is likely to be less scope for local departures from the general lockdown level across the system globally. Similar results were obtained by running similar experiments on the branching process variant of the model.

Modeling disease clusters

Since the branching version of the model represents individual cases rather than totals for the disease stage compartments they are in, we can explore the clusters that emerge as the virus moves through the population and from region to region.

Figure 10 shows disease clusters after 2 weeks of (ineffective) control at alert‐level 2 of an outbreak initialized with 100 cases. The left‐hand panel shows the clusters in their geographical context and emphasizes instances of inter‐regional transmission that have occurred. The right‐hand panel shows clusters detached from their geographical locations to emphasize the distribution of cluster sizes and structure.

FIGURE 10

Clusters that form during an outbreak, viewed in two different ways in the model interface

The views shown in Figure 10 are experimental but provide insight into the details of disease transmission that the branching process model simulates. For a more complete picture of disease clusters, we can examine the logs the model produces, to extract cluster information and reorganize it into a graph. Figure 11 shows how a single‐origin case (1923) leads to a cluster of over 200 additional cases in one simulation using the TA regions. The color of each node indicates the region where a case resides, so this cluster starts in the blue region but spreads over time to four other regions: green, orange, purple, and red. The spatial interaction in the model leads to the disease bouncing back and forth between regions: note how the disease jumps to the green region and back to the blue region at a later time, and then back to green once again (cases 2,046, 3,024, and 5,571, for example).

FIGURE 11

Graphing the paths of contagion by post‐processing log data into a D3.js radial tidy tree (https://observablehq.com/@d3/radial‐tidy‐tree). The graph shows the cluster originating from case number 1923, and how it jumps to new regions (see text for more details)

The ability to model the path of infection from any individual back to the initial infection helps us to understand the kinds of infection chains that form, and how different lockdown and suppression measures affect them.

This feature of the branching process variant of the epidemic model demonstrates the potential for more detailed simulation of contact tracing and isolation as methods of epidemic control. In the present models, the epidemic R parameter is directly modified whether system‐wide or locally as a proxy for a wide range of non‐medical interventions, such as restrictions on movement and large gatherings, working from home, school closures, masking, public education about handwashing and hygiene, and so on. Contact tracing as a public health measure has the potential to provide more targeted control that is less disruptive. With the branching process model explicitly generating an emergent sequence of connected cases, adding case progress status variables related to whether or not a case has been tested and isolated, and whether or not as‐yet‐undetected related cases have been detected, could be modeled and related to the effectiveness of the testing and contact‐tracing regime.

In the current branching process model we do this in a limited way, based on the model setting, which determines how long the delay is from infection to when a clinical case is isolated and its probability of causing subsequent exposures is reduced. We show the importance of this in the next section, but note that for now the model does not include more details of contact tracing. The current implementation thus crudely represents the effect that testing has on reducing R by slowing further spread from a known case, but does not include the possibility that a newly detected case might enable detection and isolation of previously unknown predecessor cases in disease clusters.

The impact of improved testing and contact tracing

As we have just noted, the branching process model crudely simulates the effects of contact tracing and case isolation. The graph in Figure 12 demonstrates the effect that the control has on reducing new outbreaks at various settings. This control effectively allows detected cases to be isolated after a given time delay, so they cause no further infection. These results were obtained using the branching model with the setting for the and, as above, they use 30 different simulations with the same parameter settings. We assume that there is capacity to perform contact tracing at the rate it is required by new cases being generated, but this too could be simulated with minor extensions. When the time to detection and isolation is low, simulated epidemics are more rapidly brought under control, signified by the rapid transition of the whole population to alert‐level 1 with very few recurrent outbreaks (top‐left panel). By contrast, in the top‐right and bottom‐right panels, with long detection delays, even in the best‐case scenarios it takes longer for the whole population to move into alert‐level 1, and there are numerous recurrent outbreaks, requiring frequent returns to higher alert‐levels. Even this simple example shows the importance of effective testing and contact tracing to any system of localized epidemic control.

FIGURE 12

The effects of early testing and isolation on the effectiveness of different lockdown strategies, illustrated using 30 simulations of DHB‐level data using 300 cases and a variety of time‐to‐detection measures, from 0 to 10 days

Simulating an insecure border

The arrival of exogenous cases is problematic, and a challenge for many countries and regions. But it is likely a fact of life, so it must be modeled. Figure 13 shows a simulation of New Zealand DHBs beginning with only 100 cases, but with 5 new cases arriving at random locations in the model each day. Some of these new cases will likely infect others, leading to a simulation that never terminates but oscillates around a mean daily count of around 10 new infections into the future. Note that at this point even a strict static lockdown cannot clear the infection. The control allows the number of new extraneous cases per day to form part of the simulation.

FIGURE 13

Introducing even a small number of external cases per day will usually lead to an ongoing epidemic, with several new cases reported each day into the foreseeable future

CONCLUSIONS

We have developed two spatially explicit models for simulating and addressing the evolution of the COVID‐19 pandemic on a regional basis. Results show that the application of a regionally varying series of lockdown policies is likely to be just as effective at minimizing contagion, while offering advantages of less restrictive rules for part of the population.

Trading off some additional complexity in the mitigation for a more liberal lockdown may be important. Bear in mind that to clear a large number of infections (thousands or tens of thousands) requires considerable time; our models show at least 3–4 months being needed in a high‐lockdown state for scenarios with high case numbers (see Appendix, Section  A.6) before it is safe to reduce the alert‐level and instead rely on contact tracing and isolation. This assumes that the lockdown strategy is adhered to for the duration. During this time, the population may suffer from loss of income, reduced access to healthcare and education, and the stresses of being apart from family and other social networks. All these have negative effects on wellbeing, even if we set aside all economic concerns.

At the time of writing, New Zealand sits on the cusp of an (initial) eradication of COVID‐19, with only 3 cases reported in the last 2 weeks. We do not yet know if a regional lockdown policy will be needed here, it may well be that the 8 weeks of progressive lockdown we have experienced will be sufficient, given that our case numbers have been small. If not, we may need to use this model in anger. Other places are not so fortunate, and have a pressing need now.

Our main findings are summarized as follows:

Explicitly spatial models are needed to address the complexities of managing the response to pandemics—they can offer similar performance in terms of case suppression, but with the potential for significantly less disruption to other important activities, as demonstrated in Section 4.

Our results confirm that eradicating the virus is challenging, and we see many simulations where persistent small outbreaks oscillate between regions for many months, potentially causing a lot of social and economic upheaval.

Results indicate that a regional approach could lead to less illness, because one can use a lower threshold to respond locally than makes sense globally. Although more research would be required to tune lockdown strategies to make such outcomes likely.

Using spatially explicit models of lockdown also facilitates the inclusion of other important geographical dimensions into models (see discussion in Section 5.2). We have shown how the connectivity between regions can be used to simulate likely spread of the disease and how this connectivity can be constrained to simulate various aspects of lockdown.

At lower lockdown levels, just a small number of regular external cases, say from international travel, can cause new clusters to form that need to be contained. Unless we operate regionally, the response to address this will mean imposing strict lockdown measures country‐wide, with associated negative impacts on social, economic, and medical systems. We believe this can be avoided.

Large metropolitan areas present a real challenge to manage, and are likely to experience the most disruption because they cannot easily be broken down into separately managed regions. In the New Zealand case the largest city, Auckland, with 1.5 million people, is problematic to resolve. The three DHBs that comprise Auckland effectively behave as one big region and trying to halt most travel between them would not be a realistic option.

Future work

Perhaps the most immediate need is for better configuration data for the COVID‐19 pandemic (e.g., Chen, 2020). There is still a good deal of uncertainty around the major epidemic parameters, hence they can be changed in our model as better estimates come to light. Despite these uncertainties, it is possible to still run simulations using a variety of model settings and then perform a sensitivity analysis on the results. Since all the intermediate steps in each experiment can be logged and then analyzed (we use R scripts for this), it is possible to study which results hold across a range of parameter settings and plan accordingly.

More data is urgently needed on inter‐regional movement to improve the accuracy of the model. Although we know in broad measure how personal movement has been curtailed by lockdown, we do not know important details regarding the shipping of goods and food around the country, and specifically how much interaction and related disease transmission this produces. Such data does not need to wait for the next epidemic before being gathered. We will need to modify our settings as more accurate estimates come to light. None of these parameters could have been known in advance, which is perhaps the main reason why modeling this pandemic is so difficult, although models such as these could inform planning, even without parameters tuned for a specific virus.

Extending the model into other geographic dimensions

Our models simulate the progression of COVID‐19 within New Zealand using a population of 5 million residents. However, our models do not fully represent the age–sex structure or composition of the population nationally or regionally, which would better reflect exposure, risk, and transmission patterns, particularly among at‐risk populations (James, Plank, et al., 2020; New Zealand Ministry of Health, 2020; Steyn et al., 2020; van Dorn, Cooney, & Sabin, 2020).

Within an SIR modeling framework, population structure and composition would involve the addition of more compartments, while the transition rates between I or R could be modified using accurate regional flow data. For example, Towers and Feng (2012) extended their SIR model with age heterogeneity and a contact matrix to explore the impact social distancing had on the occurrence of influenza among the elderly. Their results demonstrated that a reduction in grandparent/grandchild interactions can reduce influenza among the elderly by up to 60%. Sattenspiel and Dietz (1995) improve SIR models with regard to their ability to model population mobility, with modifications that simulate no mobility to complete, permanent migration between all regions. The variants of the SIR models mentioned above, such as those by Balcan et al. (2010) and Gatto et al. (2020), highlight more contemporary implementations that introduce the complex nature of sub‐national population structure and mobility.

Adding spatio‐demographic complexity to SIR models essentially increases the number of containers to compartmentalize the population. However, regardless of the complexities introduced, the model is dependent on an aggregation of the population, and the number of conventional assumptions remain (Getz et al., 2019; Huppert & Katriel, 2013), such as the independence of the transmission probabilities. The addition of complex population structure, or better representation of mobility at different levels of lockdown, is better suited to the stochastic branching process models, which have the potential to be extended to individual‐based models (Jacob, 2010).