A behavioural modelling approach to assess the impact of COVID-19 vaccine hesitancy.

We introduce a compartmental epidemic model to describe the spread of COVID-19 within a population, assuming that a vaccine is available, but vaccination is not mandatory. The model takes into account vaccine hesitancy and the refusal of vaccination by individuals, which take their decision on vaccination based on both the present and past information about the spread of the disease. Theoretical analysis and simulations show that voluntary vaccination can certainly reduce the impact of the disease but is unable to eliminate it. We also demonstrate how the information-related parameters affect the dynamics of the disease. In particular, vaccine hesitancy and refusal are better contained in case of widespread information coverage and short-term memory. Finally, the possible impact of seasonality on the spread of the disease is investigated.

Introduction

On 31 December 2019, Chinese public health authorities informed the WHO about a cluster of anomalous viral pneumonia cases in Wuhan city (WHO, 2019). The causal agent of the disease was identified shortly after as a new type of SARS, which was named SARS-CoV-2. Although many governments underestimated the risk of a pandemic breaking out (Magli et al., 2020), since 21 January 2020, the WHO began publishing daily situation reports on its website. The first report clearly states that the ‘WHO has issued interim guidance for countries, updated to take into account the current situation’ (WHO, 2019). Indeed, the first case outside of China was reported on 13 January 2020, and then the disease rapidly spread to other countries. Finally, it led to a devastating pandemic as we know it, causing the temporary collapse of many health systems. For example, France in the pre-COVID-19 era had about 5000 ICU beds; however, at the peak of its first wave, 7019 ICU beds were occupied by COVID-19 patients (French Public Health Agency, 2020).

During the first year of the pandemic, in the absence of a vaccine, the only possible pandemic-mitigation strategies at the local level were based on social distancing and partial and full lockdowns (Buonomo and Della Marca, 2020, Magli et al., 2020). Lockdowns were generally very effective in reducing the pressure of the pandemic on the health systems of countries but, in many countries, after a few months there was a new epidemic outbreak. Till now, most countries went through three epidemic outbreaks (also termed waves) (Center for Systems Science and Engineering at Johns Hopkins University, 2020).

Since the early stage of the pandemic, many researchers proposed models, based on traditional mathematical epidemiology, for understanding the evolution of COVID-19 and in turn controlling it (Kucharski et al., 2020, Gatto et al., 2020, Della Rossa et al., 2020, Giordano et al., 2020, Davies et al., 2020, Ngonghala et al., 2020, Dolbeault and Turinici, 2020). The early dynamics of disease transmission in Wuhan was studied by Kucharski et al. (2020) with a stochastic SEIR model using the data on the outbreak in Wuhan. Gatto et al. (2020) proposed a model to study the disease transmission across a network of 107 Italian provinces during the initial stage of the first COVID-19 wave. A network model applied to Italy was also proposed also by Della Rossa et al. (2020) to show that heterogeneity between regions plays a fundamental role in designing effective strategies to control the disease while avoiding national lockdowns. Giordano et al. (2020) introduced a model for assessing the effectiveness of testing and contact tracing combined with social distancing measures. Non-pharmaceutical interventions to fight COVID-19 in the UK and US were considered by Davies et al. (2020) and Ngonghala et al. (2020) respectively, while the effect of social distancing during lockdowns in France was studied by Dolbeault and Turinici (2020) using a variant of the SEIR model. Many other relevant studies have focused on assessing the effects of containment measures and predicting epidemic peaks and ICU accesses, see e.g., (Elie et al., 2020, Flaxman et al., 2020, Supino et al., 2020). As soon as COVID-19 vaccines became available, many compartmental models have began to appear in the literature with the specific aim of investigating vaccination’s effects on the spread of the disease as well as determining the optimal allocation of vaccine supply (Buckner et al., 2021, Choi and Shim, 2020, Mukandavire et al., 2020, Deng et al., 2021). A limitation of classical Mathematical Epidemiology (ME) is that it is built up on Statistical Mechanics; the agents are modelled as if they were molecules, and the contagion is abstracted as a chemical reaction between ‘molecules’ of the healthy species with those of the infectious species. Thus, laws like the mass action law are used in such models. The missing ingredient in ME is the behaviour of agents, specifically how people modify their contacts at risk and how their vaccine-related decisions are taken. The absence of this ingredient makes classical ME models increasingly less adaptable as a tool for Public Health. Indeed, a major challenge for global Public Health is the spread of vaccine hesitancy and refusal. This is due to the phenomenon of ‘Pseudo-Rational’ Objection to VAccination (PROVA) (Buonomo et al., 2013): people overestimate real and imaginary side effects of vaccines and underestimate real risks of the target infectious diseases (Buonomo et al., 2013, Manfredi and d’Onofrio, 2013, Wang et al., 2016). PROVA is inducing remarkable changes in the civil society’s attitude towards the prevention of infectious diseases. This increasingly severe lack of trust in vaccination is one of the many negative consequences of two distinct and synergysing phenomena of a more general nature: the post-trust society (Löfstedt, 2005) and the post-truth era (McIntyre, 2018).

The first work that explicitly modelled social distancing in ME was (Capasso and Serio, 1978), which incorporated a phenomenological behavioural response to the Kermack-McKendrick epidemic model. The spread of PROVA in the last two decades led to the birth of a new branch of ME: the Behavioural Epidemiology of Infectious Diseases (BEID) (Manfredi and d’Onofrio, 2013, Wang et al., 2016). The main aim of the BEID is to incorporate the impact of human behaviour in models on the spread and control of infectious diseases (Manfredi and d’Onofrio, 2013, Wang et al., 2016). The key role of both the present and past information in vaccination decisions and uptake as well as social distancing was first stressed, respectively, in (d’Onofrio et al., 2007, Manfredi and d’Onofrio, 2013, d’Onofrio and Manfredi, 2009) by means of phenomenological models. In a recent paper (Buonomo and Della Marca, 2020), a model describing the transmission of COVID-19 has been introduced. The model considers social distancing and quarantine as mitigation strategies implemented by the Public Health system. The model is information-dependent in that the contact and quarantine rates are assumed to be dependent on the available information and rumours about the disease status in the community. In (Buonomo and Della Marca, 2020), the model is applied to the case of the COVID-19 epidemic in Italy. The paper argues that citizen compliance along with mitigation measures played a decisive role in curbing the epidemic curve by preventing a duplication of deaths as well as about more infections.

The COVID-19 pandemic caused a worldwide research effort that resulted in the rapid development of new vaccines (Logunov et al., 2021, Knoll and Wonodi, 2021), some of which belong to the new class of mRNA vaccines (Baden et al., 2021, Polack et al., 2020). In light of the significant changes in the lives of milliards of people and the huge negative impact on the world economy, one could have expected that only a tiny proportion of people would really be hesitant towards vaccination. Unfortunately, this was not the case. As early as June 2020, Neumann-Böhme and co-workers (Neumann-Böhme et al., 2020) investigated attitudes about anti-COVID-19 vaccination from a representative sample of citizens from seven European countries. Surprisingly, although the first epidemic wave in Europe had just ended, hesitancy and opposition towards the vaccines were found among a large proportion in all classes, age groups, and sexes. In particular, of the French respondents were hesitant () or were strongly against () COVID-19 vaccines.

Before mid-December 2020, phase 3 of the experimentation of many vaccines ended, indicating that they have outstanding effectiveness in preventing COVID-19 (Baden et al., 2021, Polack et al., 2020, Logunov et al., 2021). Typically, drug regulatory agencies defined priority groups for vaccination (elderly people with serious co-morbidities, healthcare workers in senior residences, etc.). From a rational viewpoint, there were all the premises to believe that vaccine hesitancy had been strongly reduced and that mandatory vaccination campaigns could have been conducted, but this was not the case. Concerning the mandatory nature of the vaccination campaign, in many countries, the vaccines are not mandatory (Macron, 2020, La Stampa, 2020, The Guardian, 2020). As for vaccine hesitancy, an investigation conducted in October 2020 (IPSOS, 2020) suggests that of French citizens are vaccine hesitant. Other countries exhibit percentages of opposition and hesitancy that exceed in Spain and USA, in Italy, in South Africa, and in Japan and Germany. Globally, the hesitancy and objection are as large as .

Given these large percentages of hesitance and opposition to COVID-19 vaccines, we think that employing a behavioural epidemiology approach to model the implementation of a vaccination campaign for COVID-19 is appropriate. To this end, we adopt a strategy similar to the one used in (d’Onofrio et al., 2007). In other words, we assume that the vaccination rate is a phenomenological function of the present and past information that the citizens have on the spread of the epidemic. Note that, in the context of SIR and SEIR infectious diseases, more mechanistic models based on evolutionary game theories (Bauch, 2005, d’Onofrio et al., 2011, d’Onofrio et al., 2012, Wang et al., 2016) exist, but they are reduced to the approach of (d’Onofrio et al., 2007, Buonomo et al., 2013) in case of volatile opinion switching (Della Marca and d’Onofrio, 2021, d’Onofrio et al., 2011, Wang et al., 2016).

In this paper, we consider a COVID-19 affected population controlled by vaccination, where the final choice to get vaccinated or not is partially determined on a voluntary basis and dependent on publicly available information on the spreading of the disease in the community in the both the present and recent past. Our model is inspired by the compartmental epidemic model introduced in (Buonomo and Della Marca, 2020), whereby COVID-19 transmission during the 2020 lockdown in Italy was studied. An analogous situation was considered by Gumel and co-workers for the SARS epidemic when they examined a SARS model in (Gumel et al., 2004) and then considered vaccination intervention in (Gumel et al., 2006).

We performed a qualitative analysis based on stability and bifurcation theories. The analysis shows that, when the control reproduction number, , is less than 1, only a disease-free equilibrium (DFE) that is globally stable exists; when , the DFE is unstable, and an endemic equilibrium arises. The model is then parameterised based on the COVID-19 epidemic in Italy and preliminary reports on COVID-19 vaccines. In numerical simulations, we consider two possible starting times for a one-year vaccination campaign. We assess the role of critical model parameters by evaluating how they affect suitable epidemiological indicators. Finally, the effects of seasonality are investigated by making the assumption that disease transmission and severity as well as vaccination rates are lower during the warmer months.

The paper is organised as follows. In Section 2, the model is introduced, and in Section 3, the qualitative analysis is performed. Model parametrization and numerical solutions are given in Section 4 and Section 5, respectively. The case of seasonally varying parameter values is addressed in Section 6. Finally Section 7 presents concluding remarks. The paper is also complemented by a file of supplementary materials.

The model

State variables and the information index

In this section, we introduce a compartmental epidemic model to describe the spread of COVID-19 within a population, assuming that a vaccine is available, but vaccination is not mandatory.The total population at time t (say N) is divided into the following six disjoint compartments:The size of each compartment at time t represents a state variable of the mathematical model, and .

susceptible, S: healthy individuals who can contract the disease;

exposed (or latent), E: individuals who are infected by SARS-CoV-2 but not yet capable of transmitting the virus to others;

asymptomatic infectious, : this compartment includes two groups, namely the post-latent individuals, i.e., individuals who are in the incubation phase following latency, where they are infectious and asymptomatic, and the truly asymptomatic individuals, i.e., those having no symptoms throughout the course of the disease;

symptomatic infectious, : infectious individuals who show mild or severe symptoms;

vaccinated, V: individuals who received at least one dose of COVID-19 vaccine; and

recovered, R: individuals who recovered after the infectious period.

We assume that the individuals take their decision on vaccination based not only on the present information but also the past information they have on the spread of the disease, with the past information being weighted in an exponential way. Therefore, the information on the status of the disease in the community is described by means of the information index (d’Onofrio et al., 2007, Wang et al., 2016):

Such an index is an important tool in BEID (Manfredi and d’Onofrio, 2013). It is an extension of the idea of the prevalence-dependent contact rate, developed by V. Capasso in the seventies, which describes the behavioural response of individuals to prevalence (Capasso and Serio, 1978). Here, the parameter a takes the meaning of the inverse of the average time delay in the collection of the information on the disease (say ), and the parameter k is the information coverage, which summarises two opposite phenomena: the disease’s under-reporting and the level of media coverage of the disease status, which tends to amplify the social alarm. It may be assumed that ; see (Buonomo et al., 2008).

From (1), by applying the linear chain trick (MacDonald, 2008), we obtain the differential equation , determining the dynamics of M.

Assumptions regarding vaccine and immunity

Besides the role of human behaviour, the other two major hypotheses of the model are that vaccines are not perfect and both disease-induced immunity and vaccine-induced immunity do not wane. The first hypothesis is related to the scientific results from the phase 3 clinical trials as well as general knowledge concerning vaccines. The second hypothesis is stronger, and some could read it as an extreme optimistic case. This assumption is based on some very recent experimental results (Iyer et al., 2020, Wajnberg et al., 2020) and on an experimental review paper (Karlsson et al., 2020) on one of the most complex and intriguing topic concerning SARS-CoV-2: the immunological response associated with it. In particular, Iyer and co-workers (Iyer et al., 2020) showed that the igG response has practically no conversion for a long period after the onset of symptoms; in fact, only three out of 90 individuals had igG seroconversion. This very limited fraction of seroconversion can be taken into account (by means of a coefficient ; see Section 2.4) as some vaccinated individuals get infected because they had seroconversion of their vaccine-induced immune response. Moreover, in their review paper on T cell immunity to COVID-19 (Karlsson et al., 2020), Karlsson and co-workers stressed that the ‘generation of memory T cells can provide lifelong protection against pathogens. Previous studies have demonstrated that SARS-CoV- and MERS-CoV-specific T cells can be detected many years after infection. Likewise, SARS-CoV-2-specific CD4+ and CD8+ T cells are distinguished in a vast majority of convalescent donors […] Preliminary results from the two major mRNA vaccine trials in humans have demonstrated potent Th1 responses.’As far as the effectiveness of vaccine immunization is concerned, for some authorised mRNA COVID-19 vaccines, it has been estimated that the first dose induces – 14 days after inoculation – a reduction of 80% of the risk of being infected by SARS-CoV-2 (CDC, 2021) and a decreased risk of 90% in getting a very serious disease requiring hospitalization (Vasileiou et al., 2021). Thus, we approximately consider a single dose instead of two in our model.

Assumptions regarding transmission

Global research on how SARS-CoV-2 is transmitted continues to be conducted at the time of writing this paper. It is believed that infected people appear to be most infectious just before (around 1–2 days before) they develop symptoms (i.e., in the post-latency stage) and early in their illness (WHO, 2021). Recent investigations confirmed that pre-symptomatic transmission was more frequent than symptomatic transmission (Bender et al., 2021). The possibility of contagion from a truly asymptomatic COVID-19-infected person (i.e., an infected individual who does not develop symptoms) is still a controversial matter. However, it has been shown that little to no transmission may occur from truly asymptomatic patients (Bender et al., 2021).

In our model, the routes of transmission from COVID-19 patients are included in the Force of Infection (FoI) function, i.e., the per capita rate at which susceptibles contract the infection. As in (Gumel et al., 2006), the mass action incidence is considered:where .

The rationale behind this choice is that during the observed COVID-19 outbreaks, the total population remained effectively constant. For instance, in Italy (one of the countries more affected by the epidemic (Worldometer, 2020), the drop in the total population () due to COVID-19-induced deaths ( as of 19 April 2021 (Italian Ministry of Health, 2020b) is around . In this case, we expect mass action and standard incidence to give similar results.

In (2), the parameters and are modification factors that represent the level of reduced infectiousness of compartments and when compared with the subgroup of comprising post-latent individuals. Therefore, the baseline transmission rate denotes the transmission rate of post-latent individuals (see also Section 4.2, where and are estimated). For the reasons discussed above, we assume that the factor concerning the post-latent individuals is 1.

Description of the balance equations

All the state variables decrease due to natural death, at a rate . The susceptible population S increases due to the net inflow , incorporating both new births and immigration as well as decreases due to transmission and vaccination. For the time span covered in our simulations, demography could be neglected. However, including a net inflow of susceptible individuals in the model allows one to consider not only new births but also immigration, which plays an important role during the COVID-19 epidemic and can be well estimated in some cases (Buonomo and Della Marca, 2020). Therefore, since the demographic parameters can be easily obtained from data, we prefer to use an SEIR-like model with demography as successfully done for SARS models (Gumel et al., 2004, Gumel et al., 2006).

The exposed (or latent) individuals E arise as a result of new infections of susceptible and vaccinated individuals and decrease due to development at the infectious stage (at the rate ). We assume that after the latency period, the individuals enter the asymptomatic compartment , which includes post-latent and truly asymptomatic individuals, as described in Section 2.1. The number of asymptomatic individuals diminish because they enter the compartment of symptomatic individuals (at a rate of ) or they recover (at a rate of ). Mildly or severely symptomatic individuals enter from the post-latency stage and exit due to either recovery (at rate ) or disease-induced death (at rate ). Vaccinated individuals V come from the susceptible class after vaccination and their number decreases due to infections (at a reduced rate , where ).

The vaccination rate is assumed to be the sum of two terms. The first one, , is a positive constant that represents the fraction of susceptibles per time unit that decide to get vaccinated independently from the information. The motivation for this term stems from the need to model some aspects of vaccine acceptance, for instance: i) some people are strongly in favour of vaccines and act coherently by getting vaccinated; ii) some categories such as patients and healthcare workers in senior care facilities will be strongly recommended to get the vaccine (in some countries, their vaccination will be even mandatory) (La Stampa, 2020). The second term, , is a continuous increasing function of the information index M, such that and . This function represents the fraction of susceptibles per time unit that decide to get vaccinated because of the social alarm caused by the disease.

Finally, recovered individuals come from the infectious compartments and and, as discussed in Section 2.2, acquire long-lasting immunity against the disease.

The equations

According to the aforementioned description, the time evolution of the state variables is determined by the following system of balance equations: with initial conditions

Since the Eqs. (3) do not depend on R, the dynamics of the recovered compartment can possibly be studied separately, by means of the equation

The flow chart in Fig. 1 illustrates all the processes included in the model; a description of each parameter together with their baseline values is given in Table 1 (see Section 4).

Fig. 1

Flow chart for the COVID-19 model (3)–(5). The population is divided into six disjoint compartments of individuals: susceptible , exposed , asymptomatic , symptomatic , vaccinated and recovered . Blue colour indicates the information-dependent process in the model, with ruled by (3f).

Table 1

Temporal horizon, initial conditions and baseline values of parameters of model (3)–(14). The last column lists the corresponding sources for each parameter value.

Parameter	Description	Baseline value	Source
t_f	Time horizon	365-395 days	Assumed
N_0	Initial total population	6.036·10^7	ISTAT (2020)
E(0)	Initial number of exposed individuals	I_a(0)(η+ν_a+μ)/ρ	Assumed
I_a(0)	Initial number of asymptomatic infectious individuals	7,322	Estimated from ISS, 2020, Italian Ministry of Health, 2020b
I_s(0)	Initial number of symptomatic infectious individuals	7,545	Estimated from ISS, 2020, Italian Ministry of Health, 2020b
V(0)	Initial number of vaccinated individuals	0	Italian Ministry of Health (2020b)
R(0)	Initial number of recovered individuals	203,968	Italian Ministry of Health (2020b)
M(0)	Initial value of the information index	kI_s(0)	Assumed
R_0	Basic reproduction number	1.428	See (8)
R_V	Control reproduction number	0.302	See (9)
Λ	Net inflow of susceptibles	1,762 days−1	Estimated from ISTAT, 2020, Italian Ministry of Foreign Affairs and International Cooperation, 2020
μ	Natural death rate	1.07·10^-2 years−1	ISTAT (2020)
β	Baseline transmission rate	2.699·10^-8 days−1	Estimated from Italian Ministry of Health (2020b)
q	Fraction of post–latent individuals that develop symptoms	0.15	Estimated from Italian Ministry of Health (2020c)
ε_a	Modification factor concerning transmission from I_a	q+(1-q)0.033	Estimated from Gatto et al. (2020)
ε_s	Modification factor concerning transmission from I_s	0.034	Estimated from Gatto et al. (2020)
φ_0	Information-independent constant vaccination rate	0.002 days−1	Assumed
σ	Factor of vaccine ineffectiveness	0.2	CDC (2021)
ρ	Latency rate	1/5.25 days−1	Estimated from ECDC, 2020, WHO, 2020
η	Rate of onset of symptoms	0.12 days−1	Estimated from ECDC, 2020, Italian Ministry of Health, 2020c, WHO, 2020
ν_a	Recovery rate of asymptomatic infectious individuals	0.165 days−1	Estimated from ISS, 2020, Italian Ministry of Health, 2020b
ν_s	Recovery rate of symptomatic infectious individuals	0.055 days−1	Estimated from ISS, 2020, Italian Ministry of Health, 2020b
δ	Disease–induced death rate	6.248·10^-4 days−1	Estimated from ISS, 2020, Italian Ministry of Health, 2020b
D	Reactivity factor of information–dependent vaccination	500μ/Λ	Estimated from Buonomo, 2020, d’Onofrio et al., 2007
φ_max	Upper bound of overall vaccination rate	0.02 days−1	Estimated from ISS, 2020, Shim, 2019, Lee et al., 2012, Fister et al., 2016
a	Inverse of the average information delay T_a	1/3 days−1	Buonomo and Della Marca (2020)
k	Information coverage	0.8	Buonomo and Della Marca (2020)

Qualitative analysis

The following proposition ensures that the solutions of model (3) are epidemiologically and mathematically well-posed.

The regiondefined bywith initial conditions(4)is positively invariant for model (3).

See supplementary materials, Section S.1.  □

Thus, we limit our analyses to the region .

The model given by Eqs. (3) has a unique disease-free equilibrium (DFE), obtained by setting the r.h.s. of Eqs. (3) to zero, given by

To establish the local and global stability of the DFE, we compute the basic and the control reproduction numbers. The basic reproduction number, , is a frequently used indicator for measuring the potential spread of an infectious disease in a community. It is defined as the average number of secondary cases caused by one primary infection over the course of the infectious period in a fully susceptible population. If the system incorporates vaccination strategies, the corresponding quantity is then named the control reproduction number and is usually denoted by .

The reproduction number can be calculated as the spectral radius of the next generation matrix FV−1, where F and V are defined as Jacobian matrices of the new infection appearance and the other rates of transfer, respectively, calculated for infected compartments at the disease-free equilibrium (Van den Driessche and Watmough, 2002. In this specific case, if in (3), that is when a vaccination program is not in place, we obtain the expression of ; otherwise, the corresponding can be computed. We remark that by introducingas the fraction of the population vaccinated at the disease-free equilibrium (7), one can express

The basic reproduction number of model (3) is given byand the control reproduction number is given by

See supplementary materials, Section S.2.  □

Note that with equality only if (i.e., ) or . That is, despite being imperfect, the vaccine (characterised by and ) will always reduce the reproduction number of the disease.

From Proposition 2, it follows that (Van den Driessche and Watmough, 2002): As far as the global stability of the DFE is concerned, we can prove the following proposition: To stress the role of the biologically-relevant parameters and , Proposition 3, Proposition 4 can be reformulated as follows:

The DFE is locally asymptotically stable if; in contrast, if, then it is unstable.

The DFE is globally asymptotically stable if.

See supplementary materials, Section S.3, for two alternative proofs. □

Let us consider the quantityIfand, then the disease will be eliminated from the community since the DFE is globally asymptotically stable. If, then no amount of vaccination can prevent a disease outbreak in the community because the DFE is unstable. Finally, if, then the DFE is unstable, that is a disease outbreak occurs and disease elimination is not possible.

We remark that a result similar to Proposition 5 was obtained by Gumel and co-workers in (Gumel et al., 2006), where the epidemic control of SARS by means of a constant vaccination rate was studied.

Of course, the inequality can be reformulated in terms of the fraction of the population vaccinated at the disease-free equilibrium, p, given in (10). Therefore, we can also state the following:

The disease will be eliminated from the community if , with given by

As far as the existence of endemic equilibria is concerned, we can prove the following:

If, then system (3) has no endemic equilibria.

If, then system (3) has a unique endemic equilibrium.

See supplementary materials, Section S.4.  □

To derive a sufficient condition for the occurrence of a transcritical bifurcation at , we can use a bifurcation theory approach. We adopt the approach developed in (Dushoff et al., 1998, Van den Driessche and Watmough, 2002), which is based on the general centre manifold theory (Guckenheimer and Holmes, 1983). In short, it establishes that the normal form representing the dynamics of the system on the centre manifold is given bywhereand

Note that in (12), (13), has been chosen as the bifurcation parameter, is the critical value of is the state variables vector, is the right-hand side of system (3), and and denote, respectively, the left and right eigenvectors corresponding to the null eigenvalue of the Jacobian matrix evaluated at criticality (i.e., at DFE and ).

Observe that is equivalent toso that the disease-free equilibrium is stable if , and unstable when .

The direction of the bifurcation occurring at can be determined from the sign of coefficients (12), (13). More precisely, if [resp. ] and , then at , there is a backward [resp. forward] bifurcation.

For our model, we have the following:

System (3) exhibits a forward bifurcation at DFE and.

See supplementary materials, Section S.5.  □

Parametrisation

The values of demographic and epidemiological parameters are based on the COVID-19 epidemic reported in Italy since the end of February 2020 (Italian Ministry of Health, 2020b). The values of vaccine-related parameters are mainly inferred from preliminary reports about COVID-19 vaccines and from the initial trend of Italy’s immunisation campaign. A detailed derivation of such quantities is reported in the following.

Initial conditions

To provide appropriate initial conditions that mark the beginning of an epidemic wave, we make the following considerations. After the first epidemic wave (February–May 2020), Italy experienced the ‘living with the virus’ period, characterised by a relatively low level of prevalence and loosening of restrictions. However, this period ended towards the second half of August 2020 when the virus regained strength and its prevalence progressively grew, marking the arrival of the second wave.

Since data available from the beginning of the second wave are reasonably more accurate than those from the epidemic’s onset, we consider them as initial data. More specifically, we take the official national data for infectious () and recovered (R) people as of 16 August 2020 (see Italian Ministry of Health, 2020b), which is estimated as the first time after the end of the first wave when the effective reproduction number exceeded the threshold 1 (Italian Ministry of Health, 2020d). For that period, the Italian National Institute of Health estimates the fraction of asymptomatic individuals w.r.t. the total cases as about 49.25%, i.e., (ISS, 2020). As far as the initial values of exposed individuals E and the information index M are concerned, in the absence of exact data, we infer them from the corresponding expressions at the endemic state, as given in (S.7) (see supplementary materials). Hence, we get and . Finally, the initial value of susceptible individuals S is obtained from , where is the total initial population as given in (Buonomo and Della Marca, 2020).

Baseline scenario

In the absence of empirical data about attitudes towards vaccination, we follow the approach of (Buonomo and Della Marca, 2019, Buonomo and Della Marca, 2020, d’Onofrio et al., 2007) and assume that is a Michaelis–Menten function (Murray, 1989)with . Similar to what has been done in (Buonomo and Della Marca, 2019, Buonomo and Della Marca, 2020, d’Onofrio et al., 2007), we set , where . This reparametrisation means an asymptotic overall rate of days−1. The ensuing vaccination function is

As of April 2021, the rate of COVID-19 vaccination in Italy was less than 400,000 administrations per day in a population of millions of inhabitants (ISS, 2020), although acceleration plans have been laid out. Here, we take days−1, potentially implying an upper bound of 0.02 days−1 in vaccination rate under circumstances of high perceived risk. This value is in line with data concerning the 2009 H1N1 influenza pandemic, whose daily rate of vaccination has been largely investigated, and it was below 2% of the total population (see Lee et al., 2012 and references therein). Furthermore, threshold values of 1–2% per day were also considered in epidemic models of dengue (Shim, 2019) and cholera diseases (Fister et al., 2016).

To obtain a baseline value for D, we set as in (Buonomo, 2020, d’Onofrio et al., 2007), where M varied in . Here, M varies in (see (6)); hence, we expect that could be a good starting point.

Concerning the factor of vaccine ineffectiveness, , and the information-independent constant vaccination rate, , in Section 5, numerical solutions obtained by varying both and are given. As pointed out in Section 2.2, an estimate of the vaccine efficacy is 80%, which implies . Regarding the baseline vaccination rate, we take , that is, days−1.

We estimate the rate at which symptoms onset as , where represents the fraction of infected people who develop symptoms after the incubation period, and days−1 is the post-latency rate, as provided in (Buonomo and Della Marca, 2020). The fraction q is also used to infer , the modification factor concerning transmission from , so we set , where 1 [resp. 0.033] is the modification factor concerning transmission from post-latent [resp. truly asymptomatic] individuals, as considered in the models (Buonomo and Della Marca, 2020, Gatto et al., 2020).

Following the approach adopted by Gumel et al. (2004), based on the formula given by Day (2002), we estimate the disease-induced death rate aswhere is the fatality rate, and is the expected time from the onset of symptoms until death. We compute from the official national data from 16 August to 13 October 2020 (Italian Ministry of Health, 2020b) (the same period considered for the estimation of the transmission rate , as explained below), yielding %. As far as is concerned, from (ISS, 2020), we get days, providing days−1.

Similarly, the recovery rates with are estimated aswhere [resp. ] is the expected time until recovery for asymptomatic [resp. symptomatic] individuals. We assume and days on the basis of the considerations made in (Buonomo and Della Marca, 2020).

Values for and k are based on the estimates provided in (Buonomo and Della Marca, 2020). Specifically, is given by the sum of two terms: i) the inflow term due to births, , where years−1 is the country-level birth rate (ISTAT, 2020) and ii) the inflow term due to immigration, days−1, estimated from the number of citizens repatriated during about the first 11 weeks of the COVID-19 epidemic in Italy (Italian Ministry of Foreign Affairs and International Cooperation, 2020). Concerning a and k, numerical solutions by varying days−1 and are given in Section 5 (for detailed information about the ranges of values of the information parameters, see Buonomo and Della Marca, 2020).

Finally, to obtain an appropriate value for the baseline transmission rate , we consider model (3)–(14) in the absence of vaccination strategies ( days−1, D = 0) and search for the value that best fits with the initial ‘uncontrolled’ phase of the second epidemic wave in Italy. More precisely, we account for the number of COVID-19-induced deaths in Italy from 16 August, assumed as the starting date of the second wave (see Section 4.1 October 2020 (Italian Ministry of Health, 2020b). Indeed, on 13 October, the Council of Ministers approved a decree to reintroduce stricter rules to limit the spread of the disease (Italian Ministry of Health, 2020a). We use data regarding deaths which are more accurate compared with other ones. Anyway, by setting days−1, we obtain a good fit not only with the cumulative deaths (see Fig. 2 B) but also with the total infectious cases, (see Fig. 2A).

Fig. 2

Dynamics in the absence of vaccination ( days−1, D = 0). Total infectious cases (panel A) and cumulative disease-induced deaths (panel B) as predicted by model (3)–(14) (black lines) compared with Italian official data (Italian Ministry of Health, 2020b) (blue dots), in the period 16 August–13 October 2020. Initial conditions and other parameter values are given in Table 1.

All the parameters of the model as well as their baseline values are listed in Table 1.

Numerical simulations

Numerical simulations are performed in MATLAB (MATLAB, 2020). We use the 4th order Runge–Kutta method with constant step size for integrating the system and the platform-integrated functions for getting the plots.

First, we numerically investigate the impact of two vaccine-related parameters, namely the information-independent constant vaccination rate, , and the factor of vaccine ineffectiveness, , on the control reproduction number of formula (9). The corresponding contour plot of is shown in Fig. 3A , which shows the following: i) for very small values of , this parameter impacts , but days−1 suggests that depends practically only on in a linear-affine manner as shown in Fig. 3B; ii) for small values of (as those declared for some of the vaccines), the is small, for example, for , we have ; iii) for values of , comparable with those observed often for vaccines against the seasonal flu, we have ; and iv) if we define as the threshold of non-effectiveness the curve , we observe that for days−1, this threshold is reached for values of between around 0.6 and 0.7.

Fig. 3

Panel A: Contour plot of the control reproduction number (9) versus the information-independent constant vaccination rate, , and the factor of vaccine ineffectiveness, . Intersection of dotted black lines indicates the value corresponding to the baseline scenario: days−1, σ = 0.2. Panel B: Plot of versus , by setting days−1 (black line) and days−1 (blue line). Other parameters’ values are given in Table 1.

Temporal dynamics

Let us consider the time frame , where . We introduce four relevant cumulative quantities that will be used in the following: the cumulative vaccinated individuals CV, i.e., the total number of individuals who are vaccinated in ; the cumulative symptomatic cases CY, namely the number of new cases showing symptoms in ; the cumulative incidence CI, i.e., the total number of new cases in ; and the cumulative deaths CD, that is, the disease-induced deaths in . For model (3)–(14), we respectively have,

We also consider two possibilities for the time at which vaccination starts, namely.

VAX-0, that is, the baseline case that the vaccination campaign starts on day ;

VAX-30, that is, the case that the vaccination campaign starts on day .

We assume that in both cases the vaccination campaign lasts 1 year, i.e., [resp. ] days in the case VAX-0 [resp. VAX-30]. Other simulations corresponding to various starting times of the vaccination campaign (VAX-60, VAX-90, and VAX-120) are summarised in Table S.1 of supplementary materials.

Numerical simulations for the case VAX-0 are displayed in Fig. 4 . In other words, we report the temporal dynamics of three relevant state variables: susceptible individuals S (Fig. 4A), vaccinated individuals V (Fig. 4A), and symptomatic infectious individuals (Fig. 4C), as well as the cumulative number of deaths CD (Fig. 4D). We consider the following four significant scenarios (for each, we also report the observed results):

Fig. 4

VAX-0 case. Temporal dynamics of susceptible individuals S (panel A), vaccinated individuals V (panel B), symptomatic infectious individuals (panel C), and cumulative deaths CD (panel D), as predicted by model (3)–(14). Blue lines: constant vaccination with days−1, D = 0; black lines: information–dependent vaccination with days−1, D = 500μ/Λ; red lines: constant vaccination with ; green lines: constant vaccination with . Initial conditions and other parameter values are given in Table 1 and Section 5.1.

Constant vaccination (), with baseline rate days−1 (blue lines). At days, we observe the occurrence of a large peak of symptomatic cases (225,025), and a large cumulative number of deaths (28,343) at the end of simulation;

Information-dependent vaccination: days−1 and (black lines). This case is characterised by a time of peak that is halved w.r.t. the constant baseline case, at days about, and a much lower prevalence of , i.e., about one quarter w.r.t. the constant baseline case. This could be an excellent performance, but it is not the case since better performance could have led to higher vaccination rate levels;

Constant vaccination (), with rate days−1 (red lines), which is such that the peak value of is equal to the peak value observed in the case of information-dependent vaccination. One can note that, in this case, the epidemic peak occurs earlier, at days, and the final cumulative number of deaths is smaller: CD;

Constant vaccination (), with rate days−1 (green lines), where the peak of is halved w.r.t. the case of information-dependent vaccination. The epidemic peak occurs very early, at days, and the final cumulative number of deaths is relatively modest: CD.

Simulations for the case VAX-30 are, of course, graphically similar to those in Fig. 4; hence, corresponding plots are here omitted. From a quantitative point of view, to compare the results in the case VAX-30 with the case VAX-0, we focus on the scenario of information-dependent vaccination and provide in Table 2 the value of the following epidemiological indicators (not necessarily in this order): the number of susceptible and vaccinated individuals, the cumulative quantities (15) at the end of the time horizon , and the peak of symptomatic cases and its occurrence time. Comparison between the cases VAX-0 and VAX-30 is made as follows:where (see third column in Table 2).

Table 2

Information-dependent vaccination case ( days−1, ). Relevant quantities as predicted by model (3)–(14) in the case that the vaccination campaign starts at day 0, VAX-0 (first column), and in the case that it starts at day 30, VAX-30 (second column). The third column lists the differences between the values corresponding to the VAX-30 case compared with the case VAX-0. Initial conditions and other parameter values are given in Table 1.

X	X_VAX-0	X_VAX-30	X_VAX-30-X_VAX-0
S(t_f)	1.33·10^7	1.09·10^7	-2.39·10^6
V(t_f)	4.58·10^7	4.77·10^7	1.92·10^6
CV(t_f)	4.62·10^7	4.82·10^7	2.00·10^6
max(I_s)	5.76·10^4	9.14·10^4	3.38·10^4
argmax(I_s)	105.14	110.61	5.47
CY(t_f)	4.42·10^5	6.44·10^5	2.02·10^5
CI(t_f)	1.03·10^6	1.51·10^6	4.80·10^5
CD(t_f)	5.04·10^3	7.30·10^3	2.26·10^3

Note that, in both the VAX-0 and VAX-30 cases, cumulative asymptomatic people at the final time (that is, the difference CI) account for approximately 57% of the cumulative SARS-CoV-2 infections. This result is in line with the current estimates (as of April 2021) reported by the Italian National Institute of Health (ISS, 2020).

We also investigated the temporal dynamics of the ratio in the case of information-dependent vaccination. Numerical solutions are displayed in Fig. 5 for both the cases VAX-0 (black line) and VAX-30 (blue line). We note that in the case VAX-30, the ratio is larger than that in the case VAX-0 since the delay in starting the vaccination campaign induces a larger epidemic peak. In other words, in the case VAX-0, the maximum value reached by is 2.49, and the time at which it is reached is approximately days. In the case of VAX-30, this peak is reached at , i.e., 84 days after the start of VAX-30, but the peak value is much larger: .

Fig. 5

Information-dependent vaccination case ( days−1, D = 500μ/Λ). Temporal dynamics of the ratio between the information-dependent component, , and the constant component, , of the vaccination rate. Black line: VAX-0 case; blue line: VAX-30 case. Initial conditions and other parameter values are given in Table 1.

Sensitivity of epidemiological indicators to critical parameters

Here, we focus on the VAX-0 case and evaluate the sensitivity of some relevant epidemiological indicators to variations of critical parameter values. Note that for the case VAX-30, we obtain similar results, which are omitted.

Specifically, we assess how changing suitable model parameters affects the cumulative quantities (15) evaluated at the final time , the peak of symptomatic cases and its occurrence time. We anticipate here that the final cumulative incidence, CI, the final cumulative symptomatic cases, CY, and the peak of symptomatic cases, , have in all cases contour plots qualitatively similar to that of the final cumulative deaths, CD; thus, we do not plot them. Hence, the following figures display the counter plots of just three quantities:

the cumulative vaccinated individuals at days, CV;

the occurrence time of the symptomatic prevalence peak, arg;

the cumulative disease-induced deaths at days, CD.

We start by investigating how the information parameters, namely the information coverage, k, and the information delay, , may affect the epidemic’s course, see Fig. 6 . We observe that for arg and CD (as well as , CI and CY), the patterns of the contour plots are similar; in particular, for small , the range of the simulated variable when increases is large, whereas for , the range is restricted and low. The inverse phenomenon is observed for CV, where the range is restricted and small for low , while it is larger for .

Fig. 6

Impact of the information coverage, k, and the average delay, , on the VAX-0 scenario as depicted by contour plots. Panel A: cumulative vaccinated individuals at the final time days, CV. Panel B: time of symptomatic prevalence peak, arg. Panel C: cumulative deaths at the final time days, CD. The intersection of dotted white lines indicates the values corresponding to the baseline scenario: and days. Initial conditions and other parameter values are given in Table 1.

Then, we investigate how the factor of vaccine ineffectiveness, , and the information-independent constant vaccination rate, , affect the same quantities considered above. The results are shown in the contour plots in Fig. 7 for the case of constant baseline vaccination ( days−1 and ) and in Fig. 8 for the case of information-dependent vaccination ( days−1 and ). We may observe that the quantitative impact of the information-dependent vaccination is remarkable (but this was expected). As for the shapes of the plots, we note that the plots for CV (panels A) and for the time at symptomatic prevalence peaks (panels B) are remarkably different from the other plots. Moreover, the plot for CV is qualitatively different in the information-dependent vaccination case compared with the case of constant vaccination.

Fig. 7

Impact of the factor of vaccine ineffectiveness, , and the information-independent constant vaccination rate, , on the scenario VAX-0 with constant vaccination (i.e., ) as illustrated by contour plots. Panel A: cumulative vaccinated individuals at the final time days, CV. Panel B: time at symptomatic prevalence peak, arg. Panel C: cumulative deaths at the final time days, CD. The intersection of dotted white lines indicates the values corresponding to the baseline scenario: and days−1. Initial conditions and other parameter values are given in Table 1.

Fig. 8

Impact of the factor of vaccine ineffectiveness, , and the information-independent constant vaccination rate, , on the scenario VAX-0 with information–dependent vaccination (i.e., ), as shown by contour plots. Panel A: cumulative vaccinated individuals at the final time days, CV. Panel B: time when symptomatic prevalence peak is reached, arg. Panel C: cumulative deaths at the final time days, CD. The intersection of dotted white lines indicates the values corresponding to the baseline scenario: and days−1. Initial conditions and other parameter values are given in Table 1.

As discussed in Section 4.2, the baseline value of the transmission rate has been estimated using the official data from 16 August to 13 October 2020 released by the Italian authorities. Therefore, it takes into account the non-pharmaceutical interventions, such social distancing and quarantine, that were in place during that time frame. However, it is reasonable to assume that vaccination campaigns may be associated with a relaxation of non-pharmaceutical interventions.

To empirically assess the impact of this phenomenon on the epidemiological indicators, we performed additional simulations (Fig. 9, Fig. 10 ) where ranges between the value , corresponding to the threshold value , and , corresponding to , which represents the unperturbed case (this is the value of the basic reproduction number estimated at a very early stage of the epidemic in Italy (Gatto et al., 2020). Fig. 9 shows that by keeping the information coverage constant, the following is observed: i) both the cumulative number of vaccinated individuals and the time of the epidemic peak non-monotonically depend on ; ii) the cumulative number of deaths, instead, increases with . As for the interplay between and the average information delay, Fig. 9 presents features qualitatively similar to those in Fig. 9.

Fig. 9

Impact of the information coverage, k, and the transmission rate, , on the VAX-0 scenario as shown by contour plots. Panel A: cumulative vaccinated individuals at the final time days, CV. Panel B: time of symptomatic prevalence peak, arg. Panel C: cumulative deaths at the final time days, CD. The intersection of dotted white lines indicates the values corresponding to the baseline scenario: and days−1. Initial conditions and other parameter values are given in Table 1.

Fig. 10

Impact of the information delay, , and the transmission rate, , on the VAX-0 scenario as shown by contour plots. Panel A: cumulative vaccinated individuals at the final time days, CV. Panel B: time of symptomatic prevalence peak, arg. Panel C: cumulative deaths at the final time days, CD. The intersection of dotted white lines indicates the values corresponding to the baseline scenario: days and days−1. Initial conditions and other parameter values are given in Table 1.

The impact of seasonality

There is an ongoing debate on the possible effects of seasonality on the transmission and global burden of COVID-19 (Merow and Urban, 2020, Liu et al., 2021, Sajadi et al., 2020, Audi et al., 2020). Thus, for the sake of the completeness, we consider in this study the case of information-dependent vaccination and simulate the presence of seasonality on three key parameters: the transmission rate, , the rate of symptoms’ onset, , and the total rate of vaccination, , with given in (14). Indeed, for the latter the seasonality is due to a lower vaccination rate during the summer vacations.

Therefore, we use in our simulationswhere are the baseline values, and is simply two-states switch, similar to the one proposed in (Earn et al., 2000) for the transmission rate:

Since we used initial conditions corresponding to the COVID-19 data as of 16 August 2020, as officially communicated by Italian health authorities (see Section 4.1), we consider

We will denote this simulation scenario with VAX-0S.

Numerical simulations are displayed in Fig. 11 and compared with the baseline scenario, VAX-0. Corresponding relevant quantities are provided in Table 3 . Our simulation suggests that i) the impact of the summer vacation on vaccine delivery and is minimal (so they are omitted from Fig. 11); ii) despite the epidemic peak occurs many months after the summer, the peak of symptomatic cases is delayed compared with the no-seasonality case (Fig. 11A); iii) the cumulative number of deaths decreases slightly (Fig. 11B).

Fig. 11

Impact of the seasonality on the information–dependent vaccination case ( days−1, D = 500μ/Λ). Temporal dynamics of symptomatic infectious individuals (panel A), and cumulative deaths CD (panel B), as predicted by model (3)–(14). Blue lines: VAX-0S case (i.e., scenario including seasonality); black lines: VAX-0 case (i.e., no–seasonality scenario). Initial conditions and other parameter values are given in Table 1 and Section 6.

Table 3

Information–dependent vaccination case ( days−1, ). Relevant quantities as predicted by model (3)–(14) in the scenario including seasonality VAX-0S (first column). The second column lists the differences between the values corresponding to the VAX-0S case compared with the case VAX-0 (see also Table 2). Initial conditions and other parameter values are given in Table 1 and in Section 6.

X	X_VAX-0S	X_VAX-0S-X_VAX-0
S(t_f)	1.45·10^7	1.18·10^6
V(t_f)	4.47·10^7	-1.09·10^6
CV(t_f)	4.51·10^7	-1.12·10^6
max(I_s)	5.03·10^4	-7.34·10^3
argmax(I_s)	115.66	10.52
CY(t_f)	4.02·10^5	-3.97·10^4
CI(t_f)	9.44·10^5	-8.91·10^4
CD(t_f)	4.59·10^3	-444.88

Conclusions

In this paper, we introduced a mathematical model describing the transmission of COVID-19 in the presence of non-mandatory vaccination. The main novelty is that the vaccine hesitancy and refusal is taken into account. To this end, we used the information index, which mimics the idea that individuals take their decision on vaccination based on not only the present but also the past information they have on the spread of the disease.

Theoretical analysis and simulations clearly show that voluntary vaccination can certainly reduce the impact of the disease but is unable to eliminate it. Qualitatively, the time evolution of the disease remains the same, but the quantitative results are very different: an epidemic outbreak (a new epidemic wave) occurs even if (as we observed in our simulations) the information-dependent vaccination rate is, at its peak, more than three times larger than the constant baseline vaccination rate.

A key result is in particular the fact that the information-related parameters greatly affect the dynamics of the disease: large information coverage and small memory characteristic time are needed to have the best results. The different impact of behaviour and information with respect to the scenario of mandatory constant vaccination can be further appreciated by examining the contour plots in Fig. 6, Fig. 7, Fig. 8, Fig. 9, Fig. 10.

As it is reasonable, the parameter , i.e., the risk of infection for vaccinated people, has a major impact. That is, the control reproduction number essentially depends on in a linear-affine manner. This suggests to stick to vaccines that have very low , where is tiny. A very positive result is that the threshold of non-efficacy of the vaccine, which can roughly be delineated as the curve , where is located for values , i.e., for very large values of (Fig. 3A).

As far as the impact of human behaviour w.r.t. scenarios with constant vaccination rates is concerned, we found that the performances were better only w.r.t. a constant vaccination rate as low as , whereas the scenario where (see Section 5.1) led to excellent results and a substantially smaller number of deaths.

As far as the comparison between the VAX-0 and VAX-30 scenarios is concerned, we also measured the impact of VAX-30 scenario on the ratio between the information-dependent and the constant components of the vaccination rate, namely . As expected, the peak was considerably higher in the scenario VAX-30. The peaks occur in the same week if measured in the absolute time, i.e., the peak for VAX-30 occurs one month before the peaks of VAX-0 if measured in time since the start of the vaccination (see Fig. 5).

Finally, seasonality has a relative but non-negligible relevance. For example, although the decrease of the transmission rate and the onset of symptoms occurs in the summer, the predicted winter epidemic peak of symptomatic cases is decreased and delayed compared with that in the no-seasonality scenario. A small but non-negligible decrease and delay of the cumulative deaths is also observed. This overall suggests that a decrease of the transmission and of the onset of symptoms has positive impact even many months after their end (see Fig. 11).

An apparent limitation of this study is the absence of modelling for the dynamics of the transmission rate. In other words, neither the spontaneous changes of the parameter nor the imposed changes due to social distancing laws and partial/full lockdowns are taken into the account. However, these aspects are intentionally neglected here since our goal is to assess the impact of a possible voluntary vaccination campaign.

Another limitation of our work is that its aim is merely to investigate the impact of the spontaneous behavioural response of population. However, the influence of government policy, modelled as external control, may be crucial for the success of vaccination campaigns and should be the subject of future investigations. A possible approach could be the one adopted in (Buonomo et al., 2019) in the case of endemic childhood diseases where the optimal control theory has been employed.

Finally, we also plan to i) explore (mainly numerically) a realistic model of COVID-19 spread that includes the time-changes of the transmission rate and ii) assess the possibility that eradication of COVID-19 is not reached and the disease stays endemic.

CRediT authorship contribution statement

Bruno Buonomo: Conceptualization, Methodology, Supervision, Writing - original draft, Writing - review & editing. Rossella Della Marca: Methodology, Formal analysis, Software, Writing - original draft, Writing - review & editing, Funding acquisition. Alberto d’Onofrio: Conceptualization, Methodology, Supervision, Validation, Writing - original draft, Writing - review & editing. Maria Groppi: Conceptualization, Methodology, Writing - original draft, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.