We generated a synthetic community layer connecting randomly pairs of individuals living in the same district, according to the contact matrices for Spain in Prem et al.. There were contact matrices available for home, work, school, and other locations, and we chose the latter. This dataset reported the probability of connecting pairs of individuals according to their ages, in age groups of 5 years up to 75 years old. For individuals older than 75 years old, we extrapolated the data of the oldest available group.
We used the COVASIM software for modeling the spread of COVID-19. COVASIM is an open-source Python-based agent-based modeling tool. COVASIM considers a susceptible-exposed-infected-recovered or dead (SEIRD) epidemic model that includes disease parameters informed by the medical literature. The infected compartment is divided into asymptomatic and symptomatic infectious individuals, with the latter including presymptomatic, mild, severe, and critical stages. The three symptomatic stages can evolve to the recovered state, while the critical state can alternatively lead to the death of the individual (Supplementary Figure S1). The probabilities of developing symptoms, severe symptoms, a critical case, and from it the death of the individual are specified by age groups, arranged in 10-year-long age cohorts. This software has been used for studying different scenarios, for example, assessing the test-trace-quarantine strategy or quantifying the risk of outbreaks after international border opening. We modified COVASIM to include the specific details of our synthetic cities. Specifically, we included the age, sex, and contacts of the individuals in each of the considered cities. We ran independent simulations where each simulation chose one randomly infected seed as the first infected individual. Then, we kept the endemic realizations, defined as those leading to a finite number of deaths, which we set higher than 10 for the first wave. For the second wave, we also requested that there were more than 10 death observations in the last 10 days of the realization. Apart from the internal parameters of COVASIM, we considered independently the infection rate and the date of arrival of the first infected individual. We calibrated both parameters for each city analyzing the official time series of deaths (see section 2.2.2), which were more reliable than the number of cases that could suffer from high underdetection rates.
We obtained the temporal series of the number of deaths and the number of confirmed cases from the Spanish Ministry of Health at the province level, with daily resolution. Then, we multiplied these values by the fraction of the province's population living in the city. We averaged the rescaled data over a moving window of 7 days (the specific day +- 3 days) to smooth the fluctuations.
To illustrate the potential of our approach, after calibrating the first wave of COVID-19 in our model with synthetic digital cities, we implemented a case study of the second wave focusing on the city of Zaragoza. This second wave occurred between July and December 2020.
Testing and tracing. Positive tested individuals and their close contacts were isolated for 14 days until 30th September, and from then on for 10 days.
Restrictions on restaurants, cafes, and nightlife. Starting on 5th August, lifted on 4th September, and re-started on 19th October.
Opening of the schools with reduced group sizes and safety measures. Different levels started progressively, from 7th September to 17th September.
Opening of the universities with reduced group sizes. The university was opened on 14th September.
Interventions impacting the community layer. We considered the interventions related to the State of Alarm and those related to the capacity and schedules of restaurants and bars.
The non-pharmaceutical interventions that were introduced in this city to mitigate the spread of COVID-19 were the following: 
In addition to these policies, we also took into account the annual leave period of workers and considered a reduced number of interactions in the work layer starting on 15th July until the middle of September, with the maximum reduction on 1st August. In terms of the model, this implies a higher amount of time (e.g., more weight of this layer) associated with the interactions in the community layer.
We used the infection rate (beta) and date of arrival of the infected seed obtained in the calibration of the first wave and ran the model from the estimated arrival of the first case (see section 3.2) to 1st December 2020. Since we ran the simulations from the beginning of the first wave, we lifted progressively the restrictions that were active from the 14th of March in two steps: on the 1st and 20th of June, following the progressive lift of the restrictions that actually happened.
We kept 25% of the contacts in the school and university layers to simulate the small group's policy, and reduced beta by 50% in the school layer, considering the strict protocol to avoid contagion at schools.
After calibrating the model of the second wave with these non-pharmaceutical interventions and considering the results of the simulation, we introduced different counterfactual scenarios in which we quantified the effect of each of the non-pharmaceutical interventions adopted. More specifically, we computed the number of deaths and the disease prevalence by performing simulations with the same epidemiological parameters but switching on and off alternative interventions. Finally, we computed the relative change in the relevant quantity X (X = deaths or prevalence) as r = (Xcounterfactual - Xsimulated-2nd-wave)/Xsimulated-2nd-wave. Therefore, the absolute change can be obtained as Xcounterfactual = (1 + r)   Xsimulated-2nd-wave.
No testing and no contact tracing. The testing intervention was removed. Hence, as contact tracing depends on the results of the testing process, contact tracing was automatically removed.
No contact tracing. To analyze the impact of the contact tracing strategy and decouple it from the testing process, we kept the testing intervention and its parameters but removed contact tracing interventions.
Opening 100% university. We considered the opening of the university layer with 100% of the contacts, instead of the 25% contacts estimated through the small group's intervention.
Not opening university. We simulated a scenario where the university layer remained closed.
Opening all schools together (x2). The school opening was done following a staggered strategy, such that each level started on a different date. We simulated scenarios where all the levels started on the same date, either on the date of the earliest opening (7th September) or the latest opening (17th September).
100% beta in schools, whole groups. Schools were one of the sectors where strong protocols were introduced, reducing considerably the infection rate and also the group size. We simulated the absence of these protocols, keeping the same infection rate as in the rest of the layers, and considering this layer with whole groups, that is, 100% of the contacts.
Not opening schools. We quantified the changes in the outcome of the second wave if the schools had not been open.
No interventions in October. We observed that the interventions in October were key to controlling the second wave. Thus, we removed these interventions and computed this counterfactual, keeping the same final date, such that the result was comparable with the rest of the counterfactuals. However, we assumed that removing these interventions would imply that the second wave continued growing on time. To characterize this growth (in terms of both time extent and outcome), we ran additional simulations for 15 and 30 days more and compared them with extrapolations of the original second wave simulation for the same period (without including additional measurements introduced on December 2020 or calibrating the observed data in that period).
We considered the following nine different counterfactuals: