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–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₀ 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 infector-infectee 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 pre-symptomatic (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 individual-based 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.

Materials and Methods

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.

FIGURE 1

Figure 1. Reactions in the model, interactions and movements. (A) Reactions in the epidemic model. (B) Map of the Balearic Islands showing all the municipalities in the islands. (C–F) Zoom to the region in the red box in (B). (C,E) Show the interactions during the day and night, respectively. (D,F) Show the mobility to the night location (home) and day location (work), respectively. Note that the mobility patterns are reversed from (D) to (E), i.e., the same agents always do the same commute.

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),

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: β₁ from the origin of the infection on February 9th to March 15th, β₂ from March 16th to March 22nd, β₃ from March 23th to March 29th, β₄ from March 30th to April 5th, β₅ 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 χ², χ²= Σ (α Y_i - y_i)², 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 χ² for the logarithm of the cases: χ² = Σ (log(α Y_i) - log(y_i))². The minimization of χ² 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 χ² of the model values of the number of active infected cases with respect to the official values, the correction fraction α_active, and the ${\chi }_{\text{acc}}^2$ 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 χ². Once the fitting values are determined, we calculate ${\chi }_{\text{acc}}^2$ and α_acc. For each set of periods (T_lat, T_inf, T_dis), we explore the infection probabilities that minimize χ² of the number of active infected cases. The value of χ² and scale factors α^* of the best fits are shown in Supplementary Table 1 and the estimated prevalence in Supplementary Table 2.

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₀ (24), which for COVID-19 is expected to be between 29 and 74%, taking R₀ 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.

Results

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 β₁ = 0.24, β₂ = 0.12, β₃ = 0.016, and β₄ = 0.036 (Figure 2). This translates into an initial basic reproductive number R₀ = 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 χ² 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.

FIGURE 2

Figure 2. Active and cumulative infected time series in the Balearic Islands. Time evolution of the average number of (A) active infected cases (black line) and (B) cumulative number of cases with the best-fitted parameters, and confidence interval CI 95% (gray area). Confidence intervals calculated using Cox's method (26). Red dots represent the data provided by the Ministry of Health of Spain, the solid lines depict the model with the fitted parameters T_lat = 2, T_inf = 4, T_dis = 12, β₁ = 0.24, β₂ = 0.12, β₃ = 0.016, β₄ = 0.036, β₅ = 0, the orange area represents the number of undetected cases when comparing model results with the data on detected cases.

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.

FIGURE 3

Figure 3. Secondary outbreak appearing after the home-based confinement is removed. The lockdown is removed and the parameter values are as during the week of March 16th−22nd. 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.

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.

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.

FIGURE 4

Figure 4. Fraction of infected individuals (in logarithmic scale) as a function of the date of the first infection. (A) Each line corresponds to the following parameter values(T_lat, T_inf)= (1,5) (black),(2,4) (red), (3,33) (orange), (4,2) (magenta), and (5,1) (blue). (B) Fraction of infected cases as a function of the first infection and the initial number of infected individuals. The radius of the symbol is proportional to the fraction of infected cases while the color indicates the probability that a realization of the model reaches at least 50% of infected cases in the population, P₅₀.

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–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.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This manuscript has been released as a pre-print at medrxiv (44). The silhouettes in Figure 1 were designed by brgfx/Freepik.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmed.2020.563455/full#supplementary-material

References

1. Leung K, Wu JT, Liu D, Leung GM. First-wave COVID-19 transmissibility and severity in China outside Hubei after control measures, and second-wave scenario planning: a modeling impact assessment. Lancet. (2020) 395:P1382–93. doi: 10.1016/S0140-6736(20)30746-7

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Zhang L, Shen M, Ma X, Su S, Gong W, Wang J, et al. What is required to prevent a second major outbreak of SARS-CoV-2 upon lifting the quarantine of Wuhan city, China. Innovation. (2020) 1:100006. doi: 10.1016/j.xinn.2020.04.006

CrossRef Full Text | Google Scholar

3. Arenas A, Cota W, Gomez-Gardeñes J, Gómez S, Granell C, Matamalas JT, et al. A mathematical model for the spatiotemporal epidemic spreading of COVID19. medRxiv. (2020). doi: 10.1101/2020.03.21.20040022

CrossRef Full Text | Google Scholar

4. He X, Lau EHY, Wu P, Deng X, Wang J, Hao X, et al. Temporal dynamics in viral shedding and transmissibility of COVID-19. Nat Med. (2020) 26:672–5. doi: 10.1038/s41591-020-0869-5

CrossRef Full Text | Google Scholar

5. Ferguson NM, Laydon D, Nedjati-Gilani G, Imai N, Ainslie KM, Baguelin M, et al. Impact of Non-Pharmaceutical Interventions (NPIs) to Reduce COVID19 Mortality and Healthcare Demand. London: Imperial College COVID-19 Response Team Report (2020). Available online at: https://www.imperial.ac.uk/mrc-global-infectious-disease-analysis/covid-19/report-9-impact-of-npis-on-covid-19/

PubMed Abstract | Google Scholar

6. Verity R, Okell LC, Dorigatti I, Winskill P, Whittaker C, Imai N, et al. Estimates of the severity of coronavirus disease 2019: a model-based analysis. Lancet Infect Dis. (2020) 20:669–77. doi: 10.1016/S1473-3099(20)30243-7

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Sridhar D, Majumder MS. Modelling the pandemic. BMJ. (2020) 369:m1567. doi: 10.1136/bmj.m1567

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Giordano G, Blanchini F, Bruno R, Colaneri P, Di Filippo A, Di Matteo A, et al. Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy. Nat Med. (2020) 26:855–60. doi: 10.1038/s41591-020-0883-7

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Mizumoto K, Kagaya K, Zarebski A, Chowell G. Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the Diamond Princess cruise ship, Yokohama, Japan, 2020. Euro Surveill. (2020) 25:2000180. doi: 10.2807/1560-7917.ES.2020.25.10.2000180

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zhang S, Diao M, Yu W, Pei L, Lin Z, Chen D. Estimation of the reproductive number of novel coronavirus (COVID-19) and the probable outbreak size on the diamond princess Cruise ship: a data-driven analysis. Int J Infect Dis. (2020) 93:201–4. doi: 10.1016/j.ijid.2020.02.033

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Censo de Población y Viviendas. Instituto Nacional de Estadistica (Spain). (2011). Available online at: https://www.ine.es (accessed March 30, 2020).

Google Scholar

12. Available, online at: https://covid19.isciii.es/resources/serie_historica_acumulados.csv (accessed April 12, 2020).

13. Cowley G, Da Silva ET, Nabicassa M, De Barros PDP, Blif MM, Bailey R, et al. Is trachoma on track for elimination by 2020? Monitoring and surveillance after mass drug administration with azithromycin for active trachoma in Guinea Bissau. BMJ Global Health. (2017) 2:A62. doi: 10.1136/bmjgh-2016-000260.166

CrossRef Full Text | Google Scholar

14. Li Q, Guan X, Wu P, Wang X, Zhou L, Tong Y, et al. Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia. N Engl J Med. (2020) 382:1199–207. doi: 10.1056/NEJMoa2001316

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Kahn R, Peak CM, Fernández-Gracia J, Hill A, Jambai A, Ganda L, et al. Incubation periods impact the spatial predictability of cholera and Ebola outbreaks in Sierra Leone. Proc Natl Acad Sci USA. (2020) 117:5067–73. doi: 10.1073/pnas.1913052117

CrossRef Full Text | Google Scholar

16. Li R, Pei S, Chen B, Song Y, Zhang T, Yang W, et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2). Science. (2020) 368:489–93. doi: 10.1126/science.abb3221

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Lauer SA, Grantz KH, Bi Q, Jones FK, Zheng Q, Meredith H, et al. The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Ann Intern Med. (2020) 172:577–82. doi: 10.7326/M20-0504

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Fraser C, Riley S, Anderson RM, Ferguson NM. Factors that make an infectious disease outbreak controllable. Proc Natl Acad Sci USA. (2004) 101:6146–51. doi: 10.1073/pnas.0307506101

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Peak CM, Childs LM, Grad YH, Buckee CO. Comparing nonpharmaceutical interventions for containing emerging epidemics. Proc Natl Acad Sci USA. (2017) 114:4023–8. doi: 10.1073/pnas.1616438114

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wallinga J, Lipsitch M. How generation intervals shape the relationship between growth rates and reproductive numbers. Proc R Soc. (2007) 274:599–604. doi: 10.1098/rspb.2006.3754

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Ganyani T, Kremer C, Chen D, Torneri A, Faes C, Wallinga J, et al. Estimating the generation interval for coronavirus disease (COVID-19) based on symptom onset data, March 2020. Euro Surveill. (2020) 25:2000257. doi: 10.2807/1560-7917.ES.2020.25.17.2000257

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bi Q, Wu Y, Mei S, Ye C, Zou X, Zhang Z, et al. Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in Shenzhen, China: a retrospective cohort study. Lancet Infect Dis. (2020) 20: 911–919. doi: 10.1016/S1473-3099(20)30287-5

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Fine PE. Herd immunity: history, theory, practice. Epidemiol Rev. (1993) 15:265–302. doi: 10.1093/oxfordjournals.epirev.a036121

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Flaxman S, Mishra S, Gandy A, Unwin HJT, Coupland H, Mellan TA, et al. Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe. Nature. (2020) 584:257–61.

PubMed Abstract | Google Scholar

25. Riou J, Althaus CL. Pattern of early human-to-human transmission of Wuhan 2019 novel coronavirus (2019-nCoV), December 2019 to January 2020. Euro Surveill. (2020) 25:2000058. doi: 10.2807/1560-7917.ES.2020.25.4.2000058

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zhou XH, Gao S. Confidence intervals for the log-normal mean. Stat Med. (1997) 16:783–90. doi: 10.1002/(SICI)1097-0258(19970415)16:7<783::AID-SIM488>3.0.CO;2-2

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Craig AT, Heywood AE, Hall J. Risk of COVID-19 importation to the Pacific islands through global air travel. Epidemiol Infect. (2020) 148:e71. doi: 10.1017/S0950268820000710

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Mei Y, Hu J. Preparedness is essential for Western Pacific islands during the COVID-19 pandemic. Disaster Med Public Health Prep. (2020) 16:1–5. doi: 10.1017/dmp.2020.102

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kerbaj J, Cazorla C, De Greslan T, Gourinat AC, Marot B. COVID-19: the new Caledonia experience. Clin Infect Dis. (2020) 71:2279–81. doi: 10.1093/cid/ciaa600

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Perc M, Gorišek Miksić N, Slavinec M, StoŽer A. Forecasting COVID-19. Front Phys. (2020) 8:127. doi: 10.3389/fphy.2020.00127

CrossRef Full Text | Google Scholar

31. Aleta A, Moreno Y. Evaluation of the potential incidence of COVID-19 and effectiveness of contention measures in Spain: a data-driven approach. BMC Med. (2020) 18:157. doi: 10.1186/s12916-020-01619-5

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Gómez-Gardenes J, Soriano-Panos D, Arenas A. Critical regimes driven by recurrent mobility patterns of reaction-diffusion processes in networks. Nat Phys. (2018) 14:391–5. doi: 10.1038/s41567-017-0022-7

CrossRef Full Text | Google Scholar

33. Kucharski AJ, Russell TW, Diamond C, Liu Y, Edmunds J, Funk S, et al. Centre for mathematical modelling of infectious diseases COVID-19 working group. Early dynamics of transmission and control of COVID-19: a mathematical modelling study. Lancet Infect Dis. (2020) 20:553–8. doi: 10.1101/2020.01.31.20019901

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Lin Q, Zhao S, Gao D, Lou Y, Yang S, Musa SS, et al. A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action. Int J Inf Dis. (2020) 93:211–6. doi: 10.1016/j.ijid.2020.02.058

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Wu JT, Leung K, Bushman M, Kishore N, Niehus R, de Salazar PM, et al. Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan, China. Nat Med. (2020) 26:506–10. doi: 10.1038/s41591-020-0822-7

CrossRef Full Text | Google Scholar

36. Ribeiro HV, Sunahara AS, Sutton J, Perc M, Hanley QS. City size and the spreading of COVID-19 in Brazil. PLoS ONE. (2020) 15:e0239699. doi: 10.1371/journal.pone.0239699

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Stier AJ, Berman MG, Bettencourt LMA. COVID-19 attack rate increases with city size. medrxiv. (2020). doi: 10.1101/2020.03.22.20041004

CrossRef Full Text | Google Scholar

38. Bai Y, Yao L, Wei T, Tian F, Jin DY, Chen L, et al. Presumed asymptomatic carrier transmission of COVID-19. JAMA. (2020) 323:1406–7. doi: 10.1001/jama.2020.2565

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Tong ZD, Tang A, Li KF, Li P, Wang HL, Yi JP, et al. Potential presymptomatic transmission of SARS-CoV-2, Zhejiang Province, China, 2020. Emerg Infect Dis. (2020) 26:1052–54. doi: 10.3201/eid2605.200198

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Ong SWX, Tan YK, Chia PY, Lee TH, Ng OT, Wong MSY, et al. Air, surface environmental, and personal protective equipment contamination by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) from a symptomatic patient. JAMA. (2020) 323:1610–2. doi: 10.1001/jama.2020.3227

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ferretti L, Wymant C, Kendall M, Zhao L, Nurtay A, Abeler-Dörner L, et al. Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing. Science. (2020) 368:eabb6936. doi: 10.1126/science.abb6936

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Cousins S. New Zealand eliminates COVID-19. Lancet. (2020) 395:1474. doi: 10.1016/S0140-6736(20)31097-7

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Gudbjartsson DF, Helgason A, Jonsson H, Magnusson OT, Melsted P, Norddahl GL, et al. Spread of SARS-CoV-2 in the Icelandic population. N Engl J Med. (2020) 382:2302–15. doi: 10.1101/2020.03.26.20044446

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Eguíluz VM, Fernández-Gracia J, Rodríguez JP, Pericàs JM, Melián CJ. Risk of secondary infection waves of COVID-19 in an insular region: the case of the Balearic Islands, Spain. Medrxiv. (2020). doi: 10.1101/2020.05.03.20089623

CrossRef Full Text | Google Scholar