Modelling the effect of a border closure between Switzerland and Italy on the spatiotemporal spread of COVID-19 in Switzerland.

We present an approach to extend the endemic-epidemic (EE) modelling framework for the analysis of infectious disease data. In its spatiotemporal formulation, spatial dependencies have originally been captured by static neighbourhood matrices. These weight matrices are adjusted over time to reflect changes in spatial connectivity between geographical units. We illustrate this extension by modelling the spread of COVID-19 disease between Swiss and bordering Italian regions in the first wave of the COVID-19 pandemic. The spatial weights are adjusted 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. These time-dependent weights are used 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. Predictions based on a scenario where no closure of the Swiss-Italian border occurred increased the number of cumulative cases in Switzerland by a factor of 2.7 (10th to 90th percentile: 2.2 to 3.6) over the study period. Conversely, a closure of the Swiss-Italian border two weeks earlier than implemented would have resulted in only a 12% (8% to 18%) decrease in the number of cases and merely delayed the epidemic spread by a couple of weeks. Our study provides useful insight into modelling the effect of epidemic countermeasures on the spatiotemporal spread of COVID-19.

Introduction

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

One such statistical tool is the endemic–epidemic (EE) modelling framework, created for the analysis of infectious disease surveillance data (Held et al., 2005). It is a multivariate time series model that additively decomposes incidence into an endemic and an epidemic component. The epidemic component captures incidence driven by previous case counts, or force of infection, and the endemic component captures exogenous contributions to incidence (such as seasonal, socio-economic or demographic factors).

The best demonstration of the framework’s flexibility must be the wide range of infectious diseases it has been applied to. It can capture the incidence of infectious diseases of different types, with varying natural histories, vaccine preventability, and vector dynamics. In addition, the framework has already been widely used in analysis of COVID-19 outbreaks. Bekker-Nielsen Dunbar and Held (2020) give an overview of the model’s various applications. This modelling approach can handle the heterogeneities and interdependencies in high-dimensional count time series and is useful for analysing infectious disease surveillance data which is stratified by geographic regions.

The EE framework was originally designed to fit weekly surveillance data with an autoregression on counts from the previous week. The model has been recently extended to include higher-order lags (Bracher and Held, 2020a). This means that the autoregressive process can be driven by a larger time-span of past case counts. In light of the COVID-19 pandemic and the availability of daily case counts, in particular, this allows for the inclusion of infectiousness from the whole serial interval of COVID-19. Various epidemiological publications have used this feature in their analysis of COVID-19 incidence in Italy (Dickson et al., 2020, Giuliani et al., 2020), Germany (Fritz and Kauermann, 2020), England (Fronterre et al., 2020), and the African continent (Ssentongo et al., 2021) among others.

A key feature of the EE framework that will be discussed in the present paper is the possibility to include spatial dependencies, in addition to temporality. This allows for the modelling of transmission between regions over time. The spatial dependencies are included as weights in the autoregressive process of the epidemic component. The strength of spatial connectivity between regions, and the resulting transmissibility of the disease between them, can be formalised by the power law formulation introduced by Meyer and Held (2014). It reflects the decay of transmissibility between regions as the distance between them increases. Other groups have considered this approach in their endemic–epidemic modelling approaches (Fronterre et al., 2020, Ssentongo et al., 2021). It has for instance been used by Ssentongo et al. (2021) to model the effect of various socio-economic and meteorological factors, as well as containment strategies on the within- and between-country transmission in the African continent.

In this paper, we show that this approach can be extended with time-dependent weight matrices driving the epidemic component of the EE framework. We showcase how the static weights derived from the power law formulation can be adjusted over time to reflect the changing spatial relationship between geographical units as the mobility of the population shifts as a result of travel bans and border closures. This approach is of particular interest in the context of COVID-19, as in the first half of 2020 mobility-related social distancing interventions such as travel restrictions, border closures, and visit bans were imparted in many parts of the world, including within the Schengen area, in which the regions we consider here are contained. This was in spite of the World Health Organization advising against the use of such measures at the time (World Health Organization, 2020). Mobility restrictions are only considered to work if protective sequestration is feasible, i.e. autochthonous transmission has not yet occurred in the territory and so only imported cases are of concern. In our work we adjust the weight matrices according to both within-region mobility (measured via smartphones) and between-region mobility (using information from the International Organization for Migration). This provides a level of sophistication beyond simple mobility restriction/no mobility restriction indicators, as seen by e.g. Woo et al. (2020) who used a similar model to investigate the spatiotemporal spread of COVID-19.

The possibility to adjust spatial weights over time gives us the opportunity to simulate different counterfactual scenarios, representing different timelines of pandemic countermeasures targeting population mobility. This in turn allows us to retrospectively evaluate the effect of border closures as a pandemic countermeasure. We can thereby provide evidence to improve public health response and aid in decisions on optimal infectious disease control strategies, particularly when to initiate travel restrictions.

Methods

We examined the role of a border closure between Italy and Switzerland on the spatiotemporal spread of COVID-19 through the following steps:

Fit an EE model to the region-stratified time series of number of confirmed new COVID-19 cases from 24th February 2020 until 4th August 2020;

Retrospectively predict from the fitted model, given data until 2nd March 2020, under scenarios:

Closure of the Swiss-Italian border on 17th March 2020, i.e. as really implemented (base scenario (b))

No closure of the Swiss-Italian border (counterfactual scenario A)

Closure of the Swiss-Italian border on 3rd March 2020, two weeks earlier than in the base scenario (counterfactual scenario B);

Determine the total and region-specific number of cases in Switzerland based on the base and counterfactual scenarios from step 2; and

Examine the difference in total and region-specific incidence in Switzerland, to determine the usefulness of a border closure.

To fit our model, we included data from 24th February 2020, as cases have been recorded in Switzerland and Italy after that date, until 4th August 2020, as we only wanted to capture the dynamics of the first wave of the COVID-19 epidemic in Switzerland. The cutoff date of 2nd March 2020 for the (counterfactual) predictions was chosen as the date after which the three scenarios might diverge. Costantino et al. (2020) consider similar scenarios when evaluating the effect of travel bans on COVID-19 spread in Australia.

Data

We were interested in analysing the COVID-19 spread in Switzerland and in the neighbouring regions in Italy. We considered data from the regions on the second level of the Nomenclature of Territorial Units for Statistics (NUTS-2). We included the seven Swiss regions, and the four bordering Italian regions, as listed in Table 1. A map of these eleven regions with their respective NUTS-2 code can be found in the supplementary materials.

Table 1

Regions included in the analysis.

Country	Regions (NUTS-2 code)
Switzerland	Région lémanique (CH01), Espace Mittelland (CH02), Nordwestschweiz (CH03), Zürich (CH04), Ostschweiz (CH05), Zentralschweiz (CH06), Ticino (CH07)
Italy	Piemonte (ITC1), Valle d’Aosta/Vallée d’Aoste (ITC2), Lombardia (ITC4), Provincia Autonoma di Bolzano/ Bozen (ITH1)

For each of these regions, we considered daily case data, population data, daily or weekly testing data, 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.

Model

We considered the EE modelling framework. Given daily disease surveillance case counts indexed by NUTS-2 region () and day (), the region-stratified EE model was given by Held et al., 2017, Bracher and Held, 2020a

where is a region-specific overdispersion parameter, are population fractions, denotes normalised weights, is the serial interval distribution, is the transmission between regions and at day (see Section 2.3), is the maximum lag, and and are log-linear predictors. The are shifted Poisson weights as recommended by Bracher and Held (2020a) for daily disease counts: where is an unknown weighting parameter estimated from the data. Since other endemic–epidemic models of COVID-19 cases (Ssentongo et al., 2021) as well as estimates of the serial interval from the literature (Nishiura et al., 2020) indicate that its serial interval is likely to be in the range of days, rather than weeks, we considered a maximum lag of .

To determine optimal model specification, different model formulations including candidate covariates in the endemic, epidemic or both components, were compared using model selection criteria. Here we used the Bayesian information criterion (BIC). In the final model, we included an indicator for weekends and country-specific public holidays, as medical facility clinic and laboratory operating hours may be reduced and population health-seeking activities may differ on weekends or public holidays compared with weekdays. We also included a country indicator, as we expected to see differences in endemicity between regions located in Switzerland and regions located in Italy. Additionally, we included indicators for each Swiss region bordering Italy (CH01, CH05, and CH07), as potential differences in underreporting between Italy and Switzerland would make it harder to capture the whole spatiotemporal spread by the power law alone. We also included country-specific (log-transformed) testing rates, as higher testing rates were expected to increase the number of cases. In the epidemic component, we further included a gravity model to reflect how more populated regions attract more imported cases (Meyer and Held, 2014, Xia et al., 2004). As a proxy for urbanisation, we included the log-population of the largest city in each region (Kafadar and Tukey, 1993, Knorr-Held and Besag, 1998). We included log-temperature and yearly seasonality in our model, as we expected some form of seasonality in transmission related to seasonal changes in human behaviour. Finally, we included a simple linear time trend. All continuous covariates were centred.

The endemic predictor in our final model was given by: and the epidemic predictor was given by: where and are overall intercepts, is the country-specific testing rate at time , is the region-specific temperature at time , is the population of the largest city in region , is the population of , and describes the country in which region is located and takes values “Switzerland” or “Italy”. The predictor describes the gravity model (Xia et al., 2004).

Time-dependent matrix of neighbourhood order

To represent the spatial relationship between regions, driving the transmission between regions, as it changes over time a time-dependent matrix of neighbourhood order was included in our model. We started from a static matrix of neighbourhood order (see supplementary materials) and adjusted it to reflect the movement patterns and border closures between regions over time. We proceeded as follows:

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 (2014). Concretely, each entry of the matrix equals the inverse of the neighbourhood order, scaled by a decay parameter . That is, the matrix entries are given by: The decay parameter is scale-free; it is independent of the units of distance (Brockmann et al., 2006). 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 to be either. From here on, matrix rows represent “travel from” and columns represent “travel to”.

The Facebook2 Mobility data (Facebook, 2020) allowed us to perform a daily adjustment of the matrix to reflect the change in population mobility. The time-dependent entry is given by: where is the value of the Facebook movement change variable for region at time (see supplementary materials and Facebook Data for Good (2020)). The Facebook metric reports change in the quantity of movement from inhabitants of a region and not change in movement direction. Therefore we could not directly infer the impact this change in movement has for the travel between and each individually. We assumed that it is uniform across all and therefore adjusted whole matrix rows with the Facebook metric.

International Organization of Migration (IOM) data allowed us to further adjust the time-dependent matrix to reflect border closures. The time-dependent entry is given by: where is the state of restrictions to travel from region to region , at time . The IOM metric is qualitative and reports different border states as “No restrictions”, “Partial restrictions”, “Total restrictions”, and “No official restrictions reported”. We quantified these different border states as follows:

A sensitivity analysis around the quantification of the IOM metric can be found in the supplementary materials. The time-dependent matrices are visualised in Fig. 1 for the 1st day of each month, after each step of adjustment. The ordering of regions in the matrices is arbitrary. Note that the novelty of this approach to transmission weights is their time-dependent nature such that the spatial weights matrix does not just capture the relative contributions of cases from each region to each region in terms of their spatial relationship, but also how it changes over time. This also applies in particular to the case where , as we assume that a region’s contribution to its own cases also changes over time.

Fig. 1

Time-dependent weight matrices after each step of adjustment.

Counterfactual prediction of case counts

To examine the effect of border closures in place, we adjusted the neighbourhood matrices according to the three scenarios described under Section 2. For counterfactual scenario A, we implemented no border closures, which means that we only adjusted the matrix with the time-dependent Facebook metric . For counterfactual scenario B, we adjusted the matrix with and, additionally to the border closures implemented in the base scenario, we implemented travel restrictions from Italy to Switzerland () between 3rd March 2020 and 16th March 2020.

We obtained parameter estimates by fitting the EE model described in Section 2.2 to the region-stratified time series of confirmed COVID-19 cases, based on the time-varying matrix of neighbourhood order from the base scenario. Then, given data from 24th February 2020 to 2nd March 2020, we predicted from the model, based on the weight matrices adjusted for counterfactual scenarios A and B. We assessed both predictions to determine the usefulness of border closures, reflected in both the total number of cases and the cases by region.

To obtain predictions based on our EE model given by (1), (2), (3), and (4), we used the predictive moments approach (Held et al., 2017, Bracher and Held, 2020b). To compare the predicted scenarios, we considered the difference and ratio between each counterfactual scenario (A and B) and the base scenario, and thus calculated and , where is the cumulative (total or region-specific) number of cases under the counterfactual scenario at hand, and is the cumulative (total or region-specific) number of cases under the base scenario on 4th August 2020 (end of the study period).

Uncertainty

In our predictions, we take into account (1) internal uncertainty due to uncertainty in model coefficient estimates, and (2) external uncertainty caused by the decay parameter . Indeed, since the foundation of the spatial relationship between geographical units in this work is the power law application, upon which all further adjustments rely, it is crucial that it reflects the spatial connectivity correctly. In this context, it would have been important that the corresponding decay parameter be estimated from the data (Held et al., 2017). However, this parameter estimation is not (yet) compatible with the R-implementation of the higher-order lags, which are also crucial in the study of COVID-19. Therefore in this study, we sampled the decay parameter from a normal distribution motivated by Meyer and Held (2014). They estimated (95% : 1.61 to 2.01). This parameter was estimated in an analysis of the spatiotemporal spread of Influenza in Germany. We believe the Swiss setting and the German setting to be similar enough for the purpose of this analysis. We drew 100 samples from the corresponding normal distribution . Additionally, to incorporate the uncertainty around the model parameter estimates, we sampled the model parameters from a 31-dimensional multivariate normal distribution based on the estimated parameter vector and the associated variance–covariance matrix. Because the distribution of the endpoints over the 100 samples was highly skewed, the results were given as the median of the 100 predictions, with the corresponding 10th and 90th percentiles ( and ).

Implementation

This work used the EE modelling framework and associated hhh4 package ecosystem, in particular the R (version 4.0.2) packages surveillance (version 1.19.1), and hhh4addon (version 0.0.0.0.9013). Code and data for this work are available athttps://gitlab.switch.ch/suspend/COVID-19-border-CH.

Results

We fitted an EE model to the region-stratified time series of confirmed COVID-19 cases between 24th February 2020 and 4th August 2020 and the according covariates as described in Section 2. It includes a weekend and (national) public holiday indicator, a country indicator, an indicator for the Swiss border regions, testing rates, temperature, population, population of the largest city, a time trend and seasonality. The coefficients of the model and the corresponding standard errors are shown in Table 2, and the fitted values are depicted in Fig. 2. As is evident from Fig. 2, the model fit depends on time. In particular in the first months, we see a slight underestimation of cases in most regions. It is interesting to note the different proportions of endemic to epidemic incidence in the different regions. We see the highest share of endemic incidence in Italy and Ticino (CH07), followed by the other bordering regions Région lémanique (CH01) and Ostschweiz (CH05). Given that the first COVID-19 cases in Switzerland were recorded in Ticino, and imported from Italy (Bundesamt für Gesundheit BAG, 2020), this makes sense. The incidence in the French and German speaking regions of Switzerland can primarily be attributed to epidemic spread.

Table 2

Coefficient estimates and corresponding 95% confidence interval (2.5% and 97.5%) for the EE model given by (1), (2), (3), and (4).

		Estimate	2.5%	97.5%
Epidemic	β_ϕ_gravity	0.85	0.76	0.94
	β_ϕ_testing	0.36	0.28	0.44
	β_ϕ_time	0.02	0.01	0.03
	β_ϕ_WEphl	−0.48	−0.57	−0.39
	β_ϕ_country	1.98	1.77	2.18
	β_ϕ_border	0.18	0.09	0.27
	β_ϕ_city population	0.37	0.29	0.44
	β_ϕ_sin	−0.61	−0.80	−0.42
	β_ϕ_cos	1.82	1.13	2.51
	α_ϕ	−2.75	−3.01	−2.49
Endemic	β_ν_testing	0.99	0.75	1.23
	β_ν_time	0.23	0.17	0.29
	β_ν_WEphl	0.43	0.24	0.62
	β_ν_country	−1.91	−3.62	−0.20
	β_ν_border	3.14	1.63	4.64
	β_ν_temperature	0.04	0.02	0.05
	β_ν_city population	0.89	0.81	0.97
	β_ν_sin	−5.95	−7.16	−4.74
	β_ν_cos	19.62	15.10	24.15
	α_ν	0.25	0.80	1.97
Overdispersion	ψ_CH01	0.09	0.06	0.12
	ψ_CH02	0.07	0.04	0.10
	ψ_CH03	0.10	0.05	0.15
	ψ_CH04	0.21	0.13	0.30
	ψ_CH05	0.21	0.12	0.30
	ψ_CH06	0.11	0.06	0.17
	ψ_CH07	0.30	0.17	0.43
	ψ_ITC1	0.28	0.21	0.36
	ψ_ITC2	0.24	0.12	0.35
	ψ_ITC4	0.17	0.13	0.21
	ψ_ITH1	0.15	0.08	0.22

Fig. 2

Fitted values for the model given by (1), (2), (3), and (4).

The estimated lag distribution can be found in Table 3, and the corresponding plot in the supplementary materials. Our findings are consistent with the estimates obtained by Ssentongo et al. (2021).

Table 3

Estimated lag distribution for the EE model given by (1), (2), (3), and (4).

1	2	3	4	5	6	7
0.06	0.18	0.25	0.23	0.16	0.09	0.04

Given this model, data from the 8 days between 24th February 2020 and 2nd March 2020, and the counterfactual weight matrices described in Section 2.4, we predicted total and region-specific cumulative cases until 4th August 2020 under each scenario. The results are summarised in Table 4 and Fig. 3. To illustrate prediction accuracy, we plot a comparison of the daily incidence predicted under the base scenario, and observed incidence over the study period in the supplementary materials. Plots of the daily difference in new cases between the counterfactual scenarios and the base scenario are available in the supplementary materials for the complete study period.

Table 4

Counts, ratios, and differences in predicted cumulative cases for the counterfactual scenarios compared to the base scenario (b).

Territory	Base scenario (b)			Scenario A			Scenario B
	P10	Median	P90	P10	Median	P90	P10	Median	P90
Total	155818.4	174469.2	203711.6	223486.7	303558.8	461474.2	150822.9	169729.5	196765.4
Total (CH)	29560.8	36403.8	46430.7	69954.0	100472.2	166126.6	25541.0	31860.3	42519.8
Total (IT)	122976.1	138943.2	155873.2	153612.9	198950.8	287648.4	122821.3	138769.0	155614.4
CH01	12026.6	13505.6	16496.6	22447.0	30733.8	47611.2	11031.9	12396.5	15622.7
CH02	3570.7	4879.7	6346.5	9708.7	14838.6	25033.6	2947.2	4010.7	5619.0
CH03	2156.1	3177.3	4265.1	6122.5	9564.2	16704.1	1799.0	2561.6	3705.0
CH04	4813.8	7028.2	9380.3	14444.5	22120.8	37433.4	3836.2	5522.5	8185.3
CH05	2546.4	3061.3	3985.5	7069.6	10506.0	17428.7	2144.5	2582.1	3437.4
CH06	1575.0	2169.9	2801.2	4528.0	6801.3	11484.5	1307.5	1756.2	2504.4
CH07	2597.6	2898.7	3382.5	4664.6	6033.6	8600.4	2415.7	2714.2	3222.3
ITC1	27573.0	30955.6	34466.5	34895.2	45077.5	63349.2	27499.5	30912.4	34423.9
ITC2	997.6	1079.2	1181.8	1100.7	1209.2	1385.7	996.9	1078.6	1181.3
ITC4	90986.0	102897.3	117597.1	115031.4	149395.1	219346.6	90803.0	102715.9	117380.4
ITH1	3033.0	3248.5	3491.1	3408.3	3979.6	4887.8	3028.9	3245.5	3485.4
Territory	Ratio A/b			Ratio B/b			Difference A-b			Difference B-b
	P10	Median	P90	P10	Median	P90	P10	Median	P90	P10	Median	P90
Total	1.43	1.69	2.30	0.96	0.97	0.98	67089.2	124645.7	272510.7	−7585.0	−4597.6	−3172.1
Total (CH)	2.23	2.74	3.62	0.82	0.88	0.92	37961.1	64553.2	119687.8	−7257.6	−4414.9	−3067.4
Total (IT)	1.23	1.43	1.92	1.00	1.00	1.00	29372.5	58931.6	139171.6	−318.8	−190.5	−105.8
CH01	1.79	2.26	2.91	0.88	0.92	0.94	10094.1	17678.4	31455.7	−1839.0	−1139.8	−796.0
CH02	2.54	3.07	4.13	0.77	0.84	0.89	5879.5	9864.6	19687.0	−1185.5	−744.7	−500.6
CH03	2.52	3.10	4.12	0.76	0.84	0.90	3672.9	6299.2	12155.1	−830.5	−488.6	−320.1
CH04	2.68	3.29	4.33	0.75	0.83	0.89	8837.7	14984.6	29714.5	−1842.2	−1129.1	−744.5
CH05	2.71	3.37	4.36	0.80	0.86	0.90	4486.8	7428.3	13466.3	−656.0	−438.3	−316.5
CH06	2.65	3.23	4.27	0.77	0.84	0.89	2778.9	4651.8	8616.2	−537.2	−327.6	−221.2
CH07	1.71	2.14	2.66	0.91	0.94	0.96	1892.2	3240.2	5297.8	−264.3	−175.5	−126.7
ITC1	1.23	1.42	1.86	1.00	1.00	1.00	6923.0	13429.5	29314.9	−80.6	−48.1	−26.6
ITC2	1.07	1.12	1.26	1.00	1.00	1.00	67.8	125.2	268.3	−1.0	−0.5	−0.3
ITC4	1.23	1.45	1.96	1.00	1.00	1.00	22225.9	45726.1	107691.6	−234.5	−139.0	−77.4
ITH1	1.11	1.21	1.46	1.00	1.00	1.00	344.3	695.8	1546.5	−4.8	−2.7	−1.5

Fig. 3

Predictions under scenarios A, B, and base scenario b. The shaded area represents the uncertainty interval defined by and .

Under counterfactual scenario A, i.e. no border closure, we see a 2.7-fold (, ) increase in cases compared to baseline in Switzerland, with most excess cases happening in the first months. The peak number of excess cases compared to baseline would have been attained in early April with around additional daily cases (see supplementary materials for details). Overall, this scenario would have caused (, ) extra cases between 3rd March 2020 and 4th August 2020, with (, ) cases in Switzerland and (, ) in the bordering Italian regions.

Under counterfactual scenario B, i.e. earlier border closure, we see only a very small relative decrease in Switzerland (, , ). This would primarily have affected the early months, with a peak reduction of daily cases (c.a. ) around mid-March. From April onwards, no difference in daily cases would have been observed between scenario B and base scenario (b) (see supplementary materials for details). As we only implemented an internal border closure in Switzerland, i.e. only restricting entry in Switzerland from Italy, this decrease in cases naturally only affects the spread in Switzerland. Overall, this scenario would have prevented (, ) cases between 3rd March 2020 and 4th August 2020, of which the vast majority (, , ) in Switzerland.

Fig. 3 shows a timeline of predicted cases under each scenario for each of the regions. It is apparent that scenario A (blue) would have cause a marked increase in daily cases throughout the study period in all Swiss regions. The marginal effect of an earlier closure of borders (scenario B, yellow) can be clearly seen from the plot, as it remains barely distinguishable from the base scenario (pink) throughout the first wave and their uncertainty intervals overlap almost completely.

The results seem to support the base border closure scenario. Without the measures implemented, we would expect to have seen a tremendous rise in cases, whereas closing the border closures earlier would have merely delayed the spread of infection in Switzerland, and would overall have prevented only around cases.

Discussion

We highlighted the use of time-dependent weight matrices reflecting changes in spatial connectivity, to extend the traditional spatiotemporal EE framework. It is a particularly useful approach in the context of COVID-19, where pandemic countermeasures have targeted the spatial relations between geographical units with the implementation of border closures and travel restrictions, and where many data sources on patterns of population mobility have become available. This work also suitably captured most of the uncertainty inherent to our modelling approach.

Nonetheless, some limitations in the present paper need to be addressed to draw unbiased conclusions from our results. Given that the foundation of the spatial relationship between geographical units in this study is the power law application, it would have been important that the corresponding decay parameter be estimated from the data (Held et al., 2017) in order for it to reflect the spatial connectivity correctly. With the present limitations of the R-implementation in mind, in this study, we chose to sample the parameter from a distribution motivated by existing studies (Meyer and Held, 2014, Brockmann et al., 2006), but we recognise that this could have weakened our time-dependent weights. Further work could examine how to facilitate harmonising the two approaches. Similarly, the implementation of the IOM metric in this study might need further investigation. As mentioned in the methods, the IOM metric is qualitative, and the way we chose to quantify the indicators could be disputed. Indeed, we have chosen to weigh “total restrictions” as 0.1, meaning that the strength of the spatial relationship between regions separated by a border would be weakened by a factor of 10 upon closure of the latter. The sensitivity analysis exposed in the supplementary materials shows moderate robustness of our results to the choice of quantification of the multiplication factors. Regardless, the conclusions of this work remain the same in terms of support for the base border closure implementation. It is also worth noting that the model fit with the 0.1/0.5-parametrisation is better than with the other values tried in the sensitivity analysis.

It is also crucial that we address model fit, which naturally influences the subsequent predictions. In initial fittings of the EE model, we observed a particular underestimation of cases in the first wave for the regions CH01, CH05 and CH07, the three regions bordering Italy. Ideally, region-specific intercepts would have been implemented to address this, but this led to model divergence and thus unreliable parameter estimations. In this paper, we chose to address this by fitting a model with region-specific overdispersion, which significantly 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 the first months. This leads us to conclude that there is more spatiotemporal spread happening than what can be captured by the power law alone. We suspect that this is due to differing levels of underreporting, particularly in the first wave, between Italy and Switzerland but also between individual regions. In light of these limitations, we recommend that further studies take underreporting and underascertainment of case counts into consideration. Simple multiplication factors can be used (Noufaily, 2020) as a first-order approach.

Moreover, a full picture of the border situation of Switzerland would be beneficial to gain precision, as the spatiotemporal spread of COVID-19 in Switzerland arguably cannot only be explained by importation of cases from Italy and endemic factors. We intend to extend our work to additionally take the neighbouring regions in France, Germany, and Austria into consideration. Finally, the EE framework can be adjusted to incorporate pharmaceutical countermeasures such as vaccines (Herzog et al., 2011); an obvious direction to pursue were a 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 spatial weights, provides useful insight into the analysis of the spatiotemporal spread of COVID-19 and can be an important stepping stone for further research informing decision-making around COVID-19 countermeasures.

Conclusion

We used time-dependent matrices of neighbourhood order between regions, adjusted through power laws, mobility data, and border status indicators, allowing us to explore spatial interventions’ impact on early cases of COVID-19 in Switzerland and Italy. This approach has allowed us to draw counterfactual predictions based on weight matrices adjusted for counterfactual scenarios of border closures. We have observed that predictions based on counterfactual scenario A, where no closure of the Swiss-Italian border occurred over the study period, would have almost doubled the cumulative number of cases over the study period and all study regions. Additionally, based on counterfactual scenario B, where closure of the Swiss-Italian border occurred two weeks earlier than implemented in reality, predictions showed only marginally reduced numbers of cases and an epidemic spread merely delayed by a couple weeks, which is in line with WHO’s statement (World Health Organization, 2020). It is difficult to determine what an acceptable increase in cases would be, and as a result we have not considered cut-off criteria or hypothesis tests.

CRediT authorship contribution statement

Mathilde Grimée: Methodology, Software, Formal analysis, Visualisation, Writing. Maria Bekker-Nielsen Dunbar: Conceptualisation, Methodology, Software, Funding acquisition, Writing, Project administration. Felix Hofmann: Methodology, Software, Data curation, Visualisation. Leonhard Held: Conceptualisation, Methodology, Funding acquisition, Writing, Project administration, Supervision.