Finite State Graphon Games with Applications to Epidemics

We consider a game for a continuum of non-identical players evolving on a finite state space. Their heterogeneous interactions are represented with a graphon, which can be viewed as the limit of a dense random graph. A player’s transition rates between the states depend on their control and the strength of interaction with the other players. We develop a rigorous mathematical framework for the game and analyze Nash equilibria. We provide a sufficient condition for a Nash equilibrium and prove existence of solutions to a continuum of fully coupled forward-backward ordinary differential equations characterizing Nash equilibria. Moreover, we propose a numerical approach based on machine learning methods and we present experimental results on different applications to compartmental models in epidemiology.

spread mechanisms, the network of social interactions, and society-wide efforts to stop or slow the spread. As individuals, we have all made choices during the ongoing pandemic about the extent to which we minimize our personal risk of being infected. However, there is a trade-off between being careful and the pursuit of happiness. As we all have learned by now, our risk is not only determined by our own vigilance but also by others' choices and our environment.
In the framework of rational agents 1 , each individual anticipates the action of their neighbours, neighbours' neighbours, etc., and any other external influence, then selects their action as a best response to others' actions. In other words, they want optimize their private outcome while taking into account their surrounding environment, which includes other agents' actions. This type of strategic interaction is a non-cooperative game. The communication between agents in the game may be restricted by geography, social circles, and other factors. Moreover, people interact with different intensity, depending on their occupation and personality. Hence, the agents in the game, each with their own predisposition for risk, will act in a wide variety of ways and thus naturally form a heterogeneous crowd.
Consider the game discussed above where all agents anticipate the others' action and then selfishly plays a best response. Strategy profiles (the collection of all players' actions) consistent with such behavior are Nash equilibria, i.e., profiles such that no player profits from a unilateral deviation. Computing Nash equilibria in games with large number of players is particularly hard; under some specific assumptions, approximate equilibria can be found by using the mean field game approach developed independently by Lasry-Lions [45,46] and Huang-Malhamé-Caines [35]. The approach has found many practical applications, examples include the mathematical modeling of price movement in energy markets; pedestrian crowd motion and evacuation; and epidemic disease spread.
One of the fundamental assumptions in mean field game theory is that agents are indistinguishable and interact homogeneously. However, in some real-world applications such as the modeling of epidemics, the diversity of individuals and the variation of their interactions are important factors to consider. Examples of such include the effects of travel restrictions, multiple age groups with distinct social behavior and risk profiles, and the spectrum of preexisting health conditions. The aspects listed above require the game to have a non-uniform underlying network structure. Games with a large number of non-identical players can be analyzed with so-called graphon games whenever the network specifying the interactions is dense 2 .
Epidemics are driven by the spread of the disease from infected to susceptible agents. The set of susceptible agents is not necessarily the whole non-infected population, for certain diseases immunity is gained after exposure. The evidence suggests that the COVID-19 virus mainly spreads through close contact 3 . Fortunately, the disease transmission probability can be decreased by efforts of the individual. For example, an individual can choose to avoid public and closed spaces, wear a protective mask or do their shopping online. When two people meet the disease transmission probability depends on both sides' effort. The disease is less likely to occur if both parts wear protective masks than if just one part does. However, the decrease in risk of transmission is not additive in the interacting agents' efforts. Using this intuition, we assume in this paper that disease transmission likelihood depends on efforts in a multiplicative way. The effort of the agents will therefore be given the name "contact factor control", in line with the game-based model introduced in [2].
Before giving the full description of the epidemiological graphon game model, how it can be approached numerically, and all technical details, we expand upon the heuristics of graphon-based interaction with contact factor control and review the literature of related fields of research in the next sections.

The SIR Graphon Game with Contact Factor Control
The most famous compartmental model in epidemics is arguably the classical Susceptible-Infected-Removed (SIR) model 4 . In this section, we take advantage of its wide familiarity and compact formulation to further motivate for the concept of contact factor control and graphon-type interaction.
In order to give a description of the rate at which agents become infected, we need to first introduce some notation. Consider N individuals who transition between the states Susceptible (S), Infected (I), and Removed (R). An individual in state R has either gained immunity or deceased. Denote the state of agent j ∈ {1, . . . , N } at time t by X j,N t . A susceptible individual might encounter infected individuals, resulting in disease transmission. Encounters occur pairwise and randomly throughout the population with intensity β. The number of encounters with infected agents in a short time interval [t −Δt, t) is approximately proportional to the the share of the population in state I at t. Between each agent pair ( j, k), we set the interaction strength to w(x j , x k ), where w is a graphon (see Definition 1 for its definition) and x j , x k are random variables uniformly distributed on [0, 1]. Hence agent k's transition rate from state S to I is scaled by w and of the form β 1 N N k=1 w(x j , x k )1 I (X k,N t− ). Upon infection, an individual starts the path to recovery. The jump from state I to R happens after an exponentially distributed time with rate γ . The state R is absorbing.
Denoting Z j,N t := 1 N N k=1 w(x j , x k )1 I (X k,N t− ), the transition rate matrix for player j ∈ {1, . . . , N } is Here we use the order S, I, R for the columns and the rows. For instance, the term β Z j,N t encodes the rate at which an agent arrives in I coming from S. As explained above, this rate is proportional to the weighted average Z j,N t of infected agents interacting with player j. At this stage, the distinguishability of the players is seen in the aggregate variables (Z j,N t ) N j=1 which in general differ in value. As N → ∞, we expect the probability distribution flow of player j to converge to the solution of the ordinary differential equation (ODE) where p x j 0 is some given initial distribution over the states and Z x j Here p x j (t) is understood as a vector of length 3 whose coordinates correspond to the probability for player x j of being in state S, I and R, respectively at time t. Equation (2) encodes the evolution of these probabilities. Z x j t is the average number of infected agents around player x j , weighted by the pairwise interaction strength w(x j , y). Scaling p x (t), x ∈ [0, 1], by a population size N , N p x (t) =: (S x (t), I x (t), R x (t)), we retrieve a formulation of the compartmental SIR model with graphon-based interactionṡ So far we have considered a model without action. Assuming agents choose a "contact factor" in order to decrease their risk of getting infected, which enters the model in line with the discussion above, we express this new feature mathematically as follows: Given that the meeting frequency is β, pairing is random, disease spreads from infected agents to susceptible, and the spread probability is scaled by the efforts of the individuals that meet in a multiplicative way, the transition rate for individual j from S to the I is where α k t denotes the (contact factor) action of individual k ∈ {1, . . . , N } at time t, selected from set of actions A. Along the lines of the heuristics of mean field game theory, we anticipate that in an appropriate approximation of our interacting system in the limit N → ∞, agent x transitions from susceptible to infected with rate where ρ y t is the joint distribution of action and state of player y at time t and α x t is the action of individual x.

Graphon Games and Finite State Space Mean Field Games
Challenges related to large populations arise in the game theory, and mean field game (MFG) theory presents a toolbox to compute equilibria in these games with large number of players. MFGs were first developed for continuous state space models [35,45,46] and later for a finite state space models [32,33,42]. The theory for finite state MFGs has been extended in many directions, with contributions including a probabilistic approach [14], the master equation approach [4], minor-major player games [12], and extended games 5 [13]. Further, finite state mean field control, risk-sensitive control, and zero-sum games are treated in [17][18][19] which cover cases of unbounded jump intensities. Graphon games [3,[7][8][9]23,30,49] have recently been receiving an increasing research interest. The motivation is the study of strategic decision making in the face of a large non-complete dense networks of distinguishable agents. The graphon game's rising popularity stems from its ability to handle heterogeneity of agents to an extent far beyond the MFG theory.

Mean Field Games and Related Models for Epidemics
Decision making in compartmental models (e.g., the SIR model) have been studied intensively for a long time, with an increasing interest recently with the COVID-19 pandemic. In the form of games and optimal control problems, disease-combating efforts ranging from strategies for social contracts to vaccination have been analyzed in the literature. Here, we focus on work relying on the graphon game theory and the mean-field approach.
During an epidemic where the disease is prone to close-contact transmission (one example being COVID-19) control of the contact rate or other social distancing protocols are go-to solutions in the fight against the disease. Such strategies and related variations have been studied in the context of mean field games [16,28], mean field optimal control [47], and mean field type games [55]. In reality the population needs to be tested for the disease in order to accurately asses the risks in the decision making process. Two recent papers studying optimal testing policies are [15,36]. The effects of vaccination and accompanying decision making problems are not studied in this paper; they have been analyzed in MFG-related settings since before the COVID-19 pandemic [27,31,37,43,44,50]. Recently, the interplay between a population in which a disease spreads and a regulator has been studied in the form of Stackelberg games. In such models, the members in the populations are taking actions based on policies issued by the regulator, while the regulator anticipates the population's reaction and optimizes the policy. The case of a cooperative population has been studied in [36], while in [2] a population of selfish agents with contact factor control has been studied.
An example of a deterministic optimal control problem for centralized decision making during a pandemic in a society with multiple communicating subpopulations is given in [29]. The subpopulations interact over a non-uniform graph. A central planner wants to flatten the (global) curve of infections, leading to the optimal control problem. Sending the number of subpopulations to infinity, we anticipate a limit where each interaction is weighted by a graphon and the limit model would be reminiscent of the interacting system of Kolmogorov equations studied in this paper.

Networks and Graphons in Epidemiology
There is a vast body of literature on epidemiology modeling with network interactions. A review of the studies that use idealized networks 6 in epidemiology models can be found in [39]. More closely related to the ideas in this paper, there are recent contributions connecting epidemic models and graphons. In [57] a sensitivity analysis on the graphon SIS epidemic model is conducted. An infinite dimensional SIS model with application to targeted vaccination is considered in [24]. The paper [40] proposes a model with local-density dependent Markov processes interacting through a graphon structure, and considers applications to epidemiology. In a similar but more general setting to the SIR model with graphon interaction, [1] studies convergence of a stochastic particle system to an SIR-like system of PDEs with spatial interaction. We note [1] and its continuation [25,26] may be relevant for a future study of the convergence of N -player Nash equilibria to the equilibria in the finite state graphon game. The works mentioned in this section only consider the dynamics of the population without taking the agents' decision making into account.

Contributions and Paper Structure
This paper is, to the best of our knowledge, the first to address the analysis and numerical resolution of graphon games that are time-dependent and with a discrete state space. The application to epidemiology model departs from the traditional literature on epidemiology models and graphon models by the incorporation of a game theoretical aspect: here we go beyond dynamic graphon systems and find Nash equilibria for rational agents. We construct a probabilistic particle model for a continuum of interacting agents and prove that graphon aggregates must be deterministic (as in e.g. (5)) under a set of natural conditions on the strategies and transition rates. This motivates the study of the asymptotic deterministic model formulation and gives a transparent interpretation of the agent's control in the applied context. We derive theoretical results for the deterministic model: a verification theorem and an existence theorem for the coupled continuum of forward-backward ordinary differential equations (FBODEs) that characterize the finite state graphon game at equilibrium are proven. This is reminiscent of the mean field game framework, except that here the population is heterogeneous due to the graphon-based interactions. This makes the computation of solutions much more challenging. We then propose a machine learning method to solve the FBODE system. Finally, we consider a graphon game model for epidemic disease spread. Multiple test cases are solved with the proposed numerical method and the experimental results are discussed.
The outline of the rest of the paper is as follows: In Sect. 2, we introduce the model and analyze its deterministic formulation. In Sect. 3, we introduce the numerical approach and give experiment results. In Sect. 4, a theoretical framework for the model's probabilistic framework is presented and we rigorously define the graphon game. For the sake of conciseness, the proofs are postponed to the appendices.

Setup and Preliminaries
Let n ∈ N and let E be the finite set {1, . . . , n}. For each e ∈ E define the difference operator Δ e acting on functions on E by the formula [Δ e φ](e ) = φ(e ) − φ(e). We identify the set of probability measures on E, P(E), with the simplex Δ(E) := {x = (x 1 , . . . , x n ) ∈ R n + : i x i = 1} and endow it with the Euclidean distance. Throughout the paper, the notation P(·) will be used to denote the set of Borel probability measures.
Let T > 0 be a finite time horizon. A process ( f t ) t∈[0,T ] will be denoted with its bold letter symbol f . Let Let I be the unit interval equipped with the Euclidean distance. We denote by λ I and B(I ) the Lebesgue measure and Borel σ -field on I , respectively. The set I is indexing the continuum of players in the graphon game. Throughout the paper, we will employ the notation φ := (φ(x)) x∈I for functions with domain I . Furthermore, in most cases we will denote the index argument with a superscript: This paper studies games with heterogeneous interactions. When the players in the game interact, the weight they give to each others' action is parameterized by their indices. This weighted averaging is captured as the integration with respect to a graphon kernel function (see Section 1 for the discussion). Here, we give the most prevalent definition of a graphon. Definition 1 A graphon is a symmetric Borel-measurable function, w : The graphon induces an operator W from L 2 (I ) to itself: for any φ ∈ L 2 (I ),

Remark 1
The definition of a graphon varies somewhat in the literature. Some authors require neither symmetry nor a non-negative range. As mentioned in the introduction, the graphon as defined here can be used to represent the limit of a sequence of dense random graphs as the number of nodes (players) goes to infinity. For example, the constant graphon w(x, y) = p is in a sense the limit of a sequence of Erdős-Rényi graphs with parameter p ∈ [0, 1]. Conversely, random graphs can be sampled from a graphon in at least two different ways: either we sample points in I and construct a weighted graph whose weights are given by the graphon, or we sample points in I and then sample edges with probabilities given by the graphon.
A Q-matrix with real-valued entries q i, j , i, j ∈ E, is an n × n matrix with non-negative off-diagonal entries such that: In this paper, we will consider controlled Q-matrices with entries that depend on population aggregates. More specifically, we let for each We are going to work under the following assumption on the rates q x i, j : (ii) There is a finite constant C > 0, possibly depending on n, such that for p = 1, 2 and for Since the state space is finite and the rates are assumed to be uniformly bounded in Condition 1.(i), the Hölder-continuity assumption of Condition 1.(ii) is less restrictive than it otherwise would have been.

The Finite State Graphon Games Model for Epidemiology
In this section, we give a descriptive introduction to the epidemiological graphon game without going in to all the technical details. A rigorous mathematical motivation, built on the theory of Fubini extensions to accommodate for a continuum of independent jump processes, is presented in Sect. 4.
On a probability space (Ω, F , P), we consider a continuum of E-valued pure jump pro- The stochastic process X x models the state trajectory of player x. The initial state X x 0 is sampled from a distribution p x 0 ∈ P(E). Each player implements a strategy α x , a process taking values in the compact interval A ⊂ R described in more detail below. The players interact and each player's state trajectory is potentially influenced by the whole strategy profile (α x ) x∈I . To emphasize dependence we denote the state trajectory of player x as X α,x given a strategy The rate matrix is controlled by α x and influenced by the population aggregate Z α,x . The aggregates we consider are averages weighted by a graphon w, more specifically of the form In Sect. 4 we prove that the aggregate is a deterministic function of time, henceforth we write where w is a graphon. K will be called the impact function since it quantifies how much a player's joint state and control distribution impacts the aggregate variable. One example is the impact function in Sect. 1.1 where K (α, e) = α1 (e=I) , where the interpretation being that the aggregate is the averaged contact factor control of infected players. We have the following assumption on K : Condition 2 There exist finite constants L K , C K > 0 such that for all a, a ∈ A and e ∈ E: Our (for now formal) probabilistic definition of the interacting system of players is complete. The rigorous analysis of the system, the construction of a continuum of state trajectories, and conditions under which the aggregate is deterministic, i.e., of the form (9), is treated in detail in Sect. 4.
The key feature of the graphon game is that the aggregate variable is in general not the same for two distinct players. The players are therefore distinguishable and there is no "representative agent", as in MFGs 7 . As a direct consequence, there is no flow of player state distributions common to all players. Instead, each player has their private flow. Denote by p α,x (t, e) the probability that player x ∈ I that is in state e ∈ I at time t ∈ [0, T ], given that the population plays the strategy profile α. We shall argue that player x's state distribution flow p α,x solves the Kolmogorov forward equation with initial condition p α, with ρ α,y t being the joint probability law of control and state, (α y t , X α,y t− ). We turn our focus to the players' actions. We will make three standing assumptions that directly affect which strategies the players will choose. The first is that the environment that endogenously affects the players (but is known to the players) is varying smoothly over time, with no abrupt changes for example in lockdown penalties or expected recovery time. Secondly, if the players can impact their environment with their control, then the environment varies smoothly with their control too. For example, the risk of infection depends continuously on the agent's level of cautiousness. Finally, the players' strategies are decentralized, i.e., are unaffected by the transition of any agent other than themselves. Under these circumstances, the players have no apparent reason to discontinuously change their action over time except at times of transition between states. Such strategies (A-valued; decentralized; continuous in time between changes of the player's own state) will be called admissible and the set of admissible strategies denoted by A. The setting is further discussed in Sect. 4.
In this paper, we focus on the finite horizon problem where the cost is composed of two components: a running cost and a terminal cost. For each player x ∈ I , the assumptions on the conditions that the running and terminal cost functions, and g x : E × R → R, satisfy are given later in the text together with the theoretical results. The total expected cost to player x for playing the strategy σ ∈ A while the population plays the strategy profile α is As we shall see, a change in player x's control has no effect on the aggregate. Hence, the expected cost depends on the strategy profile only indirectly through the value of the aggregate variable. See Sect. 4 for the details. Therefore, hereinafter we shall use the notation J x (σ ; Z α,x ) for the right-hand side of (12). In light of this, we employ the following definition of a Nash equilibrium in the graphon game:

Definition 2
The strategy profile α is a Nash equilibrium if it is admissible and no player can gain from a unilateral deviation, i.e.,

Analysis of Finite State Graphon Games
By Definition 2, an admissible strategy profileα is a Nash equilibrium if there exists an • for all x ∈ I ,Ẑ x is the aggregate perceived by player x if the population uses strategy profileα. This alternative formulation has the advantage to split the characterization of the equilibrium into two parts and in the first part, the optimization problem faced by a single agent is performed with while the aggregate is fixed.
With a flow Z x being fixed, we define the value function of player x ∈ I as, for t ∈ [0, T ] and e ∈ E, To derive optimality conditions, we introduce the Hamiltonian of player x: where h ∈ R n and − → 1 e is the coordinate (row) vector in direction e in R n . We assume that A α → H x (t, e, z, h, α) admits a unique measurable minimizerâ x e (t, z, h) for all (t, z, h) and define the minimized Hamiltonian of player x: The dynamic programming principle of optimal control leads to the HJB equation for u x that readsu (15) can be equivalently written aṡ In the following theorem, we verify that the solution of the HJB equation indeed is the value function of the infinitesimal agent's control problem, and we provide an expression for an optimal Markovian control in terms of this value function and the aggregate.
is a continuously differentiable solution to the HJB equation (15), then u x is the value function of the optimal control problem when the flow Z x is given. Moreover, the function gives an optimal Markovian control.
Next, we prove the existence of a solution to the coupled Kolmogorov-HJB system at equilibrium. For that purpose we place the following condition: by Lâ(c), which can be bounded from above using smoothness properties of Q x and f x . More specifically, Lâ depends on the local Lipschitz coefficients of z → ∂ a Q x (a, z) and (t, z, h) → ∂ a f x (t, e, z, a), see the proof of Lemma 1). Recall that C K denotes the uniform upper bound of the impact function K guaranteed by Condition 2.

z, a) is continuously differentiable, and as a function of a it is strongly
is a probability mass function on E.

The Finite State Graphon Game for SIR with Contact Factor Control
Here, we introduce a model which we shall use as a test bed for the numerical algorithm presented in Sect. 3. It is inspired by the first example scenario in [2] and builds on the case discussed in Sect. For enforcement purposes, it also sets penalties for deviation from these levels. The cost has 3 components: The first component penalizes the agent for not following the regulator's recommended contact factor level, the second is the cost of treatment for an infected agent (this cost can be player specific due to individual differences in health care plan coverage, etc.), and the last one is the cost for being deceased. In this setting, the running cost is written as where c I and c D are nonnegative costs functions (of the player index). We set the terminal cost to be identically zero, g x (e, z) = 0 for all (e, z) ∈ E × R. The transition rate matrix for player x is given as: where β, γ , κ, ρ are nonnegative parameter functions, 0 ≤ ρ ≤ 1, determining the rates of infection, recovery, reinfection, and decease, and where in each row, the diagonal entry is the negative of the sum of the other terms on the same row. In line with the discussion in Sect. 1.1, the transition rate from state S to I depends on the agent's own decision and the aggregate variable 8 . Furthermore, when an infected agent transitions, she goes to state R with probability ρ and to state D with probability (1 − ρ). Then, for player x the optimality conditions yieldφ and the forward-backward graphon ODE system reads: We note that with a careful choice of β(x), c λ , and (λ (e) ) e∈E this system will satisfy the sufficient condition for existence from Theorem 2.

Numerical Approach
We rewrite the continuum of FBODEs (18) as the solution of a minimization problem: minimize where ( p θ , u θ ) solve the forward-forward continuous system of ODEs: 8 The rate is presented in (20) as β(x)αz. For an arbitrary impact function K this would lead to a violation of Condition 1.1. Here we are however considering contact factor control and aggregates of the form (11), which are bounded by Condition 2 and the compactness of A. The Q-matrix (20) can be modified without changing the model so that Condition 1.1 is satisfied; we defer from this for the sake of presentation.
This "shooting" strategy is reminiscent of the one used e.g. in [20,21,41] for stochastic optimal control problem and e.g. in [2,11] for mean field games in a numerical context. However, here we deal with a continuum of ODEs rather than a finite number of stochastic differential equations. Here θ is the parameter in the function ϕ replacing the initial condition of u. Typically, θ is a real-valued vector of dimension the number of degrees of freedom in the parametric function ϕ. In general, the true initial condition for u x is a nonlinear function with a potentially complicated shape. So we need to choose a rich enough class of parametric functions. In the implementation, we used a deep neural network with a feedforward architecture. See, e.g., [10] for a description of the feedforward neural network architecture we used in the implementation. Our strategy to find θ is to run a gradient-descent based method. To alleviate the computational cost and to introduce some randomness, at each iteration we replace the above cost J by an empirical average over a finite set of indices x, which is also used to approximate the value of the aggregate quantities. More precisely, for a finite set S of indices, we introduce where ( p θ,S , u θ,S ) solves the forward-forward (finite) system of ODEs:

Piecewise Constant Graphon
Let m 1 , m 2 , . . . , m K be non-negative numbers such that K k=1 m k = 1. We divide the player population into K groups, B 1 , . . . , B K , placing all players with index x ∈ [0, m 1 ) into group B 1 , etc. We assume that players belonging to the same group are indistinguishable. For example, all players within a group must have the same recovery rate γ and if x, x ∈ I are the indices of two players in the same group, then w(x, y) = w(x , y) for all y ∈ I . In this situation, we only need to specify the graphon's values on each block of indices corresponding to a group, since the graphon is a constant on each block. Let us identify the group B i with its index block (or set). We can compactly represent the interaction weights between the blocks with a connection matrix [w i j ] K i, j=1 , where w i j is the connection strength between players in block B i and players in block B j . Then, for all players x in block B i which is constant over x ∈ B i . This is a feature of that we sometimes see when the piecewise constant graphon is used. Furthermore, we assume β, γ , κ constant over each block but can differ between blocks. It opens up the possibility for us to solve the graphon game with classical numerical methods 9 and a way of evaluating the DeepGraphonGame algorithm. Turning to the remaining example set up, we specify the cost structure as a particular case of the general formulation of Sect. 2.4. We assume that the regulator has set λ (S,k) , λ (I,k) , λ (R,k) , k = 1, . . . , K , different for each block.
In this scenario, we first study the policy effects on the death ratio in the age groups. The policies compared are no lockdown (NL); quarantine for infected (QI); age specific lockdown (AL); full lockdown (FL). We can see in Fig. 1 that the death ratio decreases nearly 30% if infected individuals are quarantined, compared to no lockdown. Furthermore, if an age specific lockdown is implemented, we see that even more lives are saved while not deteriorating the economy. Zooming in on the comparison of no lockdown and quarantine  The lockdown is imposed by decreasing the λ S and λ I values of the corresponding age groups for infected (second row of Fig. 1), we note that susceptible individuals are using a smaller contact factor when there is no quarantine in place. They optimize their risk and hence want to be more cautious, since the risk of getting infected is higher (Tables 1 and 2). Secondly, in the same scenario, we model multiple cities with different attributes and study the effects of the travel restrictions. In this experiment, we compare a universal notravel policy to the policy where traveling in or out of one of the cities is restricted (totaling four policies in the comparison). City 1 is a highly populated city with a more contagious virus variant, city 2 also has this variant; however, it is a small city. City 3 is a highly populated city but with a less contagious virus variant. For visual simplification, we assume that there are no deaths (i.e. ρ = 0). In Fig. 2, we can see that the infected-density curve can be flattened the most if city 1 has travel restrictions. The reason is the existence of the more contagious variant and the large size of the city 1. We note that when this restriction is implemented the susceptible individuals feel relieved and increase their contact factor control (Tables 3 and  4). The connection matrix shown here is the case where there is a lockdown in City 1. When there is no lockdown the interaction weights are as follows: w(1, 1) = 1.0, w(1, 2) = 0.9 and w(1, 3) = 0.8 Table 4 Parameters used in the experiment with different cities Parameters All policies 40 0.1 1.0 0.9 1.0 10 1 1

Sanity Check for the Numerical Approach
Here, we test the DeepGraphonGame algorithm by comparing its solution to the solution obtained by solving the ODE system for the cities-example when city 1 has travel restrictions ( Table 4). As can be seen in Fig. 3, the DeepGraphonGame algorithm approximates the exact result well. A plot of the function x → u x (0, S), where u x is the numerically computed value function, can be seen on the right side of the bottom row in Fig. 3. We can clearly see that agents in the same block have the same u x (0, S) values. From this we infer that the DeepGraphonGame algorithm is preforming well when learning this piecewise constant function.

General Graphon
To show scalability of the proposed numerical approach now we focus on the second example in [2] with the SEIRD model where the state (E)xposed is added. An individual is in state E when infected but is not yet infectious. Hence, the agents evolve from S to E and then to I, and the infection rate from S to E depends on the proportion of the infected agents. The diagram of the dynamics can be seen in Fig. 4. The cost structure is similar to the one used in Sect. 2.4. After introducing the state E, we set In this example, we focus on an application where the agents are not homogeneous over blocks. The interaction strength between individual x and y is now given by the power law graphon: w(x, y) = (x y) −g where −∞ < g ≤ 0 is a constant 10 . Intuitively, the power law graphon models interactions in a population where a small number of individuals are responsible for a large number of the interactions. For example, a population with superspreaders 11 can be modeled with this graphon. The model with an underlying power law graphon interaction requires us to solve a continuum of coupled ODEs which is not computationally feasible. 10 We can realize that if g = 0, the setting is equivalent to a mean field game. 11 A superspreader is an infected person who is able to transmit the disease to a disproportionately high number of people. However, by using the DeepGraphonGame algorithm, the solution can be learned by using simulated particles, i.e. agents. According to CDC, COVID-19 reinfection is very rare 12 ; therefore, we assume that there is no reinfection (i.e. κ = 0). Furthermore, the recovery duration is around 10 days since the symptom onset 13 . For this reason, we assume that γ = 0.1 days −1 . According to the recent study conducted by Lauer et al. [56], an exposed person begins to show symptoms after around 5 days. Based on this observation, we choose = 0.2 days −1 . Finally, the Basic 12 https://www.cdc.gov/coronavirus/2019-ncov/your-health/reinfection.html (last accessed on 10 April, 2021) 13 https://www.cdc.gov/coronavirus/2019-ncov/hcp/duration-isolation.html (last accessed on 10 April, 2021)   Table 5). The experiment results for a sampled finite subset of the agent population are presented in Fig. 5. In the figure, each line corresponds to one agent and the color of the plot gets darker as the index of the agent (i.e. x) increases. Our first observation is that as the index of the agent increases, the aggregate also increases. In response to this high aggregate, the agent lower its contact rate, in order to protect itself. However, this protection is not enough to neutralize the effects of the high levels of the aggregate and the probability of the agent to be infected is still elevated.

The Probabilistic Approach to Finite State Graphon Games
This section contains a closer study of the continuum of interacting jump processes that constitute the graphon game dynamics. Going back to the informal discussion in the introductory Sect. 1.1, it is not clear that there would be an adequate law of large numbers so that (4) converges to (5) since the averaged random variables are dependent. A common approach in economic theory for this situation is to consider the continuum limit and average the continuum of random variables with respect to a non-atomic probability measure over I . Using the example from Sect. 1.1 again, a continuum limit of (4) is However, the integral in the expression above is ill-defined. There is an issue of constructing a continuum of independent random variables (here, that would be the driving Poisson noise) that are jointly measurable in the sample and the index. If the construction is done in the usual way via the Kolmogorov construction, then almost all random variables are essentially equal to an arbitrarily given function on the index space (i.e., as random variables they are constants). Hence, in any interesting case the function y → X y t− will not be measurable with respect to the Lebesgue measure. One solution proposed by economists is to extend the usual probability space to a so-called a Fubini extension [52], a probability space over Ω × I where Fubini's theorem holds. On the Fubini extension a continuum of random variables can be constructed that are essentially pairwise independent (e.p.i.; see Theorem 3 for the definition) and jointly measurable in sample and index. Moreover, there is hope for an exact law of large numbers [53], justifying our assumption about the determinism of the aggregate variable in the previous sections. We will construct a Fubini extension that carries a continuum of e.p.i. Poisson random measures. Then, each player path will be defined in the representation of the counting processes associated to the pure jump process as a stochastic integral with respect to a family of independent Poisson random measures with Lebesgue mean measure on R × [0, T ] as suggested in Skorokhod [51] and Grigelionis [34]. • for all G ∈ G such that G(G) < ∞, N · (G) is a Poisson random variable with rate G(G); • the random variables N · (G 1 ), . . . , N · (G n ) are mutually independent whenever the sets G 1 , . . . , G n ∈ G have finite G-measures and are disjoint; • for all ω ∈ H , N ω (·) is a measure on (Γ , G).

The Fubini Extension
In the model, state dynamics are given as E-valued jump processes. We will construct such a process from 2|E| − 1 = 2n − 1 independent Poisson random measures. The possible jumps will be those so that the state process jumps between two integers in E, at most n − 1 steps up or down. The initial state of player x ∈ I is randomly sampled from a pre-selected distribution p x 0 ∈ P(E). In order to model idiosyncratic random shocks affecting the dynamics of the individual states, we use the framework of Fubini extensions. It allows us to capture a form of independence for a continuum of random variables, while preserving joint measurability.
The following theorem summarizes the results by Sun and collaborators, see for example [53] and [54], which we use as a foundation for our model.

Theorem 3
There exists a probability space (I , I, λ) extending (I , B I , λ I ), a probability space (Ω, F , P), and a Fubini extension (Ω × I , F I, P λ) such that for any measurable mapping φ from (I , I, λ) to P(E × M 2n−1 ) there is an F I-measurable process f : Ω × I → E × M 2n−1 such that the random variables f x = f (·, x) are essentially pairwise independent (e.p.i.), i.e., for λ-a.e. x ∈ I , f x is independent of f y for λ-a.e. y ∈ I , and P • ( f x ) −1 = φ x for all x ∈ I .
where N is the probability law of the Poisson random measure introduced above, and p x 0 is the initial distribution of player x. By Theorem 3 (which holds for φ since E and M are Polish spaces) there exists a collection of random variables (ξ , N k ; k = −n + 1, . . . , n − 1) on a Fubini extension (Ω × I , F I, P λ), that are e.p.i. and φ x -distributed for all x ∈ I . With the model in Sect. 2 in mind, we assume that the mapping x → p x 0 is Lebesgue-measurable (this assumption is however not necessary for the analysis that follows).
We denote by L 2 (Ω × I ; D) the Bochner space of all (equivalence classes of) strongly (P λ, B(D))-measurable functions f : Ω × I → D for which We define L 2 (Ω × I ; C) in the same way, with C replacing D above. By e.g. For later reference, we define also the set L E as the subset of L 2 (Ω × I ; D) of P λ-a.e. D E -valued functions. One can show that L E is a closed subset of L 2 (Ω × I ; D), hence (L E , · L 2 (Ω×I ;D) ) is a complete metric space.

The Set of Admissible Strategies
We can now give a rigorous definition of the set of admissible strategy profiles. Recall that A is a compact subset of R. We will sometimes use the same notation A for the set of admissible control processes associated to an admissible strategy profile in feedback form α. At time t, for player x ∈ I , the value of such a control process is the action α(x, t, X x,α t− ) where X x,α t− is the state of player x just before time t. These control processes are predictable with respect to the filtration generated by the player's private state, and decentralized since they do not depend directly on other players' states.
The continuity in time is a strong assumption and prohibits the player from immediately reacting to abrupt changes in their environment. However, if a player transitions between two states at time t their control can be discontinuous at that time (as a multivariate function of time and state). In Sect. 2, a rationale was given for restricting our attention to such controls.

The Finite State Graphon Game in the Fubini Extension
We begin by describing an interacting system of a continuum of particles. First, we define the decoupled system where the aggregate variable vector has been "frozen". Then, we define the aggregate with a fixed point argument. Finally, we prove that the aggregate is in fact deterministic.
Consider the pure jump stochastic integral equation (here written formally) where is the rate of jumps from state i to state i + k given the action a ∈ A and the aggregate value z. The proposition below asserts that (26) has a unique solution in L E , the subset of L 2 (Ω × I ; D) defined in the end of Sect. 4.1.2.

Note that the quantity
appearing in the right-hand side of (26), is a counting process with intensity T ] so by construction the solution to (26) (granted by Proposition 1) is almost surely an E-valued pure jump process with intensity matrix . For a fixed admissible strategy profile α, consider now the coupled system The next theorem proves that (28) is well-posed with a unique solution in L 2 -sense. It further specifies the regularity of the solution: the aggregate variable Z α must P λ-a.s. be a deterministic and a continuous function of time.

Theorem 4 Let Condition 1 and 2 hold, and let α ∈ A.
(i) There exists a unique solution X α ∈ L E to (28). The corresponding aggregate Z α is a random variable in L 2 (Ω × I ; C). Theorem 4 justifies working with a model defined for all x ∈ I with a deterministic, continuous-in-time aggregate in Sect. 2. From here on, we will represent the L 2 -elements solving system (28) with the version defined for all x ∈ I and drop the check in the notation.

Remark 2
If admissible strategy profiles did not have the prescribed continuity property we could not expect the aggregate to be a continuous function of time. One example of such a case is found in [2] where a regulator imposes a penalty that is discontinuous in time, resulting in equilibrium controls and aggregates discontinuous in time. We leave the analysis of the more general case to future work.
We now turn to the notion of player costs and equilibrium. Denote by A the set of A-valued for every e ∈ E that. If the player population plays according to an admissible strategy profile α ∈ A and player x ∈ I decides to play strategy σ = (σ (t, X then (α −x , σ ) is an admissible strategy profile and the player's expected cost for using σ is In fact, A is the set of strategies a player can deviate to without destroying the admissibility of the strategy profile. Therefore, we say that ifα = (α x ) x∈I ∈ A satisfies theα is a Nash equilibrium of the graphon game. The dependence of the cost on the whole strategy profile is unnecessarily complicated, as the following reasoning shows. Notice that x ∈ I . The other players' actions appear in player x's cost indirectly, through the aggregate Z α,x , which is unaffected if one specific player changes control (it is an integral with respect to a non-atomic measure). Thus, we write J x (σ ; α) as J x (σ ; Z α,x ) a function taking an admissible strategy and an aggregate variable trajectory: In light of this, an equivalent definition of the Nash equilibrium is that a strategy profilê α = (α x ) x∈I is a Nash equilibrium in the graphon game if it satisfies further justifying the game setup in Sect. 2.

Conclusion and outlook
In this paper, we introduced stochastic graphon games in which the agents evolve in a finite state space. We provided optimality conditions in the form of a continuum of forwardbackward ODE system, for which we established existence of solutions. We proposed a numerical method based on a neural network approximation of the initial condition of the FBODE. We then applied our theoretical framework and numerical method to a class of models from epidemiology and we provided several test cases. From here, several directions can be considered for future work. An interesting aspect would be to incorporate a regulator with a Stackelberg type model as was done in [2] without graphon structure. The theoretical analysis would probably rely on a combination of the tools developed in the present work together with tools from optimal contract theory. However there would be some important challenges depending on the class of controls that are admissible for the regulator. Their controls could indeed lead to discontinuities in the incentives to the population, which would raise subtle measurability questions. This is left for future work. Another direction is to consider more realistic epidemiological models (e.g., with more compartments). Such models would be more complex and we expect our proposed machine learning numerical method to be helpful from this point of view. Furthermore, to be able to use graphon games to make epidemiological predictions, it would be interesting to investigate further how to use real data in the model and in the numerical method.
By the assumption on the strict convexity of f , we have t, e, z, a) .
So, combining the above inequalities and the property |∂ a Q x (a, z) − ∂ a Q x (a , z)| = 0, we get where C can depend on q max , C z , and C h . We conclude by using the continuity and local Lipschitz continuity properties of − → 1 e ∂ a Q x and ∂ a f x .

A.3 Theorem 2
Step 1: Definition of the solution space We start by letting, for every C 1 > 0, K C 1 be the closed ball of continuous functions (u, p) from [0, T ] into L 2 (I × E) × L 2 (I × E) such that for all (t, x) ∈ [0, T ] × I , p x (t, e) ≥ 0 for e ∈ E and e∈E p x (t, e) = 1, and such that (u, p) is bounded by C 1 in uniform norm. In other words, K C 1 is the subset of C([0, T ]; L 2 (I × E) × L 2 (I × E)) for which the second component is a probability on E and for which the uniform norm is bounded by C 1 , whose value will be fixed in Step 5 below.
Step 2: Definition of the aggregate mapping For each (u, p) ·)), e p y (t, e) dy 0≤t≤T ,x∈I .
We now prove that, if we choose the space of aggregates properly, Φ (u, p) has a unique fixed point, say (Ẑ x t ) 0≤t≤T ,x∈I , which depends continuously in (u, p) ∈ K C 1 . Indeed, let z, z ∈ Z where a.e. x ∈ I }. In light of Condition 2, Φ (u, p) (Z) ⊂ Z. Moreover, Z is a closed subset of the Banach space C([0, T ]; L 2 (I ; R)), hence a complete metric space. Using Cauchy-Schwarz inequality and the Lipschitz continuity of K andâ (given by Condition 2 and Lemma 1), we get where we used the fact that |I | = 1 and C Φ (·) := w L 2 (I ×I ) L K Lâ(·). Recall that we assume C Φ (C K ) < 1. With Banach fixed point theorem we conclude that there exists a unique fixed point in Z to Φ (u, p) . We denote itẐ (u, p) = (Ẑ (u, p),x t ) 0≤t≤T ,x∈I .
Step 3: Solving the Kolmogorov equation GivenẐ (u, p) and u, we solve the Kolmogorov equation and we get the solutionp. Existence and uniqueness of the solutionp is provided by the Cauchy-Lipschitz-Picard theorem; see, e.g., [5,Theorem 7.3] (viewing q as a linear operator acting on the Banach space L 2 (I × E)). Furthermore, given Condition 1.1, the time derivative ofp is bounded. Therefore, we conclude thatp is equicontinuous.
Step 4: Solving the HJB equation GivenẐ (u, p) andp, we solve the HJB equation and we get the solutionû. Here again, existence and uniqueness of the solutionû is provided by the Cauchy-Lipschitz-Picard theorem (viewingĤ as a Lipschitz operator acting on the Banach space L 2 (I × E)). Furthermore, there is a uniform bound on the time derivative ofû since the HamiltonianĤ is bounded given Condition 3.3 and Condition 1.1. Henceû is equicontinuous.
Step 5: Application of Schauder's theorem Let us call Ψ the mapping constructed by the above steps, namely, Ψ : K C 1 (u, p) → (û,p). By steps 3 and 4 above, if we choose C 1 large enough, Ψ maps K C 1 onto K C 1 , so it is well-defined. Furthermore, by the same steps and the Arzela-Ascoli theorem Ψ (K C 1 ) is compact. Finally we argue the continuity as follows: We first show the continuity ofẐ (u, p) in u and p. Consider a sequence (u n , p n ) n such that (u n , p n ) ∈ K C 1 for every n, and lim n→∞ (u n , p n ) = (u, p). We denote Z n := Z (u n , p n ) and prove below that lim n→∞Ẑ n =Ẑ . By Lipschitz continuity of K andâ: C p p n,x (t, ·) − p x (t, ·) 2 + C u u n,x (t, ·) − u x (t, ·) 2 dx, which tends to 0 as n → ∞. Next, we study the continuity ofp. We have, for s ∈ [0, T ]: − q x e ,e (â(t, e ,Ẑ x t , u x (t, ·)),Ẑ x t ) 2 |p n,x (t, e )| 2 + |q x e ,e (â(t, e ,Ẑ x t , u x (t, ·)),Ẑ x t )| 2 p n,x (t, e ) −p x (t, e ) 2 dtdx B Proofs for Sect. 4
Following similar lines of proof as in Proposition 1, get the initial P λ-a.e. estimates Let t j ∈ [0, T ], j ∈ N, be a sequence converging to t * ∈ [0, T ]. Without loss of generality, assume that the sequence is non-decreasing. Recall that q max denotes the uniform upper bound for the intensity rates, see Condition 1.1. For (ω, x) ∈ (Ω × I ),  . (ω, x) ∈ Ω × I . Measurability and squareintegrability of U α z follows by the same lines of proof as Lemma 1 in [3].