The epidemicity index of recurrent SARS-CoV-2 infections

Data and data processing

The modeling tools described in the following sections are applied to the Italian COVID-19 epidemic at the scale of second-level administrative divisions, i.e., provinces and metropolitan cities (as of 2020, 107 spatial units). Official data about resident population at the provincial level are produced yearly by the Italian National Institute of Statistics (Istituto Nazionale di Statistica, ISTAT; data available at http://dati.istat.it/Index.aspx?QueryId=18460). The January 2019 update has been used to inform the spatial distribution of the population.

The data to quantify nation-wide human mobility prior to the pandemic come from ISTAT (specifically, from the 2011 national census; data available online at https://www.istat.it/it/archivio/139381). Mobility fluxes, mostly reflecting commuting patterns related to work and study purposes, are provided at the scale of third-level administrative units (municipalities)^53,54. These fluxes were upscaled to the provincial level following the administrative divisions of 2019, and used to evaluate the fraction p_i of mobile people and the fraction q_ij of mobile people between i and all other administrative units j (see Supplementary Material in Gatto et al.⁷).

Airport traffic data for year 2019, used to inform the simulation shown in Fig. 4c, d, are from the Italian Airports Association (Assaeroporti; data available at http://assaeroporti.com/statistiche_201912/). Note that airports have been assigned to the main Metropolitan Area they serve, rather than to the province where they are geographically located (e.g., Malpensa Airport has been assigned to the Metropolitan City of Milano, rather than to the neighboring Varese province, where it actually lies).

Model parameters are taken from a paper by Bertuzzo et al.¹⁴, where they were inferred in a Bayesian framework on the basis of the official epidemiological bulletins released daily by Dipartimento della Protezione Civile⁵⁵ (data available online at https://github.com/pcm-dpc/COVID-19) and the bulletins of Epicentro, at ISS^51,56. The parameters estimated for the initial phase of the Italian COVID-19 epidemic¹⁴, during which SARS-CoV-2 was spreading unnoticed in the population, reflect a situation of unperturbed social mixing and human mobility, absent any effort devoted to disease control. This parameterization, in which all parameters (including the transmission rates) are spatially homogeneous, is reported in Table 2 and has been used to produce all the results presented in the main text, except for those of Fig. 6. In this case, to account for the containment measures put in place by the Italian authorities and their effects on transmission rates and mobility patterns during the first months of the pandemic, a time-varying parameterization¹⁴ for the period February 24 to May 1, 2020 has been used. In this parameterization, the transmission rates were allowed to take different values over different time windows, corresponding to the timing of the implementation of the main nation-wide restrictions, or lifting thereof. Specifically, the effect of the containment measures was parameterized by assuming that the transmission parameters had a sharp decrease after the containment measures announced at the end of February and the beginning of March, and that they were further reduced in the following weeks as the country was effectively entering full lockdown. As a by-product, these time-varying transmission rates can also at least partially account for seasonal effects on disease transmission. Due to the emerging nature of the pathogen, seasonality has not been given further consideration in this work; however, it may become a key component of future modeling efforts aimed at studying post-pandemic SARS-CoV-2 transmission dynamics³, i.e., if/when the pathogen establishes as endemic. Spatial connectivity too was modified with respect to the baseline scenario to reflect the disruption of mobility patterns induced by the pandemic and the associated containment measures¹⁴. Specifically, between-province mobility was progressively reduced as the epidemic unfolded according to estimates obtained through mobility data from mobile applications^53,57.

Spatially explicit SEPIAR with distributed controls

We consider a set of n communities connected by human mobility fluxes. In each community, the human population is subdivided according to infection status into the epidemiological compartments of susceptible, exposed (latently infected), post-latent (incubating infectious, also termed pre-symptomatic⁷), symptomatic infectious, asymptomatic infectious (including paucisymptomatic), and recovered individuals. The present model utilizes previous work aimed to describe the first wave of COVID-19 infections^7,14. In particular, it allows us to account for three widely adopted types of containment measures: reduction of local transmission (as a result of the use of personal protections, social distancing, and local mobility restriction), travel restriction, and isolation of infected individuals. To describe the effects of isolation, each infected compartment (exposed, post-latent, symptomatic and asymptomatic) is actually split into two, which allows keeping track of the abundances of infected individuals who are still in the community vs. those who are removed from it (i.e., either in isolation at a hospital, if symptomatic, or quarantined at home, if exposed, post-latent, or asymptomatic). The state variables of the model are summarized in Table 1. Supplementary Figure 1 recapitulates the structure of the model.

COVID-19 transmission dynamics are thus described by the following set of ordinary differential equations:

Susceptible individuals are recruited into community i (i = 1…n) at a constant rate μN_i, with μ and N_i being the average mortality rate of the population and the size of the community in the absence of disease, respectively, and die at rate μ. In this way, the equilibrium size of community i without disease amounts to N_i. Susceptible individuals get exposed to the pathogen at rate λ_i, corresponding to the force of infection for community i (detailed below), thus becoming latently infected (but not infectious yet). Exposed individuals die at rate μ and transition to the post-latent, infectious stage at rate δ^E. If containment measures including mass testing and preventive isolation of positive cases are in place, exposed individuals may be removed from the general population and quarantined at rate ({chi }_{i}^{E}). Post-latent individuals die at rate μ, progress to the next infectious classes at rate η^P, developing an infection that can be either symptomatic—with probability σ—or asymptomatic, including the case in which only mild symptoms are present—with probability 1 − σ, and may be tested and quarantined at rate ({chi }_{i}^{P}). Symptomatic infectious individuals die at rate μ + α, with α being an extra-mortality term associated with disease-related complications, recover from infection at rate γ^I, may spontaneously seek treatment at a hospital at rate η, and may be identified through mass screening and hospitalized at rate ({chi }_{i}^{I}). Asymptomatic individuals die at rate μ, recover at rate γ^A, and may be quarantined at rate ({chi }_{i}^{A}). Infected individuals who are either hospitalized or quarantined at home are subject to the same epidemiological dynamics as those who are still in the community, but are considered to be effectively removed from it, thus not contributing to disease transmission. Individuals who recover from the infection die at rate μ, and are assumed to have permanent immunity to reinfection. This last assumption is not fundamental, as loss of immunity can be easily included in the model. However, immunity to SARS-CoV-2 reinfection is reported to be relatively long-lasting (a few months at least), hence its loss cannot alter transmission dynamics over epidemic timescales¹⁴.

The cornerstone of model (Eq. (3)) is the force of infection, λ_i, which in a spatially explicit setting must account not only for locally acquired infections but also for the role played by human mobility. We assume that, at the spatiotemporal scales of interest for our problem, human mobility mostly depicts daily commuting flows (also coherently with the data available for parameterization; see above) and does not actually entail a permanent relocation of individuals. We thus describe human mobility (and the associated social contacts possibly conducive to disease transmission) by means of instantaneous spatial-mixing matrices ({M}_{c,ij}^{X}) (with X ∈ {S, E, P, I, A, R}), i.e.,

where p_i (0 ≤ p_i ≤ 1 for all i’s) is the fraction of mobile people in community i, q_ij (0 ≤ q_ij ≤ 1 for all i’s and j’s) represents the fraction of people moving between i and j (including j = i, (mathop{sum }nolimits_{j = 1}^{n}{q}_{ij}=1) for all i’s), r^X (0 ≤ r^X ≤ 1 for all X’s) quantifies the fraction of contacts occurring while individuals in epidemiological compartment X are traveling, and ξ_ij (0 ≤ ξ_ij ≤ 1 for all i’s and j’s) represents the effects of travel restrictions that may be imposed between any two communities i and j as a part of the containment response. Therefore, the probability that residents from i have social contacts while being in j (independently of with whom) is assumed to be proportional to the fraction r^X of the mobility-related contacts of the individuals in epidemiological compartment X, multiplied by the probability p_i that people from i travel (independently of the destination) and the probability q_ij that the travel occurs between i and j, possibly reduced by a factor 1 − ξ_ij accounting for travel restrictions. All other contacts contribute to mixing within the local community (i in this case). Note also that if ξ_ij = 0 for all i’s and j’s, then ({M}_{c,ij}^{X}) reduces to ({M}_{ij}^{X}), i.e., to the mixing matrix in the absence of disease-containment measures. In this case, (mathop{sum }nolimits_{j = 1}^{n}{M}_{ij}^{X}=1) for all i’s and X’s. It is important to remark, though, that the epidemiologically relevant contacts between the residents of two different communities, say i and j, may not necessarily occur in either i or j; in fact, they could happen anywhere else, say in community k, between residents of i and j simultaneously traveling to k. On this basis, we define the force of infection as

where the parameters ({beta }_{j}^{X}) (X ∈ {P, I, A}) are the community-dependent rates of disease transmission from the three infectious classes, ϵ_j (0 ≤ ϵ_j ≤ 1 for all j’s) represents the reduction of transmission induced by social distancing, the use of personal protective equipment, and local mobility restrictions if such containment measures are in fact in place, and the terms ({M}_{c,ij}^{X}) (with X ∈ {S, E, P, I, A, R}) describe the epidemiological effects of mobility between i and j in the presence of disease-containment measures. Note that transmission has been assumed to be frequency-dependent.

The parameters μ, δ^X (X ∈ {E, P}), σ, α, η, γ^X (X ∈ {I, A}), and r^X (X ∈ {S, E, P, I, A, R}) are assumed to be community-independent, for they pertain to population demography at the country scale or the clinical course of the disease. By contrast, the transmission rates ({beta }_{i}^{X}) (X ∈ {P, I, A}) and the control parameters, namely the isolation rates ({chi }_{i}^{X}) (X ∈ {E, P, I, A}), the reductions of transmission due to personal protection, social distancing, and local mobility restriction ϵ_i, and the travel restrictions ξ_ij, are assumed to be possibly community-dependent, thereby reflecting spatial heterogeneities in disease transmission prior to the implementation of containment measures (({beta }_{i}^{X})), testing effort and/or strategy (({chi }_{i}^{X})), local transmission reduction (ϵ_i), and travel restriction (ξ_ij).

Derivation of the basic and control reproduction numbers

Close to the DFE, a state in which all individuals are susceptible to the disease (S_i = N_i, with N_i being the baseline population size of community i) and all the other epidemiological compartments are empty (({E}_{i}={P}_{i}={I}_{i}={A}_{i}={E}_{i}^{{rm{q}}}={P}_{i}^{{rm{q}}}={I}_{i}^{{rm{h}}}={A}_{i}^{{rm{q}}}={R}_{i}=0) for all i’s), the dynamics of model (Eq. (3)) is described by the linearized system (dot{{bf{x}}}={{bf{J}}}_{{bf{c}}}{bf{x}}), where ({bf{x}}={[{S}_{i},{E}_{i},{P}_{i},{I}_{i},{A}_{i},{E}_{i}^{{rm{q}}},{P}_{i}^{{rm{q}}},{I}_{i}^{{rm{h}}},{A}_{i}^{{rm{q}}},{R}_{i}]}^{T}) (where i = 1…n and the superscript T denotes matrix transposition) and J_c is the spatial Jacobian matrix

where I and 0 are the identity and null matrices of size n, respectively, ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{X}}}) (X ∈ {E, P, I, A}) are diagonal matrices whose non-zero elements are (mu +{delta }^{E}+{chi }_{i}^{E}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}})), (mu +{delta }^{P}+{chi }_{i}^{P}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})), (mu +alpha +eta +{gamma }^{I}+{chi }_{i}^{I}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})), and (mu +{gamma }^{A}+{chi }_{i}^{A}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})), and the matrices ({{boldsymbol{theta }}}_{{bf{c}}}^{{bf{X}}}) (X ∈ {P, I, A}) are given by

where N is a diagonal matrix whose non-zero elements are the population sizes N_i, ({{bf{M}}}_{{bf{c}}}^{{bf{X}}}=[{M}_{c,ij}^{X}]) (X ∈ {S, P, I, A}) are sub-stochastic matrices representing the spatially explicit contact terms in the presence of containment measures, ϵ is a diagonal matrix whose non-zero entries are the transmission reductions ϵ_i, β^X (X ∈ {P, I, A}) are diagonal matrices whose non-zero elements are the contact rates ({beta }_{i}^{X}), and Δ_c is a diagonal matrix whose non-zero entries are the elements of vector ({bf{u}}{bf{N}}{{bf{M}}}_{{bf{c}}}^{{bf{S}}}), with u being a unitary row vector of size n.

Because of its block-triangular structure, it is immediate to see that J_c has 6n strictly negative eigenvalues, namely −μ, with multiplicity 2n, and −(μ + δ^E),−(μ + δ^P), −(μ + α + γ^I), and −(μ + γ^A), each with multiplicity n. Therefore, the asymptotic stability properties of the DFE of model (Eq. (3)), which determine whether long-term disease circulation in the presence of controls is possible, are linked to the eigenvalues of a reduced-order spatial Jacobian associated with the infection subsystem, i.e., the subset of state variables directly related to disease transmission, in this case {E₁, …, E_n, P₁, …, P_n, I₁, …, I_n, A₁, …, A_n}. Note that introducing waning immunity would not change the spectral properties of the Jacobian matrix evaluated at the DFE. The reduced-order Jacobian ({{bf{J}}}_{{bf{c}}}^{* }) thus reads

The asymptotic stability properties of the DFE can be assessed through a NGM approach^22,37. In fact, the spectral radius of the NGM provides an estimate of the so-called control reproduction number⁵⁸, ({{mathcal{R}}}_{{rm{c}}}), which can be thought of as the average number of secondary infections produced by one infected individual in a completely susceptible population in the presence of disease-containment measures. Clearly, if ({{mathcal{R}}}_{{rm{c}}},> , 1) the pathogen can invade the population in the long run, and endemic transmission will eventually be established despite the implementation of disease-containment measures. To evaluate ({{mathcal{R}}}_{{rm{c}}}) for model (Eq. (3)), the Jacobian of the infection subsystem can be decomposed into a spatial transmission matrix

and a transition matrix

so that J_c = T_c + Σ_c. The spatial NGM with large domain ({{bf{K}}}_{{bf{c}}}^{{bf{L}}}), including variables other than the states-at-infection⁵⁹ (i.e., the exposed individuals E_i) thus reads

with

Because of the peculiar block-triangular structure of ({{bf{K}}}_{{bf{c}}}^{{bf{L}}}), the spatial NGM with small domain (K_c, accounting only for E_i) is simply ({{bf{K}}}_{{bf{c}}}^{{bf{1}}}) (see again Diekmann et al.⁵⁹). The control reproduction number can thus be found as the spectral radius of the NGM (with either large or small domain), i.e.,

where

are three spatially explicit generation matrices describing the contributions of post-latent infectious people, infectious symptomatic people, and asymptomatic/paucisymptomatic infectious people to the next generation of infections in a neighborhood of the DFE in the presence of disease-containment measures.

In the absence of controls, i.e., if the isolation rates ({chi }_{i}^{X}) (X ∈ {E, P, I, A}), the transmission reductions ϵ_i, and the travel restrictions ξ_ij are equal to zero for all i’s and j’s, then the control reproduction number ({{mathcal{R}}}_{{rm{c}}}) reduces to the basic reproduction number ({{mathcal{R}}}_{0}), defined as the average number of secondary infections produced by one infected individual in a population that is completely susceptible to the disease and where no containment measures are in place. ({{mathcal{R}}}_{0}) can be evaluated as the spectral radius of matrix G^P + G^I + G^A, where

In the previous set of expressions, ϕ^X (X ∈ {E, P, I, A}) are diagonal matrices whose non-zero elements are μ + δ^E (for ϕ^E), μ + δ^P (for ϕ^P), μ + α + η + γ^I (for ϕ^I), and μ + γ^A (for ϕ^A), while matrices θ^X (X ∈ {P, I, A}) are given by ({bf{N}}{{bf{M}}}^{{bf{S}}}{{boldsymbol{beta }}}^{{bf{X}}}{({boldsymbol{Delta }})}^{-1}{({{bf{M}}}^{{bf{X}}})}^{T}), with ({{bf{M}}}^{{bf{X}}}=[{M}_{ij}^{X}]) (X ∈ {S, P, I, A}) and ({M}_{ij}^{X}={M}_{c,ij}^{X}) evaluated with ξ_ij = 0 for all i’s and j’s, and Δ is a diagonal matrix whose non-zero entries are the elements of vector uNM^S.

Derivation of basic and control epidemicity indices

The concept of epidemicity²⁶ extends previous work^24,25 where a reactivity index was defined and applied to study the transient dynamics of ecological systems characterized by steady-state behavior. To explain, in physical terms, the meaning of reactivity and of the Hermitian matrix used to derive it, consider a linear system dx/dt = Ax, where ({bf{x}}={({x}_{1},ldots ,{x}_{n})}^{T}) is the state vector and A is a n × n real state matrix. The system is subject to pulse perturbations x(0) = x₀ > 0. Reactivity is defined as the gradient of the Euclidean norm (| | {bf{x}}| | =sqrt{{x}_{1}^{2}+cdots +{x}_{n}^{2}}=sqrt{{{bf{x}}}^{T}{bf{x}}}) of the state vector, evaluated for the fastest-growing initial perturbation, and corresponds to the spectral abscissa ({{{Lambda }}}_{max }^{{rm{Re}}}(cdot )) of the Hermitian part (A + A^T)/2 of matrix A²⁴. Following Mari et al.²⁵, an asymptotically stable equilibrium is characterized by positive generalized reactivity if there exist small perturbations that can lead to a transient growth in the Euclidean norm of a suitable system output y = Wx, with matrix W describing a linear transformation of the system state.

In epidemiological applications, W should include the variables of the infection subsystem²⁶. Therefore, a suitable output transformation for the problem at hand is

where w^E, w^P, w^I, w^A are the weights assigned to the variables of the infection subsystem in the output ({bf{y}}=[{w}^{E}{E}_{1},ldots ,{w}^{E}{E}_{n},{w}^{P}{P}_{1},ldots ,{w}^{P}{P}_{n},{w}^{I}{I}_{1},ldots ,{w}^{I}{I}_{n},{w}^{A}{A}_{1},ldots ,{w}^{A}{A}_{n}]^{T}). Generalized reactivity for the DFE of system (Eq. (3)) is positive if the spectral abscissa of a suitable Hermitian matrix (either H₀ or H_c, depending on whether the spread of disease is uncontrolled or some containment measures are in place) is also positive. In SEPIAR, the expressions of matrices H₀ and H_c are far from trivial, as shown below, and the evaluation of spectral abscissae typically requires numerical techniques. Note also that, since recovered individuals are not accounted for in the system output, including waning immunity would not alter the epidemicity properties of the DFE.

Let us consider the most general case of disease-containment measures being in place (which includes as a limit case also uncontrolled pathogen spread). If we note that (ker ({bf{W}})=ker ({bf{W}}{{bf{J}}}_{{bf{c}}})), with J_c being the Jacobian of SEPIAR at the DFE in the presence of controls, matrix H_c can be defined^25,27 as the Hermitian part of WJ_c(W)⁺, i.e.,

where (W)⁺ is the right pseudo-inverse (a generalization of the concept of inverse for non-square matrices) of W, and can be evaluated as

Matrix

is Hermitian, hence real and symmetric. Therefore all eigenvalues are real and the spectral abscissa ({e}_{{rm{c}}}={{{Lambda }}}_{max }^{{rm{Re}}}({{bf{H}}}_{{bf{c}}})) coincides with the largest eigenvalue, which corresponds to the fastest-growing perturbation in the system output. Thus, e_c can be interpreted as a control epidemicity index: if e_c > 0, there must exist some small perturbations to the DFE that are temporarily amplified in the system output, thus generating a transient, subthreshold epidemic wave.

Absent any containment measures, the control epidemicity index, e_c, reduces to the basic epidemicity index, ({e}_{0}={{{Lambda }}}_{max }^{{rm{Re}}}({{bf{H}}}_{{bf{0}}})), where

and the Jacobian matrix J₀ can be obtained from J_c by setting equal to zero the isolation rates ({chi }_{i}^{X}) (X ∈ {E, P, I, A}), the transmission reductions ϵ_i, and the travel restrictions ξ_ij for all i’s and j’s.

The effective reproduction number and the effective epidemicity index

The reproduction numbers and the epidemicity indices defined above can be rigorously applied only to characterize the spread of disease in a fully naïve population (S_i = N_i ∀ i). As soon as the pathogen begins to circulate within the population, the state of the system gradually departs from the DFE. Under these circumstances, it is customary^19,21 to define a time-dependent, effective reproduction number, ({mathcal{R}}(t)), to track the number of secondary infections caused by a single infectious individual in a population in which the pool of susceptible individuals is progressively depleted, and control measures are possibly in place⁵⁸. Similarly, it is possible to define an effective epidemicity index, e(t), to evaluate the likelihood that transient epidemic waves may occur even if ({mathcal{R}}(t),<,1).

The definition of these time-dependent metrics requires to update the expressions of the spatially explicit infection matrices ({{boldsymbol{theta }}}_{{bf{c}}}^{{bf{X}}}) (X ∈ {P, I, A}) in a time-varying fashion, i.e.,

where ϵ(t) is a diagonal matrix whose non-zero elements represent the reduction of local transmission rates at time t, ϵ_i(t), ({{boldsymbol{Delta }}}_{{bf{c}}}(t)={rm{diag}}({bf{u}}{sum }_{Xin {S,E,P,I,A,R}}{bf{X}}(t){{bf{M}}}_{{bf{c}}}^{{bf{X}}}(t))), ({{bf{M}}}_{{bf{c}}}^{{bf{X}}}(t)) are spatially explicit contact matrices including time-varying travel restrictions, ξ_ij(t), and S(t), E(t), P(t), I(t), A(t), and R(t) are diagonal matrices whose non-zero elements are the time-varying abundances of susceptible (S_i(t)), exposed (E_i(t)), post-latent (_Pi(t)), symptomatic (I_i(t)), asymptomatic (A_i(t)), and recovered (R_i(t)) individuals in each community i = 1…n. The evaluation of ({mathcal{R}}(t)) and e(t) also mandates an update of the transition matrices ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{X}}}) to include time-dependent testing and isolation rates, ({chi }_{i}^{X}) (X ∈ {E, P, I, A}).

Computing the ({{boldsymbol{theta }}}_{{bf{c}}}^{{bf{X}}}) and ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{X}}}) matrices thus requires updating the state variables and epidemiological parameters of SEPIAR, i.e., to numerically solve Eq. (3) as the epidemic progresses and control strategies are put in place. The definition of time-varying transmission and transition terms gives rise to time-varying Jacobians and NGMs. In turn, the use of these matrices in the computation of reproduction numbers and epidemicity indices allows the evaluation of the time-varying threshold quantities ({mathcal{R}}(t)) and e(t). Clearly, this type of argument works best if the depletion of the susceptible pool is relatively slow, i.e., if the initial perturbation is sufficiently small and the divergence of the epidemic trajectory from the DFE is not too large^60,61.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com