key: cord-1020687-4lwzkfpf authors: Toledano, 'Oscar; Mula, Begona; Santalla, Silvia N.; Rodr'iguez-Laguna, Javier; G'alvez, 'Oscar title: Effects of confinement and vaccination on an epidemic outburst: a statistical mechanics approach date: 2021-02-22 journal: Physical review. E DOI: 10.1103/physreve.104.034310 sha: d78e44e5f0cd37f149285313378ce1f8276f97a3 doc_id: 1020687 cord_uid: 4lwzkfpf This work describes a simple agent model for the spread of an epidemic outburst, with special emphasis on mobility and geographical considerations, which we characterize via statistical mechanics and numerical simulations. As the mobility is decreased, a percolation phase transition is found separating a free-propagation phase in which the outburst spreads without finding spatial barriers and a localized phase in which the outburst dies off. Interestingly, the number of infected agents is subject to maximal fluctuations at the transition point, building upon the unpredictability of the evolution of an epidemic outburst. Our model also lends itself to testing vaccination schedules. Indeed, it has been suggested that if a vaccine is available but scarce it is convenient to carefully select the vaccination program to maximize the chances of halting the outburst. We discuss and evaluate several schemes, with special interest on how the percolation transition point can be shifted, allowing for higher mobility without epidemiological impact. Epidemic containment has been a crucial problem for humankind throughout our history, which has become of paramount importance since the COVID-19 outbreak late in 2019. The efforts to understand the spread of infectious diseases have attracted a large variety of professionals from all scientific fields, ranging from biology to sociology (see recent studies from different fields: [1] [2] [3] ). Mathematical modeling has also provided very relevant tools to analyze the stream of data concerning the infected population, and has substantially contributed to policy design (see, e.g., [4] [5] [6] ). Among the different approximations employed, agent-based models have been extensively used to provide policy recommendations even in the COVID-19 case [7] , for which simulations based on stratified population dynamics were carried out [8] . Multiscale approaches are known to improve our ability to explain the geographical expansion of the disease, and they have been also used, along with data-driven simulations, to analyze the epidemic of COVID-19 in Brazil [9] . Epidemic waves have also been considered [9, 10] , employing stratified population dynamics and non-autonomous dynamics, where mitigation effects are subsequently imposed and relaxed. During the COVID-19 epidemic and due to the scarce amount of vaccination doses in the first stages of the vaccination campaign, immunization schedules have also attracted attention of the modeling community with the aim of stifling the expansion of an epidemic burst by acting upon a number of individuals substantially below the percolation threshold [11] , e.g., through the search of certain types of motifs in the contact network. However, predictions about the evolution of an epidemic burst are known to be difficult and unreliable, because the uncertainty in the initial data propagates exponentially [12] . Yet, on occasions the inherently unpredictability of the models can be turned in our favor.As mentioned in previous works, the fluctuations in the number of infected people during the evolution of an epidemic burst can provide useful information regarding our ability to stifle an ongoing epidemic [13] . In this work, we propose a very simple agent model [14] of the susceptible-infected-recovered (SIR) type [15, 16] . Agents follow a random walk in the vicinity of their homes, which are randomly distributed on a square lattice, wandering up to a maximum distance r, which may vary from agent to agent. When r is very short, agents are effectively confined at home, and an infectious outbreak will very likely die off. As this distance is increased, outbreaks have larger chances of spreading throughout the system, until a percolation phase transition is reached, [17] in which the presence of an infinite cluster becomes certain. Further increases of mobility will have a very limited impact on the spread of the infectious disease. Interestingly, predictions about the future evolution of the epidemic are more difficult near the transition point. Indeed, above the percolation threshold the mean-field theory associated with the SIR equations provide a very accurate prediction of the evolution of the model, while below the transition point typical outbursts will not propagate beyond a certain correlation length. But as the agent mobility approaches to the percolation value, the precise geographical origin of the outbreak becomes crucial to predict the outcome. Even though an infinite connected cluster exists, the probability that the initial infected agent will be part of it is minimal at that point, leading to a maximal uncertainty. Thus, we show that the fluctuations in the number of infected people become a very useful observable in order to pinpoint the phase transition. Of course, our model is far too simple to be taken decisively in order to provide policy recommendations, which should be always undertaken with the assessment of experts from several fields. Yet, our results suggest that lock-down measures increase their effectiveness very sharply once the percolation threshold has been crossed. Finding the minimal lock-down measures that lead to the phase transition is therefore of paramount importance, and a very difficult task. We can not provide a complete recipe to establish the location of that threshold in practice, but we provide an interesting proxy as was noted before [13] : fluctuations in the infection reach. If an effective vaccine is available but scarce, a vaccination schedule becomes unavoidable, i.e.: the sanitary authorities should consider which agents must be immunized first. According to their social role, some of these agents can be selected in order to minimize the spread of the infection [18] [19] [20] . Of course, many individuals should receive the vaccine because of other considerations, such as age or health condition, and we will not discuss these very relevant aspects. Our model lends itself very easily to the evaluation of the efficiency of a vaccination program. In it, individuals are only distinguished through their relations: some homes are relatively isolated, so their inhabitants are not likely to spread the disease. Yet, a naive approach would consider that individuals with many connections should be the first candidates to be immunized. We will show that this criterion is not optimal. Indeed, some agents with few connections act as natural bridges between different clusters. Thus if they receive the vaccine the clusters will become isolated and the outburst will be halted. This article is organized as follows. Section II discusses in detail our agent model, and how our numerical simulations are performed. In Section III we describe our results regarding the percolation phase transition as the agent mobility is decreased, the secondary cases and the effects of a finite probability of recovery. Our theoretical framework is discussed in Section IV. The efficiency of several vaccination schemes is discussed in Section V, especially in connection with the shift of the percolation threshold. Our conclusions and ideas for future work are discussed in Section VI. We propose a simple model to characterize the effects of partial confinement during an epidemic expansion, which we call the confined-SIR model. Let us consider N agents moving on an L × L square lattice with periodic boundary conditions. Agent i possesses a home, determined by a fixed lattice point, H i . Agents can move freely within their wandering circles, centered at H i and with radius r i , which we will assume equal, r i = r for all i. Time advances in discrete steps ∆t = 1, in arbitrary units. At each time step, every agent performs a random step. If the step takes the agent out of the wandering circle a new step is attempted. Homes are distributed randomly over the whole lattice. Agents can be classified into three compartments: susceptible (S), infected (I) and recovered (R). Susceptible agents get infected with probability nβ when they share cell with n infected agents. Infected agents may get recovered at each time-step, with probability γ. In some cases, an effective SIR model can be written down and solved analytically. An illustration is provided in Fig. 1 . Let A i stand for the number of lattice points in the wandering circle of agent i, while A ij will stand for the number of lattice points in the intersection between the wandering circles of agents i and j. The probability that agent i and agent j will collide at time t can be estimated as Thus, if agent i ∈ S and j ∈ I, the probability per unit time that agent i will get infected by j becomes P I ij = βC ij . Thus, we can find the total probability per unit time that an agent i will get infected by just summing over all possible infection sources We will be mostly interested in the continuum limit, in which L → ∞, N → ∞ but the density ρ = N/L 2 is constant, in order to avoid lattice effects. In that limit, A i ≈ πr 2 and A ij corresponds to the overlap between two circles whose centers are a distance d ij apart, which is given by where u ij = d ij /(2r) ∈ [0, 1] is the dimensionless relative distance. Thus, we have showing that the effective infection probability does not depend on the wandering radius r. The dynamics of the model can be approximated by a Markov chain. Let us consider that, at time t, we know the epidemiological status of each agent, the location of their home H i , and the wandering radius r. The evolution of the expected number of infected agents can now be found, Notice that the C ij are constants, but the right-hand side evolves due to the migration of agents between sets I and S. In the totally mixing limit, r ∼ L and all agents wander around the whole region. Therefore, all C ij ∼ 1/L 2 . Thus, in a mean-field approximation Eq. (5) becomes and we obtain the usual SIR differential equations. Introducing fractional variables, S = S /N , I = I /N and R = R /N we reach: S ≈ −βρ SI, For early times, S ≈ 1 so the infection either expands or vanishes depending on the value of the basic reproductive number, which in this case can be defined as Our main focus will be the probability distribution for the total size of the outburst, characterized by the number of affected agents, A, that have suffered the infection at any moment in the long run, and its deviation, σ A . Of course, there are other relevant observables, such as the maximum number of infected agents, but they are not be considered in this work. Our model presents a very clear graph structure. Indeed, two agents i and j may infect each other if and only if C ij > 0. This condition determines an effective graph G, such as the one shown in Fig. 2 , in which the homes are represented as nodes and the links denote the pairs of agents that are able to meet. Both networks are obtained for N = 100 agents on 50×50 lattices, using r = 2 (top) and r = 5 (bottom). It can be readily seen that the network for r = 2 contains many disconnected clusters, while for r = 5 it contains a large cluster spanning most of the nodes. These networks are examples of random geometric graphs (RGGs) [21] [22] [23] [24] , in which N points are randomly scattered on a L × L square and joined with edges whenever their distance is less than r. RGG have proved to be a very convenient framework for the analysis of epidemic outbursts [25] , typically providing infected nodes with a certain probability per unit time of extending the infection to their neighbors. Yet, our confined-SIR model presents a substantial difference with those models. Indeed, the graph structure is equivalent to that of an RGG model, but the actual distances are still relevant in our case. Since our infected agents perform a random walk, their collisions into susceptible agents typically take place one by one, instead of through multiple infections. This difference has been shown to be relevant in the literature [26] . We have run numerical simulations of our confined-SIR model on a square lattice, placing N random homes and using the same value of the wandering radius r for all agents. Unless otherwise specified, we always average over N S = 5000 samples, with fixed values for β = 1 and γ = 0, N = 1000 agents, and ρ = 5 · 10 −6 . For large values of the wandering radius r we expect the epidemic burst to follow the mean-field approximation provided by the SIR equations for perfect mixing, Eq. (7). We have checked that conjecture numerically in Fig. 3 (top), where we show the probability distribution for the fractions S, I and R as a function of time, along with their expected values for the theoretical prediction. However, the results for smaller values of r are very different, as we can see in Fig. 3 (bottom) . Indeed, in this case the average values of S, I and R differ notably from the theoretical predictions, because the mixing assumption is inadequate in this regime. Thus, we are led to identify two different regimes. For large values of r we have nearly perfect mixing and for low r values the system is naturally divided into isolated clusters. Both are separated by a percolation transition, which we characterize in the next section. The theory of bond percolation has been one of the foremost paradigms of statistical mechanics for more than 40 years [27] [28] [29] [30] , with applications to magnetism [31] , wireless communications [32] , ecological competition [33] or sequence alignment in molecular biology [34] . Recently, a very relevant connection was described between the geodesics in strongly disordered networks and bond percolation [35] . Application of percolation theory to epidemics has been carried out by previous authors, such as those of Refs. [36, 37] , or more recently, Refs. [38, 39] . Bond percolation on a fixed lattice is characterized by a single parameter, p, the probability that a given bond will be present. Above a certain threshold value, p > p c , the probability that the system will contain an infinite connected cluster reaches one, with p c = 1/2 for the square lattice. Our system presents a strong similarity with bond percolation, but with another observable playing the role of the bond probability p. Let us consider N agents on an L × L lattice. The average distance between homes,r, can be estimated as Now, let us notice that the wandering radius r is only meaningful when compared to this average distance between homes. Thus, we introduce a mobility parameter, We will readily show that the mobility parameter ε is the only relevant variable to determine the geometry of our system. Critical points, such as the percolation transition, typically lead to large fluctuations, and can be obtained by considering the deviation of the number of affected agents, σ A . This magnitude is shown in Fig. 4 as a function of the mobility parameter for different lattice sizes (top) and densities (bottom). The validity of this strategy for characterizing the transition point has already been probed in the literature [13, 40] . Indeed, we can observe that the fluctuations in the size of the outburst present a maximum for a certain value of ε which only depends on the ratio between the recovery and infection probabilities, γ/β. Let us provide evidence that this critical value of the mobility parameter, ε c , corresponds to the position of the percolation phase transition. In the vicinity of the percolation phase-transition many observables show critical behavior in the form of power laws. The most salient of those is given by the average size of a cluster, s. Below the transition, we have with η = 43/18 ≈ 2.39 for a square lattice [17] . Figure 5 shows the average cluster size of our effective networks as a function of the mobility parameter ε, for different system sizes and populations. We can observe that, for all system sizes considered, the average size of the cluster diverges as we approach a critical value ε c , with an exponent which slightly differs from the value obtained in the square lattice, η ≈ 2.35, which seems to be robust under changes in the lattice size and the density. Epidemic phenomena are characterized by the number of secondary cases, defined as the number of susceptible agents infected by a single infected node. The average value of this number is related to R 0 , defined in Eq. (8) . Yet, the deviation of that number has also proved to be a valuable tool to understand the evolution of an epidemic outburst. Indeed, given the exponential nature of the expansion, large fluctuations in the spread rate will For comparison, we also show the connectivity histogram in the associated RGG and the fit to a Poisson law. naturally dominate the long-term evolution of the disease [41, 42] . It is thus relevant to wonder whether the spread of the infection within our model gives rise to large deviations in the number of secondary cases, which we will characterize using its full histogram. In RGG models, a proxy for that value is the number of connections for each node. Fig. 6 shows the full histogram for the number of neighbors of any given node, both before and after percolation, i.e., for r = 8 and 15 using N = 1000 and L = 447. The data fit perfectly to a Poisson distribution, corresponding to the theoretical expectation. Furthermore, the figure also shows the number of actual secondary cases, which is much more peaked, showing that the deviation is always low. The reason for this difference between RGG and our model lies in the statistical properties of time sequences of infection events involving a common infected agent. In an RGG setup infection times are uncorrelated Poisson variates, but in our model those times are not independent. Indeed, an infected agent may require a minimal time to infect two different neighbours, given by the relative positions of the overlapping areas between the three wandering circles involved. Let us discuss the effect of a non-zero recovery probability, γ > 0. In that case, infected agents may recover before they can propagate the disease, and thus we are led to compare two natural times: τ ij ≈ (βC ij ) −1 is the expected time before the infection may propagate from agent i to agent j (or viceversa), while τ R ≈ γ −1 corresponds to the expected time before recovery. Thus, if τ ij τ R we may assume that the infection will not be able to propagate from agent i to agent j, while in the opposite case, τ R τ ij , we may neglect the possibility of recovery. We may claim that a finite recovery probability provides an effective cutoff for the local infection probabilities, thus removing weak links from the graph. The top panel of Fig. 7 shows the average number of affected nodes, A , as a function of the mobility parameter as we increase the recovery probability per unit time, γ. It shows a continuous decrease, as expected. The bottom panel of Fig. 7 , on the other hand, shows the deviation in the number of affected nodes, σ A , showing that the fluctuations at the maximum do not depend strongly on γ. In addition, the value of the mobility parameter for which this maximum takes place grows rather fast with γ. Thus we are led to claim that the recovery probability per unit time strongly affects the value of the percolation threshold. Also, we can see in the bottom panel of Fig. 7 that for low γ the fluctuations show a peaked maximum, while it seems to reach a plateau for higher values of γ. Let us provide a theoretical framework in order to explain the nature of the phase transition observed in the system. A percolation transition of geometrical origin is combined with an epidemic spread transition that takes place when R 0 > 1 in a well mixed system. First, let us consider the percolation transition on an RGG of size L × L with N node and radius r. The probability that a single node will be isolated is 1 − π(2r/L) 2 [21] . The expected number of isolated nodes will be approximately given by The average connectivity of a node will be C = 4πε 2 , thus showing that the number of isolated nodes will be given by µ ∼ N e −C [23] . Interestingly, this magnitude drops quickly as r increases at the percolation transition, for large N . More rigorously, it has been proved that the longest edge of the minimal spanning tree (MST), , scales like for large N [22] . The longest edge of the MST is a good approximation for the minimal value of r for which a giant component will exist within our graph, thus signaling the percolation transition. The validity of Eq. (13) can be checked in the top panel of Fig. 8 , which shows our estimate for the critical radius r c when γ = 0 using several values of N , along with two different fits, which yields an exponent χ = 1.1 or the slightly more accurate fit to which yields χ = 0.98. The inset shows the results for the longest edge of the MST of graphs of the same type, along with both types of fits. Now, let us shift our attention to the characterization of the epidemic spread transition on generic networks [25, 43] . The standard approach to locate the spread transition on a random graph is to equate the expected number of infections and the expected number of recoveries in the unit time, which is usually expressed in the mean-field equation,β where k is the average degree of the network and β = Kβρ is the average infection rate for neighboring nodes, given by Eq. (4) . It has been shown that k can be replaced with the largest eigenvalue of the adjacency matrix Λ max , providing better results [44] . Even more, we may take into account the presence of correlations [45] . Yet, Eq. (16) is good enough as a first approximation. In our case, we will considerβ k to be the average value of the total infection rate, given by Eq. (4), but restricting ourselves to nodes belonging to the largest connected component. The probability that a node will belong to this giant component may be approximated as p C ≈ 1 − exp(−4πε 2 ), as we can readily check from Eq. (12) . Thus, we can write down a fundamental equation defining the critical mobility ε c , where γ 0 = Kρβ, and κ 1 and κ 2 are fitting parameters. Solving for ε c we obtain where γ 0 = Kβρ. Equation (18) shows that r c → ∞ as γ → γ 0 , showing that the infection will never spread beyond that point. The bottom panel of Fig. 8 shows the critical radius obtained for N = 1000 as a function of γ, along with a fit to Eq. (18), with κ 1 , κ 2 and γ 0 as fitting parameters. Interestingly, the fit is able to guess a suitable value of γ 0 ≈ 0.0008, despite the fact that the provided values were far from it. Let us consider the possibility of providing immunity through vaccination to a (small) fraction of the population, f v , for the case of zero recovery rate, γ = 0. In this section we will attempt to answer the following question: How do we select the agents that will receive the vaccine, if our only aim is to minimize the size of a future epidemic outburst? Notice that, in practice, many other issues must be considered in this situation, such as the health conditions or the age of the patients. The simplest vaccination schedule is merely to randomly select the individuals. Of course, we do not expect this schedule to be very efficient. We have considered several observables which can be employed to determine how useful it will be to provide the vaccine to a certain agent. The simplest one is the total degree, defined as the (14), while the slightly more accurate continuous line shows a fit to expression (15) . Bottom: Transition radius as a function of γ, for N = 1000, showing a fit to expression (18) . The transition value of r c /L is approximately equivalent to as shown in (13) . number of first neighbors in the effective graph. Naively, we might sort the agents by their degrees, and immunize the first f v N . Yet, it is more efficient to recompute the degree of each agent after each selection, and we will do so unless otherwise specified. The most promising observable is, nonetheless, the betweenness-centrality (BC) [46] associated to each agent, BC i , defined as the fraction of the total number of geodesics which go through agent i. This measure is global, and takes O(N 4 ) steps to compute using Dijkstra's algorithm to evaluate the geodesics, or worse when we recalculate the BC i after each vaccination [47] . In previous studies vaccination schemes based on immunizing the highest BC links have proved to bet very efficient (see e.g. [18] ). Considering the high computational cost of the BC, we are led to propose a cheaper alternative, the local betweenness-centrality (LBC) [48] , which is defined for each site as its BC corresponding to a subgraph restricted to itself and its nearest neighbors. In intuitive terms, the LBC is high for a node that connects neighbors which are otherwise disconnected among themselves, and is naturally related to the clustering coefficient. Thus we have considered the following five vaccination schedules. (0) No vaccination, considered the base case. (1) Select randomly. (2) Select the agents with highest degree (HD). (3) Select the agents with highest local betweennesscentrality (LBC). (4) Select the agents with highest betweennesscentrality (BC). The vaccination programs effectively change the topological properties of the network whenever the removed agents are not selected at random. Thus, in Fig. 9 we can observe a specific example of a network where different vaccination schemes have been performed over the same amount of agents. In all cases, let us focus on the cluster structure. For random vaccination, Fig. 9 (a), the large clusters remain untouched. For a degree-based vaccination scheme, Fig. 9 (b), we can see that links have been removed from the core of the clusters, but the clusters themselves remain connected. Indeed, immunizing the individuals with a large number of connections seems to have a low impact on the network structure in our case. The reason is that high-degree agents tend to be neighbors of other high-degree agents. Panel (c) of Fig. 9 , on the other hand, shows that removing agents with a large LBC breaks up some clusters, but some large clusters still remain active, leading to a likely propagation of an epidemic outburst to a substantial fraction of the population. Fig. 9 (d) shows the resulting network when the agents with a largest BC have been removed, and we can readily see that all large clusters have indeed disappeared. In order to describe the quantitative effect of the different vaccination schedules, the top panel of Fig. 10 depicts the expected value of the long-term number of affected agents as a function of the mobility parameter. The no-vaccination case, marked with the black curve, reaches the total number of agents, N = 1000, with the transition at around ε c ≈ 0.6. The effect of the random vaccination scheme is analogous to that obtained by performing a reduction of the population density with a factor 1−f v , and thus the mobility parameter at the transition point will follow the relation ε c = ε (0) stands for the critical value of the mobility parameter when no agents are vaccinated. Interestingly, random vaccination and highest degree vaccination provide similar outcomes, with a slightly higher critical value of the mobility parameter for the random case. Indeed, this shows that highest degree vaccination is not effective at all. Highest LBC vaccination, on the other hand, results in a substantial reduction in the number of affected agents. Yet, the most effective vaccination schedule is, Figure 9 : Effect of different vaccination schedules on a fixed network: (a) random vaccination, (b) highest degree, (c) highest local betweenness-centrality (LBC), (d) highest betweenness-centrality (BC). Black dots denote susceptible agents, and they are joined by black links. Immunized agents are colored purple, same as the dead links. The parameters for all cases are N = 1000, ε = 0.89 and a 10% vaccination fraction. Videos of the numerical simulations for these vaccinations schemes can been seen in [49] . with a large difference, the one based on the highest BC (yellow line). This result is in agreement with previous studies [18] [19] [20] . Moreover, all schemes reach a similar value of the total number of affected agents for large mobility. Thus we conclude that an effective vaccination scheme will substantially increase the mobility threshold under which an epidemic outburst will die off. Yet, the effects of vaccination can be substantially reduced above that mobility threshold. How effective are the vaccination schedules shifting the percolation transition point? In order to answer that question we have traced the percolation threshold ε c as a function of the vaccination fraction f v in the bottom panel of Fig. 10 . The base value, for no vaccination, is ε c ≈ 0.6 (as shown in Fig. 5 ). Highest BC vaccination results in a sharp increase of the percolation transition, reaching ε c ≈ 1.1 for ≈ 10% vaccination fraction. Highest LBC vaccination performs substantially worse, but still provides a significant improvement of the epidemiological situation for scarce vaccines. Random and high degree vaccinations perform similarly, and none of them are quite effective. In this work we have presented a very simple agent model for epidemic expansion, the confined-SIR model, in which the effects of partial confinements and vaccinations can be easily tested. We have developed a theoretical framework to characterize the confined-SIR model, which combines results from percolation theory in random geometric graphs (RGG) with those of mean-field theory for epidemic expansion. It is relevant to notice that our model is far too simple to lead to policy recommendations without further insight from experts from different fields, ranging from virology to sociology. The main shortcoming of our model when applied to human epidemics is that human mobility is not geographically restricted in the same way as in our model. Indeed, even under lock-down essential workers must attend their workplaces, which might be far away from their homes. Yet, we expect that some of our conclusions can be interesting for researchers in epidemic expansion and lead, in combination with other insights, to sensible policy recommendations which will help alleviate the effects of present and future epidemic outbursts. Our first conclusion is that the effects of the confinement measures vary strongly when they cross the percolation threshold. Thus, it is of paramount importance to design lock-down measures so that mobility is restricted sufficiently below the percolation threshold in order to limit the expansion of the epidemic. The determination of the percolation threshold is not easy in practice, but a good hint is provided by the fluctuations in the outburst sizes. Above the transition point, most outbursts reach a huge fraction of the population, and below it, they will only affect a few individuals. Yet, near the transition, the number of affected agents may vary enormously, depending on the location of the initial infected agents. In addition, the effect of increasing the recovery probability in our model causes a decrease of the number of final infected agents, since it provides an efficient cutoff for the infection probability between pairs of agents which seldom meet. Our second conclusion is that in order to determine an efficient vaccination schedule, disregarding other healthcare considerations, bridge individuals should be especially targeted for vaccination, i.e. individuals which move between different clusters. Their immunization will lead to an effective confinement of an epidemic outburst to its initial cluster, thus creating effective firewalls between them. Bridge individuals can be detected via betweenness-centrality (BC), which is a global and computationally expensive measure, or through easier proxies, such as the individual local betweenness-centrality (LBC), which addresses the question: Are your friends friends among them? Individuals whose friends form a clique are not good candidates for vaccination, but individuals whose friends do not know each other are. Our results in this respect are mixed: vaccinating the highest BC agents is extremely efficient, but finding those agents is very hard both in the simulations and in real life. Vaccinating the highest LBC agents makes more sense: choose the individuals whose contacts are not in contact among themselves, such as retail salespeople. Yet, more research is required in order to find optimal vaccination schedules. Confined SIR models seem to be an appropriate tool in order to improve our intuition regarding the effectiveness of different strategies to stifle an epidemic outburst and to find relevant observables in order to characterize the current situation, allowing us to make meaningful predictions. For example, in this work we have emphasized the very interesting role provided by the fluctuations in the maximal number of infected agents. Indeed, a very promising line of research is the statistical analysis of these fluctuations during real epidemic outbursts, such as COVID-19. These analysis present a very interesting challenge: fluctuations should be compared ceteris paribus, i.e. removing major differences between the different geographical areas and times. Characteristics of SARS-CoV-2 and COVID-19 Debray et al, Prediction models for diagnosis and prognosis of covid-19: systematic review and critical appraisal Spatial heterogeneity can lead to substantial local variations in COVID-19 timing and severity A hybrid strategy for network immunization Global efficiency of local immunization on complex networks Five challenges for spatial epidemic models Modelling the impact of testing, contact tracing and household quarantine on second waves of COVID-19 Steinegger, Modeling the Spatiotemporal Epidemic Spreading of COVID-19 and the Impact of Mobility and Social Distancing Interventions Outbreak diversity in epidemic waves propagating through distinct geographical scales Seoane-Sepúlveda, A SIR-type model describing the successive waves of COVID-19 Nonmassive immunization to contain spreading on complex networks The turning point and end of an expanding epidemic cannot be precisely forecast Numerical identification of epidemic thresholds for susceptible-infectedrecovered model on finite-size networks Confined SIR model for epidemic expansion, software repository A contribution to the mathematical theory of epidemics An introduction of mathematical modeling of infectious diseases An Introduction to Percolation Theory Herrmann, Suppressing epidemics with a limited amount of immunization units Effective approach to epidemic containment using link equations in complex networks Efficient vaccination strategies for epidemic control using network information Random geometric graphs The longest edge of the random minimal spanning tree Random geometric graphs Large Connectivity for Dynamic Random Geometric Graphs Epidemic spreading in random rectangular networks Discrete-time Markov chain approach to contact-based disease in complex networks First-passage percolation, subadditive processes, stochastic networks and generalized renewal theory, in Bernoulli, Bayes, Laplace anniversary volume Models of first passage percolation 50 years of FPP Kardar-Parisi-Zhang universality in first passage percolation: the role of geodesic degeneracy Ad Hoc Netw Firstpassage percolation under extreme disorder: from bondpercolation to KPZ universality Epidemics, disorder, and percolation Percolation and epidemics in random clustered networks Early epidemic spread, percolation and Covid-19 Spreading of infections on random graphs: A percolation-type model for COVID-19 Applications of the variance of final outbreak size for spreading in networks Clustering and superspreading potential of SARS-CoV-2 infections in Hong Kong Superspreading drives the COVID pandemicand could help to tame it Unification of theoretical approaches for epidemic spreading on complex networks Spectral analysis of virus spreading in random geometric networks Spectral Analysis of Virus Spreading in Random Geometric Networks High prevalence regimes in the pair-quenched mean-field theory for the susceptible-infected-susceptible model on networks troduction to algorithms Identifying central nodes for information flow in social networks using compressive sensing Epidemic Modelling UNED website We would like to acknowledge J.E. Alvarellos, P. Córdoba-Torres, R. Cuerno, E. Korutcheva and S. Ferreira for very useful discussions. This work was funded by Instituto de Salud Carlos III (Spain) through Grant No. COV20/01081, and the Spanish government through Grants No. PGC2018-094763-B-I00, No. PID2019-105182GB-I00, and No. PID2019-107514GB-100.