algorithmic Yao Xuan<EMAIL_ADDRESS>
Robert Balkin<EMAIL_ADDRESS>
Department of Mathematics, University of California, Santa Barbara, CA
93106-3080, USA and Jiequn Han<EMAIL_ADDRESS>
Department of Mathematics, Princeton University, Princeton, NJ 08544-1000, USA
and Ruimeng Hu<EMAIL_ADDRESS>
Department of Mathematics and Department of Statistics and Applied
Probability, University of California, Santa Barbara, CA 93106-3080, USA and
Hector D. Ceniceros<EMAIL_ADDRESS>
Department of Mathematics, University of California, Santa Barbara, CA
93106-3080, USA

# Optimal Policies for a Pandemic: A Stochastic Game Approach and a Deep
Learning Algorithm

###### Abstract

Game theory has been an effective tool in the control of disease spread and in
suggesting optimal policies at both individual and area levels. In this paper,
we propose a multi-region SEIR model based on stochastic differential game
theory, aiming to formulate optimal regional policies for infectious diseases.
Specifically, we enhance the standard epidemic SEIR model by taking into
account the social and health policies issued by multiple region planners.
This enhancement makes the model more realistic and powerful. However, it also
introduces a formidable computational challenge due to the high dimensionality
of the solution space brought by the presence of multiple regions. This
significant numerical difficulty of the model structure motivates us to
generalize the deep fictitious algorithm introduced in [Han and Hu, MSML2020,
pp.221–245, PMLR, 2020] and develop an improved algorithm to overcome the
curse of dimensionality. We apply the proposed model and algorithm to study
the COVID-19 pandemic in three states: New York, New Jersey and Pennsylvania.
The model parameters are estimated from real data posted by the Centers for
Disease Control and Prevention (CDC). We are able to show the effects of the
lockdown/travel ban policy on the spread of COVID-19 for each state and how
their policies affect each other.

###### keywords:

Stochastic differential game, pandemic, optimal policy, enhanced deep
fictitious play

##  1 Introduction

The pandemic of coronavirus disease 2019 (COVID-19) has brought a huge impact
on our lives. Based on the CDC Data Tracker, as of early December 2020, there
have been more than 15 million confirmed cases of infection and more than 290
thousand cases of death in the United States. Needless to say, the economic
impact has also been catastrophic, resulting in unprecedented unemployment and
the bankruptcy of many restaurants, recreation centers, shopping malls, etc.

In a classic, compartmental epidemiological model each individual is assigned
a label, e.g., Susceptible, Exposed, Infectious, Removed, Vaccinated. The
labels’ order shows the flow patterns between the compartments (SIR, SEIR,
SIRV models). Other approaches include network models, which explicitly
include the interaction of individuals, in addition to the modeling of each
individual’s dynamics, and agent-based models that are useful in informing
decision making when accurately calibrated. Moreover, the consideration of
pharmaceutical and/or non-pharmaceutical intervention policies naturally
couples game theory to epidemiological models by controlling when and how the
game is played in such models. For example, in some early studies [Bauch et
al.(2003)Bauch, Galvani, and Earn, Bauch and Earn(2004)], one can use non-
repeated games to incorporate game theory into modeling at the individual
level, where individuals (known as “players” in the game theory) maximize
their gain by weighing the costs and benefits of different strategies. We
refer to the review paper [Chang et al.(2020)Chang, Piraveenan, Pattison, and
Prokopenko] and the references therein for more details.

Differential games, initiated by [Isaacs(1965)], as an offspring of game
theory and optimal control, provide modeling and analysis of conflict in the
context of a dynamical system. They have been intensively employed across many
disciplines, including management science, economics, social science, biology,
military, etc. One of the core objectives in differential games is to compute
Nash equilibria that refer to strategies by which no player has an incentive
to deviate. However, a major bottleneck comes from the notorious
intractability of $N$-player games, and the direct computation of Nash
equilibria is extremely time-consuming and memory demanding. In a series of
recent works by [Hu(2020), Han and Hu(2020), Han et al.(2020a)Han, Hu, and
Long, Han et al.(2020b)Han, Hu, and Long], the deep fictitious play (DFP)
theory and algorithms were developed for stochastic differential games (SDG)
with a large number of heterogeneous players. The DFP framework embeds the
fictitious play idea, introduced by [Brown(1949), Brown(1951)], into designed
architectures of deep neural networks to produce accurate and parallelizable
algorithms with convergence analysis, and resolve the intractability issue
(curse of dimensionality) caused by the complex modeling and underlying high-
dimensional space in SDG.

Building from the DFP theory and algorithms for computing Nash equilibria in
SDG, we propose here to strengthen the classical SEIR model by taking into
account the social and health policies issued by multiple region planners. We
call this new model a stochastic multi-region SEIR model because it couples
the stochastic differential game theory with the SEIR model, making it more
realistic and powerful. The computational challenge introduced by the high-
dimensionality of the multi-region solution space is addressed by generalizing
the deep fictitious algorithm proposed by [Han and Hu(2020)]. This new
approach leads to an enhanced deep fictitious play algorithm to overcome the
curse of dimensionality and further reduce the computational complexity. To
showcase the performance of the proposed model and algorithm, we apply them to
a case study of the COVID-19 pandemic in three states: New York (NY), New
Jersey (NJ), and Pennsylvania (PA). We present the optimal lockdown policy
corresponding to the Nash equilibrium of the multi-region SEIR model. We
remark that our work is not to predict a pandemic, but to provide a game-
themed framework, a deep learning algorithm and possible outcomes for
competitive region planners. We hope that information can provide some
qualitative guidance for policymakers on the impact of certain policies. The
parameters used in the numerical experiments are based on the current
knowledge of the Coronavirus, which is still under development. In practice,
at the beginning of the Coronavirus, a governor might not be able to trace the
infections fully, and infected cases may not be fully identified. All these
may lead to the inaccuracy of parameter estimations despite our using of the
best available data. Therefore, it may be hard to match the predicted results
with practical observation.

The rest of the paper is organized as follows. In Section 2, we propose a
novel multi-region epidemiological model and explain how social and health
policies issued by region planners are integrated into dynamics, leading to a
game feature of the problem. We explain the numerical challenge in Section 3
and propose an enhanced deep fictitious play algorithm for memory and
computational efficiency. Section 4 focuses on a NY-NJ-PA case study with
detailed discussions on parameter choices and the resulted optimal policies.
We present some concluding remarks in Section 5 and provide technical details
in the appendices.

##  2 Mathematical modeling: A multi-region SEIR model

We consider a pandemic spreading in $N$ geographical regions, and each planner
controls the loss of her region by implementing some policies. We aim to study
how the region planners’ policies affect each other, and the equilibrium
policies.

Let us start with a modified version of the very-known epidemic SEIR model
(cf. [Liu et al.(1987)Liu, Hethcote, and Levin]), where each region’s
population is assigned to compartments with four labels: Susceptible, Exposed,
Infectious, and Removed. Individuals with different labels denote $S$: those
who are not yet infected; $E$: who have been infected but are not yet
infectious themselves; $I$: who have been infected and are capable of
spreading the disease to those in the susceptible category, and $R$: who have
been infected and then removed from the disease due to recovery or death. The
region planners can issue certain policies to mitigate the pandemic, for
instance, policies that can help reduce the transmission rates and death
rates. Mathematically, denote by $S^{n}_{t},E^{n}_{t},I^{n}_{t},R^{n}_{t}$ the
_proportion_ of population in the four compartments of region $n$ at time $t$.
We consider the following stochastic multi-region SEIR model:

$\displaystyle\,\mathrm{d}S^{n}_{t}=-\sum_{k=1}^{N}\beta^{nk}S^{n}_{t}I^{k}_{t}(1-\theta\ell^{n}_{t})(1-\theta\ell^{k}_{t})\,\mathrm{d}t-v(h^{n}_{t})S^{n}_{t}\,\mathrm{d}t-\sigma_{s_{n}}S^{n}_{t}\,\mathrm{d}W^{s_{n}}_{t},$
(1)
$\displaystyle\,\mathrm{d}E^{n}_{t}=\sum_{k=1}^{N}\beta^{nk}S^{n}_{t}I^{k}_{t}(1-\theta\ell^{n}_{t})(1-\theta\ell^{k}_{t})\,\mathrm{d}t-\gamma
E^{n}_{t}\,\mathrm{d}t+\sigma_{s_{n}}S^{n}_{t}\,\mathrm{d}W^{s_{n}}_{t}-\sigma_{e_{n}}E^{n}_{t}\,\mathrm{d}W^{e_{n}}_{t},$
(2) $\displaystyle\,\mathrm{d}I^{n}_{t}=(\gamma
E^{n}_{t}-\lambda(h^{n}_{t})I^{n}_{t})\,\mathrm{d}t+\sigma_{e_{n}}E^{n}_{t}\,\mathrm{d}W^{e_{n}}_{t},$
(3)
$\displaystyle\,\mathrm{d}R^{n}_{t}=\lambda(h^{n}_{t})I^{n}_{t}\,\mathrm{d}t+v(h^{n}_{t})S^{n}_{t}\,\mathrm{d}t,\quad
n\in\mathcal{N}:=\\{1,2,\ldots,N\\},$ (4)

where $\bm{\ell}_{t}\equiv(\ell^{1}_{t},\ldots,\ell^{N}_{t})$ and
$\bm{h}_{t}\equiv(h^{1}_{t},\ldots,h^{N}_{t})$ are policies chosen by the
region planners at time $t$. Each planner $n$ seeks to minimize its region’s
cost within a period $[0,T]$:

$J^{n}(\bm{\ell},\bm{h}):=\mathbb{E}\bigg{[}\int_{0}^{T}e^{-rt}P^{n}\big{[}(S^{n}_{t}+E^{n}_{t}+I^{n}_{t})\ell^{n}_{t}w+a(\kappa
I^{n}_{t}\chi+pI^{n}_{t}c)\big{]}+e^{-rt}\eta(h^{n}_{t})^{2}\,\mathrm{d}t\bigg{]}.$
(5)

We now give detailed description of this model (1)–(5):

  1. S:

$\beta^{nk}$ denotes the average number of contacts per person per time. The
transition rate between $S^{n}$ and $E^{n}$ due to contacting infectious
people in the region $k$ is proportional to the fraction of those contacts
between an infectious and a susceptible individual, which result in the
susceptible one becoming infected, i.e., $\beta S^{n}_{t}I^{k}_{t}$. Although
some regions may not be geographically connected, the transmission between the
two is still possible due to air travels but is less intensive than the
transmission within the region, i.e., $\beta^{nk}>0$ and
$\beta^{nn}\gg\beta^{nk}$ for all $k\neq n$.

$\ell^{n}_{t}\in[0,1]$ denotes the decision of the planner $n$ on the fraction
of population being locked down at time $t$. We assume that those in lockdown
cannot be infected. However, the policy may only be partially effective as
essential activities (food production and distribution, health, and basic
services) have to continue. Here we use $\theta\in[0,1]$ to measure this
effectiveness, and the transition rate under the policy $\bm{\ell}$ thus
become
$\beta^{nk}S^{n}_{t}I^{k}_{t}(1-\theta\ell^{n}_{t})(1-\theta\ell^{k}_{t})$.
The case $\theta=1$ means the policy is fully effective.

$h^{n}_{t}\in[0,1]$ denotes the effort the planner $n$ decide to put into the
health system, which we refer as _health policy_. It will influence the
vaccination availability $v(\cdot)$ and the recovery rate $\lambda(\cdot)$ of
this model.

$v(h^{n}_{t})$ denotes the vaccination availability of region $n$ at time $t$.
Once vaccinated, the susceptible individuals $v(h^{n}_{t})S^{n}_{t}$ become
immune to the disease, and join the removed category $R^{n}_{t}$. We model it
as an increasing function of $h^{n}_{t}$, and if the vaccine has not been
developed yet, we can define $v(x)=0$ for $x\leq\overline{h}$.

  2. E:

$\gamma$ describes the latent period when the person has been infected but not
infectious yet. It is the inverse of the average latent time, and we assume
$\gamma$ to be identical across all regions. The transition between $E^{n}$
and $I^{n}$ is proportional to the fraction of exposed, i.e., $\gamma
E^{n}_{t}$.

  3. I:

$\lambda(\cdot)$ represents the recovery rate. For the infected individuals, a
fraction $\lambda(h^{n})I^{n}$ (including both death and recovery from the
infection) joins the removed category $R^{n}$ per time unit. The rate is
determined by the average duration of infection $D$. We model the duration (so
does the recovery rate) related to the health policy $h^{n}_{t}$ decided by
its planner. The more effort put into the region (i.e., expanding hospital
capacity, creating more drive-thru testing sites), the more clinical resources
the region will have and the more resources will be accessible by patients,
which could accelerate the recovery and slow down death. The death rate,
denoted by $\kappa(\cdot)$, is crucial for computing the cost of the region
$n$; see the next item.

  4. Cost:

Each region planner faces four types of cost. One is the economic activity
loss due to the lockdown policy, where $w$ is the productivity rate per
individual, and $P^{n}$ is the population of the region $n$. The second one is
due to the death of infected individuals. Here $\kappa$ is the death rate
which we assume for simplicity to be constant, and $\chi$ denotes the cost of
each death. The hyperparameter $a$ describes how planners weigh deaths and
infections comparing to other costs. The third one is the inpatient cost,
where $p$ is the hospitalization rate, and $c$ is the cost per inpatient day.
The last term $\eta(h^{n}_{t})^{2}$ is the grants putting into the health
system. We choose a quadratic form to account for diminishing marginal utility
(view it from $\eta(h^{n}_{t})^{2}$ to $h^{n}_{t}$). All costs are discounted
by an exponential function $e^{-rt}$, where $r$ is the risk-free interest
rate, to take into account the time preference. Note that region $n$’s cost
depends on all regions’ policies $(\bm{\ell},\bm{h})$, as $\\{I^{k},k\neq
n\\}$ appearing in the dynamics of $S^{n}$. Thus we write is
$J^{n}(\bm{\ell},\bm{h})$.

The choices of epidemiological parameters will be discussed in Section 4.1.
Next, we summarize the key assumptions in the above model:

  1. 1.

The dynamics of an epidemic are much faster than the vital (birth and death)
dynamics. So vital dynamics are omitted in the above model.

  2. 2.

The planning is of a short horizon and will be adjusted frequently as the
epidemic develop. For simplicity, we assume this is no migration between
regions over the time $[0,T]$.

  3. 3.

Individuals who once recovered from the disease, are immune and free of
lockdown policy.

  4. 4.

The dynamics obeys the conservation law:
$S^{n}_{t}+E^{n}_{t}+I^{n}_{t}+R^{n}_{t}=P^{n}$. This means that the process
$R^{n}$ is redundant.

  5. 5.

The dynamics of $S$, $E$ and $I$ are subjected to random noise, to account for
the noise introduced during data recording, false-positive/negative test
results, exceptional cases when recovered individuals become susceptible
again, minor individual differences in the latent period, etc.

  6. 6.

Individuals who are not under lockdown have the same productivity, no matter
their categories. We assume this for simplicity remark that this can be
improved by assigning different productivity to individuals with or without
symptoms.

The above modeling and objectives can be viewed as a stochastic differential
game between $N$ players111Henceforth, we shall use _planner_ and _player_
interchangeably. Here, we view the whole problem as a non-cooperative game, as
many regions make decisions individually and indeed even compete for scarce
resources (frontline workers, personal protective equipment, etc.) during the
outbreak. Each player $n$ controls her states $(S^{n},E^{n},I^{n},R^{n})$
through her strategy $(\ell^{n},h^{n})$ in order to minimize the associated
cost $J^{n}$. The optimizers then are interpreted as the optimal lockdown
policy and optimal effort putting into the health system.

For a non-cooperative game, one usually refers to Nash equilibrium as a notion
of optimality. For completeness, we review the definition here.

###### Definition 2.1

A Nash equilibrium is a tuple
$(\bm{\ell}^{\ast},\bm{h}^{\ast})=(\ell^{1,\ast},h^{1,\ast},\ldots,\ell^{N,\ast},h^{N,\ast})\in\mathbb{A}^{N}$
such that

$\forall n\in\mathcal{N},\emph{ and }(\ell^{n},h^{n})\in\mathbb{A},\quad
J^{n}(\bm{\ell}^{\ast},\bm{h}^{\ast})\leq
J^{n}((\bm{\ell}^{-n,\ast},\ell^{n}),(\bm{h}^{-n,\ast},h^{n})),$ (6)

where $\bm{\ell}^{-n,\ast}$ represents strategies of players other than the
$n$-th one:

$\bm{\ell}^{-n,\ast}:=[\ell^{1,\ast},\ldots,\ell^{n-1,\ast},\ell^{n+1,\ast},\ldots,\ell^{N,\ast}]\in\mathbb{A}^{N-1},$
(7)

$\mathbb{A}$ denotes the set of admissible strategies for each player and
$\mathbb{A}^{N}$ is the produce of $N$ copies of $\mathbb{A}$. For simplicity,
we have assumed all players taking actions in the same space.

In the sequel, to fix the notations, we shall use

  * •

a regular character with a superscript $n$ for an object from player $n$;

  * •

a boldface character for a collection of objects from all players, i.e.,
$\bm{S}_{t}\equiv[S_{t}^{1},\ldots,S_{t}^{N}]^{\operatorname{T}}$;

  * •

a boldface character with a superscript $-n$ for a collection of objects from
all players except $n$, i.e.,
$\bm{S}_{t}^{-n}\equiv[S_{t}^{1},\ldots,S_{t}^{n-1},S_{t}^{n+1},\ldots,S_{t}^{N}]^{\operatorname{T}}$.

A Markovian Nash equilibrium is a Nash equilibrium defined above with
$\mathbb{A}$ being the set of Borel measurable functions:
$(\ell,h):[0,T]\times\mathbb{R}^{3N}\to[0,1]^{2}$. In other words, the
policies $(\ell_{t}^{n},h_{t}^{n})$ at time $t$ are functions of the time $t$
and the current values of all players’ state processes
$(\bm{S}_{t},\bm{E}_{t},\bm{I}_{t})$. We omit the dependence on $\bm{R}_{t}$
as it is redundant as a consequence of the conservation law.

We derive below the Hamilton-Jacobi-Bellman (HJB) equations characterizing the
Markovian Nash equilibrium. To simplify the notation, we first rewrite the
dynamics of $(\bm{S}_{t},\bm{E}_{t},\bm{I}_{t})$ defined in (1)–(3) into a
vector form
$\bm{X}_{t}\equiv[\bm{S}_{t},\bm{E}_{t},\bm{I}_{t}]^{\operatorname{T}}\equiv[S_{t}^{1},\cdots,S_{t}^{N},E_{t}^{1},\cdots,E_{t}^{N},I_{t}^{1},\cdots,I_{t}^{N}]^{\operatorname{T}}\in\mathbb{R}^{3N}$.
Again, we shall drop the redundant process $\bm{R}_{t}$. The dynamics of
$\bm{X}_{t}$ reads:

$\,\mathrm{d}\bm{X}_{t}=b(t,\bm{X}_{t},\bm{\ell}(t,\bm{X}_{t}),\bm{h}(t,\bm{X}_{t}))\,\mathrm{d}t+\Sigma(\bm{X}_{t})\,\mathrm{d}\bm{W}_{t},$
(8)

where $b,\Sigma$ are deterministic functions in $\mathbb{R}^{3N}$ and
$\mathbb{R}^{3N\times 2N}$, and $\\{\bm{W}_{t}\\}_{0\leq t\leq T}$ is a
$2N$-dimensional standard Brownian motion. Each player $n$ aims to minimize
the expected running cost

$\mathbb{E}\left[\int_{0}^{T}f^{n}(t,\bm{X}_{t},\ell^{n}(t,\bm{X}_{t}),h^{n}(t,\bm{X}_{t}))\,\mathrm{d}t\right].$
(9)

We defer the specific definitions of $b,\Sigma,\bm{W}$, and $f^{n}$ to
Appendix A.1 to facilitate the exposition. We now define the value function of
player $n$ by

$V^{n}(t,{\bm{x}})=\inf_{(\ell^{n},h^{n})\in\mathbb{A}}\mathbb{E}\left[\int_{t}^{T}f^{n}(s,\bm{X}_{s},\ell^{n}(s,\bm{X}_{s}),h^{n}(s,\bm{X}_{s}))\,\mathrm{d}s|\bm{X}_{t}={\bm{x}}\right].$
(10)

By dynamic programming, it solves the following HJB system

$\begin{dcases}\partial_{t}V^{n}+\inf_{(\ell^{n},h^{n})\in[0,1]^{2}}H^{n}(t,\bm{x},(\bm{\ell},\bm{h})(t,{\bm{x}}),\nabla_{{\bm{x}}}V^{n})+\frac{1}{2}\text{Tr}(\Sigma({\bm{x}})^{\operatorname{T}}\text{Hess}_{{\bm{x}}}V^{n}\Sigma({\bm{x}}))=0,\\\
V^{n}(T,{\bm{x}})=0,\quad n\in\mathcal{N},\end{dcases}$ (11)

where $H^{n}$ is the usual Hamiltonian defined by

$H^{n}(t,\bm{x},\bm{\ell},\bm{h},\bm{p})=b(t,\bm{x},\bm{\ell},\bm{h})\cdot\bm{p}+f^{n}(t,\bm{x},\ell^{n},h^{n}),$
(12)

$\partial_{t}$ denotes the time derivative, $\nabla_{{\bm{x}}}V$ and
$\text{Hess}_{{\bm{x}}}V$ denote the gradient and the Hessian of the function
$V$ with respect to $\bm{x}$, respectively, and Tr stands for the trace of a
matrix.

Finding the optimal policies for $N$ regions is equivalent to solving
$N$-coupled $3N+1$ dimensional nonlinear equations (11). For example, when
$N=3$, each PDE is 10-dimensional and conventional methods start to lose their
efficiency. The recently proposed deep learning algorithm in [Han and
Hu(2020)], known as deep fictitious play (DFP), has shown excellent numerical
performance in solving high-dimensional, coupled HJB equations with
convergence analysis (Han et al.(2020a)Han, Hu, and Long). In the next
section, we will first review and then propose an enhanced version of DFP to
tackle some new issues.

##  3 Numerical methodology: Enhanced deep fictitious algorithm

We first briefly review the deep fictitious play (DFP) algorithm for solving
equation (11) and we refer readers to Han and Hu(2020) for full details. With
the idea of fictitious play, DFP recasts the $N$-player game into $N$
decoupled optimization problems, which are solved repeatedly stage by stage.
Each individual problem is solved by the deep BSDE method (E et al.(2017)E,
Han, and Jentzen; Han et al.(2018)Han, Jentzen, and E). The algorithm starts
with some initial guess $(\bm{\ell}^{0},\bm{h}^{0})$, where the superscript
$0$ stands for stage $0$. At the $(m+1)^{th}$ stage, given the optimal
policies $(\bm{\ell}^{m},\bm{h}^{m})$ at the previous stage, the algorithm
solves the following PDEs

$\begin{dcases}\partial_{t}V^{n,m+1}+\inf_{(\ell^{n},h^{n})\in[0,1]^{2}}H^{n}(t,\bm{x},(\ell^{n},\bm{\ell}^{-n,m},h^{n},\bm{h}^{-n,m})(t,{\bm{x}}),\nabla_{{\bm{x}}}V^{n,m+1})\\\
\hskip
50.0pt+\frac{1}{2}\text{Tr}(\Sigma({\bm{x}})^{\operatorname{T}}\text{Hess}_{{\bm{x}}}V^{n,m+1}\Sigma({\bm{x}}))=0,\\\
V^{n,m+1}(T,{\bm{x}})=0,\quad n\in\mathcal{N},\end{dcases}$ (13)

and obtains the $(m+1)^{th}$ stage’s optimal strategy by:

$(\ell^{n,m+1},h^{n,m+1})(t,\bm{x})=\operatorname*{arg\,min}_{(\ell^{n},h^{n})\in[0,1]^{2}}H^{n}(t,\bm{x},(\ell^{n},\bm{\ell}^{-n,m},h^{n},\bm{h}^{-n,m})(t,{\bm{x}}),\nabla_{{\bm{x}}}V^{n,m+1}(t,{\bm{x}})).$
(14)

Here, $(\bm{\ell}^{-n,m},\bm{h}^{-n,m})$ stands for others’ optimal policies
from the $m^{th}$ stage and are considered to be fixed functions when solving
the PDE at the current stage. In the sequel, to simplify notations we omit the
stage label $m$ in the superscript when there is no risk of confusion. To
solve (13) at each stage, it is first rewritten in the DFP as

$\displaystyle\partial_{t}V^{n}+\frac{1}{2}\text{Tr}(\Sigma({\bm{x}})^{\operatorname{T}}\text{Hess}_{{\bm{x}}}V^{n}\Sigma({\bm{x}}))$
$\displaystyle+\mu^{n}(t,{\bm{x}};\bm{\ell}^{-n},\bm{h}^{-n})\cdot\nabla_{{\bm{x}}}V^{n}$
(15)
$\displaystyle+g^{n}(t,{\bm{x}},\Sigma({\bm{x}})^{\operatorname{T}}\nabla_{{\bm{x}}}V^{n};\bm{\ell}^{-n},\bm{h}^{-n})=0,$
(16)

with some functions $\mu^{n}$ and $g^{n}$. The solution is then approximated
by solving the equivalent BSDE
$(\bm{X}_{t}^{n},Y_{t}^{n},Z_{t}^{n})\in\mathbb{R}^{3N}\times\mathbb{R}\times\mathbb{R}^{2N}$:

$\displaystyle\bm{X}_{t}^{n}={\bm{x}}_{0}+\int_{0}^{t}\mu^{n}(s,\bm{X}_{s}^{n};(\bm{\ell}^{-n},\bm{h}^{-n})(s,\bm{X}_{s}^{n}))\,\mathrm{d}s+\int_{0}^{t}\Sigma(\bm{X}_{s}^{n})\,\mathrm{d}\bm{W}_{s},$
(17) $\displaystyle
Y_{t}^{n}=\int_{t}^{T}g^{n}(s,\bm{X}_{s}^{n},Z_{s}^{n};(\bm{\ell}^{-n},\bm{h}^{-n})(s,\bm{X}_{s}^{n}))\,\mathrm{d}s-\int_{t}^{T}(Z_{s}^{n})^{\operatorname{T}}\,\mathrm{d}\bm{W}_{s},$
(18)

in the sense of (cf. Pardoux and Peng(1992); El Karoui et al.(1997)El Karoui,
Peng, and Quenez; Pardoux and Tang(1999))

$Y_{t}^{n}=V^{n}(t,\bm{X}_{t}^{n})\quad\text{and}\quad
Z_{t}^{n}=\Sigma(\bm{X}_{t}^{n})^{\operatorname{T}}\nabla_{{\bm{x}}}V^{n}(t,\bm{X}_{t}^{n}).$
(19)

The high-dimensional BSDE (17)–(18) is tackled by the deep BSDE method
proposed in E et al.(2017)E, Han, and Jentzen; Han et al.(2018)Han, Jentzen,
and E.

In Han and Hu(2020), the algorithm solves the BSDE by parametrizing
$V^{n}(t,{\bm{x}})$ using neural networks (NN) and then obtains the
approximate optimal policy by plugging the NN outputs into (14). For memory
efficiency, the algorithm only stores the NNs’ parameters at the current and
the one-step previous stages. This strategy works well for games like the
linear-quadratic game, but it would be ineffective if $\bm{\ell}^{-n}$ or
$\bm{h}^{-n}$ explicitly appears in the minimizer in (14). In this case, when
evaluating others’ strategy $\bm{\ell}^{-n}$ or $\bm{h}^{-n}$ at stage $m$, it
does not only need NNs at stage $m$ but also at stages $m-1$, $m-2$, …, $0$.
This means one needs to store NNs’ parameters for all the previous stages from
$1,\ldots,m$, and evaluate the associated output. Therefore, the time
complexity of evaluating $(\bm{\ell}^{-n},\bm{h}^{-n})(s,\bm{X}_{s}^{n})$ up
to stage $m$ is $\mathcal{O}(m^{2})$ and the memory complexity is
$\mathcal{O}(m)$. This is infeasible in practice, as for real problems it
hundreds of stages are needed. To overcome this significant problem, we
propose an enhanced version of the original algorithm which reduces the time
complexity to $\mathcal{O}(m)$ and the memory complexity to $\mathcal{O}(1)$.
We present this new, enhanced algorithm (with pseudocode).

###  3.1 Algorithm

In order to reduce the computational complexity of evaluating
$(\bm{\ell}^{-n},\bm{h}^{-n})(s,\bm{X}_{s}^{n})$ in the situation when
$\bm{\ell}^{-n}$ or $\bm{h}^{-n}$ explicitly appears in the minimizer in (14),
we propose the Enhanced Deep Fictitious Play which parametrizes both
$V^{n}(t,{\bm{x}})$ and policy $(\ell^{n},h^{n})(t,{\bm{x}})$ by NNs. For
simplicity, we state the algorithm based on a generic stochastic differential
game, where (possibly high-dimensional) controls are denoted by
$\alpha^{n}(t,{\bm{x}})$ for player $n$.

In each stage of the Enhanced Deep Fictitious Play, for each planner $n$, the
loss that our algorithm aims to minimize consists of two parts: (1) the loss
related to solving (17)–(18) and (2) the error of approximating the optimal
strategy $\alpha^{n}$ within some hypothesis spaces. The resulted
approximation $\tilde{\alpha}^{n}$ will be used in the next stage of
fictitious play:

$\begin{split}&\inf_{Y_{0}^{n},\tilde{\alpha}^{n},\left\\{Z_{t}^{n}\right\\}_{0\leq
t\leq
T}}\mathbb{E}(\left|Y_{T}^{n}\right|^{2}+\tau\int_{0}^{T}\left\lVert\alpha^{n}(s,\bm{X}_{s}^{n})-\tilde{\alpha}^{n}(s,\bm{X}_{s}^{n})\right\rVert_{2}^{2}\,\mathrm{d}s)\\\
\text{ s.t.
}&\bm{X}_{t}^{n}={\bm{x}}_{0}+\int_{0}^{t}\mu^{n}(s,\bm{X}_{s}^{n};\tilde{\bm{\alpha}}^{-n}(s,\bm{X}_{s}^{n}))\,\mathrm{d}s+\int_{0}^{t}\Sigma(\bm{X}_{s}^{n})\,\mathrm{d}\bm{W}_{s},\\\
&Y_{t}^{n}=Y_{0}^{n}-\int_{0}^{t}g^{n}(s,\bm{X}_{s}^{n},Z_{s}^{n};\tilde{\bm{\alpha}}^{-n}(s,\bm{X}_{s}^{n}))\,\mathrm{d}s+\int_{0}^{t}(Z_{s}^{n})^{\operatorname{T}}\,\mathrm{d}\bm{W}_{s},\\\
&{\alpha}^{n}(s,\bm{X}_{s}^{n})=\operatorname*{arg\,min}_{\beta^{n}}H^{n}(s,\bm{X}_{s}^{n},(\beta^{n},\tilde{\bm{\alpha}}^{-n})(s,\bm{X}_{s}^{n}),Z_{s}^{n}),\end{split}$
(20)

where $\|\cdot\|_{2}$ denotes the 2-norm, $\tilde{\bm{\alpha}}^{-n}$ denotes
the collection of approximated optimal controls from the previous stage except
player $n$, and $\tau$ is a hyperparameter denoting the weight between two
terms in the loss function. As detailed in Section 3.2, the hypothesis space
for which we search $\tilde{\alpha}^{n}$ is characterized by another NN, in
addition to the one to approximate $Y_{0}$ and
$\left\\{Z_{t}^{n}\right\\}_{0\leq t\leq T}$. Although representing
$\tilde{\alpha}^{n}$ with a neural network introduces approximation errors, it
allows us to efficiently access the proxy of the optimal strategy
$\bm{\alpha}^{-n}$ in the last stage by calling corresponding networks,
instead of storing and calling all the previous strategies
$\bm{\alpha}^{-n,m-1},\ldots,\bm{\alpha}^{-n,1}$ due to the recursive
dependence.

Numerically we solve a discretized version of (LABEL:def:varitional_problem).
Given a partition $\pi$ of size $N_{T}$ on the time interval $[0,T]$,
$0<=t_{0}<t_{1}<....<t_{N_{T}}=T$, the algorithm reads (to ease the notation,
we replace the subscript $t_{k}$ by $k$):

$\displaystyle\inf_{\tiny\psi_{0}\in\mathcal{N}_{0}^{n^{\prime}},\left\\{\phi_{k}\in\mathcal{N}_{k}^{n},\xi_{k}\in\mathcal{N}_{k}^{n^{\prime\prime}}\right\\}_{k=0}^{N_{T}-1}}\mathbb{E}\\{|Y_{T}^{n,\pi}|^{2}+\tau\sum_{k}\left\lVert\alpha_{k}^{n,\pi}-\tilde{\alpha}_{k}^{n,\pi}(\bm{X}_{k}^{n,\pi})\right\rVert_{2}^{2}\Delta
t_{k}\\}$ (21) s.t. $\displaystyle\bm{X}_{0}^{n,\pi}=\bm{X}_{0},\quad
Y_{0}^{n,\pi}=\psi_{0}\left(\bm{X}_{0}^{n,\pi}\right),\quad
Z_{k}^{n,\pi}=\phi_{k}\left(\bm{X}_{k}^{n,\pi}\right),\quad\tilde{\alpha}^{n,\pi}_{k}(\bm{X}_{k}^{n,\pi})=\xi_{k}(\bm{X}_{k}^{n,\pi}),$
(22)
$\displaystyle{\alpha}_{k}^{n,\pi}=\operatorname*{arg\,min}_{\beta^{n}}H^{n}(t_{k},\bm{X}_{k}^{n,\pi},(\beta^{n},\tilde{\bm{\alpha}}^{-n,\pi}_{k})(\bm{X}_{k}^{n,\pi}),Z_{k}^{n,\pi}),\quad
k=0,\ldots,N_{T}-1$ (23)
$\displaystyle\bm{X}_{k+1}^{n,\pi}=\bm{X}_{k}^{n,\pi}+\mu^{n}\left(t_{k},\bm{X}_{k}^{n,\pi};\tilde{\bm{\alpha}}^{-n,\pi}_{k}(\bm{X}_{k}^{n,\pi})\right)\Delta
t_{k}+\Sigma\left(t_{k},\bm{X}_{k}^{n,\pi}\right)\Delta\bm{W}_{k},$ (24)
$\displaystyle
Y_{k+1}^{n,\pi}=Y_{k}^{n,\pi}-g^{n}\left(t_{k},\bm{X}_{k}^{n,\pi},Z_{k}^{n,\pi};\tilde{\bm{\alpha}}^{n,\pi}_{k}(\bm{X}_{k}^{n,\pi})\right)\Delta
t_{k}+\left(Z_{k}^{n,\pi}\right)^{\mathrm{T}}\Delta\bm{W}_{k},$ (25)

where $\Delta t_{k}=t_{k+1}-t_{k}$,
$\Delta\bm{W}_{k}=\bm{W}_{t_{k+1}}-\bm{W}_{t_{k}}$, and $\mathcal{N}_{0}^{n}$,
$\\{\mathcal{N}_{k}^{n^{\prime}}\\}_{k=0}^{N_{T}-1}$,
$\\{\mathcal{N}_{k}^{n^{\prime\prime}}\\}_{k=0}^{N_{T}-1}$ are hypothesis
spaces for player $n$, which will be specified later through neural network
structures.

The expectation in (21) is further approximated by Monte Carlo samples of
(24)-(25). The parameters in the hypothesis spaces are determined by
stochastic gradient descent (SGD) algorithms such that the approximated
expectation is minimized, which in turn gives the optimal deterministic
functions $(\psi_{0}^{\ast},\phi_{k}^{\ast},\xi_{k}^{\ast})$. We expect that
$(\psi_{0}^{\ast},\phi_{k}^{\ast},\xi_{k}^{\ast})$ will approximate
$(V^{n},\nabla_{\bm{x}}V^{n},\alpha^{n})$ well when this proxy of (21) is
small. Particularly, $\\{\xi^{\ast}_{k}\\}_{k=0}^{N_{T}-1}$ serves as an
efficient tool to evaluate the optimal policy at the current stage for finding
Nash equilibrium. Implementation details and the full algorithm are presented
in Section 3.2. Note that when $\tau=0$ and $\tilde{\alpha}$ are replaced by
$\alpha$ in the above algorithm, the Enhanced Deep Fictitious Play degenerates
to the Deep Fictitious Play proposed in Han and Hu(2020).

###  3.2 Implementation

Here we provide some detail to implement the methodology in Section 3.1.
First, we specify the hypothesis spaces for neural networks
$\mathcal{N}_{0}^{n^{\prime}}$,$\\{\mathcal{N}_{k}^{n}\\}_{k=0}^{N_{T}-1}$,
$\\{\mathcal{N}_{k}^{n^{\prime\prime}}\\}_{k=0}^{N_{T}-1}$, corresponding to
$V^{n},\nabla_{\bm{x}}V^{n},\alpha^{n}$ (the superscript $m$ is dropped again
for simplicity). $V^{n}(t,x)$ is parametrized directly by a neural network
$\mathrm{NN}(t,{\bm{x}})$ . Corresponding map
$\Sigma(\bm{X})^{\operatorname{T}}\nabla_{{\bm{x}}}V^{n}(t,\bm{X})$ that
defines $Z_{t}^{n}$ in the optimization problem (LABEL:def:varitional_problem)
could be parametrized by
$\Sigma({\bm{x}})\nabla_{{\bm{x}}}\mathrm{NN}(t,{\bm{x}})$. Naturally,
$\Sigma({\bm{x}})\nabla_{{\bm{x}}}\mathrm{NN}(t_{k},{\bm{x}})$ is a hypothesis
function in $\mathcal{N}_{k}^{n^{\prime}}$. Under this parametrization rule,
the hypothesis functions in $\mathcal{N}_{0}^{n^{\prime}}$ and
$\\{\mathcal{N}_{k}^{n}\\}_{k=0}^{N_{T}-1}$ share the same set of parameters.
The policy function $\alpha^{n}(t,{\bm{x}})$ is parametrized by another neural
network $\widetilde{\mathrm{NN}}(t,{\bm{x}})$ and then
$\widetilde{\mathrm{NN}}(t_{k},{\bm{x}})$ plays the role of a hypothesis
function in $\mathcal{N}_{k}^{n^{\prime\prime}}$. In other words,
$\\{\mathcal{N}_{k}^{n^{\prime\prime}}\\}_{k=0}^{N_{T}-1}$ share the same set
of neural networks. In a stochastic game of $N$ players, there are $2N$ neural
networks in total, with $N$ neural networks corresponding to
$V^{n}(t,{\bm{x}})$ and $N$ neural networks corresponding to
$\alpha^{n}(t,{\bm{x}})$. At stage $m$, the $N$ $V$-networks are trained to
approximate the solution of PDE (15) and the $N$ $\alpha$-networks are trained
to approximate the current optimal policy computed by (14) using the optimal
strategies in the last stage. The updated neural networks at stage $m$ would
be used at stage $m+1$ to simulate paths $\\{X_{k}^{n,\pi}\\}_{k=0}^{N_{T}-1}$
and optimal strategies by (14). In this work, fully connected neural networks
with three hidden layers are used.

Second, at each stage, the $2N$ neural networks could be decoupled to $N$
pairs of $V$-network and $\alpha$-network based on players. Then, the $N$
pairs of neural networks could be trained in parallel, which dramatically
reduces computational time. As Han and Hu(2020) and Seale and Burnett(2006)
pointed out, it is not necessary to solve the individual control problem
accurately in each stage; the parameters at each stage are updated starting
from the optimal parameters in the last stage without re-initialization. This
requires only a moderate number of epochs for the stochastic gradient descent
at each stage.

The full implementation of Enhanced Deep Fictitious Play is shown in Algorithm
1. For simplicity, we state the algorithm based on a generic stochastic
differential game.

Algorithm 1 Enhanced Deep Fictitious Play for Finding Markovian Nash
Equilibrium

0: $N$ = # of players, $N_{T}$ = # of subintervals on $[0,T]$, $M$ = # of
total stages in fictitious play, $N_{\text{sample}}$ = # of sample paths
generated for each player at each stage of fictitious play,
$N_{\text{SGD\\_per\\_stage}}$ = # of SGD steps for each player at each stage,
$N_{\text{batch}}$ = batch size per SGD update, $\bm{\alpha}^{0}\colon$ the
initial policies that are smooth enough

1: Initialize $N$ deep neural networks to represent $V^{n,0}$ and $N$ deep
neural networks to represent $\alpha^{n,0},n\in\mathcal{N}$

2: for $m\leftarrow 1$ to $M$ do $n\in\mathcal{N}$

3: Generate $N_{\text{sample}}$ sample paths
$\\{\bm{X}_{k}^{n,\pi}\\}_{k=0}^{N_{T}}$ according to (24) and the realized
approximate optimal policies
$\tilde{\bm{\alpha}}^{-n,m-1}(t_{k},\bm{X}_{k}^{n,\pi})$ (Remark:
$\tilde{\bm{\alpha}}$ represents social and health policies in the multi-
region SEIR model)

4: for $e\leftarrow 1$ to $N_{\text{SGD}\\_\text{per}\\_\text{stage}}$ do

5: Update the parameters of the $n^{th}$ $V$-neural network and
$\alpha$-neural network one step with $N_{\text{batch}}$ paths using the SGD
algorithm (or its variant), based on the loss function (21)

6: end for

7: Obtain the approximate optimal policy $\tilde{\alpha}^{n,m}$ represented by
the latest policy neural network

8: end for

9: Collect the approximate optimal policies at stage $m$:
$\tilde{\bm{\alpha}}^{m}\leftarrow(\tilde{\alpha}^{1,m},\ldots,\tilde{\alpha}^{N,m})$

end forThe approximate optimal policy $\tilde{\bm{\alpha}}^{M}$

Due to page limits, the exact choice of NN architectures will be detailed in
Appendix B.1. To determine the total stages of fictitious play $M$, we monitor
the relative changes of $\alpha^{n}$ and $V^{n}$, and stop the process when
the relative change from stage to stage is below a threshold. Regarding the
total number of SGD per stage, as shown in [Han and Hu(2020), Figure 1], the
original DFP is insensitive to the choice of $N_{\text{SGD\\_per\\_stage}}$.
We find the enhanced version sharing the same behavior when apply to the
COVID-19 case study. We give more details in Section 4.2, and further
experiments regarding different choices of $M$ and
$N_{\text{SGD\\_per\\_stage}}$ in Appendix B.2. For problems without
analytical solutions, one natural concern is the reliability of numerical
solutions. Theoretically, the quantity (21) serves as the indicator of the
numerical accuracy. In the original DFP where the second term in (21) does not
exist, Theorem 3 in [Han et al.(2020b)Han, Hu, and Long] ensures the
convergence to the true Nash equilibrium under technical assumptions when (21)
is small enough for each fictitious play stage and with sufficiently large $M$
and small $\Delta t_{k}$. In practice, the quantity in (21) is approximated by
its Monte Carlo counterpart, which we define as the loss function of our
algorithms. Therefore, having small training losses during all stages will
ensure convergence. Extending Theorem 3 in [Han et al.(2020b)Han, Hu, and
Long] to the current setting is beyond the scope of this paper and is left for
further work.

##  4 Application on COVID-19

Our case study is based on COVID-19. We focus mainly on the lockdown/travel
ban policy between different regions. Therefore, to simplify the presentation,
we omit the health policy $h$ in the following discussion and make
$v(\cdot)=v$, $\lambda(\cdot)=\lambda$, and $\eta=0$. Moreover, as vaccines
are not available to the population yet, we let $v(\cdot)=v=0$.

###  4.1 Parameter choices

In single-region SEIR models, the transmission rate, $\beta$, is the basic
reproductive number divided by the length of time an individual is infectious.
In our model, we assume that there is a region-independent constant $\beta$
that underlies the rate of infections for each population. The transmission
rates $\beta^{nk}$ between regions are related to the underlying transmission
rate $\beta$, and the amount of travel between regions $n$ and $k$. To
quantify the size of travel between regions, we assume there is a constant
fraction of people from region $n$ that travel to region $k$, $f^{nk}$, at any
given moment in time. We note that realistically one may expect $f^{nk}$ to
depend on time and also on the epidemic status of regions $n$ and $k$.
However, for simplicity, we will not consider these scenarios in our numerical
experiments. We assume that $f^{nn}\gg f^{nk}$, $f^{nn}\gg f^{kn}$, $\forall
k\neq n$, meaning that most of the population $n$ resides in region $n$ at any
given time, and also that most of the people in region $n$ at a given time are
from $n$ and not travelers from another region. We will see later that this
implies $\beta^{nn}\gg\beta^{nk}$ $\forall k\neq n$. To further clarify the
transmission of infection from one region due to another, we need to describe
the set of parameters $\\{\beta^{nk}:n,k\in\mathcal{N}\\}$. Here, $\beta^{nk}$
represents the rate of transmission from region $k$ to region $n$.
Specifically, $\beta^{nk}$ is the number of infected people in population $n$
per a contactable, infectious individual in population $k$ per day, assuming
that $100\%$ of population $n$ is susceptible. This definition of $\beta^{nk}$
comes from the derivation of the SDE system itself (1)–(3). Accounting for
infection in region $n$ by individuals in region $k$ from travel to both
region $n$ and $k$, we have that

$\beta^{nk}=\begin{dcases}\beta(f^{nk}f^{kk}+f^{kn}f^{nn})\frac{P^{k}}{P^{n}},&\mathrm{if}\
k\neq n\\\ \beta(f^{nn})^{2},&\mathrm{if}\ k=n.\\\ \end{dcases}$ (26)

We defer the detailed derivation of (26) to Appendix A.3. Therefore, to
specify $\beta^{nk}$, we need to provide $\beta$ and $f^{nk}$. We will specify
$f^{nk}$ for New York (NY), New Jersey (NJ), and Pennsylvania (PA) in the next
section. To estimate $\beta$, we choose a basic reproductive number
$R_{0}=2.2$, which is consistent with [Fauci et al.(2020)Fauci, Lane, and
Redfield], and assume that the length of each individual being infectious is
$13$ days. More precisely, we assume infected individuals either recover or
die in 13 days. Under these assumptions, we obtain
$\beta=\frac{2.2}{13}\approx 0.17$, consistent with [Linton et
al.(2020)Linton, Kobayashi, Yang, Hayashi, Akhmetzhanov, Jung, Yuan,
Kinoshita, and Nishiura] as a 13 day median time until death from illness
onset is used. The infection fatality rate, or IFR, is the fraction of those
infected who died from the infection. We choose the IFR to be 0.65% according
to the CDC estimate. This is also consistent with [Meyerowitz-Katz and
Merone(2020)], which suggests a point estimate of 0.68%. The assumptions of an
IFR of 0.65% and an infectious period of 13 days determines that the recovery
rate (including both recovery and death due to infection) is
$\lambda=\frac{1}{13}\approx 0.0769$, and the death rate is
$\kappa=\frac{(0.65\%)}{13}=0.0005$. We choose the latent period to be 5 days
according to [Lauer et al.(2020)Lauer, Grantz, Bi, Jones, Zheng, Meredith,
Azman, Reich, and Lessler]. This means that the we will have
$\gamma=\frac{1}{5}$. Note that this choice has also been used in other models
such as [Alvarez et al.(2020)Alvarez, Argente, and Lippi] and [Peng et
al.(2020)Peng, Yang, Zhang, Zhuge, and Hong]. We assume that the parameters
for noise-level $\sigma_{s_{n}},\sigma_{e_{n}}$, $n\in\mathcal{N}$ are all
$0.0002$, and the extent to which one adheres to the social distancing policy,
$\theta$, is either $\theta=0.9$ or $\theta=0.99$. With most of the parameters
for the SDE model (1)-(3) discussed, we now address those specific to defining
the cost. Regarding the risk-free-rate $r$, note that U.S. Treasury yields are
historically low and the uncertainty in the current level of inflation. We
choose for simplicity that $r=0$. Also, considering that we are interested in
simulations with time periods of less than a year, the discounting is
negligible. The parameter $w$ represents the dollar output per individual per
day. To estimate $w$, we use GDP per capita per day, yielding the estimate
$w=172.6$ dollars per person per day. Following [Alvarez et al.(2020)Alvarez,
Argente, and Lippi] and [Hall et al.(2020)Hall, Jones, and Klenow], we use the
value of a statistical life, $\chi$, to be 20 times GDP per capita. This
results in $\chi=1.95\cdot 10^{6}$ dollars per person. According to the CDC
summary of U.S. COVID-19 activity, the hospitalization rate was 228.7 per
100,000 population by 11/14/2020. Thus, we set $p=228.7\times 10^{-5}$. The
cost per inpatient day is $c=73300/13$ dollars, estimated according to
[Health(2020)]. The attention hyperparameter $a$ takes various values in the
case study, and will be specified in Section 4.2.

###  4.2 NY-NJ-PA COVID-19 case study

In this section, we apply our model (1)–(5) to analyze COVID-19 related policy
in three adjacent states: New York, New Jersey, and Pennsylvania. This case
study is done over 180 days starting from 03/15/2020, using the Enhanced Deep
Fictitious Play algorithm introduced in Section 3 (cf. Algorithm 3.2) and the
parameters discussed in 4.1. The exact formulas of $\mu^{n}$ and $h^{n}$ in
equation (15) are derived in Appendix A.2. We refer to New York State as
region 1, New Jersey State as region 2, and Pennsylvania State as region 3.
Their respective populations are $P^{1}=19.54$ million, $P^{2}=8.91$ million,
and $P^{3}=12.81$ million. Regarding $\beta^{nk}$, $\forall n,k=1,2,3$, we
assume that: (a) 90% of any state’s population is residing in their state at a
given time; (b) the remaining population (travelers) visit the other regions
in an equal proportion; and (c) there is no travel outside of the considered
regions, i.e., the NY-NJ-PA is a closed system. The reasoning for (c) is that,
under our model assuming that infection only occurs in the regions considered,
(c) is equivalent to allowing people traveling outside the considered regions,
but the travelers cannot be affected. For simplicity, we assume this is the
case. Under these assumptions, we will have $f^{nn}=90\%$ for $n=1,2,3$ and
$f^{nk}=5\%$ for $n\neq k$, and obtain the values of $\beta^{nk}$ through
(26). Figure 1 presents the equilibrium policy issued by the governors of NY,
NJ, and PA when the policy effectiveness is $\theta=0.99$, i.e., 99% of the
population follow the lockdown order. The hyperparameter is $a=100$, i.e.,
each governor values people’s death 100 times the lockdown cost. In this
scenario, the governors take action at an early stage and soon reach the
strictest policy. Once the disease is under control, they may relax the policy
later. The percentage of Susceptible, Exposed, Infectious, and Removed stays
almost constant in the end. As a comparison, Figure 2 illustrates how the
pandemic gets out of control if governors show inaction or issue mild lockdown
policies.

Figure 1: Plots of optimal policies (top-left), Susceptibles (top-right),
Exposed (bottom-left) and Infectious (bottom-right) for three states: New York
(blue), New Jersey (orange) and Pennsylvania (green). The shaded areas depict
the mean and 95% confidence interval over 256 sample paths. Choices of
parameters are in Section 4.1, $a=100$ and $\theta=0.99$. Figure 2: An
illustration that governors’ inaction or mild control leads to disease
spreading.

Experiment 1: dependence on $a$. We further analyze how the planners’ view on
the death of human beings changes their policies. In reality, economic loss is
not the only factor the planners concern about. It is also important to
mitigate the infections and deaths within the budget and available resources.
Different views and values from the planners will lead to different policies.
In this experiment, we consider different attitudes towards the infection,
especially death caused by COVID-19. This is reflected by the attention
hyperparameter $a$. Large $a$ implies that planners care more about human
beings and are willing to spend more effort or endure more economic loss on
lockdown to avoid further infection and death. In comparison, smaller $a$
implies that planners care less about infection and death and instead pay more
attention to minimizing the total cost.

Figure 3: Plots of optimal policies with different choice of $a$ for three
states: New York (blue), New Jersey (orange) and Pennsylvania (green), when
the lockdown efficiency is $\theta=0.9$. Figure 4: Plots of optimal policies
with different choice of $a$ for three states: New York (blue), New Jersey
(orange) and Pennsylvania (green), when the lockdown efficiency is
$\theta=0.99$.

The numerical results in Figures 3 and 4 are consistent with intuition. With a
large $a$ (top-left panels), meaning the planners give more consideration to
infection and death, they tend to issue a restrict lockdown policy, which
helps slow down the disease spreads and reduce the percentage of infected
people. As $a$ becomes smaller (top-right panels), planners weigh more the
economic loss and spend fewer efforts on lockdown. When the attention $a$ is
small enough, some states even give up controlling the disease spread due to
economic concern (bottom panels). As a result, the pandemic would get out of
control by the end of the simulation period. This mild lockdown policy leads
to a natural spread of disease (also shown in Figure 2). Experiment 2:
dependence on $\theta$. We next analyze how the residents’ willingness to
comply with the lockdown policy changes the optimal policies and the
development of a pandemic. The larger the $\theta$ is, the more likely the
residents will follow the lockdown policy, and the larger the difference the
control makes on the pandemic situation. Conversely, small $\theta$ weakens
the effect of the lockdown policy. In the extreme case of $\theta=0$, no
matter how strict the lockdown policy is, the pandemic will become a natural
spread because the control term in (1) disappears. In short, this willingness
to policy compliance should be an essential factor in decision-making. To this
end, we compare the optimal policy when $\theta=0.9$ and $\theta=0.99$ in
Figure 5. Panels (a-d) show the difference of optimal policies $\ell(t)$ and
the Susceptible $S(t)$ in the tri-state game under different $\theta$ when
$a=50$. In both situations, the pandemic is well-controlled, with the
percentage of susceptible people staying stable in the end. Moreover, in the
case of $\theta=0.99$, people are more willing to comply with the policies.
Consequently, the planners are allowed to use a less strict lockdown policy as
shown in Figure 5(b) compared to 5(a), which saves the lockdown cost. Figure 5
(e-h) shows an interesting case in the comparison of $\theta=0.9$ and
$\theta=0.99$. In this scenario, with the same attention parameter ($a=25$),
$\theta=0.9$ leads to a mild lockdown policy, see Figure 5(e), while
$\theta=0.99$ provides a possibility to stop the spread of virus, see Figure
5(f). We believe that the decision when $\theta=0.9$ is a compromise as the
lockdown is not efficient enough to reduce largely the infection and death
loss by paying lockdown cost, and also due to the limited simulation period,
i.e., the policies could have been different if we had the simulation until
the disease dies out. We also believe that the early give-up by NJ drives NY
and PA to lift lockdown policies at a later stage, because even NY and PA
issue strict policies, they are still facing severe infections from NJ due to
its high infected percentages and the existence of travel between states
whatever the policy is. So their interventions are not worth the candle.
Figure 5(f)(h) further elucidate the importance of residents’ support in
slowing down the pandemic. Further experiments based on different sets of
$(a,\theta)$ reveal the possibility of having multiple Nash equilibrium, with
more elaboration in Appendix B.3.

Figure 5: Comparison of optimal policies for three states (NY = blue, NJ =
orange, PA = green) and their susceptibles between different policy
effectivenss $\theta$ and hyperparameter $a$.

To summarize, the numerical experiments illustrate that both the balance of
economy and infection/death from the view of plan-makers and the willingness
of residents to follow the lockdown policy play an important role in decision-
making. In reality, all three states issued stay-at-home orders in March, and
attempted to reopen in June. By comparing real world policies and our
simulations of $\bm{\ell}(x)$, we may infer $\alpha$ and $\theta$ for NY, NJ,
and PA in our model, i.e., $\theta=0.99$ and $a=25$.

##  5 Conclusion

In this paper, we propose a novel multi-region SEIR model to study optimal
policies under a pandemic. Our new model, built on game theory, takes into
account how the social and health policies issued by multiple region planners
affect the progress of infectious diseases. This feature makes the model more
realistic and powerful but also introduces a formidable computational
challenge due to the high-dimensionality of the solution space and the strong
coupling of planners’ policies. We propose the enhanced deep fictitious play
algorithm to overcome the curse of dimensionality and use the model and
algorithm in a case study of the COVID-19 pandemic in three states, New York,
New Jersey, and Pennsylvania. The model parameters are estimated from real
data posted by the CDC. We are able to show the effect of lockdown/travel ban
policy on the spread of COVID-19 for each state and how people’s willingness
to comply and planners’ attitude towards deaths influence the equilibrium
strategies as a consequence of the competition between regions. We hope our
model can draw more attention to studying optimal interventions in infectious
diseases using game theory. Our numerical simulations can shed light on public
policies. In reality, during a pandemic, the planning is usually for short
periods and adjusted frequently. This can be modeled using repeated games, and
planners may infer other regions’ cost functional from past game outcomes. The
assumptions that some parameters are identical across different regions can be
relaxed and the health policy can also be added for more accurate simulations.
These will be left for future work. R.H. was partially supported by NSF grant
DMS-1953035. H.D.C. acknowledges partial support from NSF grant DMS-1818821.
This project was jointly supervised by R.H. and H.D.C.. Y.X. and R.B. have
equal contribution as the first authors. J.H. contributes on the algorithms
and early discussions of the model setup.

## References

  * [Alvarez et al.(2020)Alvarez, Argente, and Lippi] F. E. Alvarez, D. Argente, and F. Lippi.  A simple planning problem for COVID-19 lockdown.  Working Paper 26981, National Bureau of Economic Research, 2020.

  * [Bauch and Earn(2004)] C. T. Bauch and D. J. D. Earn.  Vaccination and the theory of games.  _Proceedings of the National Academy of Sciences_ , 101(36):13391–13394, 2004.
  * [Bauch et al.(2003)Bauch, Galvani, and Earn] C. T. Bauch, A. P. Galvani, and D. J. D. Earn.  Group interest versus self-interest in smallpox vaccination policy.  _Proceedings of the National Academy of Sciences_ , 100(18):10564–10567, 2003.
  * [Brown(1949)] G. W. Brown.  Some notes on computation of games solutions.  Technical report, Rand Corp Santa Monica CA, 1949.
  * [Brown(1951)] G. W. Brown.  Iterative solution of games by fictitious play.  _Activity Analysis of Production and Allocation_ , 13(1):374–376, 1951.
  * [Chang et al.(2020)Chang, Piraveenan, Pattison, and Prokopenko] S. L. Chang, M. Piraveenan, P. Pattison, and M. Prokopenko.  Game theoretic modelling of infectious disease dynamics and intervention methods: a review.  _Journal of Biological Dynamics_ , 14(1):1–33, 2020.
  * [E et al.(2017)E, Han, and Jentzen] W. E, J. Han, and A. Jentzen.  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations.  _Communications in Mathematics and Statistics_ , 5(4):349–380, 2017.
  * [El Karoui et al.(1997)El Karoui, Peng, and Quenez] N. El Karoui, S. Peng, and M. C. Quenez.  Backward stochastic differential equations in finance.  _Mathematical Finance_ , 7(1):1–71, 1997.
  * [Fauci et al.(2020)Fauci, Lane, and Redfield] A. S. Fauci, H. C. Lane, and R. R. Redfield.  COVID-19 —navigating the uncharted.  _New England Journal of Medicine_ , 382(13):1268–1269, 2020.
  * [Hall et al.(2020)Hall, Jones, and Klenow] R. E. Hall, C. I Jones, and P. J. Klenow.  Trading off consumption and COVID-19 deaths.  Working Paper 27340, National Bureau of Economic Research, 2020.
  * [Han and Hu(2020)] J. Han and R. Hu.  Deep fictitious play for finding Markovian Nash equilibrium in multi-agent games.  In _Proceedings of The First Mathematical and Scientific Machine Learning Conference (MSML)_ , volume 107, pages 221–245, 2020.
  * [Han et al.(2018)Han, Jentzen, and E] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, 2018.
  * [Han et al.(2020a)Han, Hu, and Long] J. Han, R. Hu, and J. Long.  Convergence of deep fictitious play for stochastic differential games.  _arXiv preprint arXiv:2008.05519_ , 2020a.
  * [Han et al.(2020b)Han, Hu, and Long] J. Han, R. Hu, and J. Long.  Barron metric for the convergence of empirical distribution.  _in preparation_ , 2020b.
  * [Health(2020)] Fair Health.  Costs for a Hospital Stay for COVID-19, 2020.  URL https://www.fairhealth.org/article/costs-for-a-hospital-stay-for-covid-19.
  * [Hu(2020)] R. Hu.  Deep fictitious play for stochastic differential games.  _Communications in Mathematical Sciences_ , 2020.
  * [Isaacs(1965)] R. Isaacs.  _Differential Games: A Mathematical Theory with Applications to Warfare and Pursuit, Control and Optimization_.  London: John Wiley and Sons, 1965.
  * [Lauer et al.(2020)Lauer, Grantz, Bi, Jones, Zheng, Meredith, Azman, Reich, and Lessler] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R. Meredith, A. S. Azman, N. G. Reich, and J. Lessler.  The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application.  _Annals of Internal Medicine_ , 172(9):577–582, 2020.
  * [Linton et al.(2020)Linton, Kobayashi, Yang, Hayashi, Akhmetzhanov, Jung, Yuan, Kinoshita, and Nishiura] N.M. Linton, T. Kobayashi, Y. Yang, K. Hayashi, A.R. Akhmetzhanov, S.-M. Jung, B. Yuan, R. Kinoshita, and H. Nishiura.  Incubation period and other epidemiological characteristics of 2019 novel coronavirus infections with right truncation: A statistical analysis of publicly available case data.  _Journal of Clinical Medicine_ , 538(9), 2020.
  * [Liu et al.(1987)Liu, Hethcote, and Levin] W.-M. Liu, H. W. Hethcote, and S. A. Levin.  Dynamical behavior of epidemiological models with nonlinear incidence rates.  _Journal of Mathematical Biology_ , 25(4):359–380, 1987.
  * [Meyerowitz-Katz and Merone(2020)] G. Meyerowitz-Katz and L. Merone.  A systematic review and meta-analysis of published research data on COVID-19 infection-fatality rates.  _medRxiv_ , 2020.  10.1101/2020.05.03.20089854.
  * [Pardoux and Peng(1992)] E. Pardoux and S. Peng.  Backward stochastic differential equations and quasilinear parabolic partial differential equations.  In _Stochastic Partial Differential Equations and Their Applications_ , pages 200–217. Springer, 1992.
  * [Pardoux and Tang(1999)] E. Pardoux and S. Tang.  Forward-backward stochastic differential equations and quasilinear parabolic PDEs.  _Probability Theory and Related Fields_ , 114(2):123–150, 1999.
  * [Peng et al.(2020)Peng, Yang, Zhang, Zhuge, and Hong] L. Peng, W. Yang, D. Zhang, C. Zhuge, and L. Hong.  Epidemic analysis of covid-19 in china by dynamical modeling, 2020.
  * [Seale and Burnett(2006)] D. Seale and J. Burnett.  Solving large games with simulated fictitious play.  _International Game Theory Review_ , 8(03):437–467, 2006.

##  Appendix A Technical Details to the Stochastic Multi-region SEIR Model

###  A.1 The dynamics of $\bm{X}_{t}$ in Section 2

In Section 2, for the ease of notations and clarity of the presentation, we
rewrite the dynamics of $(S_{t}^{n},E_{t}^{n},I_{t}^{n})$ defined in (1)–(3)
in the vector form

$\,\mathrm{d}\bm{X}_{t}=b(t,\bm{X}_{t},\bm{\ell}(t,\bm{X}_{t}),\bm{h}(t,\bm{X}_{t}))\,\mathrm{d}t+\Sigma(\bm{X}_{t})\,\mathrm{d}\bm{W}_{t},$
(27)

where
$\bm{X}_{t}\equiv[\bm{S}_{t},\bm{E}_{t},\bm{I}_{t}]^{\operatorname{T}}\equiv[S_{t}^{1},\cdots,S_{t}^{N},E_{t}^{1},\cdots,E_{t}^{N},I_{t}^{1},\cdots,I_{t}^{N}]^{\operatorname{T}}\in\mathbb{R}^{3N}$,
the Markovian controls $(\bm{\ell},\bm{h})$ are given by
$\bm{\ell}(t,{\bm{x}})=[\ell^{1},\ldots,\ell^{N}]^{\operatorname{T}}(t,{\bm{x}})$
and $\bm{h}(t,{\bm{x}})=[h^{1},\ldots,h^{N}]^{\operatorname{T}}(t,{\bm{x}})$.
In the sequel, for a vector
${\bm{x}}\equiv(\bm{s},\bm{e},\bm{i})\in\mathbb{R}^{3N}$, we shall index them
in two ways,

$(s_{1},\cdots,s_{N},e_{1},\cdots,e_{N},i_{1},\cdots,i_{N})\text{ or
}(x_{1},\ldots,x_{3N}).$ (28)

and use them interchangeably. The former one emphasizes the dependence on each
category, while the later notation is more condensed. Similarly, for partial
derivatives, we will have two set of notations

$(\partial_{s_{1}},\cdots,\partial_{s_{N}},\partial_{e_{1}},\cdots,\partial_{e_{N}},\partial_{i_{1}},\cdots,\partial_{i_{N}})\text{
or }(\partial_{x_{1}},\cdots,\partial_{x_{3N}}).$ (29)

We give precise definitions of (27) in this appendix.
$b(t,{\bm{x}},\bm{\ell},\bm{h})=[b_{1},\ldots,b_{3N}]^{\operatorname{T}}(t,{\bm{x}},\bm{\ell},\bm{h})$
is a deterministic vector-valued function:

$b_{j}(t,{\bm{x}},\bm{\ell},\bm{h})=\begin{dcases}-\sum_{k=1}^{N}\beta^{jk}s_{j}i_{k}(1-\theta\ell^{j}(t,{\bm{x}}))(1-\theta\ell^{k}(t,{\bm{x}}))-v(h^{j}(t,{\bm{x}}))s_{j},\qquad
j\in\mathcal{N},\\\
\sum_{k=1}^{N}\beta^{{j^{\prime}}k}s_{j^{\prime}}i_{k}(1-\theta\ell^{j^{\prime}}(t,{\bm{x}}))(1-\theta\ell^{k}(t,{\bm{x}}))-\gamma
e_{j^{\prime}},\qquad\quad j\in\mathcal{N}+N,\\\ \gamma
e_{j^{\prime}}-\lambda(h^{j^{\prime}}(t,{\bm{x}}))i_{j^{\prime}},\qquad\qquad
j\in\mathcal{N}+2N,\text{ and }{j^{\prime}}=j~{}\text{mod}~{}N.\end{dcases}$
(30)

$\Sigma({\bm{x}})=(\Sigma_{j,k}({\bm{x}}))$ is a matrix-valued deterministic
function in $\mathbb{R}^{3N\times 2N}$ with non-zero entries given below:

$\displaystyle\Sigma_{j,j}({\bm{x}})=-\sigma_{s_{j}}s_{j},\quad$
$\displaystyle\Sigma_{j+N,j}({\bm{x}})=\sigma_{s_{j}}s_{j},$ (31)
$\displaystyle\Sigma_{j+N,j+N}({\bm{x}})=-\sigma_{e_{j}}e_{j},\quad$
$\displaystyle\Sigma_{j+2N,j+N}({\bm{x}})=\sigma_{e_{j}}e_{j},\quad
j\in\mathcal{N}.$ (32)

and $\\{\bm{W}_{t}\\}_{0\leq t\leq T}$ is a $2N$-dimensional standard Brownian
motion:

$\bm{W}_{t}=[W^{s_{1}}_{t},\cdots,W^{s_{N}}_{t},W^{e_{1}}_{t},\cdots,W^{e_{N}}_{t}]^{\operatorname{T}}.$

Each region’s running cost $f^{n}$ defined in (9) is

$f^{n}(t,{\bm{x}},\bm{\ell},\bm{h})=e^{-rt}P^{n}[(s_{n}+e_{n}+i_{n})\ell^{n}(t,{\bm{x}})w+a(\kappa
i_{n}\chi+pi_{n}c)]+e^{-rt}\eta(h^{n}(t,{\bm{x}}))^{2}.$ (33)

###  A.2 The decoupled HJB equations in the form of (15)

Recall that we aim to solve (13) using the BSDE approach (nonlinear Feynman
Kac relation). To this end, in this appendix, we will rewrite it in the form
of (15) and identify $\mu^{n}$ and $g^{n}$. The first step is to identify the
minimizer in the Hamiltonian (12). Keeping in mind that our testing case is
COVID-19 where vaccines are not fully developed yet, the term $v(h_{t}^{n})=0$
is essentially zero in the past. Also, to focus on the lockdown/travel ban
policy between different regions (as we did in the case study), we shall
exclude the health policy $\bm{h}(t,{\bm{x}})$ from planners’ decision
problem, i.e., $v(\cdot)=0$, $\lambda(\cdot)=\lambda$, and $\eta=0$ in the
following derivations, and remove the dependence of $\bm{h}$ in all relevant
functions. We remark that including the health policy $\bm{h}(t,{\bm{x}})$ is
a straightforward generalization. Recall that the Hamiltonian in (13) reads:

$\displaystyle
H^{n}(t,{\bm{x}},(\ell^{n},\bm{\ell}^{-n,m}),\nabla_{\bm{x}}V^{n,m+1})$ (34)
$\displaystyle\qquad=b(t,{\bm{x}},(\ell^{n},\bm{\ell}^{-n,m}))\cdot\nabla_{\bm{x}}V^{n,m+1}+f^{n}(t,{\bm{x}},\ell^{n})$
(35)
$\displaystyle\qquad=\sum_{j=1}^{3N}b_{j}(t,{\bm{x}},(\ell^{n},\bm{\ell}^{-n,m}))\partialderivative{V^{n,m+1}}{x_{j}}+e^{-rt}P^{n}[(s_{n}+e_{n}+i_{n})\ell^{n}(t,{\bm{x}})w+a(\kappa
i_{n}\chi+pi_{n}c)],$ (36)

and recall that
$\bm{\ell}^{-n,m}=(\ell^{1,m},\ldots,\ell^{n-1,m},\ell^{n+1,m},\ldots,\ell^{N,m})$
represents the $m^{th}$ stage strategies of all players other than $n$, which
are given functions in this derivation. The first order condition requires for
$\ell^{n}$:

$\displaystyle 0=\sum_{\begin{subarray}{c}j=1\\\ j\neq
n\end{subarray}}^{N}(1-\theta\ell^{j,m})\left[{\beta^{{j}{n}}s_{j}i_{n}\left(\partialderivative{V^{{n},{m+1}}}{e_{j}}-\partialderivative{V^{{n},{m+1}}}{s_{j}}\right)}+{\beta^{{n}{j}}s_{n}i_{j}\left(\partialderivative{V^{{n},{m+1}}}{e_{n}}-\partialderivative{V^{{n},{m+1}}}{s_{n}}\right)}\right]$
(37)
$\displaystyle+2(1-\theta\ell^{n}){\beta^{{n}{n}}s_{n}i_{n}\left(\partialderivative{V^{{n},{m+1}}}{e_{n}}-\partialderivative{V^{{n},{m+1}}}{s_{n}}\right)}-\frac{1}{\theta}e^{-rt}P^{n}(s_{n}+e_{n}+i_{n})w.$
(38)

The critical point given by the above equation indeed gives a minimizer of the
Hamiltonian, as long as it is in $[0,1]$. Because we can show
$\left(\partialderivative{V^{n,m+1}}{e_{n}}-\partialderivative{V^{n,m+1}}{s_{n}}\right)>0$
by comparing $V^{n,m+1}(t,{\bm{x}}+\epsilon_{n+N})$ and
$V^{n,m+1}(t,{\bm{x}}+\epsilon_{n})$ using their definitions, where
$\epsilon_{j}$ is a $3N$-vector with only one nonzero entry $\epsilon\ll 1$ at
$j^{th}$ position. Intuitively, with all others players’ initial condition the
same, if player $n$ starts with a higher exposed proportion $e_{n}+\epsilon$,
it will produce more cost, comparing with the same increase proportion still
being susceptible $s_{n}+\epsilon$. To summarize, we deduce the optimal policy
for player $n$ at stage $m+1$ is given by:

$\displaystyle\ell^{n,m+1}(t,{\bm{x}})$
$\displaystyle=\Bigg{\\{}2{\beta^{{n}{n}}s_{n}i_{n}\left(\partialderivative{V^{{n},{m+1}}}{e_{n}}-\partialderivative{V^{{n},{m+1}}}{s_{n}}\right)}-\frac{1}{\theta}e^{-rt}P^{n}(s_{n}+e_{n}+i_{n})w$
(39) $\displaystyle+\sum_{\begin{subarray}{c}j=1\\\ j\neq n\end{subarray}}(1-$
$\displaystyle\theta\ell^{j,m})\left[{\beta^{{j}{n}}s_{j}i_{n}\left(\partialderivative{V^{{n},{m+1}}}{e_{j}}-\partialderivative{V^{{n},{m+1}}}{s_{j}}\right)}+{\beta^{{n}{j}}s_{n}i_{j}\left(\partialderivative{V^{{n},{m+1}}}{e_{n}}-\partialderivative{V^{{n},{m+1}}}{s_{n}}\right)}\right]\Bigg{\\}}$
(40)
$\displaystyle\times\left\\{2\theta{\beta^{{n}{n}}s_{n}i_{n}\left(\partialderivative{V^{{n},{m+1}}}{e_{n}}-\partialderivative{V^{{n},{m+1}}}{s_{n}}\right)}\right\\}^{-1}\wedge
1\vee 0,$ (41)

where we use the conventional notations $a\wedge b=\min\\{a,b\\}$ and $a\vee
b=\max\\{a,b\\}$. Plugging (39) into equation (13) and by straightforward
computation, one obtains for the $(m+1)^{th}$ stage,
$\mu^{n,m+1}(t,{\bm{x}};\bm{\ell}^{-n,m})=[\mu_{1}^{n,m+1},\ldots,\mu_{3N}^{n,m+1}](t,{\bm{x}};\bm{\ell}^{-n,m})^{\operatorname{T}}$
is

$\displaystyle\mu_{j}^{n,m+1}=-\beta^{jn}s_{j}i_{n}(1-\theta\ell^{j,m}(t,{\bm{x}}))-\sum_{\begin{subarray}{c}k=1\\\
k\neq
n\end{subarray}}^{N}\beta^{jk}s_{j}i_{k}(1-\theta\ell^{j,m}(t,{\bm{x}}))(1-\theta\ell^{k,m}(t,{\bm{x}})),$
(42) $\displaystyle\hskip 200.0ptj\in\mathcal{N}\setminus n$ (43)
$\displaystyle\mu_{n}^{n,m+1}=-\beta^{nn}s_{n}i_{n}-\sum_{\begin{subarray}{c}k=1\\\
k\neq n\end{subarray}}\beta^{nk}s_{n}i_{k}(1-\theta\ell^{k,m}(t,{\bm{x}}))$
(44) $\displaystyle\mu_{N+j}^{n,m+1}=-\mu_{j}^{n,m+1}-\gamma
e_{j},\quad\mu_{2N+j}^{n,m+1}=\gamma e_{j}-\lambda i_{j},\quad
j\in\mathcal{N}.$ (45)

To write $g^{n,m+1}$ as a function of $(t,{\bm{x}},z)$, we first compute

$\displaystyle\Sigma({\bm{x}})^{\operatorname{T}}\nabla_{\bm{x}}V^{n}(t,{\bm{x}})$
(46)
$\displaystyle=\left[\sigma_{s_{1}}s_{1}\big{(}\partialderivative{V^{n}}{e_{1}}-\partialderivative{V^{n}}{s_{1}}\big{)},\cdots,\sigma_{s_{N}}s_{N}\big{(}\partialderivative{V^{n}}{e_{N}}-\partialderivative{V^{n}}{s_{N}}\big{)},\sigma_{e_{1}}e_{1}\big{(}\partialderivative{V^{n}}{i_{1}}-\partialderivative{V^{n}}{e_{1}}\big{)},\cdots,\sigma_{e_{N}}e_{N}\big{(}\partialderivative{V^{n}}{i_{N}}-\partialderivative{V^{n}}{e_{N}}\big{)}\right]^{\operatorname{T}}.$
(47)

and then $g^{n,m+1}$ is given by:

$\displaystyle
g^{n,m+1}(t,{\bm{x}},z;\bm{\ell}^{-n,m})=\frac{\theta^{2}}{\sigma_{s_{n}}}\beta^{nn}z_{n}i_{n}[\ell^{n,m+1}(t,{\bm{x}})]^{2}$
(48)
$\displaystyle+\left\\{e^{-rt}P^{n}(s_{n}+e_{n}+i_{n})w-2\frac{\theta}{\sigma_{s_{n}}}\beta^{nn}z_{n}i_{n}-\sum_{\begin{subarray}{c}j=1\\\
j\neq
n\end{subarray}}^{N}\theta(1-\theta\ell^{j,m}(t,{\bm{x}}))(\frac{\beta^{nj}}{\sigma_{s_{n}}}z_{n}i_{j}+\frac{\beta^{jn}}{\sigma_{s_{j}}}z_{j}i_{n})\right\\}\ell^{n,m+1}(t,{\bm{x}})$
(49) $\displaystyle+e^{-rt}P^{n}a(\kappa i_{n}\chi+pi_{n}c).$ (50)

###  A.3 The derivation of the transmission rate $\beta^{nk}$ in (26)

We denote by $\beta$ the underlying transmission rate of the virus, which is
assumed to be region independent. This transmission rate is the average number
of people infected by an infectious person per day (assuming that the
susceptible population is 100%). Thus, if a proportion $S$ of the population
is susceptible, then $\beta S$ represents the average number of people
infected per infectious person per day. If there are a total of $P$ people in
a population with a fraction $I$ being infectious, then $\beta S(PI)$ is the
number of newly infected per day. Thus, in the context of a single-region SEIR
model, the number of newly infected (or the influx to the exposed population)
that occurs within $(t,t+\,\mathrm{d}t)$ is given by $\beta
S(t)(I(t)P)\,\mathrm{d}t$. Dividing by $P$, the influx to the scaled exposed
population in $(t,t+\,\mathrm{d}t)$ is $\beta S(t)I(t)\,\mathrm{d}t$. Now let
us consider the multi-region case and temporarily ignore the effect of
lockdown. The term from (1)–(3) that gives the influx to $E^{n}$ due to
infection from $I^{k}$ is $\beta^{nk}S^{n}(t)I^{k}(t)\,\mathrm{d}t$. To
determine $\beta^{nk}$, we build this exact influx from core assumptions.
First, we quantify the number of people from region $n$ that are infected by
those from region $k$. Specifically, the influx in the interval
$(t,t+\,\mathrm{d}t)$ to the unscaled exposed population $n$ due to
transmission from population $k$ is given by

$\sum_{\ell}(\beta f^{n\ell}S^{n}(t))(f^{k\ell}I^{k}(t)P^{k})\,\mathrm{d}t,$
(52)

where $f^{ij}$ is the (assumed constant) fraction of people from $i$ currently
in region $j$ at any moment in time. Equation (52) is obtained by summing the
number of infections in population $n$ due to population $k$ across each
region. The summand represents these infections occurring in region $\ell$.
This can be seen as the term $f^{n\ell}S^{n}(t)$ is the proportion of
population $n$ that are susceptible and within region $\ell$. Of this
population, there will be $\beta f^{n\ell}S^{n}(t)$ infections per infectious
individual per day. Since the number of infectious from $k$ that are in $\ell$
is $f^{k\ell}I^{k}(t)P^{k}$, we have that $(\beta
f^{n\ell}S^{n}(t))(f^{k\ell}I^{k}(t)P^{k})\,\mathrm{d}t$ is the number of new
infections in population $n$ due to population $k$ occurring in the region
$\ell$ within the time interval $(t,t+\,\mathrm{d}t)$. Let us assume for now
that $n\neq k$. Since we assume that $1>f^{nn}\gg f^{nk}$, the terms in (52)
besides the cases where $\ell=n$ or $\ell=k$ are negligible. Removing the
negligible terms and dividing by the population $P^{n}$, we see that the
influx to the scaled exposed population $E^{n}$ due to transmission from
population $k$ over the interval $(t,t+\,\mathrm{d}t)$ is

$\frac{P^{k}}{P^{n}}\beta(f^{nn}f^{kn}+f^{nk}f^{kk})S^{n}(t)I^{k}(t)\,\mathrm{d}t,$
(53)

which is exactly the influx represented by the model of
$\beta^{nk}S^{n}(t)I^{k}(t)\,\mathrm{d}t$. This verifies the form of
$\beta^{nk}$ for $n\neq k$ in (26). Similarly if we take $n=k$ in (52) and
ignore all terms other than $\ell=n(=k)$, then we derive the formula for
$\beta^{nn}$ in (26).

##  Appendix B Detailed discussions on the enhanced DFP algorithm

###  B.1 Implementation details

In the NY-NJ-PA case study, we choose feedforward architectures for both
$V$-networks and $\alpha$-networks. Both have three hidden layers with a width
of 40 neurons. The activation function in each hidden layer is $\tanh(x)$. We
do not apply activation function to the output layer of $V$-networks, and
choose sigmoid function $\rho_{s}(x)=\frac{1}{1+e^{-x}}$ for the
$\alpha$-networks. Other hyperparameters are summarized in the table below.

hyperparameter | $lr$ | $M$ | $N_{\text{SGD\\_per\\_stage}}$ | $N_{\text{batch}}$ | $N_{T}$ | $\tau$
---|---|---|---|---|---|---
value | 5e-4 | 250 | 100 | 256 | 40 |  $1e^{-3}$/180
Table 1: Hyperparameters in the case study: $lr$ denotes the learning rate in
stochastic gradient descent method, $M$ is the total stages of fictitious
play, $N_{\text{SGD\\_per\\_stage}}$ is the number of stochastic gradient
descent done in each minimization problem (21), $N_{\text{batch}}$ is batch
size in each stochastic gradient descent, $N_{T}$ is the discretization steps
on $[0,T]$, and $\tau$ is the weight of the control part in the loss function
(21).

###  B.2 Discussion on the choice of $M$ and $N_{\text{SGD\\_per\\_stage}}$

We provide further experiments here on various choices of $M$ and
$N_{\text{SGD\\_per\\_stage}}$. In Figure 6, we plot both validation loss and
log loss against $M$ for three states, which are produced by evaluating the
NNs using unseen data after each fictitious play stage. In each panel, loss
curves associated with different number of SGDs per stage are presented in
different colors (blue = 50, orange = 75, green = 100, red = 125, purple =
150). The numerical results show that the validation losses for all states
decrease as the number of DFP stages $M$ increases. Moreover, it shows that
different $N_{\text{SGD\\_per\\_stage}}$ generates loss curves with similar
patterns. Smaller $N_{\text{SGD\\_per\\_stage}}$ is more stable on the
validation loss of PA. This result is consistent with [Han and Hu(2020)] and
[Seale and Burnett(2006)], which convey that it is unnecessary to solve the
problem extremely accurate in each stage and that a moderate number of
$N_{\text{SGD\\_per\\_stage}}$ is sufficient. As a result, we choose
$N_{\text{SGD\\_per\\_stage}}=100$ in our case study.

Figure 6: Loss curves of each state. Left: validation losses versus rounds $M$
of the enhanced deep fictitious play; Right: $\log_{10}$ validation loss
versus rounds $M$ of the enhanced deep fictitious play. The loss curves with
respect to $N_{\text{SGD\\_per\\_stage}}=50,75,100,125,150$ are depicted in
blue, orange, green, purple and red. A smoothed moving average with window
size 3 is applied in the final plots.

###  B.3 Stability over different experiments

Here we present experiments to investigate the Nash equilibrium of the model
with different combinations of parameters. For each combination of parameters,
we use the same hyper-parameters and repeat the experiments several times. We
run the algorithm for a certain computational budget, and then filter out the
results with a fluctuating loss near the stopping and check the converged
equilibrium. In the first combination of parameters, we take $\theta=0.99$ and
$a=100$, corresponding to the case that a governor weighs the deaths much more
than the economic loss and tries to avoid it, and the residents have a strong
willingness to follow the governor’s policies. Intuitively, the pandemic is
possible to get well-controlled. Our numerical experiments confirm this
intuition: all converging trails lead to the same Nash equilibrium. A
representative plot of $X(t)=(S(t),E(t),I(t))$ is shown in Figure 7.

Figure 7: With the parameter combination $\theta=0.99$, $a=100$, the algorithm
identifies one Nash equilibrium for the NY-NJ-PA case study.

In the second batch of experiments, we take $\theta=0.9$ and $a=50$,
corresponding to the case that a governor pays less attention to the number of
death and the residents are less willing to follow the policies compared to
the first batch of experiments. The change leads to the possibility of
multiple Nash equilibrium and the pandemic being out of control. In this case,
with different NNs’ initialization, the algorithm identifies two Nash
equilibria: $75\%$ of the experiments converge to the Nash equilibrium that
the pandemic gets controlled and $25\%$ of the experiments converge to the
other Nash equilibrium where the pandemic gets out of control.

Figure 8: With the parameter combination $\theta=0.9$, $a=50$, the algorithm
identifies two possible Nash equilibria: an under control one (topic panels,
with 75% of the experiments) and on out-of-control one (bottom panels, with
25% of the experiments).

In conclusion, it is possible to have multiple Nash equilibria depending on
the parameter chosen in our stochastic multi-region SEIR model. There is
usually a single Nash equilibrium for parameters chosen at extreme values,
while for the parameters selected in the middle range, there could exist
multiple Nash equilibria. When multiple equilibria exist, we conjecture that
the possibility to reach a particular one depends on where we start the
fictitious play (the initialization of the NNs’ parameters).


* Daddi-Moussa-Ider et al. (2018) A. Daddi-Moussa-Ider, B. Rallabandi, S. Gekle, and H. A. Stone, Phys. Rev. Fluids 3, 084101 (2018).
  * Kim and Karrila (2005) S. Kim and S. J. Karrila, _Microhydrodynamics: Principles and Selected Applications_ , Butterworth - Heinemann series in chemical engineering (Dover Publications, 2005), ISBN 9780486442198\.
  * de Hoog and Lekkerkerker (2001) E. H. A. de Hoog and H. N. W. Lekkerkerker, J. Phys. Chem. B 105, 11636 (2001).
  * Song and Waugh (1990) J. Song and R. E. Waugh, J. Biomech. Eng. 112, 235 (1990), ISSN 0148-0731.
  * Song and Waugh (1993) J. Song and R. Waugh, Biophys. Journal 64, 1967 (1993), ISSN 0006-3495.
  * Style et al. (2017) R. W. Style, A. Jagota, C.-Y. Hui, and E. R. Dufresne, Annu. Rev. Cond. Mat. Phys. 8, 99 (2017).
  * Andreotti and Snoeijer (2020) B. Andreotti and J. H. Snoeijer, Annu. Rev. Flu. Mech. 52, 285 (2020).
  * Koch and Subramanian (2011) D. L. Koch and G. Subramanian, Annu. Rev. Flu. Mech. 43, 637 (2011).
  * Lauga (2016) E. Lauga, Annu. Rev. Flu. Mech. 48, 105 (2016).
  * Berg (2018) H. C. Berg, _Random walks in biology_ (Princeton University Press, 2018).
  * Takatori et al. (2014) S. C. Takatori, W. Yan, and J. F. Brady, Phys. Rev. Lett. 113, 028103 (2014).
  * Dunstan et al. (2012) J. Dunstan, G. Miño, E. Clement, and R. Soto, Phys. Fluids 24, 011901 (2012).
  * Hinch (1991) E. J. Hinch, _Perturbation Methods_ , Cambridge Texts in Applied Mathematics (Cambridge University Press, 1991).
  * Stenhammar et al. (2017) J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo, and A. Morozov, Phys. Rev. Lett. 119, 028005 (2017).
  * Strogatz (2015) S. H. Strogatz, _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (2nd ed.)_ (CRC press, 2015).
  * Guasto et al. (2012) J. S. Guasto, R. Rusconi, and R. Stocker, Annu. Rev. Flu. Mech. 44, 373 (2012).
  * Finnegan et al. (2015) S. Finnegan, S. C. Anderson, P. G. Harnik, C. Simpson, D. P. Tittensor, J. E. Byrnes, Z. V. Finkel, D. R. Lindberg, L. H. Liow, R. Lockwood, et al., Science 348, 567 (2015).
  * Simpson (2021) C. Simpson, The American Naturalist 198, 590 (2021).
  * Marcos et al. (2012) Marcos, H. C. Fu, T. R. Powers, and R. Stocker, Proc. Nat. Acad. Sci. U.S.A. 109, 4780 (2012), ISSN 0027-8424.
  * Mesibov and Adler (1972) R. Mesibov and J. Adler, J. Bacteriology 112, 315 (1972).
  * Pedley et al. (1988) T. J. Pedley, N. A. Hill, and J. O. Kessler, J. Fluid Mech. 195, 223–237 (1988).
  * Vincent and Hill (1996) R. V. Vincent and N. A. Hill, J. Fluid Mech. 327, 343–371 (1996).
  * Cisneros et al. (2010) L. H. Cisneros, R. Cortez, C. Dombrowski, R. E. Goldstein, and J. O. Kessler, _Fluid dynamics of self-propelled microorganisms, from individuals to concentrated populations_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010), pp. 99–115, ISBN 978-3-642-11633-9.
  * Painter and Hillen (2011) K. J. Painter and T. Hillen, Physica D: Nonlinear Phenomena 240, 363 (2011), ISSN 0167-2789.
  * Subramanian et al. (2011) G. Subramanian, D. L. Koch, and S. R. Fitzgibbon, Phys. Fluids 23, 041901 (2011).
  * Ezhilan et al. (2012) B. Ezhilan, A. A. Pahlavan, and D. Saintillan, Phys. Fluids 24, 091701 (2012).
  * Desai and Ardekani (2018) N. Desai and A. M. Ardekani, Phys. Rev. E 98, 012419 (2018), URL https://link.aps.org/doi/10.1103/PhysRevE.98.012419.
  * Saintillan (2018) D. Saintillan, Annu. Rev. Flu. Mech. 50, 563 (2018).
  * Bender and Orszag (2013) C. M. Bender and S. A. Orszag, _Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory_ (Springer New York, 2013), ISBN 9781475730692\.
  * Abramowitz et al. (1988) M. Abramowitz, I. A. Stegun, and R. H. Romer, _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , Applied mathematics series (U.S. Government Printing Office, 1988).