Modelling the effect of the interaction between vaccination and non-pharmaceutical measures on COVID-19 incidence

Since December 2019, the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV- 2) has spread rapidly from Wuhan (China) across the globe, affecting more than 200 countries by mid-2021, with over 190 M reported cases and around 4 M fatalities. During the first year of the pandemic, affected countries implemented a variety of non-pharmaceutical interventions to control virus transmission. In December 2020, countries started administering several authorised vaccines under a limited supply scenario. In this context, the aim of this study was to develop a SEIR-type continuous-time deterministic disease model, to determine the impact of interaction between different vaccination scenarios and levels of protection measures on disease incidence. For this, the model incorporates (i) a protection measure including low (self- protection), medium (mobility limitation), high (closure of indoor facilities) and very high (lockdown) protection levels, (ii) quarantine for confirmed cases, and (iii) vaccination rate and efficacy of four type of vaccines (Pfizer, Moderna, Astra Zeneca or Janssen). The model was verified and evaluated using the response timeline and vaccination strategies and rates in the Basque Country (N. Spain). Once the model performance was validated, different initial phase (when 30% of the population is vaccinated) vaccination scenarios were simulated, including (i) a realistic vaccine limited supply scenario, and (ii) four potential full vaccine supply scenarios where a unique vaccine type is administered. Some differences in disease incidence were found between vaccination scenarios for low and medium-level protection measures. However, regardless of the administered vaccine, a high-level protection scenario is the most effective to control the virus transmission and disease mortality in the studied initial phase of vaccination. The results obtained here may vary in further studies since there may be some unpredictable factors/covariates. With this in mind, the model here could be easily applied to other regions or countries, modifying the strategies implemented and initial conditions.

COVID-19; SARS-CoV-2; SEIR model; non-pharmaceutical measures; vaccination; simulation . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 2 1 Introduction The novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) or coronavirus disease 2019 (COVID-19), was first detected in Wuhan, China, in December 2019, as the cause of a pneumonia of unknown aetiology (Boni et al., 2020; Li et al., 2020; Mwalili et al., 2020). The SARS-CoV-2 rapidly spread all over the globe and on March 11, 2020, the World Health Organization (WHO) declared COVID-19 a global pandemic (Cucinotta & Vanelli, 2020; Mwalili et al., 2020; Rothan & Byrareddy, 2020). After more than a year and a half since the pandemic was declared, according to the WHO, there are 200 M cases reported and 4 M deaths worldwide (ECDC, 2021). The SARS-CoV-2 transmission is through exposure by (i) inhalation of very fine respiratory droplets and aerosol particles released by infected individuals, mostly between people at close range (Tang et al., 2021) (ii) deposition of respiratory droplets and particles on exposed mucous membranes (mouth, nose, or eyes) by direct splashes and sprays, and (iii) touching mucous membranes with hands that have been soiled by touching surfaces with virus ( Liu et al., 2020; Sheng, 2020; WHO, 2021a). The virus transmission in indoor settings has been the main transmission pathway when ventilation is not sufficient (Atalan, 2020; Baghat et al., 2020; Jayaweera et al., 2020; WHO, 2021a).) The main outbreaks have been related to explosive super events in indoor settings or facilities such as family gatherings, long-term health facilities, restaurants, bars and clubs (e.g. Chau et al., 2021), being these principal responsible of the dynamics and shape of the COVID-19 transmission (Althouse et al., 2020). Due to the rapid spread of the virus through these events and the lack of effective pharmaceutical treatments for the disease, particularly at the beginning of the pandemic, important control measures have been implemented worldwide: quarantine of people suspected of being exposed to COVID-19, isolation/quarantine of confirmed cases, use of face masks in public, contact tracing, social distancing, closing of indoor settings (public spaces, restaurants, bars, etc.) and schools and universities, working from home, confinement of regions with a high incidence of the virus, to total lockdown of the country to slow down the COVID-19 outbreak (Atalan, 2020; CDC (2022); ECDC, 2020; Iboi et al., 2020; MacIntyre et al., 2020; WHO, 2020). In addition to these measures, since the pandemic began, the governments of disease-impacted countries have been working on the development and implementation of strategies to return to “normal life”, including the development of several vaccines against COVID-19 (CDC, 2021; WHO, 2022). To date, a variety of vaccines have been approved by the European Medicine . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 3 Agency (EMA), under a conditional marketing authorisation due to the emergency situation, and many others are under development. These vaccines have different efficacy, understood as the percentage reduction in disease incidence, based on clinical trials (EMA, 2021). Vaccination is based on the fact that if a fraction of the population is immune to that pathogen, the susceptible host numbers decrease, so the impact of infected individuals is limited (Randolph & Barreiro, 2020; Sariol & Perlman, 2020). Herd immunity originates when a sufficiently large proportion of the population is immune to the disease (Omer et al., 2020; Randolph & Barreiro, 2020). The percentage of the population that needs to be vaccinated to achieve herd immunity varies with each disease. Several studies have concluded that to achieve COVID-19 herd immunity and relax protection measures around 50% to 70% of the population should be vaccinated (Clemente-Suárez et al., 2020; Kim et al., 2021). The immediate goal of the global COVID-19 vaccination strategy is to minimize deaths, severe disease incidence and reduce the risk of new variants. This requires fully vaccinating at least 70% of the world’s population, accounting for most adults and adolescents and for the vast majority of those at risk of serious disease (WHO, 2021b). Consequently, initial vaccination phases such as the one studied here (around 30% of the population vaccinated) may require the maintenance of some level of non- pharmaceutical protection measures to control disease transmission. In this context, the modelling approach is a determinant tool to analyse COVID-19 disease dynamics and support the development of public health policies (Wong et al., 2021). Most models for the COVID-19 pandemic are single-population continuous compartmental SEIR Kermack-McKendrick-type models, constructed using ordinary differential equation (ODE) systems (Guirao et al., 2020; Li et al., 2020; Tang, Bragazzi, et al., 2020; Wu et al., 2020). Compartmental models are a very common infectious modelling approach where the population is assigned to compartments with labels (S, Susceptible; E, Exposed; I, Infectious; R, Recovered) and individuals may progress between compartments. For COVID-19 models apart from the S and I compartments, E and R compartment are essential since: (i) there is a significant latency period during which individuals have been infected but are not yet infectious themselves, so they are exposed, and (ii) sick individuals recovered from disease they are not infectious (I) and they are immune for some months so they cannot be considered susceptible.” These population disease models may be basic in order to capture certain disease dynamic complexities. However, for any emerging pandemic, they are essential, first to develop the theoretical basis for the understanding of pathogen transmission processes and mechanisms, and second, to explore disease spread control measures. A limitation when modelling COVID-19 transmission is that only confirmed cases are known. There is a fraction of non-reported positive cases, ranging between 10-70% of the total, that correspond to people that do not get . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 4 tested or are asymptomatic to the disease. Hence, models such us the ones developed by the Imperial College of London, estimate that this fraction will show a higher number of cases compared to the reported data (Giattino, 2020). Given the uncertainty surrounding the situation after COVID-19 vaccination programmes, models estimating these unconfirmed cases, such as the one presented here, can be particularly useful for exploring different scenarios of immunisation through the vaccination effect on disease spread limitation. This work is focused on the development of a deterministic SEIR transmission model to analyse the impact of the interaction between different vaccination scenarios, regarding vaccinaton rate and efficacy, and different levels of non-pharmaceutical protection measures (from the use of mask to lockdown) on disease incidence and mortality. The model scenarios are set for the initial phase of the COVID-19 vaccination (i.e. when around 30% of the population is vaccinated) and evaluated on the response timeline of the first and second waves of the pandemic in the Basque Country (N Spain), one of the regions reporting highest disease incidence in Europe. 2 Methods 2.1 Model description and mathematical theory The model here is an extension of a Kermack-McKendrick-type model (Kermack & McKendrick, 1927). It is a deterministic SEIR transmission model that accounts for important characteristic for understanding COVID-19 disease dynamics, such as (i) incubation period, (ii) a protection measure ranging from low-level protection (self-protection; use of a mask, hygiene and social distancing), medium-level protection (mobility limitation), high-level protection (adding indoor facilities closure) and very high-level protection (lockdown), (iii) quarantine for confirmed cases, and (iv) vaccination rate and efficacy. The model is a one-population compartmental model, continuous in time, unstructured in spatial or age terms, and configured to simulate the dynamics of COVID-19 transmission processes caused by susceptible individuals contacting infected individuals or environments with infectious particles released by infected individuals. The compartmental models to describe pathogen transmission are the most frequently used class of models in epidemiology (Diekmann & Heesterbeek, 2000). Individuals can take on a finite number of discrete states, and each state is representative of a subpopulation of individuals at a given time (Table 1). These compartments and states, in consequence, are defined as the variables of the model. These . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 5 variables together with the associated parameters satisfy a system of ODEs describing the dynamics of the host-pathogen system. The model here includes seven compartments (i.e. variables or subpopulations) (Table 1) and each state assumes the following: (1) S stands for susceptible subpopulation that can become exposed to the virus by contact with an infected individual or with infectious particles released by an infected individual; (2) E represents the population exposed to the virus after being in contact with an infected individual or infected environment; (3) I represents the infected subpopulation with individuals coming from the exposed subpopulation after the corresponding incubation period of the virus (five days on average (Lauer et al. 2020; Rǎdulescu et al., 2020)) ; The I subpopulation is assumed to represent asymptomatic cases, non-confirmed and non- isolated symptomatic cases, and cases that are not yet or are not quarantined; consequently this subpopulation can be considered the source of the infection in the model. According to the WHO, although asymptomatic people can spread the virus (Rǎdulescu et al., 2020), they are most infectious in the early stages of a symptomatic stage, so that the majority of infections are caused by symptomatic individuals (WHO, 2020), (4) a fraction of this I subpopulation is the pool of infected individuals, represent confirmed and quarantined people, representing the Q subpopulation. This subpopulation of quarantined people after being diagnosed with COVID-19 includes home isolated and hospitalised patients; (5) R represents the population that has recovered from the disease and is immune to disease during a certain period of immunisation time; (6) V subpopulation represents vaccinated individuals with protection against the virus and (7) D represents individuals that die due to COVID-19, that is, this variable tracks cumulative deaths (Table 1). The variables or subpopulations of the host population are defined with respect to the number of individuals in the studied territory. Thus, the initial population N for the model is 2199711 individuals based on demographic data of the Basque Country Institute of Statistics (BIS) (BIS, 2020). The model specifies an open population where birth of new susceptible individuals is a function of the total population N. Since this is a novel coronavirus, initially everyone is susceptible to COVID-19. The model assumes some individuals were already exposed and infected at simulation day 1 (March 1) (Table 1). These values are obtained by model fitting against real cumulative mortality data (BHD, 2020) and considering that at least 33% of the cases are asymptomatic, not confirmed and able to infect, but not under quarantine (Pollán et al., 2020). . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 6 Table 1. Variable description with model assumptions and initial values obtained by model validation assuming confirmed cumulative mortality data. Initial values in this table are the ones used in the model verification/evaluation on the first wave data. For the verification/evaluation on the second wave data and simulations of vaccination scenarios, changes to these initial values are defined further on. Variable Description and modelling assumptions Initial conditions (individuals) N Population in the Basque Country 2199711 S Population of non-quarantined susceptible individuals 2199671 E Population of susceptible exposed individuals; infected individuals with no symptoms and no infectivity 30 I Population of infected individuals with infective capacity and not quarantined; asymptomatic or not reported cases 10 Q Population of quarantined infected individuals; reported cases 0 R Patients recovered from COVID-19 0 D Individuals deceased due to COVID-19 0 V Vaccinated and protected population against COVID-19 0 Another feature of the COVID-19 virus, is the incubation period, which is relatively long and an individual is able to infect others before being diagnosed (Rǎdulescu et al., 2020). In general, some model features and specific assumptions such as protection measures defined by the parameters in Table 2 may result in some predictive limitations (see section 2.7). Based on all these assumptions and simplifications, the basic model for the transmission dynamics of COVID-19 is given by the following deterministic system of nonlinear differential equations: 2.2 Model Equations The subpopulations of the model satisfy a system of ODEs describing the dynamics of the host- virus association. Variables and parameters of these equations are described in Tables 1 and 2, respectively. The numerical model for this ODE system is programmed in Matlab R2018a. The set of coupled differential equations is solved with a fourth–order predictor corrector scheme, using the Adams Bashforth predictor and the Adams-Moulton corrector. The differential equation system comprises the following differential equations: . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 7 𝑑𝑆 𝑑𝑡 = 𝑏𝑁 − (𝛼 + 𝑣)𝑆 − 𝛽 𝑆𝐼 𝑁 + 𝜎𝑅 + 𝜎𝑉 − 𝑚𝑆 (1) 𝑑𝐸 𝑑𝑡 = 𝛽 𝑆𝐼 𝑁 − 𝛾𝐸 − 𝑚𝐸 (2) 𝑑𝐼 𝑑𝑡 = 𝛾𝐸 − 𝑟1𝐼 − 𝑞𝐼 − 𝑚𝐼 (3) 𝑑𝑄 𝑑𝑡 = 𝑞𝐼 − 𝑟2𝑄 − 𝑑𝑄 − 𝑚𝑄 (4) 𝑑𝑅 𝑑𝑡 = 𝑟1𝐼 + 𝑟2𝑄 − 𝜎𝑅 − 𝑚𝑅 (5) 𝑑𝐷 𝑑𝑡 = 𝑑𝑄 (6) 𝑑𝑉 𝑑𝑡 = 𝑣𝑆 − 𝜎𝑉 − 𝑚𝑉 (7) Equation (1): The change in the number of susceptible individuals S, is a balance between (i) the loss of individuals due to protection measures (use of mask, social distancing, mobility restrictions, indoor settings closure, lockdown) and vaccination, virus transmission and

mortality and (ii) the gain of individuals from births, and recovered and vaccinated individuals who lost immunisation. Equation (2): The change in the number of exposed individuals E, is a balance between (i) the gain of individuals due to virus transmission, and the (ii) loss of individuals because of

mortality or the end of the incubation period in exposed individuals. Equation (3): The change in the number of infected individuals I, is a balance between the (ii) gain of individuals due to the end of incubation period in exposed individuals, and the (ii) loss of individuals because of background mortality, recovered individuals and the subtraction of the proportion of isolated/quarantined individuals. Equation (4): The change in the number of quarantined individuals Q, is a balance between (ii) the gain of individuals due to the proportion of the infected subpopulation who are symptomatic and theoretically recorded as confirmed cases and quarantined individuals, and (ii) the loss of individuals because of disease mortality, recovered individuals and background mortality. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 8 Table 2. Baseline model parameters with a brief description and default values used for the model. Par am eter Description Value Range Unit s

 Incubation rate 2.0× 10-1 0.1 - 0.3 day-1 Amira et al. (2020) Dhouib et al. (2021) Lauer et al. (2020) Rǎdulescu et al. (2020) r1 Recovery rate for I 7.0× 10-2 0.5 – 1 ×10-1 day-1 Rǎdulescu et al. (2020) Yang et al. (2021) r2 Recovery rate for Q 5.0× 10-2 q Quarantine rate 6.7× 10-1 0.5 -0.8 day-1 Tang, Xia, et al. (2020). Yang et al. (2021)  Infection rate 1.1 0.6 - 1.7 day-1 Lin et al. (2020) b Birth rate 2.7 × 10-5 EUSTAT (2022a) d Disease mortality rate 5.7× 10-3 0.1 - 0.6 ×10-1 day-1 López & Rodó, (2021) Rǎdulescu et al., (2020) m Natural mortality rate 2.7 × 10-5 day-1 EUSTAT (2022b)  Vaccine protection rates (Pfizer, Moderna, Astra Zeneca, Janssen) 1.3 × 10-3 day-1 Estimated in this study (Table 3) 1.8 × 10-4 3.1 × 10-5 9.8 × 10-5  Protection rate: Low (mask use, hygiene, social distancing), medium (mobility restrictions), high (mobility and indoors closure) very high (lockdown) 5.0 × 10-3 0.3 - 3 ×10-2 day-1 López & Rodó, (2021) Estimated in this study 7.5× 10-3 10 × 10-3 25 × 10-3  Immunity loss rate 1.1× 10-2 day-1 Dan et al.(2021) Equation (5): The change in the number of recovered individuals R, is a balance between the gain of recovered individuals from infection including quarantined individuals, and the loss of recovered individuals who lost immunisation and background mortality. Equation (6): The change in the number of deceased infected individuals D (cumulative mortality), is represented by the gain of individuals due to disease mortality from quarantined (and eventually hospitalised) cases. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 9 Equation (7): The change in the number of vaccinated and protected individuals V against COVID-19, is a balance between the gain of individuals due to immunisation through vaccination, and the loss of individuals due to the loss of immunisation, after a certain period of time, and background mortality. 2.3 Data Epidemiological data Data on infection-confirmed cases, cumulative deaths and vaccination in the Basque Country, were obtained from the open data system on the BHD website (BHD, 2020). The model uses realistic initial conditions for the variables and keeps parameters in the range of recent research findings (see Table 2). The new daily confirmed cases and cumulative mortality data were used for comparison with modelling results and for model verification and evaluation. Vaccination data The first mass vaccination programme in the Basque Country started on December 27, 2020, reporting by then 117000 cases and 3030 fatalities (BHD, 2020). Since then to June 8, 2021, four vaccines were administered in the Basque Country: Pfizer, Moderna, Astra Zeneca and Janssen. For two-dose vaccines, the time interval between doses is three weeks for Pfizer, four weeks for Moderna and 12 weeks for Astra Zeneca. The model uses the initially reported vaccine efficacy values (EMA, 2021) in Table 3 and 4. This efficacy values are reached after seven days from the booster dose for Pfizer and after 14 days for the Moderna and Astra Zeneca vaccines. For the single dose Janssen vaccine, the reported efficacy is achieved after 14 days. As the vaccine efficacy is very labile and dynamic depending in many factors such as virus variants, updated data on vaccine efficacy for further modelling studies can be obtained in product reports of vaccines in EMA (2022) Since the start of the vaccination programme on the 27th of December, up until the 8th of June, 643,978 complete doses of COVID-19 vaccines were administered, distributed by vaccine type as in Table 3. This administration strategy, due to vaccine supply limitations in the Basque Country, envisaged the administration of four vaccines with the following percentages in terms of population: 77.2% Pfizer, 11.2% Moderna, 3.0 % Astra Zeneca and 8.7% Janssen (BHD, 2020). . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 10 Table 3. Estimated vaccination and protection rates using the Basque Countries’ vaccination data (BHD, 2020) and the initial efficacy rates reported from EMA (EMA, 2021). Vaccine type Complete doses (CD) (individuals) Vaccination rate (VR) (individuals/day) Vaccine efficacy (VE) Vaccination protection rate (VPR) (model) Pfizer/BionTech 496894 3105 94.6% 1.3  10-3 Moderna 72164 451 93.6% 1.8  10-4 Astra Zeneca 19079 119 60.0% 3.1  10-5 Janssen 55841 349 67.0% 9.8  10-5 Table 4. Estimated vaccination and protection rates for unique vaccine administration scenarios and the efficacy reported from EMA (EMA, 2021). Vaccine type Complete doses (CD) (individuals) Vaccination rate (VR) (individuals/day) Vaccine efficacy (VE) Vaccination protection rate, (VPR) (model) Pfizer/BionTech 643978 4025 94.6% 1.7  10-3 Moderna 93.6% 1.6  10-3 Astra Zeneca 60.0% 1.0  10-3 Janssen 67.0% 1.1  10-3 The average number of vaccines administered per day (vaccination rate, VR), was estimated from data obtained from the vaccination bulletin of the BHD (BHD, 2020). According to the number of complete doses administered and the efficacy of each vaccine, a vaccine-specific vaccination protection rate was estimated for the model (Table 3) as follows: VR = 𝐶𝐷 𝑇 VPR = 𝑉𝑅 · 𝑇 𝑁 (𝑇+𝐷) · 𝑉𝐸, where T is 164 days, from December 27 to June 8, and D is the number of days needed after the second dose was administered for the vaccine to achieve peak efficacy. N is the total population in the Basque Country (N=2199711). 2.4 Descriptive analysis of the epidemic curve for the Basque Country The change in the number of daily COVID-19 new cases in the Basque Country over time, from March 1, 2020 to June 8, 2021 is described in terms of the non-pharmaceutical interventions adopted by the BHD (Figure 1). Results including simulations for model verification/validation and simulations for vaccination scenarios were discussed in terms of this initial descriptive analysis . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 11 2.5 Model verification and evaluation The model was verified and evaluated with simulations for the first and second waves of the pandemic in the Basque Country using the parameter values in Table 2. Model verification consisted of showing that the simulation model was correct, complete, and coherent. This required analysing (1) static tests involving a structured examination of the formulas, algorithms, and codes used to implement the model, by several reviewers experts in population modelling (see acknowledgements section), and (2) dynamic tests, where the computer program was run under different conditions of disease mortality and quarantine rate. This analysis were critical to ensure that results produced were correct, according to the conceptual model, and consistent with expectations of reviewers. Model evaluation was carried out against confirmed cases and cumulative mortality. Simulations and descriptions of model verification and evaluation are given in the results section. 2.6 Modelling vaccination scenarios Considering the degree of protection of each vaccine in terms of the estimated vaccination protection rate together with different nonpharmaceutical interventios (low, medium and high protection levels (Table 1) the effects of five different vaccination strategies were tested within the model. Realistic scenario with limitations on vaccine supply In the first scenario, the model tested the effect on infected cases, based on the vaccination strategy followed in the Basque Country. This realistic limited vaccine supply scenario envisages the administration of four vaccines, with the percentages in terms of vaccinated population estimated from Table 3. The model uses vaccine protection rates estimated in Table 3. Realistic scenarios with full vaccine supply Four full vaccine supply scenarios were simulated where a unique vaccine type administration is possible. For these simulations, the model uses new vaccination protection rates estimated in Table 4, as in Vaccination data section, considering that the population receives a unique vaccine. 2.7 Model limitations The four nonpharmaceutical protection levels (Table 2) are considered fully and successfully applied. However, this is a simplification of reality and may result in an underestimation of . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 12 disease incidence; since violation of low-level protection measures such as the use of mask or social distancing may be much higher than strong measured such as lockdown. For such a rapid spread of COVID-19 the population is assumed to be constant in terms of the difference between births and natural deaths, i.e. the birth rate is assumed to be the same as the natural death rate (Table 2). In the case of the studied region the birth rate is slightly lower than the mortality rate (EUSTAT, 2022a, b), however this assumption is overestimating the S population. A seroprevalence study carried out in Spain estimated 33% of the cases as asymptomatic (Pollán et al., 2020). The quarantine rate in the model (0.67 day -1, see Table 2) was estimated based on this study result and assuming that asymptomatic are non-confirmed cases and hence, non-quarantined people. This assumption could be underestimating the quarantined population particularly beyond third wave when testing was common to close contacts of Covid-19. When modelling realistic vaccine scenarios with vaccine full supply caution is required to interpret the results (see section 3.3) quantitatively; since the model is dealing with multiple dimensions, latent covariates, a local region without %100 certain data and four vaccine types. 2.8 Descriptive analysis of the curve of infection cases The curve of new daily confirmed cases by the BHD (BHD, 2020) (Figure 1) is described below in terms of the non-pharmaceutical interventions´ chronological sequence. This detailed analysis is essential for the selection of simulation cases with different protection measures and later discussion. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 13 The first epidemic wave in the Basque Country in terms of confirmed cases started with the first case of COVID-19 on February 28, 2020. This first wave reached its peak on the March 24, with 723 confirmed new cases (Figure 1). Previously, on March 13, schools were closed in all regions of Spain including the Basque Country and a nationwide lockdown (confinement) was declared, banning all public events. As the confirmed cases were decreasing, on April 13 the partial lifting of the lockdown started, allowing people to go out by age groups at certain times during the day. On May 4, the lockdown was totally lifted, starting as a phase 0, called the “new normality”. Mobility all around the country was permitted in summer and indoor public facilities were open. On July 8, mask use in public areas became mandatory for people older than six years old. By mid-July, daily cases started to increase consistently, and the second wave reached its peak on August 28, with 886 confirmed new cases (Figure 1). The peak of the third wave reached on November 5, with 1547 new confirmed cases. When cases started to increase on October 27, the activity and mobility was limited to daytime, and it was prohibited from 11pm to 6am. Also, the region of the Basque Country was closed and travel between regions was forbidden. Only groups with a maximum of six people were allowed. And on November 7, bars and restaurants were closed until December 12, likewise limiting capacity at indoor places. Figure 1. Confirmed cases in the Basque Country from March 1, 2020, to June 8, 2021 (BHD, 2020). . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 14 A fourth wave peak was recorded on January 27, 2021, with 1274 new cases. The fifth wave started in mid-March, reaching its peak on April 21, with a maximum of 1013 new cases. 2.9 Model verification and evaluation From the set of simulation results obtained during model verification and evaluation, eight simulation cases for the first two waves of the pandemic in the Basque Country are described here: Cases 1-4 (Figure 2) were simulated for the first wave and Cases 5-8 for the second wave (Figure 3). The simulations were run with the initial conditions described in Table 1 and parameter values described in Table 2. Changes to these values for specific simulations are described for each case. 2.9.1 First wave Realistic scenarios with varying disease mortality (Figure 2A) and quarantine rate (Figure 2B) were simulated to verify and evaluate the performance of the model, focusing on the dynamics of the host-pathogen system in terms of cumulative mortality. In order to simulate the protection measure of the lockdown a tangential function was considered where the protection rate changed from a low protection level (5 ×10-3 day-1, before lockdown) to a high protection level (25 ×10-3 day-1, during lockdown) (see Table 2). . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 15 Figure 2. Model verification/evaluation against the first wave of the pandemic in the Basque Country. (A) Case 1: Simulation of cumulative mortality with varying disease mortality, (B) Case 2: Simulation of cumulative mortality with varying quarantine rate (day -1), (C) Case 3: Simulation of infected, recovered and deaths (cumulative mortality), (D) Case 4: Model evaluation against new confirmed daily cases. Realistic initial conditions (Table 1) and parameter values (Table 2) were used. Simulation started on March 1 (day 1) and ended on July 8 (day 130). Increasing the disease mortality rate from low (2 ×10-3 day-1) to high (6 ×10-3 day-1) in Case 1

in an expected increased cumulative mortality conforming to expectations (Figure 2A). For this simulation case, model evaluation against real cumulative data shows that the simulation with the highest disease mortality rate (6 ×10-3 day-1) results in 1700 deaths, showing the model has the best fit to the 1687 deaths confirmed by the BHD in the first wave (BHD, 2020). Increasing the quarantine rate from low (5.7 ×10-1 day-1) to high (6.7 ×10-1 day-1) in Case 2

in an expected reduced cumulative mortality (Figure 2B) conforming to model behaviour expectations. In this Case 2, model evaluation against real cumulative data shows the following: the simulation with the higher quarantine rate results in 1700 deaths, showing the model has the best fit to the 1687 deaths confirmed by the BHD in the first wave (BHD, 2020). Once the model was verified and evaluated against mortality data, Case 3 was run to verify the model in terms of change of infected, recovered and death subpopulations with time during the first wave of the pandemic (Figure 2C). The behaviour of the model in these terms conforms to . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 16 expectations regarding the number of recovered individuals (around 21000, approximately estimated cases – confirmed deaths) and fits to deaths (1687) confirmed by the BHD in the first wave of the pandemic (BHD, 2020). Finally, Case 4 evaluated the model against new confirmed daily cases (Figure 2D). The model estimates 22230 cases as the number of true infections for the simulation period; about twice as high as the 13862 confirmed cases (BHD, 2020). 2.9.2 Second wave Similarly to the first wave, realistic scenarios with varying disease mortality (Figure 3A) and quarantine rate (Figure 2A) were simulated to verify and evaluate the performance of the model focusing on the dynamics of the host-pathogen system in terms of cumulative mortality. The simulations were run with specific initial conditions described in Table 1 and parameter values described in Table 2, except for (1) initial infected population of 40 individuals which was estimated from real data considering 33% of not confirmed infected people and (2) a low/medium-level protection rate of 6.5 × 10-3 day -1; lower than in the first wave due to the end of lockdown, opening of indoors, mobility being allowed all over the country and no restrictions on international travellers entering the country. In this second wave, simulation values for mortality rate were fitted using real deaths and quarantined individuals. These values were lower than those for the first wave since the disease incidence was much higher in young people, with a much lower fatality rate (BHD, 2020). Increasing disease mortality rate from low (1 × 10-3 day-1) to high (3 × 10-3 day-1) in Case 5

in an expected increased cumulative mortality as expected (Figure 3A). For this case, model evaluation against real cumulative data shows that the simulation with the lower disease morality rate (1 × 10-3 day-1day-1) results in 335 deaths showing the model has the best fit to the 340 deaths confirmed by the BHD in the first wave (BHD, 2020). Increasing the quarantine rate from low (5.7 × 10-1 day-1) to high (6.7 × 10-1 day-1) in Case 6

in an expected reduced cumulative mortality (Figure 3B) as in Case 2, confirming the adequate behaviour of the model. In this Case 6, the simulation with the higher quarantine rate

in 335 deaths, showing the model has the best fit to the 340 deaths confirmed by the BHD in the first wave (BHD, 2020). . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 17 Figure 3. Model verification/evaluation against the second wave of the pandemic in the Basque Country (A) Case 5: Simulation of cumulative mortality with varying disease mortality, (B) Case 6: Simulation of cumulative mortality with varying quarantine rate (day -1), (C) Case 7: Simulation of infected, recovered and deaths (cumulative mortality), (D) Case 8: Model evaluation against new confirmed daily cases. Realistic initial conditions (Table 1) and parameter values (Table 2) were used. The simulation started on July 9 (day 131) and ended on September 28 (day 220). Case 7 verified the behaviour of the model in terms of changes in the infected, recovered and death subpopulations with time during the first wave of the pandemic (Figure 3C). The number of recovered individuals (around 46200) conforms to expectations considering the estimated cases and confirmed deaths, while estimated deaths (335) fits the number of fatalities confirmed by the BHD in the second wave of the pandemic (340) (BHD, 2020). Finally, in Case 8 the model estimated 45300 cases as the number of true infections for the simulation period (Figure 3D); higher than the 31000 confirmed cases (BHD, 2020). 2.10 Vaccination scenarios Once the model was verified and evaluated, the impact of the vaccination strategy was tested considering five simulation scenarios for the fifth wave of the pandemic from March 7 to June 8, 2021. For these simulations, the varying non-pharmaceutical protection rates were from low- level to high-level protection rates as in Table 2: (i) the low-level protection rate  was estimated to be 5 × 10-3 day-1, including the use of masks, social distancing, opened indoor facilities, restaurants and bars, and easing of mobility limitations as they were initially relaxed . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 18 inside the Basque territory and eventually all around Spain to promote tourism during the Easter holidays and before summer, (i) the medium-level protection rate of =7.5 × 10-3 day-1 represented the previous scenario with no regional and national mobility restrictions, and (iii) the high-level protection rate of =1 × 10-2 day-1 adds to the previous scenario the closure of indoor public and private facilities . In addition, in these simulations the transmission rate was increased to 1.15 day-1, since the Delta variant of the virus was known to spread significantly faster than the original version of the virus (Li et al., 2021). Initial conditions are those in Table 1 with changes in susceptible (S=2000000), infected cases (I=200) and vaccinated (V=52500) subpopulations due to the course of the disease dynamics. 2.10.1 Simulation 1. Limited vaccine supply scenario: combination of Pfizer, Moderna, Astra Zeneca and Janssen In this simulation, the model tries to mirror the vaccine strategy with a combination of Pfizer, Moderna, Astra Zeneca and Janssen vaccines followed by the BHD in the Basque Country. The vaccine-specific protection rates used in the model are those estimated in Table 4. The cumulative new daily confirmed cases for the simulation period were 48577 (black dots in Figure 4A). The model estimates 57000 cases and 55000 recovered individuals. Regarding fatalities, the simulation result (460 deaths) fits confirmed deaths (462) by the BHD for the simulation period (BHD, 2020). For this vaccination scenario, simulations to test the effect of increasing the protection rate on new daily cases and cumulative mortality were run. The protection rate was increased from the adopted protection rate =5 × 10-3 day-1 to =1 × 10-3 day-2, in accordance with the strengthening of limitations in mobility and closure of indoor facilities such as restaurants and bars. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 19 Figure 4. Limited vaccine supply scenario with a combination of Pfizer, Moderna, Astra Zeneca and Janssen vaccines: (A) Model evaluation against new confirmed daily cases. Realistic initial conditions (Table 1) and parameter values (Table 2) were used. The simulation started on March 7, 2021 (day 380) and ended on June 8, 2021 (day 460), (B) Simulation of infected, recovered and deaths (cumulative mortality), during vaccine administration. Figure 5. Simulations of new daily cases (A) and cumulative mortality rate (B) from March 7, 2021 (day 380) to June 8, 2021 (day 460) with varying protection rate (day-1) in a limited vaccine supply scenario with a combination of Pfizer, Moderna, Astra Zeneca and Janssen vaccines. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 20 The highest number of new cases occurs with the lowest protection rate (=5 × 10-3 day-1) where fewer restrictions are applied (Figure 5A). This number reduces to half when the protection measures increase to = 7.5 × 10-3 day-1, decreasing even more if the protection rate is set to =1 × 10-2 day-1. The same behaviour is observed in the cumulative mortality or death cases (Figure 5B). The following scenarios (3.3.2 – 3.3.5) are full supply scenarios where the type of vaccine can be chosen. Thus, the simulations contemplate the administration of a unique vaccine type, considering (i) the number of complete doses by June 8, the same as in the real limited vaccine supply scenario, and (ii) a vaccination rate that is the same for all vaccine types (see Table 4). For these vaccination scenarios, as in the first scenario, simulations to test the effect of increasing the protection rate (from =0.004 to =0.01) on new daily cases and cumulative mortality were run. 2.10.2 Simulation 2. Pfizer scenario In this simulation (Figure 6A), Pfizer is the unique vaccine administered to the population, with a vaccine protection rate of 1.7 × 10-3 day-1 (Table 4). The model estimates 51435 cumulative new cases and 415 deaths with the lowest protection rate; a lower number of fatalities compared to that obtained in the realistic limited vaccine supply scenario (462) (Simulation 1). Similarly, to simulation 1, an increased protection rate reduces the incidence of cases and mortality (Figure 6A), particularly when the protection rate is at its maximum with strong limitations in mobility and closure of indoor facilities such as restaurants and bars. In this protection scenario, the number of cases is about five times lower than that for a low protection rate, with no cases in the last 30 days. The cumulative mortality decreases from 415 to 85 deaths. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 21 Figure 6. Simulation of new daily cases and cumulative mortality (deaths) with a varying protection rate in a unique vaccine administration scenario: A) Pfizer, B) Moderna, C) Astra Zeneca and D) Janssen. 2.10.3 Simulation 3. Moderna scenario In this simulation (Figure 6B), Moderna is the unique vaccine administered to the population with a vaccine protection rate of 1.6 × 10-3 day-1 (Table 4). The model estimates 56,500 cases in the simulation period and 458 deaths with the lowest protection rate. The impact of the high protection scenario (green line) on disease dynamics is also high; cumulative mortality decreases from 458 to 87 deaths. 2.10.4 Simulation 4. Astra Zeneca scenario This simulation represents a scenario with the Astra Zeneca vaccine as the unique vaccine administered (Figure 6C), with a vaccine protection rate of 1 × 10-3 day -1 (Table 4). The model in this case estimates 77,900 cases in the simulation period and 600 deaths with the lowest protection rate. Here, for the medium protection level the model estimates 265 deaths. However, for the higher protection rate, differences between the responses to vaccines are not significant. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 22 2.10.5 Simulation 5. Janssen scenario This simulation (Figure 6D) shows the Janssen vaccine as the unique vaccine administered with a vaccine protection rate of 1.1 × 10-3 day -1(Table 4). For the highest protection rate, new daily cases and cumulative mortality, are similar to those observed for the other vaccines. Nonetheless, for the medium and, particularly, for the lowest protection rate, the Janssen vaccine estimates 62,030 cases in the simulation period and 550 deaths. 3 Discussion and conclusions This contribution covers the theoretical and mathematical basis for modelling dynamics and epidemiology of COVID-19, specifically focusing on the effect of the interaction between the initial phase of the vaccination and non-pharmaceutical protection measures such as self- protection, mobility restrictions, closure of indoor facilities and lockdown. The Kermack and McKendrick (1927) epidemiological theory was adapted to build a SEIR deterministic model for COVID-19 to assess the impact of this interaction on infection cases and disease mortality. The model was verified and validated using the response timeline, vaccination strategies and non-pharmaceutical interventions implemented in the Basque Country (N Spain). Although robust validation of the model predictions is needed, initial results and evaluation show the potential of the model to be easily modified to match other regions’ or countries’ timelines, or the different response strategies implemented in other countries. The waves of the epidemic curve (confirmed cases) in this region can be discussed in terms of the non-pharmaceutical interventions and disease management. Fifteen days after the first case was confirmed in the Basque Country, schools were closed in all regions of Spain and a nationwide lockdown (confinement) was declared, banning all public events. In this context, during the first wave Covid-19 testing management and coverage was limited. Tests were done only, for those who presented symptoms such as fever and cough. People who did not seek medical attention were tested very rarely. The second wave peak was reached in late summer when mobility all around the country was permitted, international travellers entered the country without restrictions, and indoor public and private facilities were open. The peak of the third wave reached on November linked with the return to schools and work, and the increase of indoor activities. A fourth wave peak was recorded on January 2021 linked with Christmas holidays despite the limitations imposed by the government, limiting gatherings of people and mobility in the territory. The fifth wave was linked with the easing of mobility limitations inside the Basque territory to promote tourism during the Easter holidays. . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 23 The projections explored here for model verification and validation do not differ much from the reported cumulative mortality data and are consistent with the descriptive analysis of the epidemic curve and disease dynamics found by other similar modelling studies in the Basque Country (López & Rodó, 2020). However, estimated new daily cases are significantly higher than those reported by the BHD (BHD, 2020). This result is consistent with the fact that confirmed cases may be undercounted (Giattino, 2020) since one of the key limitations when modelling this disease is that the reported cases only become confirmed cases by a test, and there are a substantial proportion of infected people that never get tested, particularly in the first wave, because they were asymptomatic or never sought medical assistance. The model here, as well as others of this type or more sophisticated ones, use confirmed cases and deaths, testing rates, and a range of assumptions and epidemiological knowledge to estimate this proportion and consequently show a higher number of cases compared to the reported data (Giattino, 2020). This is also consistent with the seroprevalence study carried out in Spain, where around 33% of the infected cases were asymptomatic (Pollán et al., 2020). On the other hand, there is a deviation ratio and as expected if the number of cases increase the range between model prediction and real data increases. The simulation explored here for the limited vaccine supply scenario confirms that the model is correctly validated against the real data. The performance when increasing the protection measures from low-level to high-level non-pharmaceutical protection measures, follows the expected decreasing trend of COVID-19 cases and cumulative mortality. Comparing this scenario to the full vaccine supply scenarios, results suggest that the ideal scenario for limiting the impact of OCVID-19 is the one combining vaccination and high protection levels for non- pharmaceutical measures. There is not a big variation between vaccines in terms of cases and cumulative mortality when the protection rate is high (Figures 5 and 6). Differences in disease incidence response between vaccines need to be taken with caution since there may be some unpredictable latent covariates as described in section 2.7. asaas described in section 2.7 Model limitations. Overall, the results suggest that in an initial vaccination phase (30%-50% of the population is vaccinated) COVID -19 incidence, as measured on daily cases and cumulative mortality, importantly decreases when vaccination and a high level of non-pharmaceutical interventions are in place. That is, in the first vaccination phase, together with vaccination, strong mobility restrictions and closure of indoor facilities such as public spaces, restaurants and bars are critical to significantly control disease outbreaks. When the adopted measures are in the low level (no mobility restrictions and indoor facilities open), COVID-19 cases and deaths remain too high to contain the outbreak. As a positive result, it seems that the number of fatalities has decreased with respect to previous waves. This may be explained by the vaccine programme, which . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 24 prioritises the most vulnerable people by age (BHD, 2020). However, cases are still high, thus, regarding COVID-19 cases, it bears repeating that model results also suggest that initial phase vaccination may synergise with other non-pharmaceutical measures, until the proportion of the immune population increases. The relevance of these results lies in the fact that they support research regarding the relative importance of the non-pharmaceutical measures and vaccination involved in the termination of the epidemic. One limitation of the studied model approach is the assumption that most parameters, except the protection measure in the first wave, take fixed values independent of time. The assumption of constancy in time has the advantage of simplifying the models and facilitates its use. However, both the prevalence of infection and the transmission of the virus may be tied to environmental conditions (Eslami & Jalili, 2020). For the sake of simplicity and the obtention of a preliminary picture, in this model, population has not been divided by age groups. This age-structure should be considered to improve the present model as the vaccination programme has been prioritised by age groups. An age-structured version of this model would give a more accurate picture of the virus transmission (Foy et al., 2021), same as considering the new variants of the virus that could directly impact in the vaccine induced protection and the transmission rate (Moore, 2021). Given the current pandemic caused by the transmission of SARS-CoV-2, the construction of mathematical models such as this one, based on epidemiological data, has allowed us to describe the interactions, explain the dynamics of infection, as well as predict possible scenarios that may arise with the introduction of measures such as social distancing, the use of masks, mobility limitations and vaccination programmes. Mathematical models are highly relevant for making objective and effective decisions to control the disease. These models have supported and will continue to contribute to the selection and implementation of programmes and public policies that prevent associated complications, slow down the spread of the virus and minimise the appearance of severe cases of disease that may collapse health systems. 4 Acknowledgments This investigation was conducted under the framework of the Master in Public Health of the University of Basque Country (UPV/EHU) (Department of Preventive Medicine and Public Health). We appreciate this support. The COVID-19 incidence and mortality data were obtained from the Basque Health Department (BHD) (Osakidetza) (Open Data Euskadi web page: https://opendata.euskadi.eus/catalogo/-/evolucion-del-coronavirus-covid-19-en-euskadi/). The contents of this manuscript do not necessarily reflect the point of view of the BHD and in no ways anticipate the BHD´s future policy in this area. The manuscript benefited from helpful . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 25 discussions with Dr. Aitana Lertxundi and Dr. Naroa Kajarabille (Department of Preventive Medicine and Public Health, UPV/EHU) and Ane Murueta (Department of Neurosciences, UPV/EHU), and examination of the model structure by Dr. Tal Ben-Horin (Department of Clinical Sciences, North Carolina State University) and Morganne Igoe (Department of Mathematics, University of Tennessee). 5 References Althouse, B. M., Wenger, E. A., Miller, J. C., Scarpino, S. V., Allard, A., Hebert-Dufresne, L., & Hu, H. (2020). Superspreading events in the transmission dynamics of SARS-CoV-2: Opportunities for interventions and control. PLoS Biology, 18(11), e3000897. https://doi.org/10.1371/journal.pbio.3000897 Amira, F., Hamzah, B., Lau, C. H., Nazri, H., Ligot, D. V., Lee, G., Liang Tan, C., Khursani Bin, M., Shaib, M., Hasanah, U., Zaidon, B., Abdullah, A. B., Chung, M. H., Ong, C. H., Chew, P. Y., Salunga, R. E., & Hamzah, A. B. (2020). CoronaTracker: World-wide COVID-19 Outbreak Data Analysis and Prediction CoronaTracker Community Research Group Correspondence to Fairoza. https://doi.org/10.2471/BLT.20.251561 Atalan, A. (2020). Is the lockdown important to prevent the COVID-9 pandemic? Effects on psychology, environment and economy-perspective. Annals of Medicine and Surgery, 56, 38–42. https://doi.org/10.1016/j.amsu.2020.06.010 Bhagat, R., Davies Wykes, M., Dalziel, S., & Linden, P. (2020). Effects of ventilation on the indoor spread of COVID-19. Journal of Fluid Mechanics, 903, F1. doi:10.1017/jfm.2020.720 BHD. (2020). Evolución del coronavirus (COVID-19) en Euskadi - Conjunto de datos de Open Data Euskadi - Euskadi.eus. Retrieved July 6, 2021, from https://opendata.euskadi.eus/catalogo/-/evolucion-del-coronavirus-covid-19-en-euskadi/ BIS. (2020). Población. Retrieved July 6, 2021, from https://www.eustat.eus/estadisticas/tema_159/opt_0/ti_poblacion/temas.html Biswas, M. H. A., Paiva, L. T., & De Pinho, M. (2014). A seir model for control of infectious diseases with constraints. Mathematical Biosciences and Engineering, 11(4), 761–784. https://doi.org/10.3934/mbe.2014.11.761 Boni, M. F., Lemey, P., Jiang, X., Lam, T. T. Y., Perry, B. W., Castoe, T. A., Rambaut, A., & Robertson, D. L. (2020). Evolutionary origins of the SARS-CoV-2 sarbecovirus lineage responsible for the COVID-19 pandemic. Nature Microbiology, 5(11), 1408–1417. https://doi.org/10.1038/s41564-020-0771-4 CDC (2021) Developing Covid-19 Vaccines. Retrieved January 25 2022 from . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 26 https://www.cdc.gov/coronavirus/2019-ncov/vaccines/distributing/steps-ensure- safety.html CDC (2020) Nonpharmaceutical Interventions (NPIs). Retrieved January 25, 2022 from https://www.cdc.gov/nonpharmaceutical-interventions/index.html Chau, N., Hong, N., Ngoc, N. M., Thanh, T., Khanh, P., Nguyet, L….Van Tan, L. (2021). Superspreading Event of SARS-CoV-2 Infection at a Bar, Ho Chi Minh City, Vietnam. Emerging Infectious Diseases, 27(1), 310-314. https://doi.org/10.3201/eid2701.203480. Clemente-Suárez, V. J., Hormeño-Holgado, A., Jiménez, M., Benitez-Agudelo, J. C., Navarro- Jiménez, E., Perez-Palencia, N., Maestre-Serrano, R., Laborde-Cárdenas, C. C., & Tornero-Aguilera, J. F. (2020). Dynamics of population immunity due to the herd effect in the COVID-19 pandemic. In Vaccines (Vol. 8, Issue 2). MDPI AG. https://doi.org/10.3390/vaccines8020236 Cucinotta, D., & Vanelli, M. (2020). WHO declares COVID-19 a pandemic. In Acta Biomedica (Vol. 91, Issue 1, pp. 157–160). Mattioli 1885. https://doi.org/10.23750/abm.v91i1.9397 Dan, J. M., Mateus, J., Kato, Y., Hastie, K. M., Yu, E. D., Faliti, C. E., Grifoni, A., Ramirez, S. I., Haupt, S., Frazier, A., Nakao, C., Rayaprolu, V., Rawlings, S. A., Peters, B., Krammer, F., Simon, V., Saphire, E. O., Smith, D. M., Weiskopf, D., … Crotty, S. (2021). Immunological memory to SARS-CoV-2 assessed for up to 8 months after infection. Science, 371(6529). https://doi.org/10.1126/science.abf4063 Dhouib, W., Maatoug, J., Ayouni, I., Zammit, N., Ghammem, R., Fredj, S. Ben, & Ghannem, H. (2021). The incubation period during the pandemic of COVID-19: a systematic review and meta-analysis. In Systematic Reviews (Vol. 10, Issue 1, pp. 1–14). BioMed Central Ltd. https://doi.org/10.1186/s13643-021-01648-y Diekmann, O., & Heesterbeek, J. A. P. (2000). Mathematical epidemiology of infectious diseases: model building, analysis... https://books.google.es/books?hl=es&lr=&id=5VjSaAf35pMC&oi=fnd&pg=PR11&dq=di ekmann+models+in+epidemiology&ots=c3Pa0LaviK&sig=6cN0prU0MkL4WVolIXOiIo pHNRo#v=onepage&q=diekmann models in epidemiology&f=false Dong, E., Du, H., & Gardner, L. (2020). An interactive web-based dashboard to track COVID- 19 in real time. In The Lancet Infectious Diseases (Vol. 20, Issue 5, pp. 533–534). Lancet Publishing Group. https://doi.org/10.1016/S1473-3099(20)30120-1 ECDC. (2020). Guidelines for non-pharmaceutical interventions to reduce the impact of COVID-19 in the EU/EEA and the UK. European Centre for Disease Prevention and Control. Stockholm. ECDC. (2021). COVID-19 situation update worldwide, as of week 4, updated 4 February 2021. European Centre for Disease Prevention and Control. Retrieved June 6, 2021, from https://www.ecdc.europa.eu/en/geographical-distribution-2019-ncov-cases . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 27 EMA. (2021). COVID-19 vaccines: authorised | European Medicines Agency. Retrieved January 10, 2021, from https://www.ema.europa.eu/en/human-regulatory/overview/public- health-threats/coronavirus-disease-covid-19/treatments-vaccines/vaccines-covid-19/covid- 19-vaccines-authorised%0Ahttps://www.ema.europa.eu/en/human- regulatory/overview/public-health-threats/cor EMA. (2022). COVID-19 vaccines: authorised. European Medicines Agency. Retrieved January 21, 2022, from https://www.ema.europa.eu/en/human-regulatory/overview/public-health- threats/coronavirus-disease-covid-19/treatments-vaccines/covid-19-vaccines Eslami, H., & Jalili, M. (2020). The role of environmental factors to transmission of SARS- CoV-2 (COVID-19). AMB Express 2020 10:1, 10(1), 1–8. https://doi.org/10.1186/S13568- 020-01028-0 EUSTAT (2021). Euskal Estatistika Erakundea: Population. Retrieved July 6, 2021, from https://www.eustat.eus/estadisticas/tema_159/opt_0/ti_poblacion/temas.html EUSTAT (2022a) Euskal Estatistika Erakundea: Live born and rate of birth per 1,000 inhabitants by province. 1975-2020. Retrieved on January 21, 2022 https://www.eustat.eus/elementos/ele0005700/ti_Live_born_and_rate_of_birth_per_1000_ inhabitants_by_province_1975-2012/tbl0005708_i.html EUSTAT (2022b) Deaths and mortality rate of the Basque Country per 1,000 inhabitants by province. 1975-2020. Retrieved on January 21, 2022 from https://www.eustat.eus/elementos/ele0000000/defunciones-y-tasas-de-mortalidad-de-la-ca- de-euskadi-por-1000-habitantes-por-territorio-historico/tbl0000026_c.html Foy, B. H., Wahl, B., Mehta, K., Shet, A., Menon, G. I., & Britto, C. (2021). Comparing COVID-19 vaccine allocation strategies in India: A mathematical modelling study. International Journal of Infectious Diseases, 103, 431–438. https://doi.org/10.1016/J.IJID.2020.12.075 Giattino, C. (2020). How epidemiological models of COVID-19 help us estimate the true number of infections. Our World in Data. https://ourworldindata.org/covid-models Guirao, A. (2020). The Covid-19 outbreak in Spain. A simple dynamics model, some lessons, and a theoretical framework for control response. Infectious Disease Modelling, 5, 652– 669. https://doi.org/10.1016/j.idm.2020.08.010 Hu, B., Guo, H., Zhou, P., & Shi, Z. L. (2020). Characteristics of SARS-CoV-2 and COVID-19. In Nature Reviews Microbiology (Vol. 19, Issue 3, pp. 141–154). Nature Research. https://doi.org/10.1038/s41579-020-00459-7 Iboi, E. A., Ngonghala, C. N., & Gumel, A. B. (2020). Will an imperfect vaccine curtail the COVID-19 pandemic in the U.S.? Infectious Disease Modelling, 5, 510–524. https://doi.org/10.1016/j.idm.2020.07.006 Jayaweera, M., Perera, H., Gunawardana, B., & Manatunge, J. (2020). Transmission of COVID- . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 28 19 virus by droplets and aerosols: A critical review on the unresolved dichotomy. https://doi.org/10.1016/j.envres.2020.109819 Kermack, W. 0, & McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115(772), 700–721. https://doi.org/10.1098/rspa.1927.0118 Kim, J. H., Marks, F., & Clemens, J. D. (2021). Looking beyond COVID-19 vaccine phase 3 trials. In Nature Medicine (Vol. 27, Issue 2, pp. 205–211). Nature Research. https://doi.org/10.1038/s41591-021-01230-y Lauer, S. A., Grantz, K. H., Bi, Q., Jones, F. K., Zheng, Q., Meredith, H. R., Azman, A. S., Reich, N. G., & Lessler, J. (2020). The incubation period of coronavirus disease 2019 (CoVID-19) from publicly reported confirmed cases: Estimation and application. Annals of Internal Medicine, 172(9), 577–582. https://doi.org/10.7326/M20-0504 Li, B., Deng, A., Li, K., Hu, Y., Li, Z., Xiong, Q., Liu, Z., Guo, Q., Zou, L., Zhang, H., Zhang, M., Ouyang, F., Su, J., Su, W., Xu, J., Lin, H., Sun, J., Peng, J., Jiang, H., … Lu, J. (2021). Viral infection and transmission in a large well-traced outbreak caused by the Delta SARS-CoV-2 variant. https://doi.org/10.1101/2021.07.07.21260122 Li, Q., Guan, X., Wu, P., Wang, X., Zhou, L., Tong, Y., Ren, R., Leung, K. S. M., Lau, E. H. Y., Wong, J. Y., Xing, X., Xiang, N., Wu, Y., Li, C., Chen, Q., Li, D., Liu, T., Zhao, J., Liu, M., … Feng, Z. (2020). Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus–Infected Pneumonia. New England Journal of Medicine, 382(13), 1199– 1207. https://doi.org/10.1056/nejmoa2001316 Li, R., Pei, S., Chen, B., Song, Y., Zhang, T., Yang, W., & Shaman, J. (2020). Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS- CoV-2). Science, 368(6490), 489–493. https://doi.org/10.1126/science.abb3221 Lin, Q., Zhao, S., Gao, D., Lou, Y., Yang, S., Musa, S. S., Wang, M. H., Cai, Y., Wang, W., Yang, L., & He, D. (2020). A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action. International Journal of Infectious Diseases, 93, 211–216. https://doi.org/10.1016/j.ijid.2020.02.058 Liu, J., Liao, X., Qian, S., Yuan, J., Wang, F., Liu, Y., Wang, Z., Wang, F. S., Liu, L., & Zhang, Z. (2020). Community transmission of severe acute respiratory syndrome Coronavirus 2, Shenzhen, China, 2020. Emerging Infectious Diseases, 26(6), 1320–1323. https://doi.org/10.3201/eid2606.200239 López, L., & Rodó, X. (2020). The end of social confinement and COVID-19 re-emergence risk. Nature Human Behaviour, 4(7), 746–755. https://doi.org/10.1038/s41562-020-0908-8 López, L., & Rodó, X. (2021). A modified SEIR model to predict the COVID-19 outbreak in . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 29 Spain and Italy: Simulating control scenarios and multi-scale epidemics. Results in Physics, 21, 103746. https://doi.org/10.1016/j.rinp.2020.103746 MacIntyre, C. R., Costantino, V., & Trent, M. (2020). Modelling of COVID-19 vaccination strategies and herd immunity, in scenarios of limited and full vaccine supply in NSW, Australia. MedRxiv. https://doi.org/10.1101/2020.12.15.20248278 Moore, J. P. (2021). Approaches for Optimal Use of Different COVID-19 Vaccines: Issues of Viral Variants and Vaccine Efficacy. JAMA, 325(13), 1251–1252. https://doi.org/10.1001/JAMA.2021.3465 Mwalili, S., Kimathi, M., Ojiambo, V., Gathungu, D., & Mbogo, R. (2020). SEIR model for COVID-19 dynamics incorporating the environment and social distancing. BMC Research Notes, 13(1), 1–5. https://doi.org/10.1186/s13104-020-05192-1 Omer, S. B., Yildirim, I., & Forman, H. P. (2020). Herd Immunity and Implications for SARS- CoV-2 Control. In JAMA - Journal of the American Medical Association (Vol. 324, Issue 20, pp. 2095–2096). American Medical Association. https://doi.org/10.1001/jama.2020.20892 Pollán, M., Pérez-Gómez, B., Pastor-Barriuso, R., Oteo, J., Hernán, M. A., Pérez-Olmeda, M., Sanmartín, J. L., Fernández-García, A., Cruz, I., Fernández De Larrea, N., Molina, M., Rodríguez-Cabrera, F., Martín, M., Merino-Amador, P., Paniagua, J. L., Muñoz-Montalvo, J. F., Blanco, F., & Yotti, R. (2020). Prevalence of SARS-CoV-2 in Spain (ENE-COVID): a nationwide, population-based seroepidemiological study. Www.Thelancet.Com, 396, 535. https://doi.org/10.1016/S0140-6736(20)31483-5 Rǎdulescu, A., Williams, C., & Cavanagh, K. (2020). Management strategies in a SEIR-type model of COVID 19 community spread. Scientific Reports, 10(1), 1–16. https://doi.org/10.1038/s41598-020-77628-4 Randolph, H. E., & Barreiro, L. B. (2020). Herd Immunity: Understanding COVID-19. In Immunity (Vol. 52, Issue 5, pp. 737–741). Cell Press. https://doi.org/10.1016/j.immuni.2020.04.012 Rothan, H. A., & Byrareddy, S. N. (2020). The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak. Journal of Autoimmunity, 109, 102433. https://doi.org/10.1016/j.jaut.2020.102433 Rulli, M. C., D’Odorico, P., Galli, N., & Hayman, D. T. S. (2021). Land-use change and the livestock revolution increase the risk of zoonotic coronavirus transmission from rhinolophid bats. Nature Food 2021 2:6, 2(6), 409–416. https://doi.org/10.1038/s43016- 021-00285-x Sariol, A., & Perlman, S. (2020). Lessons for COVID-19 Immunity from Other Coronavirus Infections. In Immunity (Vol. 53, Issue 2, pp. 248–263). Cell Press. https://doi.org/10.1016/j.immuni.2020.07.005 . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 30 Sheng, W. H. (2020). Coronavirus disease 2019 (covid-19). Journal of Internal Medicine of Taiwan, 31(2), 61–66. https://doi.org/10.6314/JIMT.202004_31(2).01 Tang, B., Bragazzi, N. L., Li, Q., Tang, S., Xiao, Y., & Wu, J. (2020). An updated estimation of the risk of transmission of the novel coronavirus (2019-nCov). Infectious Disease Modelling, 5, 248–255. https://doi.org/10.1016/j.idm.2020.02.001 Tang, B., Xia, F., Tang, S., Bragazzi, N. L., Li, Q., Sun, X., Liang, J., Xiao, Y., & Wu, J. (2020). Erratum: The effectiveness of quarantine and isolation determine the trend of the COVID-19 epidemic in the final phase of the current outbreak in China (International Journal of Infectious Diseases (2020) 95 (288–293), (S1201971220301375), (10.1016/j.ijid.2020.03.018)). In International Journal of Infectious Diseases (Vol. 96, pp. 636–647). Elsevier B.V. https://doi.org/10.1016/j.ijid.2020.05.113 Tang. J.W, Marr, L.C, Li,Y. (2021) Dancer, S. J. Covid-19 has redefined airborne transmission. BMJ, 373 :n913 doi:10.1136/bmj.n913 Wacharapluesadee, S., Tan, C. W., Maneeorn, P., Duengkae, P., Zhu, F., Joyjinda, Y., Kaewpom, T., Chia, W. N., Ampoot, W., Lim, B. L., Worachotsueptrakun, K., Chen, V. C. W., Sirichan, N., Ruchisrisarod, C., Rodpan, A., Noradechanon, K., Phaichana, T., Jantarat, N., Thongnumchaima, B., … Wang, L. F. (2021). Evidence for SARS-CoV-2 related coronaviruses circulating in bats and pangolins in Southeast Asia. Nature Communications, 12(1), 1–9. https://doi.org/10.1038/s41467-021-21240-1 WHO. (2020), Calibrating long-term non-pharmaceutical interventions for COVID-19: principles and facilitation tools. Retrieved January 2021 from https://www.who.int/publications/i/item/WPR-DSE-2020-018 WHO. (2021a). Coronavirus Disease. Retrieved May 26, 2021, from https://www.who.int/emergencies/diseases/novel-coronavirus- 2019?gclid=CjwKCAjwlrqHBhByEiwAnLmYUHk3c3pl3RZmePd2nbVQUFdUleRglwH tHCVFuwHITtwdePAUM7nXThoC72UQAvD_BwE WHO. (2021b) Strategy to Achieve Global Covid-19 Vaccination by mid-2022 Retrieved January 20, 2022, from https://cdn.who.int/media/docs/default- source/immunization/covid-19/strategy-to-achieve-global-covid-19-vaccination-by-mid- 2022.pdf?sfvrsn=5a68433c_5 WHO (2022) Covid-19 Vaccines. Retrieved January 22 from https://www.who.int/emergencies/diseases/novel-coronavirus-2019/covid-19-vaccines Wong, W. K., Juwono, F. H., & Chua, T. H. (2021). SIR Simulation of COVID-19 Pandemic in Malaysia: Will the Vaccination Program be Effective? 1–19. Wu, J. T., Leung, K., & Leung, G. M. (2020). Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint 31 China: a modelling study. Www.Thelancet.Com, 395, 2020. https://doi.org/10.1016/S0140- 6736(20)30260-9 Yang, H. M., Junior, L. P. L., Castro, F. F. M., & Yang, A. C. (2021). Mathematical modeling of the transmission of SARS-CoV-2—Evaluating the impact of isolation in São Paulo State (Brazil) and lockdown in Spain associated with protective measures on the epidemic of CoViD-19. PLoS ONE, 16(6 June), e0252271. https://doi.org/10.1371/journal.pone.0252271 Yang, J., Tang, S., & Cheke, R. A. (2021). Impacts of varying strengths of intervention measures on secondary outbreaks of COVID-19 in two different regions. In Nonlinear Dynamics (p. 1). Springer Science and Business Media B.V. https://doi.org/10.1007/s11071-021-06294-6 Zhao, S., & Chen, H. (2020). Modeling the epidemic dynamics and control of COVID-19 outbreak in China. Quantitative Biology, 8(1), 11–19. https://doi.org/10.1007/s40484-020- 0199-0 . CC-BY-ND 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted February 9, 2022. ; https://doi.org/10.1101/2021.11.29.21266986doi: medRxiv preprint