Vaccination strategies in structured populations under partial immunity and reinfection

Optimal protocols of vaccine administration to minimize the effects of infectious diseases depend on a number of variables that admit different degrees of control. Examples include the characteristics of the disease and how it impacts on different groups of individuals as a function of sex, age or socioeconomic status, its transmission mode, or the demographic structure of the affected population. Here we introduce a compartmental model of infection propagation with vaccination and reinfection and analyze the effect that variations on the rates of these two processes have on the progression of the disease and on the number of fatalities. The population is split into two groups to highlight the overall effects on disease caused by different relationships between vaccine administration and various demographic structures. As a practical example, we study COVID-19 dynamics in various countries using real demographic data. The model can be easily applied to any other disease transmitted through direct interaction between infected and susceptible individuals, and any demographic structure, through a suitable estimation of parameter values. Two main conclusions stand out. First, the higher the fraction of reinfected individuals, the higher the likelihood that the disease becomes quasi-endemic. Second, optimal vaccine roll-out depends on demographic structure and disease fatality, so there is no unique vaccination protocol, valid for all countries, that minimizes the effects of a specific disease. Simulations of the general model can be carried out at this interactive webpage Atienza (2021 S2iyrd model simulator).

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
vaccine roll-out depends on demographic structure and disease fatality, so there is no unique vaccination protocol, valid for all countries, that minimizes the effects of a specific disease. Simulations of the general model can be carried out at this interactive webpage Atienza (2021 S 2 iyrd model simulator).
Supplementary material for this article is available online Keywords: vaccinations, epidemic modeling, SIR model, optimal strategies, population statistics, nonlinear dynamics (Some figures may appear in colour only in the online journal)

Introduction
Infectious diseases have profoundly shaped human habits, and societies at large [1]. Their tremendous effects on morbidity and mortality could only be counteracted with the advent of massive vaccination campaigns, which stand out as the most effective medical interventions ever [2]. This is so despite the fact that only smallpox has been eradicated (also rinderpest in animal hosts), while multiple other communicable diseases persist at low incidence or with significantly milder effects in infected individuals-thanks to vaccination. Eradication of pathogens is difficult if they have alternative species as reservoirs, in the absence of a vaccine, or if the disease is geographically spread or has a complicated diagnosis.
Nowadays, many infectious diseases remain endemic even if a vaccine is available and its administration is widespread. Vaccinated individuals can become infected because vaccines only confer partial protection [3]; if immunity wanes, both vaccinated and previously infected individuals can get the infection again [4]. Around 1% of infected individuals undergo reinfection by SARS-CoV-2 [5], a percentage that is higher in the case of other respiratory diseases, such as influenza, where the high variability of viral strains and continued reinfections during each flu-season have made it endemic [6]. Indeed, the emergence of mutant forms of pathogens that escape immune attack and thus limit the protective effects of past infections or vaccination is the rule [7], especially when prevalence is high [8].
The eventual fate of an emergent epidemic is uncertain as it depends on the type of immunity people acquire through infection or vaccination and on how the virus evolves [9]. The effects of disease in populations along typically lengthy transients between the emergence of a contagious disease and its dynamical stabilization depend on features of the population and of the disease, but also on controlled interventions such as non-pharmaceutical measures and, critically, vaccine roll-out [10].
It is hard to quantify all possible effects of vaccination [11], since it is ideally aimed at fulfilling several dissimilar goals, among others to prevent the disease and reduce its severity and mortality and to minimize the impact of the disease on the health care system and the economy [12]. Actually, vaccination strategies often focus on one or another priority: for instance, on vaccinating first the most vulnerable groups, or on starting with the individuals with the highest spreading potential, with the hope that the most vulnerable group will be indirectly protected [13,14]. Vaccination strategies cannot be unique and universal [15], since their effects vary as a function of the transmission mode and severity of the specific disease, of population habits (which determine contacts and therefore contagion probability), of demographic profiles and, last but not least, on the availability of vaccines [16].
In this contribution, we address the effect of different vaccination protocols in structured populations by means of a compartmental model, and assuming a non-negligible fraction of reinfections and infections of vaccinated individuals. As a measure of the effects of vaccination, we use the reduction in the number of deaths after one year, as compared with a situation of no vaccination. The model is deliberately simple to focus on the overall effects of two fixed empirical variables (the infection fatality ratio (IFR), of the disease and the demographic structure of the population) and one controllable variable: the vaccination rate. In the first part of the work, we explore the effects of these variables in synthetic populations structured into two groups where the first one has a higher number of intra-group contacts (and thus is the group that best spreads the disease) and the second one contains more vulnerable individuals (with a higher IFR). In the second part of the work, we apply the model to COVID-19 and several countries. In each particular case studied, we use available data on demography, contact matrices between age-groups [17,18], and independently evaluated, age-dependent IFR. Finally, we have built a public webpage [19] where multiple additional demographic profiles can be chosen and model parameters arbitrarily changed, also with the goal of allowing quantitative comparison with other models in the literature [13,20].

Models
Since the introduction in 1927 of the Susceptible-Infected-Removed (SIR) model [21], compartmental models have been broadly used to simulate the dynamics of epidemic propagation [22]. Major strengths of such models are their ease of use and interpretation, as well as their ability to reproduce general features of the epidemic process. Multiple extensions have been introduced through the years [23], including the addition of classes for the various stages of a disease (asymptomatic, diagnosed, hospitalized) or for individuals with variable degrees of mobility (quarantined), among others. The heterogeneous nature of populations has also been considered, thus allowing for more accurate representations of diseases that differentially affect children and adults [24], for example. Recently, compartmental models have been used to describe the current spread of SARS-CoV-2 and the effect of contention measures [25] or spatial propagation [26]. These models have been particularly useful to compare different scenarios [27] and to extract major trends under changes of parameters.

SIYRD: A compartmental model with reinfections and vaccination
The basic SIYRD model is a compartmental model with five different classes: Susceptible (S), Infected (I), Reinfected (Y), Recovered (R) and Dead (D) individuals. The recovered class merges individuals from two origins: those which were infected and have overcome the disease, thus bearing partial immunity to new infections, and susceptible individuals which have been vaccinated. Our main goal will be to study the effects of population structure and vaccination rate in the control of an epidemic. Vaccine effectiveness will play a quantitative role that is disease-and vaccine-dependent, and that is difficult to quantitatively compare with the immunity conferred by infection [28,29]. Such complexity is out of the scope of this work. As a first approximation, we assume that the degree of immunity acquired by vaccinated individuals is of the same order as that of recovered individuals after natural infection, and for simplicity make both parameters equal. Figure 1 summarizes the transitions allowed between compartments. Note that only vaccination of susceptible individuals is considered, due to the assumed equivalence between vaccination and disease overcoming.
For simplicity, we neglect waning immunity effects: Recovered individuals continue in that state unless newly infected. This situation models partial immunity acquired from infection or vaccination that lasts for the timescale of the study, which can span well over the period of . Infected individuals may recover (R, rate r) or die (D), regardless of whether they have primary infections (I, rate µ I ) or reinfections (Y, rate µ Y ). Recovered individuals can become reinfected (by class Y at rate β RY or by class I at rate β RI ), and either recover again (to class R, rate r) or die (to class D, rate µ Y ). Green dashed lines link classes whose interaction triggers transitions. (B) As in (A), for the model with two groups. Note that transitions occur only within classes of the same group, while interactions intra-and inter-groups (green dashed lines) couple the dynamics. These interactions are weighted through contact matrices and the fraction of individuals in each group, which consider the demographic structure of the population. an infection wave, given that an infection produces protection that remains high after a period of 40 weeks [30]. Our choice differs from a classic SIRS model, where recovered individuals revert to susceptible during the timescale of the study.
The equations that define the epidemic model arė Parameters are rescaled in such a way that the time unit of our simulations is one day. Vaccination is implemented in this model through a parameter v that represents the fraction of population vaccinated per time unit, i.e. per day. This rate is multiplied by a function Θ(S, θ) that takes into account its progressive slowdown and eventual halt when a fraction θ of individuals has been vaccinated. The specific functional form of Θ(S, θ) does not affect the overall dynamical properties of the system (see supplementary information).
Following widespread practice in compartmental models, we assume exponential transition rates, although more complex viral dynamics have been observed [31]. Non-exponential dynamics would be important to quantitatively predict dynamics. However, our motivation is to understand, therefore we sacrifice a quantitative for a qualitative description in order to simplify the model as much as possible.

Model parameters.
The separation between primary infections (I) and reinfections (Y) entails four different infection rates (see figure 1): β SI and β SY correspond to the infection rates of susceptible individuals by primary-infected and re-infected individuals, respectively; β RI and β RY are the equivalent infection rates for recovered individuals. Although the precise values of some of these infection rates might be difficult to estimate in general, we consider disease types where they are subject to certain relations that hold, on average, over the population. For example, reinfected individuals have a lower infective capacity in comparison to primary-infected individuals, both towards susceptible and recovered individuals. In other words, reinfection rates are smaller than primary-infection rates, β SI ⩾ β SY and β RI ⩾ β RY .
Consistently, the likelihood that a susceptible individual becomes infected is larger than that of a recovered individual, since the latter bears at least partial immunity against the disease either due to a prior infection or to vaccination. This applies both to primary infections, β SI ⩾ β RI , and to reinfections, β SY ⩾ β RY . Finally, we assume that the ratio between the infection rates of primary-infected and reinfected individuals is independent of the state of the individual that can be potentially infected, The mortality rate of primary-infected individuals can be calculated using the IFR of the disease under consideration, defined as the ratio between the number of fatalities and the number of infections in a given (sub)population (see supplementary information). The mortality rate of reinfected individuals is assumed to be significantly lower, µ I ≫ µ Y , as a result of their partial immunity against the disease. Finally, for simplicity we assume that the recovery rate r of primary-infected and reinfected individuals is identical, neglecting possible small differences.
The quantitative estimation of r depends on the IFR and on the infectious period of the disease, see supplementary information. The vaccination rate v is a variable that can be explored in a range of values depending on vaccine supply. We will consider values up to 2%, representing the fraction of population vaccinated in one day. A summary of parameters in the SIYRD can be found in table 1.

Stratified SIYRD (S 2 IYRD)
In order to consider the demographic structure of populations and how a given disease affects different groups (e.g. depending on age or sex), we extend the SIYRD model to represent two different population groups (G 1 and G 2 ) per class, each containing N 1 and N 2 individuals, respectively. Variables and parameters characteristic of each group now get the corresponding subindex: this applies to vaccination rates v i , mortality rates µ Ii and µ Yi , and recovery rates r i . Disease transmission rates β SI , β SY , β RI and β RY are independent of the group to which they apply. To be consistent, the two groups are connected to each other through a contact matrix that reflects individual connectivity patterns (see supplementary information). To simplify notation, the simple, one-class model with reinfection and vaccination will be called SIYRD all through the text, while the model with two coupled groups will be called S 2 IYRD. Equations (7)-(11) describe the dynamics of group G 1 in the S 2 IYRD model, This model has a second set of analogous equations, coupled to equations (7)-(11), obtained by simply interchanging subindexes, 1 ↔ 2.
The most important modification of the model comes from the matrix of contacts M, whose elements M ij represent the number of contacts an individual of group i has with individuals of group j, thus weighing the effect of intra-and inter-group contacts in contagion rates (see supplementary information).

S 2 IYRD model parameters.
Though mortality and recovery rates are in principle subject to relationships analogous to those described for the SIYRD model, they can be also independently estimated for each group as characteristics of the disease under consideration. Now, the values of the elements M ij synthesize the contact structure of the affected population and, in general, can be empirically obtained using actual demographic structure and surveys of contacts between different groups in a variety of daily situations (see supplementary information). Finally, the vaccination rate can be varied to represent different administration protocols, in order to test how they affect disease progression. We will consider three situations in which priority is given to G 1 , to G 2 , or both groups are simultaneously vaccinated. The result of each protocol will be compared with the baseline of no vaccination. Specifically, the prioritized population group is vaccinated at a rate v until the fraction of susceptible individuals reaches (1 − θ). When that happens, vaccination of the second, non-prioritized group begins, provided the previous threshold θ has not been yet reached through natural infections. Under simultaneous vaccination, both groups are vaccinated at rates proportional to their respective group sizes. If vaccination finishes in one of the groups first, the other group receives all doses. Note that, due to the vaccination prescription, where a fixed number of individuals (equal to the vaccine doses available) are vaccinated at each time, our simulations will be performed with finite populations of size N = 10 8 , and with the real population size in specific examples. To permit comparison of different results, however, we will use normalized variables In every specific scenario, three factors emerge as the main determinants of epidemic spread and thus of optimal vaccination protocols: the division of the population into groups and their contact habits (subsumed under the M ij terms), the etiology of the disease (IFR and recovery, death and infection rates), and the vaccination rates v i -which stand as the only variables that can be externally modified when the disease propagates freely or under constant nonpharmaceutical measures.
To account for the possibility of different timescales in the dynamics due to the different mechanisms present in the model, we have numerically integrated the equations using the variable-step, variable-order implicit solver ode15s from the MATLAB ODE suite [32]. This solver is appropriate for stiff problems with different timescales, and we have not observed unrealistic crossings of variables to negative values [33]. It also performs well for non-stiff problems, so it is a safe choice for compartmental models that tend to be well-behaved regarding numerical simulations. When vaccination is suppressed, a second steady state is also possible, (S * , 0, 0, R * , D * ), which is linearly unstable and therefore irrelevant from an epidemiological viewpoint. A neutrally stable pure endemic steady state (0, 0, Y * , R * , D * ) emerges if the mortality rate of reinfected individuals is set to zero, µ Y = 0. Further details on fixed points and their stability, together with figures illustrating the dynamics towards different final states, can be found in supplementary information. In principle, it is not possible to determine to which steady state the system will converge given the initial conditions and specific parameter values.

Long transients with a quasi-endemic disease.
Due to the introduction of reinfection, the model presents an interesting behavior in a broad range of parameters: the fraction of reinfected individuals Y and the total number of fatalities D experience a long, plateau-like transient before relaxing to their asymptotic values, see figure 2. This behavior is related to the existence of a stable endemic state for µ Y = 0. Indeed, the mechanistic origin of long transients relies on the loop that links R and Y classes, whose dynamical equations fully decouple from the rest in the limit when S = I = 0 and µ Y = 0, and present slow dynamics for small values of µ Y . In other words, these long transients are mostly controlled by the low mortality rate of secondary infections, which limits the flux of individuals out of the Y-R loop and delays the end of the epidemic. While, outside this regime, the asymptotic state is reached in about a year time, it may take over a century for the disease to disappear in the slow regime. This behavior is observed for a broad range of parameters in the relevant case when the dynamics of first infections occurs at rates several-fold faster than that of secondary infections. In practice, what the model yields in a few years time is a long plateau in epidemic incidence that can be interpreted as a quasi-endemic disease. At longer timescales, however, the model is no longer a valid description of the dynamics, since variables such as immune decay, mutations yielding new viral variants and, eventually, demographic changes, come into play. It is in these situations when diseases can become truly endemic, since the pool of susceptible individuals is replenished through one or another mechanism.

Dynamics of the S 2 IYRD model
In this section, we will discuss the dynamics of the S 2 IYRD model under the three possible vaccination strategies: priority to G 1 , to G 2 , or simultaneous vaccination. The advantage of either strategy will be evaluated under variable vaccination rates. Here, we do not make explicit the criterion to split the population into two groups, and explore instead the effect of these two groups having dissimilar sizes and being differently affected by the disease. To this end, we will study two representative values of the contact strength between groups, vary the fraction of population in either group and the IFR ratio, and illustrate how these affect the optimal vaccination strategy under increasing vaccination rates. Parameters can in principle be varied independently, but the closure relationship M 12 n 1 = M 21 n 2 (see supplementary information) has to be fulfilled for consistency. Also, we will restrict our explorations to cases where the group at higher risk is the slower spreader; otherwise, there is no conflict between vaccinating the most susceptible or the group that best spreads the disease, and the best strategy is therefore trivial. Without loss of generality, group G 1 will be the faster spreader, and group G 2 the one most affected by the disease. We will use as measure of the effectiveness of a strategy the asymptotic decrease in the number of fatalities once the epidemic halts (or, alternatively, at a sufficiently long time since vaccination started).

General dynamical features.
The behavior of the model is qualitatively robust for all scenarios explored, for different parameter values and contact matrices, as long as the initial reproductive number is slightly above the epidemic threshold (R 0 = β SI d I ≳ 1). Nevertheless, important quantitative differences arise (see examples in the next section). As time elapses, the number of susceptible individuals in both groups is reduced due to infection and vaccination. Simultaneously, there is an increase in the number of infected individuals and a subsequent increment of recovered and dead individuals. Furthermore, increases in the number of recovered individuals may cause a boost of reinfections under sufficiently high prevalence. Ultimately, when susceptible individuals become exhausted, primary infections decline and model dynamics become exclusively dependent on recovered and reinfected individuals. Infection rates and contact between groups mostly determine the speed of the epidemic, especially affecting how fast susceptible individuals move to other compartments. In turn, mortality and recovery rates play the largest effect in setting the final number of deaths. If the previous parameters are fixed, it is the vaccination rate that determines the temporal evolution of the epidemic and the optimal vaccination strategy.  Figure 3 summarizes the overall behavior and illustrates how variations in the vaccination rate and in the relative size of the groups affect optimal vaccination protocols.
Characteristic dynamics under slow and fast vaccination rates are represented in figure 3(A). If the vaccination rate is too slow, the dynamics are barely distinguishable from free propagation, and the overall reduction of the number of fatalities is very modest. For sufficiently high vaccination rates (around 0.5 and higher), the reduction of deaths becomes significant. Figure 3(B) represents the reduction obtained in number of deaths for each strategy as a function of the vaccination rate and the relative group size. At low vaccination rates, the advisable strategy is to vaccinate first group G 2 , the most vulnerable, though the optimal strategy changes to priority given to group G 1 for sufficiently high vaccination rates, the precise value of v at which the best strategy changes depending on the relative size of the groups (see also figure 3(C)). The selection of the optimal strategy is more critical at intermediate values of the vaccination rate and when the two groups have comparable sizes. At very high vaccination rates, however, vaccination strategies yield similar reductions and are weakly dependent on the relative group size. Figure 3(D) presents a quantitative summary of the effect of the best strategy. Cells to the left of the thick dark red line correspond to priority to G 2 as the advisable protocol, with a reduction in the number of deaths given for each vaccination rate and population fraction n 2 . To the right of the thick dark red line, the optimal strategy is priority vaccination of group G 1 . Reduction in the number of deaths for the optimal strategy. Panels stand for three different relative group sizes, as a function of the vaccination rate (above) and for three vaccination rates, as a function of the relative group size (below). Blue bars correspond to priority given to G 1 as the best strategy; yellow bars correspond to priority given to G 2 . (D) Quantitative effect of the best strategy as a function of the vaccination rate and the relative size of the groups. To the left of the thick red line, vaccination of the G 2 group first yields higher RD (numerically obtained values shown in each cell), and vice versa to the right of the line. Panels B, C, and D convey similar information, but underscore different effects. Parameter values common for all simulations are: N = 10 8 , M 11 = 6, M 22 = 3, M 12 = 2, IFR 1 = 0.1, IFR 2 = 0.5, d I = 15d, r 1 = 0.0667, r 2 = 0.0663, µ I1 = 6.67 · 10 −5 , µ I2 = 3.33 · 10 −4 , µ Yi = 0.01 · µ Ii , β SI = 1/d I = 0.0667, β RI = 0.1 · β SI , β RY = 0.5 · β RI , β SY = β RY · β SI /β RI and θ = 0.7.
We have explored the effect of the IFR in the optimal vaccination strategy as a function of the vaccination rate for two different group sizes (or contact matrices, recall the coupling between these two quantities through the closure relation), see figure 4. This case serves to compare the effects of diseases of different severity (as represented by the variation in the IFR) that, however, share their propagating ability and affect the same population.
At low IFR ratio, the strategy yielding the higher death reduction is vaccination of the group with the higher number of contacts, G 1 in our examples, since both groups experience similar effects of the disease. The situation is reversed as the IFR ratio increases, since as the disease impinges more severely on the G 2 group, its priority vaccination becomes more advantageous. As in the examples above, priority vaccination of the most vulnerable group is advised for low vaccination rates. Consistently, the optimal protocol changes at sufficiently high v, the threshold depending on the IFR and on the difference in contacts (c.f. group sizes) between the two groups.
As the results in this section illustrate, changes in different parameters are non-linearly related, and it does not seem possible to predict the optimal strategy without taking into account all variables and quantifying their effect through explicit simulations of the epidemic model. On the one hand, as vaccination rate increases, the proportion of vaccinated G 2 individuals grows even if this is not the priority group, since G 1 vaccination can be anyway completed earlier, and vaccination of G 2 individuals begins right after. Moreover, indirect protection of most vulnerable individuals through mitigation of G 1 -spreading capacity increases too. This is also reflected in the effect of reducing inter-group contacts. When M 12 is reduced, G 1 or simultaneous priority are more likely to improve G 2 death reduction at a given vaccination rate. Eventually, the G 1 -priority approach often becomes the optimal strategy for sufficiently high vaccination rates.

Application of S 2 IYRD to COVID-19
It has been proven that age is strongly correlated with severity and mortality of COVID-19 [34][35][36][37][38][39][40], a disease where reinfections are not rare [41][42][43][44]. Therefore, an assessment of different vaccination scenarios for SARS-CoV-2 propagation using the age-stratified S 2 IYRD model appears as a useful framework to detect major qualitative differences between vaccination strategies. In this section, we take different demographic profiles and empirically measured COVID-19 IFRs as case-examples. Populations will be divided into two groups through an age threshold. Contact matrices calculated from available data allow estimating parameters M ij for a given population split. Since there are multiple parameters that can be varied, we will discuss illustrative examples in this section and leave it to the interested reader to explore further cases through the interactive webpage [19] that we have developed.   [45,46]. IFR profiles as in figure 5(A) allow calculating the values of IFR in groups G 1 and G 2 , and the mortality and recovery rates of primary infections if the infectious period of the disease is known. For COVID-19, and in the framework of our model, it can be estimated at 13 days (see also supplementary information). The demographic composition of the population, as in the example shown in figure 5(B), is needed to calculate the values M ij . All estimated parameters vary if the threshold age that separates both groups, which affects their composition and the intra-and inter-group contacts, changes (see figure 5(C)). In the case of Spain, relative differences between the two groups in IFR and number of contacts are maximized for an age threshold at 80 years (down-right panel in figure 5(C)). We will keep this threshold fixed in this section to focus on the relationship between demographic profiles and vaccination strategies. A second case-example, with a threshold at 50 years, can be found in supplementary information (again, different thresholds and various other countries, India, Italy and Japan at the time of this writing, can be explored in the webpage [19]).
For all simulations performed, infection rates were computed as: where R 0 is the reproductive number of the epidemic disease under consideration, d I is the infectious period (the time interval during which the individual is infectious), and α i are two constants satisfying α i ⩽ 1. Since our model lacks an exposed compartment (E), we subsume under d I the average exposure period of infected individuals, estimated at 3 days [13]. Furthermore, available data suggest that patients with mild to moderate COVID-19 remain infectious no longer than 10 days after symptom onset [47]: therefore, we take d I = 13 days. The dynamics of reproductive number R(t) for COVID-19 is not a constant and depends on factors like population awareness, non-pharmacological interventions, etc [48,49]. For simplicity, we choose to fix R 0 = 1, since the COVID-19 reproductive number fluctuates all over the world around this value [50], likely pushed there by population dynamics [49]. The reinfection rate is defined through its relationship with the infection rate, β RI = α 1 β SI , and we fix α 1 = 0.01 in agreement with early estimations [44]. The last parameter is α 2 = 0.5, meaning that reinfected individuals recover twice as fast as individuals infected for the first time.
The reinfection rate is fixed at a 1% of the infection rate [44]. Finally, we consider the recovery rate of reinfected individuals to be twice as fast as that of individuals infected for the first time (see supplementary information for additional explanations). Note that the parameters here used correspond to the early propagation of COVID-19, so they would have to be reevaluated for new variants. Omicron, for example, is characterized by a different IFR (therefore different mortality and recovery rates), a different infectious period, and higher infection and reinfection rates. Preliminary evaluations of its infective rate increase it up to 4-fold with respect to previous variants [51]. Table 2 summarizes the empirical values obtained for four different countries which differ in their demographic composition and in their contacts (see supplementary information for raw data and other details). Spain has a negative growth in recent decades that yields a narrowing, cup-like demographic pyramid. It also bears the largest fraction of elder population, over six-fold that of the South African Republic (SAR). Israel's demographic pyramid has the Christmas-tree-like shape characteristic of young and rapidly growing populations, with a fraction of population above 80 years between that of Spain and the SAR. Such is the case of the elder group in the USA as well. The USA demographic pyramid displays a box-like shape that reveals a population of stable size. The IFR of all four countries is comparable for the elder group but differs for the young one. In turn, the fraction of contacts within group G 1 is high and comparable, while values of contacts between groups and within the G 2 group vary.
As previously done, we now examine three different vaccination strategies: priority to the over-80-years group G 2 , priority to the younger group G 1 , and simultaneous vaccination proportional to the population of either group. In the first two scenarios, vaccination of the corresponding second group starts once that of the priority group finishes (with a threshold θ = 0.7). Figure 6 summarizes the effect of different vaccinating protocols in the reduction of deaths as a function of the vaccination rate v for the four countries above. In all cases, vaccinating first G 2 always reduces mortality more than does any of the other two strategies, often more than doubling the effect of alternative protocols at low vaccination rates (see figure 6). Nevertheless, the difference is minor for the case of the SAR due to two main factors: the relatively high IFR of the younger group and, especially, the very low size of the elder group. As a result, most deaths come from G 1 , blurring the positive effect of protecting first the elder group. On the other hand, the vaccination of this latter group occurs early in the monitored year and, as explained, then continues with group G 1 . Therefore, its effects are not that different from simultaneous vaccination (and, actually, from priority to the G 1 group) in one-year time. Still, differences are observable at the initial stages of vaccination. Remarkably, simultaneous vaccination or priority vaccination to the G 1 group do not present major differences in any case. In the case of Israel, vaccination of the G 1 group first is slightly beneficial at low vaccination rates, while simultaneous vaccination is more efficient at high v; in the case of the USA, simultaneous vaccination is better at intermediate vaccination rates.
Vaccination rates, however, do make a significant quantitative difference in the difference of performance between the optimal strategy and the rest. In terms of reduction of fatalities, the advantage of G 1 priority barely reaches a 5% with respect to the worst-performing protocol. This even requires relatively high vaccination rates (around v = 1), and the improvement typically decreases with decreasing v.
Increasing the vaccination rate monotonically increases the reduction in the number of deaths, regardless of the protocol. The increase is almost linear for Spain and Israel, but accelerates with v for the SAR and the USA, especially. For example, if priority is given to the group G 1 , an increase in 0.05% in the vaccination rate from v = 0.2 to 0.25 rises RD in 2.3% points, while an increase from v = 0.9 to v = 0.95 causes a 4.9% increase in RD. The change in RD is also about two-fold for the other two strategies.
At sufficiently high vaccination rates, the three vaccination strategies reduce the differences in their effect. However, even in that situation, G 2 -prioritization makes it possible to achieve much earlier in time a death reduction comparable to that eventually attained with any of the other two strategies, due to the early protection of the most susceptible group. Given these observations, one may wonder: would it be possible to implement vaccination much faster than the spread of disease? For COVID-19, this hardly seems the case: even with very fast vaccination, the reduction of deaths is far from complete, figure 6. However, we dare to say that for slowly spreading diseases for which vaccines are already available, vaccination campaigns beating the spread of the disease should be possible. For example, using the default parameters of our online simulator [19], a 1% daily vaccination rate prioritizing elderly vaccination (G2 age group) results in a 67% reduction in infection and a 72% reduction in mortality using population data for the USA. Reducing the infection rate β SI from 0.075 to 0.040 and the infection rate β SY from 0.04 to 0.02 results in infection and death reductions of 99.9%, implying that very fast vaccination can really halt a disease slower than COVID-19.

Application of S 2 IYRD to influenza
To study a situation qualitatively different from the one described for COVID-19 in the last section, we obtained epidemiological data for influenza in the Australian state of New South Wales [52]. At present, ethical considerations restrict the evaluation of the effectiveness of influenza vaccine within a population to observational studies. However, during the 1991-1992 season, a single randomized placebo-controlled trial was conducted among patients aged over 60 in the Netherlands [53]. The trial revealed that the vaccine was effective by 50% against influenza confirmed by serology. However, the vaccine's efficacy was only 23% in subjects aged 70 years and above, implying a decline in effectiveness with advancing age [54]. To take this into account, we set that only 65% of vaccinated Susceptible individuals in the older age group generate defenses and pass to the Recovered class. To model a somewhat extreme situation regarding mortality, we set the age threshold between G 1 and G 2 at 20 years; doing so, we can set IFR 1 = 0, IFR 2 = 0.51 [52]: this means that the disease is not fatal for patients in the younger group. However, we have found that in this situation, for fast enough vaccination rates vaccinating the younger group first can lead to a significant reduction in mortality, figure 7, even though the disease does not pose a risk to this age group. Our analysis shows how higher vaccine efficiency and greater social contact can suggest strategies where the advantage of transmission suppression outweighs direct protection of the vulnerable.

Discussion
The goal of this work has been to analyze major effects of demographic population structure, disease fatality and vaccination rate in optimal vaccination strategies using a model with reinfections, and to provide a general scenario that can be applied to other situations. An important variable whose effect cannot be disentangled from other demographic factors is the contact between groups, which affects optimal vaccination strategies through its implicit relevance in disease spread. The dynamics have been evaluated using effective population-level contact matrices estimated in situations without mobility restrictions [17]. Contact matrices are however a key element that becomes modified in the presence of an ongoing pandemic due to non-pharmaceutical interventions that directly impact propagation dynamics [55], as implicitly shown by changes in mobility [56]. The model incorporates other simplifications, such as the implementation of a limited number of compartments and the use of only two groups.
Optimal vaccination roll-out depends on multiple variables. The precise value of the vaccination rate at which the advantage of vaccinating first the most susceptible group is lost in front of alternative strategies sensitively varies with the IFR of the specific disease. At high IFR, as for COVID-19, priority vaccination to the most connected group is only advantageous at remarkably high vaccination rates (above 1%, the precise threshold depending on the characteristics of the population). However, the relative advantage of the optimal strategy is reduced as the number of administered daily doses grows. The reduction is sensitive to demographic differences in a non-trivial way: while in the SAR and in the USA our model yields a quasi-equivalence among the three different strategies for a vaccination rate v ⩾ 1%, a significant difference persists for Spain and Israel, for example. Our results with two age thresholds (at 80 and 50 years) are consistent with the overall picture. Lowering the age threshold implies a reduction of IFR in the older group and a simultaneous increase in the average number of contacts per individual in both groups. In terms of death reduction, priority vaccination of the older group is systematically more efficient for high IFRs (as those of COVID-19). For diseases with lower IFR, as influenza, the general consensus is that children should be vaccinated first due to their high transmission capacity [16,57]: the analysis we have made using data from New South Wales strongly agrees with this picture, while highlighting the major role played by vaccine roll-out. It is likely that the optimal strategy for COVID-19 changes and resembles that for influenza if the disease becomes endemic and its effects turn milder once most of the world population has been either infected or vaccinated.
These results are consistent with those of other model-based studies showing that the optimal vaccination strategy balancing direct and indirect protection against COVID-19 is highly determined by vaccine supply and efficacy [13,14,58] and vaccination of the most susceptible group is more efficient under a broad range of situations [59,60]. Our study adds two elements to previous approaches. First, the introduction of reinfections, which are nonnegligible for COVID-19 and other corona viruses [4], shows the existence of long transients (or quasi-endemic states) that may predate the transition to a truly endemic state predicted for COVID-19 [10]. Second, the model is simple enough to allow the characterization of systematic effects due to, at least, group size, demographic composition and IFRs.
As other models used to inform optimal vaccination roll-out [13,14,58,59], our model assumes that individuals within a given group have on average the same number of contacts. More realistic models should consider heterogeneity in the number of contacts, since it has been shown that highly skewed contact distributions, with hubs (or, equivalent from a dynamic viewpoint, super-spreader individuals or events) have important effects in immunity thresholds [61] and in vaccination strategies. For example, the convenience of vaccinating network hubs first has been broadly discussed [62,63], though that strategy is hampered by the difficulty of identifying actual hubs with local information. An interesting alternative consists in vaccinating neighbors of randomly chosen individuals, since these neighbors have more contacts on average [64]. Still, this latter strategy cannot consider protective behaviors exhibited by individuals, which might be independent of their contact habits. Perhaps an improved strategy should include at once a sensible use of this last strategy coupled with individual-dependent IFR, a quantity that depends not only on age, but also on sex, comorbidities and social roles, among others. Regardless of the details of the model, the most efficient immunization strategy is conditioned by the availability of vaccines at each time, and optimization is particularly critical when availability is low-to-medium, while different strategies converge at sufficiently high vaccination rates.
We have focused on the reduction in the number of fatalities under different vaccination protocols, since this is a short-term benefit that can be easily quantified. However, different vaccination protocols also affect hospital occupation or the number of infections, and the variation can be very large depending on population structure, contact habits, and vaccination rate. Our model allows a preliminary assessment of strategies that could minimize the impact of the disease in the health care system. For example, the number of hospitalizations is proportional to the sum aI + bY and, since the second term corresponds to individuals with partial immunity, typically experiencing milder disease effects, a ≫ b. Increasing the relative size of class Y should therefore minimize the number of severely ill individuals. This minimization can be attained by increasing the vaccination rate (which depends primarily on vaccine availability and vaccine uptake), but also by increasing the maximum attainable fraction of vaccinated individuals, perhaps through devoted awareness campaigns [65]. Although quantitative benefits of reducing the number of infections are difficult to evaluate due to multiple additional variables involved, they have to be kept in mind as our knowledge of the adaptive strategies of SARS-CoV-2 improves. Beyond the burden they cause in the health care system, more infections entail an enhanced number of circulating viral variants and therefore a higher probability of emergence of strains able to escape the protective effect of vaccines. In the mid and long term, it might be advisable to seek optimal vaccination strategies that simultaneously minimize both the number of fatalities and the number of infections.