Abstract
This paper introduces high-order explicit Runge–Kutta numerical schemes in metric spaces. We show that our approach reduces to the corresponding Runge–Kutta schemes if the ambient space is Hilbert. We apply these schemes to compute the Nash equilibrium in a mean field vaccination game. Numerical simulations show improvement in the speed of convergence towards the Nash equilibrium; the numerical scheme has high order (2–4) in time.
Similar content being viewed by others
1 Introduction
1.1 State of the art and motivation
The games with a continuum of agents have been widely studied during the last decade thanks to the mean field games (MFG) theory, introduced by Lasry–Lions [34,35,36] and Huang–Caines–Malham [26,27,28]. An application of this framework is the modeling of the vaccination decision in a population where an individual has the choice between only two pure strategies: to get vaccinated or not.
This decision is a part of the behavioral epidemiology, which takes into account the interplay between human behavior and spread of infectious diseases (see [38] for a detailed presentation and [30, 40, 47] for other examples). In particular, risk perception and spread of information (on vaccination risk or probability of infection) can influence at a large scale the population behavior and thus the spread of the disease. As detailed in [38], a well-known example is Measles–Mumps–Rubella (MMR) vaccination coverage decrease between 1998 and 2003. This caused a decline in herd immunity and a measles resurgence after the publication of an hypothesis linking MMR vaccination and autism. Behavioral responses to the threat of a pandemic event can influence the spread of the disease in several ways: people can choose preventive protection (by vaccination as presented below) or limit the disease spread by isolation and quarantine (e.g., school closures). Since then, many epidemiological models have been developed to include individual response (see [47] for a review); on the contrary see [1, 6, 39, 42] for the treatment of compulsory (as opposed to voluntary) vaccination.
We use the family of compartmental models in order to describe the spread of the disease in the population. This supposes that all individuals are the same and the population is divided into several exclusive classes; the transmission between two classes is modeled through a time-dependent rate. As each compartment represents the status of an individual, it means that the structure of model with immunity is not the same as that of a disease without immunity. \({ SIR}\) terminology usually represents a disease with immunity against re-infection, indicating that the passage of an individual is from the susceptible class S to the infectious class I and then to the recovered class R. If the re-infection is possible, then the terminology used can be \({ SIS}\), which means that an individual in the infected class I goes back to the susceptible class S. If temporary immunity is given by the disease, the compartmental structure can be \({ SIRS}\), meaning that after infection, individual stays in R class before going back to the susceptible class. Other well-known compartmental models can include an exposed class (E) that represents the period between being infected and becoming contagious, like in the \({ SEIR}\) or \({ SEIS}\) models. The evolution of the size of the compartments is given by ordinary differential equations (for more details see for instance [10, 15, 25]).
In order to model individual vaccination decision in the spread of the disease, vaccination coverage can depend on available information (see [11, 18, 19]) or be determined by a vaccination game. The evolution of the vaccination coverage could be determined by an imitation dynamics depending (or not) of the available information (see [3, 16, 16, 17, 32, 41]), or be the solution of a Nash equilibrium determined at the beginning of the game, as described below. The Nash–MFG equilibrium corresponds in this case to a strategy where nobody has interest to change his own vaccination strategy. Therefore the existence and the determination of an equilibrium is central.
Before the development of the Nash–MFG theory, some work have been proposed to couple game and epidemiology. Three well-known methods that introduced specific models are presented in the next paragraph with their advantages and limitations. The compartmental model used in [5] includes five classes (susceptible, exposed, infected, recovered, vaccinated), the dynamics is presented in Eq. (1) where the population is divided into the following classes: S (susceptible), I (infected), R (recovered) and V (vaccinated). In this case preventive vaccination for the smallpox is available: the model is used to determine the probability of an individual to choose preventive vaccination. To do so, explicit individual payoffs for the two pure strategies (vaccination before or after the outbreak) are defined. Then the existence and uniqueness of the Nash equilibrium is obtained using properties of such payoffs.
Another method introduced in [5] uses the same methodology as previously but on the stationary state of the system (2) (or equivalently system (2bis) that distinguishes immunity given by the disease and immunity given by the vaccine) because of the possible non extinction of the disease. This allows to determine the probability of vaccination when disease is already present in the population with a compartmental system where vital dynamics is considered. As in previous case, the probability p to vaccinate is not time dependent and properties of payoff are used to prove existence and uniqueness of the equilibrium. Finally, solutions for other incidence rate (term \(\beta I\) is replaced by a function \(\lambda (I)\)) are given in [43] by using Markov decision process. Analytical solutions for the vaccination probability at the Nash equilibrium can be found with these approaches but it supposes constant probability of vaccination.
The third and last model being detailed here describes numerical approach when two \({ SEIR}\) models (see for example (1)) are used in order to separate two age classes (see [24]) to obtain the vaccination probability of each age class. The evolution of the probability of vaccination can be modeled by introducing a rule of thumb (see [3, 4]) to study the epidemiological stationary state and the role of parameters like vaccination and infection costs. But reaching the stationary state of the epidemiological system can take a long time as explained in [29]. To avoid this problem, [20] considers a \({ SIRS}\) model where the transmission rate \(\beta \) is a periodic function, causing the same property for the incidence. Economic approach (see [21]) based on the epidemiological model (3) (where notations are the same as above and \(\mathcal {V}(t)\) is the instantaneous vaccination rate) allows to determine the vaccination and no vaccination region (depending on the state of the disease) by equalization of the two payoffs but no proofs of optimality (existence or uniqueness) are provided.
Within the MFG framework, article [33] proposes analytic solutions in some particular cases. Furthermore, finding an equilibrium becomes more complicated if the model has specific behaviour with network interaction or epidemic spread for example [14, 23]. Thus finding the equilibrium can quickly become time consuming. One commonly used method is based on finding a fixed point of some function, which can be formulated as the attractor of an evolution system.
The existence of an equilibrium in the general case of MFG is of utmost interest (see for instance [7, 9] for an entry point to this literature). Several approaches were developed, inspired by the general framework of gradient flows (see [31]) where the procedure is iterative.
1.2 Mathematical framework
Let us recall that in a gradient flow framework (which is not necessarily related to a MFG model !) the equation to be solved is written symbolically (see [2] for an introduction):
where the function F is given but \(\nabla F\) cannot be computed. The numerical counterpart of (4) is the famous JKO numerical scheme (introduced by Jordan et al. in [31]) which, for \(\xi _k\) an approximation of \(\xi _{\tau _{k}}\), can be written as:
This scheme can be seen as the variational implicit Euler scheme and has been used for the heat equation as a gradient flow of the entropy in the Wasserstein space. Second order in time methods have been proposed very recently (see [37]).
On the other hand, a MFG equilibrium has no associated function ‘F’ and it is likely that no such function exists. For each \(\tau , \xi _\tau \) is a probability law. Rather, the search for an equilibrium takes the form of a fixed point:
In practice, inspired by the fictitious game (see [12]) and Best Reply (see [8, 9]) procedures, a new algorithm has been introduced in [44, Eq 3.2] and [46] with the following form:
In the limit \(\tau \rightarrow 0\) the Eq. (7) describes a curve which under some assumptions on \(\mathcal {P}\) will be called a solution of the evolution equation:
On the other hand, the MFG framework provides a natural function \(\mathcal {P}(\cdot , \cdot )\) even when \(\xi \) does not belong to a Hilbert space but only to a metric space. Thus we are led to consider evolution equations in metric spaces with semi-explicit numerical schemes (see [9, 46]).
1.3 Scope and structure of the paper
The purpose of this work is to introduce variational Runge–Kutta explicit methods of high order in a metric space using a generalization of linearity (presented in [13]).
These approaches are applied to the control of epidemic spread with voluntary vaccination as in [44].
The paper is structured as follows: Sect. 2 is dedicated to the Runge–Kutta method recalling its general form (Sect. 2.1) and then the introduction of variational Runge–Kutta methods (Sect. 2.2). The equivalence of the two approaches in a Hilbert space is proved in Sect. 2.3. Then Sect. 3 gives the epidemiological model used in the numerical application given in Sect. 4. An economical application is proposed in Sect. 5. Some considerations are discussed in Sect. 6.
2 Runge–Kutta methods
2.1 In a Hilbert space
Let \(\mathcal {H}\) be a Hilbert space, \(f: [0, T] \times \mathcal {H} \rightarrow \mathcal {H}\) a regular function and \(y: [0,T] \rightarrow \mathcal {H}\) satisfying (in some sense to be specified):
The Runge–Kutta approach is an iterative method using temporal discretization in order to obtain a numerical solution of an ordinary differential equation of type (9). The time horizon T is supposed to be finite and can be discretized in \((N_T+1)\) time instants \( \tau _0 = 0, \tau _1 = {\varDelta }\tau , \tau _2 = 2 {\varDelta }\tau , \ldots , \tau _{N_T} = T\) where \({\varDelta }\tau \) is the time step. Let \(y_{k} \in \mathcal {H}\) be an approximation of \(y(\tau _k)\).
Runge–Kutta method is based on the evaluation at intermediate points in time. Equation (10) presents the general form of the Runge–Kutta method for (9).
A method is defined by the values of the coefficients \(b_i, c_i\) and \(a_{i,j}\), often presented in the Butcher tableau (see table in Fig. 1). Consistency of such a method is ensured if and only if \(\sum _{i=1}^s b_i = 1\).
We will present below some particular schemes that will later be formulated in a metric space.
2.1.1 Explicit Euler scheme
The explicit Euler (denoting EE) scheme is a Runge–Kutta method with only one step. It is defined in Fig. 2 by the Butcher tableau or equivalently by Eq. (11). In this case, this scheme is of order one in time.
2.1.2 Heun scheme
Heun scheme uses an approximate value \(p_1\), in order to compute \(y_{k+1}\), see the Fig. 3. This scheme is of order two in time.
2.1.3 Explicit Runge–Kutta 3 scheme
In the same way, we introduce the Runge–Kutta 3 (denoted RK3) method. It calculates two intermediate values (\(p_2, p_3\)) to obtain \(y_{k+1}\). See Fig. 4 for the value of coefficients in the Butcher table and Eq. (13) for the explicit equations. In this case, the method is of order three in time.
2.1.4 Explicit Runge–Kutta 4 scheme
The last method is the celebrated Runge–Kutta 4 (denoted RK4) where three intermediate values (\(p_2, p_3, p_4\)) are computed to determine \(y_{k+1}\). The value of coefficients is presented in Fig. 5 and its application to Eq. (9) is given in Eq. (14). This method is of order four in time.
2.2 Variational approach
Unfortunately, none of the schemes presented in Sect. 2.1 can be used for finding a MFG equilibrium. Firstly because the problem is not really presented that way, and secondly because the space of the unknowns is, in general, not a vector space but a metric space. Therefore we need first to generalize the Runge–Kutta schemes to metric spaces and then to apply them to the MFG setting.
Let \((\mathcal {A}, d)\) be a metric space. The Definition 1 below gives the properties of a geodesic space (see [37] for additional details).
Definition 1
A curve \(\gamma : [0,1] \rightarrow \mathcal {A}\) is called a (constant speed) geodesic provided that \(d(\gamma (s), \gamma (\tau )) = |s- \tau | \cdot d(\gamma (0), \gamma (1))\) and the space is called geodesic if, for any couple of points \((X,Y) \in \mathcal {A}^2\), there exists at least a geodesic \(\gamma \) connecting them, that is, \(\gamma (0) = X, \gamma (1) = Y\).
From now on we assume that \((\mathcal {A}, d)\) is a geodesic space. Consider \(\mathcal {P}: \mathcal {A} \times \mathcal {A} \rightarrow \mathbb {R}\) a function with the following properties:
- (H1) :
-
\(\forall n \le 4, \forall r>0, \forall \xi , Y_1, \dots , Y_n \in \mathcal A\), the set of vectors \(\{ (\mathcal {P}(z, Y_1), \dots ,\mathcal {P}(z, Y_n)); z \in {\mathcal A}, d(z,\xi ) \le r \} \) is compact as a subset of \(\mathbb {R}^n\),
- (H2) :
-
for any point \(Y \in \mathcal {A}\) and any \(X_0, X_1 \in \mathcal {A}\) a constant speed geodesic \(\gamma \) exits such that \(\gamma (0) = X_0, \gamma (1) = X_1\) and the function \(\tau \mapsto \mathcal {P}(\gamma (\tau ), Y)\) from [0, 1] to \(\mathbb {R}\) is linear.
Examples:
-
If \(\mathcal {H}\) is a Hilbert space, \(f: \mathcal {H} \rightarrow \mathcal {H}\) a function then \((x,y) \rightarrow \left\langle x,f(y)\right\rangle \) satisfies hypothesis (H1) and (H2), see proof of Proposition 2.2.
-
Let \(\mathcal {A}\) be the set of probability measures on [0, T]; \(\mathcal {F}\) be the set of Lipschitz functions with constant less than 1 and the metric of weak convergence d is defined, for \(a_1, a_2 \in \mathcal {A}\) by (see [45]): \(d(a_1, a_2) = \sup _{f \in \mathcal {F}} | \int _0^T f(t) a_1(dt) - \int _0^T f(t) a_2(dt) |\). For \(Y \in \mathcal {A}, t \in [0,T]\) we define \(g_Y(t) = \left[ \sin (t) \int _0^T \cos (s) Y(ds) \right] ^2\) and for \(z, Y \in \mathcal {A}, \mathcal {Q}(z,Y) = \int _0^T g_Y(t) z(dt)\).
Let \((z_k)_{k \in \mathbb {N}}\) be a sequence of elements in \(\mathcal {A}\) converging weakly to \(z \in \mathcal {A}\) then, for \(f \in \mathcal {F}, \int _0^T f(t) z_k(dt)\) converges to \(\int _0^T f(t) z(dt)\). As \((\mathcal {A},d)\) is compact (see [45]) and the application \(z \in \mathcal {A} \rightarrow (\mathcal {Q}(z,Y_1), \dots ,\mathcal {Q}(z,Y_n)) \in \mathbb {R}^n\) continuous, \(\mathcal {Q}\) satisfies (H1).
Consider \(a_1\) and \(a_2 \in \mathcal {A}\) and for \(s \in [0,1], \gamma _s = sa_1 + (1-s) a_2\). Then \(d(\gamma _{s}, \gamma _{\tau }) = \sup _{f \in \mathcal {F}} | \int _0^T f(u) \gamma _{s}(du) - \int _0^T f(u) \gamma _{\tau }(du) | = \sup _{f \in \mathcal {F}} | \int _0^T f(u) (\tau -s)a_1(du) - \int _0^T f(u) (\tau -s)a_2(du)| = |\tau -s| \sup _{f \in \mathcal {F}} | \int _0^T f(u) a_1(du) - \int _0^T f(u) a_2(du)| = |\tau -s| \cdot d(\gamma (0),\gamma (1))\). It shows that \(\gamma \) is a (constant speed) geodesic. Then, the function \(s \rightarrow \mathcal {Q}(\gamma _s, Y)\) of the form \(s \rightarrow s \int _0^T g_Y(u) a_1(du) + (1-s) \int _0^T g_Y(u) a_1(du) \) is linear.
-
Note that previous example also works with the total variation distance.
The evolution equation that we want to solve has the form of Eq. (8) (but note that this is only a formal expression as \(\frac{d}{d \tau } \xi _{\tau }\) does not have a well-defined meaning in a general metric space, see [46]).
This section develops a variational Runge–Kutta method in a metric space. We propose a method that does not assume any vector calculus but uses the hypothesis (H2) which defines the notion of linear application in a metric space (see also [13, Section 8, Definition 8.1, p. 480] for more details). As \({\varDelta }\tau \) tends to zero Eq. (7) describes a curve. The purpose of the scheme is to describe this evolution faster and more precisely that is to say to obtain high order schemes.
2.2.1 Variational explicit Euler scheme
The first scheme of that form was introduced in [44, Eq 3.2]. Formally, we define the variational explicit Euler scheme (VEE) by:
2.2.2 Variational Heun scheme
The variational Heun scheme (VH) introduces, as in vector space, an intermediate value \(\widetilde{\xi }_{k+1}\) (Eq. 16a) and uses it to compute \(\xi _{k+1}\) (Eq. 16b).
Note that two minimizations are required in order to obtain \(\xi _{k+1}\).
2.2.3 Variational Runge–Kutta 3 scheme
The variational Runge–Kutta 3 scheme (VRK3) is defined by Eq. (17). Two intermediate values are computed in this case [\(\widetilde{\xi }^1_{k+1}\) with Eq. (17a) and \(\widetilde{\xi }^2_{k+1}\) with Eq. (17b)] in order to obtain \(\xi _{k+1}\) using Eq. (17c).
2.2.4 Variational Runge–Kutta 4 scheme
The variational Runge–Kutta 4 scheme is presented in Eq. (18). Three intermediate values [\(\widetilde{\xi }^1_{k+1}\) with Eq. (18a), \(\widetilde{\xi }^2_{k+1}\) with Eq. (18b) and \(\widetilde{\xi }^3_{k+1}\) defined by Eq. (18c)] are computed in order to obtain \(\xi _{k+1}\) by using Eq. (18d).
Note that four minimizations are needed in order to obtain \(\xi _{k+1}\).
See Sect. 4.2 for a numerical application and the study of the order in time of these schemes.
2.3 Property of variational scheme
In this subsection we prove that using hypotheses (H1) and (H2), the variational schemes introduced above are well-defined (Proposition 2.1) and correspond, when the ambient space is Hilbert, to the Runge–Kutta method presented in Sect. 2.1 (Proposition 2.2).
Proposition 2.1
Under the hypothesis (H1) and (H2), the Eqs. (15), (), () and () defining respectively the VEE, VH, VRK3 and VRK4 schemes admit a solution, i.e., the schemes are well defined.
Proof
As the argumentation is very similar for all schemes, we present the general case where \(\xi _{k+1} \in \text {argmin}_{\eta \in \mathcal {A}} \{ \mathscr {F}(\eta ) \}\) with \(\mathscr {F}: \mathcal {A}\rightarrow \mathbb {R} \), for fixed \(\xi \in \mathcal {A}, \mathscr {F}(\eta ) = \frac{d(\eta ,\xi )^2}{2 {\varDelta }\tau } + l(\eta )\) where l satisfies hypothesis (H2) and the following hypothesis (H1bis):
- (H1bis) :
-
\(\forall r>0, \forall \xi \in \mathcal A\), the set \(\{ l(z); z \in {\mathcal A}, d(z,\xi ) \le r \} \) is compact as a subset of \(\mathbb {R}\).
We will prove that the application \(\mathscr {F}\) attains its minimum in \(\mathcal {A}\).
Let \((z_k)_{k \in \mathbb {N}} \in \mathcal {A}\) be a sequence such that \(\lim _{k\rightarrow \infty } \mathscr {F} (z_k) = m \) with \( m = \inf _{\eta \in \mathcal {A}} \mathscr {F}(\eta )\). We show that \(d(z_k, \xi )\) is bounded then we find a minimizer.
Suppose \(d( z_k, \xi ) = \lambda > 1\) and let \(\gamma \) be the geodesic connecting \(\xi \) and \(z_k\) (\(\gamma (0) = \xi \) and \(\gamma (1) = z_k\)) satisfying (H2). Consider an element \(z = \gamma (t)\) on the geodesic \(\gamma \) such that \(d(z, \xi ) = 1\). By definition of a geodesic, \(d(z, \xi ) = d(\gamma (t), \gamma (0)) = |t-0|\cdot d(\gamma (0), \gamma (1)) = t d(z_k, \xi ) = \lambda \), implies necessarily \(t = 1/\lambda \) and using hypothesis (H2), we obtain:
For k large enough, \(z_k\) is such that \(\mathscr {F}(z_k) < 2m\) that is \(\frac{\lambda ^2}{2 {\varDelta }\tau } + l(z_k) < 2 m\). Rearranging Eq. (19) gives \(l(z_k) = \lambda l(z) - (\lambda -1) l(\xi )\), note that \(l(\xi )\) is constant and as \(d(z, \xi ) = 1\), hypothesis (H1bis) insures that l(z) is bounded, so \(l(z_k)\) is also bounded (by \(M \in \mathbb {R}\)). Consequently, there exists a constant C, depending on \(m, \tau \) and M such that \(\lambda < C\).
As \(d(z_k, \xi )\) is bounded, there exists a converging subsequence, \((x_k)_{k \in \mathbb {N}} \in \mathcal {A}\) such that \(d(x_k, \xi ) \rightarrow d_1\), and \(l(x_k)\) converges to \(m - \frac{d_1^2}{2 {\varDelta }\tau }\). Starting from some index \(k, x_k \in B(\xi , d_1)\) (the ball with center \(\xi \) and radius \(d_1\)) and if it is not the case, we replace \(x_k\) by the element at the intersection of the ball and the geodesic \(\gamma \) linking \(x_k\) and \(\xi \). This operation does not change the limit of \(l(z_k)\): noting temporarily \(\tilde{x}_k\) the obtained sequence, as l satisfies hypothesis (H2), we have \(l(\tilde{x}_k) = (1- \frac{d_1}{d(x_k, \xi )}) l(\xi ) + \frac{d_1}{d(x_k, \xi )}l(x_k)\) with \(d(x_k, \xi ) \rightarrow d_1\). By hypothesis (H1bis), the set \({l(x_k)}\) is closed, so there exists an element \(Z \in B(\xi , d_1)\) such that \(l(Z) = m - \frac{d_1^2}{2 {\varDelta }\tau }\) with, of course, \(d(Z, \xi ) \le d_1\). In this case, \(\mathscr {F}(Z) = \frac{d(Z, \xi )^2}{2 {\varDelta }\tau } + m - \frac{d_1^2}{2 {\varDelta }\tau } \le m\), so Z is a minimizer.
Let \(n \le 4, Y_1, \dots , Y_n \in \mathcal {A}, \lambda _1, \dots , \lambda _n \in \mathbb {R}\) and \(\mathcal {P}\) satisfying (H1) and (H2). We define the linear combination l of \(\mathcal {P}\) by: \(l(\eta ) = \sum _{i=1}^n \lambda _i \mathcal {P}(\eta , Y_i)\). It is immediate that l satisfies (H2); we prove below that if \(\mathcal {P}\) satisfies (H1) then l satisfies (H1bis). As all elements of the sum are in a compact set, the sum is bounded. Let \((z_k)_{k \in \mathbb {N}} \in \mathcal {A}\) be a sequence such that \(\lim _{k\rightarrow \infty } l(z_k) = \alpha \). Using (H1), there exists a subsequence of \(z_k\) (also noted \(z_k\)) such that the element \((\mathcal {P}(z_k, Y_1), \dots , \mathcal {P}(z_k, Y_n)) \in \mathbb {R}^n\) converges to \((\alpha _1, \dots \alpha _n)\) so \(l(z_k)\) converges to l(z) with \(l(z) = \sum _{i=1}^n \lambda _i \alpha _i = \alpha \). \(\square \)
Proposition 2.2
Let \(\mathcal {H}\) be a Hilbert space, \(f: \mathcal {H} \rightarrow \mathcal {H}\) be a function and \(\mathcal {P}(x, y) = \left\langle x,-f(y)\right\rangle \). Then the schemes VEE (Eq. 15), VH (Eq. 16), VRK3 (Eq. 17), and VRK4 (Eq. 18) correspond to Explicit Euler (Eq. 11), Heun (Eq. 12), Runge–Kutta 3 (Eq. 13) and Runge–Kutta 4 (Eq. 14) schemes respectively for the evolution Eq. (9).
Proof
We prove hypothesis (H1) in the case \(n=1\), the proof for the general case is left to the reader. We show that, for \(\xi \in \mathcal {H}, \mathcal {D} = \{ \left\langle z,-f(\xi )\right\rangle ; z \in {\mathcal H}, d(z,\xi ) \le r \} \) is bounded and closed. Denoting \(||\cdot ||\) the norm of \(\mathcal {H}\), by the triangular inequality when \(x \in B(\xi , r), ||x|| \le r + ||\xi ||\). For \(x \in \mathcal {D}, ||\left\langle x,f(\xi )\right\rangle || \le ||x||\cdot ||f(\xi )|| \le (r + ||\xi ||) \cdot ||f(\xi )|| < \infty \), so \(\mathcal {D}\) is bounded. Let z be an element of \(B(\xi ,r)\) and \(x \in \mathcal {H}\) such that \(z = \xi + x\) with \(||x|| \le r\), by linearity \(\left\langle z,-f(\xi )\right\rangle = \left\langle \xi ,-f(\xi )\right\rangle + \left\langle x,-f(\xi )\right\rangle \). Furthermore, \(\left\langle x,-f(\xi )\right\rangle \le ||x||\cdot ||f(\xi )|| \le r ||f(\xi )||\) so \(\{\left\langle z,-f(\xi )\right\rangle , z \in B(\xi , r) \}\) is included in the closed interval \([\left\langle \xi ,-f(\xi )\right\rangle - r ||f(\xi )||, \left\langle \xi ,-f(\xi )\right\rangle + r ||f(\xi )||]\). If we take \(x = - \lambda \frac{f(\xi )}{||f(\xi )||}\) with \(\lambda \in [-r,r]\) then \(\left\langle x,-f(\xi )\right\rangle = \lambda ||f(\xi )||\) which implies \(\left\langle z,-f(\xi )\right\rangle = \left\langle \xi ,-f(\xi )\right\rangle + \lambda ||f(\xi )||\). Therefore \(\{\left\langle z,f(\xi )\right\rangle , z \in B(\xi , r) \}\) is exactly \([\left\langle \xi ,f(\xi )\right\rangle - r ||f(\xi )||, \left\langle \xi ,f(\xi )\right\rangle + r ||f(\xi )||]\) and is closed as a closed interval of \(\mathbb {R}\).
To prove hypothesis (H2), consider \(x_0, x_1 \in \mathcal {H}\) and for \(t \in [0,1], \gamma (t) = t x_1 + (1-t) x_0\) the segment between \(x_0\) and \(x_1\). For \(t, s \in [0,1]\), with straightforward computation, we have \(d(\gamma (t), \gamma (s)) = |t-s| \cdot d(\gamma (0), \gamma (1))\), showing that \(\gamma \) is a (constant speed) geodesic. In this case, we can also prove the uniqueness of the geodesic. In fact, assume that there exists another geodesic \(\gamma _2\) linking \(x_0\) and \(x_1\). Let \(z = \gamma _2(t)\) be an element only on \(\gamma _2\); using the triangle inequality and that \(\gamma \) is the segment linking \(x_0\) and \(x_1\) we have \(d(x_0, z) + d(z,x_1) > d(x_0, x_1)\). Furthermore as \(\gamma _2\) is a geodesic, left member of previous equation is \(t d(x_0, x_1) + (1-t) d(x_0, x_1) = d(x_0, x_1)\). We obtain a contradiction thus there does not exist another geodesic. The function \(t \rightarrow \mathcal {P}(\gamma (t), Y)\) has the form \(t \rightarrow \left\langle t x_1 + (1-t) x_0,-f(Y)\right\rangle \), so is linear.
We now prove the equivalence of schemes. The scheme VEE is defined by the following equation:
The application \(\mathscr {F}: \mathcal {H} \rightarrow \mathbb {R}, \mathscr {F} ( \eta ) = \frac{ ||\eta - \xi _k||^2}{2 {\varDelta }\tau } + \left\langle \eta ,f(\xi _k)\right\rangle \) is differentiable and its derivative \(\mathscr {F}': \mathcal {H} \rightarrow \mathcal {H} \) is \(\mathscr {F}' ( \eta ) = \frac{ \eta - \xi _k}{{\varDelta }\tau } + f(\xi _k)\). The minimizer of \(\mathscr {F}\) is a critical point, the minimum noted \(\xi _{k+1}\), satisfies \(\mathscr {F}'(\xi _{k+1})\) = 0. After trivial rearrangement, we obtain \(\xi _{k+1} = \xi _k - {\varDelta }\tau f(\xi _k)\), that is the Explicit Euler scheme presented in (11) for Eq. (9). Arguments are similar for other schemes. \(\square \)
3 Epidemiological model
From now on, all variables can have a double dependence: on the real time t which will be subscript and on the iteration time \(\tau \) (previously k in Eq. (7)) which will be a superscript l.
This section reproduces a short presentation of the epidemiological model and of the cost computing as introduced in [44, Sect. 2]. Please refer to the original article for more details and proof.
3.1 Spread of the disease
The model simulates the dynamics of an epidemic in a population. The final time period T is supposed finite, and the time horizon can be discretized in \((N+1)\) time instants, noted \(t_0 = 0, t_1 = {\varDelta }T, \ldots , t_N = T\).
The model is compartmental, and the population is divided into several classes: susceptible individuals \((S_n)\) is the proportion of individuals susceptible to catch the disease at time instant \(t_n\); infected individuals regroups individuals who are infected, more precisely, \(I_n^{\omega }\) is the proportion of individuals who have been infected at time instant \(t_{n-\omega }\) (with \(\omega \in \{0, 1, \ldots , {\varOmega }\}\)) and \(I_n\) is the sum of all \(I_n^{\omega }\); recovered individuals \((R_n)\) is the proportion of individuals who have exited from one infected class; vaccinated individuals (\(V_n^{\theta }\) where \(\theta \) counts the time between vaccination time instant and current time instant with \(\theta \in \{ 0, 1, \ldots , {\varTheta }\}\)) is the proportion of vaccinated individuals at time instant \(t_n\) and not infected since; failed vaccinated individuals (\(F_n\)) is the proportion of individuals who were vaccinated at some time \(t \le t_n\) but whose vaccination failed and have not been infected since. As individuals can only take vaccination once, \(V^{{\varOmega }}\) represents individuals who have lost immunity. Lastly \(U_n\) represents the proportion of individuals vaccinated between \(t_n\) and \(t_{n+1}\), a proportion f of them will never develop any immunity and go to the failed vaccination class.
The specificities of disease are represented by the vector \(\gamma = (\gamma ^0, \ldots , \gamma ^{{\varOmega }}) \in (\mathbb {R}_+)^{{\varOmega }+1}\) (which reflects the recovery of an infected individual); the function \(\beta (t)\) (which characterizes the contact between an infected individual and a susceptible one at time t); and the vector \(\alpha \in [0,1]\) which quantifies the protection provided by the vaccine in terms of the probability of infection.
The model is defined by the following system of equations:
with initial conditions:
The continuous version of the previous discrete model is given by the following equations:
with initial and boundary conditions:
In this continuous time model, we use the following notation (refer to [44, Appendix] for details):
-
\(U : [0, \infty ] \rightarrow [0,1]\) represents the vaccination, U(t) is the fraction of the population that has been vaccinated by time t;
-
constants \(\beta , \gamma , f\) and function A, have the same meaning as in the discrete model;
-
\(V(t, \theta )\) is the fraction of individuals who have been vaccinated at time \(t-\theta \) and have not been infected since vaccination.
3.2 Individual vaccination: as minimization of individual cost
The system (20) presents the spread of the disease at the population level. For an individual, we consider the Markov chain \(M_n\) which describes, in terms of transition probabilities (see [44, Eq 2.10]), the state of the individual (\(M_n \in \{ S, R, F, I^j, V^l \}\) with \(j\in \{0, \dots , {\varOmega }\}\) and \(l \in \{ 0, \dots , {\varTheta }\}\)).
In this subsection, we look for the global decision of vaccination \(U_n\) which corresponds to the sum of individual decisions, as shown by Eq. (24). In order to reflect the decision of an individual, we introduce the probability law \(\mu \) defined on \(\{ t_0, \ldots , t_{N-1} \} \cup \infty \). In practice, before the epidemic starts, the individual chooses his probability of vaccination at each time step \(t_k\), assuming his non-infected status.
The collection of conditional rates \(\lambda = (\lambda _n)_{n=1}^N\) is given by the law \(\mu \):
The individual in the susceptible class at time \(t_0\) (\(M_0 = S\)) has the following cost: \(J_{indi}(\mu ; U) = \left\langle \mu ,\mathscr {C}_U\right\rangle \) where \(\mathscr {C}_U\) (the vector representing the cost of all pure strategies “vaccination sure at time t”) is defined by:
Here \(\varphi ^{V,I}_n = 1- \prod _{k=n}^{{\varTheta }} \left( 1- \beta _{{\varDelta }T}^k\alpha _{k-n-1} I_k \right) ,\) with the convention \(\alpha _{-1}=1\) and \(\varphi ^I_n= 1-\prod _{k=0}^{n} \left( 1- \beta _{{\varDelta }T}^kI_k \right) (\text {for } n=0,\ldots ,N-1), \varphi ^I_\infty = 1-\prod _{k=0}^{N} \left( 1- \beta _{{\varDelta }T}^kI_k \right) \).
The equilibrium between individual dynamics and global dynamics (20) is attained when:
The purpose of an individual is to minimize \(J_{indi}(\mu ; U) \) under the constraint \(\mu \in {\varSigma }_{N+1}\) with \({\varSigma }_{N+1} = \left\{ x \in \mathbb {R}^{N+1} \left| \sum _{k=0}^N \right. x_k =1 \text { and } x_k \ge 0, 0 \le k \le N \right\} \).
The numerical schemes give a sequence of probability laws \(\mu ^l\); here the index “l” means that \(\mu ^l\) is a numerical approximation of the solution \(\xi \) of (8) at the time \(\tau _l = l\cdot {\varDelta }\tau \). For instance, if we consider the discrete model (20), each \(\mu ^l = (\mu ^l_0, \mu ^l_1, \dots , \mu ^l_{N-1}, \mu ^l_\infty ) \in {\varSigma }_{N+1}\) is a probability law on \(\{t_0; \dots , t_{N-1}\} \cup \{\infty \}\), while if we consider the continuous time model (21), each \(\mu ^l\) is a probability law on \( [0,T] \cup \{\infty \}\). Note that here a step of the numerical schemes increments the “l” index in order to converge to a Nash–MFG equilibrium.
3.3 Numerical values
This section describes the numerical values used in the two cases presented in this paper: Short persistence, large efficacy and Long persistence, 100% efficacy. These values are sensibly similar to the ones used in [44, Subsect. 4.2 and 4.3].
For the epidemic parameters, we consider a total simulation time at 1 year (\(T = 1\)); an initial proportion of susceptibles \(S_0 = 0.94\) and infected \(I_0 = 2.0 ^{-6}\); three time instants by day (\(N = 365*3\)); a recovery rate \(\gamma ^{\omega } = 365/3.2\) (\({\varOmega }= 20\)); the reproduction number \(R_0 = 1.35\) thus \(\beta = \gamma R_0\); \(t_2^{\beta } = 1/2\) such that \(\beta (t) = \beta \) for \(t \le t_2^{\beta }\) and then \(\beta (t) = \beta _{{ min}}\) for \(t > t_2^{\beta }\) where \(\beta _{\min } = \gamma /S_0\); the relative cost of the epidemic is \(r_I= 1\).
Finally, we introduce \(t_1^{\alpha }, t_2^{\alpha }\) and set \(\alpha _{\theta } = 1 - \mathbf {1}_{[t_1^{\alpha }, t_2^{\alpha }]}\). For the case Short persistence, large efficacy, \(t_1^{\alpha } = 5/365\) and \(t_2^{\alpha }=1/12\) (\({\varOmega }= 93\)). For the case Long persistence, 100% efficacy, \(t_1^{\alpha } = 0\) and \(t_2^{\alpha }=1/2\) (\({\varOmega }= 549\)). We suppose a failure rate \(f=0\), and the relative cost of the vaccination \(r_V= 0.005\).
4 Numerical illustration
4.1 Framework
Salvarani and Turinici introduced an epidemiological model with possibility of vaccine which has imperfect efficiency and limited persistence recalled in the previous Sect. 3. To model the spread of the disease in the population, they use a compartmental model and a probability distribution to reflect the individual decisions.
To find stable individual decision they define the problem as a Nash equilibrium. They find the probability distribution \(\mu \) making an individual indifferent to change his vaccination decision if all individuals have the same \(\mu \) (see [44, Theorem 2.1] for the proof of equilibrium existence).
Let \(\mathscr {C}_{\mu }\) be the cost of pure strategies “vaccination happens at time t” when all individuals choose as vaccination strategy \(\mu \). In that case \(\left\langle \eta ,\mathscr {C}_{\mu }\right\rangle \) represents the cost of an individual with strategy \(\eta \) when others use strategy \(\mu \).
The definition of the mapping \(E(\mu )\) as introduced in [44] is the maximum gain obtained by an individual if he changes unilaterally his strategy and everybody else remains with the strategy \(\mu \); with mathematical notation:
A minimum of the mapping \(\mu \rightarrow E(\mu )\) is a Nash equilibrium. To find it they introduce an iterative method depending on a step \({\varDelta }\tau \) (algorithm 1). In this algorithm, Eq. (25) is used. The following intuitive interpretation is also provided: an individual in a population with strategy \(\mu ^l\) will, if necessary, adjust his strategy to minimize his cost function \(\eta \rightarrow \left\langle \eta ,\mathscr {C}_{\mu ^l}\right\rangle \) while at the same time keeping small the distance between the previous strategy \(\mu ^l\) and the new \(\mu ^{l+1}\). But Algorithm 1 can also be seen as numerical resolution of an evolution e.g., describing a curve in the metric space of the admissible strategies. That is why we apply variational methods in the metric space \({\varSigma }_{N+1}\) with the standard Euclidean distance to obtain faster convergence to the Nash equilibrium. Here \(\mathcal {P}(x, y) = \left\langle x,\mathscr {C}_y\right\rangle \) and for the same reasons as in the proof of Proposition 2.2, \(\mathcal {P}\) satisfies hypotheses (H1) and (H2). Note that with high order schemes VH, VRK3 and VRK4, we lose the intuitive idea but we increase the Nash equilibrium computing.
4.2 Results
In this section we test the variational schemes proposed in Sect. 2.2.
The reader is invited to refer to [44, Sect. 4.6] or Sect. 3 for the numerical values of the two studied cases Short persistence, large efficacy (with vaccine persistence at 1 month with a delay of action at 5 days) and Long persistence, 100% efficacy (with persistence of the vaccine at 6 months with no delay of action).
In order to appreciate the convergence scheme for the case Short persistence, large efficacy the evolution of the mapping \(E(\cdot )\) is presented (see Fig. 6 for VEE with \({\varDelta }\tau = 0.1\)), where the evolution for the four schemes with lower \({\varDelta }\tau \) value is added to stabilize result of mapping \(E(\cdot )\).
For the Short persistence, large efficacy case, Fig. 7 presents the numerical estimation of the scheme error and order: only the VEE is of order one, the others are of order two in \({\varDelta }\tau \). Furthermore, Fig. 8 introduces the obtained errors for each schemes depending of the number of minimizations. Recall that a high order scheme needs a very regular functional to provide high order convergence, but in our case the regularity of the function is completely unknown. However, this case remains very interesting because it shows that even if the regularity of the function is not enough to obtain order three or four, the VRK3 and VRK4 still improve the regularity of the function and have a better order than the VEE scheme.
For the other example Long persistence, 100% efficacy, the convergence is faster: Fig. 9 compares the evolution of mapping \(E(\cdot )\) for the two cases. In the second case, the graph in Fig. 10 shows that the scheme has high order as it was expected previously. VH is indeed of order 2 and VRK3 and VRK4 are respectively order 3 and 4. Furthermore, Fig. 11 introduces the error for each scheme depending on the number of minimisations: this shows the accuracy of the schemes.
5 Second application: an economical approach
5.1 Epidemiological and economic model
This numerical application is based on [22]. In order to model the spread of the disease, three classes (as explained in Sect. 1.1) are introduced: Susceptible (S), Infected (I), Recovered (R) with the following time evolution:
where S(t), I(t), R(t) represent respectively the number of susceptible, infected and recovered individuals at time t, N the total number of individuals and r(t) the rate of vaccination at time t. We also introduce the utility \(u_H\) of a healthy individual and the utility \(u_S\) when sick.
The individual cost introduced is represented by the expected loss if individual is falling ill \(\left( EL = \frac{u_H-u_S}{\gamma } \left( 1 - \frac{S(\infty )}{S(0)} \right) \right) \) and the cost of vaccination \(\theta \).
The computation that follows, taken from [22], is valid when the dynamics (26) traverses at most two regions of vaccination/non vaccination. Under this assumption, the expected payoff of the two pure strategies (vaccination or not) is equal on the frontier of the vaccination/non vaccination region, that is: \(\theta = \frac{u_H-u_S}{\gamma } \left( 1 - \frac{S(\infty )}{S(0)} \right) \).
Set \(w = \theta \gamma /(u_H-u_S)\). Rewrite previous equation under the following form: \(w = 1 - S(\infty )/S(0)\), that gives \(S(\infty ) = S(0) (1-w)\). Using that, with \(\rho = N\gamma /\beta , \frac{d}{dt} S(t)e^{R(t)\rho } = 0\) for \(r(t) = 0 \ \forall t \ge 0, S(0)e^{R(0)\rho } = S(\infty )e^{R(\infty )\rho }\) (see [1]), and after replacing of \(S(\infty )\) by \(S(0) (1-w)\), we obtain: \(R(0) = \rho \ln (1-w) + R(\infty )\). As the total number of the population is constant equal to N, we obtain Eq. (27). This last equation can be seen as a “switching curve” in the plane (S, I):
where the individual is indifferent between being vaccinating and doing nothing.
With this economical approach, there is no insurance that the curve \({\varGamma }\) generates a Nash–MFG equilibrium. In fact, [22] does not use Nash–MFG framework but assumes equality between two payoffs, false in general cases (see for example [24]).
Nevertheless we seek to replicate these results in order to benchmark the performance of our numerical schemes. For each initial point, the algorithm reports if vaccination is optimal or not, which can be compared to the region given by the switching curve \({\varGamma }\).
We use the same approach as in Sect. 3.2 to determine the global vaccination response r(t). The proportion of individuals vaccinated up to time t is the cumulative function V(t) defined by:
where r(u) is the instantaneous vaccination rate at time u with \(r(u) \in [0, r_{{ max}}]\) and \(r_{{ max}} < \infty \).
5.2 Numerical values
For the simulation, we use the following numerical values: an initial proportion of susceptibles \(S_0 = 0.75\) and infected \(I_0 = 0.1\); a cost of the vaccine \(\theta = 0.5\); an expected utility loss associated if illness \(\frac{\gamma }{u_H - u_S} = 1\); a recovery rate for an infected individual \(\gamma = 365/10\) and a transmission rate of the disease \(\beta = 73\). As in the Sect. 4, total simulation time is 1 year (\(T=1\)) with \(L = 365 \times 3\) total time steps.
5.3 Results
We obtain coherent results with the analytic solution by using the same approach as in Sect. 4.1. The equilibrium individual strategy is a mixed strategy with \(\xi _0 = 66 \%\) and \(\xi _{\infty } = 33 \%\) and its cost is 0.51. This results are coherent with the analytic result that provides a mixed strategy with \(\xi _0 = 34\%, \xi _\infty = 66\%\) and a cost at 0.5. Figure 12 presents the numerical estimation errors and order for the four schemes and Fig. 13 the errors for each schemes depending of the number of minimizations.
6 Perspectives
This work introduces three high order schemes to find a Nash equilibrium and illustrates their use in an epidemiological application.
We show numerically that the schemes VH, VRK3 and VRK4 exhibit better order than VEE and we can obtain, depending on the regularity of the function, up to order four convergence. The approach based on the Runge–Kutta method can be applied for other numerical schemes, for instance Mid Point, Leap–Frog or Adams–Bashford methods. The presented algorithms are not optimized, as minimizations can take a long time and need a high numerical precision, so the time execution can be significant. A perspective of this work is to provide an extension to bi-dimensional or tri-dimensional problems.
References
Abakuks, A.: Optimal immunisation policies for epidemics. Adv. Appl. Probab. 6, 494–511 (1974)
Ambrosio, L., Gigli, N., Savare, G.: Gradient Flows in Metric Spaces and in the Space of Probability Measures, 2nd edn. Lectures in Mathematics. Birkhäuser, Basel: ETH Zürich (2008)
Bauch, C.T.: Imitation dynamics predict vaccinating behaviour. Proc. R. Soc. Lond. B Biol. Sci. 272(1573), 1669–1675 (2005)
Bauch, C.T., Bhattacharyya, S.: Evolutionary game theory and social learning can determine how vaccine scares unfold. PLOS Comput. Biol. 8(4), 1–12 (2012)
Bauch, C.T., Galvani, A.P., Earn, D.J.D.: Group interest versus self-interest in smallpox vaccination policy. Proc. Natl. Acad. Sci. 100(18), 10564–10567 (2003)
Behncke, H.: Optimal control of deterministic epidemics. Optim. Control Appl. Methods 21(6), 269–285 (2000)
Blanchet, A., Carlier, G.: Optimal transport and Cournot–Nash equilibria. ArXiv e-prints (June 2012)
Blanchet, A., Carlier, G.: From Nash to Cournot–Nash equilibria via the Monge–Kantorovich problem. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 372(2028), 20130398 (2014)
Blanchet, A., Carlier, G.: Remarks on existence and uniqueness of Cournot-Nash equilibria in the non-potential case. ArXiv e-prints (May 2014)
Brauer, F., Van Den Driessche, P., Wu, J.: Mathematical epidemiology. No. 1945 in Lecture Notes in Mathematics/Mathematical Biosciences Subseries. Berlin: Springer (2008)
Buonomo, B., d’Onofrio, A., Lacitignola, D.: Global stability of an SIR epidemic model with information dependent vaccination. Math. Biosci. 216(1), 9–16 (2008)
Cardaliaguet, P., Hadikhanloo, S.: Learning in Mean Field Games: The Fictitious Play. ArXiv e-prints (2015)
Cheeger, J.: Differentiability of Lipschitz functions on metric measure spaces. Geom. Funct. Anal. 9(3), 428–517 (1999)
Colizza, V., Barrat, A., Barthelemy, M., Valleron, A.J., Vespignani, A.: Modeling the worldwide spread of pandemic influenza: baseline case and containment interventions. PLoS Med. 4(1), 1–16 (2007)
Diekmann, O., Heesterbeek, J.: Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation. Wiley Series in Mathematical & Computational Biology. Wiley, London (2000)
d’Onofrio, A., Manfredi, P., Poletti, P.: The impact of vaccine side effects on the natural history of immunization programmes: an imitation-game approach. J. Theor. Biol. 273(1), 63–71 (2011)
d’Onofrio, A., Manfredi, P., Poletti, P.: The interplay of public intervention and private choices in determining the outcome of vaccination programmes. PLOS ONE 7(10), 1–10 (2012)
d’Onofrio, A., Manfredi, P., Salinelli, E.: Vaccinating behaviour, information, and the dynamics of SIR vaccine preventable diseases. Theor. Popul. Biol. 71(3), 301–317 (2007)
d’Onofrio, A., Manfredi, P., Salinelli, S.: Bifurcation thresholds in an sir model with information-dependent vaccination. Math. Model. Nat. Phenom. 2(1), 2643 (2007)
Doutor, P., Rodrigues, P., Soares, M.C., Chalub, F.A.C.C.: Optimal vaccination strategies and rational behaviour in seasonal epidemics. J. Math. Biol. 73(6), 1437–1465 (2016)
Francis, P.J.: Dynamic epidemiology and the market for vaccinations. J. Public Econ. 63(3), 383–406 (1997)
Francis, P.J.: Optimal tax/subsidy combinations for the flu season. J. Econ. Dyn. Control 28(10), 2037–2054 (2004)
Fu, F., Rosenbloom, D.I., Wang, L., Nowak, M.A.: Imitation dynamics of vaccination behaviour on social networks. Proc. R. Soc. Lond. B Biol. Sci. 278(1702), 42–49 (2010)
Galvani, A.P., Reluga, T.C., Chapman, G.B.: Long-standing influenza vaccination policy is in accord with individual self-interest but not with the utilitarian optimum. Proc. Natl. Acad. Sci. 104(13), 5692–5697 (2007)
Hethcote, H.W.: The mathematics of infectious diseases. SIAM Rev. 42(4), 599–653 (2000)
Huang, M., Caines, P.E., Malhamé, R.P.: An invariance principle in large population stochastic dynamic games. J. Syst. Sci. Complex 20(2), 162–172 (2007)
Huang, M., Caines, P.E., Malhamé, R.P.: Large-population cost-coupled LQG problems with nonuniform agents: individual-mass behavior and decentralized \(\epsilon \)-Nash equilibria. IEEE Trans. Autom. Control 52(9), 1560–1571 (2007)
Huang, M., Malham, R.P., Caines, P.E.: Large population stochastic dynamic games: closed-loop McKean–Vlasov systems and the Nash certainty equivalence principle. Commun. Inf. Syst. 6(3), 221–252 (2006)
Hubert, E., Turinici, G.: Nash-MFG equilibrium in a SIR model with time dependent newborn vaccination. https://hal.archives-ouvertes.fr/hal-01389584 (2016)
Ibuka, Y., Li, M., Vietri, J., Chapman, G.B., Galvani, A.P.: Free-riding behavior in vaccination decisions: an experimental study. PLOS ONE 9(1), 1–9 (2014)
Jordan, R., Kinderlehrer, D., Otto, F.: The variational formulation of the Fokker–Planck equation. SIAM J. Math. Anal. 29(1), 1–17 (1998)
Junyuan, Y., Maia, M., Yuming, C.: Imitation dynamics of vaccine decision-making behaviours based on the game theory. J. Biol. Dyn. 10(1), 31–58 (2016). PMID: 26536171
Laguzet, L., Turinici, G.: Individual vaccination as Nash equilibrium in a SIR model with application to the 2009–2010 influenza A (H1N1) epidemic in France. Bull. Math. Biol. 77(10), 1955–1984 (2015)
Lasry, J.M., Lions, P.L.: Jeux à champ moyen. I: Le cas stationnaire. C. R. Math. Acad. Sci. Paris 343(9), 619–625 (2006)
Lasry, J.M., Lions, P.L.: Jeux à champ moyen. II: Horizon fini et contrôle optimal. C. R. Math. Acad. Sci. Paris 343(10), 679–684 (2006)
Lasry, J.M., Lions, P.L.: Mean field games. Jpn. J. Math. 2(1), 229–260 (2007)
Legendre, G., Turinici, G.: Second-order in time schemes for gradient flows in Wasserstein and geodesic metric spaces. C. R. l’Acad. Sci. Sér. I Math. 355(3), 345–353 (2017)
Manfredi, P., d’Onofrio, A. (eds.): Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases. Springer, New York (2013)
Morton, R., Wickwire, K.H.: On the optimal control of a deterministic epidemic. Adv. Appl. Probab. 6, 622–635 (1974)
Nardin, L.G., Miller, C.R., Ridenhour, B.J., Krone, S.M., Joyce, P., Baumgaertner, B.O.: Planning horizon affects prophylactic decision-making and epidemic dynamics. PeerJ 4, e2678 (2016)
Oraby, T., Thampi, V., Bauch, C.T.: The influence of social norms on the dynamics of vaccinating behaviour for paediatric infectious diseases. Proc. R. Soc. Lond. B Biol. Sci. 281, 1780 (2014)
Piunovskiy, A.B., Clancy, D.: An explicit optimal intervention policy for a deterministic epidemic model. Optim. Control Appl. Methods 29(6), 413–428 (2008)
Reluga, T.C., Galvani, A.P.: A general approach for population games with application to vaccination. Math. Biosci. 230(2), 67–78 (2011)
Salvarani, F., Turinici, G.: Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence. Working paper or preprint (2016)
Shiryaev, A.N.: Probability-1, 3 edn., vol. 95 of Graduate Texts in Mathematics. New York: Springer (2016)
Turinici, G.: Metric gradient flows with state dependent functionals: the Nash-MFG equilibrium flows and their numerical schemes. https://hal.archives-ouvertes.fr/hal-01528480 (2017)
Verelst, F., Willem, L., Beutels, P.: Behavioural change models for infectious disease transmission: a systematic review (2010–2015). J. R. Soc. Interface 13, 125 (2016)
Acknowledgements
This work was supported by the French Research Agency, program Investments for the future, project ANR-10-BINF-07 (MIHMES) and by the European Union, through the European fund for the regional development of Pays-de-la-Loire (FEDER).
Author information
Authors and Affiliations
Corresponding author
Additional information
This article belongs to the Special Issue: Demographic and temporal heterogeneity in infectious disease epidemiology.
Rights and permissions
About this article
Cite this article
Laguzet, L. High order variational numerical schemes with application to Nash–MFG vaccination games. Ricerche mat 67, 247–269 (2018). https://doi.org/10.1007/s11587-018-0366-z
Received:
Revised:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11587-018-0366-z