key: cord-0755726-psjyqlkc
authors: Grimee, M.; Bekker-Nielsen Dunbar, M.; Hofmann, F.; Held, L.
title: Modelling the effect of a border closure between Switzerland and Italy on the spatiotemporal spread of COVID-19 in Switzerland
date: 2021-05-20
journal: nan
DOI: 10.1101/2021.05.19.21257329
sha: 522dc664cd46cf08bc5ba221e59f47c1f63b30cd
doc_id: 755726
cord_uid: psjyqlkc

We present an approach to extend the Endemic-Epidemic (EE) modelling framework for the analysis of infectious disease data. In its spatiotemporal application, spatial dependencies have originally been captured by a power law applied to static neighbourhood matrices. We propose to adjust these weight matrices over time to reflect changes in spatial connectivity between geographical units. We illustrate this extension by modelling the spread of coronavirus disease 2019 (COVID-19) between Swiss and bordering Italian regions in the first wave of the COVID-19 pandemic. We adjust the spatial weights with data describing the daily changes in population mobility patterns, and indicators of border closures describing the state of travel restrictions since the beginning of the pandemic. We use these time-dependent weights to fit an EE model to the region-stratified time series of new COVID-19 cases. We then adjust the weight matrices to reflect two counterfactual scenarios of border closures and draw counterfactual predictions based on these, to retrospectively assess the usefulness of border closures. We observed that predictions based on a scenario where no closure of the Swiss-Italian border occurred, yielded double the number of cumulative cases over the study period. Conversely, a closure of the Swiss-Italian border two weeks earlier than implemented, would have only marginally reduced the number of cases and merely delayed the epidemic spread by a couple weeks. Despite limitations in the current study, we believe it provides useful insight into modelling the effect of epidemic countermeasures on the spatiotemporal spread of COVID-19.

The many types of uncertainty surrounding an ongoing emerging infectious disease outbreak with future 2 endemic potential, such as coronavirus disease 2019 (COVID- 19) , motivates the use of statistical modelling 3 to investigate current and future outbreaks of the disease. Not understanding the uncertainty of the true 4 scale of the disease, e.g. through unreported cases, leads to difficulties in assessing the impact of required 5 interventions [1, 2] . Multiple epidemic data sources provide valuable information on different aspects of such 6 an outbreak, but require appropriate statistical techniques to incorporate the associated uncertainties.

(c) Closure of the Swiss-Italian border on 3 rd March 2020, two weeks earlier than in the baseline scenario (counterfactual scenario B); 3. Determine the total and region-specific number of cases in Switzerland based on the baseline and 66 counterfactual scenarios from step 2; and 67 4. Examine the difference in total and region-specific incidence in Switzerland, to determine the usefulness 68 of a border closure. 69 To fit our model, we included data from 24 th February 2020, as cases have been recorded in Switzerland 70 and Italy after that date, until 4 th August 2020, as we only wanted to capture the dynamics of the first 71 wave of the COVID-19 epidemic in Switzerland. The cutoff date of 2 nd March 2020 for the (counterfactual) 72 predictions was chosen as the date after which the three scenarios might diverge. 73 Costantino et al. [12] consider similar scenarios when evaluating the effect of travel bans on COVID-19 74 spread in Australia. 75 76 We were interested in analyzing the COVID-19 spread in Switzerland and in the neighbouring regions in 77 Italy. We considered data from the regions on the second level of the Nomenclature of Territorial Units for 78 Statistics (NUTS-2). We included the seven Swiss regions, and the four bordering Italian regions, as listed in 79 Table 1 . A map of these eleven regions with their respective NUTS-2 code can be found in the supplementary For each of these regions, we considered daily case data, population data, daily or weekly testing data, 82 daily temperature data, data on public holidays, and data on border closures and population movement.

Information on data sources and manipulations is available in the supplementary materials. 

We considered the EE modelling framework. Given daily disease surveillance case counts Y rt indexed by 86 NUTS-2 region (r) and day (t), the region-stratified EE model was given by [13, 14] 87

where ψ r is a region-specific overdispersion parameter, e r are population fractions, ⌊·⌋ denotes normalised weights, ⌊u l ⌋ is the serial interval distribution, w rr ′ ,t−l is the transmission between regions r and r ′ at day t − l (see Section 2.4), p is the maximum lag, and ν rt and ϕ rt are log-linear predictors. The u l are shifted Poisson weights as recommended by Bracher and Held [14] for daily disease counts:

where κ is an unknown weighting parameter estimated from the data. We considered a maximum lag of p = 7

To determine optimal model specification, different model formulations were compared using model se-90 lection criteria. Here we used the Bayesian information criterion (BIC). In the final model, we included an 91 indicator for weekends and country-specific public holidays, as medical facility clinic and laboratory operating 92 hours may be reduced and population health-seeking activities may differ on weekends or public holidays 93 compared with weekdays. We also included a country indicator, as we expected to see differences in endemic-94 ity between regions located in Switzerland and regions located in Italy. Additionally, we included indicators 95 for each Swiss region bordering Italy (CH01, CH05, and CH07), as potential differences in underreporting 96 between Italy and Switzerland would make it harder to capture the whole spatiotemporal spread by the 97 power law alone. We also included country-specific (log-transformed) testing rates, as higher testing rates 98 were expected to increase the number of cases. In the epidemic component, we further included a gravity 99 model to reflect how more populated regions attract more imported cases [9, 15] . As a proxy for urbanisation, 100 we included the log-population of the largest city in each region [16, 17] . We included log-temperature and 101 yearly seasonality in our model, as we expected some form of seasonality in transmission related to seasonal 102 changes in human behaviour. Finally, we included a simple linear time trend. All continuous covariates were 103 centered.

The endemic predictor in our model was given by:

and the epidemic predictor was given by: 105 log(ϕ rt ) = α ϕ + β ϕ WEphl 1 {t is a weekend or public holiday} (t) + β ϕcountry 1 {r is in country j} (r) + β ϕ border 1 {r is CH01, CH05 or CH07} (r) + β ϕtesting log(T j(r)t ) + β ϕgravity log(P r ) + β ϕ city population log(D r ) + β ϕtime t + β ϕsin sin(2πt/365) + β ϕcos cos(2πt/365)

where α are overall intercepts, T j(r)t is the country-specific testing rate at time t, C rt is the region-specific 106 temperature at time t, D r is the population of the largest city in region r, P r is the population of r, and j(r) 107 describes the country in which region r is located and takes values "Switzerland" or "Italy". The predictor 108 β ϕgravity log(P r ) describes the gravity model [15] . . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 20, 2021.

1. To reflect the decreasing propensity of COVID-19 to spread between regions, as the distance between these regions increases, we implemented a power law model, as described in Meyer and Held [9] . Con-116 cretely, each entry of the matrix equals the inverse of the neighbourhood order, scaled by a decay 117 parameter d. That is, the matrix entries are given by:

We set the parameter d = 1.6 [9] . This parameter is scale-free; it is independent of the units of distance 119 [18] . It can thus be applied to neighbourhood order as a metric of distance.

All further manipulations were not symmetric anymore, as movement between regions is not expected 121 to be either. From here on, matrix rows represent "travel from" and columns represent "travel to". 122 2. The Facebook 2 Mobility data [19] allowed us to perform a daily adjustment of the matrix to reflect the 123 change in population mobility. The time-dependent entry is given by:

where f r,t is the value of the Facebook movement change variable for region r at time t (see supple-125 mentary materials and [20] ). The Facebook metric reports change in the quantitiy of movement from 126 inhabitants of a region r and not change in movement direction. Therefore we could not directly infer 127 the impact this change in movement has for the travel between r and each r ′ individually. We assumed 128 that it is uniform across all r ′ and therefore adjusted whole matrix rows with the Facebook metric. 

where b rr ′ ,t is the state of restrictions to travel from region r to region r ′ , at time t. The IOM metric 132 is qualitative and reports different border states as "No restrictions", "Partial restrictions", "Total 133 restrictions", and "No official restrictions reported". We quantified these different border states as 134 follows:

The time-dependent matrices are visualized in Figure 1 for the 1 st day of each month, after each step of 136 adjustment.

It should be noted that the ordering of regions in the matrices is arbitrary. 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 20, 2021. ; https://doi.org/10.1101/2021.05.19.21257329 doi: medRxiv preprint Figure 1 : Time-dependent weight matrices after each step of adjustment 6 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 20, 2021. To obtain predictions based on our EE model, we used the predictive moments approach [13, 21] . To Table 3 : estimated lag distribution for the EE model given by (1) 8 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 20, 2021 Given this model, data from the 8 days between 24 th February 2020 and 2 nd March 2020, and the 175 counterfactual weight matrices described in Section 2.5, we predicted total and region-specific cumulative is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 20, 2021. ; https://doi.org/10.1101/2021.05.19.21257329 doi: medRxiv preprint application, upon which all further adjustments rely. It is therefore crucial that it reflects the spatial con-211 nectivity correctly. In this context, it would have been important that the corresponding decay parameter 212 d be estimated from the data [13] . However, this parameter estimation is not (yet) compatible with the 213 R-implementation of the higher-order lags u l , which are also crucial in the study of COVID-19. In this study, 214 we chose to use a fixed parameter motivated by existing studies [9, 18] , but we recognise that this could 215 have weakened our time-dependent weights. Further work could examine how to facilitate harmonising the 216 two approaches. Similarly, the implementation of the IOM metric in this study might need appropriate 217 modification. As mentioned in the methods, the IOM metric is qualitative, and the way we chose to quantify 218 the indicators could be disputed. Indeed, we have chosen to weigh "total restrictions" as 0.1, meaning that 219 the strength of the spatial relationship between regions separated by a border would be weakened by a factor 220 of 10 upon closure of the latter. Whether this accurately reflects reality remains to be determined.

It is also crucial that we address model fit, which naturally influences the subsequent predictions. In initial 222 fittings of the EE model, we observed a particular underestimation of cases in the first wave for the regions 223 CH01, CH05 and CH07, the three regions bordering Italy. Ideally, region-specific intercepts would have been 224 implemented to address this, but this led to model divergence and thus unreliable parameter estimations. In 225 this paper, we chose to address this by fitting a model with region-specific overdispersion, which significantly 226 improved overall fit, and a specific indicator for the border regions, which improved fit in those regions.

However, in the present formulation of the model, we still see an underestimation of cases in most regions in 228 the first months. This leads us to conclude that there is more spatiotemporal spread happening than what 229 can be captured by the power law alone. We suspect that this is due to differing levels of underreporting, 230 particularly in the first wave, between Italy and Switzerland but also between individual regions. In light 231 of these limitations, we recommend that further studies take underreporting and underascertainment of case 232 counts into consideration. Simple multiplication factors can be used [23] as a first-order approach.

Moreover, a full picture of the border situation of Switzerland would be beneficial to gain precision, as 234 the spatiotemporal spread of COVID-19 in Switzerland arguably can not only be explained by importation 235 of cases from Italy and endemic factors. We intend to extend our work to additionally take the neighbouring 236 regions in France, Germany, and Austria into consideration. Finally, the EE framework can be adjusted 237 to incorporate pharmaceutical countermeasures such as vaccines [24] ; an obvious direction to pursue were a 238 similar analysis to be conducted for 2021.

All in all, we do believe that the present paper, extending the EE modelling approach with time-dependent 

We used time-dependent matrices of neighbourhood order between regions, adjusted through power laws, 245 mobility data, and border status indicators, allowing us to explore spatial interventions' impact on early cases implemented in reality, predictions showed only marginally reduced numbers of cases and an epidemic spread 252 merely delayed by a couple weeks, which is in line with the WHO's statement [10] . 253 We are unable to determine what an acceptable increase in cases would be as we are not policy makers, 254 and as a result have not considered cut-off criteria or hypothesis tests. However, the total number of cases 

Endemic-Epidemic framework used in COVID-19 modelling

Efficient real-time monitoring of an emerging influenza pandemic: How feasible?

A statistical framework for the analysis of multivariate infectious disease surveillance counts

Assessing the effect of containment measures on the spatio-temporal dynamic of COVID-19 in Italy

Modelling and predicting the spatio-temporal spread of coronavirus disease 2019 (COVID-19) in Italy

On the interplay of regional mobility, social connectedness, and the spread of COVID-19 in Germany. arXiv

COVID-19 in England: spatial patterns and regional outbreaks. medRxiv

Tracking and predicting the African COVID-19 pandemic. medRxiv

Power-law models for infectious disease spread

Updated WHO recommendations for international traffic in relation to COVID-19 outbreak

Spatiotemporal spread of infectious diseases and intra-or intertransmission: A gravity model approach. SSRN, 3731170

The effectiveness of full and partial travel bans against COVID-19 spread in Australia for travellers from China during and after the epidemic peak in China

Probabilistic forecasting in infectious disease epidemiology: the 13th Armitage lecture

Endemic-epidemic models with discrete-time serial interval distributions for infectious disease prediction

Measles metapopulation dynamics: A gravity model for epidemiological coupling and dynamics. the american naturalist

U.S. cancer death rates: A simple adjustment for urbanization

Modelling the risk of a disease in time and space

The scaling laws of human travel

Facebook data for good. movement range data

Protecting privacy in Facebook mobility data during the COVID-19 response

A marginal moment matching approach for fitting endemic-epidemic models to underreported disease surveillance counts

Erster bestätigter Fall in der Schweiz

Underreporting and reporting delays

Heterogeneity in vaccination coverage explains the size and occurrence of measles epidemics in German surveillance data