key: cord-0925608-8vbfqp5i
authors: Steinegger, B.; Granell, C.; Rapisardi, G.; Gomez, S.; Matamalas, J.; Soriano-Panos, D.; Gomez-Gardenes, J.; Arenas, A.
title: Retrospective study of the first wave of COVID-19 in Spain: analysis of counterfactual scenarios
date: 2021-02-18
journal: nan
DOI: 10.1101/2021.02.16.21251832
sha: e6e3d4485eb51faa9353c3c719ae2277ab46cbb1
doc_id: 925608
cord_uid: 8vbfqp5i

One of the most important questions on the COVID-19 pandemic is ascertaining the correct timing to introduce non-pharmaceutical interventions (NPIs), based mainly on mobility restrictions, to control the rising of the daily incidence in a specific territory. Here, we make a retrospective analysis of the first wave of the epidemic in Spain and provide a set of useful insights to optimize actions in the near future. We have reconstructed the exposure times, from infection to detectability, to correctly estimate the reproduction number Rt. This enables us to analyze counterfactual scenarios to understand the impact of earlier or later responses, decoupling containment measures from natural immunity. Our results quantify the differences in the number of fatalities for earlier and later responses to the epidemic in Spain.

impact of the various NPIs implemented at the timing they were implemented. The answer is essential to design the future mitigation strategies of COVID-19, as it has been done for China or Italy (2, 3) .

To assess the impact of NPIs, one can either directly analyze the epidemiological data such as case numbers, hospitalizations, and fatalities (4) or perform model-based inference (5, 6) to extract cornerstone conclusions. Here, we blend both approaches. Analyzing the epidemiological data allows us to propose a hypothesized functional form for Rt during Spain's first wave. Then, we combine this information with a model-based approach to infer the behavior of Rt. With this analysis, we can explore counterfactual scenarios.

One of the most pressing difficulties when analyzing an epidemic's impact retrospectively is estimating the unobserved past incidence of the infection (7) . Estimating the past incidence of COVID-19 relies on the accurate determination of the period's distribution between the exposure to the virus to affliction (or detectable symptomatic infection).

In Spain, the ISCIII provides a time series with symptom onset dates for reported cases (8). Using current estimates of the incubation period, we can infer the daily infections using a non-parametric back-projection method based on a maximum likelihood estimation (7) (see Materials and methods).

In Fig. 1A we show the results of the estimated infections for the Comunidades Autónomas (CCAA) in Spain, as well as for the whole country, together with the corresponding official data for symptoms onset and reporting date. The substantial delays between infection and report are evident for all regions, as corroborated in Fig. 1B showing the day when the infections and reported cases peaked, andthe difference between them. On average, infections and reported cases peaked with a 15 days difference (sd: 3). We find the maximal difference of 20 days for Catalonia, together with Castilla La-Mancha.

Then, we estimated Rt's evolution in Fig. 2A using EpiEstim (9) . In all CCAA, we observe a substantial reduction of Rt before the first confinement. Close to this first confinement, Rt is for the first time below one in almost all CCAA, and stayed mainly below one until the end of the confinement. Fig. 2B shows the day on which Rt was first below one in the CCAA. From March 13 to 21, Rt became smaller than 1 in all CCAA. On average, Rt crosses 1 on March 16 with a standard deviation of 2,1 days, which is consistent with other estimations (14, 15) . The Community of Madrid, was the first CCAA to push Rt below 1 as a consequence of the implemented containment measures on March 9, before any other CCAA. Similar results are found for other CCAA (see SM for detailed information about the evolution across CCAA).

After obtaining the retrospective values of Rt, and given that the NPIs impact starts with the behavior of people, we wonder when the epidemic response was initiated. The epidemic response should result in a breakpoint (17) in Rt after which we observe a substantial reduction of Rt. To determine this breakpoint, R1, we separate (see Materials and Methods) Rt's evolution into three linear segments: the "free" spreading phase before restrictions, the initiation of the epidemic response, and the lockdown. The result is shown in Fig. 3A . The second breakpoint pretty much coincides with the implementation of the first lockdown by the authorities in Spain. More interestingly, the reproduction number starts already to decrease between March 5 and 6 (first breakpoint). The early decrease precedes the introduction of any containment measures, also on a regional level, as well as the mobility reduction, which started between March 9 and 10 (see Fig.  S1 ).

However, this early decrease is found based on the reported case data. Therefore, the estimation of Rt is subjected to changes in the testing capacity and reporting criteria. The preceding decrease of Rt concerning mobility and containment measures could be due to a saturation of the test capacity, as illustrated in Fig. 3B for the CCAA Catalonia and Aragon. To corroborate this hypothesis, we compare the estimation with the most reliable epidemiological data in Spain, daily fatalities. In particular, we use our estimation to inform a simple, but robust, model (Methods Materials) about Rt's general form, and adjust it to the daily fatalities. The first breakpoint, BP, which indicates the start of the epidemic response, as well as the initial number of infected individuals, I0, are left as free parameters. To reduce the number of free parameters, we fix the second breakpoint on March 15, one of two break points in the time series for Rt, Fig. 3A . Additionally, we assume a constant value R2 of Rt during the lockdown. Figure S2 presents the summary statistics. The corresponding evolution of Rt is shown in Fig.4A . The estimations of R1 and R2 lead to a doubling time and half-life time of 1.84(CI: 1.69-1.99) days and 11.2(CI: 10.8-11.6) days, respectively (27) .

In Fig. 4A , we compare the Rt estimated from official incidence data, with the one fitted to reported daily fatalities. The first observation is that the actual Rt seems to have been substantially higher than what can be estimated from the reported incidence data. Additionally, we observe a reduction in Rt between March 10 and 13 substantially later than inferred from the reported data. Importantly, this indicates a decrease in Rt , preceding the implementation of the national lockdown. This seems to be the fingerprint of a spillover effect (16) from the earlier introduction of local restrictions as in Madrid, for example. These earlier restrictions raised the individual awareness of the population as indicated by an increase of COVID-19 related searches in Google from March 9 onward (see Fig. S3 ). Figure 4C shows the estimated ascertainment bias according to the ratio between the daily infections from the model and the exposure times. Our results indicate that, initially, ascertainment was only around 10%. Additionally, ascertainment continuously deteriorated towards the first confinement, when it almost reached 5%. This decreasing ascertainment indicates that the ramping up of the test capacity could not keep up with the propagation of the epidemic.

Finally, the inferred form of Rt serves us to build counterfactual scenarios. We evaluate the accumulated fatalities, as well as the attack rate, if the epidemic response (including confinement, an earlier adoption of social distancing and the specific measures taken by the CCAA) had occurred earlier or later by shifting the function of Rt.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint Figure 4B presents the total attack rate and accumulated fatalities if the entire epidemic response had occurred some days earlier/later. For a 7 days earlier response, we find an attack rate of 1.2% (CI: 1.1-1.4%) and a total number of fatalities of 4,782 (CI: 4,146-5,503). In contrast, supposing the epidemic response 7 days later, we find an attack rate of 30.8% (CI: 28.2-34.0%) and a total number of fatalities of 120,353 (CI: 110,165-132,545). Due to the high Rt in Spain, which translates into a minimal doubling time, any delay in the epidemic response would have had drastic consequences.

Our results indicate the existence of a substantial delay between exposure and notification in Spain. An average delay of 15 days proves how late individuals with symptoms were tested, the substantial test turnaround time, and the significant delay to be notified. Such pronounced delay hinders an actual evaluation of the epidemiological situation and therefore impedes effective management of the epidemic response. For illustration, as the first confinement was implemented, our results indicate that the daily infections had already peaked. Similarly, on March 28, all nonessential economic activity was shut down despite still rising case numbers. However, we find that daily infections had already peaked more than 10 days before in most CCAA.

Besides the homogeneous, impactful containment measures implemented on a national level, we found that the earlier introduction of containment measures on a regional level led to differences when the daily case numbers first stopped growing for upto a week. In general, the CCAA,which first pushed Rt below 1,coincides with the ones that implemented containment measures the earliest. The most prominent example is the Community of Madrid. Despite the regional heterogeneity, the first confinement on March 15 fell into the range when Rt became smaller than 1.

By fitting a hypothesized functional form of Rt to reproduce the evolution of daily fatalities, we found that Rt started decreasing between March 10 and 13. The first containment measures taken from March 9 onwards and the corresponding mobility reduction which started around March 10, fall into this range. This reduction of Rt, before the first confinement, highlights the role of individual awareness of the population, a factor whose consideration is crucial (16) . Nonetheless, the observation based on reconstructed infections contrasts with the reported data, which indicates a constant decrease in Rt from March 5 onwards. Available data from Aragón and Catalunya indicates that this earlier decrease can be attributed to saturation in the testing capacity. Similarly, the higher Rt value inferred from the model suggests that the expansion of test capacity could not keep up with the propagation speed of the epidemic in Spain. Accordingly, constantly fewer cases were detected towards lockdown before eventually only capturing approximately 5% of cases. In contrast, Althaus et al. find the opposite for Switzerland, where an increasing percentage of cases was detected towards lockdown (15) . A decreasing detection rate, a substantial delay between exposure and notification, together with a doubling time of about 1.8 days, hampered the evaluation of Spain's epidemiological situation and may have led to an underestimation regarding the severity of the situation.

Our counterfactual scenarios indicate that only a minimal delay in the epidemic response leads to a substantially different outcome. For an epidemic response 7 days earlier, we find approximately . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint 5,000 fatalities and an attack rate of 1%, while a 7 days later response results in 120,000 fatalities and an attack rate of 30%. In this sense, the epidemic response and NPIs prevented many fatalities. Nevertheless, an earlier introduction of these measures would have substantially reduced COVID-19 related mortality. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Solid lines indicate a centered, 7 day rolling average. In red we show the reconstructed exposure times. The shaded area in light gray indicates the lockdowns 1 and 3. In dark gray, we indicate the lockdown 2 where, in addition, all non essential economic activity was shut down. (B) Date on which the reconstructed exposure (red) times and reported cases (blue) peaked. Additionally, we explicitly indicate the delay between the two peaks. Shaded areas and bars around lines indicate 95% credible intervals in the entire manuscript. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This PDF file includes:

Figures S1 to S13 References 1 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint

To reconstruct the exposure times in the different CCAA, we use the function backprojNP (17) from the surveillance package (18, 19) in R. The method, initially proposed by Becker et al. (7), infers the expected number of exposures given the probability mass function of the incubation period through a maximum likelihood approach.

Credible intervals are calculated based on a bootstrap procedure (20) . We fixed the smoothing factor k = 6, which corresponds to a centered rolling average of 7 days. The bootstrap procedure made use of 1000 samples (B = 1000). The incubation period distribution was fixed as a gamma distribution with mean 5.2 days and standard deviation 2.8 days (21) . The time series with the symptom onset is provided by the Centro Nacional de Epidemilogía (22).

From the median incidence, obtained by the reconstruction of the exposure times, we estimate the evolution of R t using the EpiEstim package (9, 23) . Our interest lies in the impact of the NPI. Accordingly, we estimate the instantaneous R t (24) rather than the case R t (25), since the former reflects changes in the transmission dynamics more abruptly.

For the infectivity function, we choose the generation time estimated by Ganyani et al. (12) . To be more precise, we assume a generation time following a gamma distribution with mean 5.2 (CI : 3.78 − 6.78) days and standard deviation 1.72 (CI : 0.91 − 3.93) days. This generation time distribution corresponds to the estimation by Ganynai et al. for Singapore while assuming the same incubation period distribution (21) as we do in the reconstruction of the exposure times. We assume a standard deviation for the mean and the standard deviation of the generation time of 1.0 days and 1.2 days, respectively. However, we bound the values for the mean and standard deviation by the estimations of Ganyani et al..

We fixed a centered rolling average of 7 days. To bootstrap the credible intervals, we took 100 samples of the generation time distribution and considered 100 posteriors for each of these samples (n1 = 100, n2 = 100).

To identify the linear segments we use the R package segmented (13, 26) . The method proposed by Muggeo implements a maximum likelihood approach using linear predictors. The credible intervals are obtained through bootstrap. As previously pointed out, we assume three segments of R t . An initial constant value R1 as the disease was spreading freely in Spain as well as two linear evolving parts corresponding to the decrease of R t towards lockdown and the constant decrease observed during the lockdown.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02. 16.21251832 doi: medRxiv preprint Adjusting the model to the daily fatalities

To evaluate the validity of our estimation of R t , we fit a piece-wise linear R t to the epidemiological data in Spain. We then compare this fit with our estimation of R t from the exposure times. Unfortunately, the only reliable time series available are the daily fatalities (27) . Both, hospitalizations and ICU occupation, were subject to various changes of reporting criteria from March until May 2020. As an example, CCAA switched between reporting the accumulated number of admissions and the current number of ICUs occupied without announcement. For this reasons, we refrain from using hospitalization data and focus on the daily fatalities.

We assume R t to have three segments as in the previous analysis. Motivated by the identified segments and to reduce the number of parameters to fit, we fix the second breakpoint on March 15, the day the lockdown was implemented. It seems reasonable that, as the first lockdown was implemented, R t reached a relatively stable value. In contrast, the first breakpoint BP, which initiates the decrease of R t and thus the epidemic response, we leave as a free parameter. In this way, we are able to evaluate to which extent the decrease of R t before lockdown actually took place. Additionally, we assume a constant value R1 and R2 of R t before and after the second breakpoint.

We then inform an epidemic model with the assumed form of R t . Given R t , the generation time distribution w(t) and the size of the population N , the daily incidence I t on day t evolves as (9, 24) :

In the above equation, we decoupled the effect of containment measures and a reduction of R t due to natural immunity. This separation serves us later to build the counterfactual scenarios. We use the same generation time distribution as for estimating R t (12) . From the daily incidence, we propagate the symptom onset as well as the daily fatalities through convolution (28) . Given the incubation time distribution P (t) as well as the time from symptom onset to decease D(t), the daily number of individuals developing symptoms S t and the daily fatalities F t evolve as:

We fix the incubation time period distribution as for the estimation of the exposure times (21) . We set the IFR as 0.83% according to the nationwide seroprevalence study (29) . The average time from symptom onset to decease is fixed as 11 days according to the information by the public health authorities (30) . Since no information is available on the type of distribution 3 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint which fits best the data from symptom onset to decease, we make use of a gamma distribution (5, 31) . We fix the shape factor (0.63) to match the interquartile range of the empirical distribution of 7 and 17 (30) . In the Supplementary Text we perform a sensitivity analysis by considering a log-normal distribution.

Finally, in order to fit the model to the data, we make use of a MCMC approach. To be more precise, we use the DREAMzs sampler (32, 33) from the BayesianTools package (34) . We use all flat priors, namely I0 ∈ [1, 1000], R1 ∈ [3, 8] , R2 ∈ [0, 1] and BP between March 1st and 15. The number of iterations is fixed to 3 · 10 6 with three chains. We fix the snooker update probability as 0.5. We use 20% of the iterations for adaption and discard the first 30% of the iterations. The posteriors and their correlation are shown in Fig. S2 , whereas the actual fit is presented in Fig. 3C .

To analyze what would have happened if the epidemic response had occurred earlier or later, we make use of the previous fit of the piece-wise linear R t . By shifting the breaking points, we can analyze straightforwardly the impact of an earlier or later response since we decoupled the containment from the effect of natural immunity. This approach is strongly inspired by Althaus et al. (11) . We make all the simulations sampling 5000 parameter sets from the calibration's posterior to bootstrap the credible intervals.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint

Heterogeneity in the epidemic response between CCAA In the main text, Fig. 2B shows the precise day on which R t was first below one in the CCAA. Here, we describe in more detail how this day correlates with the measures implemented in the CCAA. The Community of Madrid was the first CCAA to push R t number below one. Madrid was the strongest hit CCAA early on and implemented first containment measures on March 9 before all other CCAA. Similarly, other CCAA, which first pushed R t below one, were among the first to introduce containment measures such as Asturias (35), Valencian Community (36) or Basque Country (37) . Interestingly, the last CCAA to push R t below one, are in geographical proximity to Madrid such as Aragon, Castile-La Mancha, Castile and León, and Extremadura. This is consistent with a study by Mazzoli et al. that indicates how the local incidence strongly correlated with mobility in the early phase from and to Madrid, due to multi-seeding (38) . Furthermore, Castile-La Mancha, Castile and León, and Extremadura have the highest percentages in Spain of their students attending university in the Community of Madrid; 27%, 14%, and 12%, respectively (39) . Even though it is speculative, it seems reasonable that, after the closure of the universities, the outflux of students from Madrid to their hometowns added various seeds with mostly mild symptoms, accelerating the evolution of the epidemic in these regions.

The SARS-CoV-2 epidemic led to a substantial reduction in mobility. Figure S4 presents this data from the Google mobility report. We note a strong decrease of all types of mobility before, and shortly after the lockdown. Afterwards, the mobility stays very stable, but an additional decrease can be seen coinciding with the second lockdown as all-non essential mobility was closed. This additional reduction of mobility during the closure of all non-essential is better illustrated in Fig. S5 , which shows the average mobility during the three subsequent lockdowns. We also note that the third lockdown shows a slightly stronger mobility reduction compared to the first one despite the same restrictions being applied.

Comparing the mobility reduction with the decrease in R t , we find an elevated crosscorrelation. Figure S6 presents the cross-correlation between the mobility types and R t for the CCAA. A detailed distribution for the maximal cross-correlation is shown in Fig. S7 . Additionally, with very little exceptions, the maximum cross-correlation is reached without any lag between the two time series. The exception here is the mobility related to grocery shops and pharmacies. This negative lag is due to the strong increase in visits of grocery stores and pharmacies due to fear of supply shortages as shown in Fig. S4 .

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint

In the main text we separated R t into three linear segments: the "free" spreading phase before restrictions, the initiation of the epidemic response, and the lockdown. Here, we present a more detailed analysis regarding the epidemic response and the lockdown. Therefore, we cut the time series of R t at the turning point as the decrease is initiated in the three segment analysis. Subsequently, we analyze the different segments observed until the end of lockdown. We choose the number of segments according to the Bayesian information criterion (40) and corroborated with the Akaike information criterion (41) by iterating from 1 to 15 segments.

The identified segments are shown in Fig. S8 (bottom) together with the evolution of mobility according to the COVID-19 Community Mobility report from Google (42) (top). Interestingly, between the first and second breakpoint, we observe a substantial reduction of R t before mobility starts to decrease. The second breakpoint then coincides with the initiation of mobility reduction. The third breakpoint matches the implementation of the first lockdown on March 15. At breakpoint 4, R t starts a plateau together with the mobility level. The breakpoint 6 then corresponds to the shutdown of all non-essential activity, the second lockdown. This also coincides with a further reduction in overall mobility as shown in Fig. S5 .

The lift of the second lockdown and the associated mobility increase, seem not to have caused any increase in R t . It is probable that the good weather conditions, as well as an increased individual awareness, raised by the elevated death toll and reflected in the reduced mobility (see Fig. S5 ), have contributed to a continuous decrease in R t . Finally, at breakpoint 9 the increase in R t coincides with the lift of the lockdown and corresponding increase in mobility. Various identified breakpoints cannot be associated with any nation-wide introduced measures put in place by the national authorities or changes in mobility. They may stem from measures introduced in the different CCAA.

For adjusting the model to the daily fatalities we assumed a gamma distribution. The mean of the distribution was fixed to 11 days according to the information provided by the health authorities (30) . Furthermore, we adjusted the shape factor to match the interquartile range (7 and 17 days) since this is, besides the mean, the only information we have at our disposal of the empirical distribution. Even though gamma distributions have been shown to describe particularly well the time between symptom onset and disease (31), we performed the same analysis making use of a log-normal distribution. We assumed the same priors as in the case of the gamma distribution.

Adjusting the log-normal distribution to the interquartile range, we found a logarithmic standard deviation of 0.54. The log-normal distribution is compared to the gamma distribution in CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ;  (CI: 1.75 − 2.11) days and 12.1 (CI: 11.7 − 12.5) days, respectively (12) . The results are not substantially altered with respect to a gamma distribution. The smaller R t before the epidemic leads to a higher number of initially infected individuals. Furthermore, the shape of the lognormal distribution causes R t to decrease a day earlier. The posteriors and their correlations are shown in Fig. S10 . Figure S11 shows the obtained fit between the model and the data. Figure S12 shows the inferred form of R t , while in Fig. S13 CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ; Figure S1 : First breakpoint in R t and workplace related mobility for Spain and all CCAA. Points indicate the mean, whereas horizontal bars indicate 95% credible intervals. We explicitly separated the time series into three segments, and identified them with the segmented package (26) . Similarly as in the main text, we assumed the first segment of the time series was constant, and the others linear. To discard initial fluctuations in R t due to a small incidence, we did not consider values before March 1. Additionally, we cut the series on March 28 to exclude the additional decrease in mobility due to lockdown 2. Overall, R t decreases substantially earlier than mobility. 8 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021. ; Figure S2 : Posterior distribution of the model parameters and their correlations, assuming the time between symptom onset to death follows a gamma distribution. The range of the flat priors was larger than the posteriors in all cases. For the variable BP, the integer values refer to the dates in March 2020.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021 Google trends search frequency 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021 Figure S4 : Visualization of the mobility data from the COVID-19 Community Mobility report from Google (42) . Dots indicate the data points, while solid lines represent a centered, 7 day rolling average. The shaded areas in light gray indicate the lockdowns 1 and 3. In dark gray, we indicate the lockdown 2.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 18, 2021 Figure S5 : Average mobility according to the COVID-19 Community Mobility report from Google (42) during the three lockdown stages in Spain. Mobility is substantially lower during the second lockdown since all non-essential economic activity was shutdown. Additionally, mobility was lower during lockdown 3 when compared to lockdown 1, even though restrictions were the same. This indicates an increased awareness among the population, probably due to the huge death toll to this point.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Figure S6 : Lag and value of maximal cross-correlation between different mobility types, according to the COVID-19 Community Mobility report from Google (42) and R t . Negative lag indicates that R t precedes mobility.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) Figure S7 : Maximal cross-correlation among CCAA between different mobility types, according to the COVID-19 Community Mobility report from Google (42) and R t . For the Residential mobility category we take the absolute value to facilitate visualization.

14 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) Figure S9 : Gamma and log-normal distributions used for the time between symptom onset and decease. The average of both distributions is 11 days. The gamma distribution has a shape factor of 0.63, whereas the logarithmic standard deviation of the log-normal distribution is 0.54.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint Figure S10 : Posterior distribution for the model parameters and their correlations, assuming the time between symptom onset to death follows a log-normal distribution. The range of the flat priors was larger than the posteriors in all cases. For the variable BP, the integer values refer to the dates in March 2020.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) Figure S11 : The daily fatalities in Spain from March 8 to May 8 compared to the daily fatalities for the sensitivity analysis. In the sensitivity analysis, we used a log-normal distribution describing the time between from symptom onset to decease. The shaded area represents 95% credible intervals.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) Reproduction Number R(t) estimated from the data R(t) piece-wise Figure S12 : R t estimated from the reconstructed exposure times vs. the model inferred as part of the sensitivity analysis. This is equivalent to Fig. 4A in the main text when using the log-normal distribution.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) Figure S13 : Ratio between the reconstructed exposure times and the model inferred infections as part of the sensitivity analysis. This is equivalent to Fig. 4B in the main text when using the log-normal distribution.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Shift in epidemic response (days) Number of fatalities Figure S14 : Counterfactual scenarios for the attack rate as well as the total number of fatalities. This is equivalent to Fig. 4C in the main text when using the log-normal distribution.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted February 18, 2021. ; https://doi.org/10.1101/2021.02.16.21251832 doi: medRxiv preprint

Prevalence of SARS-CoV-2 in Spain (ENE-COVID): a nationwide, population-based seroepidemiological study

Household transmission of SARS-CoV-2 and risk factors for susceptibility and infectivity in Wuhan: a retrospective observational study

Retrospective analysis of the Italian exit strategy from COVID-19 lockdown

Ranking the effectiveness of worldwide COVID-19 government interventions

Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe

Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions

Modeling the spatiotemporal epidemic spreading of COVID-19 and the impact of mobility and social distancing interventions

Time is of the essence: containment of the SARS-CoV-2 epidemic in Switzerland from

Estimating the generation interval for coronavirus disease (COVID-19) based on symptom onset data

Estimating regression models with unknown break-points

Interdependence and the cost of uncoordinated responses to COVID-19

Household transmission of SARS-CoV-2 and risk factors for susceptibility and infectivity in Wuhan: a retrospective observational study

Retrospective analysis of the italian exit strategy from COVID-19 lockdown

Ranking the effectiveness of worldwide COVID-19 government interventions

Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions

A method of non-parametric back-projection and its application to aids data

A new framework and software to estimate time-varying reproduction numbers during epidemics

Modeling the spatiotemporal epidemic spreading of COVID-19 and the impact of mobility and social distancing interventions

Time is of the essence: containment of the SARS-CoV-2 epidemic in Switzerland from

Estimating the generation interval for coronavirus disease (COVID-19) based on symptom onset data

Estimating regression models with unknown break-points

Interdependence and the cost of uncoordinated responses to covid-19

Non-parametric back-projection of incidence cases to exposure time

Monitoring count time series in R: Aberration detection in public health surveillance

Spatio-temporal analysis of epidemic phenomena using the R package surveillance

Reconstruction of the infection curve for SARS epidemic in Beijing, China using a back-projection method

Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside Hubei province, China: a descriptive and modelling study

EpiEstim: Estimate time varying reproduction numbers from epidemic curves

Estimating individual and household reproduction numbers in an emerging epidemic

Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures

segmented: an R package to fit regression models with broken-line relationships

Modeling COVID-19 dynamics in Illinois under nonpharmaceutical interventions

Infection fatality risk for SARS-CoV-2 in community dwelling population of spain: nationwide seroepidemiological study

Instituto de Salud Carlos III, COVID-19

Inference of COVID-19 epidemiological distributions from Brazilian hospital data: Inference of COVID-19 epidemiological distributions from Brazilian hospital data

Differential Evolution Markov Chain with snooker updater and fewer chains

Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling

BayesianTools: General-purpose MCMC and SMC samplers and tools for Bayesian statistics

El Principado decreta el cierre de todos los centros educativos asturianos ante la recomendación del Gobierno

Qué impacto económico tendrá la cancelación de las Fallas de València?

El gobierno vasco decreta el cierre de todos los colegios de Vitoria 15 días por el COVID-19

Effects of mobility and multi-seeding on the propagation of the COVID-19 in Spain

Datos y cifras del sistema universitario espanol

Akaike Information Criterion statistics

Madrid deja sin clases a 1,5 millones de alumnos y cancelará operaciones y citas para combatir el coronavirus