# Estimates of the reproduction ratio from epidemic surveillance may be biased in spatially structured populations

Piero Birello, Michele Re Fiorentin, Boxuan Wang, Vittoria Colizza, Eugenio Valdano

<figure><img src="image_1.png" /></figure>

## To cite this version:

Piero Birello, Michele Re Fiorentin, Boxuan Wang, Vittoria Colizza, Eugenio Valdano. Estimates of the reproduction ratio from epidemic surveillance may be biased in spatially structured populations. Nature Physics, 2024, 20 (7), pp.1204-1210. 10.1038/s41567-024-02471-7. hal-04952686

HAL Id: hal-04952686

https://hal.sorbonne-universite.fr/hal-04952686v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Estimates of the reproduction ratio from epidemic
surveillance may be biased
in spatially structured populations

Piero Birello¹, Michele Re Fiorentin², Boxuan Wang¹,
Vittoria Colizza¹, and Eugenio Valdano¹,*.

¹Sorbonne Université, INSERM, Institut Pierre Louis d'Epidémiologie et de Santé Publique, F75012, Paris, France.

²Department of Applied Science And Technology (DISAT), Politecnico di Torino, 10129, Torino, Italy.

* Corresponding author eugenio.valdano@inserm.fr

**Abstract**

Accurate estimates of the reproduction ratio are crucial to project infectious disease epidemic evolution and guide public health response. Here, we prove that estimates of the reproduction ratio based on inference from surveillance data can be inaccurate if the population comprises spatially distinct communities, as the space-mobility interplay may hide the true epidemic evolution from surveillance data. Consequently, surveillance may underestimate the reproduction ratio over long periods, even mistaking growing epidemics as subsiding. To address this, we use the spectral properties of the matrix describing the spatial epidemic spread to reweigh surveillance data. We propose a correction that removes the bias across all epidemic phases. We validate this correction against simulated epidemics and use COVID-19 as a case study. However, our results apply to any epidemic where mobility is a driver of circulation. Our findings will help improve epidemic monitoring and surveillance and inform strategies for public health response.
Main text

The reproduction ratio *R* is arguably the most used indicator to monitor the trend in the evolution of an infectious disease epidemic. *R* is the average number of secondary infections that each case generates: When it is larger than one, the epidemic wave is growing; when instead it is lower than one, it is subsiding [31, 46]. The reproduction ratio also measures the effectiveness of public health interventions, whose overarching goal is to bring an unconstrained epidemic (*R* > 1) below the epidemic threshold of *R* = 1: Accurately estimating the reproduction ratio is thus necessary to ascertain the current epidemic evolution, predict short-term trends, perform scenario analysis and plan public health action [65, 51, 60, 17, 49]. The standard way to measure *R* is to infer it from data coming from epidemiological surveillance [64, 15, 48, 14, 59]. These data may be timelines of detected cases or their proxies, like hospitalizations or deaths, and this approach applies to diseases spanning radically different epidemiology, transmission routes and burden, like influenza [6, 61], measles [24], COVID-19 [34], Ebola [38], cholera [42], dengue [11], malaria [53]. The resulting surveillance-based estimates of *R* are routinely used to design interventions [2]: Notwithstanding, we argue in this study that surveillance data may lead to biased estimates of the reproduction ratio in spatially structured populations, where geographically distinct communities (e.g., cities) are connected though human mobility. Heterogeneities in the contact network shape the way epidemics spread [50]. Here, we will show that those arising from the complex interplay between spatial transmissibility patterns and the mixing network driven by human mobility hide the true dynamic structure of the epidemic process from population-level surveillance data. This mirrors the nature of most mathematical models of epidemic spread: they integrate space and spatial data at high resolution [27, 4, 50, 57, 25, 9], but they find it harder to do the reverse, which is extracting high-resolution information from limited and coarse-grained surveillance data in the absence of knowledge of the underlying spatial dynamics [12, 54, 8]. Crucially, this means that inference on surveillance data may either overestimate or underestimate it over long periods. This is of great public health relevance: measuring for instance a reproduction ratio below one when the true value is above would falsely signal that the epidemic is under control. Here, we study this bias, identify its origin and compute its magnitude. Then, we propose a correction to case incidence data that removes this bias and ensures that surveillance-based estimates of the reproduction ratio consistently give the true reproduction ratio of the epidemic. Our theoretical findings apply to any epidemic featuring relatively short generation time and for which mobility is a contributing factor in shaping its circulation within and across communities. This covers some of the
global health threats that are being worst affected by climate change and demographic trends: viruses responsible for respiratory infections – including SARS-CoV-2 and influenza – [35], vector-borne pathogens – including the arboviruses dengue, chikungunya, Zika [41, 52], and emergence events of new viruses or new viral strains [7]. To test and illustrate our findings, we use the French COVID-19 epidemic (see Fig. 1) before the advent of vaccination as a case study.

## Theoretical formalism

The Galton-Watson branching process is a customary framework to model epidemic spread [66, 37, 26]. Let $I(0)$ be the initial number of infections, $I(1)$ the expected number of infections that the initial cases generate, and, generally, let $I(t)$ be the expected number of infections in the $t$-th generation. By definition of the reproduction ratio, we have that $I(t) = RI(t-1)$, which implies that $I(t) = R^t I(0)$. This equation means that the number of infections grows exponentially if $R > 1$. In any real outbreak other factors, like acquired immunity, seasonal effects or public-health interventions, will at some point curb this exponential growth by changing the value of $R$. Notwithstanding, we may assume $R$ to be fairly constant either in the early phase of an outbreak, when those effects have not yet kicked in, or when the timescale at which immunity and mixing change is much longer than epidemic evolution [32, 16].

In the case of a population composed of $N$ spatial communities, we may define the vector $I(t) \in \mathbb{R}^N$, whose component $I(t)_i$ is the number of infections in generation $t$ and community $i$. Likewise, the reproduction operator $\mathbf{R} \in \mathbb{R}^{N,N}$ encodes, in its component $R_{ij}$, the average number of infections generated among the residents of community $i$, by a case belonging to community $j$ [58]. Its definition thus matches that of a spatially-structured next-generation matrix [19, 47]. However, unlike the standard next-generation matrix, it can be built at any phase of the epidemic and not only close to the disease-free state, and with no need of assuming a specific transmission model. The definition of $\mathbf{R}$, and the results that we are going to derive from it, thus applies to any epidemic and disease. The specific parametrization of $\mathbf{R}$ will instead depend on the specific transmission dynamics and natural history of the disease: for directly-transmitted diseases $\mathbf{R}$ typically depends on mixing patterns among communities [39]; for vector-borne diseases the local abundance of the host vectors, modulating the effective transmissibility, needs to be factored in, too [41, 30]. The expected epidemic evolution then follows the equation

$$I(t) = \mathbf{R}^t I(0). \qquad (1)$$
$I(t)$ encodes both the total number of infections in the population in generation $t$ and its spatial distribution. We define the former as the number $I_{tot}(t) = \sum_{i} I(t)_i$ and the latter as the vector $\mathbf{x}(t) \in \mathbb{R}^N$ whose components are $x(t)_i = I(t)_i / I_{tot}(t)$.

The reproduction ratio $R$ of this process is the eigenvalue of $\mathbf{R}$ with the largest absolute value (spectral radius) [19, 47], which is real, positive and nondegenerate (see Supplementary Methods Section 1.1). Its corresponding eigenvector $\mathbf{v}$ is also strictly positive ($v_i > 0$) we normalize it so that $\sum_i v_i = 1$.

Measuring the reproduction ratio of the system thus requires knowledge of the spectral structure of $\mathbf{R}$, i.e., of the spatial structure of the epidemic. We call this reference reproduction ratio as opposed to what surveillance instead estimates from the evolution of the incidence of infections or their proxies. This may happen globally, at the level of the entire population, or locally in each community. In our framework, the population-level estimated reproduction ratio is $R^{estim}(t) = I_{tot}(t+1)/I_{tot}(t)$, i.e., the general growth rate. The local community-level estimated reproduction ratio is instead $r_i^{estim}(t) = I(t+1)_i/I(t)_i$.

A simple observation then underpins our study: in general $R^{estim}(t)$ and $r_i^{estim}(t)$ may be different from $R$, the spectral radius of $\mathbf{R}$, and, if that is the case, surveillance will not estimate the reference reproduction ratio.

To explore this, we will first determine the conditions leading to an unbiased measure of the reproduction ratio: $R^{estim}(t) = R$.

## When reference and estimated reproduction ratios match

By virtue of the Perron-Frobenius theorem, $\mathbf{R}^t \to \mathbf{R}^t \mathbf{v}\mathbf{v}^*$ asymptotically at large $t$, where $\mathbf{v}^*$ is the dual of $\mathbf{v}$ (and left Perron eigenvector of $\mathbf{R}$) and normalized so that $\mathbf{v}^*\mathbf{v} = 1$. Equation (1) becomes $\mathbf{I}(t) \to [\mathbf{v}^*\mathbf{I}(0)] \mathbf{R}^t \mathbf{v}$, which implies that $\mathbf{x}(t) \to \mathbf{v}$, $R^{estim}(t) \to R$ and $r_i^{estim}(t) \to R$. The epidemic dynamics thus brings the spatial distribution of infections toward $\mathbf{v}$, which we will refer to as the *equilibrium spatial distribution of infections*. Thus, for any epidemic dynamics, if infections are spatially distributed as the equilibrium distribution ($\mathbf{x} = \mathbf{v}$), then the error is zero and the reference reproduction ratio is measured both globally ($R^{estim} = R$) and locally ($r_i^{estim} = R$). The last equality implies the existence of a detailed balance between locally generated, imported and exported infections in every community (see Supplementary Methods Section 1.2). Fig. 1 shows evidence of the convergence to $\mathbf{v}$ during the COVID-19 epidemic in France in late 2020 and early 2021. We reconstructed $\mathbf{R}$ from mobility data from Meta [28], a multinational technology company, and $\mathbf{x}$ from hospitalization data (see Fig. 1, Reconstruction of the
reproduction operator from data and Supplementary Methods Section 1.3). In a period when R was fairly constant the angle between **x** and **v** consistently decreased, a sign that **x** was evolving towards the equilibrium distribution, as the theory predicts.

Finally, there is a class of operators $\mathbb{R}$ for which the reference and population-level estimated reproduction ratios match at all times and even out of equilibrium. Indeed, if $\rho_i \equiv \sum_j R_{ji}$, the expected total number of secondary infections, is the same across communities, then $R^{estim} = \rho_i = R \forall i$ identically (see Homogeneous transmission potential). This implies that only the combination of spatial epidemic coupling and spatial heterogeneity in the transmission potential $\rho_i$ may cause $R^{estim}$ to be different from $R$. Notably, this is true only globally and $r_i^{estim}$ may be different from $R$ even if $\rho_i = R$.

<figure><img src="image_3.png" /><figcaption>**Fig.1** | Convergence to the equilibrium spatial distribution of infections. **a**, Estimate of the daily number of new COVID-19-related infections per 100,000 residents from week 9 of 2020 to week 3 of 2022 (Feb 23, 2020 to Jan 23, 2022). See Supplementary Methods Section 1.3 for how we estimated infections from hospitalizations and the comparison with serological surveys. Each curve is that of one of the 94 departments (administrative level 2) of mainland France excluding the region of Corsica, with a 14 days window rolling average. The period from week 2 to week 9 of 2021 (Jan 11 to Mar 7), selected as a period in which **v** was fairly constant, is highlighted in grey both in **a** and **b**. The small dept. Territoire de Belfort (0.2% of the French population) is excluded due to a miscounting of hospitalizations coming from neighboring departments by the health authorities - see Supplementary Fig. S3 and S4. **b**, The orange curve (left y axis) reports the weekly cosine distance between **v** (computed from mobility data from that week - see Reconstruction of the reproduction operator from data) and its average ($\bar{v}$) on the time interval from week 51 of 2020 to week 12 of 2021 (see Supplementary Fig. S13). The cosine distance between **v** and $\bar{v}$ is defined as $1 - v \cdot \bar{v} / (\|v\|\|\bar{v}\|)$. The green curve (right y axis) reports the cosine distance between $\bar{v}$ and the weekly spatial distribution of COVID-19 estimated infections **x** in French departments. The two inset maps show the spatial distribution **x** in week 2 and week 9 of 2021, respectively (edges of the time interval in gray).</figcaption></figure>
When reference and estimated reproduction ratios do not match

We now focus on the out-of-equilibrium dynamics (x(t) ≠ v) and measure the bias on R as the relative difference between the estimated and the reference reproduction ratios

$$ \Delta(t) = \frac{R^{estim}(t) - R}{R}. \qquad (2) $$

We call $\Lambda_{\alpha}$ ($\alpha = 1, \dots, N-1$) the (possibly degenerate) eigenvalues of R other than R and, by Perron-Frobenius theorem, $|\Lambda_{\alpha}| < R$. With calculations reported in Calculation of $\Delta(t)$: proof of equation (3), we find that

$$ \Delta(t) = C(t) \sum_{\alpha} z_{\alpha} \left(1 - \frac{\Lambda_{\alpha}}{R}\right) \left(\frac{\Lambda_{\alpha}}{R}\right)^{t}, \qquad (3) $$

where $C(t)$ is positive and asymptotically constant, and $z_{\alpha}$ is a (possibly complex) number proportional to the projection of the initial condition x(0) on the $\alpha$-th mode. The modes in equation (3) for which $\Lambda_{\alpha} \approx R$, or that are almost orthogonal to the initial configuration x(0), are suppressed from the start and do not bias the estimate of the reproduction ratio. The other modes, instead, possibly do, with an effect that becomes smaller with a characteristic decay time $\tau_{\alpha} = 1/\log(R/|\Lambda_{\alpha}|)$. In addition, those modes for which $\Lambda_{\alpha}$ is not real and positive oscillate with period $T_{\alpha} = 2\pi/|\theta_{\alpha}|$, where $\theta_{\alpha} = \arg \Lambda_{\alpha}$ (with $\theta_{\alpha} \in (-\pi, \pi]$) – see Calculation of $\Delta(t)$: $\tau_{\alpha}, T_{\alpha}$. These oscillations will be visible if faster than their characteristic decay time: $T_{\alpha} \le \tau_{\alpha}$. This gives the inequality

$$ \frac{|\Lambda_{\alpha}|}{R} \ge e^{-\frac{|\theta_{\alpha}|}{2\pi}} \ge e^{-\frac{1}{2}} \approx 0.61. \qquad (4) $$

The strict equality occurs when $\Lambda_{\alpha}$ is real and negative ($\theta_{\alpha} = \pi$). Furthermore, using the same approach, we prove in Supplementary Methods Section 1.2 that oscillations may be visible in the local estimated reproduction ratios even when the global $\Delta(t)$ never changes sign (see Fig. 2e).

To test the predictions of our theory in a realistic scenario, we considered again the COVID-19 epidemic in France and built a stochastic metapopulation model using the same data as in Fig. 1b. The details of the model are reported in Epidemic simulations and Supplementary Methods Section 1.4. We measured the reference and the estimated reproduction ratios, reported in Fig. 2, which shows that surveillance-based estimates may remain consistently biased for a long period and, depending on where the epidemic wave started, they may either overestimate or underestimate the reference reproduction ratio (see Supplementary Fig. S9 for more initial conditions tested). The case depicted
in Fig. 2d is particularly concerning: during the first month of the simulated epidemic, surveillance records a lower-than-one reproduction ratio which indicates decreasing incidence (visible in Fig. 2b), mistakenly pointing to a subsiding outbreak. In reality, the reproduction ratio is well above one, and only after two months of simulated epidemic does the surveillance based estimate reach this value. Alongside the estimate of $R^{estim}$ given within the framework of the Galton-Watson process, in Fig. 2c,d we also provide an estimate of the reproduction ratio by feeding incident cases to the library *EpiEstim* [14], a popular tool to compute the reproduction ratio from surveillance data. The fact that the two measures overlap confirms that the Galton-Watson process correctly reproduces the phenomenology under study even in realistic scenarios. Notwithstanding, more detailed frameworks [19, 69, 63] could be used to study the impact of heterogeneous generation intervals.

Finally, Fig. 2e and Fig. 2f show that locally measured reproduction ratios converge to the reference value at different times and with different speeds and that, at the same moment in time, some communities overestimate $R$ and some underestimate it (analytical proof in Local over- and underestimation of $R$.

Fig. 2 shows no oscillations in the sign of $\Delta(t)$, compatible with the fact that the reconstructed operator $\mathbf{R}$ has only real and positive eigenvalues. We extended our analysis to 32 European countries and none showed oscillations in $\Delta(t)$ at the considered spatial resolution (admin-2). We discuss the details of this analysis and the possible realistic conditions leading to oscillations in Supplementary Methods Section 1.5 and Supplementary Fig. S1.

In the absence of oscillations, the estimated reproduction ratio consistently either overestimates or underestimates the reference reproduction ratio, as $\Delta(t)$ decays to zero without ever changing sign. In this case, we can determine the sign of the bias from the initial condition: $\Delta(0) = \sum_j \rho_j x(0)_j - R$. By the Perron-Frobenius theorem, $j_{min}, j_{max}$ exist so that $\rho_{j_{min}} \le R$ and $\rho_{j_{max}} \ge R$. Thus, the initial location of infections will completely determine the sign of the error that surveillance will make. If the epidemic starts $j_{max}$ – or in general in communities with high transmission potential – surveillance will consistently overestimate the reference reproduction ratio until the bias decays to zero. Conversely, if it starts in $j_{min}$ – or in communities with low transmission potential – surveillance will underestimate $R$. The spatial structure of the epidemic thus generates something similar to a *core group effect*, an epidemiological concept whereby the disease is initially confined to specific groups where strong clustering, high susceptibility or specific behavioral traits result in transmission rates higher than in the general population – like, for instance, schoolchildren and influenza [43] – and requires specific corrections to epidemic mod-
<figure><img src="image_4.png" /><figcaption>**Fig.2** | Comparison of the reference and measured reproduction ratios. We generated stochastic simulations of a COVID-19-like epidemic in France using a metapopulation model, whereby we divided the French population in 94 spatial communities corresponding to the department of mainland France minus Corsica, and modeled within- and between-community mixing using mobility data by Meta. See Epidemic simulations and Supplementary Methods Section 1.4 for details on how we built *R* from mobility data and on the epidemic spread model. Two initial conditions are tested: in **a**, **c** and **e** 100 initial infections are seeded in the department of Paris, while in **b**, **d** and **f** 100 initial infections are seeded among departments proportionally to their population. **a** and **b** show the daily number of incident infections per 100,000 individuals. *R* (yellow) is the reference reproduction ratio, *R*<sup>estim</sup> (red) is the estimated reproduction ratio as described in the paper, *r*<sub>i</sub><sup>estim</sup> are the locally estimated reproduction ratios, and *EpiEstim* (purple) is the reproduction ratio computed from incidence data using the *EpiEstim* package [14] (see Supplementary Methods Section 1.4 for details). We run 10,000 runs of the stochastic model, and plot medians with solid lines and 95% confidence intervals with shaded areas. For the *r*<sub>i</sub><sup>estim</sup>, confidence intervals are omitted for readability.</figcaption></figure>

els [5, 47, 10, 36]. In our case, areas with higher-than-average ρᵢ play the role of core
groups, skewing the estimated reproduction ratio to values higher than R. We restate,
however, that if the epidemic is instead localized in areas with lower-than-average ρᵢ the
effect is opposite and leads to an underestimate of R, something that is not typically con-
sidered in the standard core-group theory. And just as core group effects may make it
harder to assess and predict the population-level impact of an epidemic, the bias that the spatial structure induces on the reproduction ratio prevents an accurate and timely estimate of the epidemic trajectory, possibly misinforming public health response. The next section will present a method to overcome this. We remark, however, that, at present, our formalism covers only core-groups effects related to the spatial structure. Other drivers, like the presence of clusters in the contact network, require additional complexity: from models featuring an initial subexponential growth of infections [10], to some including the explicit contact network [36] inside each spatial community. Both can be future extensions of the model presented here, but they require changing the functional form of equation (1) in the former case, or forgoing the analytical approach altogether in the latter case. Another known source of bias on the reproduction ratio is the importation of infections from outside of the system, e.g, from other countries [44]. Corrections exist for the initial overestimation coming from miscounting imported infections as infections generated within the systems, and could be used in combination with ours [40]. Our correction, however, also addresses the impact that imported infections may have in changing the spatial distribution of the epidemic. Indeed if importations are distributed differently from the equilibrium distribution v, and there is no reason they should, they will enhance the bias on the reproduction ratio arising from x being different from v, which is the object of our study.

## Correction to surveillance data

So far we have proven that surveillance-based estimates of the reproduction ratio may be biased. We will now propose a way to correct for this bias. The global incidence of infections, used to estimate the reproduction ratio, is trivially proportional to the unweighted average of the incidence across communities: $I_{tot}(t) = N(\sum_{i} I_{i}(t)/N)$. We now define a modified incidence that weighs local infections with the Perron dual vector:

$$I_{tot}^{(v)}(t) = N \left( \frac{\sum_{i} v_{i}^{*} I_{i}(t)}{\sum_{i} v_{i}^{*}} \right) = N \sum_{i} v_{i}^{*} I_{i}(t) = N v^{*} I(t). \quad (5)$$

Then, we use it to get a new estimate of the reproduction ratio:

$$R^{corr}(t) = \frac{I_{tot}^{(v)}(t+1)}{I_{tot}^{(v)}(t)} = \frac{v^{*}I(t+1)}{v^{*}I(t)} = \frac{v^{*}RI(t)}{v^{*}I(t)} = R\frac{v^{*}I(t)}{v^{*}I(t)} = R, \quad (6)$$

and find it is identically equal to the reference value. The practical advantage for epidemic monitoring is clear: our correction gives an unbiased estimate of the reproduction ratio from surveillance data all along the epidemic wave, unlike traditional measures. It has,
however, two potential drawbacks. The former is that if the outbreak starts in communities where $v_i^*$ is small, then $v^*I(t)$ will be very small: stochastic fluctuations would then cause large changes in $R^{corr}$. In that case then $R^{corr}$ may well be accurate, but not precise. However, no initial condition can be orthogonal to $v^*$ whose entries are strictly positive, so even if $v^*x$ is initially small, it is likely to increase quickly and with it the precision of the measurement. In Fig. 3 we show that $R^{corr}$ accurately measures the reference reproduction ratio from the beginning of the epidemic wave, in the case of the simulated epidemics of Fig. 2. Notably, Fig. 3 also shows that if you feed $I_{tot}^{(v)}(t)$ to *EpiEstim* instead of $I_{tot}(t)$ you will also completely remove the bias on the estimate of the reproduction ratio. Our proposed modified incidence can then be readily incorporated to standard tools for public health surveillance, to improve their accuracy.

<figure><img src="image_5.png" /><figcaption>**Fig.3** | Corrected reproduction ratio. This figure uses the same stochastic epidemic model as in Fig. 2 and the same initial conditions, i.e., panels a and b correspond to panels c and d in Fig. 2 respectively. Medians and 95% confidence intervals are computed and shown over 10,000 runs. We compare the standard reproduction ratio measured from surveillance data – both in the Galton-Watson formalism (red) and with *EpiEstim* (purple) on incidence data – and the corresponding $v^*$-corrected estimates: Galton-Watson formalism (blue) and with *EpiEstim* (light blue) on modified incidence data $I_{tot}^{(v)}$.</figcaption></figure>

The latter potential drawback is that our correction requires knowing $v^*$. We argue, however, that this does not require knowing or measuring $R$ in real time and that a good estimate of $v^*$ for epidemic monitoring can be computed during *peace time*, from past population and mobility data (pre-epidemic, or from data collected during earlier epidemic phases). Indeed $v^*$ is more stable than $R$ for the fact that any change happening homogeneously across communities (e.g., changes in the rate of immunity, public health interventions) would change the latter, not the former. Fig. 4 compares the standard estimated reproduction ratio of COVID-19 in France between late 2020 and March 2021 to our correction. The former is computed with *EpiEstim* on inferred case incidence, the latter is computed with *EpiEstim* on the corrected incidence $I_{tot}^{(v)}$, with $v^*$ computed from various samples of past mobility data, all yielding stable estimates. Our correction seems
to point to the fact that traditional surveillance underestimated the reproduction ratio of
COVID-19 in France during January and February 2021. Notably, surveillance recorded
a lower-than-one reproduction ratio during more than two weeks (see also official reports
from that time [22]), indicating a subsiding epidemic wave. This is at odds with what we
know happened: a growing epidemic wave – the French *third wave* – that led to a national
lockdown, enforced on April 3 2021, i.e., immediately after the time window depicted in
Fig. 4. Our corrected reproduction ratio would have instead consistently signaled a grow-
ing epidemic wave throughout the first three months of 2021m by correctly weighing the
infections in departments with low transmission potential – see Supplementary Fig. S12
and S13. This discrepancy carries great significance when put into the context of the de-
bate over public health response at that time. In early 2021 a national curfew was in effect
but cases were rising due to the introduction and gradual takeover of the Alpha variant
of SARS-CoV-2. Authorities were wary of additional restrictions and were relying on
mass vaccination despite models suggesting that it might not be enough [18] – only 3%
of the population had received at least one dose by mid February [22] (week 6 of 2021
in Fig. 4). It is conceivable, albeit circumstantial, that the fact that surveillance underes-
timated the severity of the wave could have contributed to delaying the enforcement of
stricter movement restrictions, which became anyway inevitable later in April.

<figure><img src="image_6.png" /><figcaption>**Fig.4**|Application of the proposed correction to COVID-19 data in France. The pink curve reports the reproduction ratio of COVID-19 in France obtained from country-level incidence, in the same time interval as in Fig. 1b, using the package *EpiEstim*. The blue curves reports the corrected reproduction ratio obtained by feeding *EpiEstim* the corrected incidence as in equation (5). For this, the Perron dual vector **v**\* is reconstructed from mobility data at different times prior the observation period, specified in the legend as week ranges of year 2020. The shaded areas are the 95% credibility intervals computed by *EpiEstim*. For readability, only those of the pink curve and the darkest blue curve are shown. The different blue curves are hardly distinguishable because they overlap. See Reconstruction of the reproduction operator from data for a detailed explanation on how to compute **v**\* from data.</figcaption></figure>

We remark that our method removes a source of bias that it is inherent to the spatial
dynamics of the system and may be present even when surveillance data are accurate,
representative and unbiased. In this sense it should integrate methodologies that address
other sources of bias such as uneven susceptibility and case detection rate, reporting de-
lays, inaccurate reconstruction of infections from proxy timeseries, inaccurate estimation
of the generation time [45, 36, 23, 56, 29]: We prove the computational feasibility of
this integration in our library [1]. It is also important to stress, however, that our cor-
rection is as good as the data used to calculate it, and specifically the data used to re-
construct R. The different global health threats we identified in the introduction entail
different minimum requirements on the data: COVID-19 and influenza need data on local
and between-community mixing, arboviral diseases (dengue, chikungunya, Zika) data on
where residents of each community spend their time and the relative spatial distribution
of the vector population (Aedes mosquitoes). Finally, spreading routes and disease nat-
ural history will inform the specific data needs in the case of newly emerging pathogens
or strains. Overall, however, large-scale estimates of human mobility and mixing will
be necessary, and they increasingly come from mobile phone data. The impact of their
own limitations [68, 62, 67, 33, 55] will require case-by-case assessment when these data
are used to compute our correction. In our case study, Fig. 4 shows that the correction
is robust to choosing human mixing data from different time periods in our case study.
Also, both the data provider and we implemented adjustments to improve the accuracy
and representativeness of the colocation data (see Ref. [28] and Reconstruction of the
reproduction operator from data). Finally, our correction is currently static and works
over periods of fairly constant R: future studies should address how v* should change in
time to consistently correct for the bias on the reproduction ratio, by explicitly coupling
population mixing and disease natural history and handling exogenous shifts in R due to
seasonal trends, behavioral changes and interventions.

Our study describes a practicable way to improve the accuracy of the information
that flows from epidemiological surveillance to public health policymakers. And better
information may lead to more effective policies for preventing and controlling epidemic
threats.
# Methods

## Homogeneous transmission potential

Let us rewrite the estimated reproduction ratio in matrix form as

$$R^{estim}(t) = \frac{I_{tot}(t+1)}{I_{tot}(t)} = \mathbf{F}^T \mathbf{R} \mathbf{x}(t), \qquad (7)$$

where we introduced $\mathbf{F}$ as the unit column vector ($F_i = 1 \forall i$). If we assume that $\mathbf{v}^* = \mathbf{F}^T$ ($v_i^* = 1$) then we can apply $\mathbf{R}$ leftwards in equation (7) and get $R^{estim}(t) = R$ at any time and for any spatial distribution $\mathbf{x}$. Now, the requirement $\mathbf{v}^* = \mathbf{F}^T$ imposes that $\mathbf{R}$ is proportional to a left-stochastic matrix: indeed $\mathbf{F}^T \mathbf{R} = R \mathbf{F}^T$ means $\rho_i = \sum_j R_{ji} = R$, so that each column sums to $R$.

## Calculation of $\Delta(t)$: proof of equation (3)

Combining equation (1) and equation (7) we get the time evolution of the estimated reproduction ratio:

$$R^{estim}(t) = \frac{\mathbf{F}^T \mathbf{R}^{t+1} \mathbf{x}(0)}{\mathbf{F}^T \mathbf{R}^t \mathbf{x}(0)}. \qquad (8)$$

We insert this into equation (2) and get

$$\Delta(t) = \frac{1 \mathbf{F}^T (\mathbf{R} - \mathbf{R}) \mathbf{R}^t \mathbf{x}(0)}{\mathbf{F}^T \mathbf{R}^t \mathbf{x}(0)}. \qquad (9)$$

We introduce the eigenvectors of $\mathbf{R}$ (other than $\mathbf{v}$): $\mathbf{w}_\alpha$ eigenvector with corresponding eigenvalue $\Lambda_\alpha$. Analogously we define the corresponding dual vector $\mathbf{w}_\alpha^*$. Then, we decompose $\mathbf{F}^T$ in the dual basis: $\mathbf{F}^T = \mathbf{v}^* + \sum_\alpha (\mathbf{F}^T \mathbf{w}_\alpha) \mathbf{w}_\alpha^*$. Using this decomposition in equation (9) and applying $\mathbf{R}$ leftwards on the dual eigenvectors we get

$$\Delta(t) = - \frac{\mathbf{F}^T \left[ \sum_\alpha \left(\frac{\Lambda_\alpha}{R}\right)^t \left(1 - \frac{\Lambda_\alpha}{R}\right) \mathbf{w}_\alpha \mathbf{w}_\alpha^* \right] \mathbf{x}(0)}{\mathbf{F}^T \left[ \mathbf{v}\mathbf{v}^* + \sum_\alpha \left(\frac{\Lambda_\alpha}{R}\right)^t \mathbf{w}_\alpha \mathbf{w}_\alpha^* \right] \mathbf{x}(0)}. \qquad (10)$$

The denominator is $C(t)$ in equation (3):

$$C(t) = \frac{1}{\mathbf{F}^T \left[ \mathbf{v}\mathbf{v}^* + \sum_\alpha \left(\frac{\Lambda_\alpha}{R}\right)^t \mathbf{w}_\alpha \mathbf{w}_\alpha^* \right] \mathbf{x}(0)}. \qquad (11)$$

$C(t)$ is always strictly positive because it is proportional to $\mathbf{F}^T \mathbf{R}^t \mathbf{x}(0)$ and tends to $\mathbf{v}^* \mathbf{x}(0)$, i.e., the component of the initial condition onto the eigenspace of the Perron eigenvalue. This component is always nonzero because no $\mathbf{x}(0)$ is nonnegative (as it
is a spatial distribution of infections) and no nonnegative vector can be orthogonal to a strictly positive vector. It is thus the numerator which gives the trend and sign of $\Delta(t)$. Equation (10) then gives the value of the factors $z_\alpha$ in equation (3):

$$z_{\alpha} = -\mathbf{F}^T (\mathbf{w}_{\alpha} \mathbf{w}_{\alpha}^*) \mathbf{x}(0) \qquad (12)$$

In the case of degenerate eigenvalue one should simply replace $\mathbf{w}_\alpha \mathbf{w}_\alpha^*$ with the appropriate projector over the whole eigenspace. Note that, as discussed before, the denominator in equation (10) is always real and positive so any complex phase of $z_\alpha$ must arise from $\Lambda_\alpha$ and $\mathbf{w}_\alpha \mathbf{w}_\alpha^*$.

### Calculation of $\Delta(t)$: $\tau_\alpha, T_\alpha$

We isolate in equation (3) the contribution of each mode $M_\alpha(t)$: $\Delta(t) = C(t) \sum_\alpha M_\alpha(t)$, where

$$M_{\alpha}(t) = M_{\alpha}(0) \left(\frac{\Lambda_{\alpha}}{R}\right)^{t} = M_{\alpha}(0) \left(\frac{|\Lambda_{\alpha}|}{R}\right)^{t} e^{i\theta_{\alpha}t} = M_{\alpha}(0)e^{-t/\tau_{\alpha}}e^{i\theta_{\alpha}t}, \quad (13)$$

where we used the definition of $\tau_\alpha$ given in the main text. The decaying term with characteristic time $\tau_\alpha$ is visible.

If $\Lambda_\alpha$ is real and positive then $\theta_\alpha = 0$ and the oscillating term vanishes. If $\Lambda_\alpha$ is real and negative then $\theta_\alpha = \pi$ and the oscillating term becomes an alternating sign: $e^{i\theta_\alpha t} = (-1)^t$. This is an oscillation with period $T_\alpha = 2$, which is compatible with the definition of $T_\alpha$ given in the main text. Finally, if $\Lambda_\alpha \notin \mathbb{R}$, then then $\bar{\Lambda}_\alpha$ is also an eigenvalue, where the bar denotes complex conjugation. We will call $\bar{\alpha}$ the index corresponding to that eigenvalue: $\Lambda_{\bar{\alpha}} = \bar{\Lambda}_\alpha$. Also, the projector over the eigenspace of $\Lambda_{\bar{\alpha}}$ is the elementwise complex conjugate of the projector over the eigenspace of $\Lambda_\alpha$, meaning that $z_{\bar{\alpha}} = \bar{z}_\alpha$, and thus $M_{\bar{\alpha}} = \bar{M}_\alpha(0)$. Then $\alpha, \bar{\alpha}$ contribute in pair, as follows:

$$\begin{aligned} M_{\alpha}(t) + M_{\bar{\alpha}}(t) &= e^{-t/\tau_{\alpha}} [M_{\alpha}(0)e^{i\theta_{\alpha}t} + M_{\bar{\alpha}}(0)e^{-i\theta_{\alpha}t}] \\ &= 2e^{-t/\tau_{\alpha}} |M_{\alpha}(0)| \operatorname{Re} e^{i\theta_{\alpha}t + \phi_{\alpha}} = 2e^{-t/\tau_{\alpha}} \cos\left(\frac{2\pi}{T_{\alpha}}t + \phi_{\alpha}\right). \end{aligned} \quad (14)$$

Here we used the definition of $T_\alpha$ given in the main text, explicitly showing the emergence of the oscillating term with period $T_\alpha$.

### Local over- and underestimation of $R$

The Collatz-Wielandt inequalities tell us that, for any spatial distribution of infections $x$, $\min_{i|x_i \neq 0}(\mathbf{R}x)_i/x_i \le R$ and $\max_{i|x_i \neq 0}(\mathbf{R}x)_i/x_i \ge R$. Given that $r_i^{estim} = (\mathbf{R}x)_i/x_i$, out
of equilibrium there will always be at least one community overestimating the reference reproduction ratio ($r_i^{estim} > R$) and one underestimating it ($r_i^{estim} < R$).

## Reconstruction of the reproduction operator from data

The main data used for the reconstruction of reproduction operators for mainland France are Meta Colocation Maps, which are calculated from mobile phone data and show the proportion of time that residents of different communities spend in proximity [28]. More precisely, they give the probability $p_{ij}$ that a randomly chosen person that is resident of community $i$ and a randomly chosen person resident of community $j$ are both located in a same 600m × 600m square, during a randomly chosen five-minutes time window, in a given week. The data were provided at the resolution of departments (ADM 2). To reconstruct $\mathbf{R}$ from these data, we assumed that the expected number of secondary infections generated among the residents of community $i$, by a case who is resident of community $j$, is given by $R_{ij} = C p_{ij} n_i$, where $n_i$ is the population of spatial patch $i$, and $C$ is an overall transmissibility parameter. Notably, while the value of the spectral radius of $\mathbf{R}$ clearly depends on $C$, the left and right Perron eigenvectors $\mathbf{v}$ and $\mathbf{v}^*$ do not, and depend solely on the data. To compute $\mathbf{v}$, $\mathbf{v}^*$ from $\mathbf{R}$ we computed its right ($\mathbf{v}$) and left ($\mathbf{v}^*$) eigenvectors associated to the largest eigenvalue, using the Python function ` numpy.lINALG.eig` of ` numpy` version 1.12. The diagonal elements $p_{ii}$ of the Colocation Maps quantify the mixing within each community. In a previous study we found that they overestimate it in densely populated areas and when people spend a lot of time at home [39]. Using that same methodology, we thus corrected those entries using Meta Movement Range Maps (see Data availability). Movement Range Maps give the daily fraction $sp_i$ of residents of community $i$ that do not leave their 600m × 600m home tile. The probability $p_{ii}^{(sp)}$ to detect co-locations among people staying in their separate homes and thus not actually mixing is, in community $i$,

$$p_{ii}^{(sp)} = \frac{(sp_i d_i A)(sp_i d_i A - 1)m_i}{n_i(n_i - 1)}, \quad (15)$$

where $d_i$ is the population density in community $i$, $A = 0.36 \text{ km}^2$ is the area of a single tile, $m_i$ is the number of tiles occupied by community $i$ and $n_i$ is the resident population population. We then discount $p_{ii}^{(sp)}$ from $p_{ii}$ to remove spurious colocations. We show in Supplementary Fig. S2 the magnitude of this for the central week of the gray period of Fig. 1. For most departments it ranges between 1% and 10%, exceeding the latter value only in the densely populated area of Paris and its surroundings.
## Epidemic simulations

The model of epidemic spread used in simulations is a stochastic discrete-time metapopulation model whereby spatially distinct communities are linked through mobility [13, 3, 25, 9]. We use a synthetic population based on census data from the National Institute of Statistics and Economic Studies (INSEE) in France. We divide this population in 94 spatial communities corresponding to the departments of mainland France except Corsica. Meta Colocation Maps [28] and Movement Range Maps are used to reconstruct the coupling $p_{ij}$ between communities $i$ and $j$ and the within-community $i$ mixing $p_{ii}$. We use a SEIR compartmental model with $\epsilon$ as the rate of the E-to-I transition, $\mu$ rate of the I-to-R transition and $\beta$ overall transmission rate. The values of $\epsilon, \mu$ come from Ref. [21] (COVID-19-like), and $\beta$ is set so that $R = 1.5$. Following Ref. [28], the rate of force-of-infection to which those in community $i$ are subjected is $\beta \sum_{j} p_{ij} I_{j}$, where $I_{j}$ is the number of infected in $j$. We use a discretization time step $\Delta d = 1day$. The Supplementary Methods Section 1.4. describes in detail each step of the stochastic simulation model. There, we also compute the reproduction ratio for our model according to the next generation matrix formalism [20], obtaining

$$NGM_{ij} = \frac{\beta}{\mu} p_{ij} n_{i}, \qquad (16)$$

and $R = \text{spectral rad}(NGM)$ – see Supplementary Methods Section 1.4. Again, $n_{i}$ is the population of community $i$.

## Data availability

Meta Colocation Maps, which were used to reconstruct reproduction operators and to infer between- and within-community mixing for stochastic simulations can be requested at https://dataforgood.facebook.com/dfg/tools/colocation-maps. Meta Movement Range Maps used to correct within-community colocations are available at https://data.humdata.org/dataset/movement-range-maps. Hospital admission data in France are available at https://www.data.gouv.fr. French census data can be found at https://www.insee.fr. French departments shapefiles are available at https://www.data.gouv.fr/en/datasets/carte-des-departements-2-1/. All websites accessed November 2023.
## Code availability

The code used in this study is available here: https://github.com/ev-modelers/birello-surveillance. A continuously developed and maintained library allowing to implement our correction to estimates of the reproduction ratio from surveillance data can be found at https://github.com/ev-modelers/rt-from-surveillance.

## Acknowledgments

Colocation data were available thanks to *Data For Good at Meta*.

## Author contributions

E.V. conceived and designed the study. E.V. and M.R.F. developed the theory. P.B. developed the code and the publicly available library, performed the numerical simulations and analyzed the results. B.W. analyzed surveillance data and reconstructed infections from them. P.B., M.R.F., B.W., V.C., E.V. interpreted the results. E.V. drafted the Article. P.B., M.R.F., B.W., V.C., E.V. contributed to and approved the final version of the Article.

## Competing interests

The authors declare no competing interests.
References

[1] https://github.com/ev-modelers/rt-from-surveillance.

[2] Introducing a coherent European framework for tuning COVID-19 response measures, March 2021.

[3] Duygu Balcan, Vittoria Colizza, Bruno Gonçalves, Hao Hu, José J. Ramasco, and Alessandro Vespignani. Multiscale mobility networks and the spatial spreading of infectious diseases. *Proceedings of the National Academy of Sciences*, 106(51):21484–21489, December 2009. Publisher: Proceedings of the National Academy of Sciences.

[4] Duygu Balcan and Alessandro Vespignani. Phase transitions in contagion processes mediated by recurrent mobility patterns. *Nature Physics*, 7(7):581–586, July 2011. Number: 7 Publisher: Nature Publishing Group.

[5] Shweta Bansal, Bryan T Grenfell, and Lauren Ancel Meyers. When individual behaviour matters: homogeneous and network models in epidemiology. *Journal of The Royal Society Interface*, 4(16):879–891, July 2007. Publisher: Royal Society.

[6] Matthew Biggerstaff, Simon Cauchemez, Carrie Reed, Manoj Gambhir, and Lyn Finelli. Estimates of the reproduction number for seasonal, pandemic, and zoonotic influenza: a systematic review of the literature. *BMC Infectious Diseases*, 14(1):480, September 2014.

[7] Colin J. Carlson, Gregory F. Albery, Cory Merow, Christopher H. Trisos, Casey M. Zipfel, Evan A. Eskew, Kevin J. Olival, Noam Ross, and Shweta Bansal. Climate change increases cross-species viral transmission risk. *Nature*, 607(7919):555–562, July 2022. Number: 7919 Publisher: Nature Publishing Group.

[8] Mario Castro, Saúl Ares, José A. Cuesta, and Susanna Manrubia. The turning point and end of an expanding epidemic cannot be precisely forecast. *Proceedings of the National Academy of Sciences*, 117(42):26190–26196, October 2020. Publisher: National Academy of Sciences Section: Biological Sciences.

[9] Serina Chang, Emma Pierson, Pang Wei Koh, Jaline Gerardin, Beth Redbird, David Grusky, and Jure Leskovec. Mobility network models of COVID-19 explain inequities and inform reopening. *Nature*, 589:82–87, November 2020. Publisher: Nature Publishing Group.
[10] Gerardo Chowell, Cécile Viboud, Lone Simonsen, and Seyed M. Moghadas. Characterizing the reproduction number of epidemics with early subexponential growth dynamics. *Journal of The Royal Society Interface*, 13(123):20160659, October 2016. Publisher: Royal Society.

[11] Claudia T. Codeço, Daniel A. M. Villela, and Flavio C. Coelho. Estimating the effective reproduction number of dengue considering temperature-dependent generation intervals. *Epidemics*, 25:101-111, December 2018.

[12] Pietro Coletti, Chiara Poletto, Clément Turbelin, Thierry Blanchon, and Vittoria Colizza. Shifting patterns of seasonal influenza epidemics. *Scientific Reports*, 8(1):12786, August 2018. Number: 1 Publisher: Nature Publishing Group.

[13] Vittoria Colizza and Alessandro Vespignani. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations. *Journal of Theoretical Biology*, 251(3):450-467, 2008.

[14] Anne Cori, Neil M. Ferguson, Christophe Fraser, and Simon Cauchemez. A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics. *American Journal of Epidemiology*, 178(9):1505-1512, November 2013.

[15] Bahman Davoudi, Joel C. Miller, Rafael Meza, Lauren Ancel Meyers, David J. D. Earn, and Babak Pourbohloul. Early Real-Time Estimation of the Basic Reproduction Number of Emerging Infectious Diseases. *Physical Review X*, 2(3):031005, July 2012. Publisher: American Physical Society.

[16] Giulia de Meijere, Eugenio Valdano, Claudio Castellano, Marion Debin, Charly Kengne-Kuetche, Clément Turbelin, Harold Noël, Joshua S. Weitz, Daniela Paolotti, Lisa Hermans, Niel Hens, and Vittoria Colizza. Attitudes towards booster, testing and isolation, and their impact on COVID-19 response in winter 2022/2023 in France, Belgium, and Italy: a cross-sectional survey and modelling study. *The Lancet Regional Health - Europe*, 28:100614, May 2023.

[17] Ranu S. Dhillon, Devabhaktuni Srikrishna, and Gerardo Chowell. Getting to zero in the DR Congo Ebola outbreak. *The Lancet Infectious Diseases*, 20(4):395-397, April 2020. Publisher: Elsevier.

[18] Laura Di Domenico, Chiara E. Sabbatini, Giulia Pullano, Daniel Lévy-Bruhl, and Vittoria Colizza. Impact of January 2021 curfew measures on SARS-CoV-2 B.1.1.7
circulation in France. *Eurosurveillance*, 26(15):2100272, April 2021. Publisher: European Centre for Disease Prevention and Control.

[19] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. *Journal of Mathematical Biology*, 28(4):365–382, June 1990.

[20] Odo Diekmann, J.A.P. Heesterbeek, and M.G. Roberts. The construction of next-generation matrices for compartmental epidemic models. *Journal of the Royal Society, Interface / the Royal Society*, 7:873–85, 11 2009.

[21] Benjamin Faucher, Rania Assab, Jonathan Roux, Daniel Levy-Bruhl, Cécile Kiem, Simon Cauchemez, Laura Zanetti, Vittoria Colizza, Pierre-Yves Boelle, and Chiara Poletto. Agent-based modelling of reactive vaccination of workplaces and schools against covid-19. *Nature Communications*, 13:1414, 03 2022.

[22] Santé Publique France. COVID-19 : point épidémiologique du 11 février 2021.

[23] Katelyn M. Gostic, Lauren McGough, Edward B. Baskerville, Sam Abbott, Keya Joshi, Christine Tedijanto, Rebecca Kahn, Rene Niehus, James A. Hay, Pablo M. De Salazar, Joel Hellewell, Sophie Meakin, James D. Munday, Nikos I. Bosse, Katharine Sherrat, Robin N. Thompson, Laura F. White, Jana S. Huisman, Jérémie Scire, Sebastian Bonhoeffer, Tanja Stadler, Jacco Wallinga, Sebastian Funk, Marc Lipsitch, and Sarah Cobey. Practical considerations for measuring the effective reproductive number, Rt. *PLOS Computational Biology*, 16(12):e1008409, 2020. Publisher: Public Library of Science.

[24] Fiona M. Guerra, Shelly Bolotin, Gillian Lim, Jane Heffernan, Shelley L. Deeks, Ye Li, and Natasha S. Crowcroft. The basic reproduction number (R0) of measles: a systematic review. *The Lancet Infectious Diseases*, 17(12):e420–e428, December 2017. Publisher: Elsevier.

[25] J. Gómez-Gardeñes, D. Soriano-Paños, and A. Arenas. Critical regimes driven by recurrent mobility patterns of reaction-diffusion processes in networks. *Nature Physics*, 14(4):391–395, April 2018. Number: 4 Publisher: Nature Publishing Group.

[26] Joel Hellewell, Sam Abbott, Amy Gimma, Nikos I Bosse, Christopher I Jarvis, Timothy W Russell, James D Munday, Adam J Kucharski, W John Edmunds, Fiona Sun,
Stefan Flasche, Billy J Quilty, Nicholas Davies, Yang Liu, Samuel Clifford, Petra Klepac, Mark Jit, Charlie Diamond, Hamish Gibbs, Kevin van Zandvoort, Sebastian Funk, and Rosalind M Eggo. Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts. *The Lancet Global Health*, 8(4):e488–e496, April 2020.

[27] L. Hufnagel, D. Brockmann, and T. Geisel. Forecast and control of epidemics in a globalized world. *Proceedings of the National Academy of Sciences*, 101(42):15124–15129, October 2004. Publisher: Proceedings of the National Academy of Sciences.

[28] Shankar Iyer, Brian Karrer, Daniel T. Citron, Farshad Kooti, Paige Maas, Zeyu Wang, Eugenia Giraudy, Ahmed Medhat, P. Alex Dow, and Alex Pompe. Large-scale measurement of aggregate human colocation patterns for epidemiological modeling. *Epidemics*, 42:100663, 2023.

[29] D. C. P. Jorge, J. F. Oliveira, J. G. V. Miranda, R. F. S. Andrade, and S. T. R. Pinho. Estimating the effective reproduction number for heterogeneous models using incidence data. *Royal Society Open Science*, 9(9):220005, September 2022. Publisher: Royal Society.

[30] Frédéric Jourdain, David Roiz, Henriette de Valk, Harold Noël, Grégory L'Ambert, Florian Franke, Marie-Claire Paty, Anne Guinard, Jean-Claude Desenclos, and Benjamin Roche. From importation to autochthonous transmission: Drivers of chikungunya and dengue emergence in a temperate area. *PLOS Neglected Tropical Diseases*, 14(5):e0008320, 2020. Publisher: Public Library of Science.

[31] Matt J. Keeling and P. Rohani. *Modeling Infectious Diseases in Humans and Animals*. Princeton University Press, Princeton, NJ, USA, 2007.

[32] Adam J Kucharski, Rosalind M Eggo, Conall H Watson, Anton Camacho, Sebastian Funk, and W John Edmunds. Effectiveness of ring vaccination as control strategy for ebola virus disease. *Emerging infectious diseases*, 22(1):105, 2016.

[33] Shengjie Lai, Elisabeth zu Erbach-Schoenberg, Carla Pezzulo, Nick W. Ruktanonchai, Alessandro Sorichetta, Jessica Steele, Tracey Li, Claire A. Dooley, and Andrew J. Tatem. Exploring the use of mobile phone data for national migration statistics. *Palgrave Communications*, 5(1):1–10, March 2019. Number: 1 Publisher: Palgrave.
[34] You Li, Harry Campbell, Durga Kulkarni, Alice Harpur, Madhurima Nundy, Xin Wang, and Harish Nair. The temporal association of introducing and lifting non-pharmaceutical interventions with the time-varying reproduction number (R) of SARS-CoV-2: a modelling study across 131 countries. *The Lancet Infectious Diseases*, 21(2):193–202, February 2021. Publisher: Elsevier.

[35] You Li and Harish Nair. Trends in the global burden of lower respiratory infections: the knowns and the unknowns. *The Lancet Infectious Diseases*, 22(11):1523–1525, November 2022. Publisher: Elsevier.

[36] Quan-Hui Liu, Marco Ajelli, Alberto Aleta, Stefano Merler, Yamir Moreno, and Alessandro Vespignani. Measurability of the epidemic reproduction number in data-driven contact networks. *Proceedings of the National Academy of Sciences*, 115(50):12680–12685, December 2018. Publisher: Proceedings of the National Academy of Sciences.

[37] J. O. Lloyd-Smith, S. J. Schreiber, P. E. Kopp, and W. M. Getz. Superspreading and the effect of individual variation on disease emergence. *Nature*, 438(7066):355–359, November 2005. Number: 7066 Publisher: Nature Publishing Group.

[38] Gaël D. Maganga, Jimmy Kapetshi, Nicolas Berthet, Benoît Kebela Ilunga, Felix Kabange, Placide Mbala Kingebeni, Vital Mondonge, Jean-Jacques T. Muyembe, Eric Bertherat, Sylvie Briand, Joseph Cabore, Alain Epelboin, Pierre Formenty, Gary Kobinger, Licé González-Angulo, Ingrid Labouba, Jean-Claude Manuguerra, Jean-Marie Okwo-Bele, Christopher Dye, and Eric M. Leroy. Ebola Virus Disease in the Democratic Republic of Congo. *New England Journal of Medicine*, 371(22):2083–2091, November 2014. Publisher: Massachusetts Medical Society.

[39] Mattia Mazzoli, Eugenio Valdano, and Vittoria Colizza. Projecting the COVID-19 epidemic risk in France for the summer 2021. *Journal of Travel Medicine*, 28(7), October 2021.

[40] G. N. Mercer, K. Glass, and N. G. Becker. Effective reproduction numbers are commonly overestimated early in a disease outbreak. *Statistics in Medicine*, 30(9):984–994, 2011. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.4174.

[41] Jane P. Messina, Oliver J. Brady, Nick Golding, Moritz U. G. Kraemer, G. R. William Wint, Sarah E. Ray, David M. Pigott, Freya M. Shearer, Kimberly Johnson, Lucas Earl, Laurie B. Marczak, Shreya Shirude, Nicole Davis Weaver,
Marius Gilbert, Raman Velayudhan, Peter Jones, Thomas Jaenisch, Thomas W. Scott, Robert C. Reiner, and Simon I. Hay. The current and future global distribution and population at risk of dengue. *Nature Microbiology*, 4(9):1508–1515, September 2019. Number: 9 Publisher: Nature Publishing Group.

[42] Zindoga Mukandavire, Shu Liao, Jin Wang, Holly Gaff, David L. Smith, and J. Glenn Morris. Estimating the reproductive numbers for the 2008–2009 cholera outbreaks in Zimbabwe. *Proceedings of the National Academy of Sciences*, 108(21):8767–8772, May 2011.

[43] H. Nishiura, C. Castillo-Chavez, M. Safan, and G. Chowell. Transmission potential of the new influenza A(H1N1) virus and its age-specificity in Japan. *Eurosurveil-lance*, 14(22):19227, June 2009. Publisher: European Centre for Disease Prevention and Control.

[44] H. Nishiura and M. G. Roberts. Estimation of the reproduction number for 2009 pandemic influenza A(H1N1) in the presence of imported cases. *Eurosurveillance*, 15(29):19622, July 2010. Publisher: European Centre for Disease Prevention and Control.

[45] Hiroshi Nishiura. Correcting the Actual Reproduction Number: A Simple Method to Estimate R0 from Early Epidemic Growth Data. *International Journal of Envi-ronmental Research and Public Health*, 7(1):291–302, January 2010.

[46] Hiroshi Nishiura and Gerardo Chowell. The Effective Reproduction Number as a Prelude to Statistical Estimation of Time-Dependent Epidemic Trends. In Gerardo Chowell, James M. Hyman, Luís M. A. Bettencourt, and Carlos Castillo-Chavez, ed-itors, *Mathematical and Statistical Estimation Approaches in Epidemiology*, pages 103–121. Springer Netherlands, Dordrecht, 2009.

[47] Hiroshi Nishiura, Gerardo Chowell, Muntaser Safan, and Carlos Castillo-Chavez. Pros and cons of estimating the reproduction number from early epidemic growth rate of influenza A (H1N1) 2009. *Theoretical Biology and Medical Modelling*, 7(1):1, January 2010.

[48] Thomas Obadia, Romana Haneef, and Pierre-Yves Boëlle. The R0 package: a tool-box to estimate reproduction numbers for epidemic outbreaks. *BMC Medical Infor-matics and Decision Making*, 12(1):147, December 2012.
[49] An Pan, Li Liu, Chaolong Wang, Huan Guo, Xingjie Hao, Qi Wang, Jiao Huang, Na He, Hongjie Yu, Xihong Lin, Sheng Wei, and Tangchun Wu. Association of Public Health Interventions With the Epidemiology of the COVID-19 Outbreak in Wuhan, China. *JAMA*, 323(19):1915–1923, May 2020.

[50] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani. Epidemic processes in complex networks. *Reviews of Modern Physics*, 87(3):925–979, 2015.

[51] Benjamin Ridenhour, Jessica M. Kowalik, and David K. Shay. Unraveling R0: Considerations for Public Health Applications. *American Journal of Public Health*, 108(Suppl 6):S445–S454, December 2018.

[52] Marina Romanello, Claudia Di Napoli, Paul Drummond, Carole Green, Harry Kennard, Pete Lampard, Daniel Scamman, Nigel Arnell, Sonja Ayeb-Karlsson, Lea Berrang Ford, Kristine Belesova, Kathryn Bowen, Wenjia Cai, Max Callaghan, Diarmid Campbell-Lendrum, Jonathan Chambers, Kim R. van Daalen, Carole Dalin, Niheer Dasandi, Shouro Dasgupta, Michael Davies, Paula Dominguez-Salas, Robert Dubrow, Kristie L. Ebi, Matthew Eckelman, Paul Ekins, Luis E. Escobar, Lucien Georgeson, Hilary Graham, Samuel H. Gunther, Ian Hamilton, Yun Hang, Risto Hänninen, Stella Hartinger, Kehan He, Jeremy J. Hess, Shih-Che Hsu, Slava Jankin, Louis Jamart, Ollie Jay, Ilan Kelman, Gregor Kiesewetter, Patrick Kinney, Tord Kjellstrom, Dominic Kniveton, Jason K. W. Lee, Bruno Lemke, Yang Liu, Zhao Liu, Melissa Lott, Martin Lotto Batista, Rachel Lowe, Frances MacGuire, Maquins Odhiambo Sewe, Jaime Martinez-Urtaza, Mark Maslin, Lucy McAllister, Alice McGushin, Celia McMichael, Zhifu Mi, James Milner, Kelton Minor, Jan C. Minx, Nahid Mohajeri, Maziar Moradi-Lakeh, Karyn Morrissey, Simon Munzert, Kris A. Murray, Tara Neville, Maria Nilsson, Nick Obradovich, Megan B. O’Hare, Tadj Oreszczyn, Matthias Otto, Fereidoon Owfi, Olivia Pearman, Mahnaz Rabbaniha, Elizabeth J. Z. Robinson, Joacim Rocklöv, Renee N. Salas, Jan C. Semenza, Jodi D. Sherman, Liuhua Shi, Joy Shumake-Guillemot, Grant Silbert, Mikhail Sofiev, Marco Springmann, Jennifer Stowell, Meisam Tabatabaei, Jonathon Taylor, Joaquin Triñanes, Fabian Wagner, Paul Wilkinson, Matthew Winning, Marisol Yglesias-González, Shihui Zhang, Peng Gong, Hugh Montgomery, and Anthony Costello. The 2022 report of the Lancet Countdown on health and climate change: health at the mercy of fossil fuels. *The Lancet*, 400(10363):1619–1654, November 2022. Publisher: Elsevier.
[53] Isobel Routledge, José Eduardo Romero Chevéz, Zulma M. Cucunubá, Manuel Gomez Rodriguez, Caterina Guinovart, Kyle B. Gustafson, Kammerle Schneider, Patrick G. T. Walker, Azra C. Ghani, and Samir Bhatt. Estimating spatiotemporally varying malaria reproduction numbers in a near elimination setting. *Nature Communications*, 9(1):2476, June 2018. Number: 1 Publisher: Nature Publishing Group.

[54] Samuel V. Scarpino and Giovanni Petri. On the predictability of infectious disease outbreaks. *Nature Communications*, 10(1):898, February 2019. Number: 1 Publisher: Nature Publishing Group.

[55] Vedran Sekara, Elisa Omodei, Laura Healy, Jan Beise, Claus Hansen, Danzhen You, Saskia Blume, and Manuel Garcia-Herranz. Mobile Phone Data for Children on the Move: Challenges and Opportunities. In Albert Ali Salah, Alex Pentland, Bruno Lepri, and Emmanuel Letouzé, editors, *Guide to Mobile Data Analytics in Refugee Scenarios: The 'Data for Refugees Challenge' Study*, pages 53–66. Springer International Publishing, Cham, 2019.

[56] Katharine Sherratt, Sam Abbott, Sophie R. Meakin, Joel Hellewell, James D. Munday, Nikos Bosse, null null, Mark Jit, and Sebastian Funk. Exploring surveillance data biases when estimating the reproduction number: with insights into subpopulation transmission of COVID-19 in England. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 376(1829):20200283, May 2021. Publisher: Royal Society.

[57] D. Soriano-Paños, L. Lotero, A. Arenas, and J. Gómez-Gardeñes. Spreading Processes in Multiplex Metapopulations Containing Different Mobility Networks. *Physical Review X*, 8(3):031039, August 2018. Publisher: American Physical Society.

[58] Zachary Susswein, Eugenio Valdano, Tobias Brett, Pejman Rohani, Vittoria Colizza, and Shweta Bansal. Ignoring spatial heterogeneity in drivers of sars-cov-2 transmission in the us will impede sustained elimination. *medRxiv*, 2021.

[59] R. N. Thompson, J. E. Stockwin, R. D. van Gaalen, J. A. Polonsky, Z. N. Kamvar, P. A. Demarsh, E. Dahlqwist, S. Li, E. Miguel, T. Jombart, J. Lessler, S. Cauchemez, and A. Cori. Improved inference of time-varying reproduction numbers during infectious disease outbreaks. *Epidemics*, 29:100356, December 2019.
[60] Robin N. Thompson, Christopher A. Gilligan, and Nik J. Cunniffe. Control fast or control smart: When should invading pathogens be controlled? *PLOS Computational Biology*, 14(2):e1006014, February 2018. Publisher: Public Library of Science.

[61] Rory Thompson, James Gareth Wood, Stefano Tempia, and David John Muscatello. Global variation in early epidemic growth rates and reproduction number of seasonal influenza. *International Journal of Infectious Diseases*, 122:382–388, September 2022.

[62] Michele Tizzoni, Paolo Bajardi, Adeline Decuyper, Guillaume Kon Kam King, Christian M. Schneider, Vincent Blondel, Zbigniew Smoreda, Marta C. González, and Vittoria Colizza. On the Use of Human Mobility Proxies for Modeling Epidemics. *PLOS Computational Biology*, 10(7):e1003716, 2014. Publisher: Public Library of Science.

[63] Cristiano Trevisin, Enrico Bertuzzo, Damiano Pasetto, Lorenzo Mari, Stefano Miccoli, Renato Casagrandi, Marino Gatto, and Andrea Rinaldo. Spatially explicit effective reproduction numbers from incidence and mobility data. *Proceedings of the National Academy of Sciences*, 120(20):e2219816120, May 2023. Publisher: Proceedings of the National Academy of Sciences.

[64] J Wallinga and M Lipsitch. How generation intervals shape the relationship between growth rates and reproductive numbers. *Proceedings of the Royal Society B: Biological Sciences*, 274(1609):599–604, February 2007.

[65] Jacco Wallinga, Michiel van Boven, and Marc Lipsitch. Optimizing infectious disease interventions during an emerging epidemic. *Proceedings of the National Academy of Sciences*, 107(2):923–928, January 2010. Publisher: Proceedings of the National Academy of Sciences.

[66] H. W. Watson and Francis Galton. On the probability of the extinction of families. *The Journal of the Anthropological Institute of Great Britain and Ireland*, 4:138–144, 1875.

[67] Amy Wesolowski, Caroline O. Buckee, Kenth Engø-Monsen, and C. J. E. Metcalf. Connecting Mobility to Infectious Diseases: The Promise and Limits of Mobile Phone Data. *The Journal of Infectious Diseases*, 214(suppl 4):S414–S420, December 2016.
[68] Amy Wesolowski, Nathan Eagle, Abdisalan M. Noor, Robert W. Snow, and Caroline O. Buckee. The impact of biases in mobile phone ownership on estimates of human mobility. *Journal of The Royal Society Interface*, 10(81):20120986, April 2013. Publisher: Royal Society.

[69] Laura F. White, Brett Archer, and Marcello Pagano. Estimating the reproductive number in the presence of spatial heterogeneity of transmission patterns. *International Journal of Health Geographics*, 12(1):35, July 2013.
# Contents

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Supplementary Methods</td><td>2</td></tr><tr><td>1.1</td><td>On the applicability of the Perron-Frobenius theorem</td><td>2</td></tr><tr><td>1.2</td><td>Local estimated reproduction ratio</td><td>3</td></tr><tr><td>1.2.1</td><td>Convergence</td><td>3</td></tr><tr><td>1.2.2</td><td>Decay modes</td><td>3</td></tr><tr><td>1.3</td><td>Estimation of department-level incident infections</td><td>4</td></tr><tr><td>1.3.1</td><td>Models</td><td>4</td></tr><tr><td>1.3.2</td><td>Comparison to seroprevalence data</td><td>5</td></tr><tr><td>1.4</td><td>The stochastic epidemic model</td><td>6</td></tr><tr><td>1.4.1</td><td>Reference reproduction ratio</td><td>7</td></tr><tr><td>1.4.2</td><td>Generation interval and EpiEstim settings</td><td>8</td></tr><tr><td>1.5</td><td>Oscillations in European countries</td><td>9</td></tr><tr><td>2</td><td>Supplementary Figures</td><td>10</td></tr><tr><td>2.1</td><td>Decay time and oscillation period of the modes of Δ = R<sup>estim</sup> − R<sup>ref</sup> for selected European countries</td><td>10</td></tr><tr><td>2.2</td><td>Relative effect of the diagonal correction to co-location probabilities in French departments</td><td>11</td></tr><tr><td>2.3</td><td>Daily number of new COVID-19-related infections per 100,000 residents</td><td>12</td></tr><tr><td>2.4</td><td>National attack rate: estimates vs serological studies</td><td>14</td></tr><tr><td>2.5</td><td>Age demographic distribution in French regions</td><td>16</td></tr><tr><td>2.6</td><td>Fig. 1b and Fig. 4 of the main paper, using the age-stratified model for the estimate of infections</td><td>17</td></tr><tr><td>2.7</td><td>Comparison of the reference, estimated and corrected reproduction ratios with six different initial conditions for the stochastic model</td><td>18</td></tr><tr><td>2.8</td><td>Comparison between point-wise medians and single runs of the stochastic metapopulation model</td><td>19</td></tr><tr><td>2.9</td><td>Variation of the spatial distribution of infections between weeks 6 and 11 of 2021, and the equilibrium distribution</td><td>21</td></tr><tr><td colspan="2">Supplementary References</td><td>23</td></tr></tbody></table>
# 1 Supplementary Methods

## 1.1 On the applicability of the Perron-Frobenius theorem

The Perron-Frobenius theorem as used in the main paper requires that the matrix be strictly positive ($R_{ij} > 0$), or nonnegative ($R_{ij} \ge 0$) and irreducible. In our case $R$ may not be strictly positive if, for some $i, j$, cases from $i$ generate no cases in $j$, so we shall prove here that it is irreducible, or that it can be made irreducible. A nonnegative matrix is irreducible if and only if its associated directed graph is strongly connected [1]. The associated graph of $R$ is that which has a link between nodes $i, j$ if $R_{ij} > 0$. In general, a suitable permutation of the node indices will bring $R$ to the following form, which mirrors the general bow-tie structure of the associated directed graph:

$$R = \begin{pmatrix} T_u & 0 & 0 \\ B_1 & R_{scc} & 0 \\ B_2 & B_3 & T_d \end{pmatrix}. \qquad (S1)$$

where the blocks $T_u, T_d$ are lower diagonal and $R_{scc}$ is the adjacency submatrix of strongly connected component. The spectrum of $R$ is then the union of the diagonal elements of $T_u, T_d$ and the spectrum of $R_{scc}$. Now three options are possible. First, if $R^{ref}$, the spectral radius of $R$ and reference reproduction ratio, is among the diagonal elements of $T_d$ this means that there is one community that sustains the epidemic and at most exports cases to other sink communities (remember $T_d$ is lower diagonal), so it is a trivial case with no actual epidemic dynamics between communities.

Second, if $R^{ref}$ belongs to the spectrum of $R_{scc}$, then we shall write the Perron eigenvector in blocks as follows:

$$v = \begin{pmatrix} \frac{v_u}{v_{scc}} \\ \frac{v_{scc}}{v_d} \end{pmatrix}. \qquad (S2)$$

If we write by blocks the eigenvector equation $Rv = R^{ref}v$, on the top block we have $T_u v_u = R^{ref} v_u$, whose only solution is $v_u = 0$ as $R^{ref}$ is not an eigenvalue of $T_u$. This means $R_{scc}v_{scc} = R^{ref}v_{scc}$, and $v_d = (R^{ref} - T_d)^{-1} B_3 v_{scc}$, this being nonsingular because $R^{ref}$ is not an eigenvalue of $T_d$. The dynamics is thus completely determined by the strongly connected component, and we can restrict our study to $R_{scc}$, which represents by definition a strongly-connected graph and as such is irreducible, proving our initial claim.

Finally, if $R^{ref}$ is among the diagonal elements of $T_u$ again this means that there is one community that generates cases and exports them, possibly through several steps, to the
strongly-connected component. Again this is seeding part is trivial and underlies no actual epidemic dynamics between communities, so that again we can restrict our study to $R_{scc}$.

## 1.2 Local estimated reproduction ratio

### 1.2.1 Convergence

We prove here that when the spatial distribution of infections is at equilibrium ($x = v$) and $r_i^{estim} = R^{ref}$, this implies a detailed balance between locally generated infections, imported infections and exported infections. First, we define the three quantities as computed from the previous generation:

$$I_i^{LOC}(t) = R_{ii}I_i(t - 1); \qquad (S3)$$

$$I_i^{IN}(t) = \sum_{j \neq i} R_{ij}I_j(t - 1); \qquad (S4)$$

$$I_i^{OUT}(t) = I_i(t - 1) \sum_{j \neq i} R_{ji} = (\rho_i - R_{ii})I_i(t - 1). \qquad (S5)$$

Clearly, $I_i(t) = I_i^{LOC}(t) + I_i^{IN}(t)$. But at the equilibrium is it also true that $I_i(t) = R I_i(t - 1)$. Rewriting it as a function of $I_i^{OUT}(t)$ we can get to the detailed balance (omitting $t$):

$$I_i^{LOC} + I_i^{IN} = \frac{R^{ref}}{\rho_i - R_{ii}} I_i^{OUT}. \qquad (S6)$$

### 1.2.2 Decay modes

It is possible to get an equivalent of equation (4) of the main paper to study the deviation of the local estimated reproduction ratio ($r_i^{estim}$) from $R^{ref}$. Defining $\Delta^{(i)}(t) = (r_i^{estim} - R^{ref})/R^{ref}$ as the relative difference, and proceeding analogously to $\Delta(t)$, one can simply replace $F^T$ with $E^{(i),T}$, the canonical basis covector which has 1 in $i$ and zero otherwise. This selects the specific community $i$ instead of summing over them as in the global case. Then

$$\Delta^{(i)}(t) = \frac{1}{R^{ref}} \frac{E^{(i),T} R^t (R - R^{ref}) \mathbf{x}(0)}{E^{(i),T} R^t \mathbf{x}(0)}. \qquad (S7)$$

Proceeding, again, analogously to what described in the methods for $\Delta(t)$, one can get to

$$\Delta^{(i)}(t) = C^{(i)}(t) \sum_{\alpha} z_{\alpha}^{(i)} \left( 1 - \frac{\Lambda_{\alpha}}{R^{ref}} \right) \left( \frac{\Lambda_{\alpha}}{R^{ref}} \right)^{t}, \qquad (S8)$$

where $z_{\alpha}^{(i)} = -E^{(i),T} w_{\alpha} w_{\alpha}^{*} \mathbf{x}(0)$ encode the weight of each mode. Comparing this with equation (4) and equation (12) of the main paper one notices that the local bias has the
same decay modes as the global ones (the eigenvalues of $\mathbb{R}$), but the weight of each mode is different in each community, and different from the global weight. Also, the global weight of each mode is the sum of the local weights: $z_\alpha = \sum_i z_\alpha^{(i)}$. Importantly, given that the weights $z_\alpha^{(i)}$ may be possibly complex numbers with negative real part and/or nonzero imaginary part (see also Methods of the main paper), this means the weights of modes that are dominant in some communities may cancel out making those modes negligible in determining the speed of convergence of the global estimated reproduction ratio. This means that, in some communities, $\Delta^{(i)}(t)$ may behave very differently from $\Delta(t)$. Specifically, oscillatory modes that are suppressed globally may dominate locally: even when oscillations are not globally visible, some $\Delta^{(i)}(t)$ may oscillate.

# 1.3 Estimation of department-level incident infections

In this section we describe the models used to estimate incident infections from incident hospitalizations in each department of mainland France (excluding Corsica), and compare them to serological estimates. COVID-19-related hospitalizations were collected and made available by the French health authorities and are available at https://www.data.gouv.fr/en/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/.

## 1.3.1 Models

In the main paper we did what we previously did in Ref. [2] and we reconstructed incident infections $I(d)$ from incident hospitalizations $H(d)$ as follows: $I(d) = H(d + 7)/0.032$, where 0.032 is the infection-hospitalization rate (IHR) and 7 is the average time from infection to hospitalization [3]. $d$ indicates the physical time in days.

We also tested an alternative approach featuring age stratification. We used two datasets for hospitalizations: i) department-level, non-age-stratified; ii) region-level, age-stratified. The French health authorities did not release department-level, age-stratified hospitalization data. Note: regions are made up of various departments. In mainland France (excluding Corsica) there are 12 regions and 94 departments. Age-stratification was available at the following age classes: [0-9], [10-19], [20-29], [30-39], [40-49], [50-59], [60-69], [70-79], 80 or older. We used age-stratified IHR for the same age classes from Ref. [4]. Finally, we got the number of residents of department $i$ in age class $a$ from the French national statistical institute (INSEE - https://www.insee.fr). Using the above data, we estimated the incidence of infections in each department of each region as follows. Given that our method applies independently to each region, we assume
in the following that we are focusing on a generic region and i is the index running over the departments of that region, and we omit the index of the region. Analogously, we omit the time index d as we can do the calculation at any time and then shift back by 7 days as we did in the former methodology. So let $H_a^{(reg)}$ be the number of incident hospitalizations in the region in age class a. Let $H_i^{(dep)}$ be the number of incident hospitalizations in department i (non-age-stratified). Let $e_a$ be the IHR of age class a. Let $N_{ia}$ be the number of residents in department i and age class a. With these ingredients, we define $f_{ia}$ as the fraction of individuals in age class a and department i who get infected (incident infections). We factor $f_{ia}$ as the product of a location term $c_i \ge 0$ and an age-class term $q_a \ge 0$, i.e., $f_{ia} = c_i q_a$, with the constraint $\sum_i c_i = 1$. Then, the frequency of individuals in age class a in department i getting hospitalized is $f_{ia}^{(h)} = e_a c_i q_a$. The following consistency equations must hold:

$$ \left\{ \begin{array}{l} \sum_i N_{ia} f_{ia}^{(h)} = H_a^{(reg)} \\ \sum_a N_{ia} f_{ia}^{(h)} = H_i^{(dep)} \end{array} \right. . \qquad (S9) $$

Putting the expression of $f_{ia}^{(h)}$, we get

$$ \left\{ \begin{array}{l} e_a q_a \sum_i N_{ia} c_i = H_a^{(reg)} \\ c_i \sum_a N_{ia} e_a q_a = H_i^{(dep)} \end{array} \right. . \qquad (S10) $$

We can solve equation (S10) numerically for $q_a, c_i$, with the only requirement of parametrizing $c_i$ as a unit simplex given the constraint $\sum_i c_i = 1$. This gives the daily estimated incidence of cases in each department as $c_i \sum_a N_{ia} q_a$.

We show in Fig. S8 alternative versions of Fig. 1b and Fig. 4 of the main paper, obtained using infections estimated with this methodology, proving that choosing either model (age-stratified vs non-age-stratified) gives the same findings.

### 1.3.2 Comparison to seroprevalence data

We found three studies estimating seroprevalence overall in France and in its regions in the period from 2020-03-10 to 2020-06-30:

* SpF: from Ref. [5];

* EpiCov: from Ref. [6];

* Carrat et al.: from Ref. [7].

At the national level, we also compared our estimates (both non-age-stratified and age-stratified) with a model that we previously published [3]. We compared their estimates to the cumulative incidence estimated through our methodologies. See Fig. S5 for
the national estimate and Fig. S6 for the regional estimate. Interestingly, the reconstruc-
tion from non-age-stratified data (the one we use in the main manuscript) outperforms
that from age-stratified data at the national level both according to the available serologi-
cal estimates and our previous model. At the regional level, both perform generally well
except for two regions (Grand-Est and Île-de-France) where the non-age-stratified recon-
struction seems to overestimate infections. Also, in some regions, alternative serological
studies may also give slightly incompatible estimates making it harder to assess which of
our reconstructions schemes works better. In Provences-Alpes-Côte d'Azur, for instance,
the non-age-stratified method is a better match to EpiCov, while the age-stratified one is a
better match to Santé Publique France's. To investigate whether these regional variations
may be associated with different age profiles in the resident population we calculated the
proportion of residents older than 69 years old, and younger than 21 years old, by region:
see Fig. S7. The result is inconclusive: the discrepancy in Île-de-France - the region of
Paris, Île-de-France (IDF) in the maps below - indeed associates with the fact that the pop-
ulation there tends to be younger than the rest of France. But Grand-Est (GES in the map)
has an age profile that is not substantially different from other regions where the non-age-
stratified model does not overestimate infections. Analogoulsy, Provences-Alpes-Côte
d'Azur (PAC in the map) has, for instance, an age profile similar to Occitanie (OCC) but
the performance of the two estimation models is different: in PAC the age-stratified mod-
els seems to better reproduce SpF serological estimates, in OCC the EpiCov ones (the
opposite for the non-age-stratified model), keeping in mind that differences between the
two models are small anyways.

**1.4 The stochastic epidemic model**

We implement a discrete-time stochastic metapopulation model with an SEIR compart-
mental model, with a time step Δd = 1 day. Let S<sub>i</sub>(d), E<sub>i</sub>(d), I<sub>i</sub>(d), R<sub>i</sub>(d) be the number
of residents of department i in each compartment on day d. Being the population of each
department fixed (n<sub>i</sub>), the following continuity equation holds: S<sub>i</sub>(d) + E<sub>i</sub>(d) + I<sub>i</sub>(d) +
R<sub>i</sub>(d) = n<sub>i</sub> ∀d, where n<sub>i</sub> is the population in i. Let β, ε, μ be the rates regulating the tran-
sitions S → E, E → I, I → R. The latency and infectious periods are thus exponentially
distributed with rates ε, μ: ε is the inverse of the average latency period and μ is the inverse
of the average infectious period. One can then compute the probability of individual tran-
sition E → I in the discretization time step Δd as the probability that the latency period
is shorter than Δd: p<sub>ε</sub> = 1 − e<sup>−εΔd</sup>. Analogously, the probability of individual recovery is
p<sub>μ</sub> = 1 − e<sup>−μΔd</sup>. β is the parameter regulating transmissibility. Let p<sub>ij</sub> be the colocation
matrix, then the individual rate of becoming infected is $\beta \sum_j p_{ij} I_j(d)$ – see Ref. [8]. And as before the probability of individual infection is $p_{\beta,i}(d) = 1 - e^{-\beta \sum_j p_{ij} I_j(d) \Delta d}$.

We now describe the d-th time step of the algorithm to compute $E_i(d)$, $I_i(d)$, $R_i(d)$ from $E_i(d-1)$, $I_i(d-1)$, $R_i(d-1)$, using the probabilities just computed:

* sample $E_i^{(S \to E)} \sim \text{Binom}(n = S_i(d-1), p = p_{\beta,i}(d))$;

* sample $I_i^{(E \to I)} \sim \text{Binom}(n = E_i(d-1), p = p_{\epsilon})$;

* sample $R_i^{(I \to R)} \sim \text{Binom}(n = I_i(d-1), p = p_{\mu})$;

* update $E_i(d) = E_i(d-1) + E_i^{(S \to E)} - I_i^{(E \to I)}$;

* update $I_i(d) = I_i(d-1) + I_i^{(E \to I)} - R_i^{(I \to R)}$;

* update $R_i(d) = R_i(d-1) + R_i^{(I \to R)}$.

### 1.4.1 Reference reproduction ratio

The deterministic continuous-time ordinary differential equations of this model are

$$ 
\left\{ 
\begin{aligned}
\dot{E}_i &= \beta S_i \sum_j p_{ij} I_j - \epsilon E_i \\
\dot{I}_i &= \epsilon E_i - \mu I_i \\
\dot{R}_i &= \mu I_i
\end{aligned}
\right.
\qquad (S11) 
$$

In the notation of the next-generation-matrix formalism [9]

$$ 
\left\{ 
\begin{aligned}
\dot{E}_i &= \mathcal{F}_i^E - \mathcal{V}_i^E \\
\dot{I}_i &= \mathcal{F}_i^I - \mathcal{V}_i^I \\
\dot{R}_i &= \mathcal{F}_i^R - \mathcal{V}_i^R
\end{aligned}
\right.
\qquad (S12) 
$$

with $\mathcal{F}_i^E = \beta S_i \sum_j p_{ij} I_j$, $\mathcal{F}_i^I = \mathcal{F}_i^R = 0$, $\mathcal{V}_i^E = \epsilon E_i$, $\mathcal{V}_i^I = \mu I_i - \epsilon E_i$ and $\mathcal{V}_i^R = -\mu I_i$. Then we define the $2N \times 2N$ block matrices $F, V$:

$$ F = \left( \begin{array}{c|c} \frac{\partial \mathcal{F}_i^E}{\partial E_j} & \frac{\partial \mathcal{F}_i^E}{\partial I_i} \\ \hline \frac{\partial \mathcal{F}_i^I}{\partial E_j} & \frac{\partial \mathcal{F}_i^I}{\partial I_i} \end{array} \right)_{E=I=R=0} = \left( \begin{array}{c|c} 0 & \beta n_i p_{ij} \\ \hline 0 & 0 \end{array} \right) \qquad (S13) $$

$$ V = \left( \begin{array}{c|c} \frac{\partial \mathcal{V}_i^E}{\partial E_j} & \frac{\partial \mathcal{V}_i^E}{\partial I_i} \\ \hline \frac{\partial \mathcal{V}_i^I}{\partial E_j} & \frac{\partial \mathcal{V}_i^I}{\partial I_i} \end{array} \right)_{E=I=R=0} = \left( \begin{array}{c|c} \epsilon \delta_{ij} & 0 \\ \hline -\epsilon \delta_{ij} & \mu \delta_{ij} \end{array} \right) \qquad (S14) $$

So, by the next-generation-matrix theory:

$$ NGM_{ij} = (FV^{-1})_{ij} = \frac{\beta}{\mu} n_i p_{ij}. \qquad (S15) $$
Given our reconstruction of R from the colocation matrix (see Methods, Reconstruction of the reproduction operator from data) one can see that $R_{ij} = NGM_{ij}$ if $C = \beta/\mu$, building the final connection between our theoretical formalism and the stochastic simulations. Specifically, the equivalence between the reference reproduction ratio in our calculations and that of the stochastic model: $R^{ref} = \text{spectral rad}(NGM)$.

### 1.4.2 Generation interval and EpiEstim settings

The generation interval is the time between infection and subsequent transmission to another individual. Estimating R through the R-package EpiEstim requires feeding the generation interval distribution associated to the disease[10]. We obtain here the generation interval for our epidemic model.

The compartmental model we use (see Methods and Ref. [11]) has a rate of transition from the E to the I compartment ($\epsilon$) and a recovery rate $\mu$. We also define an effective transmissibility $\mu R^{ref}$ (see equation (15) of main paper). Let $\tau$ be the generation time. Let the probability that a transmission event occurs exactly after $\tau$ has passed since primary infection be $P(t)dt$. Also, let $\tau_E$ be the time one stays in the E compartment: $\tau_E \sim \text{Exp}(\epsilon)$, and $\tau_I$ be the time one stays in the I compartment: $\tau_I \sim \text{Exp}(\mu)$. Conditioning on $\tau_E, \tau_I$, we have

$$P(\tau|\tau_E, \tau_I)dt = \mu R^{ref}dt\theta(\tau - \tau_E)\theta(\tau_I + \tau_E - \tau), \quad (S16)$$

where $\theta$ is the Heaviside's function. So now we marginalize and get the generation time distribution $P(\tau)$:

$$P(\tau) = \int_{0}^{\infty} d\tau_{E}d\tau_{I}P(\tau|\tau_{E}, \tau_{I})P(\tau_{E})P(\tau_{I}) = \frac{\mu\epsilon R^{ref}}{\mu - \epsilon}(e^{-\epsilon\tau} - e^{-\mu\tau}). \quad (S17)$$

Note that a discrete distribution is required by the EpiEstim package. We choose to compute it over 50 bins in the interval [0, 49]. The time window over which to estimate R is set to be a week. In all simulations, we assign the estimate for R returned for the week interval $[t, t + 6]$ to the day $t + 3$. Also, given the smaller precision of early estimates, as reported in the documentation [10], we arbitrarily choose not to plot EpiEstim points associated to the first two weeks from the start of the synthetic epidemic.

For all figures shown in the main manuscript and Supplementary Information, we run multiple runs of the model and compute the reproduction ratio from incident infections produced in each of the runs (according to the framework described in the main paper, with EpiEstim, or both). We show point-wise medians computed over runs and 95%
confidence intervals as statistics of the results. In Fig. S10 and Fig. S11 we show a comparison between so obtained point-wise medians with confidence intervals, and examples of single runs of the stochastic model, in order to test their qualitative match.

## 1.5 Oscillations in European countries

We included 32 European countries: 24 members of the European Union (excluding Cyprus, Ireland and Latvia for lack of data) plus Albania, Bosnia and Herzegovina, Iceland, Montenegro, Norway, Serbia, Sweden, UK (see Data availability). For all of them we built the operator $R$ using colocation and population data at the admin-2 level, similarly to what we did for France. We found at least one real, negative eigenvalue in 11 out of 32 countries, but nowhere did they cause visible oscillations, as the oscillation period was always larger than twice the decay time (Fig. S1). We did not find nonreal eigenvalues. This begs the question whether oscillations are actually observable in real systems. To rigorously determine the conditions for a specific spectrum in a generic nonnegative matrix is not possible, except for specific or low-dimensional cases [12]. We can, however, plausibly associate the presence of an oscillating mode with period $T_α$ to the existence of a cycle of approximate length $T_α$ in the (weighted, directed) network which has $R$ as its adjacency matrix [13, 14]. Slow oscillations (large $T_α$) would then require the presence of long cycles in $R$, which are unlikely to be generated by the recurrent mobility patterns that drive the spatial spread of epidemic outbreaks following pathogen importation [15, 16]. Fast oscillations, and in particular those generated by real, negative eigenvalues, may instead be more common. They would require epidemics that are strongly coupled, i.e., where pairs of communities exist in which infected residents generate, on average, more infections in the other community than in their own, but this is not the case in the countries we examined and for the spatial resolution we considered.
# 2 Supplementary Figures

## 2.1 Decay time and oscillation period of the modes of

$\Delta = R^{estim} - R^{ref}$ for selected European countries

<figure><img src="image_7.png" /><figcaption>**Fig.S1** Decay time and oscillation period of the modes of $\Delta = R^{estim} - R^{ref}$ for selected European countries. The color shows the value, for each country, $\max_{\alpha} \tau_{\alpha}/T_{\alpha}$, where $\tau_{\alpha}, T_{\alpha}$ are defined in the main paper and are, respectively, the decay time and the oscillation period of the $\alpha$th mode. 32 countries are included: 24 members of the European Union (excluding Cyprus, Ireland and Latvia for lack of data) plus Albania, Bosnia and Herzegovina, Iceland, Montenegro, Norway, Serbia, Sweden, UK (see Data availability). Countries with real positive eigenvalues only are colored in white ($T_{\alpha} = \infty$), countries not included are in gray.</figcaption></figure>
## 2.2 Relative effect of the diagonal correction to co-location probabilities in French departments

<figure><img src="image_8.png" /><figcaption>**Fig.S2** Relative effect of the diagonal correction to co-location probabilities in French departments.</figcaption></figure>

a shows the relative effect of the diagonal corrections in the 94 departments of mainland France in week 5 of 2021, with respect to population densities. The correction is substantial in departments with high population density, consistently with the nature of the bias we want to correct. See Methods, Reconstruction of the reproduction operator from data. b map of the relative diagonal corrections.
## 2.3 Daily number of new COVID-19-related infections per 100,000 residents

<figure><img src="image_9.png" /><figcaption>**Fig.S3** Daily number of new COVID-19-related infections per 100,000 residents. Daily incident infections from week 9 of 2020 to week 3 of 2022 (Feb 23, 2020 to Jan 23, 2022), according to the estimation method we use in the main paper. Each curve is that of one of the 94 departments (administrative level 2) of mainland France excluding the region of Corsica, with a 14 days window rolling average. The area in gray goes from week 2 to week 9 of 2021 (Jan 11 to Mar 7), corresponding to the period on which we focus in Fig. 1 of the main paper. The small Dept. #90 (Territoire de Belfort), containing 0.2% of the French population, is excluded for visibility as it had an anomalously high estimated incidence peaking at around 600 due to a miscounting of hospitalizations coming from neighboring departments. This fact was known to national health authorities.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig.S4** Daily number of new COVID-19-related infections per 100,000 residents in the 12 regions of mainland France. We show new COVID-19-related daily infections estimated according to the method used in the main paper, separately for the departments in each of the twelve regions of mainland France, to facilitate readability. Names of regions (three letters codes) are reported in plot titles, while INSEE codes of departments are reported in the legend. Insets zoom in on the period from week 2 of 2021 to week 9 of 2021, the focus in Fig. 1 of the main paper. A 14 days window rolling average is used. The small Dept. #90 (Territoire de Belfort) in region Bourgogne-Franche-Comté (BFC), containing 0.2% of the French population, had an anomalously high estimated incidence peaking at around 600 due to a miscounting of hospitalizations coming from neighboring departments. This fact was known to national health authorities.</figcaption></figure>
2.4 National attack rate: estimates vs serological studies

<figure><img src="image_11.png" /><figcaption>**Fig.S5** National attack rate: estimates vs serological studies. This figure presents the cumulative attack rate across whole France, as per non-age-stratified estimation, age-stratified estimation, Pullano et al. estimation ([3]) and serological studies from SpF ([5]), EpiCov ([6]), and Carrat et al. ([7]). The estimates are plotted over time for the period from April to June 2020.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig.S6** Regional attack rate: estimates vs serological studies. This figure presents the cumulative attack rate across different French regions, as per the non-age-stratified estimation, age-stratified estimation, and serological studies from SpF ([5]), EpiCov ([6]), and Carrat et al. ([7]). The estimates are plotted over time for the period from April to June 2020.</figcaption></figure>
2.5 Age demographic distribution in French regions

<figure><img src="image_13.png" /><figcaption>**Fig.S7** | Age distribution in French regions. The top map shows the percentage of the resident population aged 70 and above, while the bottom map details the percentage of the population aged 0 to 20.</figcaption></figure>
2.6 Fig. 1b and Fig. 4 of the main paper, using the age-stratified model for the estimate of infections

<figure><img src="image_14.png" /><figcaption>**Fig.S8** Fig. 1b and Fig. 4 of the main paper, using the age-stratified model for the estimate of infections. We show an alternative version of Fig.1b and Fig.4 of the main paper, respectively in **a** and **b**. Here we consider new COVID-19 infections reconstructed as per the age-stratified model described in Estimation of department-level incident infections. Results are in agreement with the discussion related to Fig. 1b and Fig. 4 in the main paper. Also, we stress that the validity of our theoretical results is independent of the methodology employed for cases reconstruction.</figcaption></figure>
## 2.7 Comparison of the reference, estimated and corrected reproduction ratios with six different initial conditions for the stochastic model

<figure><img src="image_15.png" /><figcaption>**Fig.S9** | Comparison of the reference, estimated and corrected reproduction ratios with six different initial conditions. This figure uses the same stochastic epidemic model as in Fig. 2 and Fig. 3 of the main paper, with six different initial conditions. 500 initial cases are seeded in Alpes-Maritimes (a), Rhône (b), Bouches-du-Rhône (c), Hérault (d), Haute-Garonne (e) and Loire-Atlantique (f). The largest city in each of these departments is written in parenthesis in legend, together with the departments' population, the initial (day=0) relative difference between the estimated and reference reproduction ratios, and the time in days required for the absolute relative error to be smaller than 0.05. Both the initial measurement error and the speed of convergence depend on the initial conditions, as discussed in the main paper. We show medians and 95% confidence intervals for reproduction ratios over 1,000 runs.</figcaption></figure>
## 2.8 Comparison between point-wise medians and single runs of the stochastic metapopulation model

<figure><img src="image_16.png" /><figcaption>**Fig.S10|**Comparison between point-wise medians and single runs of the stochastic model: estimates of the reproduction number. This figure uses the same stochastic epidemic model as in Fig. 2 and Fig. 3 of the main paper, and one of their two initial conditions, i.e. 100 initial cases are seeded among departments proportionally to their population. We show the standard reproduction ratio measured from surveillance data – both in the Galton-Watson formalism (red, **a**) and with *EpiEstim* (purple, **c**) on incidence data – and the corresponding *v*\*-corrected estimates: Galton-Watson formalism (blue, **b**) and with *EpiEstim* (light blue, **d**) on modified incidence data. We show point-wise medians and confidence intervals over 10,000 runs, and we compare them with the results of single runs of the model, in order to test their qualitative match.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig.S11** | Comparison between point-wise medians and single runs of the stochastic model: locally estimated reproduction numbers. This figure uses the same stochastic epidemic model as in Fig. 2 and Fig. 3 of the main paper, and one of their two initial conditions, *i.e.* 100 initial cases are seeded among departments proportionally to their population. We show the locally estimated reproduction numbers in 6 among the 94 departments of mainland France, which we report in the legend of single plots (**a**, **b**, **c**, **d**, **e**, **f**). We show point-wise medians and confidence intervals over 10,000 runs, computed excluding points associated to locally extinct runs. Time points in which local incidence is zero are omitted.</figcaption></figure>
## 2.9 Variation of the spatial distribution of infections between weeks 6 and 11 of 2021, and the equilibrium distribution

<figure><img src="image_18.png" /><figcaption>**Fig.S12** Variation of the spatial distribution of infections between weeks 6 and 11 of 2021, against the transmission potential of communities. We show on the x axis the transmission potential of each department (rescaled so that the highest equals one) and on the y axis the difference between the fraction of national infections in that department in Week 6 and 11 of 2021. Departments with high transmission rates tended to increase their fraction of total infections, departments with low transmission rates tended to decrease it. The red line is a simple linear fit to help underline the trend. This is the overall trend, yet some departments still had their share of infections decrease even if they had a fairly high $\rho$ (the departments of the cities of Nice, Alpes-Maritimes, and Marseille, Bouches-du-Rhône, in the South of France). This is consistent with the values of these departments in the equilibrium spatial distribution of infections shown in Supplementary Fig. S13.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig.S13** Equilibrium distribution of infections computed from the average reproduction operator between week 51 of 2020 and week 12 of 2021. We show the average equilibrium distribution of infections in the time window of Fig. 1b and Fig. 4 of the main paper. Consistently with the known centrality of Paris in the French mobility network, the equilibrium distribution has high values in Île-de-France. The departments of Paris, Alpes-Maritimes and Bouches-du-Rhône are highlighted with arrows, following the discussion of Supplementary Fig. S12.</figcaption></figure>
Supplementary References

1. Horn, R. A. & Johnson, C. R. *Matrix Analysis* ISBN: 0-521-38632-2 (Cambridge University Press, 1990).

2. Mazzoli, M., Valdano, E. & Colizza, V. Projecting the COVID-19 epidemic risk in France for the summer 2021. *Journal of Travel Medicine* **28**. ISSN: 1708-8305 (Oct. 2021).

3. Pullano, G. et al. Underdetection of COVID-19 cases in France threatens epidemic control. *Nature* **590** (Feb. 2021).

4. Le Vu, S. et al. Prevalence of SARS-CoV-2 antibodies in France: results from nationwide serological surveillance. en. *Nature Communications* **12**. Number: 1 Publisher: Nature Publishing Group, 3025. ISSN: 2041-1723 (May 2021).

5. Le Vu, S. et al. Prevalence of SARS-CoV-2 antibodies in France: results from nationwide serological surveillance. *Nature communications* **12**, 3025 (2021).

6. Ministère des Solidarités et de la Santé. En mai 2020, 4,5% de la population en France métropolitaine a développé des anticorps contre le SARS-CoV-2. Premiers résultats de l'enquête nationale EpiCov Accessed: 2023-11-29. 2020.

7. Carrat, F. et al. Antibody status and cumulative incidence of SARS-CoV-2 infection among adults in three regions of France following the first lockdown and associated risk factors: a multicohort study. *International journal of epidemiology* **50**, 1458–1472 (2021).

8. Iyer, S. et al. Large-scale measurement of aggregate human colocation patterns for epidemiological modeling. *Epidemics* **42**, 100663. ISSN: 1755-4365 (2023).

9. Diekmann, O., Heesterbeek, J. & Roberts, M. The construction of next-generation matrices for compartmental epidemic models. *Journal of the Royal Society, Interface / the Royal Society* **7**, 873–85 (Nov. 2009).

10. Cori, A., Ferguson, N. M., Fraser, C. & Cauchemez, S. A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics. *American Journal of Epidemiology* **178**, 1505–1512. ISSN: 0002-9262 (Nov. 2013).

11. Faucher, B. et al. Agent-based modelling of reactive vaccination of workplaces and schools against COVID-19. *Nature Communications* **13**, 1414 (Mar. 2022).

12. Egleston, P. D., Lenker, T. D. & Narayan, S. K. The nonnegative inverse eigenvalue problem. en. *Linear Algebra and its Applications. Special Issue on the Tenth ILAS Conference (Auburn, 2002)* **379**, 475–490. ISSN: 0024-3795 (Mar. 2004).
13. Kellogg, R. B. & Stephens, A. B. Complex eigenvalues of a non-negative matrix with a specified graph. en. *Linear Algebra and its Applications* **20**, 179–187. ISSN: 0024-3795 (Jan. 1978).

14. Torre-Mayo, J., Abril-Raymundo, M. R., Alarcia-Estévez, E., Marijuán, C. & Pisonero, M. The nonnegative inverse eigenvalue problem from the coefficients of the characteristic polynomial. EBL digraphs. en. *Linear Algebra and its Applications* **426**, 729–773. ISSN: 0024-3795 (Oct. 2007).

15. Schneider, C. M., Belik, V., Couronné, T., Smoreda, Z. & González, M. C. Unraveling daily human mobility motifs. *Journal of the Royal Society Interface* **10**, ISSN: 17425662 (2013).

16. Gómez-Gardeñes, J., Soriano-Paños, D. & Arenas, A. Critical regimes driven by recurrent mobility patterns of reaction-diffusion processes in networks. en. *Nature Physics* **14**, Number: 4 Publisher: Nature Publishing Group, 391–395. ISSN: 1745-2481 (Apr. 2018).