Epidemiological impact of waning immunization on a vaccinated population

This is an epidemiological SIRV model based study that is designed to analyze the impact of vaccination in containing infection spread, in a 4-tiered population compartment comprised of susceptible, infected, recovered and vaccinated agents. While many models assume a lifelong protection through vaccination, we focus on the impact of waning immunization due to conversion of vaccinated and recovered agents back to susceptible ones. Two asymptotic states exist, the"disease-free equilibrium"and the"endemic equilibrium"; we express the transitions between these states as function of the vaccination and conversion rates using the basic reproduction number as a descriptor. We find that the vaccination of newborns and adults have different consequences in controlling epidemics. We also find that a decaying disease protection within the recovered sub-population is not sufficient to trigger an epidemic at the linear level. Our simulations focus on parameter sets that could model a disease with waning immunization like pertussis. For a diffusively coupled population, a transition to the endemic state can be initiated via the propagation of a traveling infection wave, described successfully within a Fisher-Kolmogorov framework.


Introduction
Infectious diseases have a strong impact on the dynamics of human populations and are routinely highlighted when epidemic outbreaks of deadly infections like Ebola or MERS occur. Increased human mobility, the rise of pathogens resistant to antibiotics ("Antimicrobial Resistance"), and the advent of new, so called Emerging Infectious Diseases are making infectious diseases a major health challenge of the 21st century.
To analyze transmitting diseases and epidemic outbreaks and to inform public health organisms, a wide range of mathematical models have been proposed and studied in detail. One of the simplest and well-known models illustrating the dynamics of epidemics is the SIR model, proposed by Kermack and McKendrick in 1927 [1]. The central idea of this model is to divide the entire population into three separate groups: susceptible (S) individuals that have never been infected and are not immune to the disease; infected (I) individuals, who are infectious and can spread the disease within the population; and recovered (or removed; R) individuals who have already had the disease and are therefore immune for life.
Temporal changes in the numbers of individuals in different groups of the SIR model are described by the differential equations The number of infected agents increases proportionally to the number of susceptible agents times the force of infection βI, being itself the product of the non-zero infection rate β and the number of infectious individuals. This is an example of a density dependent force of infection, being an alternative a frequency dependent force of infection βI/N [2]. Individuals recover from the infection with rate a.
Adding Eqs. (1a-1c), one can see that S + I + R is time conserved. Furthermore, there is no notion of physical space in the model, meaning that individuals are uniformly mixed in the population. Both properties are common assumptions of simple SIR models that are rarely realized in real populations.
One of the key medical advances of the 20 th century has been the proliferation of cheap and safe vaccines for a range of diseases. Vaccines are among the most important tools to prevent infections and epidemics. One particular beneficial and mathematically interesting aspect of vaccination is that not 100% of a population need to be immunized in order to prevent epidemic outbreaks, a property called herd immunity [2]. Unfortunately, vaccination encounters opposition in different socio-cultural contexts that endangers the working of the herd immunity effect as the fraction of vaccinated individuals falls below a threshold. Yet another effect that limits the efficiency of vaccines is that the protection provided has a finite lifetime: depending on the disease, the immunization effect of the vaccine fades over time and the patient has to be re-vaccinated (e.g., diphtheria, tetanus and pertussis).
In this article, we investigate an SIRV model ("V" for vaccination) that accounts for the fact that immunization of the vaccine wanes over time and that even recovered individuals can fall ill again. Furthermore, we consider the spatiotemporal dynamics of such a system. Epidemiological and SIRV models have been considered in many variants, for some reviews see Refs. [3,4,2], and some recent work consider the dynamics on networks [5,6,7], information-driven vaccination [8,9,10], or stochastic behavior [11,12]. Spatiotemporal dynamics in nonlinear systems often show traveling wave patterns or Turing-like, stationary patterns [13,14]. In the context of spatial epidemiological models [15], spatial coupling is often described by reaction-diffusion equations or networks [16,17,18,19,20,21,22].

The SIRV model with finite lifetime protection
Among the range of epidemiological models using the SIR model as a building block, we focus on those aimed to investigate the impact of vaccination, using here a modification of the SIRV models presented in Refs. [4,23]. In this model, an independent birth rate B and death rate p are taken into account. There are two types of vaccination: v 1 is the fraction of newborns being vaccinated and v 2 is the rate of vaccination of susceptible individuals. By construction, 0 ≤ v 1 ≤ 1, with these two limiting cases representing that all newborns are either susceptible or vaccinated. In contrary to the classical SIR model, protection against infection is not for life: recovered individuals become susceptible again Conversion rate from recovered to susceptible q 2 Conversion rate from vaccinated to susceptible Table 1 Parameters used in model (2). The parameters β, a, p and B are positive, the others non-negative and 0 ≤ v 1 ≤ 1. All parameters are measured per unit time, except for v 1 (unitless) and β (per unit time times population size unit).
with rate q 1 , vaccinated ones with rate q 2 , yielding the model A schematic representation of the transitions among the compartments in the SIRV model is displayed in Fig. 1 and the meaning of the parameters is found in Table 1. All parameters and variables have to be non-negative for Eqs. (2) being interpreted as a valid epidemiological model. From Eqs. (2a-2d), and defining the total population size as N (t) = S + This linear equation has the solution where N 0 is the initial population size. While the total population size can vary as a transient if N 0 = B/p, in the limit lim t→∞ N = B/p which means that asymptotically the total population size is constant. In particular, if birth and death rates are equal, then lim t→∞ N = 1. In the following, we require both birth and death rates to be strictly non-zero to exclude cases of vanishing or exponentially growing populations. Using Eq. (4), the model in Eqs. (2) can be interpreted as a time-dependent system with 3 variables (using, e.g., R = N − S − I − V ). However, in the simulations, we numerically integrate the system defined in Eq. (2) directly.

Stationary states and their stability
A key motivation of the study of population models is to find out the number of possible asymptotic states and to evaluate the relative size of the subpopulations. We find these steady state solutions in our model by setting the left-hand sides of Eqs. (2) to zero.
The first stationary state is characterized by the following solution for susceptible, infected, recovered and vaccinated agents: The number of infected individuals (and hence the number of recovered) is zero in this solution, justifying the name disease-free equilibrium (denoted with subscript "1" and abbreviated as DFE). This solution describes a population without disease where the parameters control the relative fractions of susceptible and vaccinated individuals (summing up to N = B/p). S 1 and V 1 are both independent of q 1 , meaning that in the disease-free state the finite time of disease protection of a recovered individual is irrelevant (which is consistent with R 1 = 0). Neither the infection nor the recovery rate influences the steady state values S 1 and V 1 . Note that S 1 cannot be negative since 0 ≤ v 1 ≤ 1 (see above).
The second stationary state is characterized by the following solution for susceptible, infected, recovered and vaccinated agents: Since the number of infected individuals (and hence the number of recovered) is non-zero, this state (denoted with subscript "2") is referred to as the endemic equilibrium (EE). Again, only the parameters control the relative fractions of individuals in the different compartments (summing up to N = B/p). To describe a state relevant from a population dynamics point of view, the four compartments must have non-negative population numbers, constraining the parameter values, as we see later.
To perform a standard linear stability analysis, we determine the Jacobian of Eqs. (2) and obtain where I * and S * need to be replaced by the respective stationary state solutions. Evaluated at the disease-free equilibrium, we obtain the eigenvalues as The first three eigenvalues are always negative since all rates are non-negative and p is strictly positive. The fourth eigenvalue can change sign (and therefore indicate instability of the solution), depending on the values of all parameters except q 1 , that does not influence the stability of the disease-free state at the linear level. Setting λ 4 = 0, we obtain that the DFE is unstable if the infection rate β is larger than a critical infection rate β c , given by For β = β c , one can show that S 2 = S 1 while I 2 = 0. This is schematically shown in Fig. 2. Obviously, this condition can also be expressed as critical values for the other parameters (illustrated below). The condition λ 4 = 0 coincides with the parameter value for which the DFE and the EE are identical. λ 4 can also be written as β(q 2 + p)(a + q 1 + p)(q 2 + v 2 + p) −1 I 2 which shows that the existence of the endemic equilibrium is associated with the instability of the disease-free equilibrium. The stability of the endemic equilibrium can be checked for in an analogous way but is omitted here as they lead to very lengthy expressions that have to be evaluated numerically. The existence of a stable DFE does not exclude the possibility of an appropriate initial condition mediated epidemic outbreak via a transient increase in the value of I.

Transition between states and basic reproduction number
In this section, we show how the stationary states vary as a function of some of the parameters. In particular, we consider the vaccination parameters v 1 and v 2 and the conversion rate q 2 (we have seen above that q 1 does not influence the existence or change of stability of the DFE). Figure 3 shows the stationary state solutions for all four sub-populations as a function of the fraction of vaccinated newborns (v 1 ). For the DFE, the dependence on v 1 is linear for S 1 and V 1 , showing the direct proportionality of the fraction of vaccinated people in the population on the fraction of vaccinated newborns. As v 1 is decreased. the DFE becomes unstable at a critical v 1c via a transcritical bifurcation and the EE sets in, a general feature of SIR models [3]. Then, the number of susceptible remains constant in the population while the number of infected (and also recovered) increases linearly. At the same time, the vaccinated fraction of the population decreases, and with a higher rate than when the DFE was stable.
We now consider the case of varying the vaccination rate of the susceptible individuals v 2 (Fig. 4). Considering first the EE, it can be seen that the qualitative behavior of the curves is similar to the case of varying the fraction of vaccinated newborns. However, for the DFE the fractions of susceptible and vaccinated sub-populations do not change linearly as above, see also Eq. (5). In particular, the rate of increase of V 1 as a function of v 2 starts slowing down beyond the linear regime, meaning that it becomes increasingly difficult to protect the population. Also, for v 2 = 0, V 1 = V 2 and hence if the only vaccination is taking place at birth, the fraction of vaccinated people is identical in the endemic and disease-free states.
Finally, we discuss the case of changing the conversion rate from vaccinated to susceptibles (q 2 ). The loss of protection of the vaccination plays an antagonistic role to the vaccination rate. It is not a surprise to find that for the DFE, the vaccinated fraction of the population decreases as q 2 increases. The role of q 2 in the equation for S 1 is the same as v 2 in the equation for V 1 [Eq. (5)]. For the EE, though, the situation is slightly different. While the infected fraction of the population increases with q 2 , it does so at rate that is slower than the linear growth rate at onset, showing that a waning immunization favors the endemic state, but a change in this parameter is less dangerous than a decrease of any of the vaccination parameters.
A relevant quantity in epidemiology is the basic reproduction number R 0 . It is defined as the expected number of secondary individuals infected by an individual in its lifetime (for a review see Ref. [24]). By analyzing this value it is possible to predict whether a disease present in a population will create an epidemic (if R 0 > 1).
To calculate the basic reproduction number R 0 , we use the next generation method for structured populations [24]. For that we separate the Jacobian given in Eq. (7) into a transmission part T and transition part Σ, evaluated at the DFE. We obtain: and Then, R 0 is the leading eigenvalue of the matrix [−T Σ −1 ]. It is determined as The R 0 shown in Eq. (12) above is identical to S 1 /S 2 and to β/β c , providing alternative interpretations of the onset of an epidemic. Also, R 0 = 1 + (a + p) −1 λ 4 , elucidating the relationship between the basic reproduction number and the dominant eigenvalue of the stability analysis of the DFE. Because of this link it is not surprising that for this model, the same result can be obtained by evaluating λ 4 or by setting I 2 to zero. Figure 6 shows R 0 as a function of the vaccination parameters v 1 , v 2 , and the rate of loss of protection q 2 . The panels A-C exhibit a situation involving an epidemic with low R 0 , while panels D-F use parameter values for pertussis, a disease with high R 0 . In agreement with the above figures for the stationary states, we observe that for low vaccination rates (v 1 , v 2 ) and a high rate of loss of protection (q 2 ), the endemic state is stable while the disease-free state is unstable. On the other hand, if the vaccinated fraction of the population fast loses its protection, a transition from the disease-free state to the endemic equilibrium occurs. Only the dependence of R 0 on v 1 is linear. As the infection rate β increases, the R 0 curves are shifted upwards (panels A-C), reflecting an increased tendency to destabilize the disease-free equilibrium. For the specific case of pertussis, we observe that the curve R 0 (v 1 ) is relatively flat (for three different combinations of values of v 2 and q 2 ), meaning that even complete vaccination of a newborn population is not sufficient to contain the disease if the vaccination rate v 2 is not high enough. This is confirmed by the curve R 0 (v 2 ) (Fig. 6E) which shows a sharp decay, illustrating that vaccination is an efficient way of decreasing R 0 (for three different combinations of values of v 1 and q 2 ). It is worth noting that red and blue curves differ only in v 1 and while the values are very different (0.95 and 0.05), the position of the curves are similar. In contrast to this, increasing q 2 from 0.05 to 0.3 make the control via v 2 more difficult (purple curve). Finally, the curves in Fig. 6F show that R 0 is also very sensitive to the conversion rate from vaccinated to susceptible as only small values allow that R 0 stays below 1.
As mentioned above, the stability analysis of the endemic equilibrium leads to very lengthy expressions that we exclude for the sake of brevity. We, hence, assess the stability of the endemic equilibirium numerically. In Fig. 7 we show how the four eigenvalues vary as a function of the main parameters q 1 , q 2 , v 1 and v 2 . As we know that the EE only exists if the DFE is unstable, we also plot the R 0 curve indicating the critical parameter values. The fundamental result is that where R 0 > 1, the real parts of the four eigenvalues are negative, showing that the EE is stable in these parameter regions. A particular case is the graph Fig. 7A which confirms that q 1 is not only irrelevant for the stability of the DFE, but also of the EE.
While the existence and stability of the asymptotic states are fundamental properties of any epidemiological system, an epidemic is a time-dependent process. For example, even if the DFE is stable and asymptotically obtained, an epidemic outburst can occur. Hence, in a deterministic system, the initial state of the population is fundamentally important. In the presence of birth-death processes, a high death rate can mask slow processes (if the loss of immunization is on the timescale of life expectancy) or an expanding population may require a higher vaccination rate in order to keep the population protected.

The spatiotemporal SIRV model
Epidemiological models without spatial degrees of freedom can only be applied to very well mixed populations. However, people live in confined communities that are spatially connected. As a starting approach, we assume that the population is distributed over a one-dimensional space where transport between adjacent areas is diffusive (equivalent to nearest-neighbor interactions). Therefore, we add diffusion terms to the SIRV model (2) and obtain where D F (F = S, I, R, V ) are diffusion constants for susceptible, infected, recovered and vaccinated individuals, respectively.
The two fixed points shown in Eqs. (5) and (6) of the diffusion-free system are solutions of the system with diffusion (13) in case the variables do not show any spatial dependence, i.e., represent a homogeneous solution. However, the linear stability of this homogeneous solution depends on diffusion, as we shall see now.
Perturbed around the homogeneous fixed points, in the Fourier transformed (k, t) space, the dynamics is represented through the Jacobian J k that is given by where I * and S * need to be replaced by the respective stationary state solution [Eqs. (5,6)]. Around the disease-free steady state, one can find the eigenvalues as follows The eigenvalue λ 1 (k) is a generalization of λ 2 of the ODE system and is always negative. The eigenvalue λ 2 (k) is the generalization of λ 4 of the ODE system and can therefore change sign. The eigenvalues λ 3,4 (k) depend on the sum and differences of the diffusion coefficients of the susceptible and vaccinated population fractions. It can be easily shown that λ 3,4 (k) are always negative and hence diffusion has always a stabilising effect on the DFE. The most unstable wavenumber is k = 0. Hence, whenever the DFE is unstable in the purely temporal system, it is also unstable in the spatiotemporal system. Figure 8 shows a numerical simulation of the spatiotemporal SIRV model [Eqs. (13)] in one-dimensional space. The initial condition is a disease-free state with a small nucleus of infected agents at the center of the medium. Parameters have been chosen to ensure that the disease-free state is unstable, leading to a transition to the endemic state. This can be clearly seen as a traveling wave in the space-time plot for I in Fig. 8A. Figures 8B and 8C illustrate the behavior of all variables for a fixed point in space (B) and for a fixed point in time (C). The latter displays the profile of the traveling wave front. For this set of parameters, the spatial distribution for I shows small peaks in the fronts.
The wave of infection observed in Fig. 8 can be investigated in more detail. In Fig. 9, we show how the front velocity changes with the diffusion constant D I and the infection rate β. Both functional forms follow a square root dependence reminiscent of the Fisher-Kolmogorov equation [14,25]. Indeed, for a singlespecies population model with variable u, it is known that the natural front velocity of a front triggering a transition from the unstable to the stable state is given by v = 2 f (u 1 )D, where D is the diffusion constant, f (u) describes the temporal dynamics and u 1 is the unstable steady state [25]. Applying the same rationale to Eq. (2b), we obtain The qualitative agreement between the curves is surprisingly good which is remarkable as no fitting parameters are applied and the analytic expression uses only one equation of a coupled 4-dimensional dynamical system. There is a slight quantitative difference for small velocities, as seen in Fig. 9A that could be partially explained by the fact that the simulations are performed in a finite sized system and that the calculation of the front speed from the simulation data carries an error.

Discussion
In this article, we have considered an SIRV model in the temporal and spatiotemporal domain. The model has two asymptotic states, the disease-free state and the endemic state. We have focused on the consequences of diminishing immunization, i.e., the effect when vaccinated or recovered individuals become susceptible again. The results have been obtained through bifurcation analysis of the individual solutions (for S, I, R and V ), as well as through the determination of the basic reproduction number R 0 . In the asymptotic regime the number of each sub-populations is proportional to its density in the whole population, so the results refer directly to population densities or fractions. Our exclusively temporal model shares similarities with a model studied in [23], however, the models only coincide if we set v 1 = q 1 = 0 in our model and simultaneously set µ = σ = 0 in the model discussed in Ref. [23]. However, assuming non-zero values for these parameter is crucial for both our model (possible vaccination at birth and conversion from recovered to susceptible agents) and the model discussed in Ref. [23] (variable vaccine efficacy and possibility of disease-induced deaths) and hence the interpretation and applicability of the models differ substantially. By considering the results of a linear stability analysis of the disease-free state, we have found that the loss of protection of the recovered fraction of the population (with rate q 1 ) has no influence on the onset of the endemic state. While the rate q 1 does not influence the asymptotic DFE, it can impact on the transient time to equilibrium. On the other hand, the loss of protection of the vaccinated fraction of the population (with rate q 2 ) can shift the population from a disease-free state to an endemic one. An interesting feature of this model is that the density of susceptibles in the endemic regime does not depend on q 2 . The curve of R 0 with q 2 is increasing, however, with a decreasing slope, meaning that decreasing q 2 in the epidemic regime may bring the population closer to the threshold than predicted by a linear regression. Considering the effect of the vaccination rates, we find that the fraction of vaccinated newborns v 1 changes the asymptotic fractions linearly in both stationary states, as well as R 0 . This is in contrast to v 2 where the dependence is nonlinear. There, we have found that if v 2 is decreased in an disease-free state, the basic reproduction number increases more strongly than predicted by a linear regression. This implies that the critical R 0 = 1 may be reached for higher v 2 than assumed.
Our results show that in the diffusion-free state the dominant eigenvalue of the disease-free state λ 4 = β(q 2 + p)(a + q 1 + p) q 2 + v 2 + p I 2 . This means that if the endemic state exists, the disease-free equilibrium is unstable that is associated with positive λ 4 values and R 0 > 1. In the complete absence of adult vaccination, implying v 2 = 0, Eqs. (5) and (6) show that the vaccinated number density is the same for the two asymptotic states Bv 1 q 2 + p . This then implies that one cannot predict the actual epidemic state from the proportion of the vaccinated population alone. We have presented a numerical solution for the stability problem of the endemic equilibrium. It indicates that while the EE exists, it is stable.
The features obtained from a study of this model can be put in the appropriate context of epidemiological data. Diphtheria and Pertussis (whooping cough) are amongst the diseases that are associated with waning immunization. Repeated vaccinations ("boosts") are needed to prevent the spreading of such diseases. Due to the high R 0 values of these diseases, children are vaccinated at early ages. Without epidemiological control, the R 0 of pertussis has been estimated at 16-18 [26], a value lowered to 11-15 [27] later. In the presence of vaccination, the value could be lowered to around 5.5 [28]. The incidence among adults are explained by waning immunization and the possibility of evolving subclinical strains that are held responsible for persistence of pertussis in vaccinated populations [28]. In Fig. 10, we show a short time series of a population suffering from pertussis infection and for which the endemic equilibrium is stable. The initial state consists of a population with very few infected agents. We clearly notice some outbursts of infection, with a characteristic time gap between 1 and 2 years. This timescale is not far off from known deterministic models of pertussis which consistently predict annual epidemics [29]. Note, however, that detailed and more realistic models for pertussis rely on an SEIR mechanism, with an exposed/latent phase and/or age structure, and possibly term-time forcing. Furthermore, stochastic effects are also known to be crucial in the disease dynamics [30]. A recent work compares the different classes of models including reinfection of recovered and loss of infection-derived immunity and subsequent reinfection [31]. In the context of this article, we simply want to illustrate an example of a specific disease for our model.
For all realistic epidemiological models, spatial interactions have to be considered. In our model, we have assumed a nearest-neighbor interaction, modeled by diffusion terms. Based on the definition of a spatial basic reproduction number R 0k , our linear stability analysis of the DFE confirms that the most unstable wavenumber is k = 0, and that the disease-free equilibrium cannot be destablized by controlling the diffusion rates. A discussion of the spatial stability problem of the endemic equilibrium is beyond the scope of the present work.
A well-known feature of infection models with diffusion is that they are able to describe the propagation of waves, of particular interest being waves that represent the onset of an epidemic. We have shown that in spite of the comparatively high complexity of the model (4 coupled equations), the wave speed still approximately follows the one-species Fisher-Kolmogorov model, similar to what has been observed for a different model [16].
Temporal and spatiotemporal epidemiological models have been studied in many variants. A series of recent works try to find optimal vaccination strategies, for example by a probabilistic modelling of infection in networks [32], by minimizing the number of infected and susceptibles [33], by a Poisson distributed vaccination schedule on networks [6], by an information (and time) dependent vaccination rate [10], or by optimizing the vaccination rate through a stochastic maximum principle [12]. In contrast to these articles, we analyze the front speed of a general SIRV model, similar to the approaches of [18] for a stochastic SIR model and [20] for an SIR model with non-smooth treatment (vaccination) functions.
As an outlook to further work in the model, we mention the spatial stability analysis of the endemic equilibrium, an analytical and numerical investigation of fronts in two space dimensions and the incorporation of social effects. Fig. 2 Illustration of the stability of the disease-free state (I 1 = 0). Note that I 2 < 0 corresponds to a non-physical solution and hence absence of the endemic equilibrium state (population numbers have to be positive).