Abstract
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 nonpharmaceutical 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 types of vaccines (Pfizer, Moderna, Astra Zeneca, and 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. Significant differences in disease prevalence and cumulative mortality were found between vaccination scenarios for low and medium-level protection measures. For high-level protection measures, any vaccine scenario is effective at limiting the virus transmission and disease mortality. 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.
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 etiology [1–3]. 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 [3–5]. After more than a year and a half since it was declared a pandemic, according to the WHO, there are 200 M cases reported and 4 M deaths worldwide [6].
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 [7], (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 [8–10]. The virus transmission in indoor settings has been the main transmission pathway when ventilation is not sufficient [10–13]. 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 [14], being these events principal drivers determining the dynamics of COVID-19 transmission [15].
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 [11]; CDC [16–20].
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 [21, 22]. 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 [23]. 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 [24, 25]. Herd immunity originates when a sufficiently large proportion of the population is immune to the disease [24, 26]. 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 [27, 28]. 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 [29]. Consequently, initial vaccination phases such as the one studied here (around 30% of the population vaccinated) may require the maintenance of some level of nonpharmaceutical 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 [30]. Most models for the COVID-19 pandemic are single-population continuous compartmental SEIR Kermack-McKendrick-type models, constructed using ordinary differential equation (ODE) systems [2, 31–33]. 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 compartments 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 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 nonreported positive cases, ranging between 10 and 70% of the total, that correspond to people that do not get tested or are asymptomatic to the disease. Hence, models such as 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 [34]. Given the uncertainty surrounding the situation after COVID-19 vaccination programmes, models estimating these unconfirmed cases, such as the ones 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 vaccination rate and efficacy, and different levels of nonpharmaceutical 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 the 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 [35]. It is a deterministic SEIR transmission model that accounts for important characteristics 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 [36]. 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 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 [37, 38]); The I subpopulation is assumed to represent asymptomatic cases, nonconfirmed, and nonisolated symptomatic cases, and cases that are not yet 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 [38], they are most infectious in the early stages of a symptomatic stage, so the majority of infections are caused by symptomatic individuals [20], (4) a fraction of this I subpopulation is the pool of infected individuals, representing 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 2,199,711 individuals based on demographic data of the Basque Country Institute of Statistics (BIS) [39]. 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 [40] and considering that at least 33% of the cases are asymptomatic, not confirmed and able to infect, but not under quarantine [41].
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 [38]. 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:
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 background 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 (ii) the loss of individuals because of background 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 (i) the gain of individuals due to the end of incubation period in exposed individuals and (ii) the 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 (i) 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.
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.
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
2.3.1. 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 [40]. 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.
2.3.2. Vaccination Data
The first mass vaccination programme in the Basque Country started on December 27, 2020, reporting by then 117000 cases and 3030 fatalities [40]. 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 [23] in Tables 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 on 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 [50].
Since the start of the vaccination programme on December 27 up until June 8, 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 [40].
The average number of vaccines administered per day (vaccination rate, VR) was estimated from data obtained from the vaccination bulletin of the BHD [40]. 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 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 nonpharmaceutical 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.
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 programme was run under different conditions of disease mortality and quarantine rate. This analysis was 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 interventions (low, medium, and high protection levels (Table 1)), the effects of five different vaccination strategies were tested within the model.
2.6.1. 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.
2.6.2. 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 disease incidence; since infringement of low-level protection restrictions such as the use of mask or social distancing may be much more common than strong restrictions 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 [46, 48]; however, this assumption is overestimating the S population.
A seroprevalence study carried out in Spain estimated 33% of the cases as asymptomatic [41]. 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 nonconfirmed cases and hence, nonquarantined people. This assumption could be underestimating the quarantined population particularly beyond the third wave when testing was common to close contacts of COVID-19.
For scenarios with vaccine full supply caution is required to interpret the modelling results (see Section 3) quantitatively; since the model is dealing with multiple dimensions, latent covariates, a local region without %100 certain data and four vaccine types.
3. Results
3.1. Descriptive Analysis of the Curve of Infection Cases
The curve of new daily confirmed cases by the BHD [40] (Figure 1) is described below in terms of the nonpharmaceutical interventions’ chronological sequence. This detailed analysis is essential for the selection of simulation cases with different protection measures and later discussion.
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. 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 was 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 11 pm to 6 am. 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. The fifth wave started in mid-March, reaching its peak on April 21, with a maximum of 1013 new cases.
3.2. 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
3.2.1. First Wave
Realistic scenarios with varying disease mortality (Figure 2(a)) and quarantine rate (Figure 2(b)) 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).
Increasing the disease mortality rate from low (2 × 10−3 day−1) to high (6 × 10−3 day−1) in Case 1 results in an expected increased cumulative mortality conforming to expectations (Figure 2(a)). 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 that the model has the best fit to the 1687 deaths confirmed by the BHD in the first wave [40].
Increasing the quarantine rate from low (5.7 × 10−1 day−1) to high (6.7 × 10−1 day−1) in Case 2 results in an expected reduced cumulative mortality (Figure 2(b)) 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 [40].
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 2(c)). The behaviour of the model in these terms conforms to 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 [40].
Finally, Case 4 evaluated the model against new confirmed daily cases (Figure 2(d)). The model estimates 22230 cases as the number of true infections for the simulation period; about twice as high as the 13862 confirmed cases [40].
3.2.2. Second Wave
Similar to the first wave, realistic scenarios with varying disease mortality (Figure 3(a)) and quarantine rate (Figure 2(a)) 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 [40]. Increasing disease mortality rate from low (1 × 10−3 day−1) to high (3 × 10−3 day−1) in Case 5 results in an expected increased cumulative mortality as expected (Figure 3(a)). For this case, model evaluation against real cumulative data shows that the simulation with the lower disease morality rate (1 × 10−3 day−1 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 [40].
Increasing the quarantine rate from low (5.7 × 10−1 day−1) to high (6.7 × 10−1 day−1) in Case 6 results in an expected reduced cumulative mortality (Figure 3(b)) as in Case 2, confirming the adequate behaviour of the model. In 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 [40].
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 3(c)). 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) [40]. Finally, in Case 8, the model estimated 45300 cases as the number of true infections for the simulation period (Figure 3(d)); higher than the 31000 confirmed cases [40].
3.3. 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 nonpharmaceutical 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 inside the Basque territory and eventually all around Spain to promote tourism during the Easter holidays and before summer, (ii) 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 of 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 [51]. Initial conditions are those in Table 1 with changes in susceptible (S = 2000000), infected (I = 200), and vaccinated (V = 52500) subpopulations due to the course of the disease dynamics.
3.3.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 4(a)). 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 [40].
(a)
(b)
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.
The highest number of new cases occurs with the lowest protection rate (α = 5 × 10−3 day−1) where fewer restrictions are applied (Figure 5(a)). 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 5(b)).
(a)
(b)
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.
3.3.2. Simulation
(2) Pfizer Scenario. In this simulation (Figure 6(a)), 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). Similar to simulation 1, an increased protection rate reduces the incidence of cases and mortality (Figure 6(a)), 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.
(a)
(b)
(c)
(d)
3.3.3. Simulation
(3) Moderna Scenario. In this simulation (Figure 6(b)), 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.
3.3.4. Simulation
(4) Astra Zeneca Scenario. This simulation represents a scenario with the Astra Zeneca vaccine as the unique vaccine administered (Figure 6(c)), 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.
3.3.5. Simulation
(5) Janssen Scenario. This simulation (Figure 6(d)) 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.
4. 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 nonpharmaceutical protection measures such as self-protection, mobility restrictions, closure of indoor facilities, and lockdown. The Kermack and McKendrick [35] 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 nonpharmaceutical 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 nonpharmaceutical 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 in November linked with the return to schools and work, and the increase of indoor activities. A fourth wave peak was recorded in 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.
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 [52]. However, estimated new daily cases are significantly higher than those reported by the BHD [40]. This result is consistent with the fact that confirmed cases may be undercounted [34] 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, uses 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 [34]. This is also consistent with the seroprevalence study carried out in Spain, where around 33% of the infected cases were asymptomatic [41]. On the other hand, there is a deviation ratio and as expected if the number of cases increases 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 nonpharmaceutical 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 COVID-19 is the one combining vaccination and high protection levels for nonpharmaceutical 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,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 nonpharmaceutical 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 [40]. 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 nonpharmaceutical 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 nonpharmaceutical 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 [53]. 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 [54], same as considering the new variants of the virus that could directly impact on the vaccine-induced protection and the transmission rate [55].
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 minimize the appearance of severe cases of disease that may collapse health systems.
Data Availability
The data used to support the findings of this study were obtained from the open data system on the Basque Health Department (Osakidetza) website [40].
Disclosure
The contents of this manuscript do not necessarily reflect the point of view of the BHD and in no way presents the BHD´s future policy in this area. The preprint of this article was posted on medRxiv (Canga et al. [56]).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
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). The authors appreciate this support. The manuscript benefited from helpful discussion with Dr. Aitana Lertxundi and Dr. Naroa Kajarabille (Department of Preventive Medicine and Public Health, UPV/EHU), 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).