Abstract
The promise of efficacious vaccines against SARS-CoV-2 is fulfilled and vaccination campaigns have started worldwide. However, the fight against the pandemic is far from over. Here, we propose an age-structured compartmental model to study the interplay of disease transmission, vaccines rollout, and behavioural dynamics. We investigate, via in-silico simulations, individual and societal behavioural changes, possibly induced by the start of the vaccination campaigns, and manifested as a relaxation in the adoption of non-pharmaceutical interventions. We explore different vaccine efficacy, vaccination rollout speeds, prioritization strategies, as well as multiple behavioural responses. We apply our model to six countries worldwide (Egypt, Peru, Serbia, Ukraine, Canada, and Italy) selected to sample diverse socio-demographic and socio-economic contexts. To isolate the effects of age-structures and contacts patterns from the particular pandemic history of each location, we first study the model considering the same hypothetical initial epidemic scenario in all countries. We then calibrate the model using real epidemiological and mobility data for the different countries. Our findings suggest that early relaxation of safe behaviours can jeopardize the benefits brought by the vaccine in the short term: a fast vaccine distribution and policies aimed at keeping high compliance of individual safe behaviours are key to mitigate disease resurgence.
Introduction
The COVID-19 pandemic has been largely fought with non-pharmaceutical interventions (NPIs). Bans of events and social gatherings, limitations in national and international travels, school closures, shifts towards remote working, curfews, closure of pubs and restaurants, cordon sanitaires, national and regional lockdowns are examples of governmental interventions implemented around the world to curb the spreading of SARS-CoV-2 [1–5]. While extremely effective, such top-down NPIs induce profound behavioural changes, bring many social activities to a halt, and thus have huge socio-economic costs. Hence, alongside these measures, governments nudged and/or mandated populations to adopt another set of NPIs. Social distancing, face masks, and increased hygiene are examples [6]. Although far from being cost-free, they are more feasible, sustainable, and allow for higher levels of socio-economic activity. As such, they have been the leitmotif of the post first COVID-19 wave in many countries. Unfortunately, awareness, adoption, and compliance with these NPIs have been spotty [6]. Furthermore, they have not been complemented with sufficiently aggressive test and trace programs. As result, many countries experienced marked disease resurgences after the summer 2020 and some had to resort to new lockdowns [7]. As we write, we have turned a crucial corner in hampering the resurgence, diffusion and severe outcomes of the disease. Several vaccines have shown great results of their phase 3 trials, and have been authorized for emergency use by several regulating agencies and tens of others are in the pipeline [8–10]. As result, we have witnessed the start of vaccination campaigns around the world. However, the logistical issues linked to the production, delivery, and administration of billions of doses on a global scale impose caution when evaluating the impact vaccines will have on the pandemic in the short and medium term. They will be a scarce resource and it will take time to vaccinate the fraction of the population necessary for herd immunity [11, 12]. Furthermore, vaccines are not perfect and there are many other unknowns [13]. The results show high levels, around 95% [14, 15], of direct protection against the disease. Despite the first signs of the vaccine impact in the real world are extremely encouraging [16, 17], these figures might end up to be lower, it is still unclear to what extent vaccines limit further transmission, and how long the immunity lasts [13]. Finally, vaccines’ acceptance is a complex challenge. A recent survey among 13, 426 participants in 19 countries shows that, while 71.5% of the sample is very or somewhat likely to take the vaccine, there are large heterogeneities [18]. Acceptance rates vary from 90% in China to less than 55% in Russia. Furthermore, they are linked to socio-economic, socio-demographic features, and education attainment [18]. Arguably, vaccines alone will not be able to contain the spreading of the virus, at least in the short term [19]. Social distancing, face masks, hygiene measures, and other NPIs will be still key during the delivery of such vaccination programs, especially considering the emergence and fast spread of variants of concern [20].
In this context, an important question emerges. What will happen to adoption and compliance to NPIs as vaccination campaigns progress? Their arrival and delivery might induce individual and collective behavioural changes. Some might see this milestone as the official end of the emergency and as result relax their COVID-safe behaviours. Somehow paradoxically, vaccines might have, at least initially, a net negative effect. According to the health-belief model, one of the most commonly used psychological theories to characterize health-related behaviours, beliefs, perceptions, barriers to take action, and other modifying variables such as socio-demographic and socio-economic factors are key ingredients driving behavioural changes [21–23]. A recent study may offer an empirical evidence of this tendency [24]. By means of surveys delivered before the vaccination rollout in December 2020, authors found that vaccine information reduce adherence to social distancing, hygiene measures, and the willingness to stay at home. Several surveys conducted during the COVID-19 pandemic, well before any concrete hope for a vaccine, confirm this picture, provide hints of how the arrival of vaccines might corrode even more adoption, highlight how compliance is a complex multi-faced problem [5] and that risk-perception as well as NPIs adoption are indeed associated to several socio-economic determinants such as age, gender, wealth, urban-rural divide [25–40].
The literature aimed at estimating the epidemiological and societal impact of COVID-19 vaccines has been focused mainly on two very important points. The first line of research has been devoted to quantifying the effects of a vaccine on the evolution of the pandemic, considering different efficacy and coverage levels [41]. The second instead tackled the issue of vaccine allocation investigating strategies that target first different groups (i.e., age brackets, high-risk individuals) or particular occupations (i.e., doctors, nurses) [10, 42–44]. Very recently, the intuition that social distancing remains key during vaccination rollout stimulated few studies on the effects of a vaccine on the adoption of NPIs in specific settings [45–48], but the impact on real-world scenarios is still largely unexplored.
To tackle such limitation, we introduce an age-structured compartmental epidemic model capturing the possible relaxation of NPIs adoption in response to the vaccine rollout. We model different compliance levels as distinct compartments and consider different behavioural dynamics driving the relaxation of NPIs. We test extensively the effects of behaviour change on disease spreading for different prioritization strategies, vaccine efficacy, and vaccination rollout speeds, using real demographic data and contacts matrices for six countries: Egypt, Peru, Serbia, Ukraine, Canada, and Italy. We choose these countries sampling levels of economic development. Indeed, in the World Economic Situation and Prospects 2020 issued by the United Nations, Egypt and Peru are classified as Developing Economies [49], Serbia and Ukraine as Economies in Transition, Canada and Italy as Developed Economies. Furthermore, considering dissimilar countries allows us to explore the interplay between vaccination and behaviours also as a function of population pyramids and intra/inter-generational mixing observed around the world. High income countries are typically characterized by higher average age, but lower inter-generational interactions respects to mid/low income countries [50]. These observations, together with the dependence of COVID-19 fatality rates on age, point to the possibility of non-trivial dependencies which we aim to explore here. As a way to realistically account for the different epidemic trajectories, we also explore the model after calibrating it on COVID-19 weekly deaths in the period 2020/09/01 − 2020/12/31 in these countries. To this end, we incorporate the timeline and effects of governmental restrictions on social contacts.
Our results provide quantitative insights on the interaction between sustained NPIs and an effective vaccination program. We show that an early relaxation of COVID-safe behaviours may lower, and even nullify, the advantages brought by the vaccine in the short term. Overall, the picture that emerges from the analysis of the different countries is consistent: a high level of compliance towards NPIs such as mask-wearing, social distancing, and avoidance of large gatherings, will be needed in the next months in order to avoid spoiling the great effort of the vaccination campaigns.
Results
Vaccine-Behaviour Model
We consider an age-structured epidemic model based on a Susceptible-Latent-Infectious-Recovered compartmentalization with the addition of pre-symptomatic and asymptomatic infection stages. Susceptible individuals (S), interacting with the infectious compartments, may enter the Latent stage (L). Disease progression let latent individuals evolve in the pre-symptomatic condition (P) and afterwards either in the asymptomatic (A) or symptomatic (I) stage. Both A and I individuals finally recover and enter the Recovered compartments (RA, RI). We also compute the number of daily deaths according to the age stratified Infection Fatality Rate (IFR) on the symptomatic infections. The individuals that are responsible for the new infections are I, A, and P, with the latest two characterized by a lower infectiousness. On top of the disease dynamics, we model both the vaccination rollout and the behavioural responses linked to it. After the start of the vaccination rollout, on each day a fraction of the susceptible population receives a vaccine that reduces the probability of infection proportional to vaccine efficacy V E and enters in the compartment V. We introduce the rollout speed rV as the number of daily administered vaccine doses expressed as a percentage of the total population. We study and compare three different strategies of vaccine prioritization. The first considers distributing vaccines in decreasing order of age. Previous research on COVID-19 vaccines allocation has shown that this strategy is the most effective in reducing the number of deaths and overall severity [10, 42, 43]. It is also the main strategy currently deployed in the different vaccination campaigns across the globe. In the second and third strategy vaccines are either distributed homogeneously or first to the individuals in age brackets 20-49 [51] and then to the rest of the population. The last strategy seeks to reduce symptomatic transmission targeting the most social active part of the population [10]. For simplicity, we will refer to the three rollout approaches as vaccine strategy 1, vaccine strategy 2, and vaccine strategy 3. In parallel to the vaccination, individuals (both vaccinated and not) may start giving up safe behaviours and expose themselves to higher infection risks. We introduce the compartments SNC and V NC for susceptible and vaccinated individuals respectively, to describe different behavioural classes (for convenience “NC” refers to “non-compliant” to COVID-safe behaviours). The increased infection risk for these individuals is captured by the parameter r > 1 (for example, r = 1.3 indicates an increased risk of 30%). The model accounts also for possible transitions towards and from the non-compliance compartments. In particular, we model the transition from S and V towards riskier behavioural classes as a function of the fraction of the vaccinated population and the parameter α. In turn, we imagine that a worsening of the epidemiological conditions may push non-compliant individuals back to safer behaviours. Here, we consider the number of fatalities per 100, 000 individuals in the previous time step (i.e., day) and the parameter γ to control the second behavioural transition. This modeling choice aims to depict the adaptive nature of human behaviour where individual choices are influenced by vaccination and pandemic progression [52]. We refer the reader to the Materials and Methods section and the Supplementary Information for more details on the model. In the Supplementary Information, we also test another (simpler) mechanism where individuals relax their behaviours (i.e., transit to SNC and V NC) at a constant rate α. Similarly, non-compliant embrace COVID-safe behaviours again (i.e., transit back to S and V) at a constant rate γ. The overall picture emerging adopting this mechanism does not change significantly.
Interplay between NPIs adoption and vaccination campaign
We use data for six different countries: Egypt, Peru, Serbia, Ukraine, Canada, and Italy. As mentioned above, these countries have been selected to sample a range of demographics and socio-economic contexts. In Figure 1A, we represent some key characteristics of the demographic and the mixing patterns between age groups of these countries. First, we consider the fraction of people aged over 65. This is a measure of the epidemic fragility of a country: indeed 65+ individuals are at particularly high risk of mortality from COVID-19 disease [53]. Italy is the country showing the highest fraction of 65+ people (24%), while the developing economies show the lowest fraction (Egypt 5.3% and Peru 8.7%). Indeed, low/middle-income countries tend to have a younger population with respect to the high-income ones. Second, we present also the contacts intensity of different age groups. This is defined as the total number of contacts that an individual in a certain age group has, on average, with all the others in a day. We observe a typical decreasing trend, with younger people that tend to have more contacts. In particular, developing economies show a much higher number of daily contacts for individuals aged under 30. This can play an important role in the spreading even in developed economies. According to a recent study the resurgence of the COVID-19 Pandemic in the US after Summer 2020 was mainly sustained by younger people [51]. Finally, we represent a measure of inter-generational mixing. We define it simply as the number of daily contacts that an individual in the age groups at high mortality risk from COVID-19 (65+) receives from individuals in the age brackets 0 − 49. We observe that Egypt is the country showing the highest inter-generational mixing, followed closely by Peru.
It is important to notice how heterogeneities in health infrastructures, access to health care, and comorbities are expected across the six countries under examination [50]. These features might induce variations in the IFR. For simplicity however, in the following, we use the same fatality rates for all countries. We leave the exploration of this important aspect to future work.
Prioritization strategies and vaccine efficacy
The proposed model allows to investigate the impact of the behavioural response under different conditions. We consider six populations matching the characteristics of the countries under consideration and explore the model with the same epidemic initial conditions. The experimental setup considers that each population has already experienced a previous wave of infections and that restrictive measures are in place to mitigate the spreading. Therefore, we set the basic reproductive number R0 = 1.15, 1% of initially infected individuals, and 10% of immune individuals. In Figure 1B we show the relative deaths difference for the three vaccine prioritization strategies and three vaccine efficacy V E = 50%, 70%, 90%. The relative deaths difference is defined as the fraction of deaths averted thanks to the vaccine with respect to a baseline simulation without vaccine and, therefore, no behavioural response. Moreover, we explore a range of behavioural responses running the simulations for a range of values of the parameter α. Starting from α = 0 (which implies no behavioural response), we perform simulations with increasing α values (implying stronger reactions), while keeping constant the other behavioural parameter (γ = 0.5).
As a first observation, we note that, across the different countries and vaccine efficacy considered, the strategy aiming to reduce the severity of the pandemic (i.e., vaccine strategy 1) is indeed the most effective in reducing the number of deaths, followed by the homogeneous (i.e., vaccine strategy 2) and then the strategy targeting a reduction in transmission (i.e., vaccine strategy 3). This is consistent with the findings of previous studies focused on the allocation of vaccines among groups [42, 43]. As illustrative and concrete example let us consider the case of Canada. When α = 0 (i.e., no behaviour response) and V E = 90% with this strategy the fraction of averted deaths is 0.68 with respect to a baseline without vaccine. This fraction reduces to 0.62 with the homogeneous, and to 0.59 with the strategy prioritizing younger individuals. Furthermore, the difference among strategies diminishes when lower vaccine efficacy are considered. In the previous example, when V E = 70% the fractions of averted deaths for the three strategies are 0.62, 0.58, 0.57, and when V E = 50% are 0.54, 0.53, 0.52.
When the behavioural response to the vaccine rollout (i.e., α > 0) is considered, we note a consistent decreasing trend of the relative deaths difference for increasing values of α. This shows that the behavioural response impacts the unfolding of the epidemic, and that a relaxation of NPIs leads to a smaller fraction of averted deaths. Interestingly, a concerning effect also emerges: in some conditions, as non-compliance becomes larger, the benefit brought by the vaccine is nullified and the number of observed deaths increases with respect to the no-vaccine and no-behaviour change scenario (i.e., the relative deaths difference become negative). This is solely attributable to the behavioural reaction to the vaccination campaign which in turn is not efficient enough to balance behaviour relaxation. Indeed, this phenomenon is observed in particular when strategies that do not target a reduction in severity are employed, and when the vaccine efficacy is low. In the case of Serbia, with V E = 90% and strategy aimed at reducing severity, the fraction of averted deaths goes from 0.77 when α = 0 to 0.56 when α = 10, with a potential loss of 0.21 attributable to the NPIs relaxation. This effect is much more pronounced in the case of the vaccination strategy 2 (or 3): the fraction of averted deaths goes in this case from 0.67 (0.63) when α = 0 to 0.34 (0.20) when α = 10, with a potential loss of 0.33 (0.43) attributable to the NPIs relaxation. Similarly, lower vaccine efficacy are impacted more significantly by the behaviour. In the previous example, vaccines with an efficacy of 90%, 70%, and 50% would avert, in the case of the vaccination strategy 1, respectively 0.77, 0.69, 0.61 of deaths when α = 0. When instead α = 10, these figures drop to 0.56, 0.24, − 0.19 with a drop of 0.21 in the first, of 0.45 in the second, and of 0.80 in the third case.
In the Supplementary Information we repeat the analysis considering the fraction of averted symptomatic infections instead of deaths. We find that the ranking of allocation strategies is inverted. Indeed, when considering averted infections the most efficient strategy is the one targeting the younger population first, followed by the homogeneous, and the strategy aimed at curbing severity. This is in line with previous findings [42]. Similarly to the case of deaths, also in this additional analysis we find that early behavioral relaxation reduces the fraction of averted cases, influences heterogeneously the prioritization strategies, and impacts more significantly lower vaccine efficacy.
In the results presented in Figure 1B, we observe that the most disadvantaged countries are the developing economies, Egypt and Peru. At first, this may seem counter-intuitive. Indeed, these countries have a lower fraction of individuals aged over 65. A possible explanation is given by the high activity of young individuals combined with the high inter-generational mixing. Moreover, it is worth stressing that these results are to be intended in relative terms: a relative worst performance in averting deaths does not necessarily imply a worst absolute performance. We observe also that the difference between countries is more profound for the vaccination strategy aimed at reducing transmission. Furthermore, behavioural relaxation widens the distance between the countries. With V E = 90% and the rollout strategy targeting a reduction in transmission, the gap between the relative deaths difference between Italy and Egypt is 0.22 when α = 0. This figure increases to 0.53 when α = 10.
Vaccine rollout speed
In Figure 2, we investigate the interplay between NPIs relaxation and vaccine rollout speed, keeping the vaccine efficacy constant (V E = 90%) and varying rV between 0.1% and 1%. We choose these values to cover the spectrum of real vaccination rollout speeds of the vaccination campaigns across the globe. Peru, for example, administered on average 0.05 daily doses per 100 people in the week commencing on the 8th of March, Italy administered 0.30, and Serbia 0.70 [54]. We observe that, besides allocation strategy and vaccine efficacy, rollout speed is a fundamental factor of the vaccination campaign. In the case of Egypt and a vaccination strategy that targets a reduction of severity, when α = 0, the fraction of averted deaths spans between 0.54 with the fastest rollout considered (rV = 1%) and 0.23 with the slowest one (rV = 0.1%). This implies a drop of 0.31. When, instead, we consider also behavioural response, this gap widens significantly. Indeed, for α = 10, we go from a fraction of averted deaths of 0.31 with rV = 1% to − 1.12 with the ten times slower rollout rV = 0.1%, marking an absolute drop of 1.43. This pattern is observed consistently also for other countries, and the effect is more intense for higher values of α. Similarly to previous results, we find that Egypt and Peru are generally more disadvantaged while the difference between countries reduce when the homogeneous or the strategy that reduces severity are employed.
To further explore how to mitigate the impact of behavioural relaxation once that the vaccination campaign started, we systematically explored the interplay between vaccine efficacy and vaccination rollout speed in Figure 3. The red dashed lines highlight the combinations of vaccine efficacy vs. rollout speed to achieve a 30% drop in observed deaths (taken as reference value), in the absence of a behavioural response (i.e., α = 0). In most countries, it can be achieved with a vaccination rollout speed ranging from 0.15% to 0.25% when the available vaccine has an efficacy ranging widely from 90% to 60%. On the contrary, when even a mild behavioural response is considered (α = 1), the rollout speed has to increase greatly when vaccine efficacy diminishes: for the case of Italy, a 30% drop in deaths would be achieved with a rollout speed of 0.47% when the vaccine efficacy is 90%, but in the case of V E = 60% the rollout speed rV has to grow up to 1.1% to achieve the same result.
Epidemic trajectories and impact of behavioural relaxation
As a way to ground the model on more realistic epidemic scenarios, we calibrate it to the real epidemic trajectories of the six countries considered. The calibration is performed via an Approximate Bayesian Computation technique (ABC) [55] on weekly deaths for the period 2020/09/01 − 2020/12/31. In doing so, we estimate some key epidemiological parameters matching the context of each country, thus defining realistic initial conditions for the model. To capture the effects of non-pharmaceutical interventions on the contacts between people we consider data from the Google Mobility Report [56] and the Oxford Coronavirus Government Response Tracker [57] until week 11 of 2021. After the calibration step, the model evolves from 2021/01/01 up to 2021/06/01 with the vaccination rollout for each of the countries. Note that in some of these countries, the vaccination campaigns officially started in the second half of December, though mostly symbolically.
In Figure 4A we represent key indicators resulting from the calibration. In particular, for the six countries we report the boxplot for the calibrated infection parameter β, and also the projected number of symptomatic infections per 100, 000 and the fraction of recovered individuals on the 2021/01/01, which is the start of the vaccination campaign in our simulation. We acknowledge a significantly high acquired immunity in the case of Peru, where the estimated attack rate as of 1st January, 2021, is around 38%. Nonetheless, this figure is line with other available estimates [54]. We observe how the calibration allows for an heterogeneous representation of the epidemiological conditions of the different countries. In Figure 4A we also represents the effects of restrictive measures on contacts. As a proxy, we consider the ratio between the leading eigenvalue of the contacts matrix considering the restrictions and the leading eigenvalue of the baseline contacts matrix without restrictions. The leading eigenvalue of the contacts matrix influences the reproductive number of the disease (see Materials and Methods), and more broadly it measures the intensity of contacts among people. By normalizing with respect to the baseline contacts matrix (with no restrictions) we can grasp the strictness of the measures in place. For example, a ratio of 0.3 would imply, in our simulations, a 70% reduction of the reproductive number with respect to the baseline without restrictions. The effect of restrictive measures on contacts varies over time up to week 11, 2021. In the case of Italy, for example, we observe the partial ease of the measures during January and February (i.e., ρ(Ct)/ρ(Cb) increases) followed by the tightening of measures in March aimed at curbing the third wave of infections. Afterwards, we keep the mixing levels as those observed for week 11 (more details in the Materials and Methods section). This is a conservative assumption as we can imagine that the seasonality and impact of the vaccination campaign might induce a relaxation of NPIs as the number of cases and deaths, hopefully, will go down. In this scenario, the countries with the strictest measures in place after week 11, are Peru and Italy, while Egypt and Ukraine show the most permissive measures. In the Supplementary Information, we repeat the analysis considering scenarios in which we partially relax some of the measures. The results are in line with the observations presented below.
In Figure 4B we represent the relative deaths difference for the calibrated model in the six countries. We study two vaccination rates. A fast deployment that manages to vaccinate rV = 1% of the population daily (aligned with the preparedness plan for Influenza pandemic [58] and vaccination rates achieved by countries like Israel or Chile), and a slower yet more realistic rollout rV = 0.25% (aligned with current rates in European countries). Worryingly, in the different countries, even in a scenario with several restrictions and a successful vaccination campaign in place, it is possible to witness a significant increase (factors reach a level of 2) in deaths exclusively due to a relaxation of COVID-safe behaviours. As a result of the calibration step, the countries under investigation face different epidemiological conditions (e.g., effective reproductive number, immunity from previous waves, initial number of infected, estimated effects of restrictions). These translate into phenomenological differences with respect to the results of the previous analyses. In particular, Ukraine and Egypt are significantly impacted by behavioural relaxation in the case of the slower rollout rV = 0.25%. From Figure 4A we see that in both countries the estimated effect of restrictions on contacts (especially in the last weeks) is smaller with respect to the others. This underlines how, besides the characteristics of the vaccination campaign, also the epidemiological conditions and the measures in place are important factors influencing the behaviour-vaccine interplay. In the case of the faster rollout rV = 1%, the ranking of the countries is similar to the one obtained previously (Peru and Egypt are generally more disadvantaged), when we imposed the same conditions for all countries. A possible explanation is that when a fast rollout is employed the influence of epidemiological conditions and restrictions become less important (because of the efficiency of the rollout), and the effects observed in the exploratory analysis emerge more easily.
Figure 4B shows that while the vaccination strategy targeting a reduction of severity performs slightly better when no and mild behavioural responses are in place, it appears to be much more robust when strong behavioural responses are considered. Also, while in the case of a faster rollout (rV = 1%) behavioural changes impact mildly the number of avoided deaths, with a 4 times slower rollout (rV = 0.25%) we find that the number of deaths consistently increases due to a relaxation of COVID-safe behaviours. In fact, a relaxation of protective behaviors not backed by a fast enough rollout, might induce an increase of deaths with respect to scenarios without vaccines. Finally, in case of homogeneous vaccination or prioritization aimed at reducing transmission overconfident behaviours further increase such death toll. These findings are consistent across the different countries considered.
In order to provide a measure about vaccination strategy robustness against behavioural relaxation, we sum the relative deaths difference (or relative symptomatic cases difference) across the spectrum of considered α values. In Figure 5 we show for every country how these two key metrics would be affected by different vaccination strategies when a wide range of behavioural responses are considered.
The plots show that vaccination strategy 1 is the most effective in reducing the number of deaths, while vaccination strategy 3 is the most effective in reducing the number of symptomatic cases, thus confirming that a potential behavioural response would not alter the already known vaccination prioritization strategies. Nonetheless, because of the unfolding of the epidemic in Egypt and Peru and the high intergenerational mixing, the vaccination strategy 2 (homogeneous across age classes) appears to be slightly better to minimize symptomatic cases when COVID-safe behaviours are lifted too soon. Information about vaccination strategy robustness against behavioural changes might be used to tune resilient roll-out campaigns. Policy-makers might also consider to optimize vaccination strategies with respect to a combination of multiple metrics. In this direction, our analysis provide additional insights. Interestingly, in Italy, Canada, Serbia and Ukraine, strategy 1 slightly outperforms strategy 3 in terms of robustness against behavioural changes to decrease the number of deaths, but strategy 3 greatly outperform strategy 1 in terms of robustness towards avoided infections, thus suggesting once again that behavioural changes play an important role that might modify our understanding of prioritization strategies.
Discussion
For almost a year, in the midst of a global pandemic, policymakers struggled to implement sustainable restrictions to slow SARS-CoV-2 spreading. Every non-pharmaceutical intervention was aimed at flat-tening the incidence curve and buying time for the development, test, production and distribution of vaccines that might ultimately protect the population. With an impressive scientific endeavor, several vaccines have been developed and an early distribution campaign was rolled out by the last days of December 2020. Besides the potential novel threats emerging from new virus strains [59], the current vaccination campaign represents the beginning of a new normal and a gigantic step towards complete virus suppression. However, we demonstrated that if the growth in vaccination uptake would lead to overconfident conducts inducing relaxation of COVID-safe behaviours, additional avoidable deaths will occur. We extended the literature proposing a mechanistic compartmental model able to simulate the unfolding of COVID-19, the vaccination dynamics and the compliance/non-compliance transition modulated by different behavioural mechanisms. Performing in-silico simulations allowed us to explore, from a theoretical standpoint, the interplay among different vaccination strategies, vaccine efficacy, rollout speeds, and behavioural responses. We found that behaviour impacts non-linearly vaccine effectiveness. Indeed, early NPIs relaxation affects more significantly lower vaccine efficacy, slower rollouts, and allocation strategies that target reduction in transmission rather in severity. We included in our analysis six different countries (Egypt, Peru, Serbia, Ukraine, Canada, and Italy), representatives of various points of the spectrum of world economies. This allowed us to observe the effects of behaviour-vaccine relationships for various population pyramids and mixing patterns. We observed that the developing economies, characterized by a younger population, but higher contacts activity and inter-generational mixing, are generally more affected by behaviour change with respect to developed economies. Then, as a way to ground the model on more realistic epidemiological conditions, we calibrated it using real epidemic and mobility data for the six countries considered and simulated the unfolding of the first months of the vaccination campaign. In such realistic scenario, we observed that with restrictive measures in place and a successful vaccination campaign, it is possible to witness to non-negligible increases in deaths attributable to an early relaxation of COVID-safe behaviours. The calibration step allowed us to highlight that also the epidemiological conditions related to the country-specific unfolding of the disease are an important factor influencing the interplay behaviour-vaccine.
We acknowledge some limitations in the present study. First, the vaccine is modeled as protective towards infection, without modeling the protection against severe COVID-19 complications and potential reduction of viral shedding in infected vaccinated individuals. At the time of writing, such information are still matters of study though initial reports support high level of protection from severe outcomes and a significant reduction in transmission [16, 17, 60]. For the sake of simplicity, we considered the vaccines fully working immediately after the first dose and we neglected that vaccination campaigns are using a portfolio of vaccines rather than a single one. We have also studied three simple vaccination strategies that neglect the complexities of an unprecedented mass vaccination. As result, both the vaccination priorities, vaccine effects and vaccination rates are an approximation of reality. In the Supplementary Information, we studied a data-driven vaccination rollout for Italy, where vaccines are distributed to the various age brackets following the real daily administration data. The results we obtained are qualitatively similar to those presented in the main text. While the model calibration suggests that our approach can nicely capture national trends, the model is not meant to provide accurate forecasts of the local unfolding of the disease, but rather to test what-if scenarios in a comparative fashion. We have considered a simple age-structure compartmental model that does not capture spatio-temporal heterogeneity both in terms of spreading and of NPIs implementation which have instead been observed in the countries under investigation. Our model does not include new, more transmissible, variants of concern and assumes the same IFR across the six countries. Finally, we propose and model two potential mechanisms leading to behavioural changes, but data are not available to perform a quantitative validation of the behavioural components of the model.
Implementation of individual protective behaviours and adherence to NPIs have been vital in order to reduce the transmission of SARS-CoV-2 leading to substantial population-level effects [5, 61–72]. Behavioural science can provide valuable insights for managing policies, incentives, communication strategies and can help coordinate efforts to control threats and evaluate such interventions [73]. As during the first waves of COVID-19, when NPIs were the only available mitigation measures [74], the results of our paper call for adequate strategies to keep high the attention and compliance towards individual COVID-safe behaviours, such as mask-wearing, social distancing, and avoidance of large gatherings now that vaccines are finally available. Communication strategies and policies should keep targeting such non-pharmaceutical intervention to avoid frustrating the immense effort of the vaccination campaigns.
Materials and Methods
Epidemic Model
Model Definition
We propose a compartmental model that incorporates both the vaccination process and the behaviour dynamics possibly linked to it. See Figure 6 for a schematic representation. The virus transmission is modeled using the following approach. Susceptible individuals (S compartment) in contact with infectious compartments become latent (L compartment). After the latent period (ϵ −1), L individuals enter the pre-symptomatic stage of the infection (P). Then, they can transition either in the asymptomatic (A) or the symptomatic stage (I) at rate ω (the length of time spent in the L and P compartments is the incubation period). The probability of being asymptomatic is f. After the infectious period (µ−1), both I and A individuals enter the Recovered compartment (R). Alternatively, I individuals can also die and transit to the D compartment. Note how we include a delay of Δ days from the time individuals enter the compartment D and die (Do compartment, the superscript stands for “observed”). The infectious compartments are P, A, I. The transmission rate is β and the force of infection is dependent on the age-stratified contact matrix C ∈ ℝ K×K, whose element Cij represents the average number of contacts that an individual in age group i make with individuals in j per day. The matrix C has four location-specific contributions: contacts at home, workplace, school, and other locations. We adopt country-specific contacts matrices provided in Ref. [75]. We assume that P and I have lower infectiousness with respect to symptomatic I (βχ, with χ < 1). Similar approaches have been used in several modeling studies in the context of the current pandemic [67, 76].
On top of the disease dynamics, we model both the vaccination process and the behavioural change that is possibly coupled to the vaccination. More in detail, after the start of the vaccination campaign at tV, at each time step, a fraction rV of the susceptible population receive a vaccine and transit to compartment V. As mentioned above, we consider three vaccination strategies: one in which the vaccine is given in decreasing order of age (vaccination strategy 1 aimed at reducing the severity), one in which it is given homogeneously to the population (vaccination strategy 2), and one in which it is first given to individuals in age brackets 20-49 and then to the rest of the population (vaccination strategy 3 aimed at reducing the transmission of the virus) [10]. For convenience, we derive the model’s equations in the homogeneous case, we refer the reader to the following sections for more details. We consider a “leaky” vaccine that reduces susceptibility with a certain efficacy V E. In other words, the infection rate for V individuals is β(1 − V E). In parallel, we imagine a behavioural dynamics triggered by the presence of the vaccine. Indeed, susceptible individuals (both vaccinated and not) may start adopting less safe behaviours because reassured by the presence of an effective vaccine. This is encoded in the model with a transition from the compartment S (V) to a new compartment SNC (V NC) of individuals that protect less themselves and as a result get infected at a higher rate (rβ, with r > 1). Note that NC stands for non-compliant. We propose two mechanisms to model this transition. In the first mechanism, the transition towards non-compliance happens at a rate α and it is catalyzed by the cumulative fraction of vaccinated individuals vt. The opposite transition from SNC (V NC) to S (V) happens at rate γ and is catalyzed by the number of observed deaths per 100, 000 individuals in the previous time step . Indeed, an increase in deaths is frequently used - especially by media - as an indicator of the severity of the current epidemiological situation. Existing literature suggests that risk perception (in the form of number of infected individuals or deaths) and communication of such risk significantly affect adherence to personal mitigation strategies such as social distancing and wearing face masks [5, 77]. Expressing the number of deaths in proportion to the population allows us to compare countries of different size. In the second mechanism, S (V) individuals transit to the non-compliant compartment SNC (V NC) at a constant rate α. We also account for the possibility of non-compliant going back to safer behaviours again at a constant rate γ. To simplify the narrative, we present the results considering this second mechanism only in the Supplementary Information. The overall picture discussed above does not change significantly. Note that in order to avoid issues with transition probabilities large than one we model the rates as with X = [S, V] and where the specific expression of the exponent depend on the two mechanisms described above. Note how for small values of g(α), h(γ) the rate converges to the usual mass-action law. The model can be written down as the following system of differential equations for individuals in age group k: Where δ[t − tV] is a step function which is equal to 1 if t > tV and to 0 otherwise. The basic reproduction number is , where is the contacts matrix weighted by the relative population in different age groups. In the Supplementary Information we provide details on its derivation. We adopt the model integrating numerically the equations, thus it is deterministic. However, it is worth stressing that when we calibrate the model to the real epidemic trajectory in the six countries, we use a probabilistic framework through an Approximate Bayesian Computation technique. Said differently, the calibrated parameters are characterized by posterior probability distributions rather than exact values. For this reason, the results presented in Figure 4 (i.e., the median of model’s projections) are to be intended as an ensemble of multiple trajectories generated sampling from the posterior distribution of the fitted parameters. The model is implemented in the programming language python qwith the use of the libraries scipy [78], numpy [79], and numba [80]. The visualization of the results is realized with the library matplotlib [81].
Phase space exploration
In Figure 7 we explore the phase space of the behavioural parameters (α, γ) of the model. We set for the six countries R0 = 1.15, 1% of initially infected individuals, and 10% of immune individuals, rV = 1%, and V E = 90%. After setting the initial conditions, we let the model evolve, individually for each country, for one year exploring a grid of the parameters α and γ. This allows us to observe the phase space of parameters regulating the behavioural transitions. In particular, for each (α, γ) pair we compute the relative deaths difference due to vaccines and behaviour change. As stated previously, the relative deaths difference is the fraction of deaths averted with respect to a baseline without vaccine and behavioural response. For the countries of focus, we consider two different values of the parameter r (r = 1.3, 1.5) which defines the increase in infection risk for individuals relaxing preventive behaviours. The obtained relative deaths difference varies from a maximum of 0.77 to a minimum of 0.26. This indicates that, in our simulations in the best scenario about 80% of deaths are averted thanks to the vaccine rollout. However, this potential gain reduces to around 25% in the worst case, with a potential waste of 55% of the benefit brought by the vaccine in terms of reduced mortality. Since α and γ are the only varying parameter in these simulations, such a reduction is only attributable to the relaxation of NPIs from the individuals. More in detail, across the different settings considered, a common pattern emerges. For α fixed, as γ grows we observe progressively an increase in the relative deaths difference. Indeed, if the population reacts promptly, non-compliant individuals turn back to COVID-safe behaviours, and the fraction of averted deaths benefits from this behaviour. On contrary, for a fixed γ an increase of α induces a stronger behavioural response causing more deaths otherwise avoided thanks to the vaccination. Furthermore, we observe that, for a given pair of the behavioural parameters, the fraction of averted deaths is lower when r = 1.5 with respect to r = 1.3. Indeed, in this case non-compliant individuals expose themselves to a higher risk of infection.
In the Results section, we kept for simplicity γ = 0.5 and we let vary α. This choice is informed by the maximum number of deaths observed in the countries of focus. In Italy, for example, the maximum number of deaths reported on a single day is around 1, 000. Therefore, this value of γ is such that, in a similar situation, non-compliant individuals would likely return to COVID-safe behaviours. As described in the previous paragraph, in the dynamic-rate behavioural mechanism we model the transition rate as , where is the number of deaths per 100, 000 individuals observed at time t − 1. Hence at the peak of deaths the transition rate towards compliance is . The sensitivity to this choice is discussed in the Supplementary Information.
Vaccination Strategies
As mentioned above, we consider three vaccination strategies. In the first one, the vaccine is given in decreasing order of age. Since the IFR of COVID-19 strongly correlates with age, many countries world-wide are adopting similar strategies prioritizing the vaccination among the elderly. Previous modeling works in the context of COVID-19 showed that this is the preferable strategy when considering the number of averted deaths [10, 42, 43]. In practice, this means that in our simulations we start giving the vaccine to the 75+ age bracket and we proceed in decreasing order only when everyone in this group is vaccinated. In the second strategy the vaccine is given homogeneously to the population respecting the age distribution. This means that, of the X vaccines available at step t, the fraction Nk/N is given to age group k (where Nk is the number of individuals in age bracket k, and N is the total number of individuals). In the third strategy, we prioritize age groups 20 − 49. Indeed, individuals at high mortality risk from COVID-19 disease may be protected indirectly by vaccinating age brackets that sustain the transmission [82, 83]. In practice, in our simulations we start giving the vaccine homogeneously to the 20 − 49 age brackets, and when they are all vaccinated we give it homogeneously to all other groups. In the previous two strategies, if, on a given time step, the number of people remained in the age group which is currently being vaccinated is smaller than the total number of vaccines available, the exceeding part is given to the next age group in line. This implies that the number of vaccines given in different time steps is always constant. In all cases, the vaccines are given only to the susceptible population.
Model Calibration
In this section we describe the methods used to calibrate the model to the real epidemic trajectories in the countries considered. We use an Approximate Bayesian Computation technique to calibrate the model on weekly deaths during the period 2020/09/01 − 2020/12/31. We account for government-mandated restrictions and their effect on the contacts between individuals modifying the contacts matrices using data from the Google Mobility Report [56] and the Oxford Coronavirus Government Response Tracker [57].
Approximate Bayesian Computation
We calibrate the model for each country using the Approxi-mate Bayesian Computation (ABC) rejection method [55, 84]. At each step of the rejection algorithm, a set of parameters θ is sampled from a prior distribution and an instance of the model is generated using θ. Then, an output quantity E′ of the model is compared to the corresponding real quantity E using a distance measure s(E′, E). If this distance is greater (smaller) than a predefined tolerance ψ, then the sampled set of parameters is discarded (retained). After accepting N sets, the iteration stops and the distribution of accepted parameters is an approximation of the real posterior distribution P (θ, E). The free parameters of our model and the related prior distributions are:
The transmission rate β. We explore uniformly values of β such that the related R0 is between 0.8 and 2.2. The basic reproduction number of SARS-CoV-2 is higher [85], but we consider lower values since our calibration starts on 2020/09/01 when restrictions were in place to mitigate the spreading.
The delay in deaths Δ ∼ [14, 25]. Indeed, for COVID-19 the average time between symptoms onset and death is about 2 weeks [86] and we also account for possible additional delays in death reporting.
The initial number of infected individuals. We explore uniformly values between 0.5 and 15 times the number of reported cases in the week before the start of the simulation. We then assign these individuals to the infected compartments (L, P, A, I) proportionally to the time spent there by individuals (ϵ−1 for L, ω−1 for P, and µ−1 for I, A), and we split between I and A individuals considering the fraction of asymptomatic f.
The model is calibrated on the period 2020/09/01 − 2020/12/31 using the weighted mean absolute percentage error on weekly deaths as an error metric with a tolerance ψ = 0.3 (ψ = 0.4 for Egypt for convergence issues) and 5, 000 accepted parameters set.
Model Initialization
The number of individuals in different age groups is initialized considering the 2019 United Nation World Population Prospects [87]. We consider 16 age brackets of five-years, except for the last one that includes 75+ individuals. As specified above, the initial number of infected individuals is calibrated considering the total number of confirmed cases in the week before the start of the simulation according to the data issued from the European Centre for Disease Prevention and Control [88]. To initialize the number of non-susceptible individuals (placed in the R compartment) we compute the average of several publicly available projections of total COVID-19 infections up to 2020/09/01 (i.e., start of the simulation) from different modeling approaches [54]. Both the initial number of infected and of non-susceptibles is assigned homogeneously across age groups. Other parameters used are E−1 = 3.7 days, µ−1 = 2.5 days, ω−1 = 1.5 days χ = 0.55, f = 0.35 in line with current estimates of COVID-19 infection dynamics parameters [89–93]. We use the age-stratified Infection Fatality Rate (IFR) from Ref. [53].
Modeling the Effects of NPIs on Contact Matrices
We incorporate in our model the implementation of top-down NPIs by changing the contacts patterns defined by the contacts matrix C. As mentioned, we use the country-specific contacts matrices provided in Ref. [75]. These are made by four contributions: contacts that happen at school, workplace, home, and other locations. In a baseline scenario, the overall contacts matrix C is simply the sum of these four contributions. Following Ref. [94], we implement the reductions in contacts, due to the restrictions, multiplying the single contribution by a reduction factor ηi(t). Thus, in general the overall contacts matrix at time t become: For simplicity, we assume no changes to contacts at home, though lockdowns tend to increase them [61]. For the contacts locations work and other locations we use data from the Google Mobility Report [56]. In detail, we consider the fields workplaces percent change from baseline to model reduction in the contacts location work and the average of the fields retail and recreation percent change from baseline and transit stations percent change from baseline for other locations. A general entry of the report pl(t) represents the percentage change on day t of total visitors to a specific location l with respect to a pre-pandemic baseline. From pl(t) we derive a contacts reduction coefficient as follows: ηl(t) = (1 + pl(t)/100)2. Indeed, the number of potential contacts in a specific location is proportional to the square of the the number of visitors. We model contacts reduction in school considering instead the Oxford Coronavirus Government Response Tracker [57]. More in detail, we use the ordinal index C1 School closing. The index ranges from a minimum of 0 (no measures) to a maximum of 3 (require closing all levels). We turn this quantity into the contacts reduction coefficient in school as follows: ηs(t) = (3 − C1 School closing)/3. We use these datasets to inform contacts reductions up to week 11, 2021. After, we assume that contacts remain at the same level of week 11. Note how we let vary both the transmission rate β and the contact matrix. The former describes the risk of infection given contacts with infectious individuals. This is function of the disease (which is assumed to be the same, we don’t consider multiple or emergent new strains possibly more transmissible) and of the protective behaviours such as social distancing and use of face masks. The latter describes variations to the number and types of contacts induced by top-down NPIs as for example remote working, schools closure and lockdowns. By splitting the contributions to the force of infection of transmission rate and contact matrix we are able to take into consideration different behavioural attitudes which, given the same number of contacts, might lead to higher or lower risks of infection. This allows us to consider explicitly both top-down and bottom-up NPIs.
Calibration Result
In Figure 8, we report the results of the calibration. It is important to stress how our goal is not to develop a predictive model aimed at forecasting the pandemic trajectory. The fit is used to ground the model and to define the epidemic conditions at the start of the vaccination campaign in the six countries. In fact, our aim is to understand the possible interplay between behaviours and vaccine rollout which is also function of the epidemic progression. In the Figure we report the official and simulated weekly number of deaths (median, 50% and 95% confidence interval). Despite its simplicity and approximations, the model is able to reproduce the evolution of the pandemic in the six countries capturing well its progression after the summer. In the Supplementary Information we also report the posterior distributions for the parameters.
Data Availability
All the data and code used in this paper are publicly available at https://github.com/ngozzi/vaccine-behaviour
Author Contributions
NG, PB, NP conceived the research and designed the analyses. NG conducted the numerical simulations. All authors wrote the manuscript.
Competing Interests
The authors declare that they have no competing financial interests.
Data Availability
All the data and code used in this paper are publicly available at https://github.com/ngozzi/vaccine-behaviour.
Acknowledgements
All authors thank the High Performance Computing facilities at Greenwich University. NG acknowledges support from the Doctoral Training Alliance. PB acknowledges support from Intesa Sanpaolo Innovation Center. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Footnotes
↵* paolo.bajardi{at}isi.it
References
- [1].↵
- [2].
- [3].
- [4].
- [5].↵
- [6].↵
- [7].↵
- [8].↵
- [9].
- [10].↵
- [11].↵
- [12].↵
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].
- [23].↵
- [24].↵
- [25].↵
- [26].
- [27].
- [28].
- [29].
- [30].
- [31].
- [32].
- [33].
- [34].
- [35].
- [36].
- [37].
- [38].
- [39].
- [40].↵
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].
- [47].
- [48].↵
- [49].↵
- [50].↵
- [51].↵
- [52].↵
- [53].↵
- [54].↵
- [55].↵
- [56].↵
- [57].↵
- [58].↵
- [59].↵
- [60].↵
- [61].↵
- [62].
- [63].
- [64].
- [65].
- [66].
- [67].↵
- [68].
- [69].
- [70].
- [71].
- [72].↵
- [73].↵
- [74].↵
- [75].↵
- [76].↵
- [77].↵
- [78].↵
- [79].↵
- [80].↵
- [81].↵
- [82].↵
- [83].↵
- [84].↵
- [85].↵
- [86].↵
- [87].↵
- [88].↵
- [89].↵
- [90].
- [91].
- [92].
- [93].↵
- [94].↵