A number of projects for the development of a vaccine for the SARS-cov2 virus are in advanced stages of development and clinical trials [24, 8, 19, 13, 16, 2, 11] , and it seems possible that by the first or second trimester of 2021, a tested, off-the-lab vaccine could be available. However, off-the-lab availability is not enough: the vaccine will have to be produced and distributed in great quantities, and this will take time. The dire reality of world power makes it quite likely that the early production will the distributed mainly in the richest and most powerful countries, while many countries will initially receive an amount insufficient to vaccinate more than a fraction of the population. This state of affairs poses the problem of designing a vaccination strategy: who should be vaccinated first so that the limited amount of vaccine will have a maximum impact, at least while the country is waiting for more substantial amounts to be delivered? A general consensus is that medical and paramedical personnel, crucial in an epidemic, should be the first to receive the vaccine. However, assuming that once the essential personnel has been vaccinated there are still doses available, one faces the problem of how to distribute this surplus among the general population to minimize the number of casualties while the country waits for more vaccine to arrive.
We develop here an infection spread model suitable for studying vaccination strategies in the presence of a vulnerable portion of the population, and we use it to begin a study of vaccination schedules under various hypotheses regarding the constraints that the limited availability of vaccine or the limited vaccination capacity may impose.
The most common models in the study of epidemic diffusion are derivations of the standard SIR model [10, 1, 25, 23] , which makes strong hypotheses about the distribution of the population and its interactions. In all these models, if S(t) and I(t) are the numbers of susceptible (healthy) and infected people at time t, the probability that a healthy person enters in contact with an infected one is proportional to S(t)I(t). Because of this, the model doesn't account for the different sizes of social circles that different people may have [6, 7, 9] , thereby not allowing the modeling of targeted vaccination strategies.
We define a richer model of social interaction, one that takes into account the variability of social circles, and use it to simulate and evaluate various vaccination strategies. We consider a general model and ground it using the social and epidemiological data from Spain as a test-bed 1 .
The model divides the population in two groups: a more vulnerable group (conventionally referred to here as the elderly) and a less vulnerable one (the young). The parameters for the two groups have been set considering the Spanish population age 20-64 for the young, and age ≥65 for the elderly. The social connections among the young have been estimated using average social circle data [26] . The social connections within the group of elderly, as well as between elderly and young, are parameters that we vary in our tests.
Using this model, we test various scenarios that differ in the fraction of vaccine that is administered to the elderly as well as in the policy that we use to decide which young should be vaccinated first. We also consider the case in which the vaccination capacity is limited and the vaccine must be administered over an extended period of time, not negligible with respect to the spread of the epidemic. In this case, we use a genetic optimization algorithm to design an optimal vaccination schedule.
Social groups are customarily modeled as random graphs in which nodes represent individuals and edges represent social interactions between individuals [6, 20, 21] . Let G = (V, E) be such a graph, with V = {1, . . . , N } and E ⊆ V × V . We assume that the graph is undirected without self-loops ((u, v) ∈ E ⇔ (v, u) ∈ E and (u, u) ∈ E for all u, v), and let m = |E|/2 be the number of edges 2 . The neighbor of a node is
1 The epidemiological data are derived from the daily reports from the Ministerio de Sanidad, published in the web page www.mscbs.gob.es/en/profesionales/saludPublica/ccayes/alertasActual/nCov/home.htm 2 The factor 2 is due to the fact that each arc between u and v is represented twice in E, once as (u, v) and once as (v, u). and its degree is d(u) = |N (u)|. We indicate withd = k d(k)/N the average number of neighbors of nodes in the graph. Let P k = P u d(u) = k the probability that a random node, selected with uniform probability, had degree k. One crucial observation in most social networks is the existence of relatively few people with many contacts, and a "long tail" of many people with less and less social contacts [5, 27, 14, 22] . Specifically, many a social network exhibit a power law or scale-free distribution of the type P k ∼ k −γ where γ > 1 is an exponent that determines how "fat" the log tail is, and that depends on the specific social network [18, 4] .
Several random graph models have been defined to mimic this behavior. The one that we use here was described in [5, 15] . The algorithm creates the graph by adding one node at the time, and is controlled by two parameters, N and q, where N is the number of node of the final graph and q a parameter that determines connectivity. The graph is built in N steps, indexed by a parameter t = 1, . . . , N . At step t, a new node is created, node number t. At this point the graph already has t − 1 nodes. Let d(t, i), i = 1, . . . , t − 1 be the number of neighbors of node i at step t. From the newly created node t, q edges are drawn using a preferential attachment (or "rich gets richer") strategy: node t is connected to node i with probability
The algorithm that generates the graph is shown in Figure 1 . The function prefRnd(V, L) generates a random element of the set V using the list L = [d(1), . . . , d(t − 1)] to determine the probabilities of generating elements of V as in (2) . It can be shown that this algorithm generates a power law graph with P k ∼ k −γ , where the exponent γ depends on the parameter q [15] .
(V, E) = random network(n, q) 1. Figure 2 : Algorithm that joins two graphs through preferential attachment. It assumes that V y ∩ V e = ∅, that is, the identifiers of the nodes for the elderly and those for the young are different, Both V y and V e are ordered list, allowing to maintain the correspondence between the nodes and the values in d. The function unif(V ) returns and element of V chosen randomly with uniform probability.
The complete model is composed of two such graphs interconnected through preferential attachment. First, two power-law graphs are generated. The first is a graph G y = (V y , E y ) of young people generated with parameters N y and q y (|V y | = N y ), resulting in a graph with exponent γ y andd y average neighbors per node. The second graph, of elderly people, G e = (V e , E e ) is generated with parameter N e and q e (|V e | = N e ), resulting in a graph with exponent γ e andd e average neighbors per node.
Finally, the two graphs are connected: a numberd ey is selected, which will be the average number of young relations that an elderly person has. Then N e ·d ey edges are generated between the two graphs, selecting an elderly person at random with uniform probability, a young person using preferential attachment, and connecting them. The algorithm that joins the two graphs is shown in Figure 2 The rationale for using preferential attachment here is that certain young people (e.g., social workers, cleaning personnel, delivery people, etc.) have contact with many elders, while the majority of young people have little contact with them. The procedure allows us to create a model in which the interrelations within the groups and the relations between groups can be independently controlled.
Each node u of the graph represents an individual which, at any given time t, can be in one of four states, stored in the parameter state t [u]: susceptible (S), infected (I), recovered (R), or victim (V ). People transition from one state to another depending on their contacts and their recovery probability according to the following diagram:
The only transition that depends on the interaction between nodes is S → I, whereby a susceptible individual becomes infected. The probability of this transition depends on the parameter β as well as on the number of infected neighbors of a given individual. The model assumes that an individual interacts with their neighbors uniformly 3 . So, given a node u in state state t−1 [u] = S at time t − 1, a node v ∈ N (u) is picked randomly with uniform probability. If state t−1 [v] = I then, with probability β the node goes to state t [u] = I and, with probability 1 − β, it remains in state t [u] = S. The other states evolve independently of the state of the neighbors of u, and do so according to the schema of Table 1 . These transitions refer to non-vaccinated individuals. A separate indicator for each individual indicates whether the individual was vaccinated. Individuals can be vaccinated only if they are in state S (viz., only healthy individuals are vaccinated), and once vaccinated they will permanently stay in state S (viz., the vaccine is 100% effective). of the graph, while epidemiological parameters (those of (3) for young and elderly, that is, β y , φ y , ρ y , y , β e , φ e , ρ e , e ) determine the spread of the infection in each neighborhood of the graph. The epidemic parameters are not observable, but they can be derived based on observable values, as shown in Section 3.2. The observables used here are relative to the early spread of covid-19 in Spain, specifically the period up to March 25, when the effects of the lock-down decreed by the government on March 14 began to have a measurable effect effect on the spread of the epidemic.
The tests use graphs with N = 5, 000 nodes, with 80% of the population being young (N y = 4, 000) and 20% being elderly (N e = 1, 000). We executed preliminary tests with various values of N , observing that for N > 1, 000 the results were stable and their qualitative characteristics did not change. The social graph of the young is fixed, generated with q = 20, resulting in γ y = 1.77 andd y = 16.9. These values were chosen to be in compliance with the model in [26] . Then, four different graphs are generated, with different values of the social connection between the elderly (d e ) and of the connection between young and elderly (d ey ) 4 . These data result in four graphs, identified by the codes in Figure 2 . Graphs of type A A A are characterized by sparse relations between the elderly, that is, each elderly person interacts with few other elderly people, while in graph of type B B B there are many more elderly-to-elderly interactions. Graphs of type α α α have few cross-age interactions: each elderly person is in contact, on average, with only a few young people, while β β β graphs have denser interaction. These distinctions will allow us to analyze the effects of different vaccination strategies when used in conjunction with other sanitary measures such as the (relative) isolation of categories at risk or the reduction of contacts between these categories and the general population.
In order to determine the hidden epidemiological parameters φ, ρ, and from observable data, we use an SIR model as a first approximation. At this time, very few cases of reinfection have been reported [3, 12] so, for the time being, one can assume that the immunity is permanent (at least in the time span considered in the simulation), and set y = e = 0.
Assume a population of I(t) infected people maintained in isolation. In the one-compartment SIR model, their evolution is described by the equations
The solution for I(t) is given by
where I 0 is the initial number of infected individuals. The average time a person stays infected before they either die or recover is 1/(φ + ρ). Assuming the day as the time unit, and an average D days duration of the disease, independently of age, one obtains:
The values of the individual parameters are determined using the measured letality, that is, the number of victims divided by the number of cases. Let L y and L o be the letality for young and old people, respectively. From (5) and the first of (4) we have
For large t,
Therefore the letality in terms of the φs and the ρs can be expressed as
These equations, together with (4) yield
The value of β is derived from its relation to the unitary infection factor R 0 , which is given by Table 3 : Parameter values used in our simulation. The observable parameters in the first part are derived from the epidemic data for Spain in the period prior to March 30, before the lock-down showed its effects. The unobservable parameters in the second part are derived from those, and are those used in the model.
The values used in the simulations of Section 4 are given in Table 3 .
In the first set of simulations, we assume that all vaccination is done on a healthy population before the onset of the epidemic. In the case of an ongoing epidemic such as the one of covid-19, this corresponds to a scenario in which only the still healthy population is considered, keeping it effectively isolated from the people who were already infected. We simulate the different vaccination scenarios given by the parameters specified below and then, once the suitable number of people has been vaccinated, we simulate the unfolding of the epidemic, keeping track of the victims as a function of time. The simulations are characterized by three parameters: The parameters are chosen so that F · ψ is in any case no greater than the fraction of elderly people in the population. If F · ψ is greater than this value then, once all the elderly have been vaccinated, the notion of "vaccination strategy" would cease to be useful: all the vaccine that is left after vaccinating all the elderly would be administered to the young.
mode mode mode: The criterion with which the vaccine is administered to the young. The fraction of the vaccine allotted to the elderly is always distributed uniformly inside the group. The fraction that is used to vaccinate the young is allocated according to one of the following modes: popular popular popular: the vaccine is administered in order of "popularity", that is, the nodes in the graphs with the most neighbors are vaccinated first;
connected connected connected: the vaccine is administered in order of connections with the elderly: the young people with the most connections to the elderly are vaccinated first.
The rationale for the popular mode is that by vaccinating the people with the most connections will slow down the spread of the infection by removing dangerous "hubs" 5 . The rationale for the connected mode is to protect the vulnerable population by vaccinating the people that can spread the infection from the highly socialized young people (where the spread is rapid) to the less socialized but vulnerable elderly.
The duration of the simulations is fixed and is set to T = 100 days.
Mode: popular popular popular; vaccination level: 5% 5% 5% (F = 0.05) The second series of simulations assumes a limited capacity to administer the vaccine, so that at most p people can be vaccinated on any given day. The total amount of vaccine available in the given period is still assumed to be sufficient to vaccinate a fraction F of the population (here F = 0.2), so the daily vaccination capacity is p = F/T (expressed as the fraction of the population that is vaccinated each day).
In this case the controlled variable is a list ψ[t], t = 1, . . . , T , which determines, for each day, the fraction of daily vaccination capacity that is administered to the elderly. That is, on day t, p · ψ[t] elderly people and p · (1 − ψ[t]) young people are vaccinated (both these values are expressed as a fraction of the total population). As in the previous case, the elderly to be vaccinated are chosen at random with uniform probability, while the young are chosen according to the selected mode. For T = 100, ψ is a list of 100 values, so it clearly unfeasible and uninformative to sample it and present all the possible sequences in a graph. Instead, we use an optimization algorithm to determine the sequence ψ[t] that minimizes the number of victims during the vaccination period. The optimization is done using a genetic algorithm whose details are given in Appendix A.
In this section, the results are presented in the hypothesis that the delivery capacity is infinite, the only limitation being the amount of vaccine available. The three parameters that affect the results are those in Section 4. On day 0, a number of young and elderly people are vaccinated and from that point on we track the number of fatality for T = 100 days. The infection begins at t = 0 with one randomly chosen young person infected (the infection never starts with the elderly). Figures 3 and 4 show the data in the hypothesis that the amount of vaccine available is enough to vaccinate 5% of the population. In Figure 3 , the young have been vaccinated using the popular mode, while in Figure 4 the connected mode was used.
Each series has four graphs, one for each of the graph types of Table 3 , and the four curves in each graph are the victims for different values of ψ, the percentage of doses that are given to the elderly (ψ = 0.25, 0.5, 0.75, 1). In almost all cases, at the end of the 100 days, the result is the expected: the number of victims decreases as the fraction of vaccine administered to the elderly increases. The only exception is Figure 4 .a, in which the differences are not statistically significant.
The smallest number of victims is achieved in Figure 4 .a. This is the graph Aα, with little contact among the elderly and little contact between the elderly and the young. The value is lower in Figure 4 .a, the connected mode, suggesting that if the vulnerable population is relatively isolated, vaccinating the (relatively few) contacts is a good strategy.
Mode: connected connected connected; vaccination level: 10% 10% 10% (F = 0.1) A peculiar phenomenon, which we call risk inversion takes place in the popular mode, especially in the Aβ and Bβ graphs: if all doses are given to the elderly (ψ = 10) the number of victim in the first 20-30 days is higher than if some of the doses are given to young people. In this case, vaccinating only the elderly does offer some protection to the most vulnerable (albeit an incomplete one: there are, in the case F = 0.05, not enough doses to vaccinate all the elderly), but leaves the more heavily connected young completely exposed, leading to a rapid expansion of the epidemic. One the one hand this leads, in the β scenarios in which young and old are more connected, to more victims among the unprotected elderly and, on the other hand, to relatively many victims among the young. On a longer time span, the lack of protection of the elderly for low ψ compensates the effects of slow spreading, and the victims increase when ψ is small. Figure 6 .a with respect to Figure 4 .a, which suggests that relative social distancing is key to obtain the best results in a situation of scarcity of vaccine. The risk inversion is present here as in the previous cases.
One striking result is that of Figure 6 .a and 6.b, in which the mortality behaves opposite as in other cases: giving only 25% of the doses to the elderly results in a significant decrease in mortality. These figures correspond to the connected mode of the α graphs, with scarce connection between the elderly and the young. In this situation of high isolation and relatively high availability of doses, using 75% of the vaccine for the connected young effectively isolates the elderly while reducing the victims among the young. Compare this with Figures 8.a and 8 .b. In this case, the number of doses is sufficient to provide a better isolation to the vulnerable group, so that in the α graphs the number of victims is low regardless the value of ψ. In the β graph of Figure 8 .b, the stronger connection between elderly and young balances out the higher number of doses available, and the phenomenon of Figure 6 .a is repeated. The convenience of targeting the "connection points" (the young people that have connections with the vulnerable population) depends on the amount of vaccine available and on the density of the relations between elderly and young.
The risk inversion in the early days of vaccination is present in all cases, although not with the same prominence. This phenomenon suggests that if the vaccine is available in batches, the vaccination policy might be different for different batches: a first vaccination in part to the young to limit the speed of the infection and, if a second batch is available within 20 or 30 days, a more massive vaccination to the elderly to protect them. Figures 9 and 10 show the optimal vaccination schedule as determined by the optimization algorithm. The vaccination schedule is assumed to extend over a period of 100 days, and the total number of doses is assumed to be sufficient to vaccinate 20% of the population. The results were obtained using 500 generations of 800 genes each (see Appendix A). The results of the best genes of several generations were checked to verify that the solutions were stable. In the figures, the black crosses mark the output of the optimizer, that is, they represent the optimal fraction of vaccine to be administered to the elderly each day. The red line is a polynomial least squares approximation of these data using a polynomial of degree 10, which allows us to detect the general trends of the solution without the abrupt changes: our considerations will be based on these curves. The blue and green lines are the victims among the young and old people respectively, expressed as a fraction of the respective population. In Figure 9 the vaccination of the young is done using the popular mode, while in Figure 10 the connected mode was used. The four graphs in each figure are relative to the four graph types of Table 3 .
Mode: popular popular popular One common feature of all the strategies is the high initial fraction of vaccine administered to the elderly followed by an increase in vaccine given to the young when the number of casualties begins to grow. The peaks are more pronounced in the popular mode: in this case, after an initial short vaccination period to protect part of the vulnerable population, the focus moves to the young in order to stop the spread of infection by cutting the social hubs through vaccination. The situation is less clear in the β graphs (with many connections between the two groups). In this case, especially in the Aβ graph, the optimal strategy seems to involve a more balanced approach, with alternating focus on the elderly and the young.
For the connected mode, in the α graphs, the peak of vaccination of the elderly is wider. Given the low value ofd ey , relatively few vaccinations of the young are sufficient to block the spread of the infection to the vulnerable population, so the initial vaccination effort can concentrate on the vulnerable for a longer period of time. In the Aβ graph, with few elderly-to-elderly relations and strong connections with the young, the focus shifts on a much earlier vaccination of the young.
Mode: connected connected connected 
We have developed a model that, with respect to the standard SIR model, allows us to take into account the variability of the social relations of people in a social scenario. Given the existence of social groups especially vulnerable to the disease (herein referred to as "the elderly"), two factors appear crucial when defining a vaccination strategy: the social relations among the elderly, and the connections between the elderly and the young. The most general result that we can draw from the model is that a prudent vaccination strategy should balance two requirements: protect the elderly, and vaccinate enough young people (with intense social relations) to slow down the spread of the infection. The exact nature of this balance depends on the social connections and on who, among the young people, is vaccinated. Overall, the most effective strategy seems to be to give priority, among the young, to those who are more heavily connected to the vulnerable people (caregivers, cleaning personnel, etc.). If the amount of available vaccine is not too scarce (at least 10% of the population in our simulations), the results are good even if only 25% of the doses is administered to the vulnerable. If additional measures are taken to keep the vulnerable isolated, the 25% solution is actually the one that gives the best results. If the elderly are not isolated, then the strategy of splitting the vaccine between them and the young has positive results in a short (20-30 days) period but, in a longer period, it is necessary to move to a heavier vaccination to the vulnerable.
A more detailed strategy can be devised if we take into account a limited delivery capacity and the distribution of the vaccination over the whole duration of the epidemic. In this case, the optimal strategy seems to begin with a higher fraction of elderly vaccinated then, when the number of victims begins to grow, with a more intense vaccination of young people to keep the spread under control. The strategy that follows this first phase is less clear, but it appear to always consist in an oscillation with periods of more intense vaccination of the elderly alternated with periods of more intense vaccination of the young, all with a general tendency to increase the fraction of vaccine administered to the young.
A few words of caution are in order. All these results were obtained through a mathematical model in which, necessarily, many of the details of the actual situation are abstracted. While the general conclusions that we derived have undoubtedly a validity, applying them in the field requires extreme caution and a deeper analysis of the evolution of the epidemic. Also, our model has assumed that the "most popular" and "most connected" people can be identified and summoned for vaccination. This hypothesis might be reasonably realistic for the connected mode, in which the focus can often be identified by profession than with the more elusive popular mode. The distribution of the number of neighbors of the node one arrives to when crossing a random edge [21] may help implementing this mode.
Finally, we have ignored here all the problems of public perception, and the media backlash that may be caused by a strategy that, rational as it might be, does not follow what the people (or, more importantly, the media) perceive as "common sense". used as a vaccination schedule (fraction of the vaccination that is given to the elderly each day). The value of G is not a trivial function of ψ. In particular, we cannot assume any of the standard properties (continuity, derivability, analiticity, etc.) that most optimization algorithms assume. Because of the nature of the function, a genetic algorithm [17] is a reasonable option to solve this optimization problem.
We transform the problem into a discrete one by restricting ψ[t] to values of the type ψ[t] = 1 n[t] + 1 (13) with n[t] ∈ [0, 15]. Each n[t] can be represented as a four-bit number, so the sequence n[1], . . . , n[T ] that determines the solution can be represented as an integer of 4 · T bits. This is the "gene" of our algorithm. A generation is composed of a collection of G genes γ i , i = 1, . . . , G. Given a gene, we derive from it the vaccination schedule using (13) and execute the simulation. This will result in a score s i . Standard genetic algorithms are function maximizers, so we consider as score the inverse of the number of fatalities at time T (s i = −V γi (T )). There are several methods to create the following generation of individuals. Since the performance of the algorithms seems to have little dependence in the specific method used, we use one of the simplest, based on the creation of an intermediate gene pool. The gene pool is a set P of P individuals possibly replicated (generally |P | = G: the pool has the same size as the generations) such that the number of "copies" of an individual in the pool is proportional to its score. An easy algorithm for generating a pool is the tournament: we do P comparisons of pairs of individuals taken at random from the generation: the individual with the highest score goes into the pool: P ← ∅ for k=1 to P do i ← rnd(1,G) j ← rnd(1,G) if s i ≥ s j then P ← P ∪ {γ i } else P ← P ∪ {γ j } fi od
In order to build the next generation, pairs of genes are taken at random from the pool (with uniform distribution) and crossed to create two new individuals that will go into the next generation (this requires that G be even). We use the method of the double cut to cross the genes. Two values a, b ∈ [0, 4T ] are chosen randomly. The two offspring are then generated as in Figure 11 , in which we assume a < b. We also define a small mutation probability: for each new gene, with a (small) probability p, we pick a random bit and flip it. Note that this method doesn't guarantee that the best individual of a generation will pass unchanged to the next, so we actually use the crossing to create G − 2 individuals to which we add the two best performers of the previous generation.