key: cord-0886554-yxk9su29
authors: Ferreyra, Emanuel Javier; Jonckheere, Matthieu; Pinasco, Juan Pablo
title: SIR Dynamics with Vaccination in a Large Configuration Model
date: 2021-07-24
journal: Appl Math Optim
DOI: 10.1007/s00245-021-09810-7
sha: 97beac425d5f6efaf676d1ad952acc3e372585ad
doc_id: 886554
cord_uid: yxk9su29

We consider an SIR model with vaccination strategy on a sparse configuration model random graph. We show the convergence of the system when the number of nodes grows and characterize the scaling limits. Then, we prove the existence of optimal controls for the limiting equations formulated in the framework of game theory, both in the centralized and decentralized setting. We show how the characteristics of the graph (degree distribution) influence the vaccination efficiency for optimal strategies, and we compute the limiting final size of the epidemic depending on the degree distribution of the graph and the parameters of infection, recovery and vaccination. We also present several simulations for two types of vaccination, showing how the optimal controls allow to decrease the number of infections and underlining the crucial role of the network characteristics in the propagation of the disease and the vaccination program.

While epidemic dynamics have been studied extensively in the context of mean-field models where each individual potentially interacts with every other individual, there has been more recently a research effort to include the effect of local interactions using sparse random graphs, see [30, 35, 40] . The non-homogeneity in the inter-individual interactions can be reflected in the model by defining a degree-distribution describing the statistics of the interactions, see also the complete and pedagogical review [44] which contains both historical and modern references on epidemic processes on complex networks.

A rigorous mathematical description of epidemics on configuration models is a challenging problem, since a mean field approach implies to consider a system of infinitely many differential equations of SIR (susceptible-infectious-removed/recovered) or SIS (susceptible-infectious-susceptible) type. These systems are coupled by the distribution of links between nodes with k and j neighbors, for any pair of values k, j ∈ N. A seminal work in that direction was [14] , where a set of finitely many differential equations describing the asymptotic dynamics of a SIR on a sparse configuration model has been rigorously derived as projections of the infinite dimensional system. This is a remarkable result as it allows the model to grasp both the specifics of the interaction graph and the epidemic dynamics in a simple finite dimensional deterministic dynamical system. Their work agrees with other approximations in the literature, like the work of Volz which describe a Poissonian SIR epidemics using coupled non-linear ordinary differential equations [53] , and shows that these equations are indeed verified in the thermodynamic limit (i.e., when the number of nodes tends to infinity in a sparse configuration model). See also [29] for more on the SIR dynamics on the configuration model and [38] for equivalent formulation of the ODE dynamics.

On the other hand, vaccination processes were extensively studied in last years, since the anti-vaccine movements threaten social health programs by playing a common goods dilemma: they try to avoid the individual costs of vaccination and simultaneously pretend to enjoy the advantages of herd immunity. By modeling vaccination as a game, we are faced with the classical difference between individual and social optima, and worse equilibria are reached due to individual actions, than the ones obtained by a centralized planner.

So, the optimal vaccination problem has been studied using control and game theory tools, see [2, 22, 34, 49, 55] . Two main points of view are considered: a rational individual immersed in the population who maximize its own benefit, or a centralized agent who takes decisions for the overall population, for example a government.

The optimization problem for a centralized agent can be thought of as to minimize the costs of a vaccination program while preventing the epidemic spread of a disease, for example, achieving herd immunity. In [26] , the authors consider a deterministic epidemic compartmental model and use a time-dependent vaccination rate and linear costs to get the optimal vaccination strategy.

Let us observe that the authors in [21] propose an evolutionary game-theoretic problem, where individuals use evidence to estimate costs of vaccination, the model being based on the agent point of view. Vaccination strategies can also be influenced by a neighborhood behavior and may depend on the individual's beliefs about their neighborhood vaccination strategies [36, 45] . Another approach can be found in [22] , where they study how the psychology of individuals intervenes in their perception of their risk, susceptibility or mortality rates.

Regarding a network background, the authors in [54] show that the vaccination is most effective when the full network structure is known by the individual agents although in the real world the decisions are based on partial information of the contact underlying graph. This could justify that individuals decide to get vaccinated with rates that depend on their degree in the network. Highly-connected individuals (hubs) have a high incentive to vaccinate, whereas individuals with few contacts have less incentive to vaccinate as exposed for example in [13, 36, 49] . We refer the interested reader to [9] for an extensive review of compartmental models for epidemic modeling, both in mean field setting and in networks mode, the discussion on the trade-off between simple vaccination models which miss a lot of details but are very useful to reach a general qualitative analysis, and more detailed models usually designed for quite specific diseases and populations.

In [3, 11] , the authors considered different local vaccination strategies, where randomly chosen neighbors of the selected individuals are vaccinated, prior to the introduction of the disease. They obtained formulas for the final size of the epidemic as function of the vaccination intensity using branching processes techniques and generating functions.

In this work, we study degree related vaccination strategy including the case of the "proportional to degree" strategy, which is very similar to the acquaintance vaccination model, originally proposed in [12] . Let us remark however that a major difference in our model is that the vaccination process and the epidemic occur simultaneously. We consider a Markov-variant of the SIR model, where infectious periods are exponentially distributed, see for details and more general models [4] . We adapt the techniques developed in [14] and [29] to show the convergence of the degree measures that describe the Markovian dynamics of propagation on the random graph in this context and obtain a generic and deterministic description of the epidemic evolution. As in the case without vaccination, we are able to derive a finite-dimensional differential system describing the evolution of the quantities that usually describe the epidemics, namely the number of individuals in each compartment (Susceptible, Infected, Recovered and Vaccinated) and the basic and effective reproduction number, whose definitions and description are in Sect. 5.

The vaccination model deals with a context where individuals do not know if their contacts are infected or not (otherwise, they might actually suppress them), and could have an incentive of vaccination based on their social interactions. Usually, macroscopic indicators at the level of city or state give a clue of the severity of the epidemic, although as in the case of the actual Covid-19 pandemic, an exact knowledge of the state of your contacts is challenging due to the presence of asymptomatic individuals.

Then, we study the optimal controls for the vaccination formulated as a game both in a decentralized and centralized setting following ideas developed in [16] (in a purely mean-field setting). We show in particular that the optimal vaccination strategy boils down to a bang-bang control, i.e., the optimal solution consists in vaccinating with the highest possible rate until some fixed time-threshold depending on the connectivity of the network and the costs, and then not to vaccinate at all anymore. On the other hand, using techniques from continuous optimization for systems with restricted controls, we first define a general assumption that allows proving uniqueness and existence of optimal control in the sense of viscosity solutions following [10] and [51] . We also show that the optimal centralized strategy must be developed at the highest rate possible when deployed. We could not however prove that there is only one phase of vaccination.

Finally, we consider four network examples, compute the optimal vaccination strategy and simulate the propagation of the disease and the vaccination program. We observe that the optimal strategy has the desirable feature to reduce considerably the total number of infected individuals, while a more conservative vaccination program would increase costs without having significant effects on the epidemic. We also relate our results to graph measures as centrality coefficients and to the classical indicators associated to the epidemic.

In Sect. 2, we introduce the necessary notation and the epidemic model. For the sake of completeness, we add a short description of the configuration model random graph, together with the relevant measure spaces considered.

In Sect. 3 we generalize the results of [14] describing the propagation of an epidemic on a configuration model random graph by incorporating a strategy of vaccination for the susceptible population depending arbitrarily on time and out-degrees.

We first present our main scaling limit result: we obtain an infinite system of measure valued differential equations that describe the Markovian dynamics of propagation of a disease on the random graph, subject to a vaccination strategy, see Theorem 1. Then, we derive a finite-dimensional differential system describing the evolution of the main variables that describe the epidemic, namely the number of individuals in each compartment (S, I , R, V ), and the probability of interaction between the susceptible population with agents in different compartments.

In this edge-based system, p X will denote the probability that an edge connects a susceptible with a node in state X , for X = S, I , R, or V . We use a generalization of the probability generating function g, associated to the initial degree distribution and depending on two variables that describe the probability that an individual remains susceptible: α, the probability that an edge has not transmitted the disease to a given node (which help us to describe the probability that the node is not infected), and θ to describe the probability of no vaccination. The closed system we get is the following (see the next section for the precise definitions and notation):

where r , γ and π t are the contagion, recovery and vaccination rates. The heterogeneity of the population connections, which is an important feature in many propagation models [2, 40, 44, 52] , can be grasped here through the function g of the degree distribution; in our case, through the expression α t ∂ αα g(α t ,θ t ) ∂ α g(α t ,θ t ) . We also write the system in the cases where the vaccination rate depends linearly on the degree of a node or is constant. Here the system relies on the probability generating function of the initial degree distribution ψ and the mean excess degree of the nodes generalizing known results [42, 53] .

Let us observe that this heterogeneity is not present in classical mean field models which rely on the assumption of a completely homogeneous mixed population. However, we will show that the system (1) in a configuration model with a Poisson degree distribution converges to the classical SIR model when the mean degree goes to infinity. The system of equations (1) allows us to compare agent based simulations or real data with the curves obtained by numerical integration [19] .

In Sect. 4, we study the optimal control problems associated with vaccination. In order to consider the effect of the graph structure, we define a vaccination strategy as a bounded and measurable time dependent function, followed uniformly by all the individuals in the population but depending on the connectivity of an individual. We assume that the rate of vaccination is not decreasing in the degree, in agreement with the literature where highly connected individuals have more incentive to get vaccinated whereas individuals with few contacts have less incentive. On the other hand, this is already a large family of controls, needing a quite general theoretical treatment, involving in particular weak viscosity solutions.

The optimal vaccination policy will vary if we consider the individual point of view. This decentralized case is approached using the theory of mean field games, by considering the perspective of a single rational individual added to an infinite population with a given arbitrary vaccination policy. We show in particular that the optimal vaccination strategy boils down in that case to a bang-bang control, i.e., the optimal solution consists in vaccinating with the highest possible rate until some fixed stopping time depending on the connectivity of the network and the costs, and then not to vaccinate at all anymore,

where ν is the rate of vaccination and 1 [0,τ ] is the characteristic function of the set [0, τ ]. Then, we consider the social optimum, and we introduce a particular cost functional in order to find the optimal centralized strategy. As before, we show that the optimal strategy π t is of threshold type, taking either the maximum value or zero.

In Sect. 5 we analyze the theoretical results and we present numerical computations. We obtain a modified Basic Reproduction Number for the epidemic on the graph before the vaccination process started, namely

Therefore, an epidemic outbreak will occur with strictly positive probability if R 0 > 1, and corresponds to the already known critical threshold stated in [41, 44] . We compute the optimal vaccination strategy and we simulate both the propagation of the disease and the vaccination program. We observe that the optimal strategy has the desirable feature to reduce considerably the final size of the epidemic, depending on the contagion and recovery rates r and γ , the function g, the threshold τ and the rate of vaccination ν, while a more conservative vaccination program would increase costs without having significant effects on the epidemic. Finally, we consider different networks generated by four degree distributions, all with the same mean degree first, and with the same R 0 later: (a) Poisson, which corresponds with mean field models with homogeneous mixing, as we show in Sect. 3.2.1; (b) Bimodal, where a fraction of the population has few links distributed Poisson, while the rest of the nodes has many links, as proposed in [9] . (c) a Regular graph, as an example of networks where all nodes are equals on its rates; and (d) Power Law, a classical model of social networks [42] .

We conclude in Sect. 6, and the full proofs of the different theorems can be found in Sect. 7.

Let us introduce the stochastic environment of the epidemic process. We use the configuration model random graph introduced by Bollobás [7] which can be constructed as follows. First denote N n 0 = {0, . . . , n} and suppose we have n nodes, and a sequence of degrees k 1 , . . . , k n independent and identically distributed according to p (n) = ( p (n) k ) k=1,...,n−1 1 such that the sum of the degrees is even. One initially assigns a quantity k i of half-edges to the ith-node and then choose two of them uniformly from the unmatched ones, establishing the connection between the nodes, until all the half-edges are matched. Nodes will represent individuals, thus we will use both terminology throughout the paper.

Under the assumption that p (n) converges in probability to p = ( p k ) k∈N 0 and E[ p (n) 2 ] converges to E[ p 2 ] < ∞ when n goes to infinity, there is asymptotically a probability bounded away from 0 of obtaining a simple graph as showed in [28] . Thus, we may repeat the matching procedure until we obtain a simple graph, i.e., until the resultant graph does not contain self-loops nor multiedges [18] .

As a consequence, the degree of a randomly chosen node is distributed according to p (n) . Now, the probability that a neighbor has degree k is kp k j=1,...,n j p j , which is the so-called size-biased degree distribution. This distribution will greatly influence the dynamics on configuration models as opposed to the mean field point of view, where the underlying graph of potential connections is complete.

Given a degree distribution p = ( p k ) k∈N , the associated probability generating function is defined by ψ(z) = k∈N p k z k .

We now describe the dynamics of the epidemic with vaccination. For a given Susceptible node (i.e., not having contracted the illness nor vaccinated) we consider several independent exponential clocks with parameter r , one for each edge connecting an Infected node (i.e., potential encounters). It will describe the contact process: if this clock rings, the Susceptible makes a transition to state Infected and remain infectious during an exponential time with mean 1/γ , whereupon it will not longer infect any node, going to Recovered state. Also for the Susceptible, we consider a nonhomogeneous Poisson Process with rate π t (k) depending on the degree k of the node, which is the time dependent control variable and represents the rate at which the individual becomes Vaccinated, indicated by the jumps of this process.

Both Recovered and Vaccinated states are absorbing states of the resulting continuous time Markov chain, and we differentiate between both to keep track of the epidemic characteristics. Also, they have a different impact on the total costs of an epidemic.

In the sequel, we suppose that π t (k) = ξ(k)π t , where π t is a measurable bounded function on [0, T ] and ξ : N → R. For the existence of the fluid limit, some technical assumption on ξ is needed, and we suppose that

where μ S 0 is the initial degree distribution of the susceptible population. The cases ξ(k) = 1 and ξ(k) = k appear more frequently in the literature of local vaccination strategies. Constant rate indicate homogeneity of the agents on the immunization effort; and a rate proportional to the degree, as in the case of acquaintance vaccination [3, 11, 12] , takes into account that vaccination of hubs and higher degree nodes is convenient to stop the transmission of a disease. This can also be achieved using a targeted vaccination where higher degree nodes are vaccinated, which corresponds to ξ(k) = 1 κ≤k (k) for some threshold κ. In the acquaintance vaccination, one samples a fraction of the nodes, and for each sampled node, one of their neighbors uniformly chosen is vaccinated. Hence, the probability of a node to be selected by any neighbor is proportional to their degree since it has the size biased distribution.

We denote S n t , I n t , R n t and V n t the total number of Susceptible, Infected, Recovered and Vaccinated nodes respectively, and S n t , I n t , R n t and V n t the sets containing the nodes in each state. These quantities are of central interest in the literature and they are the main variables describing the dynamics (we refer the reader to [30] for an informative review on epidemics dynamics). Nevertheless, let us remark that our model being over a Configuration Model random graph, computing its dynamics is in principle very demanding, as we should study a stochastic process in (growing) dimension n, the number of nodes.

The principal signification of part (ii) in Theorem 1 below is that the limit behavior is a good approximation of the case n large.

Instead, we study the dynamics of four measures in M F (N 0 ), the set of finite measures on N 0 embedded with the topology of weak convergence (we refer the reader to [6] for a complete description of topological properties and results), describing the connection between the susceptible population and the rest. For this purpose, we resort to the so-called principle of deferred decisions, revealing the graph simultaneously with the propagation of the disease, regarding the types of the edges connecting the different states of the individuals. This trick is possible since the random environment for the epidemic dynamics acts over a configuration model which is constructed using a uniform matching. We follow here the ideas developed in [14] and [5] .

We now describe the quantities involved in our formulation of the dynamics precisely. For the ith-node, we denote k S i the random number of edges connecting the ith-node to a susceptible individual, and δ k the Dirac measure on N 0 . The empirical measure μ S,n t ∈ M F (N 0 ) describes the degree of the susceptible individuals: for each k ∈ N 0 , μ S,n t (k) denotes the number of susceptible nodes with degree k at time t,

Similarly,

. Instead of considering the proportion of the population in each state, our work describes the dynamics of the following edge-based quantities:

the number of semi-edges connecting a susceptible node; and, analogously, N I S,n t , N RS,n t , N V S,n t are the number of edges linking a susceptible node with an infected, recovered or vaccinated one, when the size of the population is n.

We also consider the proportions of edges associated to those quantities:

Our main result consists in the convergence of the normalized empirical measures in the Skorokhod space embedded with the weak topology. We write the semi-martingale decomposition stated in Proposition 5 and a description of the Markovian process through a system of stochastic differential equations derived form Poisson point measures, see Sect. 7.1.

For each μ ∈ M F (N 0 ), the set of finite measures on the natural numbers including zero, and f ∈ B b (N 0 ), the set of bounded real functions on N 0 , we write

In order to show that the normalized degree empirical measures of each type converge to the solutions of an infinite system of differential equations as the population size tends to infinity, we scale the measures in the following way: for n ∈ N, we set

The suppression of n will denote the limit measures associated to the fluid limit

In the limit we also define:

that allow us to describe the probability that a node with degree k has not contracted the disease, α k t , and the probability of no vaccination at time t, θ ξ(k) t .

We now state the convergence of the degree empirical measures when the size of the population tends to infinity. The corresponding deterministic solution of the measure-valued system will in turn give interesting insights on the effect of the vaccination in the propagation of the epidemic for large populations. The proof is similar to the proof of the main theorem in [14] , for a SIR without vaccination. We use strongly the fact that π t is uniformly bounded on time, and this introduces several differences in steps 2 and 5 (compared to their original proof), where we need to invoke results from [15] in order to ensure uniqueness of weak solutions to the transport equation involved. We add it in the last section for completeness. (ii) the sequence (μ (n) ) n∈N converges in distribution to μ in the Skorokhod space when n goes to infinity.

The derivation of the system of equations (3) can be explained by the underlying Markovian process describing an epidemic on a configuration model graph, which is also useful for a possible simulation. The edges based measures must be updated when any of the three possible events occur: the infection or vaccination of a susceptible individual (with a rate that depends linearly on their degree k), and the removal of an infected one (according to exponential clocks of parameter γ ).

Let us recall that probability that an individual of degree k that started susceptible remains in that state at time t is θ ξ(k) t α k t , which gives the first equation. In the second equation, the first integral accounts for the removing of infectious individuals; the second one corresponds to the addition of the newly infected, rkp I is the rate of infection for a susceptible of degree k and the multinomial part is the random draw of the neighbors according to the probabilities of the edges. The third and fourth integrals correspond to the infection or vaccination of the neighbors of each infected individual, whose degree is chosen according

, the size biased distribution. The number of edges to the susceptible population decreases by one in both events, infection and vaccination.

Now (3) we obtain a countable system of ordinary differential equations that allows us to describe the infection propagation in terms of the measures:

In Proposition 1, we now derive a generalization of the equations proposed by Volz [53] including a generic vaccination strategy. If the vaccination function was constant or bounded and continuous, then existence and uniqueness would be trivial. However, we are allowing the control to be bounded and measurable implying the necessity of a general treatment as the one studied in [10] or [51] which states that the functional describing the dynamics needs to be Lipschitz. This hypothesis is straightforward using classical computations and bounds together with the properties of the involved functions.

We define the function

and denote with a subscript the partial derivation of g, namely:

and analogously for the rest. This function will be key to reduce the number of equations to seven, obtaining in this way a tractable description of both the epidemic dynamics and the optimal vaccination strategy.

In the limit of infinitely many nodes, the system (3) can be reduced to the following system of differential equations,

Moreover, the right-hand side of the system is Lipschitz and uniformly bounded, and the problem (5) hence admits a unique solution for any initial datum and any measurable π

Proof We use the function g to compute a closed expression for N S t , N I S t , N RS t and N V S t and its derivatives. We denote 1(k) := 1 and χ(k) = k.

is the proportion of susceptible individuals at time t. Similarly,

The next step is to find the dynamics for p S t and the other edges probabilities. Before computing it, let us note that:

Using the definition of p I t andα t = −r p I t α t , we obtain:

.

We replace f by χ in (3), and after some basic computations, by rearranging terms using the multinomial theorem, we havė

Reasoning similarly with the other probabilities, and putting all the equations together, we have, for each control π : [0, T ] → [0, ν], the closed system of equations (5) .

This finishes the proof.

The particular case ξ(k) = ak+b which includes constant and acquaintance vaccination can be studied using only the generating function of the initial degree ψ(z) = k∈N μ S 0 (k)z k of the network. To this end, we define the quantities β t = e − t 0 r p I s +aπ s ds and φ t = e − t 0 bπ s ds , and rewrite the system in terms of these new variables. Let us note that the derivation follows from system (3) . Noticing that S t = φ t ψ(β t ), the dynamics are described by:

Observe first that the dynamics of the epidemic depend strongly on the degree distribution, since the expression β t ψ (β t ) ψ (β t ) governs the expected number of susceptible individuals connected with a neighbor of a given degree at time t.

Let us quickly clarify the differences between this model with a linear vaccination rate and the usual mean field model. In the sparse Erdös-Renyi model, the number of neighbors in the graph follows a binomial distribution, which can be approximated in a large population by a Poisson distribution. On the other hand, when the graph is fully connected and the contact process is determined by a Poisson process, the number of neighbors with whom each node effectively connects is also Poisson distributed.

In the Mean Field model with vaccination [16] , an individual of an homogeneous population encounters others following a Markov process in continuous time with rate r . So, individuals can be in three states: susceptible, infected and recovered or vaccinated; and we denote S t , I t and R t , its respective proportions of the total population. Here, vaccinated individuals are treated as recovered, since the influence in the propagation of the disease is the same. If the initial individual of the contact process is susceptible and the encountered one is infected, the first one becomes infected. An infected individual recovers at rate γ , and a susceptible can choose its own vaccination rate π , going to recovery state. The optimal strategy π played for all the players is called a mean-field equilibrium, defined as a fixed point of the best response functional, this is, π ∈ B R(π ), which minimizes a properly defined cost functional [17] . A mean eld equilibrium consists in a strategy where no player has an incentive to deviate from the common strategy for their own benefit. When the size of the population goes to infinity, the dynamics of the population where all players use the vaccination strategy π is described by the following system of equations:

The main difference between the two systems lies on the term associated to the infection process in the mean field case, r S t I t , which is now r p I t φ t ψ (β t ) = r N I S t and the neighbors are chosen according to the size biased distribution. In the particular case of Poisson distribution, the size biased distribution is also Poisson, and, as we show in this section, both models are asymptotically similar.

Here, we are considering an edge-based dynamics instead of an individual-based one. In the propagation of the disease, exponential clocks of the contact process are assigned to edges connecting susceptible with infected nodes, the mechanism is not to choose uniformly between all the individuals in the population but in a neighborhood, which modifies quantitative and qualitatively the generator of the Markov process and therefore the limit equation.

Though the dynamics are not the same, (even in the case of a Poisson distribution for the configuration model), one can actually show that they are asymptotically equivalent, when the mean number of connections between individuals grows large.

If we consider a population of size N in which every individual has C possible contacts and scale it such thatr = rC remains constant, the mean field equation of these dynamics is described by (9) replacing r and π byr andπ respectively. Then our limiting dynamics on a configuration model with Poisson degree distribution, not fully connected but uniformly linked, solve the MF equations when C goes to infinity taking ψ(z) = Ce C(z−1) . Being S = φψ(β), we have: which is asymptotically the first equation in (9) . The third equation follows from a similar reasoning, and the second one from the fact that S = 1 − I − R.

As underlined before, we aim in this section at describing the optimal strategy which will consist in vaccinating at the maximum rate as soon as possible up to some critical time (not necessarily deterministic) depending on the parameters of the model.

We define a generic vaccination strategy, which is time dependent and also depends on the individual connection in the network environment. Using the theory developed in [10, 51] , we characterize under mild assumptions the existence and uniqueness of a viscosity solution for our optimization problem. Of course, the model is a crude simplification of a very complex reality, but we believe it still allows to grasp an interesting phenomenology (see Sect. 5.1). We do not involve explicitly in our setting more individual features like age, risk, or beliefs, nor susceptibility to the disease. However, many of these characteristics may be modeled through adequate parameters and cost values.

In this section, we analyze the optimal control problem from an individual point of view. We focus on the perspective of a particular individual immersed in the population, who will take decisions in order to minimize their cost in a game against the whole population. This rational individual can be considered as a player seeking for the best response to a fixed strategy followed by the population, and we are therefore in the context of the theory of mean field games, which shall provide our context, definitions of equilibrium and existence results.

Suppose we add an individual to a population that evolves according to a vaccination strategy π . Since the population is infinite, the behavior of this new individual will not affect the evolution of the whole population, hence its dynamics will be described by the already stated equations (5), which we represent in the formẋ = ϕ(x, π). We denote byπ the vaccination strategy for the new individual andx = (S t ,Ĩ t ,R t ,Ṽ t ) their probability distribution over the four possible states.

Finally, let us remark that despite the mean field game denomination, the dynamics of the population follows the network based system that we described in Sect. 3.2.

We suppose that the new individual has degree k, therefore, we have that

SinceS t = α k tθ ξ(k) as before, but depending on the dynamics of the population infected edges and on their own vaccination rate, their state is determined by the following system of the formẋ = f 0 (x,x, π,π):

We do not write the dependence of the variables on time nor on the degree of the node in order to avoid heavy notations. We consider a cost function similar to the one proposed in [16, 26, 31] , which contains a linear term in the vaccination rate c V π t where c V may depend on the cost of the vaccine and its possible side effects; and a term modeling the cost incurred by an infected patient per time unit, which could include the possible loss generated by being unable to attend work, the costs of treatment and medical consultation, and could consider the severity of an illness such as its sequels or even death. Hence, the new individual wants to minimize their cost defined by:

So, the new individual looks at the best response to strategy π , i.e., they want to play B R(π ) ∈ arg minπC t (π,π). The minimum is taken over

which is a compact set for the weak topology. This implies that B R(π ) is not empty, since any minimizing sequence has a limit.

In Theorem 2 of [17] , the authors show that if the dynamics and the costs are continuous functions of the involved variables there always exists a mean-field equilibrium in such games; thus the existence of a solution for our problem follows from their result, since the cost is linear in I , the function g is analytic, and the rates of transition for the new individual depend linearly on p I or are constant for a fixed π .

Although this result guarantees the existence of a mean field equilibrium, we will compute the best response strategy for any π played by the population analyzing our problem as a continuous time Markov Decision Problem with finite horizon.

Denote J S (t), J I (t) the optimal cost starting at time t in states susceptible and infected, respectively. The optimal cost J and the strategyπ * that realize it, satisfy the following Hamilton-Jacobi-Bellman optimality equation [46] :

The equation explains how state transitions and their rates impact on the individual expected cost.

Proposition 2 Letπ * be the strategy of a node with degree k that realizes the optimal cost J . Then,π * is threshold, that is,π * = ν1 [0,τ ] (t) for some τ ∈ [0, T ].

Proof We will prove that the optimal strategy is constantly the maximum rate of vaccination until some time τ , and after that instant, the optimal strategy boils down to zero. Let us remark that the costs associated to two different strategies that differ in a null measure set stay the same. So, we have uniqueness up to a zero measure set. We can see from the third equation in (12) that

Hence J I decreases from J I (0) = c I γ (1 − e −γ T ) to J I (T ) = 0. Let us also observe that if J S (t) > c V thenπ(t) = 0. Since J S (T ) = 0 and the costs are continuous,J S (T ) = 0 therefore, if we call τ the first instant at which J S is below c V , we have J S (t) ≤ J I (t) for all τ ≤ t ≤ T , such that the second term in the second equation in (12) is non-negative. If J S does not cross c V , then we take τ = 0.

Moreover, if the cost of being susceptible is bigger than the vaccine cost, the derivative of J S will be even smaller. Hence, J S (t) ≤ J I (t) for all 0 ≤ t ≤ T and therefore J S is always decreasing before τ . This concludes the proof.

Let us remark that the optimal cost depends on the degree k of the individual, thus π * t =π * k t . Moreover, if we suppose that ξ(k) is increasing in k, meaning that the nodes with more contacts have more incentive to vaccinate, the derivative of J S (t, k) is always bigger than the derivative of J S (t, j) when k is bigger than j, and both J S and J I are equal at the beginning of the propagation of the disease if the initial proportion of infected is small enough. Thus, the threshold τ k will be bigger than τ j . We have:

be the optimal strategy for an individual of degree k. Then τ k > τ j if k > j.

Now we consider a centralized planner trying to optimize from the point of view of the total population. We first consider the control system of the formẋ = ϕ(x, π) like (1) where the set of admissible controls is compact, and the family of admissible control functions π is only restricted by its measurability. Given the initial data x(0) = x 0 the Cauchy problem has a unique solution, as we stated in Proposition 1.

Given an initial data (s, y) we consider the general optimization problem:

where L is the instant cost functional, is the final cost, and the state variable x depends not only on time but on the control and the initial data. The optimal π is taken over the set of measurable functions π : [0, T ] → [0, ν].

As stated by the dynamic programming method, the optimal control can be characterized by the value function V (s, y) := inf π ∈ J (s, y, π), but the classical point of view does not allow discontinuous control functions. Hence, we start by verifying the hypothesis that our setting must satisfy in order to ensure existence and uniqueness of the optimal control, based on more general results on viscosity solutions theory [10, 51] . In the case of measurable control, we can also apply the Pontryagin's Maximum Principle with less restrictive assumptions.

According to Lemma 9.2 in [10] , the functionals involved must satisfy

for all x 1 , x 2 ∈ R 7 , and π ∈ , for some constant C. Under these assumptions, the value V is a bounded, Lipschitz continuous function, and it can be characterized as the unique viscosity solution to a Hamilton-Jacobi equation.

Since the epidemic dynamics satisfy these assumptions on ϕ, to the best of our knowledge, the most general condition on the cost functions to admit a solution is to be Lipschitz in the state variable and bounded, which allows us to model a wide range of real situations.

As a particular case inspired by the individual optimization problem exposed above and in agreement with [16, 22, 26, 31] , we define the cost:

We can easily check, by basic computations and bounding the second term using the regularity of g and the Mean Value Theorem, that our setting satisfies hypotheses (14) .

Further, given the data x(0) = x 0 , let t → x * (t) = x(t, π * ) be an optimal trajectory corresponding to the optimal control π * . Following Theorems 7.18 and 11.27 in [51] , there exists an absolutely continuous application t → p(t) ∈ R 7 called the adjoint vector, and a real number p 0 ≥ 0, such that ( p, p 0 ) is non trivial, and such that for almost every t ∈ [0, T ]

x * = f (x * , π * ),

where the Hamiltonian of the system is H = p 0 L 1 + p f . Summarizing, we have the following result.

Let π * the strategy that minimizes (13) for the cost functions defined above. Then π * is threshold.

Proof Writing the equation for π * we get

which is a linear function of π with principal coefficient ρ * . Since we are minimizing over π ∈ [0, ν] we can conclude that

and the proof is finished.

Since it is impossible to solve analytically the system (16), the method of Forward-Backward Sweep presented in [37] can be useful to understand the behavior of ρ * . From (14) we can deduce that ρ * (T ) > 0 indicating that the vaccination rate must be zero after some time. As in the preceding section, the optimal vaccination strategy is of threshold type, in the sense that vaccination must be intended with the maximum effort (maximum rate) and otherwise discouraged. This result may be explained by the intuition that the individuals would vaccinate if the vaccination cost is lower than the potential cost associated to the infection, and the probability of ever contracting the disease decreases from the start of the epidemic. However, remark that a quantitative characterization of the change of regime depends on the complete behavior of the epidemic. This is solved by calculating Bellman equations and using the theory of Dynamic Programming. lation and mild economical cost. Note that it is natural to suppose that there exists an upper bound on the vaccination rates which takes into account both economical and organizational considerations.

Knowing the effective contact rate r of the disease in consideration and its recovering rate γ , the vaccination budget represented in ν, and taking in account the connectivity of the population, we aim at finding this optimal vaccination effort.

In usual SIR mean field models, there exist two regimes: above a threshold in the ratio r γ the epidemic will propagate and the final number of infected reach a fraction of the population even when the number of initially infected ε is arbitrary small. Below this threshold the final size of the epidemic will be proportional to ε.

Here, this threshold can be described in terms of the connectivity of the graph which depends on the average number of contacts for the size biased distribution. In our model, the number of newly infected in a small time interval is proportional to the quantity p I . The threshold is then determined by the following equation:

If we take an ε-proportion of initially infected individuals, then we have the initial conditions:

Assuming ε 1, after a simple computation we get the inequality:

As we mentioned before, ∂ αα g(1,1) ∂ α g (1, 1) is the expectation of the size biased distribution, therefore it represents the mean degree of a random neighbor (a new infected agent). By subtracting one (the spread will not return to the infecting agent), we get the expected out degree, which represent the possible number of contacts that the newly infected individual has. Therefore, Eq. (20) states that the epidemic outbreak will occur if the total rate of infection is bigger than the rate of recovering. We can rewrite it and find back the already known critical threshold stated in Prop 6.1 of [44] :

Written in terms of the transmissibility T = r r +γ (i.e., probability of infection given contact between a susceptible and infected individual) we recover the formula stated in [41] for the epidemic threshold.

Finally, for the optimal strategies that we found in the previous section, we compute the impact of the vaccination process in terms of epidemic sizes, showing that the decrease of the susceptible population can be described through the generating function of the network, the proportion of IS edges and the vaccination rate, through a simple formula.

Following Miller [38] , we can compute the final epidemic size in our model which allows to measure the impact of the vaccination process. To that end, we compute α k ∞ the probability that a randomly chosen susceptible node u with degree k is never infected in terms of the transmissibility T = r r +γ , interpreted here as the probability that one of the edges of u with an infected neighbor transmit the disease to u. The probability that this neighbor is never infected, given that u does not transmitted the disease, is (1, 1) because the neighbor degrees follow the size biased distribution. Thus, the probability that the edge is not an infected contact for u solves the fixed point equation:

we have:

Proof After some simple computations we get

Now, we can also compute the probability that an initially susceptible node remains susceptible at the end of the epidemic, corresponding to never be infected nor vaccinated:

With this,

In order to compute the size of the epidemic (total number of recovered agents), first we compute the final number of vaccinated agents. This can be thought of as an exponential race with non-homogeneous rate. Thus, the probability that a node of degree k will be vaccinated is:

where λ π k (t) = t 0 ξ(k)π s ds and λ I k (t) = t 0 rkp I s ds. With this, V ∞ = k μ S 0 (k)V ∞ (k), and therefore the final epidemic size is R ∞ = 1 − S ∞ − V ∞ . The case of the degreeproportional vaccination, ξ(k) = k, can be described by the generating function of the initial degree of the susceptible. By taking

, one can write:

On the other hand, α ∞ is the fixed point of α ∞ = 1 − T + T ψ (α ∞ e −τ ν )/ψ (1) and S ∞ = ψ(α ∞ e −τ ν ). This finishes the proof.

Here, we can see the strong dependence on the connectivity model of the network through the generating function g (or ψ) and on the maximum rate of vaccination in order to reduce, (exponentially in τ ν and through the functions g and ψ respectively), the propagation of the epidemic. Additionally, the maximum rate of vaccination can be translated in the budget of the decision maker, because it may indicate how effective the decision to vaccinate may be.

Currently, the development of vaccines is usually subsequent to the appearance of a possible epidemic. In many of them, it is through the value R 0 that the epidemiological parameter r (or β for homogeneous compartmental models) can be estimated. It contains the probability that in an encounter between a susceptible and an infected person, the former becomes ill, and the number of encounters in a unit of time; and γ which is related to the average time that a person is infectious. Knowing these parameters, the quotient ∂ αα g(1,1) ∂ α g (1, 1) can be inferred from the formula (1.1). In this section we fix the parameters of the epidemic and simulate the propagation of the disease and the vaccination process in the cases ξ(k) = ak +b solving numerically the system of equations (8) .

We suppose a threshold vaccination strategy of the form π t = ν 1 [0,τ ] in the period [0, T ] thus the global optimization cost is

We do not write explicitly the dependence of the variables on τ to reduce notation. Let us remember that an agent of degree k will vaccinate according to the rate ak + bν during the vaccination period [0, τ ], which is hidden in the probability generating function g.

Hence, the optimization problem reduces to find τ * = arg min τ C(τ ). Thus, for each τ we run the numerical integration and select the optimal threshold that minimize this cost.

We do this for two vaccination policies, constant and degree-proportional. In order to compare both strategies, we set a = 1 and b = ψ (1), the mean number of neighbors, obtaining the same vaccination rate.

As we stated in Sect. 3.2.1, the mean field model corresponds to taking a Poisson degree distribution. In [9] , the authors establish that the Poisson distribution is not realistic for the modeling of contacts, and propose to use a Bimodal distribution in which a proportion p of the population follows a Poisson law of mean C and the rest a delta-L distribution. If the whole population follows a delta distribution, the associated graph is called regular. Other distribution were also proposed in [25, 39, 41] , in particular the power law, found in several social networks [1, 32, 33, 43] .

Thus, we consider four different networks of size N = 10000, associated to a correspondent generating function of the degree distributions, choosing the parameters in order to get the four with mean degree ψ (1) = 5: Here we can see that the graph structure makes a difference for the propagation or eradication of the disease through the existence of hubs that could play the role of super-spreaders, or almost disconnect the networks due to vaccination at higher rate. We summarize some network measures of the generated graphs, and show the numerical solutions of system (5) in each of them in order to compare the spread of the epidemic. We use the python package networkx for the graph creation and the calculus of the coefficients [24] , and we solve the system (5) using a Runge-Kutta 4.

We can see in the plot of Fig. 1 , that the Bimodal network has typically more nodes of high degree than a Poisson (with same mean), and less than a Power Law network, according to the results in [9] .

For the epidemic evolution we take r = 3 and γ = 1, and in Table 1 we can see that networks (b) and (d) have a bigger basic reproduction number due to the size biased distribution and the expected excess degree.

Typical parameters of the nodes of social networks are the clustering coefficient, which indicate the tendency of a node to form a cluster, and the betweenness centrality, measuring the presence of the node in the shortest paths between different nodes, see [27] . Both play an important role in the propagation of a disease, the former due to potential contagion to larger groups of nodes, and the later because imply a quick connection between different parts of the network. Also, the closeness coefficient measures how far apart are the nodes, see [8] . Proper definitions of these coefficients could be:

• Betweenness centrality of a node v is the sum of the fraction of all-pairs shortest paths that pass through v. • The density for undirected graphs is 2m n(n−1) , where n is the number of nodes and m is the number of edges.

• The clustering of a node v is the fraction of triangles passing through that node that effectively exist. • Closeness centrality of a node v is the reciprocal of the average shortest path distance to v over all n − 1 reachable nodes.

See [42] for more details. In Table 2 we show the average number for these measures after several experiments with networks of 10 4 nodes. We observe that density, related to the mean degree of the four networks, is similar. Nevertheless, there are important differences between the averages of the clustering coefficients.

In Fig. 2 , we present three simulations for a Bimodal graph for different optimization costs and vaccination rates, for the constant (abbreviated Const) and degree proportional (abbreviated Deg) vaccination strategies. We plot the total cost, the final number of vaccinated and the epidemic size as a function of the vaccination threshold τ , obtaining the same behavior and shape in all the studied cases, namely, the final vaccinated population is increasing in τ , and the costs decrease up to the optimum.

We can see from the simulation of the optimization in Fig. 2 that the epidemic size decreases notably with respect to the case without vaccination which corresponds to the vaccination threshold τ = 0. On the other hand, we observe that the epidemic size curve flattens around the value τ * that minimizes the cost. Therefore, it makes no sense (neither financially nor epidemiologically) to continue a vaccination program after that time, because the final cost will be bigger and the population infected will not decrease further with more vaccinated individuals.

In Fig. 3 , we present the result of the simulation for the variables S, I , R and V when ν = 0.2, c V = 10 and c I = 50, and the vaccination threshold is obtained by solving the optimization problem.

We expect a faster spread of the disease in the networks with a bigger clustering and closeness coefficients and lower betwenness, from the ideas exposed in the last paragraph.

Comparing Fig. 3 with Tables 1 and 2 we can check that the velocity of propagation of the disease and the maximum number of infected agents is bigger in networks (b) and (d), where the basic reproduction number and the clustering and closeness coefficients are bigger and the betwenness is smaller. Moreover, the size of the epidemic in (d) is greater than the size for networks (a) and (c).

On the other hand, the final number of vaccinated is bigger in networks (a) and (c), almost doubling the number of vaccinations in network (d). The former are the most homogeneous in terms of degree and mixing modeling, meaning that the vaccination rate of different individuals are similar, while (d) have a big number of degree 1 individuals and a considerable proportion of high-degree nodes, whose vaccination rates have great variability. The Bimodal and Power Law networks seem to best reflect the interactions between people [1, 32, 33, 43] . In those cases, we observe that the spreading of the disease is faster and infect a bigger proportion of the network, due the presence of hubs, and the degree-dependent vaccination effect is not enough to contain the epidemic. Let us note that for networks (b) and (d) we have R 0 ∼ 10, and the herd immunity is reached when about 50-60% of the population is infected or vaccinated. In a complete network, we need about 90% of the population recovered or vaccinated. Also, for R 0 ∼ 3-4, we need about of 70-75% in the complete network, while in networks (a) and (c) the vaccinated and recovered agents surpass the 90% of the population.

Finally, we present the results of the numerical integration for four networks with the same R 0 = 3.75. The degree distribution is the same as the one used in the previous plots, but with a fine-tuned selection of parameters: λ = 5 for the Poisson, λ = 3, L = 8, p = 0.73 for the Bimodal, the degree of the Regular network is 6, and for the Power Law we take α = 2 with exponential cutoff κ = 20.

In this case, the Power Law (d) network is very different to the other three. Its coefficients of centrality indicate that a more rapid spread of the disease is expected. Here, the role of the hubs is more relevant because the quotient ψ (1)/ψ (1) is the same. This is shown in the Table 3 . We can also expect from Fig. 4 that vaccinating the hubs, which will occur at high rate, will somehow cleave the network, stopping the propagation of the disease.

Effectively, in Fig. 5 we observe that the final size of the epidemic is substantially lower than in the other networks, and the final number of vaccinated is also much Fig. 4 Instances of graphs when N=80 for the same R 0 lower, indicating that in this case the network used to model the reality plays a central role.

From the presented instances, we can also observe that acquaintance vaccination is more effective in all the cases, arriving to lower numbers of infected and with a less number of vaccinated, except for the case of the Regular network (c), where both vaccination strategies perform the same because the degree is equal for all the nodes.

We considered a generalization of the SIR model on a large configuration model adding a vaccination strategy. The variables describing the dynamics were determined by a Markovian contact process on a diverse population where individuals have different degrees in a network of random connectivity, the vaccination mechanism depending on the number of possible contacts of an individual.

We derived large graph limits for the evolution of the epidemic in this context leading to a system of seven differential equation that describe the SIRV dynamics.

We then solved the associated optimal control problems for a general vaccination rate, proving existence and uniqueness of the solutions under mild assumptions on cost functions. We also characterized the optimal solutions as threshold type.

For this type of strategies, we computed the impact of the vaccination process in terms of epidemic sizes, finding a fixed point equation that describes the impact of the vaccination intensity and the time on which it is developed through the function g (or ψ) and depending on the proportion of I S edges and the vaccination rate.

Finally, we studied four particular networks and their centrality coefficients, relating it with the epidemic indicator R 0 . Given a maximum vaccination rate and for fixed disease parameters, we solved numerically the system of equations for the cases of degree-proportional and constant vaccination, on networks with the same mean degree first and then with the same mean excess degree. We observed that in the networks that describe better the interaction of the individuals, the epidemic level of infection can be significantly different from in homogeneous models.

In this section we describe the theoretical background of the proofs and modeling methods. We first introduce the stochastic differential equations that define the process and then explain the scaling used in order to get the fluid limit result. We apply a tightness-uniqueness strategy to prove convergence and characterize the fluid limit thanks to a martingale representation.

Inspired by [14] and [20] we will represent the behavior of our dynamic as a process which is solution of a system of a stochastic differential equations derived from Poisson Point Measures (PPM).

We will use three different PPM for each event which modifies the quantities we are interested in: an infection, a recovery, or a vaccination. We need to identify the rates of these events and how to update the measures on the graph.

Suppose an event occur at time T , and let us analyze the first case, an infection. For that, it is convenient first to consider

the rate of infection of a given k-degree individual at time T . They will have their halfedges connected according to the quantities μ T and distributed following a multivariate hypergeometric distribution. We denote

Finally, given k, j, l and m, we have to update the measures μ I S T , μ RS T and μ V S T choosing the infected, recovered and vaccinated individuals who will be connected to the newly infected. In order to do that, we draw three vectors u = (u 1 , . . . ,

indicating how many links each I , R or V node has with the newly infected. We consider U = n∈N (N 0 ) n and for each μ ∈ M F (N 0 ) and n ∈ N we define U ⊇ U(μ, n) := u = (u 1 , . . . , u μ,1 ) :

where ζ i (μ) := F −1 μ (i) is the degree according μ of the i-th node. Similarly, we define v, w ∈ U depending on the measures μ RS T and μ V S T .

Thus, the number of edges of type I S, RS or V S will be given respectively by

We define

and

Then, we update our measures as follows, introducing some notation:

Another event in consideration is a recovering. Here we choose uniformly an infected i and set:

This happens with probability 1/I T − . The last event is vaccination. The corresponding rate is π t N S t . We remark the strong dependence on the degree of the individual, because it is more probable that a higher connectivity node to be vaccinated first. More precisely, the probability that the new vaccinated has degree k is

. Once we draw the vaccinated individual, and supposing their degree is k, we update the measures as follows:

Now we introduce three Poisson Point Measures that will be very useful to describe the M F (N 0 )-valued stochastic process (μ t ) t≥0 . For a similar point of view, see [14] or [20] .

The first one will provide us the possible instant in which an infection occurs. We define d N 1 (s, k, θ 1 , j, l, m, θ 2 , u, θ 3 , v, θ 4 , w, θ 5 ) as a product measure on 3 , where ds and dθ are Lebesgue measures and dn are counting measures on N 0 or U, accordingly.

The degree k infected agent will be connected with j infected, l recovered and m vaccinated agents, drawn according u, v and w as we explained above.

We also have d N 2 (s, i) on E 2 = R + × N a PPM with intensity γ for the recovering process. That is, for each atom we have associated a possible recovering time s and the identification number i of the new recovered.

The last PPM, d N 3 (s, k, θ 1 , j, l, m, θ 2 , u, θ 3 , v, θ 4 , w, θ 5 ) is defined in R + × E 3 where E 3 = E 1 and it is very similar to the first one. It assigns a mass to each possible time s where a degree k vaccinated agent is connected with j infected, l recovered and m vaccinated agents, drawn according u, v and w.

In all the cases, the auxiliary variables θ are useful to take into account the rates in this integral representation. In order to simplify notation we will not write the dependency on the variables, and consider the following indicator functions to represent the rates:

Now that we have clarified the evolution of the measures according to the events that may occur, we are ready to write an integral form for this evolution in terms of the Poisson Point Measures, for example for the second coordinate

Doing the same for the four coordinates, we can write the system of Stochastic Differential Equations: N 1 , N 2 , N 3 there exists a unique strong solution to the system (30) in the Skorokhod space D(R + , (M F (N 0 )) 4 ).

Proof First note that all the measures are dominated by the expectation of μ S 0 + μ I S 0 + μ RS 0 + μ V S 0 and the supports are bounded on the positive integers. The proof can be completed in the same way as in [50] .

Inspired by the techniques developed in [14] and [20] we write a renormalization of the system when the number of individuals is n and the number of edges is proportional to n. We observe that the intensity of the jump process has the same order, and deduce the scaling for the fluid limit renormalization. We prove the convergence of the solution of the finite case system of equations to the solution of (30) in the weak sense of the Skorokhod space [23] .

Let us consider four sequences of measures indexed by n ∈ N, (μ n,S ), (μ n,I S ), (μ n,RS ) and (μ n,V S ) satisfying the system of equations (30) 

The proof is divided in five steps. The first one consists in the tightness of the renormalized process regarding several criteria of convergence in the Skorokhod space. The second step is the uniqueness of the solution. The third step prove that the renormalized process satisfies asymptotically the deterministic system 1, showing that the limit is a good approximation for large populations. In step 4 we prove that the limit satisfy this equation. This is done for the measure of the susceptibleinfected edges. In step 5 we prove the convergence of the other measures regarding the solution of a transport equation. Let us consider, for any ε ≤ 0 and A > 0, the closed set of M F (N 0 ),

and M 0+,A = ε>0 M ε,A . Some topological properties of this can be found in the appendix of [14] or [20] . We suppose that μ (n) 0 converges to μ 0 and that μ (n) 0 ∈ M 4 0,A for any n, with μ I S 0 , χ > 0. We also make the assumption π t μ t , 1 ≤ ν k μ S 0 (k)ξ(k) < C μ S 0 , χ j for j = 3 just in order to keep the hypothesis of finite 5-th. moment, although we can change it to finite j + 2-th. moment if necessary imposing more restrictive conditions but deriving the result making the changes needed.

In order to prove (ii), since lim ε →0 t ε = ∞, is enough to prove the result in D([0, t ε ], M 4 0,A ) for ε sufficiently small. From now on, we take 0 < ε < ε < μ I S 0 , χ .

Step 1: Tightness of the renormalization. Take (μ n ) n∈N , t ∈ R >0 and n ∈ N. By assumptions, we have: 

This implies that the sequence μ (n) t is tight for each t. By the criterion of convergence of measure valued processes proposed by Roelly [48] we have to prove that, for each test function f ∈ B b (N), μ (n),S , f , μ (n),I S , f , μ (n),RS , f , μ (n),V S , f n∈N is tight in D(R >0 , R 4 ). We present here the calculations only for μ (n),I S , f because the others are similar or simpler. Since we have a semimartingale decomposition, applying the Rebolledo criterion for weak convergence of sequences of semimartingales, we have to prove that both the finite variation part, and the quadratic variation satisfy the Aldous criterion. We want to prove that, for all θ > 0 and η > 0 there exist n 0 ∈ N and δ > 0 such that for all n > n 0 and for all stopping times S n and T n with S n < T n < S n + δ we have P(|A 

For the finite variation condition (34), we take the following bound:

E |A Since j+l+m≤k p n s ( j, l, m|k − 1)2 j is twice the mean number of edges with the infected population conditioned to having degree k, this number is bounded by k, and using the definitions of λ n , π and p we have that:

Then, applying Markov's inequality:

which is smaller than θ if δ is small enough. We bound the quadratic variation of the martingale reasoning analogously, 

In order to prove C 

as long as s ≤ τ n ε and n ≥ 1/ε. Let us remember that k S n t i is the number of edges of the i-th node connecting with the susceptible population at time t, for a population of size n.

Additionally, for all the possible draws u ∈ U( j, μ n,I S s ) we have 

This last expression tends to zero because of the weak convergence of μ 

Classes of small-world networks

Infectious Diseases of Humans: Dynamics and Control

Acquaintance vaccination in an epidemic on a random graph with specified degree distribution

Approximating the epidemic curve

The jamming constant of uniform random graphs

Convergence of Probability Measures

Random graphs

Centrality and network flow

Viscosity solutions of Hamilton-Jacobi equations and optimal control problems

Graphs with specified degree distributions, simple epidemics, and local vaccination strategies

Efficient immunization strategies for computer networks and populations

Erratic flu vaccination emerges from short-sighted behavior in contact networks

Large graph limit for an SIR process in random network with heterogeneous connectivity

Ordinary differential equations, transport theory and sobolev spaces

A mean-field game analysis of SIR dynamics with vaccination

Discrete mean field games: Existence of equilibria and convergence

Random Graph Dynamics

Eb-devs: A formal framework for modeling and simulation of emergent behavior in dynamic complex systems

A microscopic probabilistic description of a locally regulated population and macroscopic approximations

Imitation dynamics of vaccination behaviour on social networks

Long-standing influenza vaccination policy is in accord with individual self-interest but not with the utilitarian optimum

The Theory of Stochastic Processes I

Exploring network structure, dynamics, and function using networkx

Using network science to propose strategies for effectively dealing with pandemics: The covid-19 example

Optimal vaccination schedules in a deterministic epidemic model

Transitivity in structural models of small groups

The probability that a random multigraph is simple

Law of large numbers for the SIR epidemic on a random graph with given degrees

Mathematics of Epidemics on Networks: From Exact to Approximate Models

Global optimal vaccination in the SIR model: Properties of the value function and application to cost-effectiveness analysis

The web of human sexual contacts

How viruses spread among computers and people

The effect of constant and pulse vaccination on SIR epidemic model with horizontal and vertical transmission

Infection dynamics on scale-free networks

The impact of imitation on vaccination behavior in social contact networks

Convergence of the forward-backward sweep method in optimal control

A note on a paper by Erik Volz: SIR dynamics in random networks

A primer on the use of probability generating functions in infectious disease modeling

Epidemic outbreaks in complex heterogeneous networks

Spread of epidemic disease on networks

The structure and function of complex networks

Epidemic spreading in scale-free networks

Epidemic processes in complex networks

Dynamics of epidemic spreading with vaccination: impact of social pressure and engagement

Markov Decision Processes: Discrete Stochastic Dynamic Programming

Continuous Martingales and Brownian Motion

A criterion of convergence of measure-valued processes: application to measure branching processes

Effectiveness of realistic vaccination strategies for contact networks of various degree distributions

Modèles particulaires stochastiques pour des problèmes d'évolution adaptative et pour l'approximation de solutions statistiques

Contrôle optimal: théorie & applications

An epidemic model to evaluate the homogeneous mixing assumption

SIR dynamics in random networks with heterogeneous connectivity

Efficient vaccination strategies for epidemic control using network information

Stability analysis and optimal vaccination of an SIR epidemic model

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

nodes susceptible, infected, recovered and vaccinated, respectively. We take the scaling μ We assume that the sequences of initial conditions converge weakly in M F (N 0 ) to μ S 0 , μ I S 0 , μ RS 0 and μ V S 0 when n goes to infinity. We obtain the renormalized system:Let us define where the finite variation is given byand the associated martingale is square integrable with quadratic variation,

We first calculate the infinitesimal generator L of our process, and we write the Levy's martingale with φ = μ, f and φ 2 . Then we apply the integration by parts formula [47] , and identifying the martingales in the expression, we rearrange the terms in order to get the quadratic variation. For a detailed proof see [20] .Our fluid limit result may be proved in the same way as the proof of the main theorem of [14] but we add it for completeness.Let μ be a limit point in C(R + , M 4 0,4 ) of the sequence of stopped processes and let (μ (n) ) n∈N be a subsequence that converges to μ, denoted with only one over-script to simplify notation. Since the limit is continuous, the convergence is uniform over compact sets of the positive reals.Define, for all t ∈ R + and f ∈ C b (N) the applicationssuch that (3) can be read asStep 2: Uniqueness of the solution inThe second step consists in proving the limit values are the unique solution of (3). The strategy will be to prove that the total measure and the first and second moments of two solutions are equals and then prove that the generating functions of those measures satisfies a partial differential equation that admits a unique solution in a weak sense.Due to extension by regularity, it is enough to prove the uniqueness inTake μ i = (μ S,i , μ I S,i , μ RS,i , μ V S,i ) for i = 1, 2 two solutions of (3) in this space with the same initial condition and defineNote that, for all t ∈ [0, T ) and i = 1, 2, we have N S,i t ≥ N I S,i t > ε and thereforeSince μ i are solutions of (3), we have, for j = 0, . . . , 3 and taking f = χ jOne can reproduce similar computations for the other quantities:With this, ϒ satisfies a Gronwall's type inequality which implies that is identically 0 for all t ≤ T . Then, for all t < T and for j = 1, 2, 3, we have It remains to prove the uniqueness for the other 3 measures. The method that we will use to prove μ I S,1 = μ I S,2 can be used for the rest.We consider the generating functions:for any t ∈ R + , i = 1, 2 and η ∈ [0, 1). Let us defineUsing f (k) = η k in the second equation of (3) and after some basic computations we getNow, H (t, η) is continuously differentiable with respect to time and it is well defined and bounded in [0, T ]; and K t is piece-wise continuous in L 1 and also it is well defined and bounded on [0, T ]. Further, H and K do not depend on the solution we choose, because we already have μ S,1 = μ S,2 and p I ,In view of the regularity of H and K it is known that this equation admits only one solution in a weak sense (see last section in [15] and for all η ∈ [0, 1). Since both measures have the same mass, we have μ I S,1 = μ I S,2 .We can use similar arguments in order to prove that μ V S,1 = μ V S,2 and μ RS,1 = μ RS,2 .Step 3: μ (n) satisfies asymptotically the deterministic system (3). Let us remember that, for each f ∈ C b (N), we can write:In order to characterize the limiting values, for each n ∈ N and for all t ≥ 0, we havewhere n, f ·∧τ n ε vanishes in probability and uniformly in t over compact time intervals. We can take bounds in a similar way as in Step 1 in order to get that:which implies the sequence (M (n),I S, f t ) n∈N vanishes in probability and in L 2 , and therefore in L 1 by Cauchy-Schwartz.On the other hand, the finite variation part can be split in two: one considering the simple edges between the newly infected node and the infected population, and a second part regarding multiple edges, that we know is expected to vanish as the size of the population grows. Formally,where τ j f (k) := f (k − j) for all f : N → R and ∀k ∈ N. Now we introduce some notation for the proportions of edges the newly infected agent has, discarding the edge involved in the infection process. It is important here to make a difference between the term that comes from the infection from the one that comes from the PPM modeling the vaccination process, because in this case we do not assume a priori that there is at least one infected neighbor. We define, for each t > 0 and n ∈ N 0 ,Let us remember that In the case of the infection process, we also define, for all the j, l, m such that j + l + m ≤ k − 1, and for all n ∈ N,and for the vaccinations,the probabilities of the multinomial variables counting the quantities of each type of neighbors that will has the newly infected or vaccinated, respectively.We can write whereThus, if we consider the differenceswe can bound:Since the multinomial term is a good approximation of the multivariate hypergeometric as n goes to infinity, the last expression tends to zero due to dominated convergence. On the other hand,Putting all the bounds together, we can conclude that μ (n),I S , f converges in probability uniformly over compact intervals.We are considering the sequence (μ (n) ·∧τ n ε ) n∈N and we already proved that its limit in the closed set M 4 0,A is μ, we want to prove the same for the non stopped sequence. According to the Skorokhod representation theorem there exists a subsequence on the same probability space of μ whose marginal probability distributions are the same as those of the original sequence such that μ is the almost sure limit. With an abuse of notation, we will denote (μ (n) ·∧τ n ε ) n∈N this subsequence. The mappingsAccording to Lemma 1 we have that, for p ≤ 5, p : D(R + , M ε,A ) → D(R + , R) which assigns ν · → ν · , χ p is continuous.Using this, and that the quotient Let us define t ε = inf{t ∈ R + : N I S t ≤ ε }. We do not know this number to be deterministic, but we can say that:Therefore, splitting in the following way:we have, from the bounds and estimations we made in Step 3, that I S, f ·∧τ n ε ∧T (μ (n) ) is bounded for the fourth moment of μ (n) . Since μ (n) 0 → μ 0 and using (55) , the first term in (56) converges in L 1 and in probability to zero. On the other hand, the continuity ofconverges to I S, f (μ) and therefore,·∧t ε ∧T (μ). So, this convergence and (55) implies that, the second term converges toStep 3 again, and the estimations done in it, we can conclude this sequence also converges in probability to zero. Therefore, we have that μ I S is a solution to the system (3) on the interval [0, t ε ∧ T ].If either μ RS 0 , χ > 0 or μ V S 0 , χ > 0, then we could apply similar techniques with both. If not, the result can be immediately deduced because for all t ∈ [0, t ε ∧ T ], μ (n),I S t , χ > ε and the terms p n t ( j, l, m | k − 1) and q n t ( j, l, m | k) are negligible when l or m are positives.So, μ is almost surely the unique continuous solution of the deterministic system (3) in [0, t ε ∧ T ], which implies t ε = t ε and the convergence in probability of (μ (n) ·∧τ n ε ) n∈N to μ holds, uniformly on the interval [0, t ε ], due to the continuity of μ. In order to prove the convergence in the Skorokhod space, for η > 0, we write:Using the continuity of f and the uniform convergence in probability that we have proved, the first term in the last expression converges to zero. In order to show that the second term vanish, we can reproduce the bounds taken in Step 2 of this proof and apply Doob's inequality. Finally, since P(τ n ε > T ∧ t ε ) → 1 we have that the three terms goes to zero.Hence, due to uniqueness proved in Step 2, the original sequence (μ (n) ) n∈N converges.Step 5: The convergence of the other measures What we have done for the infectedsusceptible connectivity measure can be also done for the recovered and vaccinated measures in much the same way. For the susceptible connectivity measure, one can reason in the following way. If we consider the renormalized equation 31 and we take limit in n, the sequence (μ (n),S ) n∈N converges in D(R + , M 0,A ) to the solution to the transport equationthat can be solved as a function of p I and π , for any test function f ∈ C 0,1 b (N×R + , R) with bounded derivative respect time variable [15] .If we take f (k, s) = ϕ(k)e − t−s 0 rkp I u +π u (k)du we obtainas the first equation of (3) establishes. The proof is finished.Lemma 1 For any p ≤ 5, the map p : D(R + , M ε,A ) → D(R + , R) that assigns (ν . ) → ν . , χ p is continuous.Proof The proof can be obtained by following the steps of Lemmas 1-5 in the appendix of [14] .