Where the virus travelled

In the spring of 2020 our research group in Ghent built a model of the COVID-19 epidemic in Belgium, to tell policymakers where the hospital numbers were heading. It followed people through ten stages of infection and nine age groups, but it treated the country as one well-stirred pot. For China, France or Italy that would have been a serious omission: there, how people moved between regions explained much of how the first wave spread. Belgium is different. It is small, it is one of the most densely populated countries in the world, and with five kilometres of road for every square kilometre of land it is also one of the best connected. Does it matter where Belgians travel, or does the virus reach every corner at once anyway?

Before building a spatial model we went looking for an answer in the data, in a paper with Gisele Miranda, Jenna Vergeynst, Joris Meys, Tijs Alleman and Jan Baetens. Its answer, in one sentence:

In the early stage of the epidemic, the movement of people and the spread of the virus were strongly correlated; once the virus had spread across the country, the correlation weakened.

Arrondissements well connected to the first outbreaks saw the most excess deaths five to six weeks later; by May, where the virus had started no longer predicted where it struck. Over all pairs of arrondissements, the better connected ones had hospital curves of more similar shape, running roughly in step, though that link is weak. The tabs below build this up the way the paper does, on Belgium's 43 arrondissements: did people move less when the virus was around, did the virus follow the connections to the first outbreaks, and do the local epidemic curves show an order and a shape that follow the connections too?

A word on the data. The paper rests on phone data from Proximus, hospital admissions per arrondissement from Sciensano and deaths per arrondissement from Statbel. The first two are confidential, so this page runs on public stand-ins, named where they are used: commuting from the 2011 census for the phone data, Sciensano's admissions per province for those per arrondissement, and, as a second view, its confirmed cases per municipality. The deaths are the real ones. Everything is computed in your browser; where I quote the paper, I say so.

Two waves

Belgium's first confirmed case was a person repatriated from Wuhan in early February 2020, who infected nobody. The virus arrived for real with holidaymakers coming home from northern Italy after the half-term break at the end of February. By early March it was spreading, on 18 March the country went into lockdown, and two days later the borders closed to all but essential travel. A second wave followed in the autumn.

New hospital admissions per day in Belgium, 2020, 7-day centred mean (Sciensano, public, summed over the provinces). Shaded: school holidays. Dotted: events, named when the pointer is on them. Above: the waves as the paper cuts them, and the two lockdowns.

The paper cuts 2020 into two waves with the same dates for every arrondissement. The first starts on 1 March, with the first confirmed cases, rises until the national peak of admissions on 28 March and falls until the summer holidays begin on 1 July. The second starts when the schools reopen on 1 September, rises until 30 October and falls until the end of the year. Each wave has a rising and a falling part. With the borders closed for much of this time Belgium was close to an island, which lets the analysis ignore travel abroad.

Fewer trips

The phone data come as a staytime: $P^{gh}(t)$ is the total time the residents of arrondissement $g$ spent in arrondissement $h$ on day $t$ . Add up the time spent outside the home arrondissement over a period of $\Delta t$ days and divide by the number of residents $N^g$ to get the outward mobility of $g$ :

m^g = \frac{1}{\Delta t}\,\frac{1}{N^g} \sum_{t} \sum_{h \neq g} P^{gh}(t).

Before the pandemic, from 20 to 28 February, the typical arrondissement spent about two and a half hours per resident per day elsewhere. While admissions rose in the first wave that fell to about two; while they fell, under lockdown, to about one and a half. The second wave started close to the old level and dropped again once the second lockdown came. (These are the medians of the paper's figure 5, read off by eye.) People moved less when the virus was around, which surprises nobody, but it is a useful check that the phone data measure what they should.

To compare arrondissements pair by pair, the paper condenses the staytimes into a connectivity index:

\mathrm{CI}^{gh} = \ln\left[\frac{1}{\Delta t} \sum_{t \,\in\, \text{wave}} \bigl(P^{gh}(t) + P^{hg}(t)\bigr)\right],

the average daily time the residents of $g$ spend in $h$ plus the reverse, on a log scale. It is symmetric, and on purpose not per head, so an arrondissement with many people is well connected. The analyses that follow all use rank correlations, so only the order of the index matters. Here the phone data are replaced by commuting. With $M^{gh}$ the number of people who live in $g$ and work in $h$ , from the 2011 census, this page uses

\mathrm{CI}^{gh} = \ln\bigl(M^{gh} + M^{hg}\bigr).

Commuting is the backbone of the phone data too, but it knows nothing about 2020: it stays the same before, during and after a lockdown.

Connected to the first outbreaks

Sciensano placed the first clusters of infections in the southwest of Limburg, in the arrondissements of Tongeren and Hasselt, and in the Borinage around Mons. If the virus travelled with people, then arrondissements strongly connected to those clusters should have been hit harder, some weeks later, than poorly connected ones. The paper measures "hit" by excess deaths: the deaths $D^g(w)$ in arrondissement $g$ in week $w$ of 2020, compared with the average over the same week of the ten years before,

\begin{gathered} \mathrm{ED}^g(w) = \frac{D^g_{2020}(w) - \bar D^g(w)}{\bar D^g(w)} \times 100\%, \\ \bar D^g(w) = \frac{1}{10} \sum_{y=2010}^{2019} D^g_y(w). \end{gathered}

Excess deaths need no test and no convention about which death counts as a COVID-19 death, which makes them about the least arguable measure there is. How well connection to an outbreak arrondissement $s$ predicts them is measured with Spearman's rank correlation over the 42 other arrondissements, week by week:

\rho_s(w) = \operatorname{corr}_{g \neq s}\bigl(\operatorname{rank} \mathrm{CI}^{sg},\ \operatorname{rank} \mathrm{ED}^g(w)\bigr).

It is 1 when the two rankings agree, 0 when they are unrelated and −1 when one is the reverse of the other.

connected to first outbreaks

outcome

week of

On the left, how many people commute between Tongeren (hatched) and each other arrondissement; on the right, the excess deaths in the week of 15 April. The two maps share a pattern: the northeast of the country, around and west of Tongeren, is dark on both. The scatter puts one against the other. For Tongeren in mid-April $\rho$ is 0.73 on these public data; the paper found 0.63 with the phone data, and 0.59 for Hasselt. The curve on the right shows $\rho$ week by week. It climbs above the chance band two to three weeks into March, stays there through April and falls away in May. From then on the virus is everywhere, and where it started no longer predicts where people die.

Press play and watch the excess deaths spread over the country and fade. Then try Mons: its connections predict much less, in the paper as here, although Sciensano counted it among the first clusters. Its admissions tell the opposite story: in the paper they run ahead of nearly every other arrondissement, by up to 6.3 ± 1.3 days on Namur, where Tongeren and Hasselt do not stand out. Any other arrondissement can be picked too, by clicking it on either map. With 53 weeks on the curve, the odd week crosses the band by chance alone; what counts is a run of weeks in a row, in spring.

In France the corresponding effect of the first large outbreak, in the Haut-Rhin, faded within two weeks of the lockdown. In Belgium it lasted about four weeks longer. The paper suggests two reasons: no single large spreading event started the Belgian epidemic, so the excess deaths built up gradually, and in a country this connected the best- and worst-connected arrondissements differ less than they do in France.

Hospital admissions per arrondissement are not public, but confirmed cases per municipality are: switch the outcome to see them. Sciensano publishes a day with one to four cases only as "<5", but its exact totals per province and day are public too, so the page spreads what a province's "<5" days add up to evenly over them. The chart of $\rho$ always shows both outcomes, the one on the map as the solid line. The paper left cases out because they depend on how much was tested, and the two curves show both why and what cases would have added. Until mid-March, cases per head were highest in the arrondissements best connected to Tongeren and Hasselt, with $\rho$ about 0.6, some three weeks before the deaths show the same pattern. On 11 March testing was restricted to hospital patients and health workers with symptoms. From then on a case said less about where the virus circulated than about who was admitted or worked in care, and $\rho$ collapses; after that it wanders in and out of the chance band. For Tongeren it even turns negative in the autumn: the second wave hit Wallonia and Brussels hardest, and Limburg least.

Reading a noisy curve

The second half of the paper compares the epidemic curves of the arrondissements themselves: which ran ahead, and which had the same shape. Daily admissions in an arrondissement are noisy, from weekday effects and from plain chance, and a small arrondissement may see nobody for days and then three patients at once. Comparing such curves directly would mostly compare noise. So each curve is first smoothed with a generalised additive mixed model. The logarithm of the expected number of admissions $y(t)$ on day $t$ is a sum of $k$ smooth bumps, B-splines $B_j(t)$ , about one per week, and the counts may scatter more widely around it than pure chance would make them (overdispersion, $\phi > 1$ ):

\begin{gathered} \log \mathbb{E}\,y(t) = \sum_{j=1}^{k} \beta_j B_j(t), \\ \operatorname{Var} y(t) = \phi\,\mathbb{E}\,y(t). \end{gathered}

The coefficients are fitted with a penalty against wiggling: with $\mathcal{D}(\beta)$ the deviance, the misfit to the counts,

\hat\beta = \arg\min_{\beta} \Bigl[\, \mathcal{D}(\beta) + \lambda \sum_{j} \bigl(\beta_{j+1} - 2\beta_j + \beta_{j-1}\bigr)^2 \Bigr].

The best fit $\hat\beta$ is one curve. But it comes with an uncertainty, a covariance matrix $\hat V$ that grows with $\phi$ , and drawing

\beta^{(r)} \sim \mathcal{N}\bigl(\hat\beta,\ \hat V\bigr), \qquad r = 1, \dots, 100,

gives 100 curves: a fan as wide as the data are unsure. Every result from here on is computed on all 100 and reported as a mean and a spread.

show

On top, Brussels-Capital: 1.2 million people and a clean curve, where the 100 fits are almost one line. Below, Diksmuide: 52 000 people. Its admissions are confidential, so the figure keeps each patient admitted in West-Vlaanderen with the probability that they lived in Diksmuide, which gives counts with the noise of a population that size. The dots are now mostly zeros and single patients, the fan opens up, and the date of the peak is uncertain by days. Near the ends of the window, where a stray patient or two carry the whole fit, the fan can flare up: that is the model admitting how little it knows there. "New sample" draws another Diksmuide; any province or arrondissement can go in either panel. This is why the paper stops at arrondissements: for the 581 municipalities most curves would look like this one or worse.

The fit here is the paper's without one refinement. Its model also lets the random deviations of neighbouring days correlate, a first-order autoregressive term, which makes the fans somewhat wider than these. The spline has about one basis function per week, as in the paper, and the weight $\lambda$ of the penalty is chosen by generalised cross-validation.

Ahead or behind

Two tools compare the smoothed curves pair by pair. Both first rescale the curves to zero mean and unit spread, $s^g(t)$ and $s^h(t)$ . Time-lagged cross-correlation then slides one along the other and asks at which shift $k$ , in days, they overlap best:

\begin{gathered} C^{gh}(k) = \sum_{t \,\in\, \text{wave}} s^g(t + k)\, s^h(t), \\ \mathrm{TLCC}^{gh} = \arg\max_{k}\, C^{gh}(k). \end{gathered}

A negative lag means $g$ runs ahead of $h$ . Dynamic time warping then lines the two curves up at that shift and measures the difference in shape that is left. It may stretch and squeeze one curve in time to match the other as closely as possible, and reports the average gap that remains:

\mathrm{DTW}^{gh} = \frac{1}{2n} \min_{\pi} \sum_{(i,j) \,\in\, \pi} w_{ij}\, \bigl|s^g(i) - s^h(j)\bigr|.

Here $n$ is the number of days the shifted curves overlap, and $\pi$ runs over the ways of matching the days of one curve to those of the other in order, from the first day to the last; a step forward in both curves at once counts twice ( $w_{ij} = 2$ ), a step in one of them once. The lag is a number of days; the DTW distance has no unit, and smaller means more alike.

colour by

shift

With the paper's data these are heatmaps of 43 by 43 arrondissements; with the public data they are 11 by 11, one row and column per province, with Brussels as the eleventh. Click a cell to see its pair. Drag along the cross-correlation strip or use the slider to slide the second curve over the first, "best shift" jumps to the top of $C(k)$ , and "warping" draws which day of one curve the time warping matches to which day of the other.

The paper's findings on timing survive the coarser data. In the first wave the provinces run within about a day of each other, with Namur trailing, as it does in the paper. In the second the lags roughly double, with Brussels and Liège clearly ahead and West-Vlaanderen and Limburg behind, much as the paper found for their arrondissements; it confirmed that the lags grew with a one-sided Wilcoxon signed-rank test (p < 0.001). The shapes, too, differ more in the second wave: switch the colouring to DTW and flip between the waves.

The line under the figure brings the two halves together, as the paper's figure 11 does: over all pairs, does a stronger connection go with a shorter lag and a more similar shape? Over its 903 pairs of arrondissements the paper found both correlations negative, as expected, but weak: −0.15 for the lags and −0.31 for the shapes in the first wave. The provinces tell the same story, if less sharply for the lags: between provinces those are small, often under a day, and leave little to rank. Connection is one influence on how the virus moved, but far from the only one: why people travel, how many stops they make and the short trips within an arrondissement, none of which these data see, all matter too.

What it adds up to

The three steps build up to the paper's main conclusion:

In the early stage of the epidemic, the movement of people and the spread of the virus were strongly correlated; once the virus had spread across the country, the correlation weakened.

People moved less while the virus was around, so the phone data measure what they should. Connection to the first outbreaks, in Tongeren and Hasselt, predicted where the excess deaths would be highest five to six weeks later, and stopped predicting anything by May. Over all pairs of arrondissements, a stronger connection went with hospital curves of more similar shape that ran more nearly in step, but only weakly: other factors, beyond what these data see, cloud a clean correlation. None of this proves that travel carried the virus; it shows that the data look the way they would if it did, early on.

It also says what a model of the epidemic in Belgium should contain: a notion of mobility, above all for the early stage, before the virus has spread evenly through the country. So we built one: a spatially explicit version of our model, run on the same phone data, for the provinces and the arrondissements. Its verdict is modest in the way this analysis foreshadows: once the virus had spread, travel between provinces barely changed the course of the epidemic, and cutting social contact did far more than restricting movement.

Two caveats about this page. Commuting in 2011 is not movement in 2020: it misses shopping, visits and the school run, and it cannot drop during a lockdown the way the phone data did. And the provinces average away most of the geography the paper was about. Where the public data agree with the paper, as they mostly do, that is reassuring; where they differ, these are the data I would doubt first.

Mobility and the spatial spread of SARS-CoV-2 in Belgium
M. Rollier, G. H. B. Miranda, J. Vergeynst, J. Meys, T. W. Alleman, the Belgian Collaborative Group on COVID-19 Hospital Surveillance, J. M. Baetens
Mathematical Biosciences 360, 108957, 2023 · preprint
A stochastic mobility-driven spatially explicit SEIQRD COVID-19 model with VOCs, seasonality, and vaccines
T. W. Alleman, M. Rollier, J. Vergeynst, J. M. Baetens, Applied Mathematical Modelling 123, 2023. The spatial model this analysis led to.
Assessing the effects of non-pharmaceutical interventions on SARS-CoV-2 transmission in Belgium
T. W. Alleman et al., Epidemics 37, 2021. The national model that came first.
How human mobility explains the initial spread of COVID-19
S. M. Iacus et al., Joint Research Centre, 2020. The French and Italian analysis the connectivity index comes from.
UGentBiomath/COVID19-Model
The code of the group's models, with the public data used on this page: the census commuting matrix, the weekly excess deaths of 2020 and the arrondissement borders.

Public data: hospital admissions per province and confirmed cases per municipality from Sciensano (epistat.sciensano.be), weekly excess deaths per arrondissement and the 2011 commuting matrix as prepared in the repository above, borders simplified with mapshaper. Keys once a figure has focus: space plays the weeks, the arrows step through them or shift the curve, m switches between deaths and cases, r resets, n draws a new sample, 1 and 2 pick a wave, b jumps to the best shift, w shows the warping, f goes fullscreen. The address bar follows the settings, so a view can be linked.

Where the virus travelled

Two waves

Fewer trips

Connected to the first outbreaks

Reading a noisy curve

Ahead or behind

What it adds up to

Read more