Rebound in epidemic control: how misaligned vaccination timing amplifies infection peaks

Global pandemics of infectious diseases are one of the greatest threats that humanity faces and, as a consequence, the issue of how to control them becomes increasingly more urgent as time passes. To rise up to the challenge, traditional epidemic control offers two alternatives: non-pharmaceutical interventions (NPIs) such as school closures, quarantine, or mask mandates on the one hand, and pharmaceutical interventions such as vaccination, drugs, or treatments on the other. However, in the case of NPIs¹ as in that of vaccination^2,3, apparently similar implementations can lead to different outcomes, hinting at the fact that epidemic control is, at its heart, a complex problem that requires a detailed understanding of the underlying mechanisms of the epidemic.

Mathematical models have long been used to provide a framework to gather such insights⁴. In the context of COVID-19, for instance, they have proved to be a powerful tool for the description, prediction and prevention of the ongoing pandemic^{5,6,7,8,9,10,11,12,13,14,15}. Therefore, it comes as no surprise that the search for an optimal strategy for the minimization of the cases is a well-researched topic, in particular in the context of NPIs^{16,17,18,19,20}. The reason behind this interest is obvious: non-pharmaceutical interventions come with high social and economic costs attached and it is therefore crucial to make them last as little as possible, possibly to contain the epidemic until pharmaceutical treatments are available. The same kind of optimization problem has been studied in relation to vaccination, for example in connection to issues of economical costs²¹, vaccination rates²², age prioritization²³, boosters’ distribution^13,24,25 and NPIs-vaccination synergy²⁶. Following this literature, in this paper, we will also focus on the optimization of vaccination strategies, especially from the point of view of timing and duration. Unlike the previous works on this subject^27,28,29, we employ a SIRS model with a time-limited vaccination campaign, which stands as the simplest model incorporating waning immunity (a feature shared by vaccines of several respiratory illnesses^30,31). Additionally, we assess the effectiveness of a vaccination strategy by measuring the peak prevalence observed after the conclusion of the vaccination campaign.

From our analysis, a counter-intuitive rebound effect is observed: the timing of vaccination and subsequent immunity waning can synchronize with the increase in susceptible population, potentially leading to an infection peak larger than what would be expected in the absence of vaccination. Earlier attempts to model reinfections and waning immunity employing a SIRS model exist^32,33, but they primarily deal with the stability of such systems under the assumption of an unlimited supply of vaccines, and therefore no rebound was observed. The resurgence of infected individuals after the start of a mass vaccination campaign has been documented and modeled in various papers^34,35,36, some of them pinpointing the gradual accumulation of the susceptible compartment as a key driver³⁷. Additional factors include the significance of transient phases in the dynamics³⁸ and the impact of declining immunity within hosts³⁹. In contrast, investigations into time-limited vaccination strategies are relatively rare, often focusing on specific scenarios such as pathogens with seasonal infectivity⁴⁰ or vaccinations administered in pulsatile campaigns⁴¹. However, to the best of our knowledge, the effect has not yet been mechanistically described in a simple compartmental model depending only on the timing of the vaccination campaign.

Finally, we are interested in investigating whether the rebound effect seen in simpler models persists when examining a specific disease and a more advanced model. For this purpose, we expand upon the model by Arenas et al.⁹, originally developed for COVID-19, to include vaccination campaigns, focusing on the effective management of vaccination timing to minimize the second peak of infections and reduce cumulative hospitalizations and deaths. The age-stratified nature of this model also allows us to introduce an additional variable alongside the timing and duration: the age priority. Using data from a recent wave of an infectious disease, our findings provide key insights for optimizing vaccination distribution, aiming to mitigate the effects of future epidemic waves and inform public health policies.

In the following, we outline our findings about what separates a good vaccination strategy from a bad one, and in particular on how to get to an optimal one. This section is divided into four parts: The first part introduces the SIRS model with vaccination and its basic features. The second part focuses on the rebound effect (i.e., the unexpected growth of cases after a vaccination campaign) and gives a mechanistic explanation on why it happens. The third part outlines the importance of timing and duration of a campaign to achieve an optimal result. Lastly, the fourth part shifts to a more complex epidemiological model and deals with the problem of identifying the Pareto-optimal solutions that simultaneously minimize both the second wave of cases/hospitalizations and overall deaths.

The model

Let us consider a continuous-time mean-field SIRS model with vaccination, a compartmental model that divides the population into three groups: susceptible (S), infected (I) and recovered (R) individuals. In this model, susceptible individuals can become infected through contact with an infected person at a rate of β. Alternatively, they may gain immunity to the disease through vaccination, occurring at a rate of α(t), and consequently move into the recovered category. Infected individuals transition to the recovered compartment at a rate μ, while recovered individuals lose their immunity at a rate δ, becoming susceptible again. Note that the vaccination process in this model can be dynamically managed; it can be turned on and off at any moment according to the employed vaccination strategy. The dynamics of the model is sketched in Fig. 1.

A graphical representation of the SIRS compartmental model with reinfections and vaccination described in Eqs. (1)–(3).

Full size image

The equations that describe the model are:

Here, variables S(t), I(t) and R(t) represent the fraction of individuals in their corresponding compartments, which fulfill S(t) + I(t) + R(t) = 1; we also write the vaccination rate depending explicitly on time t, while the rest of the parameters are constant. We begin by studying the behavior of the system of Eqs. (1)–(3) when α(t) = α, to understand what different levels of vaccination have on the epidemic in the long term. We opted for a constant vaccination rate because, in many real-world scenarios, the pace of vaccination rollout is constrained by logistical factors and is independent of the number of susceptible individuals. First of all, we recall that, if β > μ, the stationary states S^*, I^*, R^* of the standard SIRS model with no vaccination (α = 0) are⁴²:

The basic reproduction number of this model is ({R}_{0}=frac{beta }{mu }) and does not change once we add a non-zero vaccination term. In the case β < μ the system rapidly reaches the “disease-free equilibrium”, where S = 1 and I = R = 0, which we consider uninteresting for the scope of this paper. Therefore, from now on we will always assume that β > μ.

Depending on the rate of vaccination α, we can identify three regimes: the coexistence regime, the eradication regime, and the total immunity regime, illustrated in Fig. 2. In the first regime, which occurs when the vaccination rate α ∈ [0, α_er], all compartments are present and the qualitative behavior of the system remains the same as the standard SIR but with stationary values ({S}_{{rm{co}}}^{* },{I}_{{rm{co}}}^{* },{R}_{{rm{co}}}^{* }):

The critical vaccination rate for eradication α_er = (1 − 1/R₀)δ is defined as the value of the vaccination rate for which the endemic fraction of infected individuals I^* goes to zero. Thus, the eradication regime, where no infective individuals are present, happens when α ∈ [α_er, α_ti], where the stationary values become

The critical vaccination rate required to achieve total immunity, denoted α_ti = δ, is defined as the vaccination rate at which the fraction of susceptible individuals also vanishes. If α exceeds this threshold, the system transitions into the third regime–the total immunity regime–characterized by the presence of only recovered individuals. Finally, during all of our numerical simulations we never allow the fraction of susceptibles to drop lower than 0.005 due to vaccination, in order to avoid pathological behaviors with unrealistically small numbers for the fraction of individuals in any compartment.

Illustration of the different vaccination regimes based on the value of the vaccination rate α. Please note that the three regimes are shown with identical lengths in this diagram, though their actual extents are governed by the values of α_er and α_ti.

Full size image

The rebound effect of misaligned vaccination

In the previous section, we assumed that the rate of vaccination of the SIRS model is constant over time; however, this is far from a realistic scenario. In many real-world emergency situations, vaccine supplies are constrained and distributed within specific timeframes. Consequently, vaccination efforts extend over defined periods known as vaccination campaigns. These campaigns have a designated start date and follow a daily schedule for distributing vaccine doses.

Building upon this premise, to incorporate a vaccination campaign into the previous model, we propose to use a piecewise function for the vaccination rate α(t). This function is defined for the duration of the campaign, with a starting time denoted as t_start and an ending time denoted as t_stop, as follows:

Let us begin our analysis by simulating the previous epidemic model using two different values for the vaccination rate α, thereby placing the model in two distinct vaccination regimes (Throughout all our numerical simulations, we maintain the fraction of susceptibles above 0.005 following vaccination to prevent pathological behaviors associated with unrealistically small numbers in any compartment.). Concerning timing, we will let the model reach a stable equilibrium before initiating the corresponding vaccination campaign. We simulate the system of Eqs. (1)–(3) with a value of α = 0.003 (corresponding to the coexistence regime) and α = 0.008 (corresponding to the eradication regime), see Fig. 3.

Illustrative example of the effects of limited-time vaccination in two different vaccination regimes. On the left, the vaccination rate α falls in the coexistence regime (α = 0.003), while on the right it is in the eradication regime (α = 0.008). The shaded area denotes the time period when the vaccination strategy is active. In both regimes, when the vaccination stops, a resurgence wave is observed, but with a different magnitude: in the coexistence regime, the wave is relatively modest, whereas in the second scenario, it is more pronounced, comparable in size to the initial wave. This observation suggests that a higher intensity of vaccination efforts may paradoxically exacerbate the rebound effect. Simulation with β = 1, μ = 0.5, δ = 0.01 and a vaccination campaign taking place between t = 400 and t = 600. The initial conditions are I(0) = 0.01 and S(0) = 0.99.

Full size image

In the figure, we observe that initially, without any vaccination strategy in place, the system experiences fluctuations until reaching a stable equilibrium. Following this, the initiation of the vaccination campaign leads to a reduction in the number of infected individuals. This reduction also involves the susceptibles if the strategy falls under the eradication regime.

Then, after ending the vaccination campaigns, we observe the anticipated “rebound effect”, characterized by a significant peak in the number of infected individuals, which is particularly evident in the eradication regime. This phenomenon hints that vaccination, despite its benefits, can have complex and sometimes unexpected effects on disease spread dynamics, and deserves a deeper investigation into its origins and underlying mechanisms.

This counter-intuitive phenomenon, like many aspects of epidemic modeling, can be understood through the availability of susceptible individuals. A disease can only spread if there are susceptible individuals to fuel its transmission. In this light, the efficacy of a vaccination campaign lies in its ability to transfer people directly from the susceptible category to the recovered one, thereby reducing the pool of individuals available for contagion. However, once the campaign concludes and assuming vaccine-induced immunity wanes over time, vaccinated individuals gradually revert to susceptibility. If the campaign was highly effective, a significant portion of the population might become susceptible simultaneously, lacking natural immunity from recent infections and becoming available for new infections. This synchrony can lead to an increase in infection numbers. Analogous to pulling and releasing a spring, the vaccination campaign initially lowers infection rates, but their subsequent increase can lead to a sharp rise in cases.

It is important to note that the extent of this rebound varies significantly depending on the vaccination strategy employed. In the co-existence regime, the rebound is small and can simply be seen as the system settling from one equilibrium (with vaccines, Eq. (5)) to another (without vaccines, Eq. (4)). On the other hand, in the eradication and total-immunity regimes, the rebound can be very large (depending on the duration of the vaccination campaign).

To mathematically understand the mechanism behind the rebound in the latter case, we must return to the modeling approach. As we mentioned, the eradication regime starts at the value of α_er where the number of infected individuals goes to zero. This allows us to develop a perturbation theory framework that gives us an approximated solution for the unfolding of the epidemic just after the end of the vaccination campaigns. At zeroth order in this approximation the contagion term βSI and the recovery term μI vanish so the system of Eqs. (1)–(3) can be approximated by:

whose analytical solutions are governed by exponential functions. The full calculations at first order together with all the definitions can be found in the Methods section. The full solution for t ⩾ t_stop then reads:

where (widetilde{Delta S}(I(t),t),widetilde{Delta R}(I(t),t)) and I(t) are the first-order contributions, and S(t_stop), R(t_stop) and I(t_stop) are the value for each compartment right after the end of the vaccination process (the agreement between these equations and the numerical solution is shown in Fig. S1 of the Supplementary Material). It is important to note that these equations are applicable only for a certain period after t_stop. Indeed, as long as the number of infected I(t) stays low, the contributions of (widetilde{Delta S}(I(t),t)) and (widetilde{Delta R}(I(t),t)) also remain small, allowing for the zeroth order of Eqs. (10)–(11) to dominate and the susceptible compartment to approach its limit S → 1.

Despite their limitations, these equations provide valuable insights and allow us to draw a causal connection between the strength of the vaccination campaign and the height of the rebound. For instance, we know that the approximation stops working at the time t_p where the number of infected, Eq. (12), reaches its theoretical upper bound, i.e., I(t_p) = 1 − μ/β. A closed, although complicated, expression for t_p exists and thanks to that it is still possible to demonstrate that (frac{partial {t}_{p}}{partial S({t}_{{rm{stop}}})} ,<, 0) and (frac{partial {t}_{p}}{partial I({t}_{{rm{stop}}})} ,<, 0) (see proof in Section S1 of the Supplementary Material). This implies that the smaller the number of susceptible and infective individuals remaining after vaccination, the larger t_p will be; thus, the validity of the exponential solutions to Eqs. (10)–(11) extends over a longer time period. Since t_p is connected to the growing phase of the zeroth order susceptible in Eq. (10), a larger t_p also means a larger build-up in the susceptible compartment and therefore a larger peak of susceptibles ({S}_{max }). Finally, one can prove, albeit based on other simplifying assumptions (see Section S2 of the Supplementary Material), that a larger ({S}_{max }) results in a taller epidemic peak (i.e., the height of the first peak of infectives following vaccination). This, therefore, demonstrates that a robust vaccination campaign can cause the rebound peak to be taller.

Timing as a fundamental feature to avoid the rebound

In the preceding section, we exclusively examined scenarios where the vaccination campaign commenced after the standard SIRS model had reached its equilibrium. Furthermore, the campaign had a designated duration intended to allow the system enough time to stabilize into a new equilibrium state under the influence of vaccination before returning to its initial condition upon the campaign’s completion. However, vaccination campaigns frequently overlap with epidemic waves that are still far from equilibrium, and their durations may differ. This highlights the importance of identifying the optimal timing and duration of vaccination campaigns to reduce the rebound effect while the epidemic is still ongoing.

To address this challenge, we have developed a simulation framework to systematically evaluate a broad spectrum of vaccination strategies, spanning various starting times and durations. Our objective is to find the strategy that most effectively counters the rebound effect by specifically aiming to minimize the height of the next largest peak (which usually coincides with the second peak, but that is not always the case, as shown in Fig. S4 in the Supplementary Material.), which defines our optimal criterion in this context.

To achieve this, we initially conduct a simulation without any vaccination to establish a baseline that serves as a reference point for comparing all subsequent vaccination scenarios. In simulations that incorporate vaccination, we preserve the baseline scenario’s parameters and initial conditions. A vaccination campaign is characterized by a specific starting day and duration. Note that the duration of the campaign is inversely proportional to the rate of vaccination α. We perform a separate simulation for each potential start day, ranging from day 1 to the day when the baseline scenario reaches its second peak (indicated by the second vertical line in Fig. 4(a)).

a Illustrates the proportion of infectious individuals over time in the baseline scenario, i.e., without vaccination. The vertical dashed lines indicate the peaks of the epidemic and act as reference points in the subsequent panels. b The variation in the height of the subsequent largest peak relative to the baseline (vertical axis) due to the most effective vaccination strategy initiated on each day (horizontal axis). Here, the term “best” strategy refers to the approach that results in the greatest reduction in the height of the subsequent largest peak. The color of each bar indicates the vaccination rate of that strategy. c The same quantities as (b) but, instead of using a constant vaccination rate α, we have used a vaccination rate α(t, S) ∝ Θ(S − S^*), i.e., a vaccination that is only active when the susceptible compartment exceeds its stationary value. The horizontal black dashed lines in panels b) and c) correspond to the maximum possible reduction, obtained when the second peak is just as tall as the equilibrium value. In this case, the strategy that comes closer to achieve the optimal value is the one in panel b) that starts on day 89 with speed α = 0.002.

Full size image

Additionally, we vary the duration of the vaccination campaign, exploring periods between 30 and 150 days. In all simulations, the total number of vaccines is fixed, equal to 30% of the population. This quantity ensures that, under the fastest rollout scenario of 30 days, the daily vaccination rate reaches 1% of the population. We select this rate as the maximum feasible daily coverage in a realistic scenario. The decision to use a fixed total quantity of vaccines was driven by our aim to optimize the use of a resource we consider to be finite.

Our goal is to identify the strategy that maximizes the benefits derived from the limited number of vaccines available. This approach is based on the understanding that having access to additional vaccines would be advantageous, as a greater quantity of vaccines generally leads to more favorable outcomes. This assumption is corroborated by the data presented in Fig. S2 in the Supplementary Material, which illustrates the inverse relationship between the total number of vaccines and the cumulative number of cases. However, it is critical to recognize that merely increasing the vaccine supply does not unequivocally resolve all challenges. Specifically, as the vaccine allocation escalates, the rebound effect may become more pronounced, as evidenced in Fig. S3 of the Supplementary Material. This observation highlights the need to develop strategies that not only make the most of the existing vaccine supply but also minimize the risk of a rebound effect. By doing so, we can improve the overall effectiveness of vaccination campaigns in controlling epidemic outbreaks.

Our findings are illustrated in Fig. 4b, where we present the outcomes of vaccination campaigns that start on different days. Each bar in these plots corresponds to a specific start day for the vaccination, with the bar’s color indicating the campaign’s duration and the height indicating the relative reduction or increase in the size of the subsequent peak compared to the baseline scenario. A bar extending upwards signifies an increase in the next peak, indicating a rebound effect, whereas a downward extension suggests a reduction in peak size, indicating no rebound effect. It is important to note that although multiple scenarios with varying durations were analyzed for each start day, the plots only display the outcome of the optimal strategy, defined as the one resulting in the smallest subsequent peak.

From our analysis, we find that if the vaccination roll-out starts around the time of the first peak, even the most successful vaccination strategies will suffer from the rebound effect. This can be seen in Fig. 4b, where we see that those simulations that consider a starting time around the first peak experience the rebound effect, indicated by bars with positive heights, signifying an increase in the epidemic’s peak magnitude. In turn, strategies that start after the first peak have a better chance of avoiding any resurgence, particularly when vaccines are administered over an extended period (i.e., with vaccination within the coexistence regime). The optimal strategy observed occurs just a few days following the first peak and is the only approach that achieves a complete flattening of the infectious curve. This is evident from its proximity to the dashed line in Fig. 4b, which indicates a reduction making the second peak equivalent to the equilibrium value. Outside of this optimal strategy, determining a straightforward rule for how soon before the second peak to start and the ideal length of the campaign proves challenging.

To test the idea introduced in the previous section, i.e., that the rebound effect is caused by an excessive departure of the susceptible and infected from their equilibrium values S^* and I^*, we introduce a new kind of vaccination campaign. Previously, our model used constant, nonzero vaccination rates α(t) during a specific time frame, with rates dropping to zero outside of this interval. To avoid pulling the number of susceptible too far away from their equilibrium and thus cause the rebound effect, we test a new vaccination strategy. Here we introduce a new formula for the vaccination rate α(t, S) ∝ Θ(S − S^*), where Θ indicates the Heaviside function. This functional form ensures that the vaccination will only be active when the susceptibles are above their equilibrium number. The vaccination finally stops at the moment in time when all the doses at our disposal (30% of the population in this case) run out. Please note that this strategy is characterized by intermittent activation, as the susceptible population’s numbers oscillate around the equilibrium point. This fluctuation causes the vaccination campaign to switch on and off, until it finally ceases when the supply of vaccines is fully depleted. It can be seen in Fig. 4(c) that this “fine-tuned” strategy works wonderfully, resulting in a consistent relative reduction across a wide range of starting times and almost always outperforming the constant strategy. The only exception to this rule are the truly optimal solutions leading to the total flattening of the second peak, that can only be obtained with the constant strategy.

Adding complexity: rebound effect in an age-stratified metapopulation model

In the preceding section, we explored how the rebound effect is influenced by two factors: the vaccination start day and its duration. To verify the rebound effect’s consistency in more complex scenarios, we transition to the model by Arenas et al.⁹, which was developed to simulate the spread of COVID-19 during its initial outbreak. A brief overview of this model’s key components is provided here, with a comprehensive explanation available in the Methods section and further details in Sections S3 and S4 of the Supplementary Material, while a discussion of the rebound effect in this model can be found in Fig. S9.

Building upon the conventional SIRS framework, this model integrates further phases of the infection cycle, notably the exposed (E) and asymptomatic (A) stages, in addition to hospitalization (H) and death (D). Additional compartments were introduced to accurately represent the timeframes leading to hospitalization and death. An illustration of the compartmental model is provided in Fig. 5. The model employs a metapopulation approach, incorporating the movement-interaction-return mechanism⁴³. This implies that the population is distributed across various patches, with individuals moving to adjacent areas and then returning to their original locations, without permanent migration. Moreover, the model incorporates age stratification, categorizing individuals as young, adult, or old, to account for the infection’s outcome varying significantly with age. However, the original model in⁹ does not account for vaccination and waning immunity. Therefore, we adapted the model to include these aspects, as shown in Fig. 5.

The acronyms correspond to susceptible (({S}_{v}^{g})), exposed (({E}_{v}^{g})), asymptomatic infectious (({A}_{v}^{g})), symptomatic infectious (({I}_{v}^{g})), pre-hospitalized in ICU (({P}_{Hv}^{g})), pre-deceased (({P}_{Dv}^{g})), in ICU before recovery (({H}_{Rv}^{g})), in ICU before death (({H}_{Dv}^{g})), deceased (({D}_{v}^{g})), and recovered (({R}_{v}^{g})), where g denotes the age stratum of all compartments and v the vaccination status. The arrows indicate the transition probabilities.

Full size image

As with the SIRS model, we start by running a baseline simulation without vaccination (see Fig. S5 in the Supplementary Material). Then, we run an extensive number of simulations where we systematically explored various combinations of durations, starting days, and age priorities, all while maintaining a constant total number of vaccines. We vary the starting day of the vaccination, considering all days between the beginning of the simulation and the midpoint between the first and second peak. For the duration of the vaccination, we use 30, 60, 90, and 120 days and we keep the total number of doses fixed (so that all strategies deliver the same amount of vaccines). To explore the age group priority, we used all the possible combinations of strictly positive values that are multiples of 10% and that sum up to 100%. The resulting strategies are then categorized as young, adult, old, or mixed based on the priority levels of each age group. For instance, a strategy labeled “young” means that the priority given to the young group is higher than any of the priorities given to the two remaining groups. The mixed strategies are those that do not fulfill the previous condition.

This time, to determine the optimal strategy, we employ two metrics: the relative variation of the height of the second peak, as in the previous section, and the relative variation of cumulative deaths compared to the baseline simulation. Our objective is to minimize both metrics simultaneously by identifying Pareto-optimal solutions. These are solutions where improvements in one metric would lead to worsening the other. The latter metric, introduced due to the inclusion of age structure in our model, accounts for varying risk levels across different age groups, making it crucial for evaluating the impact of prioritizing one age group over another in vaccination plans. However, this metric is significantly affected by the total runtime of the simulation. To circumvent arbitrary decision-making, we conclude the simulation after the epidemic’s transient phase.

The results are shown in Fig. 6. It can be seen that most of the simulations (each represented by a point) fall on positive values of the difference in cases, meaning that they perform poorly or are counterproductive, so we restrict our attention to those solutions that are able to improve the baseline scenario. Among those, the optimal ones are the ones that satisfy the aforementioned Pareto condition. Such points lay at the left-most side of the scatter plot and taken together lay on a line called the Pareto front, which we highlight in Fig. 7, where we also show a detailed description of the features of each of these optimal solutions. From this figure, it emerges that there are two kinds of optimal strategies, each distinctly different from the other. The first type, that we mark as Type 1, contains all those strategies that are better at reducing the number of deaths while overlooking the height of the next largest peak. In order to achieve this the vaccination has to start before the first peak is reached, it must be fast and must prioritize old people. On the other hand, the second type performs better at lowering the height of the next largest peak while allowing for a larger number of deaths. In this second case the vaccination has to start well after the first peak, must have a longer duration, and should also prioritize old people.

In this scatter plot, each point corresponds to a unique simulation based on a distinct vaccination strategy. The horizontal axis displays the relative variation in the number of cases at the next largest peak relative to the baseline case, while the vertical axis represents the relative variation in the number of cumulative deaths, both with respect to the baseline. Our objective is to minimize these two metrics in order to identify strategies that outperform the baseline scenario. Each plot point conveys three distinct variables: color, shape, and size. The color signifies the strategy’s implementation time (t_start), with brighter shades indicating later starting time. The point’s shape indicates the primary age group targeted by the strategy, and the point’s size reflects the campaign’s duration, where larger sizes denote longer duration. It is important to note that the total number of vaccines administered remains constant, irrespective of the campaign’s duration. As we can see, many of the analyzed strategies yield negative outcomes. Points on the positive side of the horizontal axis indicate an increase in the number of cases during the second peak. However, it is important to note that none of the simulations results in an increase in the number of deaths. Finally, the black solid line marks the Pareto front, which indicates the optimal strategies.

Full size image

This table shows the characteristics of Pareto-optimal strategies alongside a scatter plot that includes only those data points belonging to the Pareto front. The columns of the table provide information on the duration and starting time of each strategy, the priority vector that determines the vaccines distribution by age group (ordered as [Young, Adult, Old]), as well as the age stratum that is most prioritized for each strategy. The asterisk indicates simulations that are optimal for reducing both the number of cases and the number of hospitalizations (see Fig. S7). It is possible to identify two different sets of strategies, highlighted by yellow and violet ovals and labeled Type 1 and Type 2. In Type 1 strategies, the primary objective is to minimize the number of deaths, indicated by the lower values on the vertical axis. This is best achieved by vaccinating early, before the first peak, and using short-duration campaigns. Conversely, Type 2 strategies aim to reduce the number of cases, indicated by the lower values on the horizontal axis. To achieve this goal, it is optimal to vaccinate later, before the second peak, and use longer duration campaigns. Both types suggest to strongly prioritize elderly people over everyone else.

Full size image

An analogous analysis optimizing the number of hospitalizations instead of cases can be found in Figs. S6 and S7. This analysis demonstrates similar behavior, as the strategies are divided into the same two groups. However, the transition between Type 1 and Type 2 strategies is much smoother and does not present the characteristic right angle that we see in Fig. 6. This can be attributed to the fact that hospitalizations’ peak and cumulative deaths are closely connected and therefore the strategies to optimally reduce the first are bound to have a deep impact on the second as well. Furthermore, in the case of daily hospitalization, the rebound effect disappears almost entirely, due to the fact that the vaccines are modeled to give long-lasting protection against the risk of hospitalization. Different strategies still produce widely heterogeneous outcomes, but in almost all of them, the second peak is lower than the one in the no-vaccines scenario. It is interesting to note that strategies that fall on the Pareto front in Fig. 6 are not the same strategies that constitute the Pareto front in Fig. S6. In other words, the strategies that minimize cases and deaths are not the same than those that minimize hospitalizations and deaths, although some overlap can be found (see the asterisks in Fig. 7). This means that, perhaps unsurprisingly, cases and hospitalizations are not interchangeable as metrics of the outcome of an epidemic, especially when age-specific risk factors are included in a model.

From our analysis, it emerges that the outcome of a vaccination campaign, assuming vaccine waning and immunity decay, is more strongly influenced by the timing of its initiation than by any other factor. This finding is in agreement with previous studies^29,40. Even when we have a large number of vaccines at our disposal, a poor choice in the vaccination timing could lead to counterproductive effects. This happens because two forces are at play: on the one hand, vaccines have a high impact in the short term, helping to reduce the number of cases and fatalities; on the other hand, if the vaccination reduces the infections well below their equilibrium value, that can provoke a rise in the susceptible population and a subsequent increase in the severity of subsequent waves. This rebound effect is the main obstacle to an efficient vaccination strategy and it is therefore a phenomenon worth studying. For this reason in this work we investigated this effect by giving it a mechanistic explanation based on an approximation of the SIRS model. In this regard, all the simulations hint to the fact that the best time to start vaccinating is before the cases start rising, i.e., either before or after the epidemic peak. Given the non-linear origin of the effect, it could not be eliminated even by using vaccination boosters, since every subsequent campaign would obey the same trade-off of the one that came before. The best one can do in order to avoid the rebound completely is using a fined-tuned strategy specifically created to keep the number of susceptible people close to their equilibrium value. Finally, in a more complex scenario, different strategies are available depending on the aim of the campaign and of time at our disposal: to minimize deaths, the priority should be given to fast and early strategies; if instead our goal is to lower the peak prevalence, slower and later strategies are the better solution.

Unfortunately, the rebound effect we identified and its dependency on timing is hard to observe empirically, even in principle. The main problem to overcome would be the lack of a counterfactual scenario (similar to our baseline simulation) to compare the outcome of the any vaccination strategy against. The ideal setup to deal with this is a natural experiment in which two countries, as similar as possible to each other in health standards, age distribution and vaccination statuses, react to the same epidemic wave with differently timed vaccination campaigns; any difference in the resulting outcome could then be interpret in light of such a difference in timing. Situations like these are hard to achieve in a real-world scenario, first and foremost because countries are complex systems, containing a multitude of potentially alternative explanatory variables to the one we are trying to test. A more indirect way to demonstrate the potential existence of a rebound phenomenon would be to show that the hidden mechanism behind it (i.e., the slow build-up of the susceptible population) does in fact produce long periods of calm followed by large outbreaks, even if said mechanism has a different origin than the one discussed in this paper. An example of such an empirical confirmation comes from all the instances of re-emergences of diseases in highly vaccinated populations: measles in Japan (2007)⁴⁴, France (2011)⁴⁵ and various countries in Southern Africa (2009-2010)⁴⁶. In all these examples the effect is not as pronounced as the one we have shown here, and the reason for that is probably that in none of these cases the vaccination coverage ever drops to zero.

Finally, we recognize that from a public health perspective, the suggestion of delaying the beginning of a vaccination campaign right when cases are ramping up might be hard to apply, given the strong pressure that policymakers would likely be facing at such a moment. However, it is worth to point out that our study offers some insights into how to handle that situation as well, and it can help to decide whether a slow or a fast roll-out is preferable depending on where we are, relative to the epidemic peak. This task can be performed with the forecasting methods available to us (the discussion of which is beyond the scope of this paper) whose uncertainties are usually small enough to provide us with a satisfactory knowledge of the state of the epidemic. Finally, our work not only offers insights on how to react during the rapid unfolding of events, but also within the context of proactive and preventive planning. It is in fact clear from our analysis that, no matter how small the number of infected individuals is, the susceptible population should not be left free to grow uncontrollably, but should be kept under control with timely interventions.

No-infected approximation of the SIRS model

As already mentioned, once a vaccination campaign in the eradication regime comes to an end, it is possible to find an approximate solution for the evolution of the system. That happens because in that timeframe the number of infectious individuals is so small that the system starts evolving as if there were no infections nor recoveries and the whole dynamics was dominatated by the waning immunity. That means that the system is evolving according to the no-infected equations in (8) and (9), whose solution are

where S₀(t) + R₀(t) = 1, since we explicitly assume that I₀(t) = 0. The subscript 0 has been introduced to differentiate the solution of this simplified system from those of the full SIRS model. Furthermore, we can introduce the quantities ΔS(t) ≡ S(t) − S₀(t) and ΔR(t) ≡ R(t) − R₀(t), which are the differences between the solutions of the full SIRS model and the solutions of the infected-less system. Both of these terms are of order O(I), therefore as long as I(t) stays small they will also remain small. Finally it should be noted that since we know that S(t) + I(t) + R(t) = 1 and S₀(t) + R₀(t) = 1 it follows that ΔS(t) + I(t) + ΔR(t) = 0.

By rewriting Eq. (2) using the definition of ΔS(t) and dropping the term which is of order O(I²), we can obtain an approximated equation for the infected compartment I(t)

The solution to this equation is found through the method of separation of variables and is given in Eq. (12). By repeating the same process with Eq. (3) we can then use the newly found I(t) to find a first order estimation of the function ΔR(t), which we call (widetilde{Delta R}(I(t),t))

This equation can be solved by recognizing that it is a linear equation where I(t) plays the role of a forcing term. The solution is:

Finally, a first order estimate of (widetilde{Delta S}(I(t),t)) can be found by exploiting the relation ΔS(t) + I(t) + ΔR(t) = 0. This, in turn, allows us to numerically estimate both the magnitude of the maximum value reached by the susceptible compartment ({S}_{max }) as well as the time at which it is reached ({t}_{max }) by looking for the point in which the condition (frac{{rm{d}}}{{rm{d}}t}({S}_{0}(t)+widetilde{Delta S}(I(t),t))=0) is satisfied, see the Supplementary Material.

Details of the COVID-19 model

As already mentioned, the COVID-19 model that we used in this work is an extension of the model introduced by Arenas et al.⁹ and, therefore, similarly to that model, it has a compartmental dynamics that includes susceptible, exposed, asymptomatic, infected, hospitalized, recovered and deceased, together with some additional compartments to regulate the latency periods. All of these groups are arranged as shown in Fig. 5. Similarly to the original model, the current version also takes place on a metapopulation network where the mobility is recurrent and modeled through the MIR (Movement-Interaction-Return) framework⁴³.

However, differently from the original, the version presented here introduces some features (such as vaccinations and reinfections) that make this model more useful to describe the late stages of the pandemic. In particular, vaccinations have been designed to realistically resemble the mass vaccination campaigns that started in 2021. First of all, unlike all the other transition terms (see Eqs. (S.28)–(S39) in the Supplementary Material), the rate of vaccination does not depend on the number of unvaccinated individuals, but rather involves a fixed number of vaccines that are distributed among the population according to two rules. The first rule prioritizes patches with a high population density, while the second rule assigns different levels of priority to different age groups according to ε^g. Using these rules, we obtain ({epsilon }_{i}^{g}(t)), which represents the absolute number of vaccine doses per age g and location i, as follows:

where the normalization is chosen so that the number of doses per day is fixed. Additionally, we introduced a mechanism to limit and redistribute the doses whenever the number of vaccinated people outnumbers the total number of people in each patch.

From the modeling perspective, being in a vaccinated compartment signifies reduced probabilities of initial infection or, if already infected, of transmitting the disease, being hospitalized, or dying, compared to those in the non-vaccinated compartment. After a certain latency period, individuals lose their vaccine-acquired immunity. However, they transition to a “post-vaccinated” status rather than reverting to full susceptibility. In this status, their risk of infection and transmission aligns with that of unvaccinated individuals, but their chances of hospitalization or death remain as low as those who are vaccinated.

Several technical details of our simulations were shared across all simulated scenarios. The specific parameter values used are compiled in Tables S1 to S4 in the Supplementary Material, along with the rationale behind their selection. Some parameters were based on our own assumptions or those of previous studies⁹, while others were estimated using scientific reports^47,48 and through a calibration process. The initial conditions for our simulations included 10000 individuals in each of the exposed, asymptomatic, and infected compartments. The rest of the population was assumed to be in the susceptible non-vaccinated compartment. The simulations were run for a duration of 600 days, which was chosen to ensure that both the first two infectious peaks were included in the analysis, independently of the vaccination strategy.