key: cord-0968722-7cl1d663
authors: Canga, A.; Bidegain, G.
title: Modelling the effect of the interaction between vaccination and non-pharmaceutical measures on COVID-19 incidence
date: 2021-11-29
journal: nan
DOI: 10.1101/2021.11.29.21266986
sha: 02b19dfeefa8c70591ae80c0b296ce8f75defb86
doc_id: 968722
cord_uid: 7cl1d663

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, a SEIR-type continuous-time deterministic disease model was developed to explore the effect of vaccination in terms of vaccination rate and efficacy, together with varying non-pharmaceutical protection measures, on disease incidence in the initial phase of vaccination. 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 veri[fi]ed 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. The Pfizer scenario resulted in the lowest prevalence of infection and cumulative mortality, particularly for low- and medium-level protection rates. 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 model here, which is based on this example, could be easily applied to other regions or countries, modifying the strategies implemented and initial conditions.

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

To date, a variety of vaccines have been approved by the European Medicine 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) . The reported efficacy rates for vaccines requiring two doses are 94.6 and 93.6, respectively, for the Corminaty (Pfizer/BioNTech) and the COVID 19 Moderna (ModernaTX) vaccines and up to 60% for the Vaxzevria (Oxford-AstraZeneca AZD1222) vaccine (EMA, 2021) . The efficacy rate for the monodose Janssen vaccine (Johnson & Johnson) is up to 67% (EMA, 2021) .

In the last few decades, mathematical models have been used to analyse the spread of infectious diseases and explore control strategies (Chowell et al., 2016) . 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 singlepopulation continuous compartmental SEIR Kermack-McKendrick-type models, constructed using ordinary differential equation (ODE) systems (R. Wu et al., 2020) . 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. 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 this this version posted November 29, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 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 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 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 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 an extension of a Kermack-McKendrick-type model (Kermack & McKendrick, 1927) .

Its main objective is to explore the impact of the interaction between different vaccination scenarios and levels of protection measures on disease incidence in the initial phase of the COVID-19 vaccination (i.e. when around 30% of the population is vaccinated). For this, as a first step, the model is verified and evaluated on the response timeline of the first and second waves of the pandemic in the Basque Country (N. Spain) using data from the Basque Health Department (BHD, 2020). This example represents a region with a very high COVID-19 incidence, since Spain has been one of the countries most affected by the pandemic (Guirao, 2020) becoming the country with the sixth highest excess mortality in the world (Dong et al., 2020) , and the Basque Country has been for many months, one of the regions reporting the highest numbers of positive cases and deaths, within the country. The first mass vaccination programme in the Basque Country started on December 27, 2020, reporting by then 117000 cases and 3030 fatalities (BHD, 2020) . This model could be easily applied to other regions or countries with different disease control strategies and initial conditions.

The model is a one-population deterministic compartmental model, continuous in time, unstructured in spatial or age terms, and configured to simulate the dynamics of COVID-19 . CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 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 variables together with the associated parameters satisfy a system of ODEs describing the dynamics of the host-pathogen system.

The model here is an adapted SEIR model with 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); 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) . 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 . CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint 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.

Description and modelling assumptions Initial conditions The model specifies an open population where birth of new susceptible individuals is a function of the total population N. For this analysis, and for such a rapid spread of disease 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) .

One of the features 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) . There is also a group of asymptomatic patients, who are infected and able to infect, but not under quarantine. Because of these factors it is very difficult to measure the real dynamic of the virus (Zhao & Chen, 2020). The model here assumes as valid the result of the seroprevalence study carried out in Spain where around 33% of the cases were asymptomatic (Pollán et al., 2020) .

The quarantine rate was estimated considering this result as a rate for the non-quarantined and not confirmed infectious cases. 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 this this version posted November 29, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 In general, the model contemplates a set of specific assumptions, such as some protection measures defined by the parameters in Table 2 . 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.

The subpopulations of the model satisfy a system of ODEs describing the dynamics of the hostvirus 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:

Equation (1) 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 this this version posted November 29, 2021. ; https://doi.org/10. 1101 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 background mortality or the end of the incubation period in exposed individuals.

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.

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. 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.003-0.03 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) . CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10. 1101 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.

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.

Equation (7): The change in the number of vaccinated and protected individuals 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.

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.

From December 27, 2020, 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 previously mentioned efficacy rates 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 vaccine, the efficacy is achieved after 14 days (EMA, 2021).

Since the start of the vaccination programme on the 27 th of December, up until the 8 th 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.16% Pfizer, 11.21% Moderna, 2.96% Astra Zeneca and 8.67% Janssen (BHD, 2020), . CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint Table 3 . Estimated vaccination and protection rates using the Basque Countries' vaccination data (BHD, 2020) The average number of vaccines administered per day (vaccination rate), 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:

where T is 160 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).

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 license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint

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) 

Considering the degree of protection of each vaccine in terms of the estimated vaccination protection rate, the effects of five different vaccination strategies were tested within the model.

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 .

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.

The curve of new daily confirmed cases by the BHD (BHD, 2020) ( Figure 1 ) is described below in terms of the non-pharmaceutical interventions:

. CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10. 1101 On February 28, 2020 the first case of COVID-19 was confirmed in the Basque Country. On

March 13, schools were closed in all regions of Spain including the Basque Country and on March 14 a nationwide lockdown (confinement) was declared, banning all public events.

During this period the first wave reached its peak on the March 24, with 723 confirmed new cases. At this point it must be considered that tests were conducted only for those who presented symptoms such as fever and cough. Thus, people who did not seek medical attention were tested very rarely.

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". However, on July 8, mask use in public areas became mandatory for people older than six years old to limit the virus transmission. By mid-July, daily cases started to increase consistently, and the second wave reached its peak on August 28, with 886 confirmed new cases. This peak was reached when mobility all around the country was permitted, international travellers entered the country without restrictions, and indoor public and private facilities were open. This wave was linked to the increase of people gathering for leisure purposes, coinciding with summer holidays (BHD, 2020).

The peak of the third wave reached on November 5, with 1547 new confirmed cases. This wave could be linked to the change of the season, the return to schools and work, and the 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint increase of indoor activities. 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.

A fourth wave peak was recorded on January 27, 2021, with 1274 new cases. This wave was linked with Christmas holidays despite the limitations imposed by the government, limiting gatherings of people and mobility in the territory.

The fifth wave started in mid-March, reaching its peak on April 21, with a maximum of 1013 new cases. This wave was linked with the easing of mobility limitations inside the Basque territory to promote tourism during the Easter holidays.

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.

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 (0.005 day -1 , before lockdown) to a high protection level (0.025 day -1 , during lockdown) (see Table 2 ).

. CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint 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 (0.002 day -1 ) to high (0.006 day -1 ) in Case 1 results 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 morality rate (0.006 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 (0.57 day -1 ) to high (0.67 day -1 ) in Case 2 results 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). 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 this this version posted November 29, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 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 expectations regarding the number of recovered individuals (around 21000, approximately estimated casesconfirmed 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).

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 lowmedium-level protection rate of 0.065 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 (0.001 day -1 ) to high (0.003 day -1 ) in Case 5 results 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 (0.001 day -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 (0.57 day -1 ) to high (0.67 day -1 ) in Case 6 results 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 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).

. CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 (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).

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 lowlevel to high-level protection rates as in Table 2 : (i) the low-level protection rate  was estimated to be 0.005 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 license It 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 this this version posted November 29, 2021. ; https://doi.org/10. 1101 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 = day -1 represented the previous scenario with no regional and national mobility restrictions, and (iii) the high-level protection rate of = 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 (B. 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.

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 =0.005 to =0.01, 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 license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint (Table 1) and parameter values (Table 2) 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 this this version posted November 29, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 The highest number of new cases occurs with the lowest protection rate (=0.005 day -1 ) where fewer restrictions are applied ( Figure 5A ). This number reduces to half when the protection measures increase to = 0.0075 day -1 , decreasing even more if the protection rate is set to =0.01 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.

In this simulation ( Figure 6A ), Pfizer is the unique vaccine administered to the population, with a vaccine protection rate of 0.0017 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 license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 

In this simulation ( Figure 6B ), Moderna is the unique vaccine administered to the population with a vaccine protection rate of 0.0016 day -1 (Table 4 ). The model estimates 56,500 cases in the simulation period and 458 deaths with the lowest protection rate. The effect of the high protection scenario (green line) on disease dynamics is similar to that for the Pfizer scenario, with a higher impact in cumulative mortality, which decreases from 458 to 87 deaths.

This simulation represents a scenario with the Astra Zeneca vaccine as the unique vaccine administered ( Figure 6C) , with a vaccine protection rate of 0.001 day -1 (Table 4 ). Here, for the medium and low protection level, the numbers of cases and deaths are notably higher compared to the Pfizer and Moderna vaccines. Nevertheless, for the higher protection rate, differences between the responses to vaccines are not significant.

. CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint

This simulation ( Figure 6D ) shows the Janssen vaccine as the unique vaccine administered with a vaccine protection rate of 0.0011 day -1 (Table 4) . Similarly to the Astra Zeneca scenario, 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 is less effective in reducing cases and deaths compared to Pfizer and Astra Zeneca.

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 selfprotection, 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 projections explored here for model verification and validation do not differ much from the reported cumulative mortality data and are consistent with 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 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) .

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 . CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10. 1101 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 Pfizer vaccine scenario is the one that maintains daily positive cases and deaths at the lowest levels even when the protection rate is low. Nevertheless, 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 ). The limited vaccine scenario ( Figure 5 ) and the Moderna scenario are very similar in terms of disease incidence, regardless of the protection level. The Astra Zeneca and Janssen vaccine scenarios are the ones with the worst results in COVID-19 cases and cumulative mortality for medium and low protection levels.

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 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) . 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint . CC-BY-ND 4.0 International license It 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 this this version posted November 29, 2021. ; https://doi.org/10.1101/2021.11.29.21266986 doi: medRxiv preprint . CC-BY-ND 4.0 International license It 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 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 this this version posted November 29, 2021. ; https://doi.org/10. 1101 

These values are obtained by model fitting against real cumulative mortality data (BHD, 2020) and considering that at least 33% of the cases are not confirmed

Effects of ventilation on the indoor spread of COVID-19

Evolución del coronavirus (COVID-19) en Euskadi -Conjunto de datos de Open Data Euskadi -Euskadi.eus

A seir model for control of infectious diseases with constraints

Evolutionary origins of the SARS-CoV-2 sarbecovirus lineage responsible for the COVID-19 pandemic

Symptoms of COVID-19 | CDC

Superspreading Event of SARS-CoV-2 Infection at a Bar

Mathematical models to characterize early epidemic growth: A review

Dynamics of population immunity due to the herd effect in the COVID-19 pandemic

WHO declares COVID-19 a pandemic

Immunological memory to SARS-CoV-2 assessed for up to 8 months after infection. and meta-analysis

Mathematical epidemiology of infectious diseases: model building, analysis

An interactive web-based dashboard to track COVID-19 in real time

COVID-19 situation update worldwide, as of week 4

COVID-19 vaccines: authorised | European Medicines Agency

The role of environmental factors to transmission of SARS-CoV-2 (COVID-19)

Comparing COVID-19 vaccine allocation strategies in India: A mathematical modelling study

How epidemiological models of COVID-19 help us estimate the true number of infections

The Covid-19 outbreak in Spain. A simple dynamics model, some lessons, and a theoretical framework for control response

Characteristics of SARS-CoV-2 and COVID-19

Will an imperfect vaccine curtail the COVID-19 pandemic in the

A contribution to the mathematical theory of epidemics

Looking beyond COVID-19 vaccine phase 3 trials

The incubation period of coronavirus disease 2019 (CoVID-19) from publicly reported confirmed cases: Estimation and application

Viral infection and transmission in a large well-traced outbreak caused by the Delta SARS-CoV-2 variant

Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus-Infected Pneumonia

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

Community transmission of severe acute respiratory syndrome Coronavirus 2

The end of social confinement and COVID-19 re-emergence Spain and Italy: Simulating control scenarios and multi-scale epidemics

Modelling of COVID-19 vaccination strategies and herd immunity, in scenarios of limited and full vaccine supply in

Approaches for Optimal Use of Different COVID-19 Vaccines: Issues of Viral Variants and Vaccine Efficacy

SEIR model for COVID-19 dynamics incorporating the environment and social distancing

Herd Immunity and Implications for SARS-CoV-2 Control

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

Management strategies in a SEIR-type model of COVID 19 community spread

Herd Immunity: Understanding COVID-19

The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak

Land-use change and the livestock revolution increase the risk of zoonotic coronavirus transmission from rhinolophid bats

Lessons for COVID-19 Immunity from Other Coronavirus Infections

Coronavirus disease 2019 (covid-19)

An updated estimation of the risk of transmission of the novel coronavirus (2019-nCov)

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

Covid-19 has redefined airborne transmission

Evidence for SARS-CoV-2 related coronaviruses circulating in bats and pangolins in Southeast Asia

Coronavirus Disease

SIR Simulation of COVID-19 Pandemic in Malaysia: Will the Vaccination Program be Effective?

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

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

Impacts of varying strengths of intervention