Optimal COVID-19 vaccine prioritization by age depends critically on inter-group contacts and vaccination rates
Abstract
The limited availability of COVID-19 vaccines has prompted extensive research on optimal vaccination strategies. Previous studies have considered various non-pharmaceutical interventions, vaccine efficacy and distribution strategies. In this work, we address the combined effects of inter-group contacts and vaccination rates under contact reduction, analysing the Spanish population’s demographic and age group contact patterns and incorporating reinfection dynamics. We conduct an exhaustive analysis, evaluating 362 880 permutations of nine age groups across six vaccination rates and two distinct, empirically quantified scenarios for social contacts. Our results show that at intermediate-to-high vaccination rates with unrestricted social contacts, optimal age-based vaccination strategies only slightly deviate from older-to-younger prioritization, yielding marginal reductions in deaths and infections. However, when significant reductions in social contacts are enforced—similar to the lockdowns in 2020—there are substantial improvements, particularly at moderate vaccination rates. These restrictions lead to a transition where infection propagation is halted, a scenario that became achievable during the pandemic with the observed vaccination rates. Our findings emphasize the importance of combining appropriate social contact reductions with vaccination to optimize age-based vaccination strategies, underscoring the complex, nonlinear dynamics involved in pandemic dynamics and the necessity for tailored context-specific interventions.
1. Introduction
Compartmental models have been broadly used to characterize the dynamics of infection propagation in a variety of scenarios. Though most studies have addressed the impact of refining basic susceptible–infected–recovered dynamics [1] to include a multiplicity of additional states [2], others have thoroughly explored the effect of various external interventions. The use of such models has become widespread, addressing issues such as the effect of contention measures [3–5], the immune state of the population [6], population density [7] or vaccine rollout [8], among others. One of the major strengths of modelling approaches is that they offer the possibility of comparing various scenarios under similar conditions. Altogether, such studies have refined our understanding of epidemic dynamics and the role played by multiple relevant variables.
The availability of a vaccine against SARS-CoV-2 largely before the infected population attained a sufficiently high immunity level to halt propagation prompted the study of optimal protocols to administer available vaccine doses [9–11]. Due to the high, positive correlation of SARS-CoV-2 mortality with the age of the infected individual, there was a largely preferred protocol of administration where vaccination of the elderly was prioritized, to continue in decreasing age order once most vulnerable groups had been protected. Results obtained through models incorporating a variety of scenarios, but implementing COVID-19 features, systematically supported that mortality and years of life lost were minimized if the population over 60–70 years old was given priority [11–15]. Some of these studies used data from sociological studies elaborated in the absence of social contact limitations or, in the best case, implemented such limitations in a phenomenological way, using for instance effective values of the reproduction number in the simulations [12,14]. However, modifications in social contact patterns during an epidemic event can change due to both non-pharmaceutical interventions and as a result of unsupervised responses of individuals [16,17] in a highly heterogeneous way [18]. During the COVID-19 pandemic, patterns of contact between individuals were severely modified [19], differently affecting citizens involved in public services, children or elders, for example.
In this contribution, we carry out an exhaustive exploration of age-based vaccination protocols by measuring the reduction in deaths and infections obtained under various scenarios. Our baseline situation considers an age-structured population with nine age groups and a specific demographic structure: the Spanish population. In our first scenario, no restrictions to social contacts are implemented, though empirical data regarding inter- and intra-group contacts are used [20]. Under the previous conditions, we carry out simulations with a compartmental model that takes into account partial immunity and, therefore, reinfection of individuals, as well as vaccination at various rates. All possible orderings of the population, attending to the age of the different groups defined in our study, are simulated. This is a total of possibilities for each vaccination rate. As a result, we are able to identify the best possible protocol and quantify its benefits with respect to the baseline case. The whole study is repeated when lockdown measures are implemented, using the results of two studies [21,22] that have measured contact reduction in different age groups. Our findings align with results in the literature under no contact restrictions, showing clearly that the largest benefit in terms of mortality reduction is obtained when vaccination prioritizes elderly people. However, we obtain a substantial advantage of vaccinating first people in their thirties under lockdown conditions. In this second situation, reductions in mortality and the number of infections can be simultaneously optimized under similar vaccinating strategies. We also show that there is a non-trivial interaction between the degree of social contact and the vaccination rate and that only a suitable combination of both variables can yield the largest beneficial effects.
2. Model and data
2.1. SIYRD: a model with reinfection and vaccination
In a previous study [23], we introduced SIYRD, a compartmental model with five different classes: Susceptible (S), Infected (I), Reinfected (Y), Recovered (R) and Dead (D) individuals. The model implemented two main novelties: vaccination (S individuals move to class R upon vaccination) and reinfection (R individuals can move to the reinfected class Y). Only vaccination of susceptible individuals is considered, and an equivalence between vaccination and disease overcoming regarding the immune state is assumed. Empirical studies have found no significant differences in recovered individuals between the immune effects elicited by vaccination or by infection [24,25], thus supporting this hypothesis.
Numerical and analytical studies of the SIYRD model showed that the introduction of reinfections causes the appearance of very long transients where the disease behaves as quasi-endemic. Extending the model to include a synthetic population divided into two age groups revealed that at high vaccination rates—whenever disease impact intensifies with age—the optimal vaccination strategy shifts from prioritizing the elderly to prioritizing younger individuals. Here, our main goal is the exhaustive exploration and quantification of this observed effect in realistic scenarios.
To this end, we extend the SIYRD model to encompass nine different age groups. This allows us to implement the demographic structure of the population and inter- and intra-group empirical contact matrices in two situations: without and with limitations to social contacts. The numerical analysis uses data from Spanish demography and COVID-19 data, though it is easily applicable to other diseases, countries or group-contact scenarios.
2.2. SIYRD: a model implementing demographic and group-contact structure
Let us define as the number of susceptible individuals in age group , with including from newborns to individuals of age 9, individuals between ages 10 and 19 and so on until . The last group, , includes all individuals of age 80 and older. The number of individuals for other classes and all age groups is analogously described by , , and , .
Equations (2.1)–(2.5) describe the dynamics of each of the nine groups in the SIYRD model, as schematically shown in figure 1:

Figure 1. Schematic of the SIYRD model. Five different compartmental classes are considered: susceptible, S; first infected, I; recovered (either recovered from I or vaccinated from class S), R; reinfected (infected after vaccination or after recovering from a first infection), Y; dead, D. Individuals are further classified in one out of nine different age groups. Contacts between infected individuals (I or Y) and those that can be infected (S or R) occur between all possible age pairs with a weight depending on contagion parameters, demographic structure and empirical contact matrices (orange arrows in the schematic).
Parameters are rescaled such that the time unit of numerical simulations is one day. The meaning of parameters and the estimated values based on empirical data are discussed in the next subsection.
2.3. SIYRD model parameters
Parameter values chosen for our numerical simulations consider the early propagation of COVID-19 and the Spanish population as an example. A summary of the main characteristics of the latter used to feed model parameters is represented in figure 2. Though the qualitative results obtained are robust to variations in virus and population characteristics, quantitative results would need to be re-evaluated for different situations. Omicron, for instance, is characterized by a milder impact fatality rate (thus yielding different mortality and recovery rates) [27,28], a different infectious period and higher infection and reinfection rates than early COVID-19 circulating strains. Also, epidemic dynamics should differ between countries with different social contact habits or regions with an expansive population pyramid—unlike the constrictive Spanish population pyramid.

Figure 2. SIYRD parameters for COVID-19 in Spain. (a) Demographic pyramid for the Spanish population in 2020; blue corresponds to male population, and dark red to female population; in the vertical axis, populations are grouped in 5 year intervals, starting with 0–4 years at the bottom. (b) Population size considering nine age groups in 10 year intervals. (c) Effective values of the number of daily contacts for each group, considering the original (grey) and reduced (purple) contact matrix. In gold (log scale, axis on the right-hand side), we represent the infection fatality rate for each age group in Spain [26]. (d) Heatmap of the contact matrix representing the logarithm of the average number of contacts between individuals [20] for the nine-group-stratified Spanish population. (e) Reduced contact matrix, as in (d), but in a scenario with limitations to the number of social contacts (see text for details).
2.3.1. Infection rates
The separation between primary infections (I) and reinfections (Y) entails four different infection rates: and correspond to the infection rates of susceptible individuals due to contacts with primary-infected (I) and reinfected (Y) individuals, respectively; and are the analogous rates for recovered individuals. Note that since we do not contemplate an exposed class in the model, incubation periods are effectively included in the infection and recovery rates. Although the precise values of some of these infection rates might be difficult to estimate in general, there are sensible relationships among the parameters that should hold, on average, over the population. We assume individuals in the reinfected class Y bear a lower viral load due to their previous exposure to the virus and, therefore, are less infective than primary-infected individuals in class I, both towards susceptible S and towards recovered R individuals. This implies that reinfection rates are smaller than primary infection rates: and .
In the same vein, 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 prior infection or to vaccination. This applies both to primary infections, , and to reinfections, . 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,
Lacking specific data that suggest otherwise, we assume that COVID-19 transmission rates , , and are independent of the age group considered. Precise values can depend on different COVID-19 strains but are largely country-independent.
Specifically, infection rates were estimated as follows:
where parametrizes the transmissibility of the epidemic disease under consideration, is the infectious period (the time interval during which the individual is infectious), and are two constants satisfying .
In general, the reproductive number varies along epidemic propagation, usually starting at high values at the beginning of an epidemic wave and decaying towards one as the epidemic progresses. Such has been the case with the COVID-19 reproductive number, which has been observed to fluctuate all over the world around its critical transmissibility value , with occasional departures at the emergence of new strains [29] or when social contact restrictions are put in place [30]. Under free epidemic propagation, has consistently tended towards 1 [31], probably due to changes in individual risk aversion in response to the epidemic state [17]. In view of this evidence, for simplicity, we choose to fix , as representing the baseline long-term trend in .
The infectious period results from the sum of two contributions: the average exposure period of infected individuals, estimated at days [12], and the infectious period proper. Data suggest that patients with mild to moderate early COVID-19 remain infectious no longer than 10 days after symptom onset [32]. Therefore, we take days.
Finally, we have chosen , modelling the first COVID-19 wave [33] before significantly mutated strains appeared, and , implying that reinfected individuals recover twice as fast as primary-infected individuals.
2.3.2. Mortality and recovery rates
The quantitative estimation of the primary mortality rate depends on the infection fatality rate (IFR) and on the infectious period of a disease. The values of the IFR for COVID-19 in Spain [26] are represented in figure 2c, though measures in different world countries return comparable values [34,35]. For primary-infected individuals, can be calculated from the IFR in each age group, , and the infectious period as
where is given by the ratio between the number of fatalities and the number of infections in a given group . The recovery rate from a primary infection is analogously estimated from
Mortality rates of secondary infections, , are set to for all nine groups. This, by definition, establishes an endemic disease where, asymptotically, all individuals that either have been vaccinated or have survived a primary infection belong to R or Y classes. Considering that this model does not incorporate the possibility of mutations in the circulating strain, assuming that secondary infections do not cause fatalities seems sensible. Still, this limits the applicability of the model to long time intervals, where changes in the health state of individuals or in the circulating virus might take place. Since the goal of the model is to establish optimal vaccination protocols in a mostly susceptible population (at the onset of the disease propagation), a timescale of about 1 year suffices, thus justifying the study of the limit case . Finally, we assume that the recovery rate of primary-infected and reinfected individuals is identical, neglecting possible small differences [15,36].
2.3.3. Vaccination rates
Vaccination is implemented through a parameter that represents the fraction of susceptible individuals in group vaccinated per time unit. The rate is multiplied by a function to indicate that vaccination is achieved only in a fraction of individuals in that group. Since it has been shown that qualitative results are independent of the specific functional form of [23], we assume that vaccination takes place at a constant rate until the fraction , with for all groups , is achieved.
We explore priority vaccination protocols determined by a specific order of the nine age groups. Specifically, the prioritized population group is vaccinated at a rate until the fraction of susceptible individuals reaches in that group. When that happens, vaccination of the next prioritized group begins—provided the previous threshold has not yet been reached through natural infections, a situation that may hold for too low values of . Vaccination continues with other groups until the vaccination protocol is completed. The vaccination rate is a variable that can be explored in a range of values depending on vaccine supply. We consider six different vaccination rates: , , , , and representing the fraction of population vaccinated per day.
In the case of the Spanish population, with a total around of 47 million individuals, slow vaccination rates of and would correspond to 23 500 and 47 000 daily vaccines administered, while the fastest rate considered in the simulations, , implies that 940 000 individuals would be vaccinated per day. The highest 7 day average daily vaccine administration actually achieved—during the Omicron wave in late —was [37].
2.3.4. Demographic structure and social contacts
The demographic structure of the population (figure 2a) sets the size of each age group (figure 2b), a quantity that serves to normalize the value of . The elements of the contact matrix are the number of contacts an individual of group has with individuals of group , thus weighing the effect of intra-group () and inter-group () contacts in contagion rates.
We have used different sources to quantify contact matrices in a situation where contacts are not limited and when there are contact restrictions. In the former case, with no restrictions to contacts, data for Spain could be obtained and used as found in a previously published comprehensive study [20] (see also figure 2d). Data on contact reduction under lockdown for Spain were not available, however, so we resorted to two independent studies for The Netherlands [21] and England [22] to estimate a realistic matrix under contact reduction. Both studies compared the reduced matrix to a baseline, pre-pandemic situation. By using the difference between contacts under lockdown versus contacts without restrictions, we illustrate how lockdown measures impact depends on the age group, affecting in turn optimal vaccination protocols.
The study with the Dutch population [21] is a cross-sectional survey in which participants reported the number and age of their contacts the previous day. It was conducted at three different times. In this work, we are interested in two of these times, namely data collected from February 2016 to October 2017 (pre-pandemic) and data collected in April 2020, after strict physical distancing measures were implemented. These included closing daycare centres, schools, universities, cafés, pubs, restaurants, theatres, cinemas and sports clubs, as well as cancelling events with more than 10 participants. The advice to citizens was to work from home whenever possible and to maintain a 1.5 m distance from others outside the household. Overall, the effective number of contacts per person was found to be reduced by 69%. The distancing measures taken are very similar, only slightly more permissive than those enforced in Spain in April 2020. Data in [21] can be directly applied to our case once we have taken care of the differences in age grouping, which simply implies a re-calculation of the number of contacts per age group using weights given by the size of each group.
The baseline for the effect of lockdown measures in the English study [22] was a previous population-based prospective survey of mixing patterns carried out in eight European countries [38]. Interestingly, this study showed that mixing patterns and contact characteristics were remarkably similar across countries (Spain was not included in the study). From March to June 2020, the lockdown in England entailed work from home, schools closed, restaurants closed and mandatory masks in some areas. The mean number of contacts was about 75% fewer than at pre-pandemic time [22]. Data had to be processed again just to meet the age grouping used in our study.
Reduced matrices obtained from the Dutch and the English studies show qualitatively similar patterns and are quantitatively comparable. They can be found in the file SuppFile_ConstrainedMatrix.xlsx, available at [39]. We have averaged the two values obtained for each of the matrix elements, resulting in values that we have used in our simulations as a case example of contacts under restrictions (see appendix A). The final matrix that we implemented is represented in figure 2e.
2.4. SIYRD numerical implementation
To account for the possibility of different emerging timescales in the dynamics due to the action of various mechanisms in the model, we have numerically integrated the equations using the variable-step, variable-order implicit solver ode15s from the MATLAB ODE suite [40]. This solver is appropriate for stiff problems with different timescales. We have not observed unrealistic crossings of variables to negative values [41]. 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. Data and relevant code for this research work are stored in GitHub [42] and have been archived within the Zenodo repository [39].
Given the fixed parameters (all rates specified in previous sections), a vaccination rate and a social scenario (unrestricted or reduced physical contact), simulations are run for a time interval corresponding to 1 year using all possible permutations for vaccination of age groups. Each vaccination protocol (permutation of the nine age groups) is written as a vector with a specific group order. For example, indicates that group 7 (individuals from 60 to 69 years old) is vaccinated in the first place, followed by group 8, then 9 and so on until either the simulation finishes or 70% (or over) of individuals in a given group have been infected. This latter case is common at low vaccination rates; for example, with the parameters chosen and for , and if vaccination starts with the oldest group, all other groups reach a fraction of over 70% infected individuals before their vaccination turn arrives. We represent this situation as to indicate that only group 9 could be vaccinated, while the others reached their ‘immunity thresholds’ due to natural infections. There are different orderings, which are exhaustively explored in this contribution. The efficacy of each ordering (each vaccination protocol) is quantified through the percentage reduction RD% in the total number of deaths and the percentage reduction RI% in the total number of infections, as compared to the baseline of no vaccination.
Figure 3 shows the dynamics of the SIYRD model in four representative scenarios, using the parameters described earlier. This figure demonstrates how a well-coordinated combination of social contact reduction and an optimal age-based vaccination protocol can significantly reduce the total number of fatalities and infections (compare figure 3a,d). Additionally, the number of simultaneously infected individuals is much lower in the latter scenario, thereby lessening the strain on the healthcare system. These dynamics will be systematically explored in the following section.

Figure 3. Dynamics of the SIYRD model in representative scenarios. The -axes represent time in days since the start of epidemic propagation. The inset values show the percentage reduction in deaths (RD) and infections (RI) at the simulation endpoint (1 year), calculated relative to the baseline scenario of no vaccination. In this figure, the vaccination rate is . The two upper plots illustrate the dynamics under an older-to-younger vaccine administration protocol, while the two lower plots correspond to dynamics using a representative protocol that maximizes RD at intermediate vaccination rates (see table 2). The scenarios depicted are as follows: (a) unconstrained social contacts with an older-to-younger vaccination protocol; (b) constrained social contacts with the same vaccination protocol as in (a); (c) unconstrained social contacts with the protocol ; and (d) constrained social contacts with the same vaccination protocol as in (c).
3. Results
3.1. Optimal ordering under unconstrained contacts
We have first explored the situation of epidemic propagation using a contact matrix that reflects the behaviour of a population when no limitations to social contacts are imposed [20]. This scenario has been explored in various studies [12,13], despite the fact that assuming no changes in contact behaviour is unrealistic when a threatening contagious disease is spreading: even in the absence of imposed institutional measures, individual behavioural changes typically arise as a result of the perceived risk of contagion [43–45]. Still, the case of no changes in measures that hinder contagion, spontaneous or imposed, serves as a baseline to compare different scenarios under otherwise identical conditions.
3.1.1. Reduction in fatalities and infections
The percentage reduction in the number of fatalities and infected individuals depends significantly on the vaccination rate and on the vaccination ordering of age groups. We present a statistical summary of these reductions in figure 4 through a violin plot of RD% and RI% values, for the ensemble of all possible orderings. As could have been expected, increasing the vaccination rate improves, on average, the effect of vaccination, diminishing the total number of both deaths and infections. It is important to note that a bad choice of the vaccination protocol can, however, suppress the advantage of increasing the vaccination rate, affecting especially the total number of deaths.

Figure 4. Distribution of the percentage reduction in deaths (a) and infections (b) for all possible permutations of the vaccination protocol, for each vaccination rate, in the case with unconstrained contacts. Dots represent average values, and black bars span from maximum to minimum effect.
Too low vaccination rates, , change only slightly the course of the epidemic, which is dominated by its unconstrained propagation. However, the effect of vaccination changes drastically for larger rates: for and higher, the order of the administration protocol can cause benefits that vary from about only 10% improvement (in both deaths and infections) to over 50% decrease in total number of fatalities when the optimal ordering is chosen.
3.1.2. Optimal protocols differ from elder-to-younger ordering
While for an insufficiently high vaccination rate the specific ordering of vaccination of age groups does not significantly affect the total number of deaths and infections, the situation changes as increases. For average vaccination rates, and , there is an optimal ordering to minimize the number of deaths that coincides with the strict age ordering from older-to-younger groups , widely favoured for COVID-19. However, for the higher rates tested, and , there are protocols that outperform the oldest-to-youngest (strict age order (SAO) strategy) ordering.
Optimal strategies imply an important reduction in infections with respect to strict age ordering, although the results are similar in terms of RD% (see table 1).
protocol | RD% | RI% | SAO RD% | SAO RI% | ||
---|---|---|---|---|---|---|
best RD | 0.05 | [9 0 0 0 0 0 0 0 0] | 9.8 | 1.0 | 9.8 | 1.0 |
0.1 | [9 8 0 0 0 0 0 0 0] | 18.5 | 1.9 | 18.5 | 1.9 | |
0.5 | [9 8 7 0 0 0 0 0 0] | 54.6 | 9.0 | 54.6 | 9.0 | |
1.0 | [9 8 7 6 0 0 0 0 0] | 65.6 | 17.3 | 65.6 | 17.3 | |
1.5 | [7 9 8 6 5 0 0 0 0] | 69.4 | 25.5 | 69.2 | 25.3 | |
2.0 | [3 2 4 9 8 7 6 5 0] | 71.8 | 51.2 | 70.7 | 33.2 | |
best RI | 0.05 | [3 2 1 5 6 7 8 9 0] | 1.3 | 1.1 | 9.8 | 1.0 |
0.1 | [3 2 4 6 7 8 9 0 0] | 2.5 | 2.2 | 18.5 | 1.9 | |
0.5 | [3 4 5 6 7 8 9 0 0] | 8.7 | 11.5 | 54.6 | 9.0 | |
1.0 | [3 4 2 1 5 6 7 8 9] | 14.6 | 25.4 | 65.6 | 17.3 | |
1.5 | [3 2 4 1 5 6 8 7 9] | 28.0 | 43.5 | 69.2 | 25.3 | |
2.0 | [3 2 4 1 5 6 7 8 9] | 52.9 | 62.3 | 70.7 | 33.2 |
At any vaccination rate, we observe that the protocol that minimizes the number of fatalities does not coincide with the protocol that minimizes the number of infections, as can be seen by comparing the upper (RD) and lower (RI) halves of table 1.
3.1.3. Many age-ordering protocols yield similar advantages
As the spread of RD% and RI% values in figure 4 shows, there are significant quantitative differences among protocols. In this section, we explore in further detail how similar are orderings that perform comparably regarding the reduction caused in the two former variables. We focus our analysis on the subset of protocols that cause a reduction of 95% or higher than the optimal protocol, in both RD% and RI%. We recall here that, for each vaccination rate, we have analysed 9! = 362 880 different orderings; the top 5% performance includes several thousand different possibilities. It is important to clarify that not all these possibilities are different in practice. For low vaccination rates, in particular, epidemic propagation proceeds faster than vaccination, so that just one or two groups can be vaccinated (before all the groups reach the 70% threshold). Therefore, what we call protocol (as in table 1, best RD for ) actually contains 8! = 40 320 protocols (corresponding to all possible orderings of groups 8 to 1 from the second to ninth position) that, in practice, become indistinguishable.
Given the ensemble of ordering protocols performing at 95% the optimal one, we calculate the fraction of age groups that appear in the first, second and third vaccinating positions. For example, if the ensemble contains only two protocols, say and , groups 9 and 8 appear once in the first and once in the second position, and group 7 appears twice in the third position. Figure 5 summarizes group prioritization for these three positions and each vaccination rate. Let us examine each vaccination rate in turn.

Figure 5. First, second and third position distribution for each vaccination rate considering top protocols with a level of performance higher than 95% of the best one in terms of (a) RD and (b) RI. The colours indicated are for the first position; the second has a lighter shade, and the third is even lighter.
For , as described, epidemic dynamics are dominated by free contagion from infected to susceptible individuals, and the flux from susceptible to recovered due to vaccination is mostly inconsequential: contacts do not play any major role; reduction of casualties is maximal if the prioritized group is the eldest, and reduction in the number of infections (obtained through direct vaccination) does not depend significantly on which group is vaccinated first. As shown in table 1, the RI of the best protocol is only 1.1%, while the RI of the SAO protocol achieves 1.0%, representing a minor difference.
However, prioritizing vaccination for the eldest results in the smallest reduction in infection numbers. This is because this group has the highest mortality rate, making vaccination less effective at lowering infection rates; it is also the group with fewer contacts. Such minimal influence of the eldest group in the RI is reflected in figure 5. For , there is no strategy prioritizing this group included among the best-performing RI protocols (RI within the optimal, , and 5% variation, 1.045).
This explains why in figure 5b group 9 is absent from the first position, while the advantage of beginning with any other group yields comparable reductions in RI.
Quantitative changes appear for (figure 4) since the reduction in the fraction of casualties reaches 18.5% for the best protocol, which is the SAO protocol . When all protocols causing reductions of 17.58% or higher are considered, up to three groups can be vaccinated, but a strong preference for protocol is observed: only strategies vaccinating in this age order are able to perform among the best 5% in RD. As for RI, the situation is very similar to that for the previous vaccination rate. Contacts do not affect the dynamics in any significant way, which continues to be dominated by free propagation plus the direct effect of vaccination: twice the number of vaccines doubles the reduction in the number of infections, with an irrelevant dependence on the vaccination order.
The situation changes qualitatively for , where signs of a nonlinear dependence between vaccination and observed benefits show up. The difference in performance among protocols becomes larger for RD, while all of them still perform similarly for RI (figure 4). For RD, the SAO protocol continues to be the preferred option, now allowing the vaccination of three groups. Protocols with groups 7, 6 and 5 in the third position are found among the best 95% orderings (that is, causing a reduction of at least 51.87% in RD). Interestingly, some protocols that begin with the vaccination of the second-oldest group, as , perform now better, regarding RD, than some protocols that begin with the oldest group, as .
Together with the previous rate, is the rate where the effect of the order of vaccination has the strongest effect. The difference between the optimal and the worst protocol rises to about 50% for RD; still, the dispersion of RI values remains small. At the same time, protocols different from SAO offer comparable advantages, and even protocols that begin with vaccination of groups 1, 2 or 3 reach RD values over 60%. Optimal reductions in the number of infections are achieved when vaccination starts with groups of intermediate ages, followed by children and leaving the eldest groups at the end of the protocol. These orderings, however, perform very poorly regarding RD, so they are not an option for COVID-19 at this vaccination rate.
When increases further, the effects of contacts, IFR and vaccination rate intermingle in a complex way, eventually leading to orderings that turn out to be optimal both for RD and RI when vaccination begins with groups 3, 2 and 4. At this point, the spread of RD values diminishes, and vaccination speed overcomes the basal dynamics of contagion, yielding an advantage when more connected groups are vaccinated first. This trend continues for values of the vaccination rate above , the largest value used here in our simulations.
3.2. Optimal ordering under social contact constraints
Limitations in social contacts, reducing the average number of contacts an individual has, should impact disease transmission. As explained in §2.3.4, we have used two empirical studies carried out under lockdown measures to estimate contact reduction between all age-group pairs. This allows us to quantify the joint effect of vaccination and social contact restrictions and to compare the result with the baseline case in the previous section. In the case of COVID-19, we observe that the effect the union of the two dissimilar mechanisms may have had in its spread is not only quantitative. According to our model and as we show in the following, the combination of an optimal vaccination protocol with a sufficient reduction in the number of effective contacts may have caused locally a halt in epidemic propagation. An illustration of this effect through epidemic dynamics of the model can be seen in electronic supplementary material, figures S1–S4.
3.2.1. Reduction in fatalities and infections
As in the previous case, we analyse the distribution of RD% and RI% considering all possible orderings at each vaccination rate. The results are represented as violin plots in figure 6 and should be compared with those in figure 4.

Figure 6. Distribution of the percentage reduction in (a) deaths and (b) infections for all possible permutations of the vaccination protocol, for each vaccination rate, in the case with restricted contacts. Dots represent average values, and black bars span from maximum to minimum effect.
Though general considerations are similar to those discussed for the case with unconstrained contacts, two main differences emerge. First, a significant dispersion in RD% values occurs even at the lowest vaccination rates, highlighting on the one hand a nonlinear interaction among adjustable variables (contact limitation and vaccination rate) and, on the other hand, the relevance of properly selecting the optimal vaccination order. With reduced contacts, infection progresses slower, and therefore there is more ‘time to act’, explaining why there can be bigger differences among strategies. Second, even at intermediate vaccination rates (), RD% and RI% values attained are above the fraction of vaccinated individuals; for , both infection propagation and mortality are fully suppressed. This indicates that, for , epidemic propagation can be halted through a suitable selection of vaccination ordering and measures to limit social contact. Interestingly, nonlinear effects produce a relatively fast transition from vaccination rates at which there is a direct proportionality between the number of vaccinated individuals and the total reduction in fatalities plus infections, and vaccination rates able to fully inhibit infection propagation.
3.2.2. Optimal protocols substantially differ from older-to-younger ordering
A significant reduction of contacts between individuals at intermediate vaccination rates causes a full inhibition of infection propagation due to nonlinear effects between contact reduction and vaccination. The selection of the optimal vaccination order, once the two latter variables are fixed, becomes essential to maximize their synergistic effect. Figure 7 and table 2 summarize the results regarding optimal orderings. Remarkably, even at intermediate vaccination rates the optimal ordering substantially differs from the baseline, SAO ordering. At , vaccination should ideally start with the youngest group, followed by groups 4 and 3, to yield an improvement in RD% of about 10% with respect to SAO ordering; under the same protocol, the improvement in RI% more than doubles. Remarkably, both variables are significantly improved under protocols that coincide with the first four groups to be vaccinated (see table 2). The epidemic is practically fully contained at vaccination rates and higher.

Figure 7. First, second and third position distribution for each vaccination rate considering top protocols with a level of performance higher than 95% of the best one in terms of (a) RD% and (b) RI% under restricted contacts. The colours indicated are for the first position; the second has a lighter shade, and the third is even lighter.
v | protocol | RD% | RI% | SAO RD% | SAO RI% | |
---|---|---|---|---|---|---|
best RD | 0.05 | [9 8 0 0 0 0 0 0 0] | 30.27 | 3.01 | 30.27 | 3.01 |
0.1 | [9 8 7 0 0 0 0 0 0] | 47.33 | 6.30 | 47.33 | 6.30 | |
0.5 | [1 4 3 2 9 8 6 7 5] | 86.11 | 74.67 | 77.42 | 33.62 | |
1.0 | [1 4 3 2 5 6 8 7 9] | 96.80 | 92.06 | 91.71 | 80.03 | |
1.5 | [4 5 3 2 8 9 7 6 1] | 99.99 | 95.72 | 95.18 | 89.35 | |
2.0 | [3 1 2 6 8 4 9 7 5] | 99.99 | 98.46 | 96.37 | 91.29 | |
best RI | 0.05 | [4 3 5 6 7 8 9 0 0] | 6.34 | 5.59 | 30.27 | 3.01 |
0.1 | [4 3 5 6 7 8 9 0 0] | 12.16 | 11.57 | 47.33 | 6.30 | |
0.5 | [1 4 3 2 5 6 7 8 9] | 80.70 | 81.11 | 77.42 | 33.62 | |
1.0 | [1 3 4 2 5 6 7 8 9] | 96.79 | 92.07 | 91.71 | 80.03 | |
1.5 | [2 6 5 3 4 8 1 7 9] | 99.99 | 99.98 | 95.18 | 89.35 | |
2.0 | [3 1 2 6 8 4 7 9 5] | 99.99 | 99.99 | 96.37 | 91.29 |
Note that the response of the different groups and of the whole population at a given vaccination rate under limited social contacts is not a simple translation of that obtained at higher vaccination rates under unlimited contacts, as comparison of figures 5 and 7 shows (see also electronic supplementary material, figures S1–S4). One reason is the changes in the contact matrix, with different group ages leading to the number of contacts under the two different conditions analysed. In particular, it is only in a situation with limitations in social contacts that the vaccination of the youngest and intermediate age groups yields high advantages in both RD% and RI% and is therefore justified.
4. Discussion
In this study, we conducted an exhaustive exploration of all possible vaccine administration protocols across nine age groups, six vaccination rates and two scenarios informed by field data, with and without social contact reductions. Using the demographic structure of the Spanish population, we simulated epidemic dynamics through a model that includes reinfections, potentially leading to long-term endemic states [23], as observed in COVID-19. Our findings highlight that the effectiveness of vaccination protocols is highly dependent on vaccination rates and the implementation of measures to limit social contact.
With unrestricted social contacts, a vaccination rate of 1.0% of the population per day makes a strict older-to-younger vaccination order the optimal strategy to minimize casualties. However, if the vaccination rate increases to 2%, prioritizing the 20–29 age group, followed by the 10–19 and 30–39 age groups, becomes more effective. Still, even at the highest explored vaccination rates (), prioritizing younger groups never yields more than a 4% reduction in RD without limiting social contacts. This quantitative advantage, while notable, should be interpreted cautiously, given that our model excludes other potentially relevant factors, such as sex-dependent responses to infection, age-dependent vaccination efficacy, values different from 70% maximum vaccination within each age group, or variations in the value of transmissibility. Nonetheless, the fact that some countries or regions experienced vaccination rates exceeding 2% per day underscores the potential of alternative protocols for minimizing fatalities, even if the improvements are limited. Other studies have presented results qualitatively consistent with our findings when differential allocation of vaccine doses to age groups is considered. For instance, a study on the Delta variant’s spread in Australia demonstrated the benefits of prioritizing socially active groups for vaccination over vulnerable groups, showing greater benefits as vaccination coverage increased [46]. These results align with findings for other contagious diseases; in the case of influenza, vaccinating schoolchildren and adults aged 30–39 offers significant protection by reducing transmission from the most active spreaders to the wider population [47]. However, caution is advised, since in our work we have not considered the existence of asymptomatic individuals for which the infection goes undetected. Some studies taking this aspect into account have highlighted the risks of high vaccination rates in diseases with a high proportion of asymptomatic individuals, who can unknowingly contribute to the spread of infection [48].
The situation changes when social contacts are restricted. Our results demonstrate a complex interaction between social contact patterns and vaccination rates, which directly influences the optimal vaccine rollout strategy. Preference for an older-to-younger vaccination order shifts under a lockdown scenario: at a 1.0% vaccination rate, prioritizing the youngest age group (0–9 years old), followed by the 30–39 and 20–29 groups, with the oldest group (over 80) last, is optimal. If the vaccination rate increases to 2%, the recommended order changes to 20–29, followed by 0–9 and 10–19, leaving the 40–49 group in the last position.
Increasing vaccination rates and limiting social contacts are distinct and non-equivalent strategies for controlling the spread of contagious diseases. The effectiveness of each strategy varies both quantitatively and qualitatively across different age groups. This finding aligns with a mathematical model analysing the reopening of New Zealand’s borders, which indicates that even with a highly effective vaccine, additional public health measures are necessary to manage COVID-19 risks [49]. Similarly, work using agent-based models has shown a positive synergistic effect between mild and social restrictions and intermediate vaccination rates in controlling epidemics [50]. Another study, using COVID-19 dynamics with demographics similar to the general US population, found that if no measures to limit social contacts are implemented, prioritizing older age groups is advisable when vaccine effectiveness or supply is low. However, when both are high, prioritizing younger age groups is more effective due to their role in driving transmission [9].
For diseases with IFRs varying by age group, incorporating both vaccination strategies and contact reduction measures simultaneously could be beneficial, although direct evidence is still limited. A modelling study comparing the propagation of COVID-19 and influenza using demographic data from Brazil, Uganda and Germany found that early and widespread vaccination, combined with measures to reduce social contacts, was crucial for reducing deaths. This approach proved to be more important than the specific vaccination strategy employed [10].
While our study primarily examines the Spanish population, its analysis readily extends to diverse countries or regions. Critical vaccination rates, triggering a shift in age-group prioritization, depend on demographic distribution and contact habits. Prior findings [23] suggest that countries or regions with younger populations or increased contact with average-age groups may benefit from vaccinating younger groups at lower rates. It is crucial to emphasize that any quantitative advantage hinges on demographics, vaccination rates, vaccine effectiveness, disease-specific IFR and additional epidemiological variables. Consequently, no absolute recommendation can be made for optimal vaccination protocols, and each case needs a tailored analysis considering its unique characteristics, as has been already pointed out [10,11].
Under social contact restrictions, the benefit of prioritizing vaccination for groups with higher social contacts becomes more pronounced, highlighting the importance of considering protocols that focus on younger or more socially active groups. It is crucial to recognize that the value of these alternative vaccination orderings extends beyond just the reduction in deaths and infections. During severe lockdowns, vaccination at intermediate rates—achieved globally during the COVID-19 pandemic—can significantly slow the spread of the virus. This not only reduces the total number of infections but also flattens the curve of active cases, thereby easing the burden on healthcare systems. While a primary goal of COVID-19 vaccination is to minimize fatalities [51,52], reducing infections is equally important, particularly in preventing the emergence of new variants [53–55]. The optimal vaccination strategy differs depending on the desired outcome: for instance, different strategies minimize fatalities or infections. Although an approach minimizing infections may not produce the most desirable results in terms of public health for some diseases where some groups have high fatality rates under unrestricted social contact, it is effective for illnesses like flu, where children are often vaccinated first [23,47,56,57]. For any contagious disease, making informed decisions about vaccination prioritization is crucial, as it impacts both mortality and infection rates, as demonstrated by the models used in this study.
Ethics
This work did not require ethical approval from a human subject or animal welfare committee.
Data accessibility
Data and relevant code for this research work are stored in GitHub [42] and have been archived within the Zenodo repository [39]. Appropriate references to the different files are provided in the paper.
Supplementary material is available online [58].
Declaration of AI use
We have not used AI-assisted technologies in creating this article.
Authors’ contributions
I.A.-D.: conceptualization, data curation, formal analysis, investigation, methodology, resources, software, validation, visualization, writing—review and editing; G.R.-M.: conceptualization, data curation, formal analysis, investigation, methodology, resources, software, validation, visualization, writing—review and editing; S.A.: conceptualization, investigation, methodology, project administration, resources, supervision, writing—review and editing, funding acquisition, validation; S.M.: conceptualization, funding acquisition, investigation, methodology, project administration, resources, supervision, validation, writing—original draft, writing—review and editing.
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.
Conflict of interests
We declare we have no competing interests.
Funding
This research was supported through grants PID2020-113284GB-C21 (I.A.-D. and S.M.), PID2019-109320GB-100, (BADS, G.R.-M. and S.A.) and PID2022-142185NB-C21 (PGE, S.A.), funded by MICIU/AEI/10.13039/501100011033 and by ‘ERDF/EU’ (PGE). MICIU/AEI/10.13039/501100011033 has also funded the ‘Severo Ochoa’ Centers of Excellence to CNB, CEX2023-001386-S and the special grant PIE2020-20E079 (CNB) entitled ‘Development of protection strategies against SARS-CoV-2’.
Acknowledgements
The authors are indebted to A. Vespignani for support with the use of their data.
Appendix A
A.1. Calculation of contact matrices
As introduced in the main text, our main objective was to study the effects of vaccination by comparing a ‘normal’ situation, where contacts are not limited, with a ‘confinement’ scenario in which contact restrictions exist. Data on contact reduction under lockdown for Spain were not available, however, so we resorted to two independent studies for The Netherlands [21] and England [22] to estimate a realistic matrix under contact reduction. We used both studies as a guideline in our particular case in Spain, as they compared the reduced matrix to a baseline, pre-pandemic situation. However, data cannot be directly applied to our case. We need to process their original contact matrices in order to meet the age grouping used in our study. Reduced contact matrices can be found in the file SuppFile_ConstrainedMatrix.xlsx, available at [39].
A.2. Dutch population study (NLD_Backer2021)
We use two of their contact matrices corresponding to data collected from February 2016 to October 2017 (Baseline) and April 2020 (Lockdown), respectively. Both matrices consider the same age group, so they are comparable with each other. However, one of our original 10 year groups is missing, G1 [0,10) years, which is divided into two smaller subgroups, [0,4] and [5,9]. We need to group all the information of their contacts, either among themselves or with the rest of the groups, in order to compare it later with the original matrix without contact restrictions [20]. Considering the structure of the matrix, correspond to the contacts that G1 individuals have towards the other groups. We calculate these values considering the smaller subgroups of the study as follows:
where and are the contacts of the original subgroups with the group in the matrix of the study and , the number of participants in the survey for each subgroup.
In the particular case, , we re-calculate the contacts within the group:
which accounts for all contacts that individuals from groups [0,4] and [5,9] have with individuals from the other two groups. This is because both subgroups are now part of the same ‘major’ group [0,10).
The elements therefore represent the contacts that other groups have with the individuals of the new group [0,10). In this case, we just have to add the contacts with the original two subgroups:
We apply this transformation for both matrices in the study obtaining the corresponding matrices in the pre-pandemic (Baseline, ) and lockdown (April 20, ) situations; matrices are available in the file SuppFile_ConstrainedMatrix.xlsx at [39]. The green area in both matrices highlights the region that remains constant with respect to the original matrix. Once their matrices are updated to consider our 10 year age groups we calculate the reduction factor matrix, , as the ratio between lockdown and pre-pandemic contacts:
A.3. English population study (GBR_Gimma2022)
In this study, the transformation of the data is more complex because the age groups defined are different when comparing the baseline and lockdown matrices.
— | Baseline: The POLYMOD [38] matrix is , as we have 5 year intervals for individuals aged and a final group with contact information for individuals over . In this case, the transformation is simple: we apply the process described in equations (A 1), (A 2) and (A 3), considering 5 year age groups two by two, to compute the outcoming, internal and incoming contacts for each original 10 year group. | ||||
— | Lockdown: The distribution of age groups in this matrix is different compared to the baseline POLYMOD matrix. Four groups are defined for individuals aged less than : [0,4], [5,11], [12,17] and [18,29]. The problem here is that we cannot combine their contact information into any of the original groups [0,9], [10,19] or [20,29] as they are not well defined by intervals. To solve this issue, we have to define some intermediate age groups in order to redistribute contacts: the [5,11] group is divided into [5,9] and [10,11] and the [18,29] into [18,19] and [20–29] (see Lockdown1_inter matrix in SuppFile_ConstrainedMatrix.xlsx, available at [39]). Since we lack specific data on the number of individuals of these intermediates, we assume that all ages were equally represented within the groups in the survey. Therefore, their contact contributions are divided proportionally. Considering [5,11] as an example, the intermediate group [5,9] takes the information of 5 years within the original group, so it takes of their contacts. Analogous to that, the [10,11] subgroup takes the remaining . Following this process, the two intermediates, [18,19] and [20,29], receive and of the [18,29] contacts, respectively. In the end, we recover the original 10 year age groups by adding the corresponding contact information as previously indicated: G1 [0–10)—[0,4] and [5,9]; G2 [10–20)—[10–11], [12-17] and [18,19]; and G3 [20,30) (directly obtained from the second split). |
In both matrices, the oldest age group comprises individuals aged 70 years and older. Thus, our two original older age groups, G8 [70,80) and G9 80+, are within it. Since we lack further information on the group over 80 years of age, we assume that the reduction in contacts is the same. Again, we apply equation (A 4) to obtain the reduction factor matrix.
A.4. Spanish matrix reduction
Finally, we average the two reduction factors ( elements) obtained from both studies to multiply the original elements in the pre-pandemic situation in Spain (Unconstrained SIYRD_9G, in SuppFile_ConstrainedMatrix.xlsx, available at [39]). Thus, we obtain the reduced contact matrix under lockdown in Spain as a case example of contacts under restrictions (Constrained SIYRD_9G in SuppFile_ConstrainedMatrix.xlsx). As expected, there is a significant average reduction of in the number of daily contacts in Spain. As table 3 illustrates, there is a certain degree of variability between age groups.
group | G1 | G2 | G3 | G4 | G5 | G6 | G7 | G8 | G9 |
---|---|---|---|---|---|---|---|---|---|
reduction (%) | 59.1 | 66.6 | 73.0 | 66.8 | 66.0 | 63.7 | 68.0 | 68.8 | 70.0 |