Risk of Secondary Infection Waves of COVID-19 in an Insular Region: The Case of the Balearic Islands, Spain

The Spanish government declared the lockdown on March 14th, 2020 to tackle the fast-spreading of COVID-19. As a consequence, the Balearic Islands remained almost fully isolated due to the closing of airports and ports, these isolation measures and the home-based confinement have led to a low prevalence of COVID-19 in this region. We propose a compartmental model for the spread of COVID-19 including five compartments (Susceptible, Exposed, Presymptomatic Infective, Diseased, and Recovered), and the mobility between municipalities. The model parameters are calibrated with the temporal series of confirmed cases provided by the Spanish Ministry of Health. After calibration, the proposed model captures the trend of the official confirmed cases before and after the lockdown. We show that the estimated number of cases depends strongly on the initial dates of the local outbreak onset and the number of imported cases before the lockdown. Our estimations indicate that the population has not reached the level of herd immunization necessary to prevent future outbreaks. While the low prevalence, in comparison to mainland Spain, has prevented the saturation of the health system, this low prevalence translates into low immunization rates, therefore facilitating the propagation of new outbreaks that could lead to secondary waves of COVID-19 in the region. These findings warn about scenarios regarding after-lockdown-policies and the risk of second outbreaks, emphasize the need for widespread testing, and could potentially be extrapolated to other insular and continental regions.


INTRODUCTION
The rapid propagation of the new COVID-19 pandemic requires timely responses, including the alignment of evidence generation by scientists and decision-making by policy stakeholders. As of the current date, several mathematical models have been developed to help policy-making in a wide arrange of interventions in various countries, encompassing from testing strategies to lockdown measures (1)(2)(3)(4)(5)(6). Although modeling pandemics is not without flaws, and its predicted scenarios cannot be uncritically adopted and therefore directly translated into policy (7), modeling can be a valuable support tool to guide policy when assessed in an integrated way.
Recent studies have dealt with the possibility of a second-wave of COVID-19 after the retirement of lockdown and confinement measures in China (1,2). Recently, the value of restrictive social distancing measures has been largely proved since the first outbreak of the pandemic in European countries such as Italy (8). The analysis of data from closed confinements such as sea cruises allows us to address some key questions regarding the risk of second waves in an environment without external perturbations (9,10). The study of the evolution of the pathogen in islands offers an opportunity to learn how the propagation occurs, and how the mobility restrictions are shaping the propagation in relatively isolated areas, either due to transport lockdowns implemented to contain COVID-19 dissemination or because of their geographical conditions.
The Balearic Islands archipelago is composed of four inhabited islands in the Mediterranean Sea, i.e., Majorca, Menorca, Ibiza, and Formentera, with a total population of 1,095,426, as per 2011 data (11). The main economic activity is tourism with principal connections to the UK and Germany. The first reported case in Spain was identified in the Canary Islands on January 31st, while in the Balearic Islands the first case (second in Spain) was confirmed on February 9th. He was a British citizen resident in Majorca who had been in contact with an infected person with SARS-CoV-2 during a stay in France from January 25-29. In Spain, the schools were closed on March 16th and the lockdown was implemented at the national scale from March 17th. As of April 11th, the number of infected cases in the Balearic Islands was 140 per 100,000 inhabitants (1,507 confirmed cases) to be compared with 345 in Spain (161,852 confirmed cases; data updated with values of April 11th, 2020) (12). The lockdown of the Balearic Islands includes the closing of airports and ports for passengers, rendering the archipelago a virtually closed system. In this regard, archipelagos are "living laboratories" suggesting insights about the ecology and evolution of infectious diseases and offering unique experimental testing protocols to reduce or eliminate the diseases not only in the islands but potentially across the world (13). Thus, the Balearic Islands present an opportunity to be used as a benchmark to explore how isolation and after-lockdown measures impact secondary COVID-19 waves.
COVID-19 has a particular structure in the timings of the disease that make it particularly dangerous in terms of a silent spreading potential. First, the incubation period, i.e., the time since infection to symptom onset, is relatively large around 5.2 days [95% confidence interval (CI), 4.1-7.0] (14). This itself is a driver of the predictability of the spatiotemporal patterns to expect from this disease (15). Furthermore, the latent period, i.e., the time from infection to the start of being infectious, does not align completely with the incubation period (4). Although latency periods of median 3.69 days (95% CI, 3.30-3.96) have been reported (16) infections can occur days 1 and 2 after exposure (4,17). This leaves a period of presymptomatic infectivity, that increases R 0 through silent spreading, as not even the carrier might be aware of its own infectivity (18). The relative effectiveness of different non-pharmaceutical interventions will depend critically on the relation of those times (incubation and latent period) (19). Other related periods that shape the dynamics of the outbreaks are the generation interval (time between infection of infectorinfectee pairs) (20) [example of mean values of generation intervals are 5.20 days for Singapore data, and 3.95 days for Tianjin, (21)] and the serial interval (the time between symptom onsets of an infector-infectee pair), which has also been used to estimate viral shedding dynamics for COVID-19 (example of serial interval values are characterized by a mean of 5.8 days [95% confidence interval (CI), 4.8-6.8 days] (4).
We aimed to study the dissemination of COVID-19 in an isolated system through a compartmental model that included, besides the susceptible (S), diseased (D) and recovered (R) compartments, an exposed (E) compartment, and a presymptomatic (I) infective compartment to account for the incubation period, as the times of transit between the latter two compartments are crucial for the modeling of COVID-19 (3,4). Due to population size, we can implement an individualbased model where we consider each inhabitant as an individual in the model. In particular, first, we compare the results of an individual-based model tailored for the Balearic Islands and identify the parameter values that best fit the data. Second, we explore the likelihood of a second-wave scenario as a function of the initial date of the first imported case and the number of imported cases before the lockdown.

Data
Population data for the 67 municipalities in the Balearic Islands were taken from the Instituto Nacional de Estadística (INE, Spain), which gathers all the census data (11). The census also provides the commuting flows for people that, according to the registry, are living in one municipality and work in another. This allows assigning a living location and working location to each individual. For small municipalities, these commuting fluxes are not included. We avoid the isolation of these municipalities (to be specific, Formentera with a population of 12,111 inhabitants, Escorca 280, Estellencs 389, Banyalbufar 605, and Deyá 755) considering commuting flows of 10 people toward each of the neighboring municipalities and Palma, the capital of Balearic Islands.
Data for the active infected and cumulative infected cases are obtained from the Ministry of Health (12). In particular, the official reports provide data on the cumulative number of infected, recovered, and deaths. The number of active infected cases is taken as the cumulative number of infected cases and subtracting deaths and recovered. Unfortunately, the values for recovered cases are only reported from March 22nd. Thus, for the fitting, we considered all the historical series for the cumulative number of cases while only values starting from March 22nd for the number of active cases.

Model
The relatively small population size of the Balearic Islands allows us to develop an individual-based model. Each individual is placed in one of 67 municipalities according to the census (Supplementary Figure 1). The mobility between municipalities is considered with commuting data from the 2011 census provided by the INE (11). For each simulation day, we consider two steps, one where each individual is located in its residence municipality, and a second step where each individual is placed in the working place. At each step, individuals can interact with any of the individuals placed at the same location (Figure 1). The number of individuals with residence in location i and commuting with location j is denoted by N ij . Thus, the population at location i is given by N i = j N ij , and the population at location i during working time is given by N i ′ = j N ji . The location assignment is done sequentially starting from N 1,1 , then N 1,2 , . . . until N 67,67 . This assignment is done initially and such positions remain unchanged during the time evolution of the model.
The states of the individuals correspond to a SEIDR model: S, susceptible; E, exposed, corresponding to the latent period; I, infectious, corresponding to the presymptomatic infective period; D, diseased, corresponding to be infective with or without symptoms; and R, recovered. The transitions between these states are as follows, S becomes E in contact with an infected individual (I or D) with probability β. After T lat (latent period) days, E becomes I; after T inf (presymptomatic infective period) days I becomes D, and after T dis (disease), D becomes R ( Figure 1A).
The values of T lat , T inf , and T dis were obtained from the time evolution of the number of active infected and cumulative infected cases in the Balearic Islands. The lockdown was imposed in Spain on March 16th and the effect of the mobility restrictions can be identified on March 22nd. The 6 days in this period are reflected in the incubation period, T lat + T inf ( Supplementary Figures 2A,C), and agrees with recent estimations (22). Finally, from the data on the cumulative number of infected cases, the change in slope is observed on April 2nd, that is, T dis = 12 days (Supplementary Figures 2B,C).
To implement the mobility restrictions, we observe from the data that the cumulative number of infected cases shows a bending every 7 days approximately, which is in accordance with the beginning of the lockdown, and the restriction imposed on March 15th, and later corrected on March 22nd and March 29th (Supplementary Figure 2D). Thus, the model has the freedom to adjust the infection probability every week after March 15th.
For a single day, the modeling proceeds as follows (Figure 1),

First, it considers the population in their residence location,
for each municipality pairs of individuals in the same municipality are selected, say i and j. Then, i updates his/her state according to the dynamic rules. For each municipality p, N p pairs are chosen randomly where N p is the population size of the municipality p. 2. Second, we consider the individuals distributed in the municipalities of work. For each municipality p' , we chose N p ' pairs of random individuals working in the same municipality p'. 3. Resume from 1.
Thus, on average, in a day, each individual is updated twice.
For calibration, the model is run exploring all the parameters: T lat + T inf = 6 and T dis = 12; and β is explored in the range [0, 1] in the following periods: β 1 from the origin of the infection on February 9th to March 15th, β 2 from March 16th to March 22nd, β 3 from March 23th to March 29th, β 4 from March 30th to April 5th, β 5 from April 7th to April 11th The total number of infected cases depends on the date of the first infection. Models assume that the beginning of the outbreak is typically 30 days before the day when 10 infections are recorded (23). In the case of the Balearic Islands, on March 8th, 11 confirmed cases were reported. The first case reported in the Balearic Islands corresponds to an imported infection notified on February 9th. Consequently, the beginning of the outbreak was set on February 7th, 2 days before the first infected case was identified. Thus, we explore the date of the beginning of the infection between Jan 28th and Feb 7th.

Model Validation
The results of the model are validated with the official number of active infected and the cumulative number of infected cases between March 15th and April 11th. As the official values do not take into account the non-tested asymptomatic and the diseased not consulting to the healthcare systems, we assume that the reported values are a proportion of the values obtained from the model. Then, to validate the model parameters we minimize χ 2 , where α is a scale factor, that is, the ratio between estimated and confirmed cases, Y i is the value obtained from the model in day i, and y i is the official value in day i. Due to the initial exponential growth of the epidemics, we calculate χ 2 for the logarithm of the cases: The minimization of χ 2 leads to an optimal scale factor log(α * ) = 1/n (log(y i / Y i )), where n is the number of observation days. For this value of α * , we finally calculate the optimal values. Our assumption implies that the scale factor should be similar to both the active and cumulative infected cases.
For each set of parameters, we report the χ 2 of the model values of the number of active infected cases with respect to the official values, the correction fraction α active , and the χ 2 acc of the model values of the number of cumulative infected cases with respect to the official vales and the correction fraction α acc . For each set of parameters, the best fit is considered as the one leading to the minimum χ 2 . Once the fitting values are determined, we calculate χ 2 acc and α acc . For each set of periods (T lat , T inf , T dis ), we explore the infection probabilities that minimize χ 2 of the number of active infected cases. The value of χ 2 and scale factors α * of the best fits are shown in Supplementary

Quarantining Mechanism
In order to explore the effect of quarantines on the impact of the second wave we implemented a stochastic quarantine mechanism. When an individual reaches the symptom onset, i.e., when they reach the D state, with probability q the individual is isolated from further interaction with other individuals. The probability q can be understood as the fraction of cases that will be quarantined or isolated.

Herd Immunity Assumptions
An approximation to the herd immunity threshold is given by 1-1/R 0 (24), which for COVID-19 is expected to be between 29 and 74%, taking R 0 between 1.4 and 3.9 (14,25). To explore whether the number of cumulative infections reach the herd immunity threshold and therefore avoidance of potential second waves is to be expected, we run the model for the same parameters leading to the best fit (24). After the system has relaxed to zero infection, we select a random susceptible from the populations and infected her. As we are interested in whether the epidemics will spread again, we use the transmission rate obtained at the beginning of the epidemics in the Balearic Islands, that is, before any restriction on mobility had been applied. We can expect that once the mobility restrictions have been removed, the transmission will be reduced in comparison to the initial values, especially due to an improvement in the hygiene of the population. This will affect how fast the COVID-19 will spread and the intensity of the wave. If the estimated number of infected cases is lower than the threshold for herd immunization, we assume that SARS-CoV-2 will spread. In the Supplementary Materials, we show how the data can be collapsed using a proper combination of the initial number of infected cases and the time of the first infection.

Number of Active Infected Cases
The best fit of the model to the confirmed cases, allows us to extract the transmission probabilities and also the scaling factor that captures the ratio between the estimated and the confirmed cases. For the scenario where the initial date was on Feb 7th, and a latent period of 2 days (T lat = 2), an infective period of 4 days (T inf = 4), and disease period of 12 days (T dis = 12), the values of the infection probability leading to the best fit are β 1 = 0.24, β 2 = 0.12, β 3 = 0.016, and β 4 = 0.036 (Figure 2). This translates into an initial basic reproductive number R 0 = 3.84. The fitting of the data also informs us that the correction factor is 0.054, that is, that the confirmed cases are 5.4% of the model estimates. At the same time, we obtain that the percentage of added recovered cases and fatalities according to the official sources the model estimate are 5.3% of the confirmed cases. For the other values of the latent, infective, and disease periods, we obtain similar accuracy, given by χ 2 and similar scaling factors (Supplementary Table 1). The scaling factor, which gives the fraction of the model estimates, that corresponds to the confirmed case, increases to 10% in the case of T lat = 5, T inf = 1, which is the case with less infected individuals in the model.
The introduction of a single imported infected case after the first wave has expired produces a secondary wave that strongly depends on the first one (Figure 3). The intensity and duration of the second wave depend on specific values capturing the conditions applicable when newly infected cases appear, e.g., the transmission probability, which depends on the habits of the population, hygiene, and social distancing, and mobility restriction. Qualitatively similar results were obtained for the other set of values of the characteristic periods. The peak of the second wave is very sensitive to the date of the first exposure. If it happened on January 28th, the intensity of the second peak is less pronounced and similar to the one for the first peak, in contrast to the case of a more recent exposure, when the second peak can be more than one order of magnitude larger the first peak.

Effect of Quarantines on the Second Wave
We show the results of applying different levels of the stochastic quarantine mechanism in Figure 3. The greater the percentage of quarantined cases, the smaller is the peak of the second wave, and the longer is it in duration. We also observe that, depending on the fraction of cases that occurred during the first wave, the effect of quarantining varies. So for the case when the first infection occurred on Feb 7th (Figure 3A), the fraction of quarantined cases that we need to reduce the peak of the second wave by one order of magnitude is around 60%, while for a stronger first wave ( Figure 3B) the reduction of one order of magnitude of the second wave is obtained with 40% quarantine. We simulated different degrees of quarantines to explore its effect on the second wave. The greater the percentage of infected cases that are quarantined at symptom onset, the smaller the incidence of the second wave. The shaded area represents the 95% CI (for clarity, only shown for the case of no quarantine). Confidence intervals calculated using Cox's method (26). (A) First infected case is on February 7th, the number of infected in the first wave has not reached herd immunization and a second wave is triggered by a single infected case. The intensity of the second wave is one order of magnitude larger than the first. (B) First infected case is on January 28th, the number of infected cases in the first wave has reached a larger fraction of the population and the intensity of the second wave is, in the scenario, smaller than the first. Average over 100 realizations.

Herd Immunization Estimates
Assuming recovered individuals get immunity, to estimate whether the Balearic Islands have reached herd immunization, we explored the estimated number of infected cases under two immunization scenarios based on the date of the first infection ( Figure 4A) and the number of imported cases before air and maritime transport lockdown ( Figure 4B). Firstly, we analyzed how the estimation of infected individuals is sensitive to the date of the first infection. We explored the time range of the first infection from January 28th (which corresponds to the stay in France before returning to Majorca) to Feb 7th (which corresponds to 30 days before the 10th confirmed case). Secondly, we explored the estimates under the assumption that more than one imported case could have gone unnoticed into the Balearic Islands before the closing of the airports. Depending on these two parameters, the range of immunization spans from 3% (for one initial infected on February 7th) to 64% (for 20 initial infected cases on January 28th). With the assumption of immunity after recovery, the achievement of herd immunization in the population is very sensitive to the date of the first infection and the number of imported cases before air and maritime transport lockdown. The interpretation of herd immunization indicates that if infected individuals become immune, then a percentage r of herd immunization prevents the spreading of reproductive numbers smaller than 1/(1-r). Assuming that the first case was exposed to the infection during his stay in France in the last days of January, the percentage of the population that was infected can be as high as 50%, which could prevent a high second peak, for values of the basic reproductive number below 2. Conversely, if the first case was infected 30 days before 11 confirmed cases were reported in the Balearic Islands, the percentage of infected individuals could lower to <10%, therefore falling far from potential herd immunity (only for values of the basic reproductive number below 1.1). The relation between the number of initial infected cases, the date of the first infection, the number of cases, and the number of confirmed cases is further explored in the Supplementary Materials.

CONCLUSIONS
Our study shows that a model including five compartments together with information on mobility between municipalities can be used to capture the spread of the epidemics in a closed community. The validation of the model with the official data allowed us to obtain the parameters that best fitted the data. Once the model was validated, we extracted an estimation of the number of the total infected in the Balearic Islands that indicates, assuming immunization after recovery, that these figures would reach the herd immunization threshold depending critically on the date of the first infection and the initial number of seeds, being herd immunization achievement more likely for an initial date before January 31st and number of initial infected above 10. Our exploration of the forecasted scenario of a newly infected individual entering the community after the lockdown confirmed that the number of potential cases widely varies according to the initial date of infection, which correlates with the percentages of immunity. Although we cannot determine with precision the start of the infection in the Balearic Islands, the model suggests that the Balearic Islands population is below the herd immunization threshold and thus, also susceptible to new outbreaks depending on how immunity is acquired and how the mobility restrictions are further implemented. In particular mobility and transmission probability, which depends on the general use of masks and hygiene protocols by the population, might alter the attack rate.
Focusing on second waves in insulated areas during the COVID-19 pandemic is of great value to analyze the spreading and containment of infectious diseases, where the lockdown of islands constitutes a paradigmatic scenario, with the potential to be applied to continental regions (13). For example, the risk of COVID-19 import to the Pacific Islands had been assessed and analyzed in the early stages (27)(28)(29). The use of modeling tools is a complement to field studies that can be used to anticipate the progress of a pandemic and thus help health authorities to make decisions. In the case of the Balearic Islands, there are two foremost advantages in terms of model precision. First, since the incidence during the first peak was relatively low and hospital capacity including ICU beds was not overpassed, the forecasted scenario of a second wave presenting with more intensity is more feasible than in other areas. Second, the relatively small size of the Balearic Islands and the organization of health and epidemiological surveillance systems make the official accounts of reported cases more reliable than in other areas were due to low rates of testing, overloaded hospitals, and lack of centralized data collection hampered the initial estimates.
The implications of the forecasted potential second wave yielded by our model for an insular territory can be useful also for other areas that are either naturally geographically isolated or closed to external perturbations due to strict lockdowns. According to our results, the date of the first infection and the import of cases while the airports and ports were open appear to be key to assess the likelihood and intensity of future waves and outbreaks. Knowing the approximate date of infection of the first reported case in an outbreak proved critical to estimate the current and foreseen number of cases. Whether a second wave occurs and the intensity of the peak strongly depended on the date of the first infection, as the number of infected cases grows exponentially, but also on the number of imported cases, which contribute additively to the number of cases, and also on the real herd immunization. Our estimates rely on calculations assuming conditions far from the behavior of the population, and on the habits, for example, regarding hygiene, the use of masks, and social distancing, of the population after the lockdown is relaxed. We further show that the effect of quarantining measures strongly depends on the level of immunization reached by the population during the first wave.
Our model is an individual-based model for which, due to the population size, we identify each inhabitant with an individual in the simulation model. This approach is different from other models considering pan-mixing and ordinary differential equations (8), and forecasts based on iteration methods (24,30). Other approaches implement recurrent mobility (3,31,32), which selects the individuals that perform the commuting randomly at each iteration step, thus increasing the mixing in the complete population. Our approach assigns a residence and a working municipality as initial condition and these locations remain unchanged during the time evolution for the model. Our implementation assumes that the same person commutes between two locations and thus it has to be fixed initially in the model. A random selection at each day will increase the number of effective connections, which could be compensated by a reduction in the transmission probability. We believe this approach is more comprehensive and better captures the reality of commuting under home confinement conditions, which essentially limit mobility from households to workplaces for those individuals that cannot work remotely or are not exempted from any work under the regulations of each country, while the rest of the population are not supposed to move from the vicinity of their households and even then only for justified reasons such as basic food supply. We use here a stochastic approach similar to other works (22,33) which lets us compute confidence intervals even for single combinations of the parameters instead of deterministic ordinary differential equations (8,34,35) or discrete-time dynamic equations (3). We also use a fixed time for the transit through the E, I, and D compartments, T lat , T inf , and T dis , respectively. We believe this approach is more realistic than an approach based on rates, where individuals transit the compartments at a given rate, giving rise to exponentially distributed times of transit through the compartments. In this case, infected cases will have the opportunity to be infectious immediately, or to transit the I compartment also immediately, bringing the start of secondary infections closer to the time they were infected for many individuals. A similar effect happens with the length of the disease (time in the D compartment), having then individuals that immediately recover.
The model also has several limitations. First, as it is constructed for fitting the global numbers of infected patients, it is missing finer structure, needed for the evaluation of risks of subpopulations that are differently exposed to the virus or have different outcomes, such as the population of elderly people or health workers, and the effect of city size (36,37). Second, for COVID-19 there is evidence of three main transmission channels, namely direct contact with an infected individual with symptoms (14), contacts with an asymptomatic individual (38,39), and environmental transmission (40). The present model takes into account the first two modes of transmission, but not the environmental one explicitly, although probably the fitting has assigned part of this transmission to the processes included in our model. Therefore, there is not a direct way of measuring the effect of interventions to reduce environmental transmission. Third, the model also considers asymptomatic and symptomatic individuals to be infectious in the same way, although the viral shedding in asymptomatic individuals is indeed lower (5). This can have an impact on the number of infected individuals and deserves future research. Fourth, the model assumes that the mobility restrictions are applied in the same way to all of the agents in the system and thus is lacking the fact that symptomatic infected individuals will modify their mobility drastically, either if they are quarantined at home or admitted to a hospital. We are therefore overestimating mobility, but this is probably passed to the infectivity in the fitting procedure. Fifth, the model also takes fixed times to transit through the E, I, and D compartments (T lat , T inf , and T dis , respectively), which is artificial. More refined models would take these transit times from specified distributions matching the parameters of the disease (33,35). While this will render the model more realistic, we believe that fixing the times is a good compromise between using rates for transiting the compartments and implementing distributions for those times, as it already captures the delays induced by these particular timings of the infection. Finally, the model also assumes that individuals are granted immunity to the virus, at least for the timescales explored here.
In conclusion, the risk of secondary infection waves should be comprehensively and cautiously addressed before removing confinement measures. Our study provides several relevant findings that could be useful to support policy design at avoiding second waves once measures to return to the societal usual activities start to be applied. First, the isolation of asymptomatic individuals that tested positive for SARS-CoV-2 and close contacts to infected individuals during the prior 2 weeks might reduce the number of new infections after the establishment of the usual activity by preventing dissemination from asymptomatic carriers during the incubation period. This requires proper testing strategies tailored according to the estimated prevalence of infection, population density, the openness of the community, and other relevant factors. Second, contact tracing measures are crucial, and digital tools might enhance the identification of high-risk individuals to be tested or preemptively isolated (41). Yet, data privacy and other relevant ethical considerations should be carefully balanced when designing contact tracing in the community. Third, progressive return to the normal activity instead of an abrupt change will facilitate the monitoring of new cases and may avoid a sharp growth of the number of infected individuals, which is expected when herd immunity has not been reached. Experience from other archipelagos illustrates the potential of pursuing an elimination strategy (including a full lockdown) together with quarantine of travelers from abroad (42,43). Further modeling studies on second-waves of COVID-19 are warranted to strengthen the knowledge on the best theoretical assumptions and data to be used to increase forecasting precision. In addition, these models should be validated through real-world data as these are collected during and after the pandemic.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. We obtained the data from official sources of the Spanish government, including the Ministry of Health and the Instituto de Salud Carlos III for the account of official cases (https://covid19.isciii.es/resources/serie_historica_acumulados. csv) and the Instituto Nacional de Estadística for the mobility data (https://www.ine.es). Computing codes are available at https://github.com/juanfernandezgracia/Balearic-epi.

AUTHOR CONTRIBUTIONS
VE and JF-G designed the work. VE performed the analysis. VE, JF-G, JR, JP, and CM prepared the figures, tables, wrote the first draft, and provided final approval to the manuscript.
All authors contributed to the article and approved the submitted version.
FUNDING VE and JF-G acknowledge funding from the Ministry of Science and Innovation (Spain) and FEDER through project SPASIMM [FIS2016-80067-P (AEI/FEDER, UE)]. JF-G acknowledges funding from the Vicerrectorado de Investigación e Internacionalización of the University of the Balearic Islands and Campus de Excelencia Internacional CEI15-09 (Ministerio de Educación, Cultura y Deporte, Spain) through its talent attraction program.