Multiple outbreaks in epidemic spreading with local vaccination and limited vaccines

How to prevent the spread of human diseases is a great challenge for the scientific community and so far there are many studies in which immunization strategies have been developed. However, these kind of strategies usually do not consider that medical institutes may have limited vaccine resources available. In this manuscript, we explore the Susceptible-Infected-Recovered (SIR) model with local dynamic vaccination, and considering limited vaccines. In this model, susceptibles in contact with an infected individual, are vaccinated -with probability $\omega$- and then get infected -with probability $\beta$. However, when the fraction of immunized individuals reaches a threshold $V_L$, the vaccination stops, after which only the infection is possible. In the steady state, besides the critical points $\beta_c$ and $\omega_c$ that separate a non-epidemic from an epidemic phase, we find for a range of $V_L$ another transition points, $\beta^*>\beta_c$ and $\omega^*<\omega_c$, which correspond to a novel discontinuous phase transition. This critical value separates a phase where the amount of vaccines is sufficient, from a phase where the disease is strong enough to exhaust all the vaccination units. For a disease with fixed $\beta$, the vaccination probability $\omega$ can be controlled in order to drastically reduce the number of infected individuals, using efficiently the available vaccines. Furthermore, the temporal evolution of the system close to $\beta^*$ or $\omega^*$, shows that after a peak of infection the system enters into a quasi-stationary state, with only a few infected cases. But if there are no more vaccines, these few infected individuals could originate a second outbreak, represented by a second peak of infection. This state of apparent calm, could be dangerous since it may lead to misleading conclusions and to an abandon of the strategies to control the disease.


I. INTRODUCTION
Human interactions have a structure that can be well described in the form of a complex network [1][2][3][4]. In the last few years, new technologies allowed us to record large amount of data of contact patterns [5][6][7][8]. This data has become accessible to researches that use data-driven network modeling approaches to analyze and understand spreading in social systems, for example, how epidemic and even rumors spread in real populations.
Scientists have focused [9][10][11][12][13], on modeling and analyzing disease spreading since it can lead to catastrophic health consequences as well as large economic losses. Several mathematical approaches have been developed and used to study different epidemic models, improving the understanding of disease spreading on complex networks [14,15] (and references therein).
Since one of the goals of health authorities is to minimize health catastrophes and economic impact of health policies, many studies have focused on establishing immunization and mitigation strategies for enhancing the functionality of a society and reduce the economic cost [16,17]. For example, vaccination programs [18] are very efficient in providing immunity to individuals and as a consequence, the final number of infected people decreases considerably. However, these strategies are usually very expensive and unrealistic, because vaccines against new strains are usually not available during the initial propagation stage.
As a consequence, non-pharmaceutical interventions are needed to protect the society. One of the most effective and studied strategy to slow an epidemic is quarantine. However, it has the disadvantage that full isolation has a negative impact on the economy of the region and it is difficult to implement it in a large population. Thus, it is important to find a balance between these two strategies. Another policy, such as social distancing strategies, have been modeled and implemented in order to reduce the average contact time between individuals [19,20]. This kind of strategies, usually include closing schools, cough etiquette, travel restrictions, intermittent connections, etc. Unfortunately, in most cases, these strategies do not prevent a pandemic, but only delays its spread.
One of the most remarkable cases of a disease spreading was the pandemic occurred by the H1N1 strain in 2009, which caused about 15.000 deaths. Initially, the disease propagated over the network of close contacts, and then through the airline network, transporting infected individuals to different cities, thus spreading the disease all over the world. One of the most used models to mimic these kind of epidemics is the Susceptible-Infected-Recovered (SIR) model [21][22][23][24]. In this model an individual can be in only one of three possible states: Susceptible (S), Infected (I) or Recovered (R). An individual in state S in contact with an I, changes to an I state with probability β. After a period of time t r the infected individual changes to an R state and stops transmitting the disease. This model presents, in the steady state, two regimes governed by an effective probability of contagion T = T β,tr = 1−(1−β) tr , such that for T ≤ T c the system is in an epidemic-free phase and for T > T c it is in an epidemic phase, where the disease reaches a high fraction of the population. The SIR model has been also successfully applied to model the case of SARS and other diseases of influenza type [25,26]. For decades, researchers have studied different scenarios of the SIR epidemic model [27][28][29][30] and develop mitigation strategies to prevent the epidemic [18,[31][32][33][34], such as isolation, quarantine and random and targeted vaccination. Particularly, Valdez et al. [35] studied, using the SIR model, the effect of an intermittent social distancing strategy on the propagation of epidemics in adaptive complex networks. Based on local information, a susceptible individual interrupts the contact with an infected individual with a certain probability and restores it after a fixed period of time. In a similar way, Ref. [36] extended the model and was able to successfully predict the date of extinction of the Ebola's outbreak in Liberia in 2014. Ebola outbreaks have been studied by the scientific community due to the high impact of this epidemic on certain regions of southwest Africa, mainly in Guinea, Sierra Leona and Liberia. Fortunately, there exists available data for the scientific community enabling to study more accurately the behavior and propagation of this disease [36,37].
During this Ebola outbreak, a vaccine trial has been performed and used in Guinea in 2015 in the capital city. It has been found that the strategy of the vaccine trial applied to mitigate the transmission of Ebola-Virus-Disease (EV D) in Guinea in 2015, was very efficient [38]. The vaccine trial tested the efficacy of an experimental vaccine against Ebola.
The trial used a "ring" vaccination strategy based on the approach that was used to eradicate the smallpox [39]. This involves the identification of a newly diagnosed Ebola case, and then the vaccination of all his contacts and the contacts of those contacts, which are usually their family members, neighbors, co-workers and friends. In practical terms, the close contacts of a newly identified Ebola case have been vaccinated, if they consent to it. This is the basis of the motivation of the present study where we immunize the neighbors of an infected individual in network models. Moreover, the amount of vaccines, sometime due to economic restrictions, may not be enough to protect all the susceptible neighbors of the infected individuals during the whole spreading process. Thus, we propose and study here models with the scenario of limited vaccines availability, i.e. not enough for all the vulnerable population. Generally speaking, the available resources to control, avoid, or maybe enhance a spreading process are limited, and many studies have focused in how to optimally use these scarce resources against a disease [40], or even to deal with illicit drug usage [41].
We are interested in testing how this limited amount of accessible vaccines affects the spread of the epidemic. Motivated by this, we present here a model of localized vaccination that mimics the vaccine trial in Guinea, in which only neighbors of infected individuals could be immunized using limited vaccination units. This could help to understand how the propagation of a disease is affected by the localized and limited vaccination. Therefore, the question we wish to answer is how to use, in an efficiently way, this limited amount of vaccines with the aim of reducing the propagation of an epidemic.
In our model, the state of an individual can be Susceptible (S), Vaccinated (V), Infected (I) or Recovered (R). A susceptible individual in contact with an infected one will become vaccinated with probability ω, until the total number of vaccines V L is used. If he does not become vaccinated, then with probability β he will become infected. Then, after a certain period of time t r , the infected individuals will recover.
Using an edge-based compartmental model and the generating functions theory [24,42], we obtain and study the evolution equations for the fraction of S, I, R and V individuals and find a perfect agreement between theoretical and simulation results. Then, we study the steady state of the epidemic process, for which there is no more infected individuals.
We find two different phases, an epidemic and a non-epidemic phase, separated by a critical threshold β c , which depend on ω. Below β c the disease can not spread and the fraction of recovered individuals R approaches to zero. On the other hand if β is fixed, there is a critical vaccination probability ω c , above which the epidemic can not develop, since the disease is successfully blocked by immunized individuals.
Above β c , depending on the parameters, we find either a continuous phase transition or a discontinuous phase transition for a second epidemic break at higher critical threshold β * , which depends on V L . Similarly for fixed β, below ω c we find a discontinuous transition for ω = ω * in the fraction of recovered individuals. In Refs. [43][44][45] the authors also found a discontinuous transition but in the density of infected individuals, and using an endemic epidemic model (susceptible-infected-susceptible), where the recovery of sick individuals depends on the availability of healing resources.
We also find, that depending on the parameters there are values of β = β † > β * and ω = ω † < ω * that characterize a crossover between two regimes, one for which the available amount of vaccines is sufficient to immunize the population during the whole process and the other in which it is not.

II. THE MODEL
At the initial state, all the individuals in the network are susceptible except one, the patient zero, which is in the infected compartment or state. Before spreading the disease, all the susceptible neighbors of this individual receive a vaccine with probability ω. Notice that this vaccination is local and dynamic, and is done through the links that connect infected to susceptible nodes. Then the patient zero will try to infect all its neighbors that have not been vaccinated, and this event will occur with probability (1 − ω)β. In the next time step of the process, all the susceptible neighbors of the infected nodes will be vaccinated with probability ω again or will be infected with probability (1 − ω)β. The infected individuals will move to the recovered state R after t r units of time since they become infected, and the vaccinated or immunized individuals will remain in state V. While the disease spreads through the population, the number of vaccinated people increases, until the health institutes run out of vaccines. We define V L as the fraction of available vaccines over the entire population.
When this limit is reached, no more individuals can be immunized and hence those in the infected state will infect their neighbors with probability β. In the steady state the epidemic is over when the fraction of infected individuals is zero, thus the individuals can only be in state S, R or V. For demonstration of the model see

III. THEORETICAL FORMALISM
The edge-based compartmental model (EBCM) [8,24], was applied to model the SIR and was adapted by Valdez et al. for discrete time, and a fixed recovery time t r [27]. We can solve theoretically the evolution and the steady state of this model with unlimited vaccines [46] and also adapt it here for the limited case. The EBCM is based on a generating function formalism, implemented in branching and percolation processes on complex networks.
This approach allows to study not only the steady state but also the temporal evolution of the process. First, we derive the general equations for the case of unlimited vaccines and then we explain the effect of the depletion of the vaccines. Denoting the fraction of susceptible, infected, vaccinated and recovered individuals at time t by S(t), I(t), V (t) and R(t), respectively, the EBCM approach lies on describing the evolution of the probability that a randomly chosen node is susceptible. In order to compute S(t), a link is randomly chosen and then a direction is given, in which the node in the target of the arrow is called the root node, and the base is its neighbor, called base node. We denote θ t to the probability that at time t, the base node does not transmit the disease to the root node and neither induces the immunization of the root node. In this approach, the state of the base node can not be affected by the root node, so that we can treat the state of the roots neighbors as independent [23,24,27]. A node remains as susceptible if none of its k neighbors cause its infection or immunization, then the fraction of individuals in the susceptible state at time t is given by where G 0 (x) = kmax k=k min P (k)x k is the generating function of the degree distribution, P (k), of the network [47]. To compute θ t we have to take into account all the possible states of the base node. Suppose an edge that connects the root node and the base node. Then, this edge has not been used yet to infect or vaccinate the root node if the base node is • in state S, with probability Φ S .
• infected but did not spread the disease to the root node, nor induced the immunization of the root node, which is expressed by Φ I .
• in state R but during the time it was infected, it did not propagate the disease to the root node, nor induced vaccination to the root node. This probability is denoted by • vaccinated or immunized, with probability Φ V . We summarize these probabilities in Table I.
In Fig. 2 a) we demonstrate the configurations of the root and the base node.
Thus, accounting all these cases (see Fig. 2 b), θ t is given by Similar to S(t) in Eq. (1), we can write an expression for Φ S (t) (see Fig. 2

b). The
neighbor of the root node, with degree k, is in state S if the disease does not spread through its k − 1 links and if none of its k − 1 neighbors, omitting the root node, do not induce its vaccination. Recall that the edge coming from the root node is not considered. Hence the probability that the base node is susceptible at time t is θ k−1 t and thus, where G 1 (x) = kmax k=k min kP (k)/ k x k−1 is the generating function of the excess degree distribution of the network and k is the average degree [48]. The evolution equations that describe the process for unlimited vaccination, i.e. V L = 1, are (see Appendix A for detailed derivation), . if none of its k partners induce its vaccination. We consider that the root node can not change the state of the base node, thus the latter node is susceptible if it does not get infected through its k − 1 links, and if none of its k − 1 neighbors other than the root, cause its immunization.
In the first equation, θ t decreases if the base node is infected at time t, and induces the vaccination of the root node with probability ω, or if it spreads the disease to the root, with probability (1 − ω)β. Note that Φ I takes into account that the base node and the root node had no prior interaction. The second equation represents the evolution of the probability that the base node is in state S, which is the finite difference of Φ S (see Eq. (3)).
The third equation is a bit more complicated. The root node has, with probability Φ I (t), an infected neighbor at time t that did not induce its immunization or caused its infection.
This probability changes if the infected base node causes the immunization of the root node or infects it, which is reflected in the first term. In the second term we have the susceptible individuals that become infected at time t and will be in state I in the next time step.
Unlike [27], where the authors used the EBCM to solve the classical SIR model, this term does not account all the variation of Φ S (t), since a fraction of the susceptible individuals go to state V. The fraction of the nodes in state S that go to state I is weighted with the factor and (1 − ω) tr (1 − β) tr is the probability that during the period t r an infected base node did not infect nor induced the vaccination to the root node. Finally, the last equation takes into account the immunized neighbors of the root node. The variation of Φ V (t) increases with time and is proportional to the negative change of Φ S . In this case the factor C ω = ω/ ω + (1 − ω)β , is the probability that the link between the root node and the base node is used to immunize. We explain these additional probabilities in Table II.
After computing θ t using Eqs.(4), we can compute the evolution of the fraction of susceptible, infected and vaccinated individuals at time t by, Notice than using these magnitudes we can compute the fraction of recovered individuals, The base node did not infect nor did it induce the vaccination of the root node change in the fraction of infected individuals is also proportional to the variation of the susceptible individuals, but here the factor C β is related to the transition from state S to I.
Hence, −C β ∆S(t) is the fraction of new infected individuals at time t. On the other hand, is the fraction of individuals that got infected t r temporal units earlier. Thus, this fraction represents the individuals that move to state R at time t, and hence contribute negatively to the fraction of infected individuals.
The set of equations (4) and (6) describes the temporal evolution of the process with unlimited vaccines (V L = 1). Now we assume that we have a limited amount of vaccination units, lower than the number of individuals in the system. Thus we impose a limit V L as the maximal fraction of vaccinated individuals.
The evolution of the system is the same as the unlimited case until V (t) reaches the vaccination limit, V L . At this point there is no available vaccines and thus the vaccination probability becomes zero, allowing the disease to spread without barriers. Hence the equations should be iterated normally until V (t) = V L , and then setting ω = 0 for the rest of the process. Nevertheless, since in Eq. (6), V (t) changes by finite increments, it is unlikely that V (t) matches with V L exactly. Thus it is not clear when iterating the equations, the precise moment at which the immunization process has to be stopped. We explain in Appendix B the procedure that has been performed to solve this problem and to reproduce exactly the results from the computational simulations.

IV. RESULTS
We perform stochastic simulations of the localized and limited vaccination SIR model over single networks with N = 10 6 nodes, whose degree distribution is Erdős Rényi (ER) with an average degree k = 10. The networks are built using the Molloy-Reed algorithm [49].

A. Temporal Evolution
To demonstrate the validity of the theoretical formalism, in Fig. 3 we show simulations and theory of the temporal evolution of the process for an infection probability β = 0.168, a vaccination probability ω = 0.45, a vaccination limit V L = 0.5, and a recovery time t r = 3.
The dashed lines are the theoretical results from the EBCM described in the previous section, while the solid lines represent different realizations of the stochastic simulations. We can see the excellent agreement between the theoretical equations (6) and the simulation results.
It can be seen from Fig. 3  the fraction of susceptible nodes shows a plateau, since as seen in Fig. 3(b), the epidemic almost vanished. and there are only few infected individuals that can infect the susceptible people. This plateau, which is also observed for the recovered individuals (Fig. 3(c)), seems to indicate that the system begins to stabilize, since the magnitudes change slowly with time. However, when suddenly the vaccines are exhausted, the fraction of infected individuals starts to increase again, reaching a second peak ( Fig. 3(b)) that could be even higher than the first one. This increase obviously cause a further decrease and increase of the susceptible and recovered individuals respectively. Finally the disease starts to fade away as the fraction of infected individuals decreases, then the system reaches the steady state and all the magnitudes stabilize.
From now on, in the rest of the manuscript, our results will be based only from the theoretical equations, since we find excellent agreement (Fig. 3) with the stochastic simulations.
Next we will show the temporal behavior of the process for different values of the infection probability, β. In the standard SIR model, without vaccination, there is a critical value β c below which there is no epidemic. This value satisfies the equality T c = 1/κ [50], where T = 1 −(1 −β) tr is the transmissibility, the effective probability of contagion, and κ = k 2 / k is the branching factor of the degree distribution of the network. k and k 2 are the first and second moments of the degree distribution respectively. From this relation, in which T c ≡ T (β c ), the critical infection probability can be obtained as β c = 1 − (κ − 2)/(κ − 1) [21]. In the present model of limited dynamical vaccination, the relation between T c and κ holds but in this case the transmissibility depends also on the vaccination probability ω (see Appendix D for the expression of T ). Closed expressions of β c for t r > 1 are quite complicated, however for t r = 1 the critical infection probability is simply β c = 1−(κ−2)/ (1−ω)(κ−1) . We see that this probability depends on ω but not on the vaccination limit V L . Besides β c , in our model there are also specific values of β associated with dramatic changes in the behavior of the magnitudes at the steady state. Unlike β c , these values depend on the number of immunization units. Next we will show how the magnitudes evolve with time when β is one of these specific values. In Figs. 4 (c) and (d), for larger values of β we observe that the two peaks start to get closer, until eventually they start to fuse together forming a single peak as in Fig. 4 (a), but much higher. On the other hand, for larger β values the curves become smoother since there is only a single outbreak.  [46]. It is important to point out that the existence of this maximum is highly influenced by the vaccination probability ω and the topology of the network. In Appendix E (Fig. 10) we show the fraction of vaccinated and recovered individuals at the steady state for networks with a heterogeneous power law degree distribution, and for different values of ω.
Next we observe what happens if we impose a limit on the fraction of available vaccines.
For V L = 0.5, represented by the red squares, we see in Fig. 5 (a) a plateau between two values of β. This happens since V can not surpass the vaccination limit. The lower of these values is β * , which we introduce in Fig. 3, and we call the other β † , which is greater than β * (see Fig. 5 (a)). Between these two values is the range of β for which V reaches its limit value, V L . Now we ask what is the effect of the vaccination limit on the fraction of recovered individuals. In Fig. 5  shows that there is no β † , and thus for any β > β * the fraction of infected is significantly higher, compared to the case of unlimited vaccines. The dashed lines that denote these points indicate the emergence of a second peak of infection, as can be seen in Figs. 4 (b) and (c).
In Fig. 5 (c) we show the fraction of susceptible individuals, which decrease with β and also show a discontinuous jump at β * , associated with V L .
In Fig. 5 (d) we compute the time it takes the process to reach the steady state, as seen in for β * = 0.2581 ( Fig. 4 (b)) the process takes much longer time to reach the steady state compared to β = 0.25809 (Fig. 4 (a)).
Thus, we can infer that similar to β c , the probability of infection β * denotes a transition point, that separates a region in which the immunization strategy stops the spreading of the disease, and another in which the vaccination units are insufficient to stop it. On the other hand, we do not observe a peak at β † , indicating that this is not a transition point.
Instead, this point denotes a crossover between the regime of insufficient vaccines and a regime in which the immunization units can not be used completely, since the probability of infection is too high. Moreover we recall that the existence of β † is related to the topology of the network, see Appendix E ( Fig. 10). To see also how the curves of Fig. 5 behave for a different recovery time t r , see Fig. 9 in appendix E where we show the steady state for t r = 3 and for V L = 0.4, the vaccination limit used in Fig. 3.
Next we fix the infection probability β and analyze how the magnitudes at the steady state change with the vaccination probability ω. In Fig. 6  This discontinuous jump is analogous to the behavior observed in Fig. 5 at β * , and one can easily relate ω † to β † .
Thus, for a disease with β fixed and for a fixed number of available vaccines V L , based on the vaccination rate we can predict the number of infected individuals in the system when the epidemic comes to an end. Furthermore and very importantly, for a given β and V L , we can chose the optimal rate of vaccination, ω, such that the fraction of infected be minimal or even zero.
The numerical values of β * , β † , ω * and ω † can be calculated theoretically using the generating functions formalism and branching theory [47,53] (see Appendix D for the derivation of the formula).
Finally in Fig. 7 we show the model phase diagrams, which exhibit different regions depending on the parameters. In Fig. 7 (a) In Fig. 7 (b) we fix β = 0.5 and show the phase diagram in the plane V L − ω. Here the solid and dashed curves represent w * and w † respectively. Here we see that when ω increases the depletion region becomes broader because more vaccines are applied. However when ω further increases, the immunization strategy gets more effective against the disease, and then a smaller amount of vaccines is required to control the epidemic. Furthermore when ω ≥ ω c , the disease can not propagate since all the paths are blocked by immunized individuals.
Using these phase diagrams we can learn how the regions of insufficient vaccines change with the available immunization resources in medical institutes, the infection probability, which depends on the disease, and the vaccination probability, which may depend on the on the vaccination limit V L . We also find that β * is a transition point, at which the curve of recovered individuals has a discontinuous jump, whose height depends on V L . This type of behavior, in which a discontinuous transition is observed, has been seen when the dynamics of propagation of epidemics is coupled with social processes [54,55]. Furthermore, we analyze the temporal evolution of the process close to β * . We find that for β β * , the temporal evolution of the fraction of infected individuals presents two peaks. When the disease is about to vanish the vaccines are exhausted, and then the infection probability β * is sufficiently large for an extremely small fraction of infected individuals to cause a sudden second outbreak. On the contrary, we observed that β † is not a transition point but a crossover, and that its existence depends on the topology of the network.
On the other hand we analyze the steady state of the process as a function of ω, finding other points of interest. One of them is ω † , below which the vaccination probability is too low to use all the vaccines, thus the immunization units are not exhausted but the epidemic is not effectively halted. Another point is ω * , above which the vaccination probability is high enough to control the epidemic with the available immunization resources, and shows a discontinuous transition in the fraction of recovered individuals. These results are of significant importance since the vaccination probability is one of the few parameters that can be controlled by the health institutes. Thus, ω can be chosen to minimize the number of infected individuals or even halt the epidemic in the primary stages, according to the available resources.
We solved the model using an EBCM, finding an excellent agreement with the stochastic simulations. Also, we used the branching theory to find the values of β * , β † , ω * and ω † .
In future studies we will analyze different features of the local vaccination model, as the immunization of first and second neighbors of infected individuals. Thus, intending imitate more accurately the ring vaccination strategy used against the Ebola outbreak in Guinea during 2015.
On the other hand, a node is in state V if is not susceptible and if was more likely to receive a vaccine rather than be infected, hence Similar to Eq.(A2) we can writeΦ Next we study the variation of Φ R , the probability that the base node is recovered and also, that during the time it was infected did not cause the infection or the vaccination of the root node. Since individuals recover with rate γ hencė To obtain Φ R first we have to rewrite Eq.(A1). The probability that the disease or the vaccination spread through at least one link to the root node is 1 − θ and thus Now combining Eqs. (A5) and (A6) Now integrating this equation and considering that 1 − θ and Φ R are negligible at the beginning of the process, then simply Combining Eqs. (2), (A3) and (A8) then and finally using Eq.(A1) we can write a single differential equation for θ , and also the rates r β and r ω become probabilities β and ω respectively, while γ is replaced by the recovery time t r .
Appendix B: Temporal evolution of the discrete-time equations close to the threshold V L When iterating Eqs. (4), one approach is to set ω = 0 in the temporal step that V (t) would surpass V L , ensuring that V (t) ≤ V L , but not that V (t) = V L at the steady state.
This would cause many fluctuations when computing V (t) as a function of β at the steady state. Thus, to reproduce exactly the results from the computational simulations another approach should be used. Next we detail the procedure we use to avoid this fluctuations.
At time n * we calculate ∆V (n * ), and if it turns out that V (n * ) = V (n * − 1) + ∆V (n * ) is greater than V L , then we calculate what is the value of ω * that satisfies V (n * ) = V L . Thus, instead of setting ω = 0, we use a smaller probability ω * < ω 0 , where ω 0 is the vaccination probability at the beginning of the process. This adjustment of ω may have to be performed a couple of times until finally ω = 0, but also V (t) = V L .
Thus in order deal with the limit of the vaccine units while iterating the equations, we have to use a vaccination probability that has a slight dependence on time. Furthermore, we have to take into account how this procedure affects Ω, since it depends on ω, as we can see in Eq. (5). Suppose that we choose a set of parameters for which at some point the population runs out of vaccines and we iterate the theoretical equations. Then, at the beginning of the process Ω(t) is given by Eq. (5), and when the process ends Ω = 1−(1−β) tr , which is simply the transmissibility of the SIR model [21]. Recall that Ω is the probability that a node in state I infects one of its neighbors or induces its vaccination during the time that this node remains in this state. Consider that during the time that a node is infected the probability of vaccination changes. Thus in this case the effective probability of infection or immunization Ω is lower than the one described in Eq. (5) but higher than the transmissibility of the SIR model. Considering the different probabilities of vaccination that may have to be used during the process we can write a general expression for Ω at time t, This expression takes in account that ω depends on time and is proved in detail in Appendix C. If ω t = ω for all t, then Eq. (B1) leads to Eq. (5): Appendix C: Derivation of Ω when ω depends on time Ω is the probability that a node infects one of its neighbors or induces its immunization during the time that it is infected, which is the recovery time t r . In the standard SIR model this probability is known as the transmissibility T and is calculated as follows: At n = 1 we simply consider the probability of infection β. At n = 2 we have to consider that the infection did not occur at n = 1, which happens with probability 1 − β. Next for n = 3, now we have to consider that there was no infection at n = 1 and n = 2, which happens with probability (1 − β) 2 , and so on.
For the vaccination model we have to include the immunization probability ω. For simplicity we define Ω 1 and Ω 2 as the effective probabilities of infection and immunization respectively, during t r units of time. Thus, similar to Eq.(C1) Note that if ω = 0, then Ω 1 is the same as Eq. (C1). Thus Ω 1 plays the role of the transmissibility in the dynamical vaccination model. Furthermore, we can write Eq. (C2) in a closed form and thus which is Eq. (5). When ω depends on time, a closed expression can not be found. Suppose that ω = ω t , and consider a node that was infected at t = 0, then Eq. (C2) takes the form .
which can be summarized in the following expressions (1 − ω j ).
Eq. (C6) represents the immunization and infection transmission for t = t r , since we are considering a single node infected at t = 0. We can generalize these equations for any time (1 − ω t−tr+j ), and finally adding Ω 1 (t) and Ω 1 (t) leads to Eq. (B1).
Appendix D: Computation of β * and β † In the steady state of our model, the fraction of vaccinated individuals when there is no limit in the number of vaccines is [46]: where C β is the same factor used in Eq. (4) and f ∞ satisfies the transcendental equation f ∞ is the probability that the branches of infection expand indefinitely, and T is the probability that an infected node spreads the disease through a link, also known as transmissibility [46] Note that this expression is the same as Ω 1 in Eq. (C3). If we compute the process for fixed ω and start increasing β, when β = β * the fraction of vaccinated nodes reach the limit V L . Then for greater values of β the number of vaccinated individuals in the steady state is equal to V L , but beyond β = β † if exists, the vaccination threshold is no more reached. If we observe closely Fig. 5, we see that for these values of β, the unlimited-vaccines curve equals the vaccination limit V L . Thus we can find β * and β † by solving the following system: where T β * ,ω is the transmissibility for β = β * and fixed ω. The same applies for β † and T β † ,ω . For a Poisson network G 0 (y) ≡ G 1 (y) and hence we can write a single transcendental equation to find β * or β † : This equation has one, two, or none solutions between 0 and 1 depending on the parameters. Two different solutions correspond to β * and β † while a single solution means that β † does not exist. Moreover if V L is large enough, there is no solution, which means that there are always available vaccines when needed.
Based on a similar reasoning we can find w * and w † for fixed β. From Eqs. (D3) we can derive a similar expression to Eq. (D4) for a Poisson network: which has two or none solutions depending on V L .
In Fig. 8  now we have two intersection points, which correspond to ω * and ω † . For V L = 0.606 the curve is tangential to the identity and hence there is a single solution, which means ω * = ω † .
As we can see beyond this point there is no solution. In Fig. 9 we show how the steady state changes when the recovery time t r is greater than 1. Note that this figure is similar to Fig. 5, with the same vaccination probability ω = 0.45, vaccination limits V L = 1 and V L = 0.5, and a different recovery time t r = 3. We see that β * , unlike β † which barely changes, is lower than for t r = 1. Furthermore, because the individuals are infected during a larger period of time, the fraction of recovered is larger and so is the discontinuous jump.
On the other hand in Fig. 10 (a) we show the fracion of vaccinated and recovered individuals in the steady state as a function of the infection probability β, for a network with a power law degree distribution, and for different vaccination probabilities ω. The vaccination limit and the recovery time are fixed, V L = 1 and t r = 1. We observe that for low values of ω the curves of vaccinated individuals exhibit a maximum while, for larger values they are monotonically increasing. As explained in the main text, β † exists as long as the curve of vaccinated individuals has a maximum, hence in this case for ω > 0.6 only β * exists. In addition, in Fig. 10 (b) we show the fraction of recovered individuals, which as expected decrease when the vaccination probability is higher.