Delay in booster schedule as a control parameter in vaccination dynamics

The use of multiple vaccine doses has proven to be essential in providing high levels of protection against a number of vaccine-preventable diseases at the individual level. However, the effectiveness of vaccination at the population level depends on several key factors, including the dose-dependent protection efficacy of vaccine, coverage of primary and booster doses, and in particular, the timing of a booster dose. For vaccines that provide transient protection, the optimal scheduling of a booster dose remains an important component of immunization programs and could significantly affect the long-term disease dynamics. In this study, we developed a vaccination model as a system of delay differential equations to investigate the effect of booster schedule using a control parameter represented by a fixed time-delay. By exploring the stability analysis of the model based on its reproduction number, we show the disease persistence in scenarios where the booster dose is sub-optimally scheduled. The findings indicate that, depending on the protection efficacy of primary vaccine series and the coverage of booster vaccination, the time-delay in a booster schedule can be a determining factor in disease persistence or elimination. We present model results with simulations for a vaccine-preventable bacterial disease, Heamophilus influenzae serotype b, using parameter estimates from the previous literature. Our study highlights the importance of timelines for multiple-dose vaccination in order to enhance the population-wide benefits of herd immunity.


Introduction
Vaccination remains the most effective intervention measure in preventing many infectious diseases (Ehreth 2003). Conferring high levels of protection against a number of vaccine-preventable diseases requires more than one dose of vaccine that may be offered at different ages according to specific schedules set by vaccination programs (Jackson et al. 2012;Riolo et al. 2013;Riolo and Rohani 2015). For instance, vaccine schedules against Haemophilus influenzae serotype b (Hib) recommended for infants includes either 3 primary doses without a booster, or 2 to 3 primary doses plus a booster given at least 6 months after completing the primary series (World Health Organization et al. 2016). However, even in the presence of booster doses, resurgence and outbreaks of some vaccine-preventable diseases still occur, notwithstanding substantial levels of routine primary vaccine series (Jackson et al. 2012;Riolo et al. 2013;Riolo and Rohani 2015). Reduced effectiveness of vaccination has been explicated for the occurrence of such outbreaks due to factors associated with incomplete protection efficacy of primary vaccine series, inadequate coverage of booster doses, waning immunity over time, and the duration of vaccine-induced protection that may be significantly shorter than the average lifetime of the population (Alexander et al. 2006;Riolo and Rohani 2015).
While the importance of age-at-vaccination and booster doses has been documented, the optimal vaccine schedules remains unclear for several vaccine-preventable diseases and scheduling is mainly determined based on epidemiological context in individual settings (Low et al. 2013;Jackson et al. 2013;Riolo and Rohani 2015). The diversity of booster dose schedules observed in immunization programs worldwide could have a significant impact on disease elimination, since the scheduling may also affect the uptake rates of booster vaccination (Fitzwater et al. 2010). This poses a particular challenge for public health immunization programs in the context of deferral and subsequent refusal of booster doses that diminish the herd immunity (Omer et al. 2009;Dubé et al. 2013;Briere et al. 2014), and could lead to disease resurgence. Identification of the optimal booster schedule is therefore an important component of vaccination policies.
Despite the importance of dosing interval between primary and booster vaccination, a theoretical framework to investigate the impact of such interval and varying vaccination schedules on disease dynamics in the population is currently lacking. In this study, we aimed to establish this framework by developing a vaccination model, represented by a system of delay differential equations that describe the dynamics of disease transmission. Using this system, we evaluated the effect of delay in booster dose after primary vaccination on the long-term disease prevalence. We incorporated a number of key parameters into the model including the protection efficacy of primary vaccination, duration of vaccine-induced protection, and the coverages of primary and booster vaccination. We considered the delay as a control parameter, and analyzed the transient and steady-state behaviours of the system, in addition to determining the effect of time interval between primary and booster doses on disease elimination and persistence. We show, by means of simulations, that the threshold of disease control depends critically on the parameter of delay in booster dose for a given protection efficacy of primary vaccination.
To study the dynamics of our vaccination model, we first propose a basic framework without vaccination, derive the basic reproduction number (R 0 ), and prove a threshold result for disease elimination in terms of R 0 . We then construct the general model by incorporating primary and booster vaccination into the basic framework, and analyze its behaviour. Stability of the disease-free equilibrium is investigated, and represented in terms of the control reproduction number (R c ). When R c > 1, we show the uniform persistence, indicating that the disease elimination is infeasible. Finally, we use parameter values estimated for a bacterial disease, Heamophilus influenzae serotype b (Hib), and perform simulations to illustrate the model results by varying time interval between primary and booster doses, represented by the delay parameter.

The basic framework
To develop the basic framework, we divided a population of constant size N into several compartments to represent the epidemiological statuses of individuals as susceptible (S), infectious (I ), recovered and fully protected (R), and partially protected (W , and R w ). The distinction between the two classes W and R w is based on the consideration that partial protection following natural infection may last for a certain period of time (on average) before declining towards negligible levels. We assumed a fixed duration of full protection following recovery from infection. Once this period has elapsed, individuals will have only partial protection, and may become infected at a reduced rate compared with fully susceptible individuals. The duration of partial protection is divided into a period of fixed length (for those in the W class), followed by an exponentially distributed period of waning immunity (for those in the R w class) leading to full susceptibility. Due to the fixed periods of full and partial protection, the dynamics of disease transmission can be expressed by the following set of equations: where β is the baseline transmission rate; γ is the recovery rate of infectious individuals, η is the reduction of susceptibility to infection due to partial protection, μ is the natural death rate (assumed to be the same as the birth rate), τ r represents the fixed duration of full protection, τ w represents the fixed duration of partial protection; and θ is the rate of loss of immunity in partially protected individuals in the R w class. A schematic diagram of the transitions between different classes of individuals in this basic framework is represented in Fig. 1.
Let a, referred to as 'age', be the time elapsed since individuals enter each class. Thus, r (t, a) and w(t, a) represent the density, with respect to age a at time t, of recovered individuals having full protection, and those having partial protection during the fixed period before moving to the exponentially distributed period, respectively. Using the above notation, we have the following equations: Solving along the characteristics gives: where, for simplicity we used the notationβ = β/N . Note that system (1) includes integral and differential equations. By differentiating the R and W equations in (1) with respect to t, we obtain the following model, which together with (3), constitutes a closed system of delay differential equations: One can easily check that the population Let τ = τ r + τ w , and denote by C the Banach space C([−τ, 0], R 5 ) of continuous functions mapping the interval [−τ, 0] into R 5 equipped with the norm: where φ ∈ C and | · | is a norm in R 5 . For a continuous function u : We choose the initial conditions for system (4) from the set Ω ⊆ C, defined by: The following result shows that system (4) is well-posed in Ω, and the solution semiflow admits a global attractor on Ω.

Proof
We start with the last assertion. From the third equation of (4), we have: Integrating both sides yields: Similarly, by integrating and using the initial condition (4) gives: For a given φ ∈ C, we define Thus, system (4) can be written as u (t) = G(u t ). We show that G(φ) is Lipschitzian in φ within each compact set in C, that is for all M > 0 there is a K > 0 such that for all φ, ψ ∈ C with φ ≤ M and ψ ≤ M, the inequality |G(φ) − G(ψ)| ≤ K φ − ψ holds. We note that there are terms of linear, quadratic and exponential types in G. Quadratic terms are all Lipschitzian, and one can see that: The Lipschitzian property for the most involved term can also be seen from: where we used the mean value theorem e x − e y = e ξ (x − y). Hence, there is a unique solution of the system through (0, φ) on its maximal interval of existence [0, σ φ ). We note from (4) that any solution satisfies: and hence if I (0) ≥ 0, then I (t) ≥ 0 for all t ∈ (0, σ φ ). From (6) and (7), we find that R and W are non-negative. Then the non-negativity of S and R w follows from the inequalities The non-negativity and relations (6) and (7) ensure that Ω is forward invariant. Since the total population is constant, it follows that S(t) and I (t) are bounded by N and the solutions exist globally. Therefore, the solution semiflow Φ(t) = u t (·) : Ω → Ω is point dissipative. By Theorem 3.6.1 in Hale (1977), Φ(t) is compact for any t > τ. Thus, from Theorem 3.4.8 in Hale (1988), it follows that Φ(t) has a compact global attractor in Ω.

Basic reproduction number
The basic reproduction number (R 0 ) is the average number of new infected individuals generated by a single infected individual introduced into an entirely susceptible population, during the course of infection. According to the theory of epidemics, we expect that the disease will vanish if R 0 < 1, while it will persist in the population if R 0 > 1. New infections occur only in the S class with the rateβ I . In a fully susceptible population, S/N ≈ 1 and the average length of infection is (μ + γ ) −1 , and therefore we define the basic reproduction number as R 0 = β/(γ + μ). We proceed by presenting the threshold dynamics for system (4), which determines whether the disease dies out or persists.

Threshold dynamics
It is clear that system (4) has a unique disease-free equilibrium E 0 = (N , 0, 0, 0, 0). We first show that E 0 is globally asymptotically stable when R 0 < 1. Then we establish the uniform persistence of the disease when R 0 > 1 using techniques of persistence theory (Smith and Thieme 2011).
Theorem 2.2 If R 0 < 1, then the disease-free equilibrium E 0 of system (4) is globally asymptotically stable in Ω.
Proof Linearizing system (4) at E 0 , we obtain the following system: where and all other components are zero. The characteristic equation of system (9) has the form: which implies that E 0 is asymptotically stable. Now we prove that E 0 is globally attractive in Ω. We consider the following inequality: (6) and (7), we find that R(t) and W (t) are also converging to zero. For any > 0, and a sufficiently large t > 0, we have w(t, τ w ) < , and the following inequality holds: This means that lim sup t→∞ R w (t) ≤ μ + θ , and therefore R w → 0 as t → ∞. Since the total population is constant, we obtain that S(t) → N as t → ∞. Thus, lim t→∞ u(t, φ) = (N , 0, 0, 0, 0), and we conclude that E 0 is globally asymptotically stable.
We now prove the persistence of disease for R 0 > 1. Considering the semiflow Φ on Ω, we define the persistence function by: where Ω 0 is called the extinction space corresponding to P (that is the collection of states where the disease is not present). From the relation (8), it follows that the sets Ω 0 and Ω + are forward invariant under the semiflow Φ. We now introduce some terminology from persistence theory (Smith and Thieme 2011, Chapters 3.1 and 8.3).

Definition 2.3
Let X be a nonempty set and P : X → R + .

A semiflow
For such a solution, which contradicts the convergence of I to zero. Hence, E 0 is weakly P-repelling.
Theorems 2.2 and 2.4 indicate that the disease dynamics are completely determined by the basic reproduction number R 0 . In the following, we extend our model by introducing the primary and booster vaccination, and analyze the persistence dynamics of the resulting system.

The general vaccination model
The vaccination model includes additional classes of individuals who are vaccinated with primary series; partially protected following primary vaccination; and fully protected following booster vaccination (See Table 1). A schematic diagram for timelines of primary and booster vaccination with delays is represented in Fig. 2. We assume that a fraction p of newborns will receive primary vaccine series in the first τ s period of their life. The remaining fraction of newborns will be recruited to the susceptible class and can therefore become infected through contacts with infectious individuals. The primary vaccination is assumed to provide partial protection for a certain period of time during which infection can occur with a lower rate compared with fully susceptible individuals. We also assume that partial protection induced by primary vaccine gradually wanes over time, and individuals who forgo booster vaccination will eventually become susceptible. Similar to the basic framework, we consider a fixed duration of partial protection after primary vaccination, followed by an exponentially distributed period of waning immunity leading to full susceptibility. Those who have received primary vaccination may also receive booster dose. We assume that, similar to recovery from infection, booster vaccination provides a fixed duration of full protection, followed by fixed and exponentially distributed durations of partial protection.  In order to mathematically express the model, we let a represent the age since the individuals enter each class, and define v s (t, a), v p (t, a), and v d (t, a) to represent the density, with respect to age a at current time t, of primary vaccinated individuals, those who are partially protected as a result of primary vaccination and eligible to receive booster dose, and those who are partially protected by primary vaccination and will not receive booster dose, respectively. Let τ p (0 ≤ τ p ≤ τ w ) represent the delay in receiving booster vaccination within the fixed period of partial protection following primary vaccination. For those who will not receive the booster, we denote by τ d the remaining period of time within the fixed duration of partial protection, i.e., τ p + τ d = τ w . Description of all model parameters are provided in Table 2. With the above notation, the model can be expressed by the following system of integro-differential equations: with population densities expressed by: In this model, we considered a single class V w for partially protected individuals with exponentially distributed duration of protection, regardless of whether the immunity was conferred by vaccination or natural infection. Thus, the R w class from the basic framework is included in the V w class in the vaccination model. Solving along the characteristics gives: Differentiating the integral equations V s , V p , V d , V b , W and R in (11) and substituting the density functions, we obtain the following system of differential equations with delays: where

The total population S(t)+V s (t)+V p (t)+V d (t)+V w (t)+V b (t)+W (t)+I (t)+R(t) =
N . Note that the equations of V b and R are decoupled from the rest of the model. In the following, we show that (19) is well-posed, and further define the diseasefree equilibrium and the control reproduction number. Let τ c = max{τ s + τ p + τ b + τ w , τ r + τ w }. We choose the initial conditions for system (19) from the set X , which is defined by Theorem 3.1 For any φ ∈ X , system (19) has a unique non-negative solution U (t, φ) satisfying U 0 = φ and U t ∈ X , and the solution semiflow Φ(t) = U t (·) : X → X has a compact global attractor. Moreover, the solutions of system (19) with initial conditions in X satisfy the integro-differential equations system (11).
Proof The proof is similar to Theorem 2.1.

Reproduction number
Recall that in the absence of vaccination, system (19) reduces to (4) and the basic reproduction number is given by R 0 = β/(μ + γ ). Letting I (t) ≡ 0, we obtain the unique disease-free equilibrium of system (19), Linearizing system (19) at E * 0 , we obtain the following equation for the infection class: We now introduce the reproduction number following the idea in (Xu and Zhao 2012).
Denote by x 0 the number of infectious individuals at time t = 0, and x 1 (t) be the remaining population at time t. Thus, Thus, from (21), the total number of newly infected cases is Therefore, we define the reproduction number in the presence of vaccination by which can be interpreted as the total number of new infections generated by a single infectious individual in all non-infection classes during the average infectious period 1/(μ + γ ).

Local stability
Here, we show disease elimination for sufficiently small I (corresponding to solutions in a small neighbourhood of E * 0 ) if R c < 1. Notice that V b and R in system (19) are independent of other state variables. In the following, we consider (19) with the additional equations: Linearizing (19) at E * 0 , we have where

t), I (t), A(t), B(t), C(t) T , and
This is a block triangular matrix with the zero block O 4×6 , and All matrices D j , j = 1, . . . 9, that appear in (23) can also be derived from (19) and (22).

Theorem 3.2
If R c < 1, then the disease-free equilibrium E * 0 of system (19) is locally asymptotically stable.

Remark 3.3
We have not been able to establish the global stability of E * 0 when R c < 1. There are some vaccination models that exhibit the phenomenon of backward bifurcation, where a stable endemic equilibrium co-exists with the stable disease-free equilibrium Gumel (2002). However, based on the simulation results presented in the next section, we conjecture that Theorem 3.2 holds for the entire domain of system (19) and E * 0 is globally stable.

Uniform persistence
In this section, we prove the disease persistence when R c > 1. Consider the semiflow Φ(t) in X , defined by the unique global solutions. We define the persistence function: where X 0 is the extinction space corresponding to P (i.e., X 0 is the collection of states without disease presence). From Theorem 3.1, it follows that X 0 and X + are forward invariant under the semiflow Φ.
Theorem 3.4 If R c > 1, then the semiflow Φ is uniformly P-persistent, i.e. there is a δ > 0 such that for any solution lim inf t→∞ I (t) ≥ δ.
Proof First we show that E * 0 is weakly P-repelling. Suppose that there exists ψ 0 ∈ X such that P(ψ 0 ) > 0 with For such a solution, I (0) > 0 and lim t→∞ I (t) = 0. For sufficiently small > 0, we have R c − 6 γ + μ > 1. Hence, for sufficiently large t, we get Φ(t, ψ 0 )− E * 0 < , and which contradicts the convergence of I to zero. Thus, E * 0 is weakly P-repelling. We also note that whenever I (t) ≡ 0, from the R equation in (11) and Theorem 3.1, it follows that 0 }, and one can see that Φ is uniformly weakly P-persistent (Smith and Thieme 2011, Theorem 8.17). Since Φ has a compact global attractor on X , the application of Theorem 4.5 in Smith and Thieme (2011) guarantees that Φ is uniformly P-persistent.

Remark 3.5
It is expected that an endemic equilibrium exists when the disease is uniformly persistent. Due to the exponential terms in the model, it is not possible to derive an explicit formula for the components of an endemic equilibrium, and we could not prove the existence or uniqueness of such equilibrium using established methods (such as fixed point arguments). However, our numerical experiments suggest that there is a unique endemic equilibrium that emerges as R c increases and passes the threshold of one. Regarding its stability, it is known that SIRS models with delay can exhibit periodic oscillations Hethcote et al. (1981), and their endemic equilibria can either be stable or unstable. For our model, a linear stability analysis seems very difficult to conduct due to the various delay terms. In numerical simulations, however, we can readily find a combination of parameter values for which the vaccination model (19) exhibits sustained oscillations in the disease prevalence. This typically occurs for a small vaccination coverage, and increasing this coverage first stabilizes the endemic equilibrium, and then can lead to disease elimination when it is sufficiently high to bring R c less than one.

Simulation results
To illustrate the effect of booster schedule on the dynamics of disease spread in the population, we simulated the model while varying the protection efficacy of primary vaccination and the coverage of booster vaccination. For the simulation results presented here, we used parameter values estimated for Haemophilus influenzae serotype b (Hib) in the published literature. Primary vaccination for Hib in most routine infant immunization programs includes 2 to 3 doses of vaccine offered between 2 to 6 months of age (World Health Organization et al. 2016), and we therefore assumed τ s = 6 months for completion of primary series. Primary vaccination is estimated to provide partial protection for a fixed duration of τ w = 4 years, followed by an exponentially distributed time period with the average of 1/θ = 6 years (Konini and Moghadas 2015;Jackson et al. 2012). Booster vaccination provides full protection for a fixed period of τ b = 6 years (Konini et al. 2016;Leino et al. 2000). We assumed that, after the period of full protection has elapsed, partial protection follows the same timelines as primary vaccination. Similar to booster vaccination, we assumed that recovery from infection provides full protection for a fixed period of τ r = 4 years (Konini and Moghadas 2015). Infection in the form of carriage (i.e., asymptomatic without showing clinical symptoms) contributes more significantly to the incidence of Hib compared to symptomatic disease, and has a prolonged infectious period from several days to several weeks (Leino et al. 2000;Jackson et al. 2012). We therefore assumed an average infectious period of 1/γ = 50 days.
For the purpose of simulations, we used a population of size N = 100, 000 with an average lifetime of 1/μ = 70 years. The transmission parameter β was calculated based on a given basic reproduction number, while fixing other parameters of the model. We assumed R 0 = 1.2 in the range 1.1-1.4 estimated in studies of Hib (Farrington et al. 2001) and other pathogens that cause bacterial meningitis, such as Neisseria meningitidis serotype C (Stephens 2011). We fixed the coverage of primary vaccine at p = 0.9, and varied the protection efficacy of primary vaccination, reflected in the reduction of susceptibility to infection. We ran the simulations while changing the time for booster vaccination within the fixed period of partial protection following primary vaccination, that is, 0 < τ p ≤ τ w = 4 years.
We also simulated R c as a function of two model parameters, namely the protection efficacy of primary vaccination (η), and the time for booster vaccination following primary series (τ p ). For these simulations, we considered the coverage of booster vaccination as a function of τ p in three different scenarios: (i) Fixed coverage of booster vaccination: ρ = 0.95 (Fig. 3, solid line). (ii) Exponentially declining coverage of booster vaccination: This coverage reduces as the time delay τ p in booster vaccination increases (Fig. 3, dashed line). (iii) Inverted logistic declining coverage of booster vaccination: This coverage reduces with time delay τ p in booster vaccination in a functional form similar to van Genuchten-Gupta model (Fig. 3, dotted line). Figure 4 shows the variation in R c corresponding to the scenarios of booster coverage. For a fixed coverage of booster vaccination (ρ = 0.95), Fig. 4a shows that when the protection efficacy of primary vaccination is sufficiently high (approximately above 70%), R c decreases with increasing delay in booster schedule following primary vaccination, and the disease can be eliminated if R c < 1 (in the region to the left side of the white line). When primary vaccination provides a protection efficacy that is nearly as good as that conferred by the booster dose (i.e., η > 0.9), the disease can be eliminated regardless of the time for booster schedule. However, for a moderate to low protection efficacy of primary vaccination, the delay in booster vaccination has little or no effect in reducing R c and the disease persists in the population.
When the coverage of booster vaccination declines exponentially, we observed lower R c for early booster schedule, regardless of the protection efficacy of primary vaccination (Fig. 4b). In our simulations, disease elimination can occur with a protection efficacy above 90%, but requires booster vaccination within 12 months following the primary vaccination (i.e., the region below the white line in Fig. 4b). In contrast to the scenario for a fixed coverage of booster vaccination, these simulations suggest that an early booster dose may be essential in curtailing disease spread if the coverage of booster vaccine is expected to decline (exponentially) over time. This scenario may correspond to vaccine refusal in the contexts of booster deferral (Centers for Disease Control and Prevention et al. 2009).
Further simulations indicate that functional form of the decline in booster coverage can also play an important role in determining the optimal timing of booster vaccination. With the inverted logistic functional form of ρ represented by the dotted curve in Fig. 3, we observed that the protection efficacy of primary vaccination can influence the magnitude of reduction in R c with the time delay in booster schedule (Fig. 4c). For a moderate to low protection efficacy of primary vaccination, early booster (similar to the scenario of exponential decline in ρ) leads to the maximum reduction in R c . However, as the protection efficacy of primary vaccination increases (approximately above 50% in these simulations), the maximum reduction of R c corresponds to an intermediate time-interval for booster vaccination. Figure 4c indicates that an optimal timing of booster schedule may lead to disease elimination, while the disease can persist in the population if the booster dose is offered too early or too late following primary vaccination.
To further illustrate our findings in terms of disease prevalence, we simulated the model for the scenarios of booster coverage, while fixing the protection efficacy of the primary vaccination. Figure 5a shows that for fixed ρ = 0.95 and η = 0.8, the disease will be eliminated if the booster dose is offered 30 months after the primary vaccination. However, an earlier schedule of a booster dose 6 months after the primary vaccination leads to the disease persistence in the population. This situation is reversed for the scenario of booster coverage that declines exponentially with time delay in booster vaccination. Figure 5b shows the disease elimination and persistence for η = 0.95, with the booster dose schedules of 2 and 24 months after the primary vaccination, respectively. When the booster coverage declines in a functional form similar to the inverted logistic function, Fig. 5c shows the disease persistence for early and late booster schedules of 1 and 24 months after the primary vaccination with η = 0.88. However, for an intermediate delay of 9 months in booster vaccination, the disease is eliminated over time.  Fig. 5b). Colour curves represent the disease prevalence over time with changing population size. Parameter values are the same as Fig. 5b with τ p = 24 months, while the birth rate changes by Δ Remark 4.1 For our theoretical results and simulations presented in Fig. 5, we assumed that the population size (N ) is constant. To illustrate the effect of a changing population size on the disease dynamics, we considered the parameter setting of Fig. 5b with τ p = 24 months, and modified the birth rate by some constant value Δ. Hence, Δ > 0 corresponds to a growing population size, while Δ < 0 represents a declining population. The results are illustrated in Fig. 6, showing the change in disease prevalence for different values of Δ as the population size changes.

Discussion
In this study, we investigated the role of booster schedule on the long-term disease dynamics. We developed a system of delay differential equations to include several key parameters describing the protection efficacy of primary vaccine series, durations of partial and full protection following vaccination, and coverage of primary and booster doses. In addition to investigating its dynamics, we simulated the model with a delay in booster dose after completing primary series using parameter values estimated for Hib. Simulation results indicate that, for a given protection efficacy of primary vaccination, the reduction of disease transmissibility, reflected in the reproduction number (R c ), depends critically on the timing of a booster dose. However, the coverage of booster vaccination remains a key parameter influencing the optimal timing of a booster dose. When the uptake of a booster dose is expected to remain high, a delay in booster vaccination (within the expected duration of protection induced by primary vaccine series) may be beneficial in reducing R c , and could lead to disease elimination for a sufficiently high protection efficacy of primary vaccination. This is particularly important if the booster dose provides only a relatively short period of full protection compared with the average lifetime. However, vaccination programs may contend with the possible drop-out and decrease in the coverage of booster vaccination, whether due to acquiring infection after receiving the primary series, or simply due to individuals voluntarily forgoing (e.g., refusal of) the booster dose (Omer et al. 2009;Dubé et al. 2013;Briere et al. 2014). In this case, our simulations illustrate that if the coverage of booster vaccination decreases over time, then the timing of a booster dose can be essential to achieve the greatest reduction of disease prevalence over time.
This study has important implications for public health vaccination policies. First, for routine infant immunization programs with high primary and booster coverages, deferral of a booster dose within the average duration of protection induced by the primary series may be beneficial. However, the protection efficacy of the primary vaccine series remains an important parameter in determining the optimal dosing interval between primary and booster vaccination (Charania and Moghadas 2016). Second, in the absence of efforts to achieve an optimal schedule, having a booster program does not necessarily guarantee the elimination of disease, even though the incidence may be reduced as has been observed for Hib. Given the high protection efficacy of primary vaccine series (> 85%) against Hib disease (Jackson et al. 2013), and weak evidence of additional protection from booster within one year following complete primary vaccine series, our results suggest that immunization programs should consider a longer time interval between primary and booster doses. Furthermore, the sensitivity of longterm disease outcomes to the booster schedule underscores the importance of targeted efforts towards improving uptake rates of both primary series and booster vaccination. It is also important to note that our results herein apply to vaccines that confer only a temporary protection. Conspicuously, for a vaccine that provides a long-term full protection comparable to the average life-time, the best outcomes are achieved with the shortest time interval between the primary and booster vaccination.
Our study has several limitations that merit further investigation. Our model is based on the assumption of homogeneous mixing in the population dynamics of disease spread. It is well documented that heterogeneities and contact patterns can influence vaccination dynamics at both the individual and population levels (Metcalf et al. 2015). We assumed a uniform protection efficacy of primary and booster vaccination without considering immunological characteristics of individuals that affect the within-host immune dynamics. Our simulation results are based on the assumption that an anti-Hib polysaccharide conjugate vaccine provides stronger immune protection due to effects of carrier protein on stimulation and proliferation of immune responses. This is an important consideration in the development of conjugate vaccines for T-cell independent pathogens (such as Hib) in order to enhance immunogenicity in infants and young children (Goldblatt 2000). We therefore assumed a shorter period of full protection following recovery from natural infection compared to booster vaccination. However, in older individuals with competent immune system, natural infection can also lead to the development of adaptive immune memory and therefore induce strong protection effects with timelines similar to those conferred by conjugate vaccines (Goldblatt 2000). These considerations can be included in advanced computational frameworks, such as agent-based modelling (Laskowski and Moghadas 2014;Shoukat et al. 2018), in order to evaluate the effect of individual level characteristics on the population dynamics of disease spread and control in the presence of vaccination. Despite these limitations, our study provides a theoretical foundation for future studies involving more detailed computational and quantitative models to help improve vaccination programs and booster schedules against vaccine-preventable diseases that require multiple vaccine doses.