The honeymoon period after mass vaccination

1.   Introduction

Vaccines are biological preparations that induce a host's immunity against a disease without causing the actual illness. Their importance has been highlighted recently by the demand for the rapid development of a vaccine for COVID-19. Models help us determine how vaccination programs can change disease dynamics, depending on many factors, including the effectiveness of the vaccine and how it is distributed. Here we consider how different types of individual-level protection afforded by a vaccine can result in different population-level dynamics of a disease, with a focus on "honeymoon periods".

McLean and Anderson ^[2] introduced the term honeymoon period to denote the "long period of low incidence that immediately follows the initiation of a mass vaccination program". There are many papers that note the occurrence of honeymoon periods in both real data and epidemiological model trajectories. For example, Hefferman and Keeling ^[3] showed that a waning vaccine model with moderate waning times (between 40 to 80 years) and high vaccine coverage (greater than 70%) can lead to considerable oscillations in the infected class, with a long low incidence period preceding possibly large-scale epidemics. Mossong and Muller ^[4] examined the time to re-emergence of measles in different situations using age-structured models. In Figure 1 we present an example of a possible honeymoon period in pertussis data from Massachusetts.

Most research on mathematical modeling of infectious diseases focuses on the equilibrium dynamics of epidemiological systems. However, transient dynamics can also have some very unique and informative characteristics, complementing the steady state analysis of systems. The honeymoon period is a transient feature that can be an important consideration for public health agencies.

In this paper, we conduct a study of honeymoon periods arising from a simple model of an imperfect vaccine examined in ^[1]. This model allows for three different modes of vaccine failure: "leakiness", "all-or-nothingness" and "waning." We reviewed existing linearization methods for measuring transient dynamics in ecology, such as reactivity and amplification envelope ^[5], and conclude that these are not enough to determine the characteristics of honeymoon periods. We present proof of global asymptotic stability of the endemic equilibrium when the disease-free equilibrium is unstable and give a technical definition of the honeymoon period for these models. We also provide some bounds for the lengths of these periods and the size of the drop in the "effective susceptible" fraction of the population.

This paper is structured as follows: In Section 2 we review the imperfect vaccine model analyzed in ^[1] and review measures of transient dynamics defined by ^[5] that are based on linearization about the equilibrium. In Section 3 we demonstrate that the measures based on linearizaton are not helpful in distinguishing between models that display honeymoon periods and those that do not. We also define the effective susceptible class and present a proof of global stability of the endemic equilibrium, in order to give a technical definition of the honeymoon period. In Section 4 we prove some basic properties of the solution during the honeymoon period. Finally, in Section 5 we give a summary of the results and discuss future directions.

2.   Materials and method

2.1.   VSIR model

The imperfect vaccine model in Magpantay et al. ^[1] is an extension of the vaccinated-susceptible-infected-recovered (VSIR) model with unstructured population described by a system of ordinary differential equations (ODEs) similar to that presented in ^[6]. The three different modes of vaccine failure considered in the model are:

● Leakiness ("failure in degree")

The probability of getting infected upon exposure is decreased but not eliminated by vaccination.

● All-or-nothingness or primary vaccine failure ("failure in take")

Lifelong perfect protection is only provided to a fraction of the vaccinated population while the remaining fraction remains as susceptible as if they were never vaccinated.

● Waning ("failure in duration")

Perfect protection is provided to all vaccinees for an exponentially distributed period of time after which each vaccinee becomes as a susceptible as if they were never vaccinated.

It was shown that these different moves of protection have different epidemiological consequences ^[1]. Here we present a short description of the model: Let $S(t)$, $I(t)$ and $R(t)$ be the fractions of susceptible, infected and recovered individuals respectively in the population at time $t$. Let $V(t)$ be the proportion of individuals that has immunity from vaccination that has not waned yet. We assume a constant fraction $p$ of newborns are vaccinated soon after birth but the vaccination fails in a fraction $\varepsilon_A$ of the vaccinees. Thus $(1-p)+\varepsilon_A p = 1-(1-\varepsilon_A)p$ of newborns remain susceptible while $(1-\varepsilon_A)p$ go to the vaccinated compartment. We assume that individuals in the vaccinated compartment can still become infected, although at a lower rate than susceptible individuals. This "leakiness" of the vaccine, denoted by $\varepsilon_L$, reflects the reduction in probability of infection upon exposure of individuals in the $V$ class relative to those in the $S$ class. Finally, we also assume that vaccine protection wanes at a rate of $\alpha$. As shown in Table 1, we can re-parametrize from the waning rate $\alpha$ to a waning probability $\varepsilon_W$ by setting $\varepsilon_W = \frac{\alpha }{\alpha+\mu}$. The full set of model equations is given by (2.1):

with the requirement (2.2) on the initial conditions at time $t = 0$,

Here $\mu \geq 0$, $\gamma > 0$, $\beta > 0$, $p \in [0, 1]$, and $(\varepsilon_L, \varepsilon_A, \varepsilon_W) \in [0, 1)^3$. The model is illustrated in Figure 2 and the parameters are described in Table 1. If $(\varepsilon_L, \varepsilon_A, \varepsilon_W) = (0, 0, 0)$ then there is no way for a vaccinated individual to become infected so this special case is when we have a vaccine that induces perfect immunity.

2.2.   Disease-free equilibrium

To find the disease-free equilibrium, we set all the equations in (2.1) equal to zero with $I = 0$. Solving this yields the disease-free equilibrium given by

As in ^[1], we can show that the basic reproduction number with vaccination is given by

where $R_0 = \frac{\beta}{\gamma+\mu}$ is the basic reproduction number with no vaccination ($p = 0$). The vaccine impact $\varphi$ is defined to be

The equation for the vaccine impact above shows that it can be thought of as the relative change in the reproduction number per unit of vaccination. This can be used to derive the critical vaccination coverage $p_c$ (the minimum fraction that needs to be vaccinated in order to drive the disease to extinction) given by

if this quantity is less than one.

2.3.   Endemic equilibrium

To obtain the endemic equilibrium of the system we again set the derivatives in (2.1) to zero but assume this time that the infected class is not zero. This yields

where

and $a = \varepsilon_L-(1-\varepsilon_W)(R_0-1)+1$ and $R_0 = \frac{\beta}{\gamma+\mu}$. In ^[1], it was shown that this equilibrium is in the positive cone and is locally asymptotically stable when $R_p > 1$. We prove in the next section that the endemic equilibrium is globally stable when $R_p > 1$.

It was also shown in ^[1] that for a fixed vaccine impact parameter $\varphi$ such that $R_p > 1$, the vaccine failure parameter set $(\varepsilon_L, \varepsilon_A, \varepsilon_W)$ that leads to the highest value of the endemic steady state value of the infected class is $(1- \varphi, 0, 0)$; that is the purely leaky model. Additionally, the endemic equilibrium of $(\varepsilon_L, \varepsilon_A, \varepsilon_W) = (0, 1- \varphi, 0)$ is the same as that of $(\varepsilon_L, \varepsilon_A, \varepsilon_W) = (0, 0, 1- \varphi)$. Thus, we cannot distinguish between these purely all-or-nothing and purely waning models with the same vaccine impact by just looking at the equilibrium values of the components. However, these models can be differentiated if we look at their transient dynamics, in particular what if we initialize the models at pre-vaccine era equilibrium (with $p = 0$) then we increase the value of $p$ (refer to Figure 3).

2.4.   Measures of transient dynamics

Suppose that we have an ordinary differential equation given by

with an asymptotically stable equilibrium point at $x = \tilde{x}$. Let $z = x-\tilde{x}$ be the deviation of the system from equilibrium and $z_0 = z(0) = x(0)-\tilde{x}$ be the deviation at time zero. Let $J = \left.\frac{df}{dx}\right|_{x = \tilde{x}}$ be the Jacobian of the system (2.7) evaluated at $\tilde{x}$. The linearization of (2.7) about the equilibrium $x = \tilde{x}$ is given by,

Since all solutions of (2.1) are invariant in the positive cone with $V(t) + S(t) + I(t) + R(t) = 1$ for $t\geq0$ then $R(t)$ is determined by the values of the three other states. Thus, when we look at the Jacobian of the system, we can just consider the matrix of partial derivatives of the first three equations in (2.1) with respect to the first three states. The general form of the Jacobian at $\tilde{x}$ is,

To evaluate the Jacobian at the disease free equilibrium we set $\tilde x = x^0 = (V^0, S^0, I^0, R^0)$ and get,

To evaluate it at the endemic equilibrium we set $\tilde x = x^* = (V^*, S^*, I^*, R^*)$ and obtain,

We note that the element in the third row and the third column simplifies to zero after substituting in the values of $V^*$ and $S^*$ from (2.3).

Next, we review quantities that measure how sensitive a system is to deviations from equilibrium. Later on we will evaluate how these quantities depend on the parameter values of the imperfect vaccine model. We begin by discussing the reactivity and amplification envelope, two quantities defined by Neubert and Caswell ^[5] which have been used to measure transient dynamics when approaching a bifurcation point ^[7].

Definition 1 (Reactivity, from ^[5]).

The reactivity of an asymptotically stable equilibrium is the maximum, over all perturbations, of the rate at which the trajectory departs from the equilibrium in the linearized system. It measures the maximum instantaneous amplification of perturbations of that equilibrium. Using $|z| = \sqrt{z^T z}$ in (2.10), we derive

Let $H(J) = (J+J^T)/2$. This is the Hermitian part of $J$ and therefore its eigenvalues are all real-valued. The eigenvalues of $H(J)$, which do not necessarily have the same sign as the eigenvalues of $J$, determine the instantaneous behavior of the linearized system (2.8). Let $\lambda_H$ be the largest eigenvalue of $H(J)$. Then we derive the following expression for reactivity in terms of $J$,

Definition 2 (Amplification envelope, from ^[5]).

The amplification envelope is the maximum possible deviation from the equilibrium at any time following a perturbation. We note that the solutions $z(t)$ of the linearized system (2.8) are given by $z(t) = e^{J t} z_0$, where $e^{J t} = \sum_{n = 0}^{\infty} {\frac{(J t)^n}{n!}}$. Using this in (2.13) we derive,

where $\sigma(e^{J t})$ denotes the square root of the largest eigenvalue of $(e^{J t})^T(e^{J t})$. Note that the amplification envelope changes with time.

2.5.   Numerical solutions of the model equations

We used numerical methods to determine the size and extent of the honeymoon periods we obtain for the models evaluated at a range of different parameter values. MATLAB's $\mathtt{ode45}$ was used for the numerical integration of the systems of differential equations.

3.   Results

In Figure 3, the dynamics of the infected class for the leaky, all-or-nothing (AON) and waning models with the same vaccine impact value are shown after the start of mass vaccination. We see clearly different transient dynamics between the three. All three initially start off at the pre-vaccine ($p = 0$) endemic steady state and then vaccination coverage is set to $p = 0.9$. The waning model displays a dip below the new equilibrium that lasts for about 50 years and then the solution displays large oscillations as it approaches its new equilibrium. In contrast, the all-or-nothing model undergoes smaller oscillations about its new equilibrium. As expected from the previous chapter, both the all-or-nothing and waning models approach the same new endemic equilibrium value despite very different transient dynamics before reaching it. The leaky model approaches a different new endemic equilibrium value which has a higher infected population at its new endemic equilibrium. It also approaches this equilibrium in a different manner, after dropping below it at first and then rebounding, taking many years of small oscillations before reaching values close to the new equilibrium value.

We see from Figure 3 that the initial reduction in incidence after the start of mass vaccination is often much larger than the long-term reduction in incidence (once the system reaches its new equilibrium). We also observe that the waning model exhibits a long, pronounced honeymoon period.

3.1.   Plots of reactivity and amplification envelope

We derived the reactivity about both the endemic and disease-free equilibrium point of the three different vaccine models with the same vaccine impact. These were plotted in Figure 4(a)-(b) against the vaccine failure parameters. We see that there is only a small difference between the values of reactivity for the all-or-nothing and waning models despite Figure 3 showing a large difference in their transient dynamics. This indicates that reactivity is not a good measure for determining whether or not a model will display a honeymoon period. Interestingly, it also illustrates that the maximum possible instantaneous rate of amplification in the leaky model is less than the two other models.

We derived the amplification envelope of the three different vaccine models and plotted the results in Figure 4(c)-(d). Here we can observe again that there is not much of a difference between the results for the all-or-nothing and waning models, and that of the leaky model is also only slightly different from the other two. We find that the amplification envelope is also not a good indicator of the existence of significant honeymoon periods.

When observing honeymoon periods, we set the initial conditions of the models to be the pre-vaccine era equilibrium and consider how it approaches the new vaccine era equilibrium. Usually these two equilibria are significantly different from each other. This is why the reactivity and amplification envelope, which are based on linearization about equilibria and are more appropriate for measuring the behavior of small perturbations from such equilibria, may not be good indicators of which models will display substantial honeymoon periods.

3.2.   The effective susceptible class

It turns out that to analyze the honeymoon period, it will be important to estimate how the susceptibility of the population changes over time. Since individuals in the vaccinated compartment are still partially susceptible to the disease if $\varepsilon_L > 0$, we can think of the "effective" susceptible number of individuals in a population to include both the susceptible individuals and a fraction $\varepsilon_L$ of the vaccinated compartment.

Definition 3. We define the "effective susceptible class" $S_E(t)$ by

This effective susceptible class is related to the effective reproduction number $R_\text{eff}(t)$ of system (2.1) via the relationship $R_\text{eff}(t) = R_0 (S(t) + \varepsilon_L V(t)) = R_0 S_E(t)$. In Figure 3, we see that the infected class does not change noticeably during the honeymoon period. This is usually the only compartment that is observed and recorded. However other important changes are occurring in the other components of the system during the honeymoon period, causing the resurgence of the disease later on. In Figure 5, we present the trajectory of the effective susceptible class that corresponds to the trajectory shown in Figure 3. We observe that the leaky and all-or-nothing models show very similar behavior whereas the waning model displays very different dynamics. We also note from the trajectories that the equilibrium value of the effective susceptible class is independent of the type of vaccine failure or vaccine coverage. This is proven in Proposition 4.

Proposition 4. Suppose that $R_0 > 1$. For any $p\in[0, 1]$ such that $R_p > 1$, the value of the effective susceptible class at the endemic equilibrium is given by

Proof. If $\varepsilon_L = 0$ then from (2.3), $S_E^* = S^* = \frac{1}{R_0}$. If $\varepsilon_L > 0$ then again from (2.3), $S_E^* = S^*+\varepsilon_L V^* = \frac{1}{R_0}-\varepsilon_L V^*+\varepsilon_L V^* = \frac{1}{R_0}$.

3.3.   Global stability of endemic equilibrium when $R_p > 1$

Before we formally define different measures of the honeymoon period, we first prove the global stability (for all non-negative initial conditions) of the endemic equilibrium when $R_p > 1$. Here we closely follow the work of Li et al. ^[8] using the LaSalle Invariance Principle.

Proposition 5 (Global asymptotic stability of endemic equilibrium). Let $R_p > 1$. For any set of positive initial conditions such that $V(0)+S(0)+I(0)+R(0) = 1$, the trajectory of the solution to (2.1) approaches the endemic equilibrium $(V^*, S^*, I^*, R^*)$ given by (2.3)-(2.6).

Proof. If $R_p > 1$ then $(V^*, S^*, I^*, R^*)$ is in the positive cone ^[1]. Let

The system of Eqs (2.1) is equivalent to the system given by,

The stability of the endemic equilibrium is maintained after transformation and is given by $(x^*, y^*, z^*) = (1, 1, 1)$. Define the function $L(x, y, z)$ by,

where the positive numbers $a_1, a_2$, and $a_3$ will be found later. Differentiating $L$ with respect to $t$ along solutions of (3.3) yields,

Following Li et al. ^[8], we determine the constants $a_1, a_2$, and $a_3$ and rearrange $F (x, y, z)$ in four steps:

Step 1. We construct Table 2 with the terms $x$, $\frac{1}{x}$, $y$, $\frac{1}{y}$, $\frac{y}{x}$, $z$, $xz$ and $yz$ (which all appear in $F(x, y, z)$) as the column names. Then, we identify groups of terms that multiply to one and indicate each group by a row with a check mark on the relevant terms. We also write $p_k : = b_k(n_k - h_{k, 1} - h_{k, 2} - \dots - h_{k, n_k})$ in the last column, where $h_{k, i}$ are the terms that multiply to 1 (i.e., $\prod_{i = 1}^{n_k} h_{k, i} = 1$), $n_i$ is the number of terms and $b_k$ is a non-negative number to be found later.

We now rewrite $F(x, y, z)$ as,

Step 2. Set the terms multiplying the unmarked terms in Table 2 for function $F (x, y, z)$ to zero. This means we set the terms multiplying $z$, $xz$, $yz$ all equal to zero:

A solution to this system is $a_1 = a_2 = a_3 = 1$. Using this, the function $L(x, y, z)$ is positive definite. Substituting these values into $F (x, y, z)$ and using the fact that $(V^*, S^*, I^*, R^*)$ is an endemic equilibrium (substituting them into the right-hand-side of (2.1) yields zeros) to simplify terms we derive,

Step 3. Let the coefficients for the same terms between $F (x, y, z)$ and $\sum_{k = 1}^{5} p_k$ be equal, which yields the following equations:

We note that this is a linear system of five equations with five unknowns and rank of four. We can solve for the $b_1$ through $b_4$ in terms of $b_5$,

Note that we did not use $b_1 + b_2 = \beta S^* I^* + \mu S^*$ in (3.6) to get (3.7). Substituting in the expressions for $b_1$ and $b_2$ to this equation yields,

This equation indeed holds from endemic equilibrium, so this verifies that the system of Eqs (3.6) is consistent. To guarantee that all the $b_i$'s are nonnegative, $b_5$ should satisfy the following inequality:

We can further simplify this by noticing from the endemic equilibrium equations that

which implies that $\varepsilon_L\beta V^*I^* +\mu V^* < (1-\varepsilon_A)p\mu$. Additionally, we can rearrange the equation to yield,

which implies that $(1-\varepsilon_A)p\mu-\alpha V^* = \varepsilon_L\beta V^*I^* +\mu V^* \geq\varepsilon_L\beta V^*I^*+\mu V^*-(1-(1-\varepsilon_A)p)\mu$. Thus, the inequality (3.8) simplifies to $\max\Big\{0, \varepsilon_L\beta V^*I^*+\mu V^* \Big\}\leq b_5 \leq \varepsilon_L\beta V^*I^*+\mu V^*$ which yields,

Step 4. From (3.7),

The last line comes from substituting in $\varepsilon_L \beta V^* I^* +\mu V^* = (1 - \varepsilon_A) p\mu-\alpha V^*$ which comes from the equations for the endemic equilibria. This step verifies that the constant terms in $F (x, y, z)$ and $\sum_{k = 1}^{5} p_k$ are equal.

Finally, by putting together Steps 1-4 we find that the function $F (x, y, z)$ can be rewritten in the form of $\sum_{k = 1}^{5}p_k = \sum_{k = 1}^5 b_k(n_k - h_{k, 1} - h_{k, 2} - \dots - h_{k, n_k})$ with

Thus,

The property that the arithmetic mean is greater than or equal to the geometric mean implies that $F(x, y, z) \le 0$ for all positive values of $x$, $y$ and $z$, and that the equality only holds only for $x = y = 1$. This means that $\{(x, y, z) \in\mathbb{R}^3_+: F (x, y, z) = 0\} = \{(x, y, z) : x = y = 1, z > 0\}$. The maximum invariant subset of this set only contains $(1, 1, 1)$. By LaSalle's Invariance Principle, this equilibrium is globally asymptotically stable for all initializations in the positive cone. It follows that the endemic equilibrium $(V^*, S^*, I^*, R^*)$ is also globally stable for all non-negative initial conditions.

3.4.   Technical definition of the honeymoon period

When $R_p > 1$, the global stability of the endemic equilibrium (Proposition 5) and the fact that $S_E^* = \frac{1}{R_0}$ independent of the value of $p$ (Proposition 4) motivates us to introduce the following technical definition for the honeymoon period in Definition 6.

Definition 6. Let $p > 0$ and $R_p > 1$. Suppose that the initial conditions of the model Eqs (2.1) are the pre-vaccine era endemic equilibrium (the endemic equilibrium of the system when $p$ is set to zero). This means,

Define the length of the honeymoon period of the system to be

if this exists and $\mathcal{T} = \infty$ otherwise. The relative drop in the effective susceptible class is defined to be

if this exists.

Under the conditions given by Definition 6, we say that the honeymoon period is the time interval $(0, \mathcal{T})$. During this time, we expect that the effective number of susceptibles in the population is below its equilibrium value of $\frac{1}{R_0}$ and that the infectious fraction of the population is decreasing. These statements are proven later in Proposition 7.

We plotted the length of the honeymoon period in Figure 6 for different values of transmission rate $\beta$ and vaccine failure parameter ($\varepsilon_L$, $\varepsilon_A$ and $\varepsilon_W$ for the purely leaky, all-or-nothing and waning models respectively). The length of the honeymoon period is shown using a color bar in a log-scale. This plot shows that the waning model generally has longer honeymoon periods over the two other models.

We also plotted the relative drop in the effective susceptible class ($S_E$ in Definition 6) after the start of mass vaccination in Figure 7. We note that the waning model generally displays a larger relative drop in effective susceptibles than the other two models. The border between the colored and white area corresponds to the region where $R_p = 1$ and a transcritical bifurcation occurs with the stability switching between the endemic and disease-free equilibria.

4.   Discussion

The results in Figures 4 show that reactivity and amplification envelope ^[5], two measures of transient dynamics based on linearization techniques, may not clearly distinguish which types of vaccine models will display significant honeymoon periods. This is shown when we computed these quantities for models and parameter regimes that exhibit very clearly different honeymoon periods in their numerically computed trajectories shown in Figures 3 and 5. In both the all-or-nothing and waning models we have $\varepsilon_{L} = 0$ which sets several elements in the Jacobian (2.9) to zero. Additionally, these two models have the same values at endemic equilibrium resulting in very similar values for their reactivity and amplification envelopes. We also note that the drop in the reactivity for the leaky model in Figure 4(a) occurs near the so-called reinfection threshold where $R_0 = \frac{1}{\varepsilon_L}$ ^[10]. This is a subject we are looking into for further study.

In Figures 3 and 5 the drop in incidence that occur during honeymoon periods follow an initial drop in the effective number of susceptibles in the population. This drop in incidence eventually results in the effective number of susceptibles slowly building up again. At the end of the honeymoon period, disease outbreaks may occur due to a buildup in effectively susceptible individuals during the low incidence period. Eventually the effective susceptible population recover to its equilibrium value of $\frac{1}{R_0}$ which is the same value independent of the vaccine coverage and type of vaccine failure (Proposition 4). This allowed us to give definitions for two features of the honeymoon period: its length $\mathcal{T}$ and the relative drop in the effective susceptible population $S_E$ (Definition 6). These quantities are measured in Figures 6, 7. From these figures we clearly see that the waning model has a propensity for longer and more pronounced honeymoon periods than the leaky and all-or-nothing models. The length of the honeymoon period depends on the model type and parameter values, as we can see in the Figure 6. The longest honeymoon periods seem to occur close to the boundary where $R_p = 1$.

As discussed earlier, we expect that during the honeymoon period the effective number of susceptibles in the population is below its equilibrium value and that the infectious fraction of the population is decreasing. We establish these basic properties of the honeymoon period and the behavior of solutions to (2.1) during this period in Proposition 7.

Proposition 7. Let $p > 0$ and $R_p > 1$.

1. $S_E(t)\in(0, \frac{1}{R_0})$ for $t\in(0, \mathcal{T})$.

2. $\frac{dI(t)}{dt} < 0$ for $t\in(0, \mathcal{T})$.

Proof. To prove Part 1 we first solve for $\frac{dS_E}{dt} = \frac{dS}{dt}+\varepsilon_L\frac{dV}{dt}$ and group some terms together to get,

We find that this derivative must be negative at $t = 0$ by using the values of the prevaccine era endemic equilibrium as initial conditions (Definition 6).

Thus there exists $\delta > 0$ such that $\frac{dS_E}{dt} < 0$ for all $t\in[0, \delta)$, implying that $\mathcal{T} > 0$ and the honeymoon period interval is nonempty. Since we initialize the system with $S_E(0) = \frac{1}{R_0}$, by the definition of $\mathcal{T}$ we must have $S_E(t) < \frac{1}{R_0}$ for $t\in(0, \mathcal{T})$. The invariance of the system (2.1) in the positive cone implies that $S_E(t)\in(0, \frac{1}{R_0})$ for $t\in(0, \mathcal{T})$.

We note also here that from Proposition 4 the endemic equilibrium with or without vaccination is $S_E^* = \frac{1}{R_0}$. Thus $S_E(t)\to \frac{1}{R_0}$ by the global stability of the endemic equilibrium (Proposition 5). We note that in many numerically generated trajectories, we have found that $S_E(t)$ approaches $\frac{1}{R_0}$ in an oscillatory manner, implying that $\mathcal{T}$ is finite. Looking at the properties of how the vaccine era endemic equilibrium is approached is a topic for future study.

To prove Part 2, we rewrite the equation for the derivative of $I$ in (2.1) as,

Then $\frac{dI}{dt} < 0$ for $t\in(0, \mathcal{T})$ follows from Part 1.

Part 2 of Proposition 7 imply that if the disease is going to rebound, this can only occur after the honeymoon period. The $I$ class is continuously decreasing during the entire honeymoon period.

5.   Conclusion

It is difficult to model childhood diseases and vaccination ^[9]. There is a lot of interest in modeling these correctly because vaccination is one of the most powerful and cost-effective interventions for diseases. While most studies focus on the equilibrium behavior of these systems, measuring and understanding the transient behavior is also important as they can last for very long periods of time. Transient behaviors, such as honeymoon periods can also provide clues into the type of imperfect vaccine that is being used, and that can be very important for the design and cost-benefit evaluation of the rollout of vaccination for some diseases.

In this manuscript we reviewed the VSIR model analyzed in ^[1] which allows for three modes of vaccine failure: leaky, all-or-nothing, and waning. It has previously been observed that, given the same value of vaccine impact, the all-or-nothing and waning models have the same endemic equilibrium, so it is not possible to distinguish between them at equilibrium. However, by studying their transient dynamics numerically, especially the dynamics after the start of mass vaccination, we often observe very clear differences in how these models approach the vaccine era endemic equilibrium. For some parameter values, the waning model may display a very pronounced honeymoon period while the all-or-nothing model with the same vaccine impact displays a less eventful drop in incidence to its new equilibrium value.

We attempted to use linearization methods in ecology to characterize which models can display honeymoon periods. We reviewed the reactivity and amplification envelope ^[5], however the results show that these cannot distinguish which of the all-or-nothing and waning models will exhibit significant honeymoon periods. So far, we can only confirm if a model displays honeymoon periods when we numerically solve for the trajectory of the models.

As a basis for future study on honeymoon periods, we presented a technical definition of some features of this transient period in Definition 6 and proves some properties of the solutions during this time interval in Proposition 7 (using the global stability of the endemic equilibrium in Proposition 5). These features were (1) the length of the honeymoon period (plotted in Figure 6 for different values of transmission rate and vaccine failure parameter) and (2) the relative drop in the susceptible class (plotted in Figure 7). From this we saw that the waning model usually has the largest relative drop in the effective susceptible class after the start of mass vaccination, and the longest honeymoon period. Interestingly, the trajectory of the effective number of susceptibles for the leaky model is similar to that of the all-or-nothing model.

There is still a lot to learn about honeymoon periods. The results of this study indicate that they may not be detected by the reactivity and amplification envelope, measures of transient dynamics that are based on linearization techniques. Further study of the honeymoon periods will probably require nonlinear techniques.

Acknowledgments

We would like to thank Queen's University, the University of Manitoba and NSERC for funding this research. We would also like to thank J. Arino, S. Portet, P. Taylor, T. Day, P. Rohani and A.A. King for their helpful comments and suggestions that have helped us with the development of ideas in this manuscript.

Conflict of interest

We have no conflicts of interest to disclose.