ABSTRACT
Throughout the last two centuries, vaccines have been helpful in mitigating numerous epidemic diseases. However, vaccine hesitancy has been identified as a substantial obstacle in healthcare management. We examined the epidemiological dynamics of an emerging infection under vaccination using an SVEIR model with differential morbidity. We mathematically analyzed the model, derived , and provided a complete analysis of the bifurcation at . Sensitivity analysis and numerical simulations were used to quantify the tradeoffs between vaccine efficacy and vaccine hesitancy on reducing the disease burden. Our results indicated that if the percentage of the population hesitant about taking the vaccine is 10%, then a vaccine with 94% efficacy is required to reduce the peak of infections by 40%. If 60% of the population is reluctant about being vaccinated, then even a perfect vaccine will not be able to reduce the peak of infections by 40%.
1. Introduction
Vaccination against infectious diseases is considered one of the most successful achievements of public health management. The world's first vaccine was demonstrated in 1796 in the UK by Dr. Edward Jenner against smallpox [Citation1]. Since then, vaccines have been developed and successfully used in controlling and eliminating various communicable diseases including Polio, Rubella, Measles, Mumps, and Chickenpox [Citation2]. Some vaccines are effective at blocking infection whereas other vaccines reduce the development of disease [Citation3]. According to the CDC, mRNA vaccines for SARS-CoV-2 could potentially reduce the risk of severe symptoms by 90% [Citation4].
Despite the successes of vaccination, however, public reluctance to inoculation has become a considerable barrier to disease control. The WHO defines vaccine hesitancy as the reluctance or refusal to vaccinate despite the availability of vaccines [Citation5]. This phenomenon can be particularly acute in the case of emergence of a new infection, as was observed during the COVID-19 pandemic [Citation6]. Potential reasons for vaccine hesitancy include mistrust in the health care system, religious or cultural beliefs, concerns about vaccine safety, costs, and logistical barriers, among others [Citation7,Citation8]. While vaccine hesitancy could lead to decreased levels of vaccine coverage and a heightened risk of disease outbreaks, experts worldwide agree that vaccine hesitancy exhibits an increasing trend [Citation7]. However, vaccine efficacy has proven to have had a substantial effect on disease control in the past [Citation4,Citation9]. Therefore, this study seeks to examine the impact of vaccine efficacy and vaccine hesitancy, as well as their interplay in the control of an emerging infectious disease.
Despite the importance of vaccination in disease control, the number of ODE-based mathematical studies that consider the impact of vaccine hesitancy is modest [Citation6,Citation10–15]. Using a compartmental model with ODEs, Mesa et al. examine how the refusal of measles vaccines could lead to outbreaks and higher societal costs [Citation10]. They use data on measles in England as a case study and find that even low levels of vaccine reluctance could cause a significant societal burden.
Buonomo et al. use an extended SVEIR model of COVID-19 to explore individuals' decision making on taking vaccines [Citation6]. They use numerical simulations to demonstrate how present and past vaccine-related information influences vaccination decisions, and in turn, the dynamics of disease spread. Leung et al. use an SIRSV model of COVID-19 to explore the impact of vaccine hesitancy on disease dynamics [Citation11]. Analytical and numerical methods demonstrate that the inclusion of vaccine hesitance in the model avoids significant biases in the process of determining the stability condition, maximum peak infection time, and threshold behaviours of their system.
In the current study, we aim to examine the impact of vaccine hesitancy on infections and deaths due to the emergence of a novel pathogen, and we ask what vaccine efficacy would be required to overcome the brunt of vaccine hesitancy. Furthermore, we investigate whether infected individuals with mild, moderate, or severe symptoms drive the disease dynamics.
The model we developed is a modified SEIR compartmental model that includes vaccination and differential morbidity and is comprised of seven compartments: susceptible, vaccinated, exposed, infected (mild symptoms), infected (moderate symptoms), infected (severe symptoms), and recovered. Our model is unique in that we consider three infected categories (mild, moderate, and severe symptoms), modeled with different transmission coefficients and recovery rates, and we include the consideration of vaccine hesitancy. The effects of vaccine hesitancy are explored by expressing vaccine coverage as a function of vaccine hesitancy, vaccine availability, and the fraction of susceptible individuals in the community. In this study, we perform a detailed mathematical analysis and use numerical simulations to answer the following research questions:
How does increasing the vaccine efficacy and vaccine coverage help reduce cases and deaths due to the emerging infection?
What is the impact of vaccine hesitancy on infections and deaths, and what efficacy of a vaccine would be required to overcome the brunt of vaccine hesitancy?
Which model parameters are crucial in controlling the disease burden, and how does change with variation in these parameters?
Which class of infected individuals (i.e. mild, moderate, or severe) contributes most to the overall infection burden?
The organization of this paper is given below. An introduction to the study along with a brief review of literature is provided in Section 1 (current section). The model development is comprehensively explained in Section 2. Section 3 delves into the mathematical analysis of the model. A detailed numerical analysis of the model is presented in Section 4. Finally, a discussion of the results is given in Section 5.
2. Model development
In this study, we use a deterministic compartmental model in order to examine the dynamics of a respiratory disease caused by an emerging infectious pathogen, in the presence of vaccination that prevents infection. The model divides the total population in a hypothetical community into seven mutually exclusive compartments based on the disease status of each individual at any time t.
Unvaccinated susceptible (S)
Vaccinated susceptible (V)
Exposed (E)
Infected with mild symptoms ()
Infected with moderate symptoms ()
Infected with severe symptoms ()
Recovered (R)
This compartmental model tracks the flow of individuals into and out of the seven states over time. The schematic diagram pertaining to the compartmental model is given in Figure .
The following system of non-linear autonomous ordinary differential equations represents the inflow and outflow of individuals to and from the compartments and community over time. (1) (1) The descriptions of the parameters used in the model are presented in Table . The parameter values used in the numerical simulations are given in Table .
A few examples of previous models with differential morbidity in the context of COVID-19 can be found in the literature [Citation3,Citation20,Citation22]. In particular, Šušteršić et al. [Citation20] includes a classification of infected individuals on the spectrum of mild, severe, and critical symptoms. Moreover, some vaccines may only block infections, whereas others would only provide protection against severe disease, but not infection [Citation3,Citation23,Citation24]. In particular, Beukenhorst et al. [Citation3] focus on understanding which type of vaccination strategies would be more effective at a population level.
Key assumptions used in developing the model are as follows:
Demography. Susceptible (S) individuals are recruited at a rate θ. The natural community mortality rate due to factors excluding the disease is denoted by μ.
Infection. The per-capita force of infection is the sum , where through are the infected compartments, and through are transmission coefficients. Transmission occurs by mass action; the per-time rate of influx to the exposed (E) compartment due to infection of susceptible individuals is .
Leaky vaccine. The vaccine is imperfect (or leaky), with efficacy , understood as the percentage chance that the vaccine will successfully protect an individual from infection. Susceptible individuals (S) are vaccinated at per capita vaccine coverage rate p, at which point they enter the vaccinated (V) compartment. Vaccinated individuals can become exposed (E), in the event of vaccine failure, at per-time rate . A leaky vaccine with similar traits can be found in Esteban & Almodovar-Abreu [Citation17].
Disease progression. Exposed individuals (E) become infected with mild symptoms () after an incubation period of days. Individuals infected with mild symptoms () can progress to the state of moderate symptoms () according to the rate or recover at the rate . Once acquiring moderate symptoms, individuals can progress to severe symptoms () at the rate , recover at the rate , or die due to disease at the rate . Individuals with severe symptoms can recover at the rate or die due to disease at the rate .
Severity of infection. Individuals infected with mild symptoms are more likely to be in the community spreading the infection as opposed to those infected with moderate and severe symptoms, who may self-isolate or be quarantined or hospitalized. We therefore assume . We also assume that more severe infections are more deadly, and therefore . Finally, we assume that more severe infections require longer to clear, and so we assume .
Single viral strain. Due to the complexities of disease spread caused by the coexistence of multiple viral variants, the current model is limited to a single variant. Since we do not consider the generation of mutant variants that are most likely responsible for reinfection, we assume that recovery from the disease induces long-term immunity, which excludes the possibility of reinfections. Previous studies exclude the waning of disease-induced and/or vaccine-induced immunity for similar reasons [Citation17,Citation18,Citation20,Citation25–27].
2.1. Functional form of vaccine coverage in consideration of vaccine hesitancy
In this work, we consider the model (Equation1(1) (1) ) with two different choices for the vaccination coverage rate p. The simpler of the two is when it is treated as constant; see Section 3 for mathematical analysis and Sections 4.1 and 4.2 for numerical analysis. In the alternative, we focus on the impact of vaccine hesitancy and define a time-dependent vaccine coverage rate. We derive the functional form for the coverage rate as follows.
We assume that the per-capita vaccination rate is proportional, with proportionality (rate) constant k, to the probability that a randomly selected member of the population is both susceptible and vaccine-approving (i.e. not hesitant). Assume the probability of a randomly sampled member of the population being vaccine-hesitant is ψ, so that is the probability that the individual is vaccine-approving. The probability of a randomly sampled member of the population being susceptible is , where is the total community population. Assuming these sampling events are independent, the probability of a random individual being susceptible and vaccine-approving is . Using the definition of the proportionality (rate) constant k, we arrive at the formula (2) (2) The simulations in Sections 4.3 and 4.4 are based on the functional form of p given in (Equation2(2) (2) ).
3. Mathematical analysis of the model
This section is organized as follows. In Section 3.1, we present some elementary mathematical properties of the model. The disease-free equilibrium, basic reproduction number, and its threshold properties are proven in Section 3.2. In Section 3.3, we discuss the details surrounding the endemic equilibrium. The bifurcation at is analyzed in Section 3.4.
We will frequently make use of a collection of aggregate loss rate parameters: which represent the total loss rates from the infectious compartments and , respectively. They will allow for more compact expressions in the subsequent sections.
3.1. Preliminary results
Since the model (Equation1(1) (1) ) is defined by a polynomial vector field, local existence and uniqueness of solutions are guaranteed. We will therefore only focus on non-negativity and boundedness properties.
Proposition 3.1
Non-negativity and Boundedness of Solutions
All solutions of system (Equation1(1) (1) ) with non-negative initial conditions remain bounded and non-negative for all .
Proof.
We show that no component of the model can leave the positive orthant. First, suppose that , but there exists a time such that . Let be minimal; that is, . Then , which implies for all s>0 sufficiently small. This contradicts the definition of . Therefore, for all where it is defined. Similarly, if , then for all t>0.
Given the above analysis, the second equation of system (Equation1(1) (1) ) leads to as for all . Solving this differential inequality, we get Similarly, by the third equation of system (Equation1(1) (1) ), we obtain as . Solving this differential inequality yields . Similar calculations show that and R remain non-negative for all where they are defined.
Adding all seven equations of system (Equation1(1) (1) ) together and using the non-negative nature of each component, we obtain the inequality . It follows that and . Therefore, solutions to system (Equation1(1) (1) ) are bounded.
Remark: It is important to verify whether the state variables of system (Equation1(1) (1) ) remain positive and bounded for all time in order for the system to be epidemiologically relevant.
Proposition 3.2
Positive Invariance
The closed set is a positively invariant region for system (Equation1(1) (1) ).
Proof.
Let denote the sum of all components, and suppose so that the initial condition is in Ω. From the proof of Proposition 3.1, we know that
3.2. Disease-free equilibrium and basic reproduction number
The disease-free equilibrium (DFE), also referred to as the infection-free equilibrium, of an epidemiological ODE system is the equilibrium solution (i.e. stationary point) where all the disease related variables are equal to 0. A routine calculation yields (3) (3) The basic reproduction number or basic reproductive ratio, usually denoted by , gives valuable insights into the dynamics of an epidemic, as it delineates if an outbreak will occur or if the disease will die out. is biologically defined as the expected number of secondary infections produced by a single infectious case in an otherwise completely susceptible population [Citation28]. Mathematically, it is a threshold quantity that determines the stability of the disease-free equilibrium [Citation28,Citation29]. The disease spreads (i.e. an epidemic occurs) if , and it dies out if . The following theorem proves exactly this result for the case of a small outbreak. The stability of the DFE when is treated in Lemma 3.5.
Theorem 3.3
The basic reproduction number for model (Equation1(1) (1) ) is (4) (4) The disease-free equilibrium is locally asymptotically stable if and only if , and unstable if and only if .
Proof.
To this end, we use the next-generation matrix (NGM) method developed by Diekmann et al. as well as van den Driessche and Watmough [Citation29,Citation30]. First, we divide the model compartments into two main categories as uninfected, , and infected, . We rewrite as , where represents all rates from X to Y, and represents all other rates.
The Jacobians of and evaluated at the DFE given in (Equation3(3) (3) ) are and is the spectral radius of the next-generation matrix , and is given by (5) (5) Then, by simplifying further, we obtain a simple expression, which we readily find is equal to (Equation4(4) (4) ). The stability assertions follow from Theorem 2 of [Citation29].
Remark 3.1
The expression (Equation5(5) (5) ) has a succinct biological interpretation. It is the sum of rates at which individuals become exposed (not yet infectious) by individuals in each of the three infectious states, weighted by the amount of time a typical individual spends in each of these states as they progress through their infection. Specifically, the unconditional (i.e. without assuming transfers between states) average length of time an individual spends with mild symptoms (), moderate symptoms (), and severe symptoms () is , , and , respectively. The probability of the transition is the ratio , so that the average time spent with mild symptoms after exposure is . The probability of transition from is . Therefore, the probability of the transition is , so that the average amount of time spent with moderate symptoms (after exposure and mild symptoms) is . Similarly, the average amount of time spent with severe symptoms is
Remark 3.2
If , the above theorem guarantees only that a sufficiently small outbreak will die out. While numerical simulations strongly suggest that global stability of the DFE () follows from , we make no attempt to mathematically prove that in this work.
3.3. Endemic equilibrium
The endemic equilibrium of an epidemiological ODE system is an equilibrium solution to the system where the disease-related variable(s) are not necessarily equal to 0. For this model, an explicit expression is available for the endemic equilibrium (it can be shown to be equivalent to one of the zeroes of a second degree polynomial), but the resulting formula is not especially useful due to the introduction of several very complicated bulk parameters. We instead analyze the endemic equilibrium using a more geometric approach, which can be connected to .
Lemma 3.4
The model (Equation1(1) (1) ) has at most one (positive) endemic equilibrium. In particular, the endemic equilibrium exists if and only if .
Proof.
Notice that any steady state of system (Equation1(1) (1) ) must satisfy the following equations: Recall that . Hence, we obtain where . Substituting this into the expression for E, we get Suppose , so that we eliminate the disease-free equilibrium. Canceling E on both sides and substituting in S and V, we have Now define by . Thus, the endemic equilibria are positive solutions of the non-linear equation . Clearly, the function is continuous and non-negative on its domain. Moreover, we have for all . Therefore, is a strictly decreasing function on its domain. In addition, we have . Thus, a positive solution of the equation exists if and only if . Furthermore, this solution will be unique due to the monotonicity of the function . Also, no endemic equilibrium exists if . Straightforward algebraic simplifications can be used to show that .
Lemma 3.4 characterizes the existence of a unique endemic equilibrium in terms of . Figure shows a plot of the endemic equilibrium curve of model (Equation1(1) (1) ).
3.4. Bifurcation at
Bifurcation analysis concerns the qualitative changes to the stability of a dynamical system's equilibria with respect to the variation of a system parameter. We conduct an analysis based on the centre manifold theory in order to investigate the nature of the bifurcation at in system (Equation1(1) (1) ). The theory of bifurcations in an epidemiological context has been extensively discussed [Citation31–33].
Lemma 3.5
The model (Equation1(1) (1) ) has a transcritical bifurcation with respect to ϵ at . Namely, the following points hold for small:
When , there are no (positive) endemic equilibria near the disease-free equilibrium, and the latter is locally asymptotically stable;
When , there is a unique endemic equilibrium near the disease-free equilibrium, and it is locally asymptotically stable. The disease-free equilibrium is unstable.
Proof.
Note that can be written as for a constant K>0 that does not depend on ϵ. We will take advantage of this fact and define a change of parameters. Namely, solving the above expression for ϵ, we get The basic reproduction number is now an explicit parameter in the ODEs, and the impact is that the differential equations for E and now read as In this new parameterization, the Jacobian at the disease-free equilibrium has a single zero eigenvalue when , and all other eigenvalues have a negative real part.
In what follows, we will neglect the recovered (R) component from the model, since it is not relevant to the qualitative dynamics and represents a terminal state in the flow diagram. The Jacobian of system (Equation1(1) (1) ) at and is Following the framework of Theorem 4.1 of Castillo-Chavez and Song [Citation34], we require left and right eigenvectors v and w of the Jacobian. When , one can verify that suitable choices are where is a dimensionless quantity, denotes the row i column k entry of the Jacobian, and is a normalizing constant that ensures . Note that since , we have , so is positive, and it follows that ρ is well-defined and positive.
The next step is to form the sums where are the components of the vector field (Equation1(1) (1) ), and all partial derivatives are evaluated at and . The vector field is quadratic with few non-linear terms in the variables . The nonzero partial derivatives required to calculate a are for j = 1, 2, 3, and their symmetric versions. Taking this into account, together with the fact that , we have and we conclude that a<0, given that , and the other terms in parentheses are all positive.
For the b sum, the nonzero partial derivatives are Then, we have Since , , and , we conclude that b>0.
The dynamics on the centre manifold attached to the disease-free equilibrium are equivalent to which has a transcritical (i.e. forward) bifurcation with a nontrivial equilibrium . The first-order Taylor expansion of the resulting non-trivial (i.e. not disease-free) equilibrium is therefore The conclusions of the theorem follow by taking into account the sign pattern of the vector w and the signs of a<0 and b>0.
Remark 3.3
The previous lemma was proven by formally introducing into the differential equations as a parameter, thereby eliminating the incubation rate ϵ. However, similar (in fact, easier) proofs can be used to show that the conclusions of the lemma are the same if instead, one varies the vaccine efficacy ϕ, coverage rate p, or any of the transmission coefficients for i = 1, 2, 3. Establishing that the endemic equilibrium is stable for in full generality (i.e. not just for small) is beyond the scope of this work.
4. Numerical results
This section is organized as follows. In Section 4.1, we present a sensitivity analysis of with respect to the model parameters. Computational simulations of the model that address our primary research questions are provided in Sections 4.2–4.4. The software packages used for the simulations include R, MATLAB, and Mathematica.
4.1. Sensitivity analysis
In this section, we provide an analysis of the local sensitivity of the model to its parameters. We investigate the effect that a small change in a given parameter has on . In epidemiological modelling, sensitivity analysis of with respect to the parameters appearing in its expression provides important insights into disease control. For example, if a parameter has a high positive (or alternatively, negative) sensitivity index with respect to , reducing (or alternatively, increasing) the factors contributing to that parameter could help control the disease.
We conduct a sensitivity analysis of with respect to the model parameters using differential sensitivity analysis (i.e. direct method), which uses point values of the parameters. Given a model parameter ω, the normalized forward sensitivity index of with respect to ω is given by [Citation28]. Observe that we can rewrite the general sensitivity index formula as . From the latter structure, the sensitivity index formula essentially accounts for the ratio of the corresponding percentage change in resulting from a percentage change in the parameter ω. Given below are the sensitivity index formulae of the parameters listed in Table .
where and .
The sensitivity index of a parameter may depend on other parameters of the model, or it may be a constant value with no dependence on other parameters. For example, the sensitivity index of θ is independent of any parameter values. The expressions and are all unconditionally positive, so an increase in the parameters and always increase the value of . Conversely, the expressions and are all negative for all applicable values of the relevant parameters. Thus, an increase in the parameters and always decrease the value of . However, the sensitivity index expressions and have an ambiguous structure, so those parameters' correlation (positive or negative) on generally depends on the values of some or all other parameters.
The aforementioned ambiguity can be discerned by evaluating those sensitivity indices at the numerical values of the parameters, and then examining their signs. The sensitivity indices are given in Table and illustrated graphically in Figure . Notice that means that if the value of is increased (or decreased) by 10%, then the value of will consequently increase (or decrease) by 7.114%. Moreover, indicates that if the value of is increased (or decreased) by 10%, then it will result in a 7.686% decrease (or increase) in the value of . We note that has a remarkably high negative value compared to the other sensitivity indices. Defining a vaccine failure parameter, ζ, as , however, yields , which is in . Thus, we use the vaccine failure parameter (ζ) for the graphical illustration of sensitivity indices in Figure .
The parameters and ϵ have a positive influence on , while and have a negative influence on . Further, θ, μ, ϕ, , and are the most sensitive parameters in the model. However, θ and μ are demography parameters, and are not easily amenable to control measures. Thus, the remaining parameters ϕ, , and are of greatest interest in our study.
Some of the results from our sensitivity analysis were expected, since they follow from and are consistent with the model formulation. For example, we expected that an increase in recovery rates would lead to a decrease in infections, and hence , which is indeed the case according to the sensitivity analysis. However, our sensitivity analysis also highlighted a less obvious result, namely that takes a positive value, despite the expectation that a higher progression rate from mild to moderate symptoms would lead to a reduction of disease burden as .
4.2. Numerical simulations on vaccination (with constant coverage)
Numerical simulations were performed on the model in order to answer our research questions about disease dynamics. The initial values used for the state variables are . The basic reproduction number of the system, evaluated at the parameter values given in Table , is . Since , this indicates that an outbreak occurs thus causing an epidemic in the community.
Figure shows how the value of changes with different pairs of parameters. As expected, varies as transmission coefficient and the vaccine coverage rate p change, but it is much more sensitive to transmission than coverage, if the coverage is low (i.e. between 0.005 and 0.05) (Figure (a)). As long as the vaccine coverage is high enough, increases with increasing transmission, but the vaccine is still sufficient to prevent an epidemic that has a low to moderate transmission coefficient (for mild cases). If vaccine coverage is too low, the transmission coefficient does not matter and any value for (across the range shown) gives rise to an epidemic.
The sensitivity of to transmission coefficient and vaccine efficacy ϕ is shown in Figure (b); the greatest variation in results is observed for vaccine efficacies between 90% and 100%. Notice how stronger vaccine efficacies are able to increasingly prevent epidemics with higher values of , but lower efficacies are insufficient to control the spread of infection with more moderate or low coefficients of transmission. Looking at vaccine efficacy together with vaccine coverage (Figure (c)), notice that efficacies above 96% would be needed to prevent an epidemic across the entire displayed range of vaccine coverage rates (i.e. 0 to 20%).
Figure shows the cumulative number of disease-induced deaths over 200 days for different values of vaccine efficacy (. Notice that the number of deaths after 200 days could be reduced drastically by increasing vaccine efficacy (ϕ) from 90% to 100%. A vaccine with a 97.5% efficacy would save the lives of nearly twice as many infected individuals compared to a vaccine with an efficacy of 90%. Furthermore, it can be observed that a perfect vaccine would reduce the cumulative number of deaths by about 65% as compared to a vaccine which is 90% effective. Therefore, it becomes evident from Figure that high efficacy vaccines considerably reduce disease mortality, even when the vaccine is designed to block infection rather than to reduce disease severity.
4.3. Numerical simulations on vaccine hesitancy
We next use the functional form (Equation2(2) (2) ) of vaccine coverage to assess the impact of vaccine hesitancy on the spread of infection and the number of disease-induced deaths per day. Figure shows the number of infections (with mild, moderate, and severe symptoms) over time for values of vaccine hesitancy ranging from 0 to 100%. Infections decrease with declining vaccine hesitancy; the peak in mild infections, for example, drops in half (approximately) without vaccine hesitancy () as compared to 100% vaccine hesitancy (, i.e. no vaccination). The peaks of the infection curves shift very slightly to the right as ψ decreases, showing that the time to infection peak is not as rapid without vaccine hesitancy. Notice in Figure that individuals with mild symptoms contribute most to the infection burden. Due to exposed individuals' progression through the three infected stages, with a certain percentage of each infected class being removed from infection, the peak of infections decreases through the progression (which follows from the model formulation).
Figure illustrates the cumulative number of disease-induced deaths over 200 days for different values of vaccine hesitancy (ψ). Notice how the predicted cumulative deaths increase evenly as ψ increases from 0 to 100% by quarters. Lower vaccine hesitancy in the community can give rise to epidemics with fewer deaths; the simulation shows that three-quarters as many deaths result from the epidemic with as compared to that with .
4.4. Quantification of results on vaccine hesitancy
In this section, we use simulations to determine the vaccine efficacy needed to reduce the peak of infections by 40% (as compared to the peak in the absence of vaccination), given various levels of vaccine hesitancy. The peak is taken as the maximum value of , or (i.e. ), and the absence of vaccination is the case where , and hence p = 0. The value of vaccine efficacy ϕ is then calibrated accordingly (i.e. giving 40% of ) with levels of vaccine hesitancy of 10%, 25%, and 60%.
Results are shown in Table , revealing the considerable impact of vaccine hesitancy. If 60% of the population is reluctant about being vaccinated, then even a perfect vaccine will not be able to reduce the peaks of infections by 40%. It is interesting to note that even if the vaccine hesitancy is small (10%), the efficacy needs to be above 94% to reduce the peaks of infections by 40%.
Figure shows how the percentage reduction in the peak of total infections () changes as vaccine hesitancy varies from 0 to 100% and vaccine efficacy varies from 70% to 100%. These results suggest that if and , the peak of the total infections could be reduced by more than a quarter of the peak infections in the absence of vaccination. It is hence evident that considerable reductions in peak infections could be achieved with lower vaccine hesitancies and much higher vaccine efficacies.
5. Conclusions
In this study, we developed an SEIR model with vaccination and differential morbidity to gain insights into the epidemiological dynamics of an emerging infectious disease. We also constructed a functional form for vaccine coverage in terms of vaccine hesitancy. The resulting compartmental model of seven state variables expressed as a system of non-linear ODEs was analyzed, both mathematically and numerically. The objectives of the study were to identify how vaccine efficacy and vaccine hesitancy affect disease dynamics; to estimate the efficacy of a vaccine needed to overcome the impact of vaccine hesitancy; to explore how differential morbidity impacts disease dynamics; and to determine which parameters play a pivotal role in an outbreak of disease.
Firstly, we proved some important results on the qualitative nature of the mathematical model. We computed the basic reproduction number () of the model using the NGM method. We also proved that the model has a unique endemic equilibrium when and no endemic equilibria when , despite the quadratic equation of endemic equilibria. Then, we performed a bifurcation analysis in order to show that there exists a transcritical bifurcation with respect to ϵ at .
The numerical analysis of the model involves a sensitivity analysis, graphical simulations, and a quantification of results in order to address the research questions, including the tradeoff between vaccine efficacy and vaccine hesitancy. The sensitivity analysis showed that the parameters that are most sensitive to and easily amenable to control measures are the vaccine efficacy (ϕ) (or vaccine failure (ζ)), transmission coefficient for individuals infected with mild symptoms (), and recovery rate for individuals infected with mild symptoms (). The sensitivity analysis and Figure indicate that milder infections contribute more to the disease burden; this result follows from the model formulation. Figures and reinforce that more efficacious vaccines could help control the disease burden more efficiently.
In order to examine the impact of vaccine hesitancy on the spread of disease, we performed numerical simulations based on the functional form of p. Figures – indicate how the disease burden and disease-induced deaths could be decreased with lower vaccine hesitancy. Through the quantification of our results, we established that if the vaccine is not perfect, a lower vaccine hesitancy would be needed to reduce the peaks of infections to the same degree. Thus, it could be inferred that taking steps to overcome vaccine hesitancy could help substantially reduce the disease burden. These results could prove to be useful for healthcare management.
In summary, we used a compartmental epidemiological model of an emerging pathogen to show how the disease burden could be reduced substantially through effective vaccines and reduced vaccine hesitancy levels. Future research in this direction could focus on considering the effects of multiple mutant variants of the pathogen, which could lead to the waning of both vaccine-induced and infection-acquired immunity. Moreover, considering the vaccinated but infected individuals as a separate compartment in the model would enable evaluation of the vaccine efficacy against transmission separately from the vaccine efficacy against severe infections, hospitalizations, and deaths.
Acknowledgments
We would like to thank two anonymous reviewers for feedback that greatly improved our work.
Disclosure statement
No potential conflict of interest was reported by the authors.
Additional information
Funding
References
- World Health Organization. History of smallpox vaccination; n.d. Available from: https://www.who.int/news-room/spotlight/history-of-vaccination/history-of-smallpox-vaccination.
- Centers for Disease Control & Prevention. 14 Diseases you almost forgot about (thanks to vaccines); 2022 Sep 15. Available from: https://www.cdc.gov/vaccines/parents/diseases/forgot-14-diseases.html.
- A. Beukenhorst, C.M. Koch, C. Hadjichrysanthou, et al. SARS-CoV-2 elicits non-Sterilizing immunity and evades vaccine-Induced immunity: implications for future vaccination strategies. Eur J Epidemiol. 2023;38(3):237. doi: 10.1007/s10654-023-00965-x
- Centers for Disease Control & Prevention. Coronavirus disease 2019 (COVID-19); 2020 Feb 11. Available from: https://www.cdc.gov/coronavirus/2019-ncov/science/data-review/vaccines.html.
- World Health Organization. Ten health issues WHO will tackle this year; 2019. Available from: https://www.who.int/news-room/spotlight/ten-threats-to-global-health-in-2019.
- B. Buonomo, R. Della Marca, M. Groppi. A behavioural modelling approach to assess the impact of COVID-19 vaccine hesitancy. J Theor Biol. 2022;534:Article ID 110973. doi: 10.1016/j.jtbi.2021.110973
- E. Dubé, C. Laberge, M. Guay, et al. Vaccine hesitancy. Hum Vaccin Immunother. 2013;9(8):1763–1773. doi: 10.4161/hv.24657
- World Health Organization. Vaccine hesitancy: a growing challenge for immunization programmes; 2015. Available from: https://www.who.int/news/item/18-08-2015-vaccine-hesitancy-a-growing-challenge-for-immunization-programmes.
- World Health Organization. Coronavirus disease (COVID-19): vaccines; n.d. Available from: https://www.who.int/news-room/questions-and-answers/item/coronavirus-disease-(covid-19)-vaccines.
- D. Olivera Mesa, P. Winskill, A.C. Ghani, et al. The societal cost of vaccine refusal: a modelling study using measles vaccination as a case study. Vaccine. 2023;41(28):4129–4137. doi: 10.1016/j.vaccine.2023.05.039
- C.H. Leung, M.E. Gibbs, P.E. Pare. The impact of vaccine hesitancy on epidemic spreading. In: American Control Conference (ACC); 2022, p. 3656–3661. doi: 10.23919/ACC53348.2022.9867327
- B. Oduro, A. Miloua, O.O. Apenteng, et al. COVID-19 vaccination hesitancy model: the impact of vaccine education on controlling the outbreak. Math Model Appl. 2021;6(4):81–91. doi: 10.11648/j.mma.20210604.11
- M. Angeli, G. Neofotistos, M. Mattheakis, et al. Modeling the effect of the vaccination campaign on the COVID-19 pandemic. Chaos Solit Fractals. 2022;154:Article ID 111621. doi: 10.1016/j.chaos.2021.111621
- G. Ledder. Incorporating mass vaccination into compartment models for infectious diseases. medRxiv; 2022. p. 2022.04.26.22274335. doi: 10.1101/2022.04.26.22274335
- A. de Miguel-Arribas, A. Aleta, Y. Moreno. Impact of vaccine hesitancy on secondary COVID-19 outbreaks in the US: an age-structured SIR model. BMC Infect Dis. 2022;22(1):511. doi: 10.1186/s12879-022-07486-0
- Centers for Disease Control & Prevention. COVID-19 and your health. Centers for disease control and prevention; 2020 Feb. Available from: https://www.cdc.gov/coronavirus/2019-ncov/your-health/about-covid-19/basics-covid-19.html.
- E.P. Esteban, L. Almodovar-Abreu. Assessing the impact of vaccination in a COVID-19 compartmental model. Inform Med Unlocked. 2021;27:Article ID 100795. doi: 10.1016/j.imu.2021.100795
- I. Ahmed, G.U. Modu, A. Yusuf, et al. A mathematical model of coronavirus disease (COVID-19) containing asymptomatic and symptomatic classes. Results Phys. 2021 Feb;21:Article ID 103776. doi: 10.1016/j.rinp.2020.103776
- National Institute of Health. Clinical spectrum of SARS-CoV-2 infection; n.d. Available from: https://www.covid19treatmentguidelines.nih.gov/overview/clinical-spectrum/.
- T. Šušteršič, A. Blagojević, D. Cvetković, et al. Epidemiological predictive modeling of COVID-19 infection: development, testing, and implementation on the population of the benelux union. Front Public Health 9 (2021), pp. 1–24. doi:10.3389/fpubh.2021.727274.
- Johns Hopkins Coronavirus Resource Center. Mortality analyses; n.d. Available from: https://coronavirus.jhu.edu/data/mortality.
- C. Kerr, R. Stuart, D. Mistry, et al. Covasim: an agent-Based model of COVID-19 dynamics and interventions. PLoS Comput Biol. 2021 Jul;17(7):e1009149. PLoS Journals. doi: 10.1371/journal.pcbi.1009149
- A. Wilder-Smith. What is the vaccine effect on reducing transmission in the context of the SARS-CoV-2 delta variant? Lancet Infect Dis. 2022 Feb;22(2):152–153. doi: 10.1016/S1473-3099(21)00690-3
- D. Swan, C. Bracis, H. Janes, et al. COVID-19 vaccines that reduce symptoms but do not block infection need higher coverage and faster rollout to achieve population impact. Sci Rep. 2021 Jul 1;11(1):Article ID 15531. doi: 10.1038/s41598-021-94719-y
- Z.S. Kifle, L.L. Obsu. Mathematical modeling for COVID-19 transmission dynamics: a case study in Ethiopia. Results Phys. 2022 Mar;34:Article ID 105191. doi: 10.1016/j.rinp.2022.105191
- J. Mondal, S. Khajanchi. Mathematical modeling and optimal intervention strategies of the COVID-19 outbreak. Nonlinear Dyn. 2022;109(1):177–202. doi: 10.1007/s11071-022-07235-7
- P. Samui, J. Mondal, S. Khajanchi. A mathematical model for COVID-19 transmission dynamics with a case study of India. Chaos Solitons Fractals. 2020 Nov;140:Article ID 110173. doi: 10.1016/j.chaos.2020.110173
- C.W. Castillo-Garsow, C. Castillo-Chavez, A tour of the basic reproductive number and the next generation of researchers, in An Introduction to Undergraduate Research in Computational and Mathematical Biology, H. Callender Highlander, A. Capaldi, and C. Diaz Eaton, ed., Foundations for Undergraduate Research in Mathematics, Birkhäuser, Cham, 2020. doi:10.1007/978-3-030-33645-5_2
- P. van den Driessche, J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002 Nov;180(1):29–48. doi: 10.1016/S0025-5564(02)00108-6
- O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol. 1990 Jun;28(4):365–382. doi: 10.1007/BF00178324
- A.B. Gumel. Causes of backward bifurcations in some epidemiological models. J Math Anal Appl. 2012 Nov;395(1):355–365. doi: 10.1016/j.jmaa.2012.04.077
- B. Song. Basic reinfection number and backward bifurcation. Math Biosci Eng. 2021;18(6):8064–8083. mbe-18-06-400. doi: 10.3934/mbe.2021400
- M. Martcheva. Methods for deriving necessary and sufficient conditions for backward bifurcation. J Biol Dyn. 2019;13(1):538–566. doi: 10.1080/17513758.2019.1647359
- C. Castillo-Chavez, B. Song. Dynamical models of tuberculosis and their applications. Math Biosci Eng. 2004;1(2):361–404. doi: 10.3934/mbe.2004.1.361