The epidemicity index of recurrent SARS-CoV-2 infections
Data and data processingThe 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 pi of mobile people and the fraction qij of mobile people between i and all other administrative units j (see Supplementary Material in Gatto et al.7).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.14, where they were inferred in a Bayesian framework on the basis of the official epidemiological bulletins released daily by Dipartimento della Protezione Civile55 (data available online at https://github.com/pcm-dpc/COVID-19) and the bulletins of Epicentro, at ISS51,56. The parameters estimated for the initial phase of the Italian COVID-19 epidemic14, 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 parameterization14 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 dynamics3, 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 measures14. Specifically, between-province mobility was progressively reduced as the epidemic unfolded according to estimates obtained through mobility data from mobile applications53,57.Spatially explicit SEPIAR with distributed controlsWe 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-symptomatic7), symptomatic infectious, asymptomatic infectious (including paucisymptomatic), and recovered individuals. The present model utilizes previous work aimed to describe the first wave of COVID-19 infections7,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:$${dot{S}}_{i} =mu ({N}_{i}-{S}_{i})-{lambda }_{i}{S}_{i}\ {dot{E}}_{i} ={lambda }_{i}{S}_{i}-(mu +{delta }^{E}+{chi }_{i}^{E}){E}_{i}\ {dot{P}}_{i} ={delta }^{E}{E}_{i}-(mu +{delta }^{P}+{chi }_{i}^{P}){P}_{i}\ {dot{I}}_{i} =sigma {delta }^{P}{P}_{i}-(mu +alpha +{gamma }^{I}+eta +{chi }_{i}^{I}){I}_{i}\ {dot{A}}_{i} =(1-sigma ){delta }^{P}{P}_{i}-(mu +{gamma }^{A}+{chi }_{i}^{A}){A}_{i}\ {dot{E}}_{i}^{{rm{q}}} ={chi }_{i}^{E}{E}_{i}-(mu +{delta }^{E}){E}_{i}^{{rm{q}}}\ {dot{P}}_{i}^{{rm{q}}} ={chi }_{i}^{P}{P}_{i}+{delta }^{E}{E}_{i}^{{rm{q}}}-(mu +{delta }^{P}){P}_{i}^{{rm{q}}}\ {dot{I}}_{i}^{{rm{h}}} =(eta +{chi }_{i}^{I}){I}_{i}+sigma {delta }^{P}{P}_{i}^{{rm{q}}}-(mu +alpha +{gamma }^{I}){I}_{i}^{{rm{h}}}\ {dot{A}}_{i}^{{rm{q}}} ={chi }_{i}^{A}{A}_{i}+(1-sigma ){delta }^{P}{P}_{i}^{{rm{q}}}-(mu +{gamma }^{A}){A}_{i}^{{rm{q}}}\ {dot{R}}_{i} ={gamma }^{I}({I}_{i}+{I}_{i}^{{rm{h}}})+{gamma }^{A}({A}_{i}+{A}_{i}^{{rm{q}}})-mu {R}_{i}.$$
(3)
Susceptible individuals are recruited into community i (i = 1…n) at a constant rate μNi, with μ and Ni 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 Ni. 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 timescales14.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.,$${M}_{c,ij}^{X}=left{begin{array}{ll}{r}^{X}{p}_{i}{q}_{ij}(1-{xi }_{ij})hfill&,{text{if}},i,ne, jhfill\ (1-{p}_{i})+(1-{r}^{X}){p}_{i}+{r}^{X}{p}_{i}{q}_{ij}(1-{xi }_{ij})&,{text{if}},i=j,end{array}right.$$
(4)
where pi (0 ≤ pi ≤ 1 for all i’s) is the fraction of mobile people in community i, qij (0 ≤ qij ≤ 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), rX (0 ≤ rX ≤ 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 rX of the mobility-related contacts of the individuals in epidemiological compartment X, multiplied by the probability pi that people from i travel (independently of the destination) and the probability qij 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$${lambda }_{i}=mathop{sum }limits_{j=1}^{n}{M}_{c,ij}^{S}frac{(1-{epsilon }_{j})left({beta }_{j}^{P}mathop{sum }nolimits_{k = 1}^{n}{M}_{c,kj}^{P}{P}_{k}+{beta }_{j}^{I}mathop{sum }nolimits_{k = 1}^{n}{M}_{c,kj}^{I}{I}_{k}+{beta }_{j}^{A}mathop{sum }nolimits_{k = 1}^{n}{M}_{c,kj}^{A}{A}_{k}right)}{mathop{sum }nolimits_{k = 1}^{n}left({M}_{c,kj}^{S}{S}_{k}+{M}_{c,kj}^{E}{E}_{k}+{M}_{c,kj}^{P}{P}_{k}+{M}_{c,kj}^{I}{I}_{k}+{M}_{c,kj}^{A}{A}_{k}+{M}_{c,kj}^{R}{R}_{k}right)},$$
(5)
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 rX (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 numbersClose to the DFE, a state in which all individuals are susceptible to the disease (Si = Ni, with Ni 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 Jc is the spatial Jacobian matrix$${{bf{J}}}_{{bf{c}}}=left[begin{array}{llllllllll}-mu {bf{I}}&{bf{0}}&-{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&-{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&-{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{delta }^{E}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{{boldsymbol{chi }}}^{{bf{E}}}&{bf{0}}&{bf{0}}&{bf{0}}&-(mu +{delta }^{E}){bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{{boldsymbol{chi }}}^{{bf{P}}}&{bf{0}}&{bf{0}}&{delta }^{E}{bf{I}}&-(mu +{delta }^{P}){bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&eta {bf{I}}+{{boldsymbol{chi }}}^{{bf{I}}}&{bf{0}}&{bf{0}}&sigma {delta }^{P}{bf{I}}&-(mu +alpha +{gamma }^{I}){bf{I}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{{boldsymbol{chi }}}^{{bf{A}}}&{bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-(mu +{gamma }^{A}){bf{I}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{gamma }^{I}{bf{I}}&{gamma }^{A}{bf{I}}&{bf{0}}&{bf{0}}&{gamma }^{I}{bf{I}}&{gamma }^{A}{bf{I}}&-mu {bf{I}}end{array}right],$$
(6)
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$${{boldsymbol{theta }}}_{{bf{c}}}^{{bf{X}}}={bf{N}}{{bf{M}}}_{{bf{c}}}^{{bf{S}}}({bf{I}}-{boldsymbol{epsilon }}){{boldsymbol{beta }}}^{{bf{X}}}{({{boldsymbol{Delta }}}_{{bf{c}}})}^{-1}{({{bf{M}}}_{{bf{c}}}^{{bf{X}}})}^{T},$$
(7)
where N is a diagonal matrix whose non-zero elements are the population sizes Ni, ({{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 Jc 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 {E1, …, En, P1, …, Pn, I1, …, In, A1, …, An}. 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$${{bf{J}}}_{{bf{c}}}^{* }=left[begin{array}{llll}-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}\ {delta }^{E}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&{bf{0}}&{bf{0}}\ {bf{0}}&sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}\ {bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}end{array}right].$$
(8)
The asymptotic stability properties of the DFE can be assessed through a NGM approach22,37. In fact, the spectral radius of the NGM provides an estimate of the so-called control reproduction number58, ({{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$${{bf{T}}}_{{bf{c}}}=left[begin{array}{llll}{bf{0}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}end{array}right],$$
(9)
and a transition matrix$${{boldsymbol{Sigma }}}_{{bf{c}}}=left[begin{array}{llll}-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&{bf{0}}&{bf{0}}&{bf{0}}\ {delta }^{E}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&{bf{0}}&{bf{0}}\ {bf{0}}&sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}\ {bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}end{array}right],$$
(10)
so that Jc = Tc + Σc. The spatial NGM with large domain ({{bf{K}}}_{{bf{c}}}^{{bf{L}}}), including variables other than the states-at-infection59 (i.e., the exposed individuals Ei) thus reads$${{bf{K}}}_{{bf{c}}}^{{bf{L}}}=-{{bf{T}}}_{{bf{c}}}{({{mathbf{Sigma }}}_{{bf{c}}})}^{-1}=left[begin{array}{llll}{{bf{K}}}_{{bf{c}}}^{{bf{1}}}&{{bf{K}}}_{{bf{c}}}^{{bf{2}}}&{{bf{K}}}_{{bf{c}}}^{{bf{3}}}&{{bf{K}}}_{{bf{c}}}^{{bf{4}}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}end{array}right],$$
(11)
with$${{bf{K}}}_{{bf{c}}}^{{bf{1}}} ={delta }^{E}left[{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}+sigma {delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}+(1-sigma ){delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}right]{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}})}^{-1}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})}^{-1}\ {{bf{K}}}_{{bf{c}}}^{{bf{2}}} =left[{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}+sigma {delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}+(1-sigma ){delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}right]{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})}^{-1}\ {{bf{K}}}_{{bf{c}}}^{{bf{3}}} ={{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}\ {{bf{K}}}_{{bf{c}}}^{{bf{4}}} ={{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}.$$
(12)
Because of the peculiar block-triangular structure of ({{bf{K}}}_{{bf{c}}}^{{bf{L}}}), the spatial NGM with small domain (Kc, accounting only for Ei) is simply ({{bf{K}}}_{{bf{c}}}^{{bf{1}}}) (see again Diekmann et al.59). The control reproduction number can thus be found as the spectral radius of the NGM (with either large or small domain), i.e.,$${{mathcal{R}}}_{{rm{c}}}=rho ({{bf{K}}}_{{bf{c}}}^{{bf{L}}})=rho ({{bf{K}}}_{{bf{c}}})=rho ({{bf{G}}}_{{bf{c}}}^{{bf{P}}}+{{bf{G}}}_{{bf{c}}}^{{bf{I}}}+{{bf{G}}}_{{bf{c}}}^{{bf{A}}}),$$
(13)
where$${{bf{G}}}_{{bf{c}}}^{{bf{P}}} ={delta }^{E}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})}^{-1}\ {{bf{G}}}_{{bf{c}}}^{{bf{I}}} =sigma {delta }^{E}{delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}\ {{bf{G}}}_{{bf{c}}}^{{bf{A}}} =(1-sigma ){delta }^{E}{delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}$$
(14)
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 GP + GI + GA, where$${{bf{G}}}^{{bf{P}}} ={delta }^{E}{{boldsymbol{theta }}}^{{bf{P}}}{({{boldsymbol{phi }}}^{{bf{E}}}{{boldsymbol{phi }}}^{{bf{P}}})}^{-1}\ {{bf{G}}}^{{bf{I}}} =sigma {delta }^{E}{delta }^{P}{{boldsymbol{theta }}}^{{bf{I}}}{({{boldsymbol{phi }}}^{{bf{E}}}{{boldsymbol{phi }}}^{{bf{P}}}{{boldsymbol{phi }}}^{{bf{I}}})}^{-1}\ {{bf{G}}}^{{bf{A}}} =(1-sigma ){delta }^{E}{delta }^{P}{{boldsymbol{theta }}}^{{bf{A}}}{({{boldsymbol{phi }}}^{{bf{E}}}{{boldsymbol{phi }}}^{{bf{P}}}{{boldsymbol{phi }}}^{{bf{A}}})}^{-1}.$$
(15)
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 uNMS.Derivation of basic and control epidemicity indicesThe concept of epidemicity26 extends previous work24,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) = x0 > 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 + AT)/2 of matrix A24. Following Mari et al.25, 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 subsystem26. Therefore, a suitable output transformation for the problem at hand is$${bf{W}}=left[begin{array}{llllllllll}{bf{0}}&{w}^{E}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{w}^{P}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{w}^{I}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{w}^{A}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}end{array}right],$$
(16)
where wE, wP, wI, wA 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 H0 or Hc, 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 H0 and Hc 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 Jc being the Jacobian of SEPIAR at the DFE in the presence of controls, matrix Hc can be defined25,27 as the Hermitian part of WJc(W)+, i.e.,$${{bf{H}}}_{{bf{c}}}=H({bf{W}}{{bf{J}}}_{{bf{c}}}{({bf{W}})}^{+})=frac{1}{2}left{{bf{W}}{{bf{J}}}_{{bf{c}}}{({bf{W}})}^{+}+{[{({bf{W}})}^{+}]}^{T}{({{bf{J}}}_{{bf{c}}})}^{T}{({bf{W}})}^{T}right},$$
(17)
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$${({bf{W}})}^{+}={({bf{W}})}^{T}{[{bf{W}}{({bf{W}})}^{T}]}^{-1}.$$
(18)
Matrix$${{bf{H}}}_{{bf{c}}}=left[begin{array}{llll}-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&frac{{w}^{P}}{2{w}^{E}}{delta }^{E}{bf{I}}+frac{{w}^{E}}{2{w}^{P}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&frac{{w}^{E}}{2{w}^{I}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&frac{{w}^{E}}{2{w}^{A}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}\ frac{{w}^{P}}{2{w}^{E}}{delta }^{E}{bf{I}}+frac{{w}^{E}}{2{w}^{P}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&frac{{w}^{I}}{2{w}^{P}}sigma {delta }^{P}{bf{I}}&frac{{w}^{A}}{2{w}^{P}}(1-sigma ){delta }^{P}{bf{I}}\ frac{{w}^{E}}{2{w}^{I}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&frac{{w}^{I}}{2{w}^{P}}sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}\ frac{{w}^{E}}{2{w}^{A}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}&frac{{w}^{A}}{2{w}^{P}}(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}end{array}right]$$
(19)
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, ec can be interpreted as a control epidemicity index: if ec > 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, ec, reduces to the basic epidemicity index, ({e}_{0}={{{Lambda }}}_{max }^{{rm{Re}}}({{bf{H}}}_{{bf{0}}})), where$${{bf{H}}}_{{bf{0}}}=H({bf{W}}{{bf{J}}}_{{bf{0}}}{({bf{W}})}^{+})=frac{1}{2}left{{bf{W}}{{bf{J}}}_{{bf{0}}}{({bf{W}})}^{+}+{[{({bf{W}})}^{+}]}^{T}{({{bf{J}}}_{{bf{0}}})}^{T}{({bf{W}})}^{T}right}$$
(20)
and the Jacobian matrix J0 can be obtained from Jc 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 indexThe 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 (Si = Ni ∀ 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 customary19,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 place58. 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), More
