Derivation of the effective reproduction number ℛ for COVID-19 in relation to mobility restrictions and confinement

Supplementary Note 1. Epidemic spreading model

We propose a tailored model for the epidemic spread of COVID-19. We use a previous framework for the study of epidemics in structured metapopulations, with heterogeneous agents, subjected to recurrent mobility patterns [7, 17, 10, 18].To understand the geographical diffusion of the disease, as a result of human-human interactions in small geographical patches, one has to combine the contagion process with the long-range disease propagation due to human mobility across different spatial scales. For the case of epidemic modeling, the metapopulation scenario is as follows. A population is distributed in a set of patches, being the size (number of individuals) of each patch in principle different. The individuals within each patch are well-mixed, i.e., pathogens can be transmitted from an infected host to any of the healthy agents placed in the same patch with the same probability. The second aspect of our metapopulation model concerns the mobility of agents. Each host is allowed to change its current location and occupy another patch, thus fostering the spread of pathogens at the system level. Mobility of agents between different patches is usually represented in terms of a network where nodes are locations while a link between two patches represents the possibility of moving between them.

We introduce a set of modifications to the standard metapopulation model to account for the different states relevant for the description of COVID-19, and also to substitute the well-mixing with a more realistic set of contacts. Another key point is the introduction of a differentiation of the course of the epidemics that depends on the demographic ages of the population. This differentiation is very relevant in light of the observation of a scarcely set of infected individuals at ages (< 25), and also because of the severe situations reported for people at older ages (> 65). Our model is composed of the following epidemiological compartments: susceptible (S), exposed (E), asymptomatic infectious (A), infected (I), hospitalized to ICU (H), dead (D), and recovered (R). Additionally, we divide the individuals in N_G age strata, and suppose the geographical area is divided in N regions or patches. Although we present the model in general form, its application to COVID-19 only makes use of the three age groups mentioned above (N_G = 3): young people (Y), with age up to 25; adults (M), with age between 26 and 65; and elderly people (O), with age larger than 65.

Let us suppose we have a population of N individuals distributed in N_P regions, with n_i individuals residing in region i. We also consider that individuals belong to one of N_G different age strata, in such a way that individuals of age group g live in region i. Thus, where n^g is the total population of group age g.

Our system is completely characterized with the fraction of individuals in state m, where m ∈{S, E, A, I, H, D, R }, for each age stratum g, and associated with each patch i, that we denote in the following as . The temporal evolution of these quantities is given by:

These equations correspond to a discrete-time dynamics, in which each time-step represents a day. They are built upon previous work on Microscopic Markov-Chain Approach (MMCA) modelization of epidemic spreading dynamics [6], but which has also been applied to other types of processes, e.g., information spreading and traffic congestion [8, 9, 16].

The rationale of the model is the following. Susceptible individuals get infected by contacts I with asymptomatic and infected agents, with a probability , becoming exposed. Exposed individuals turn into asymptomatic at a certain rate η^g, which in turn become infected at a rate α^g. Once infected, two paths emerge, which are reached at an escape rate µ^g. The first option is requiring hospitalization in an ICU, with a certain probability γ^g; otherwise, the individuals become recovered. While being at ICU, individuals have a death probability ω^g, which is reached at a rate ψ^g. Finally, ICUs discharge at a rate χ^g, leading to the recovered compartment. See in Supplementary Figure 1 a sketch of the compartmental epidemic model with all the transitions, and Supplementary Table 1 for a list of the parameters and their values to simulate the spreading of COVID-19 in Spain, which will be discussed in Supplementary Note 6.

The value of encodes the probability that a susceptible agent belonging to age group g and patch i contracts the disease. Under the model assumptions, this probability is given by: where p^g denotes the degree of mobility of individuals within age group is the mobility matrix (fraction of individuals of group age g that choose destination j while living in region i), and denotes the probability that those agents get infected by the pathogen inside patch i. This way, the first term in the r.h.s. of Eq. (S9) denotes the probability of contracting the disease inside the residence patch, whereas the second term contains those contagions taking place in any of the neighboring areas. Furthermore, we assume that the number of contacts increases with the density of each area according to a monotonously increasing function f . Finally, we introduce an age-specific contact matrix, C, whose elements C^gh define the fraction of contacts that individuals of age group g perform with individuals belonging to age group h. With the above definitions, reads

Parameters β_A and β_I correspond to the infection probabilities for contacts of a susceptible individual with asymptomatic and symptomatic infected individuals, respectively. The exponents represent the number of contacts made by an agent of age group g in patch i with infectious individuals —compartments A and I, respectively — of age group h residing at patch j. Accordingly, the double product expresses the probability for an individual belonging to age group g not being infected while staying in patch i. Finally, parameter ν is introduced to account for the self-isolation of symptomatic infected individuals, either because of awareness or health state, which effectively reduces their appearance in the patches.

The term z^g⟨k^g⟩f (ñ_i/s_i) in Eq. (S10) represents the overall number of contacts (infectious or non infectious), which increases with the density of patch i following function f, where s_i is the area of patch i measured in km², and also accounts for the normalization factor z^g, which is calculated as: where the effective population at patch i, i.e., the number of people present at patch i when commuting takes place, is given by which is distributed in age groups of size

For convenience, we have defined the effective mobility matrices , that take into account both the degree of mobility of the population, and the transition probabilities to the neighboring patches.

From now on, we will use the tilde to indicate variables measured while the commuting is active, to distinguish them from variables when all the population is in its residence patch.

The function f(x) governing the influence of population density has been selected, following [11], as:

The last term of the exponents in Eq. (S10) contains the probability that these contacts are contagious, which is proportional to , the expected number of individuals of age group h in the given infectious state m (either A or I) which have moved from region k to region i:

The discrete time nature of this model allows for an easy computation of the time evolution of all the relevant variables, providing information at the regional level. See Supplementary Note 7 for the details of its application to the COVID-19 outbreak in Spain. Additionally, the model is amenable for analytical inspection, which has allowed us to find the epidemic threshold, see Supplementary Note 3.

Supplementary Note 2. Modeling containment measures

Here we assess the performance of different containment measures to reduce the impact of COVID-19 using the mathematical model. To incorporate containment policies in our formalism, we consider that a given fraction of the adult population κ₀ is isolated at home whereas both young and elderly people are assume to stay at home . In this sense, let us remark that parameter κ₀ allows us to tune the strength of those containment mesaures proposed to contain COVID-19. In this sense, κ₀ = 1 reflects a total lockdown of the population which isolates the households from each other, thus constraining the transmission dynamics at the level of household rather than municipality. From the former assumptions, we compute the average number of contacts of agents belonging to each group g ∈ {Y, M, O} as where σ encodes the average household size.

In this scenario, a relevant indicator to quantify the efficiency of the policy is the probability of one individual living in a household, inside a given municipality i, without any infected individual. Assuming that containment is implemented at time t_c, this quantity, denoted in the following as CH _i(t_c), is given by

This way, Eqs. (S17) and (S18) become time-dependent: where Θ(x) is the Heaviside function, that is 1 if x ≥ 0 and 0 otherwise. Accordingly, the mobility parameters p^g change as which make and z^g also dependent on time, see Eqs. (S11)–(S13). Unfortunately, complete isolation of households is impossible, thus we introduce a permeability factor ϕ that reduces the effective isolation.

This containment strategy is introduced in the dynamical Eqs. (S2)–(S7) by modifying Eqs. (S2) and (S3) for the time after t_c: where we have added a new compartment CH to hold the individuals under household isolation after applying containment κ₀, and δ_a,b is the Kronecker function, which is 1 if a = b and 0 otherwise. Containment also affects the average number of contacts, thus we must also update Eq. (S10):

In summary, we have introduced a containment strategy characterized by the isolation of a certain fraction of the population, which in turn reduces the average number of contacts, the mobility and, as a result, the probability of becoming infected, thus reducing the overall prevalence of the disease. The containment is parameterized with the containment factor κ₀, and the permeability factor ϕ.

Supplementary Note 3. Calculation of ℛ₀ using the next generation matrix

The model is amenable for analytical calculations. We calculate the basic reproduction number ℛ₀ using the next generation matrix (NGM) approach [4]. Accordingly, we need to analyze the stability of the disease free equilibrium. We do so by making a first order expansion of the model equations for small values (ϵ) of the non-susceptible states m: and , where m ∈ {E, A, I, H, D, R}. The expansion allows us to transform our discrete time Markov Chain into a system of continuous time differential equations. We start by expanding the infection probabilities : where we have defined and

We then insert the above expression into , leading to: where it is convenient to separate tensor J if four terms:

With

These tensors encode the four different ways in which the epidemic interactions may take place: individuals belonging to the same patch i = j and not moving (J₁); interaction in the patch of i with individuals coming from patch j (J₂); interaction in the patch of j with individuals coming from patch i (J₃); and individuals from i and j interacting at any other patch k (J₄).

In the next generation matrix framework, we only need to consider the epidemic compartments. Additionally, the framework requires to work in absolute numbers of infected individuals. Accordingly, we define the absolute number of individuals in the respective compartments as

Making use of the above definitions, the corresponding differential equations take the form:

Defining the vector N^g =(N^E,g, N^A,g, N^I,g)T, the above system of differential equations can be rewritten as: where we have defined with, and

In the above tensor we have defined

With the above differential equation, the basic reproduction number is given by:

We can calculate the inverse of the tensor V as

The inverse of the matrix V ^g is given by:

Accordingly, we have:

As we are looking for the eigenvectors of the tensor FV ⁻¹, we note that their components associated to the compartments A and I —rows 2 and 3— must be zero, since the associated rows in the above matrix are zero. To be more precise, we have for a = 2N_P +1, … 3N_P, which are the elements associated to the compartments A and I. Accordingly, we can restrict the above matrix only to the vector space associated to the compartment E and the eigenvalues will be equivalent, which gives us the basic reproduction number ℛ₀:

Where and Λ_max(Z) denotes the spectral radius of tensor Z.

Supplementary Note 4. Effective reproduction number and critical containment

The effective reproduction number ℛ is given by where ρ^S(t) is the fraction of reachable susceptible individuals at time t and ℛ₀ is the basic reproduction number [3]. In our case, the containment has an effect on both terms. First, κ₀ reduces the pool of individuals susceptible of contracting the disease, i.e. it reduces the term ρ^S. And, second, κ₀ reduces the average degree of the individuals and thus it reduces ℛ₀. Eq. (S50) allows us to accurately determine the effective reproduction number and capture the impact of human mobility and containment measures on the spread of the disease.

Computing Eq. (S50) by means of the NGM approach involves the calculation of the spectral radius of tensor Z, which hinders an understanding of the roots behind the transition from flattening to bending the epidemic curve triggered by an efficient confinement, as illustrated in the main text. In order to have an analytical estimation of the critical value of the needed confinement to observe this transition, , let us consider a single well mixed population. In this case, we have and

If we remove the age strata, the basic reproduction number is reduced to where we have approximated which is valid for low values of the infection probabilities β_m. The same result for ℛ₀, Eq. (S53), is obtained if all age groups share the values of population n^g, average number of contacts ⟨k^g⟩, asymptomatic infectious rate ⟨α^g⟩, and escape rate ⟨µ^g⟩, since Λ_max(C) = 1 due to the contact-by-age matrix C being a stochastic matrix.

Accordingly, after introducing containment, we get for the effective reproduction number: where the temporal dependence has been omitted. The parameter ϕ encodes the permeability (household isolation level), ranging from ϕ = 0 for which households are perfectly isolated, to ϕ = 1 which recovers the well-mixed scenario before the intervention, but now with a lower average number of contacts.

In order to bend the epidemic curve and not to flatten it we impose the condition ℛ ≤ 1. Accordingly, the critical containment value must fulfill the equation:

Defining the following coefficients: we can write Eq. (S56) as whose solution is given by:

Note that only the physically meaningful solution, , must be retained.

Supplementary Note 5. Dynamical approximation of the effective reproduction number

A naive approach to computing the effective reproduction number ℛ is considering the mean field scenario, in which an infected subject i contacts ⟨k⟩ individuals each time step. Assuming an infection probability β, the expected number of individuals infected by i at each time step is ρ_S ⟨k⟩β, where ρ_S is the fraction of susceptible individuals on the population. Taking into account the duration of the infectious period τ, we can estimate how many individuals have been infected by subject i over time as: ℛ = ρ_S ⟨k⟩βτ .

We can extend this approach to include both the temporal dependence and the heterogeneity of contacts induced by daily commuting patterns, making use of all the probabilities provided by our model equations. Moreover, here we aim at studying the temporal evolution of the effective reproduction number beyond the early stage of the disease. In this sense, we define the reproduction number of patch i and age group g as the number of infections observed if we seed an infected individual in the aforementioned patch and age group [14]. This involves considering that this infected individual, with residence in patch i, may commute to patch j, where it will be able to contact and infect susceptible individuals coming from any other patch k.

First, we calculate the expected number of susceptible individuals of age group h which have moved from region k to region j as:

The fraction of susceptible individuals of each age group g in patch j at time t is expressed by:

Next, we compute the number of susceptible contacts made by an individual of age group g and patch i, which can be expressed as:

Using this expression, and assuming that we have an individual from patch i and age group g that has become infectious, we can compute the effective reproduction number for each patch and age group as:

For estimating the global effective reproduction number, we then make use of a weighted average taking into account the distribution of the population across patches and age groups:

In Eq. (S65), we assume that the fraction of susceptible individuals as well as the average number of contacts is constant during the infectious period. However, the fraction of susceptible individuals obviously changes during that time. More importantly, as containment measures are put in place, the number of contacts varies as well. To account for the temporal variability of these quantities, we calculate the contributions to of the infections produced since the individual became asymptomatically infected at time t. Since the infectious period is divided in two phases, one asymptomatic with infectious rate α^g, and another symptomatic with infectious rate µ^g, what we have is a coupled Poisson process. The evolution of the probabilities to be asymptomatic, ζ^A,g, or symptomatic, ζ^I,g, are given by:

Since the seed is initially placed as asymptomatic, the initial conditions are given by ζ^A,g(0) = 1 and ζ^I,g(0) = 0. Accordingly, solving explicitly the recurrent equations, we find:

Incorporating these probabilities into Eq. (S65) we obtain

Supplementary Note 6. Parameters for the modelization of the spreading of COVID-19

In this subsection, we detail our parameters choice to study the current epidemic outbreak in Spain. Regarding epidemiological parameters, the incubation period has been reported to be η⁻¹ + α⁻¹ = 5.2 days [12] in average which, in our formalism, must be distributed into the exposed and asymptomatic compartments. In principle, if one does not expect asymptomatic transmissions, most of this time should be spent inside the exposed compartment, thus being the asymptomatic infectious compartment totally irrelevant for disease spreading. However, along the line of several recent works [13, 15, 5] we have found that the unfolding of COVID-19 cannot be explained without accounting for infections from individuals not developing any symptoms previously. In particular, our best fit to reproduce the evolution of the real cases reported so far in Spain yields α⁻¹ = 2.86 days as asymptomatic infectious period. In turn, the infectious period while being symptomatic is established as µ⁻¹ = 7 days [1, 20], except for the young strata, for which we have reduced it to 1 day, assigning the remaining 6 days as asymptomatic; this is due to the reported mild symptoms in young individuals, which may become inadvertent [2]. Furthermore, we have included an isolation factor ν = 0.6 which reduces the infectiousness of symptomatic agents as a result of their self-isolation. We fix the fatality rate ω = 42% of ICU patients by studying historical records of dead individuals as a function of those requiring intensive care. In turn, we estimate the period from ICU admission to death as ψ⁻¹ = 7 days [19] and the stay in ICU for those overcoming the disease as χ⁻¹ = 20 days [1].

Supplementary Note 7. Calibration and validation of the model for the evolution of COVID-19 in Spain

Equations (S2)–(S8) enable to monitor the spatio-temporal propagation of COVID-19 across Spain. To check the validity of our formalism, we aggregate the number of cases predicted for each municipality at the level of autonomous regions (comunidades autónomas), which is a first-level political and administrative division, and compare them with the number of cases daily reported by the Spanish Health Ministry. In this sense, we compute the number of cases predicted for each municipality i at each time step t as:

As our model is designed to predict the emergence of autochthonous cases triggered by local contagions and commuting patterns, those imported infected individuals corresponding to the first reported cases in Spain are initially plugged into our model as asymptomatic infectious agents. In addition, small infectious seeds should be also placed in those areas where anomalous outbreaks have occurred due to singular events such as one funeral in Vitoria leading to more than 60 contagions. Overall, the total number of infectious seeds is 47 individuals which represents 0.2 % of the number of cases reported by March 20, 2020.

Supplementary Figure 2 shows that our model is able to accurately predict not only the overall evolution of the total number of cases at the national scale but also their spatial distribution across the different autonomous regions. Moreover, the most typical trend observed so far is an exponential growth of the number of cases, thus clearly suggesting that the disease was spreading freely in most of the territories when the model was calibrated. Note, however, that there are some exceptions such as La Rioja or País Vasco in which some strong policies targeting the most affected areas were previously promoted to slow down COVID-19 propagation.

We also show in Supplementary Figure 3 the distribution of at three different stages of the epidemic spreading in Spain, according to our model: free spreading without containment (March 14, 2020); 10 days after the implementation of partial confinement, showing the initial effects of containment (March 24, 2020); and 21 days after the implementation of partial confinement (April 4, 2020).

Supplementary Note 8. Sensitivity Analysis

We have performed a sensitivity analysis to check the robustness of the critical value of the effective reproduction number ℛ = 1, identifying the transition from supercritical (flattening) regime to the subcritical (bending) one when applying confinement. For this purpose, we report here the results obtained by tuning some of the parameters of the model under the same assumptions and show that the obtained results are consistent with the results reported in the main text of the manuscript and do not alter the conclusions of our study.

In particular, we fix the permeability parameter to ϕ = 0.20 and the confinement to κ₀ = 0.66, and sample a region of the parameters space by varying the average household size, σ ∈ [2, 3], the asymptomatic infectious rate of each strata g, α^g ∈ [0.25, 0.5], and the self-isolation parameter of symptomatic infected individuals, ν ∈ [0.5, 0.7].

Results of this sensitivity analysis are depicted in Figure 4, where we show that the effective reproduction number remains almost unaltered while exploring different points of the parameters space.