# Epidemiological indices with multiple circulating pathogen strains

Cristiano Trevisin, Lorenzo Mari, Marino Gatto, Vittoria Colizza, Andrea Rinaldo

► To cite this version:

Cristiano Trevisin, Lorenzo Mari, Marino Gatto, Vittoria Colizza, Andrea Rinaldo. Epidemiological indices with multiple circulating pathogen strains. Infectious Disease Modelling, 2025, 10 (3), pp.802-812. 10.1016/j.idm.2025.03.006. hal-05058136

HAL Id: hal-05058136

https://hal.science/hal-05058136v1

Submitted on 6 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Epidemiological indices with multiple circulating pathogen strains

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

Cristiano Trevisin<sup>a, b, *</sup>, Lorenzo Mari<sup>c</sup>, Marino Gatto<sup>c</sup>, Vittoria Colizza<sup>b</sup>, Andrea Rinaldo<sup>a, d</sup>

## ARTICLE INFO

**Article history:**

Received 15 January 2025

Received in revised form 12 March 2025

Accepted 12 March 2025

Available online 17 March 2025

Handling Editor: Dr Daihai He

**Keywords:**

COVID-19

Epidemic control

Variants

Reproduction numbers

Epidemicity indices

Infectious diseases

## ABSTRACT

Epidemiological indicators (e.g. reproduction numbers and epidemicity indices) describe long- and short-term behaviour of ongoing epidemics. Their evolving values provide context for designing control measures because maintaining both indices below suitable thresholds warrants waning infection numbers. However, current models for the computation of epidemiological metrics do not consider the stratification of the pathogen into variants endowed with different infectivity and epidemiological severity. This is the case, in particular, with SARS-CoV-2 infections. Failing to account for the variety of epidemiological features of emerging variants prevents epidemiological indices from spotting the possible onset of uncontrolled growth of specific variants, thus significantly limiting the prognostic value of the indicators. Here, we expand an existing framework for the computation of spatially explicit reproduction numbers and epidemicity indices to account for arising variants. By analysing the data of the COVID-19 pandemic in Italy, we show that embedding additional layers of complexity in the mathematical descriptions of unfolding epidemics reveals new angles. In particular, we find epidemiological metrics significantly exceeding their thresholds at the emergence of new variants. Such values foresee a recrudescence in new infections that only becomes evident after emerging new variants have effectively replaced the previous active strains. The demography of the variant composition flags the presence of specific strains growing more rapidly than the total number of infections generated by all variants combined. Variant-aware epidemiological indicators thus allow to engineer better control measures tailored to the shifting patterns of severity and evolving features of infectious disease epidemics.

© 2025 The Authors. Publishing services by Elsevier B.V. on behalf of KeAi Communications Co. Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

<footnote><sup>a</sup> Laboratory of Ecohydrology, Ecole Polytechnique Fédérale de Lausanne, Station 2, Lausanne, 1015, Switzerland</footnote>

<footnote><sup>b</sup> Sorbonne Université, INSERM, Pierre Louis Institute of Epidemiology and Public Health, 27 rue Chaligny, Paris, 75012, France</footnote>

<footnote><sup>c</sup> Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Via Ponzio 34/5, Milano, 20133, Italy</footnote>

<footnote><sup>d</sup> Dipartimento di Ingegneria Civile, Edile e Ambientale (ICEA), Università di Padova, Via Marzolo 9, Padova, 35131, Italy</footnote>

<footnote>* Corresponding author. Sorbonne Université, INSERM, Pierre Louis Institute of Epidemiology and Public Health, 27 rue Chaligny, Paris, 75012, France. E-mail address: cristiano.trevisin@inserm.fr (C. Trevisin).</footnote>

<footnote>Peer review under the responsibility of KeAi Communications Co., Ltd.</footnote>
# 1. Introduction

After its emergence in late 2019, the virus responsible for the COVID-19 pandemic, SARS-CoV-2, has undergone many mutations. A major effort has been devoted to stemming the spread of the virus after the pandemic outbreak—via both vaccines and non-pharmaceutical interventions (Bertuzzo et al., 2020; Colizza et al., 2021; Lemaitre et al., 2022; Leung et al., 2023)—generated strong evolutionary pressure on the virus to evade immunity conferred by prior vaccination and infection (Carabelli et al., 2023; Grabowski et al., 2022; Koelle et al., 2022). In particular, newly emerging strains achieved higher reproduction numbers, shorter doubling times (Goldberg et al., 2023; Koelle et al., 2022) and reduced their mean generation time (Kim et al., 2022; Zhang et al., 2021). This pattern has not only affected significantly the course of SARS-CoV-2, but also that of other infectious diseases such as dengue fever and influenza (Guzman & Harris, 2015; Lin et al., 1999). Several studies have included variant stratification in their compartmental models to understand the effect of control measures (Saviano et al., 2023; Cuevas-Maraver et al., 2024) or to predict the course of an epidemic (Baccega et al., 2024) in the presence of co-circulating strains. Another study (Sabbatini et al., 2024) accounted for the effect of the variant prevalence to estimate the exposure and the hospitalization rates.

The fate of an epidemic is recapitulated by the computation of effective prognostic indicators. Foremost among them is the reproduction number, $R$, which measures the number of secondary infections produced by each primary infector within a given community. Values of $R$ above the unity flag an accelerating spread of infections, whereas subunit values characterize an asymptotic waning epidemic because no endemic equilibrium may be established. In practice, reproduction numbers command stricter or milder imposed restrictions via the enforcement of control measures aimed at achieving subunit values of $R$. Depending on the chosen epidemiological modelling framework, $R$ may be computed via next-generation matrices, whenever a SIR-like compartmental model is adopted (Diekmann et al., 2010; Van Den Driessche, 2017), or via renewal equations (Cori et al., 2013) that allow one to compute reproduction numbers directly from data.

Another key indicator is the epidemiicity index (Mari et al., 2018; Neubert & Caswell, 1997; Trevisin et al., 2022, 2024) pinpointing the largest amplification attained in the short-term by any impulsive perturbation to a disease-free equilibrium (Mari et al., 2017, 2018, 2021, 2024). Notably, the asymptotic stability condition $R < 1$ does not prevent local transient outbursts of infections (Mari et al., 2024). It has been formally demonstrated that excluding the possibility of short-term flare-ups requires $R$ substantially lower than its endemicity threshold ($R \ge 1$) (Mari et al., 2024; Trevisin et al., 2022, 2024). This questions the prognostic power of reproduction numbers alone to sort out the most effective control measures during in times of spreading infections.

Within a metapopulation, infections may be transmitted among spatially connected communities via e.g. contacts driven by human mobility. The COVID-19 pandemic highlighted the relevance of spatial connectivity in both our understanding of the spread of SARS-CoV-2 towards proper identification of control measures. Spatial formulations of the reproduction number (Birello et al., 2024; Trevisin et al., 2023) and the epidemiicity index (Trevisin et al., 2024) were derived. Results highlighted how incorporating a spatial stratification improves the estimates of $R$, in particular because they account for export/import of infections within connected communities, allowing a more precise identification of appropriate control measures.

Another stratification, pertaining to the emergence of different variants of the virus of interest, is still missing to date. We posit that embedding strain diversification in a framework for the computation of epidemiological indices and looking at anomalies in their temporal patterns sheds light on the behaviour of new variants. To test this hypothesis, we propose to expand the existing spatially explicit model for the computation of epidemiological metrics (Mari et al., 2024; Trevisin et al., 2024) to provide a comprehensive framework that seamlessly derives, in real-time, space- and variant-specific values of both the reproduction number and the epidemiicity index. This is claimed to allow us to single out the fastest-growing strain within the total number of new infections, thus anticipating possible recrudescence in new infections. Our approach differs from previous attempts at incorporating variant stratification into compartmental models as we do not calibrate a compartmental model but rather infer time-series of epidemiological indices that may serve as a basis for control measures. While our framework requires one to know the variants' generation times in advance, it does not require specific assumptions for direct transmission, and may be computationally lighter than compartmental approaches.

# 2. Methods

We propose a framework based on a Leslie projection matrix (Mari et al., 2024; Trevisin et al., 2024), here expanded to account for variant emergence. A key assumption is that once a variant has emerged, it is assumed to act independently of the other circulating strains. Moreover, any interaction among variants (notably competition) is subsumed by their respective $R$ and generation time distributions ($\beta(\tau)$ i.e., the time elapsed between the infector's exposure and that of their infectee(s)).

Let us consider a metapopulation composed of $N$ interconnected human communities (or nodes), where several variant genotypes ($v = 1, ..., V$) of a generic infectious disease are circulating. The compartment infected by strain $v$ in an arbitrary community $j$ ($j = 1, ..., N$) is indicated as $I_j^v(t, \tau)$, a variable of both time $t$ and age of infection $\tau$. Also, let $C_{ij}(t)$ be the proportion of residents of the community $j$ who commute daily to $l$ for their daily activities (Trevisin et al., 2023). The following differential problem holds (Mari et al., 2024; Trevisin et al., 2023):
$$ \frac{\partial I_j^v(t, \tau)}{\partial t} + \frac{\partial I_j^v(t, \tau)}{\partial \tau} = -\gamma^v(\tau)I_j^v(t, \tau) \quad (1a) $$

$$ I_j^v(t, 0) = \sum_{k=1}^{N} z_{jk}(t) \mathcal{R}_k^v(t) \int_0^\infty \phi^v(\tau) I_k^v(t, \tau) d\tau \quad \forall t \ (\text{B.C.}), \quad (1b) $$

where $\gamma^v(\tau)$ represents the instantaneous rate of exit from the infected compartment (because of recovery, death, quarantine, or isolation), and $\phi^v(\tau)$ is the rate of secondary transmission per single infectious case. Regardless of the variant being considered, the elements $z_{jk}$ are defined as

$$ z_{jk}(t) = \sum_{l=1}^{N} \frac{C_{lj}(t)n_j}{\sum_{m=1}^{N} C_{lm}(t)n_m} C_{lk}(t), \quad (2) $$

denoting each node's exposure to infectious individuals residing in a different one via several spatial routes of infection. Both matrices $\mathbf{C} = [C_{ij}]$ and $\mathbf{Z} = [z_{ij}]$ are column-stochastic. The rate of secondary transmission $\phi^v(\tau)$ relates to the distribution of generation times, $\beta^v(\tau)$, via the following relation:

$$ \beta^v(\tau) = p^v(\tau)\phi^v(\tau), \quad (3) $$

where $p^v(\tau) = \exp(-\int_0^\tau \gamma^v(\xi)d\xi)$ represents the fraction of infected individuals that are still infectious $\tau$ days after infection.

To comply with a daily or weekly collection of epidemiological data, one may discretize the above process with respect to both the time and the age of infection, which leads to the following system of finite-difference equations:

$$ I_{j,1}^v(t) = \sigma_0^v \sum_{k=1}^{N} z_{jk}(t) \mathcal{R}_k^v(t) \sum_{i=1}^{q} f_i^v I_{k,i}^v(t-1) \quad (4a) $$

$$ I_{j,i}^v(t) = \sigma_{i-1}^v I_{j,i-1}^v(t-1) \quad \text{for } i \ge 2, \quad (4b) $$

where $f_i^v = \int_{i-1}^i \phi^v(\tau) d\tau$ is a discretization in time of function $\phi^v(\tau)$, and $\sigma_{i-1}^v = \frac{p^v(i)}{p^v(i-1)}$ represents the proportion of infectious individuals with age of infection $i-1$ still infectious on the next day. Inference of terms $\mathcal{R}_k^v(t)$ is possible via sequential Monte Carlo algorithms (Trevisin et al., 2023).

## 2.1. The multi-strain leslie projection matrix

Let us rewrite process (4) in matrix form and define a column vector $\mathbf{I}(t)$ collecting all the infectious compartments at time $t$, that is:

$$ \mathbf{I}(t) = \bigcup_{v=1}^{V} \mathbf{I}^v(t), \quad (5) $$

with:

$$ \mathbf{I}^v(t) = \{I_{1,1}^v, \dots, I_{N,1}^v, \dots, I_{1,q}^v, \dots, I_{N,q}^v\}^T, \quad (6) $$

where $\mathbf{T}$ indicates matrix transposition.

The projection matrix associated with model (4) is defined such that $\mathbf{I}^v(t+1) = \mathbf{L}^v(t)\mathbf{I}^v(t)$. The non-zero block $\mathbf{L}^v(t)$ for a given variant strain $v$ has the structure of a so-called Leslie projection matrix and reads

$$ \mathbf{L}^v(t) = \begin{pmatrix} \mathbf{T}_1^v(t) & \mathbf{T}_2^v(t) & \cdots & \mathbf{T}_{q-1}^v(t) & \mathbf{T}_q^v(t) \\ \sigma_1^v \mathbf{1}_N & \mathbf{0}_N & \cdots & \mathbf{0}_N & \mathbf{0}_N \\ \mathbf{0}_N & \sigma_2^v \mathbf{1}_N & \cdots & \mathbf{0}_N & \mathbf{0}_N \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0}_N & \mathbf{0}_N & \cdots & \sigma_{q-1}^v \mathbf{1}_N & \mathbf{0}_N \end{pmatrix}, \quad (7) $$

where $\mathbf{1}_N$ and $\mathbf{0}_N$ are the identity and null matrix of order $N$, respectively, and the sub-matrices $\mathbf{T}_i^v(t)$ are defined as:
$$ \mathbf{T}_i^v(t) = \sigma_0^v \mathbf{f}_i^v \begin{pmatrix} \mathcal{R}_1^v(t) z_{11}(t) & \mathcal{R}_2^v(t) z_{12}(t) & \dots & \mathcal{R}_N^v(t) z_{1N}(t) \\ \mathcal{R}_1^v(t) z_{21}(t) & \mathcal{R}_2^v(t) z_{22}(t) & \dots & \mathcal{R}_N^v(t) z_{2N}(t) \\ \vdots & \vdots & \ddots & \vdots \\ \mathcal{R}_1^v(t) z_{N1}(t) & \mathcal{R}_2^v(t) z_{N2}(t) & \dots & \mathcal{R}_N^v(t) z_{NN}(t) \end{pmatrix}. \quad (8) $$

Therefore, the first N rows of $\mathbf{L}^v(t)$ represent the appearance of new infections. In contrast, each of the $\sigma_i^v \mathbf{I}_N$ blocks corresponds to the transition of infectious individuals with age of infection i to their i + 1-th day of infection.

For multiple circulating variants, the global projection matrix $\mathbf{M}(t)$ is block-diagonal and reads:

$$ \mathbf{M}(t) = \begin{pmatrix} \mathbf{L}^1(t) & \mathbf{0}_{Nq} & \cdots & \mathbf{0}_{Nq} \\ \mathbf{0}_{Nq} & \mathbf{L}^2(t) & \cdots & \mathbf{0}_{Nq} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0}_{Nq} & \mathbf{0}_{Nq} & \cdots & \mathbf{L}^v(t) \end{pmatrix}, \quad (9) $$

with $\mathbf{0}_{Nq}$ being a null matrix of order $Nq$.

## 2.2. Long-term analysis

To assess the long-term behaviour of the system and the feasibility of a disease-free equilibrium, one first has to extract the next-generation matrix (van den Driessche & Watmough, 2002) and its spectral radius from the relevant Leslie projection matrix. The latter is defined as global effective reproduction number (Mari et al., 2021; Trevisin et al., 2024). To yield the next-generation matrix, each diagonal block $\mathbf{L}^v(t)$ of the general projection matrix $\mathbf{M}(t)$ must be split into two matrices, a transmission matrix $\mathbf{T}^v(t)$ and a transition matrix $\mathbf{\Sigma}^v(t)$, such that $\mathbf{T}^v(t) + \mathbf{\Sigma}^v(t) = \mathbf{L}^v(t)$. The variant-specific next-generation matrix (Allen & van den Driessche, 2008) then reads

$$ \mathbf{K}^v(t) = \mathbf{T}^v(t) (\mathbf{1}_{N \times q} - \mathbf{\Sigma}^v(t))^{-1}, \quad (10) $$

where $\mathbf{1}_{Nq}$ is the identity matrix of order $Nq$. It is quite straightforward to prove that the following relation holds:

$$ \mathbf{K}(t) = \begin{pmatrix} \mathbf{K}^1(t) & \mathbf{0}_{Nq} & \cdots & \mathbf{0}_{Nq} \\ \mathbf{0}_{Nq} & \mathbf{K}^2(t) & \cdots & \mathbf{0}_{Nq} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0}_{Nq} & \mathbf{0}_{Nq} & \cdots & \mathbf{K}^v(t) \end{pmatrix}. \quad (11) $$

For each variant, its national $\mathcal{R}_v^G$ (i.e., the indicator that predicts the long-term fate of that specific variant in the considered metapopulation) is defined as the spectral radius of the variant-specific next-generation matrix, $\mathbf{K}^v(t)$:

$$ \mathcal{R}_v^G(t) = \rho(\mathbf{K}^v(t)), \quad (12) $$

with $\rho(\cdot)$ indicating the spectral radius operator. In particular, the condition $\mathcal{R}_v^G(t) < 1$, if maintained for all t's, guarantees that strain v cannot establish in the community, thus preventing endemic transmission. Because of the property (11), one can easily conclude that the global $\mathcal{R}^G$, which encompasses all variants and predicts the fate of the disease as a whole, is the largest of the variant-specific global $\mathcal{R}$:

$$ \mathcal{R}^G(t) = \max_v \mathcal{R}_v^G(t). \quad (13) $$

If $\mathcal{R}^G(t) < 1$ for all t's, no variant can ever establish in the community.

## 2.3. Short-term reactivity analysis

One may also analytically assess whether a short-term flare-up is possible under given epidemiological conditions. Similarly to the long-term analysis, the block-diagonal nature of the Leslie projection matrix brings a substantial simplification of the relevant algebra. As such, one may define the epidemiicity index for each variant subset v as the maximum one-step growth of the compartments infected with variants v, that is:
$$ \mathcal{E}_v(t) = \max_{\|\mathbf{Y}^v(t+1)\|_p \neq 0} \frac{\|\mathbf{Y}^v(t+1)\|_p}{\|\mathbf{Y}^v(t)\|_p} = \max_{\|\mathbf{I}^v(t)\|_p \neq 0} \frac{\|\mathbf{W}\mathbf{I}^v(t)\|_p}{\|\mathbf{W}\mathbf{I}^v(t)\|_p}, \quad (14) $$

where $\|x\| = \sum_i |x_i|$ is the $\ell^1$-norm of vector $x$, $\mathbf{I}^v(t)$ is the set of infectious classes (in space and age of infection) for strain $v$ at time $t$, $\mathbf{Y}^v$ is a customary observed variable based on $\mathbf{I}^v(t)$, and $\mathbf{W}$ is a suitable transformation matrix. The condition $\mathcal{E}_v(t) < 1$, if maintained for all $t$'s, prevents the abundance of individuals infected with strain $v$ from ever increasing in the short term, even if infected people enter the community. Other algebraic norms (see Trevisin et al. (2024)) are possible, but here we only focus on the $\ell^1$-norm because of its straightforward interpretation. The block-diagonal property of the general projection matrix $\mathbf{M}(t)$, allows one to define the global epidemiicity index as:

$$ \mathcal{E}^G(t) = \max_v \mathcal{E}_v(t). \quad (15) $$

If $\mathcal{E}^G(t) < 1$ for all $t$'s, no variant can ever undergo a temporary prevalence increase in the population.

In the following, we discuss the short-term epidemiicity of the system through the observed variable $\mathbf{Y}$, defined as the vector of length $NV$ that contains as elements the absolute abundances of infectious individuals per each variant strain and spatial node. This vector reads

$$ \mathbf{Y} = \bigcup_{v=1}^{N} \mathbf{Y}^v, \quad \text{with } \mathbf{Y}^v = \{Y_j^v, \text{ for } j=1, \dots, N\}, \quad (16) $$

where

$$ Y_j^v = \sum_{i=1}^{q} I_{j,i}^v. \quad (17) $$

The corresponding transformation matrix $\mathbf{W}$ such that $\mathbf{Y}(t) = \mathbf{W}\mathbf{I}(t)$ takes the form

$$ \mathbf{W} = \underbrace{(\mathbf{1}_q \quad \mathbf{1}_q \quad \cdots \quad \mathbf{1}_q)}_{N \times V \text{ times}}, \quad (18) $$

with $\mathbf{1}_q$ being a diagonal matrix of order $q$. Computing the epidemiicity index applying the $\ell^1$-norm to this observed variable highlights the variant-node pairs that could achieve the largest short-term growth in the total number of infected individuals. Other state transformations could be applied as well (see Trevisin et al., 2024) but here we limit our attention to this one out of simplicity.

Given the column-stochasticity of matrix $\mathbf{Z}(t) \forall t$, it is possible to derive an explicit formula for the computation of the epidemiicity index for variant $v$, namely

$$ \mathcal{E}^v(t) = \sigma_0^v \max_j \mathcal{R}_j^v(t) \max_i f_i^v + \sigma^v, \quad (19) $$

with the additional assumptions that $\sigma_i^v = \sigma^v \forall i \in \{1, \dots, q-1\}$, indicating a constant rate of exit from the infectious compartments, and $f_q^v = 0$, which underlines negligible infectiousness at the age of infection $\tau \ge q$.

Inverting this formula allows one to compute the threshold reproduction number $\mathcal{R}_\star^v < 1$ that should be achieved to guarantee that a specific variant strain is non-reactive (Mari et al., 2024):

$$ \mathcal{R}_\star^v = \frac{1 - \sigma^v}{\sigma_0^v \max_i f_i^v}. \quad (20) $$

Thanks to the block-diagonal property of the Leslie projection matrix, one may easily compute the global epidemiicity index as the largest variant-specific one:

$$ \mathcal{E}^G(t) = \max_v \mathcal{E}^v(t) = \max_{v,k} \left[ \sigma_0^v \mathcal{R}_k^v(t) \max_i f_i^v + \sigma^v \right]. \quad (21) $$

## 2.4. Application to COVID-19

While this study can be implemented for any directly transmitted infectious disease, we tailor it to the COVID-19 pandemic in Italy. The distributions of generation time for the considered SARS-CoV-2 variants are assumed to be gamma-shaped, with
details shown in Table 1 and shown in Fig. 1. The depletion of the infectious compartment, regardless of the considered variant, follows a negative exponential distribution, namely $p(\tau) = \exp(-0.068\tau)$, which yields $\sigma_j = \exp(-0.068) \forall j \ge 1$, corresponding to a mean residence time within the infectious compartment of approximately 11.7 days (Gatto et al., 2020). Whenever estimates of the generation times' distribution are not available, the serial interval is used instead (Chen et al., 2022; Svensson, 2007; Torneri et al., 2021).

Daily reported infections in the 20 Italian regions are gathered from the Italian Department of Public Protection (github.com/pcm-dpc/COVID-19). As for human mobility is considered, the connectivity matrix C as in (Trevisin et al., 2024) through mobility data collected from the Italian Institute of Statistics (istat.it/it/archivio/139381) and updated via the “Workplace mobility” time series of the Google Community Mobility Reports (google.com/covid19/mobility/). The simulations shown here extend from February 23, 2020, through October 15, 2022. The data on the prevalence of variants in Italy are obtained from the data portal of the European Centre for Disease Prevention and Control (GISAID, see Khare et al. (2021); Elbe and Buckland-Merrett (2017); Shu and McCauley (2017), and TESSy frameworks) and is shown in Fig. 1. The following variants are considered for this study: the wild type, Alpha (B.1.1.117), Gamma (P.1), Delta (B.1.157.2), and Omicron subtypes (BA.1, BA.2, BA.4, BA.5+BQ.1). To contain the computational burden, all other variants subtypes are not considered, as their prevalence among sequenced samples was low at all times. In addition, isolated spurs in the reporting of cases associated with a given variant occurring more than a month apart from that variant's main transmission period are removed to reduce the computational burden.

## 2.5. Numerical experiments

In the following, we compare a variant-specific scenario and a non-stratified one. Concerning the non-stratified scenario, the vector $Y_1$ has length N, and its transformation matrix has size (q, N). Model (4) loses its dependency on the variant and takes the form (Trevisin et al., 2024)

<caption>**Table 1**</caption>

Metrics of the generation time distribution of the considered variants in the multi-variant setup. The mean of the generation time $\mu_{\delta}$ is given. The standard deviation $\sigma_{\delta}$ is obtained from the confidence intervals. $R_{\star}^{v}$ computed using Equation (20) and assuming a daily discretization of the generation time distribution are also shown.

<table><thead><tr><th>Common name</th><th>Lineage</th><th>μ<sub>δ</sub> (95 % CI) [days]</th><th>σ<sub>δ</sub> [days]</th><th>R<sub>★</sub><sup>v</sup></th><th>Ref.</th></tr></thead><tbody><tr><td>Wild type</td><td>—</td><td>5.20 (4.40–6.00)</td><td>1.87</td><td>0.23</td><td>†</td></tr><tr><td>Alpha</td><td>B.1.1.117</td><td>5.20 (4.87–5.47)</td><td>0.70</td><td>0.09</td><td>‡</td></tr><tr><td>Gamma</td><td>P.1</td><td>5.18 (4.93–5.43)</td><td>0.58</td><td>0.07</td><td>§</td></tr><tr><td>Delta</td><td>B.1.157.2</td><td>4.43 (4.36–4.49)</td><td>0.15</td><td>0.06</td><td>§</td></tr><tr><td rowspan="4">Omicron</td><td>BA.1</td><td>3.30 (2.80–3.70)</td><td>1.05</td><td>0.14</td><td>¶</td></tr><tr><td>BA.2</td><td>2.90 (2.70–3.10)</td><td>0.47</td><td>0.07</td><td>¶</td></tr><tr><td>BA.4</td><td>2.80 (1.50–6.70)</td><td>6.08</td><td>0.22</td><td>¶</td></tr><tr><td>BA.5+BQ.1</td><td>2.30 (1.60–3.10)</td><td>1.75</td><td>0.19</td><td>¶</td></tr></tbody></table>

**References:** †Chen et al. (2022), ‡Jusot (2022), §Galmiche et al. (2023), *Madewell et al. (2023).

<figure><img src="image_3.png" /><figcaption>**Fig. 1. A:**Variant prevalence for the COVID-19 pandemic in Italy; **B:** estimated number of infections attributable to each variant; **C:** same as in **B**, zoom on wild type and B.1.1.117 strains between end of September 2020 and mid-April 2021; **D:** visual representation of each variant's generation interval distribution.</figcaption></figure>
$$I_{j,1}(t) = \sigma_0 \sum_{k=1}^{N} z_{jk}(t) \mathcal{R}_k(t) \sum_{i=1}^{q} \hat{f}_i I_{k,i}(t-1) \quad (22a)$$

$$I_{j,i}(t) = \sigma_{i-1} I_{j,i-1}(t-1) \quad \text{for } i \ge 2, \quad (22b)$$

where, to allow inter-scenario comparisons, we compute each $\hat{f}_{k,i}$ term as the weighted average of the variant-specific $f_i^v$, with weights equal to the active infections (per node), that is

$$\hat{f}_{k,i} = \frac{\sum_{v=1}^{V} f_i^v I_{k,i}^v (t-1)}{\sum_{v=1}^{V} I_{k,i}^v (t-1)}. \quad (23)$$

# 3. Results

## 3.1. Computation of the variant-specific effective reproduction numbers

The spatially connected effective reproduction numbers ($\mathcal{R}$) of the considered variant strains are computed by running a sequential Monte Carlo algorithm independently for each circulating variant (Trevisin et al., 2023). We show in Fig. 2 both the local variant-specific $\mathcal{R}$ for each considered variant and, for the sake of completeness, the local variant-unstratified $\mathcal{R}$ (Trevisin et al., 2024). At the beginning of the spread of each strain, the variant-specific $\mathcal{R}$ take a value that is substantially larger than the unit threshold—the divide between endemic spread and disease-free conditions—thus indicating a rapid initial spread, where the number of new infections caused by a given variant undergoes exponential growth. After the variant has consolidated, its $\mathcal{R}$ decreases and starts oscillating around the unit threshold, a consequence of immunity build-up in the population and control measures. Interestingly, yet unsurprisingly, one may observe that, as soon as a new, fitter variant emerges, the $\mathcal{R}$ of the previous ones tend to steadily remain below the unit threshold until the old variants disappear or become undetectable, a purely mathematical consequence of the strain waning. In addition, one may note a larger variability among the local variant-specific values of $\mathcal{R}$ for the first strains (wild type, Alpha and Delta) than for later strains. These differences, however, remain small, and all nodes seem to react quite similarly to the circulation of SARS-CoV-2. The variant-unstratified local $\mathcal{R}$, on the other hand, appears much less prone to temporal fluctuations. Specifically, the initial jump observed in the variant-specific reproduction number is not matched by a substantial increase of the variant-unstratified index above the endemicity threshold $\mathcal{R} = 1$. Still, the rise in new recorded infections leads to a temporary growth of $\mathcal{R}$ above the endemicity threshold (e.g., before January 2021, or after July 2021).

## 3.2. Computation of the national epidemiological indices

The national $\mathcal{R}^G$ and epidemiicity index $\mathcal{E}^G$ are shown in Fig. 3 for both the variant-specific and variant-unstratified scenarios. As previously observed (Fig. 2), whenever a new variant emerges, its local variant-specific $\mathcal{R}$ values are characterized by an initial spike. This sudden jump is also seen if the global variant-specific $\mathcal{R}^G$ is computed, essentially the largest among the national variant-specific values of $\mathcal{R}$ (Fig. 3A).

The emergence of a new variant often brings the two time series (variant-specific and -unstratified global reproduction numbers) to the two different sides of the unit threshold, the variant-unstratified denoting disease-free equilibrium conditions. Out of the considered simulation span of 965 days, this occurred for 317 days (33 % of the time). In fact, while the variant-unstratified $\mathcal{R}^G$ only focuses on the total number of recorded infections, the variant-specific one detects whether any specific strain is growing. To better disentangle the behaviour of the variant-specific $\mathcal{R}^G$, we investigated the period during which the Alpha variant has replaced the wild type (see Fig. 3B) by computing the national $\mathcal{R}$ for both strains. At the emergence of the Alpha variant in October 2020, the new infections caused by this strain followed a steady exponential growth. Nonetheless, these additional infections did not shift the total number of recorded infections until January 2021 (Fig. 1C). This is reflected in the national $\mathcal{R}$ of the Alpha variant, which immediately diverged from the variant-unstratified national $\mathcal{R}$, which remained closer (and below the endemicity threshold) to the national $\mathcal{R}$ of the wild type, a waning strain, and only started growing above the endemicity threshold around March 2021.

One may observe larger increases for the first new variants (Alpha and Delta), at which emergence the $\mathcal{R}^G$ skyrockets to values up to three times higher. Conversely, later variants are characterized by minor, but still relevant spikes. The variant-unstratified $\mathcal{R}^G$ does not display this behaviour and settles around the endemicity threshold with only modest increases, a trend also shown by the local variant-unstratified $\mathcal{R}$. The two time series (variant-specific and unstratified), however, converge as soon as a variant establishes as the only circulating strain, as seen during the first phase of the epidemic up to October 2020 and, briefly, around autumn 2021.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Local effective reproduction numbers for different SARS-CoV-2 strains and different model set-ups during the COVID-19 pandemic in Italy. **A–H**: Variant-specific *R*; **I**: Variant-unstratified *R*.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** General epidemiological indices evaluated for different model set-ups. **A**: national effective reproduction number; dashed vertical lines highlight the temporal span shown in panel **B**, while filled diamonds mark the emergence (first recorded case) of a new variant and are colour-coded according to Fig. 1. **B**: as in **A**, showing national *R* for SARS-CoV-2 wild type and B.1.1.117 (here highlighted in red) between the end of September 2020 and mid-April 2021 against the variant-unstratified national *R*. **C**: as in **A**, but showing the epidemicity index; **D**: as in **C**, but between the end of September 2020 and mid-April 2021.</figcaption></figure>
The epidemiicity index (Fig. 3C) also holds a similar behaviour, with sudden spikes at the emergence of a new variant. Actually, the divergence between the variant-specific epidemiicity index and the variant-unstratified one appears even wider than the discrepancy observed for $\mathcal{R}$. This is apparent throughout the investigated timespan, but especially during the replacement of the wild type with the Alpha variant (Fig. 3D). Once again, the gap between the epidemiicity index of the newer strain and that of the one being replaced proves wider than the difference between the national reproduction numbers of the two variants. A higher epidemiicity index of the Alpha strain results from a lower epidemiicity threshold (see Table 1). In addition, $\mathcal{R}$ takes higher values (Fig. 3B). Owing to a different generation time's distribution, the threshold reproduction number $\mathcal{R}^*$ of the Alpha variant is much lower than that of the wild type, meaning that the control measure should target a much lower $\mathcal{R}$ to achieve unreactive conditions.

# 4. Discussion and conclusions

The emergence of new pathogenic variants during an epidemic induces shifts in the generation time distributions and in the values of reproduction numbers currently not accounted for by algorithms for the computation of epidemiological indices from data. Here, we have extended an existing sequential Monte Carlo algorithm for the real-time estimation of epidemiological indicators inclusive of variant-induced effects. Our main result is that adding such an additional layer of complexity is not a technicality nor a negligible and incremental effect, but rather a tool to provide significant corrections to the estimates of prognostic indices like e.g. the evolving values of the reproduction numbers from data. The epidemiicity index is likewise susceptible to major corrections owing to the emergence of viral variants. The two metrics are traditionally used to assess the fate of an outbreak of infectious disease in the long vs. short term, and in the immediate aftermath of a variant's emergence. This has significant implications, in the authors' view, on effective epidemic control.

To that end, we computed the variant-specific local $\mathcal{R}$ (Trevisin et al., 2023). Our computations suggest that the $\mathcal{R}$ values of an emerging variant typically overshoot well above the endemicity threshold, while both $\mathcal{R}$ of the circulating strain and of the variant-unstratified one merely fluctuate around (or below) the critical threshold. This is the expected behaviour when the strain is replaced, as its disappearance is directly related to the condition $\mathcal{R} < 1$ being verified for a prolonged period. On the other hand, the behaviour of the variant-unstratified $\mathcal{R}$ only focuses on the total number of recorded infections. Thus, $\mathcal{R}$ proves insensitive to whether one strain has begun to grow until the emergent strain overcomes the previously circulating ones, and drives and controls the total number of infections.

This behaviour also emerged with the global index $\mathcal{R}^G$, which accounts for all the variants spreading within a connected metapopulation. While the variant-unstratified $\mathcal{R}^G$ showcases the same behaviour of the local $\mathcal{R}$ (being only dependent on the total number of infections), the variant-specific $\mathcal{R}^G$ displays sudden spikes at the appearance of a new and evolutionary fitter strain. These jumps underline the presence of a variant growing more rapidly than the total number of infections, as was the case around November 2020 when the Alpha strain started replacing the wild type of SARS-CoV-2. This replacement went unnoticed in terms of the effective reproduction number until March 2021, when Alpha became the dominant strain and the variant-unstratified $\mathcal{R}^G$ finally started growing above the epidemic threshold. The variant-specific $\mathcal{R}^G$ anticipates this growth and highlights a potentially concerning situation well before a dominant strain is replaced. Failing to consider variants' stratification may thus lead to serious underestimations of the severity of an epidemic.

The epidemiicity index follows a qualitatively similar trend, and the variant-specific time series also diverges from the variant-unstratified at the emergence of newer and fitter strains. Similar spikes may be observed throughout the whole period studied. In addition to such spikes, the divergence between variant-specific and variant-unstratified time series appears substantially larger than the one observed with the reproduction numbers. This is a direct consequence of the higher reactivity of several strains than that of the wild type, characterized as it is by lower values of $\mathcal{R}^*$ (Table 1). For this reason, epidemic controls should ideally not simply aim at achieving $\mathcal{R}^G < 1$, but rather $\mathcal{R}^G < \mathcal{R}^*.$

Our computational experiments are not devoid of limitations. First, to implement the framework introduced here, one has to gather data on, or analytically derive, the number of new infections caused by each relevant strain. While, in practice, the chance for a variant to be sampled might be slightly different from its actual prevalence among the sequenced samples, this probability is unlikely to swiftly change over time and impact the computation of $\mathcal{R}$. Secondly, our framework requires knowledge of the variants' generation time distributions. This might prove challenging and in need of major assumptions at the emergence of a new variant (or at the emergence of the wild type of the pathogen) as reliable estimations of the serial interval and the generation time distribution usually require at least a few weeks of observations to be sorted out. As a result of its underlying algebra, the variant-stratified framework always provides more conservative estimates of the global epidemiological indices compared to the variant-unstratified ones. As such, even tentative estimates of the relevant generation time distributions can provide informative estimates of the values of the epidemiological metrics and help understand an evolving situation. Furthermore, strain replacement typically occurs in a few months, as suggested by the data shown in Fig. 1, and updated estimates of the generation time distributions may be used to adjust the inferred epidemiological numbers progressively. This implies that variant-specific frameworks point to the radically different description of the status of an outbreak – and its containment measures – than that of variant-unstratified ones. Inference of the variants' reproduction numbers may prove difficult during lull phases of the outbreak when the number of new infections for a specific variant is very low. While this is an established challenge for inference frameworks (Cori et al., 2013; Parag, 2021; Trevisin et al., 2023), working in discrete time allows one to exclude variants that do not achieve a steady prevalence or to aggregate data by adopting a larger time step, provided it is smaller than the mean of the generation time distribution. Working with a larger
time step may thus increase the prevalence over a given time window and improve our sequential Monte Carlo's estimates of inferred reproduction number time series. However, we maintain that the framework put together improves the state of the art in this field.

In conclusion, our results show that whenever variants' emergence is considered, embedding the pathogens' stratification into variants in algorithms for the computation of prognostic epidemiological indices produces substantially different results. In turn, this yields broad ramifications for optimal control measures. Our improved framework detects whether one specific strain grows more rapidly than the total number of infections and therefore foresees possible growths of new infections once the emerging variant becomes the dominant strain. Because our framework combines stratification in both space and variants, it provides a more precise estimate of the epidemiological metrics in a comprehensive environment. This is particularly relevant when outbreaks involve several communities and co-circulating pathogen strains. Finally, the sequential Monte Carlo process used here makes it possible to update the estimates of effective reproduction numbers and epidemicity indices as new surveillance data becomes available, with notable expediency. This framework may be applied to several practical scenarios, including when two or more airborne infectious diseases with similar infection routes (e.g., COVID-19, flu, etc.) circulate at the same time. It may be applied to any geographical setting, provided the data are available. Coupling this comprehensive framework with public health surveillance systems would create a pipeline that extends from data collection to the computation of the epidemiological indices, thereby securing updated indices that may be considered to plan control measures. Such integration would also make it possible to correct the values of these indices when needed, especially when data on variant prevalence is scarce, not collected daily, or needs frequent reassessments.

**CRediT authorship contribution statement**

**Cristiano Trevisin**: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Lorenzo Mari**: Writing – review & editing, Visualization, Validation, Investigation, Formal analysis, Data curation. **Marino Gatto**: Writing – review & editing, Validation, Methodology, Investigation, Formal analysis. **Vittoria Colizza**: Writing – review & editing, Visualization, Validation, Investigation, Formal analysis. **Andrea Rinaldo**: Writing – review & editing, Validation, Investigation, Funding acquisition, Formal analysis, Conceptualization.

**Data availability statement**

All codes and data involved in this study are available at the following Zenodo database https://doi.org/10.5281/zenodo.15011376 and are also available at the following public GitHub repository: https://github.com/cristianotrevisin/epidemiological-indices-metapopulation.

**Declaration**

All authors declare no competing interests.

**Declaration of competing interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Acknowledgements**

C.T. and A.R. acknowledge funding from the Swiss National Science Foundation via the project 'Optimal control of intervention strategies for waterborne disease epidemics' (grant number 200021–172578). L.M. acknowledges funding from the Italian Ministry of University and Research through the project 'Epidemiological data assimilation and optimal control for short-term forecasting and emergency management of COVID-19 in Italy' (FISR_2020IP_04249).

**References**

Allen, L. J., & van den Driessche, P. (2008). The basic reproduction number in some discrete-time epidemic models. *Journal of Difference Equations and Applications*, 14(10–11), 1127–1147.

Baccega, D., Castagno, P., Fernández Anta, A., & Sereno, M. (2024). Enhancing COVID-19 forecasting precision through the integration of compartmental models, machine learning and variants. *Scientific Reports*, 14(1), 1–17.

Bertuzzo, E., Mari, L., Pasetto, D., Miccoli, S., Casagrandi, R., Gatto, M., & Rinaldo, A. (2020). The geography of COVID-19 spread in Italy and implications for the relaxation of confinement measures. *Nature Communications*, 11(1), 4264.

Birello, P., Re Fiorentin, M., Wang, B., Colizza, V., & Valdano, E. (2024). Estimates of the reproduction ratio from epidemic surveillance may be biased in spatially structured populations. *Nature Physics*, 20(7), 1204–1210.

Carabelli, A. M., Peacock, T. P., Thorne, L. G., Harvey, W. T., Hughes, J., Peacock, S. J., Barclay, W. S., de Silva, T. I., Towers, G. J., & Robertson, D. L. (2023). SARS-CoV-2 variant biology: Immune escape, transmission and fitness. *Nature Reviews Microbiology*, 21(3), 162–177.
Chen, D., Lau, Y.-C., Xu, X.-K., Wang, L., Du, Z., Tsang, T. K., Wu, P., Lau, E. H. Y., Wallinga, J., Cowling, B. J., & Ali, S. T. (2022). Inferring time-varying generation time, serial interval, and incubation period distributions for COVID-19. *Nature Communications*, **13**(1), 1–12.

Colizza, V., Grill, E., Mikolajczyk, R., Cattuto, C., Kucharski, A., Riley, S., Kendall, M., Lythgoe, K., Bonsall, D., Wymant, C., Abeler-Dörner, L., Ferretti, L., & Fraser, C. (2021). Time to evaluate COVID-19 contact-tracing apps. *Nature Medicine*, **27**(3), 361–362.

Cori, A., Ferguson, N. M., Fraser, C., & Cauchemez, S. (2013). A new framework and software to estimate time-varying reproduction numbers during epidemics. *American Journal of Epidemiology*, **178**(9), 1505–1512.

Cuevas-Maraver, J., Kevrekidis, P., Chen, Q., Kevrekidis, G., & Drossinos, Y. (2024). Vaccination compartmental epidemiological models for the delta and omicron SARS-CoV-2 variants. *Mathematical Biosciences*, **367**, Article 109109.

Diekmann, O., Heesterbeek, J. a. P., & Roberts, M. G. (2010). The construction of next-generation matrices for compartmental epidemic models. *Journal of The Royal Society Interface*, **7**(47), 873–885.

Elbe, S., & Buckland-Merrett, G. (2017). Data, disease and diplomacy: GISAID's innovative contribution to global health. *Global Challenges*, **1**(1), 33–46.

Galmiche, S., Cortier, T., Charmet, T., Schaeffer, L., Chény, O., Platen, C. v., Lévy, A., Martin, S., Omar, F., David, C., Mailles, A., Carrat, F., Cauchemez, S., & Fontanet, A. (2023). SARS-CoV-2 incubation period across variants of concern, individual factors, and circumstances of infection in France: A case series analysis from the ComCor study. *The Lancet Microbe*, **4**(6), e409–e417.

Gatto, M., Bertuzzo, E., Mari, L., Miccoli, S., Carraro, L., Casagrandi, R., & Rinaldo, A. (2020). Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures. *Proceedings of the National Academy of Sciences*, **117**(19), 10484–10491.

Goldberg, E. E., Lin, Q., Romero-Severson, E. O., & Ke, R. (2023). Swift and extensive Omicron outbreak in China after sudden exit from 'zero-COVID' policy. *Nature Communications*, **14**(1), 3888.

Grabowski, F., Kochanćzyk, M., & Lipniacki, T. (2022). The spread of SARS-CoV-2 variant omicron with a doubling time of 2.0–3.3 Days can Be explained by immune evasion. *Viruses*, **14**(2), 294.

Guzman, M. G., & Harris, E. (2015). Dengue. *The Lancet*, **385**(9966), 453–465.

Jusot, J.-F. (2022). An update of serial interval estimates for COVID-19: A meta-analysis. *40pen*, **5**, 16.

Khare, S., Gurry, C., Freitas, L., Schultz, M. B., Bach, G., Diallo, A., Akite, N., Ho, J., Lee, R. T., Yeo, W., Team, G. C. C., & Maurer-Stroh, S. (2021). GISAID's role in pandemic response. *China CDC Weekly*, **3**(49), 1049–1051.

Kim, D., Ali, S. T., Kim, S., Jo, J., Lim, J.-S., Lee, S., & Ryu, S. (2022). Estimation of serial interval and reproduction number to quantify the transmissibility of SARS-CoV-2 omicron variant in South Korea. *Viruses*, **14**(3), 533.

Koelle, K., Martin, M. A., Antia, R., Lopman, B., & Dean, N. E. (2022). The changing epidemiology of SARS-CoV-2. *Science*, **375**(6585), 1116–1121.

Lemaitre, J. C., Pasetto, D., Zanon, M., Bertuzzo, E., Mari, L., Miccoli, S., Casagrandi, R., Gatto, M., & Rinaldo, A. (2022). Optimal control of the spatial allocation of COVID-19 vaccines: Italy as a case study. *PLoS Computational Biology*, **18**(7), Article e1010237.

Leung, X. Y., Islam, R. M., Adhami, M., Ilic, D., McDonald, L., Palawaththa, S., Diug, B., Munshi, S. U., & Karim, M. N. (2023). A systematic review of dengue outbreak prediction models: Current scenario and future directions. *PLoS Neglected Tropical Diseases*, **17**(2), Article e0010631.

Lin, J., Andreasen, V., & Levin, S. A. (1999). Dynamics of influenza A drift: The linear three-strain model. *Mathematical Biosciences*, **162**(1–2), 33–51.

Madewell, Z. J., Yang, Y., Longini, I. M., Halloran, M. E., Vespignani, A., & Dean, N. E. (2023). Rapid review and meta-analysis of serial intervals for SARS-CoV-2 Delta and Omicron variants. *BMC Infectious Diseases*, **23**(1), 1–9.

Mari, L., Casagrandi, R., Bertuzzo, E., Pasetto, D., Miccoli, S., Rinaldo, A., & Gatto, M. (2021). The epidemicity index of recurrent SARS-CoV-2 infections. *Nature Communications*, **12**(1), 2752.

Mari, L., Casagrandi, R., Rinaldo, A., & Gatto, M. (2017). A generalized definition of reactivity for ecological systems and the problem of transient species dynamics. *Methods in Ecology and Evolution*, **8**(11), 1574–1584.

Mari, L., Casagrandi, R., Rinaldo, A., & Gatto, M. (2018). Epidemicity thresholds for water-borne and water-related diseases. *Journal of Theoretical Biology*, **447**, 126–138.

Mari, L., Trevisin, C., Rinaldo, A., & Gatto, M. (2024). Sufficient reproduction numbers to prevent recurrent epidemics. *medRxiv*, page 2024.02.15.24302902.

Neubert, M. G., & Caswell, H. (1997). Alternatives to resilience for measuring the responses of ecological systems to perturbations. *Ecology*, **78**(3), 653–665.

Parag, K. V. (2021). Improved estimation of time-varying reproduction numbers at low case incidence and between epidemic waves. *PLoS Computational Biology*, **17**(9), Article e1009347.

Sabbatini, C. E., Pullano, G., Di Domenico, L., Rubrichi, S., Bansal, S., & Colizza, V. (2024). The impact of spatial connectivity on NPIs effectiveness. *BMC Infectious Diseases*, **24**(1), 21.

Saviano, M., Fierro, A., & Liccardo, A. (2023). A deterministic compartmental model for the transition between variants in the spread of Covid-19 in Italy. *PLoS One*, **18**(11), Article e0293416.

Shu, Y., & McAuley, J. (2017). Gisaid: Global initiative on sharing all influenza data – from vision to reality. *Euro Surveillance*, **22**(13), Article 30494.

Svensson, A. (2007). A note on generation times in epidemic models. *Mathematical Biosciences*, **208**(1), 300–311.

Torneri, A., Libin, P., Tomba, G. S., Faes, C., Wood, J. G., & Hens, N. (2021). On realized serial and generation intervals given control measures: The COVID-19 pandemic case. *PLoS Computational Biology*, **17**(3), Article e1008892.

Trevisin, C., Bertuzzo, E., Pasetto, D., Mari, L., Miccoli, S., Casagrandi, R., Gatto, M., & Rinaldo, A. (2023). Spatially explicit effective reproduction numbers from incidence and mobility data. *Proceedings of the National Academy of Sciences*, **120**(20), Article e2219816120.

Trevisin, C., Lemaitre, J. C., Mari, L., Pasetto, D., Gatto, M., & Rinaldo, A. (2022). Epidemicity of cholera spread and the fate of infection control measures. *Journal of The Royal Society Interface*, **19**(188), Article 20210844.

Trevisin, C., Mari, L., & Rinaldo, A. (2024). Epidemicity indices and reproduction numbers from infectious disease data in connected human populations. *Infectious Disease Modelling*, **9**(3), 875–891.

Van Den Driessche, P. (2017). Reproduction numbers of infectious disease models. *Infectious Disease Modelling*, **2**(3), 288–303.

van den Driessche, P., & Wmatough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences*, **180**(1), 29–48.

Zhang, M., Xiao, J., Deng, A., Zhang, Y., Zhuang, Y., Hu, T., Li, J., Tu, H., Li, B., Zhou, Y., Yuan, J., Luo, L., Liang, Z., Huang, Y., Ye, G., Cai, M., Li, G., Yang, B., Xu, B., ... Li, Y. (2021). Transmission dynamics of an outbreak of the COVID-19 delta variant B.1.617.2 — guangdong province, China, may–june 2021. *China CDC Weekly*, **3**(27), 584–586.