key: cord-0876378-ncnaast9
authors: Petrakova, Viktoriya; Krivorotko, Olga
title: Mean Field Game for Modeling of COVID-19 Spread()
date: 2022-04-19
journal: J Math Anal Appl
DOI: 10.1016/j.jmaa.2022.126271
sha: c75d0fdfbd1387855ccb708da3d01f49ea576dcb
doc_id: 876378
cord_uid: ncnaast9

The paper presents one of the possible approaches to pandemic spread modelling. The proposed model is based on the mean-field control inside separate groups of population, namely, suspectable (S), infected (I), removed (R) and cross-immune (C) ones. The numerical algorithm to solve this problem ensures conservation of the total population mass during timeline. The numerical experiments demonstrate modelling results for COVID-19 spread in Novosibirsk (Russia) for two 100-day periods.

The COVID-19 pandemic has not only caused hundreds of thousands of deaths around the world, but also had a huge impact on the global economy, and everyone's daily life. As a rule, SIR-type compartmental differential models are used to describe the dynamics of transmission of infectious diseases in a 5 population. However, a SIR model stops working if it is necessary to account for population heterogeneity or some stochastic phenomena that are significant Mean Field Game for Modeling of COVID- 19 Spread for small populations at the initial phase of disease spread. Since the spread of COVID-19 has a significant spatial characteristic, various restrictive measures have been used to slow down disease propagation, exerting an additional in- 10 fluence on the dynamics of population behavior (by "spatial" we mean certain local features characteristic of a particular population, determined, for example, by location, introduced policies, mass-media coverage, etc.). Thus, to study an infectious-disease spread, it is also necessary to take into account the stochastic parameters of the system and its heterogeneity, which, in the classical sense, 15 leads to systems with a large number of differential equations, even if it is possible to cluster the population somehow. Computationally, the approach to describing population dynamics can be simplified by using the mean-field game (MFG) theory, which enables one to describe the behavioral dynamics of a system with a small number of partial differential equations, taking into account 20 population heterogeneity.

The paper presents a SIRC model for COVID-19 spread, based on a system of ODEs using the MFG approach for which optimality conditions for mean-field system have been obtained (see 2). In Section 3, the optimal control problem with external influence has been solved. In Section 4, a modified computational 25 scheme for the MFG system has been proposed and the results of numerical experiments for COVID-19 spread in the Novosibirsk Region have been discussed.

There are a lot of mathematical models of infectious disease spread that are based on the mass balance law and described by systems of nonlinear or-30 dinary differential equations [1] . The basic compartmental model to describe transmission of infectious diseases was proposed by Kermack and McKendrick in 1927 [2] where the population was divided into three groups: susceptible, 

In this subsection, we will focus on the macroscopic description of the spread without taking into account the individual effect of an agent on the whole population. For these purposes epidemiological models with different characteristics 45 are generally used. Most of these models are based on the susceptible-infectedremoved (SIR) separation. In Casagrandiet et al. [3] , a SIRC model was presented for describing the dynamic behavior of influenza A by introducing a new component, namely, cross-immunity (C) of a population for people who has recovered from infection with different strains of the same virus subtype at pre-50 vious time intervals. In [4, 5, 6] this idea was developed and it was shown in [7] that this model can be used to account for an infected but asymptomatic part of the population. Now, let us consider the dynamics of an epidemic following a simplified SIRC model. Assume that the epidemic occurs in a short time period if compared to 55 the population dynamics (birth and death), so the last one can be neglected.

The flow diagram of the simplified SIRC model is presented in fig. 1 . From the point of modeling, the population at any time moment is divided into four compartments: susceptible individuals S(t), i.e., those who do not have specific immune defenses against the infection; infected individuals I(t); and two com-60 partments of those who are totally or partially immune i.e., recovered R(t) and

cross-immune C(t).

Formally, in accordance with the diagram shown in Fig. 1 and the mass balance law, the SIRC model can be represented by the following system of four 

with initial conditions:

Note, that S(t)+I(t)+R(t)+C(t) = 1 represents the total population. Here the following parameters are used: μ is the rate at which the cross-immune population becomes susceptible again; β is the contact/infection transmission rate; ε is the average reinfection probability of a cross-immune 70 individual; γ is the recovery rate of the infected population; δ is the rate at which the recovered population becomes cross-immune and moves from total to partial immunity. Parameter μ is connected with the average number of days neutralizing antibodies disappear, and δ is the estimated number of days when CD4+ and CD8+ T-cells decrease [8] . Note that in the absence of cross-75 immunity, i.e. (1 − ε) = 0, the SIRC model becomes similar to the SIRS model for the fractions S and C become immunologically indistinguishable [6] .

In fact, greater consistency with some real-world phenomena can be achieved if stochastic changes in a system are taken into account. SIR-based epidemic models assume that observed dynamics are determined by the deterministic cases assuming the population is homogenic and neglecting the characteristics of individual agents. At the same time, accounting for the rationality of an individual who is a part of the population, i.e., for their ability to be in a state X(t) and be able to change it, leads to management problems due to a large number of participants in the system. Taking into account the stochasticity of the process, we assume that the dynamics of a single individual is subject to theÎto differential equation

where i ∈ 1, .., N ; W N i is an independent standart Wiener processes; α i (t) is the i-th agent's strategy, and θ N (t) is the empirical measure of distribution of agents (individuals, players) in the system at time t [9]. Let function b, σ be continuous in time and same for all players. In [9] it was shown that when the number of agents in system rises extremely (N → ∞), the mass of single individuals can be substituted by a representative agent whose state is determined by the following control equation:

Here,

is agents distribution over state space Ω at time t, and α(t) is the control process (or, in other words, strategy the representative agent's strategy) to ensure the Nash equilibrium of a system of interacting agents and minimize the function:

where f and G are Lipschitz functions. This approach to controlling a population with a large number of interacting agents is known as Mean Field Game 80 (MFG). For the more detailed description, refer to original works [9, 10, 11].

For our goals, a population of "atomized" agents is characterized by their states (or positions) x ∈ Ω at each time moment t ∈ [0, T ]. The term "atomized" means that each agent of an infinite set has no influence on the situation (because of its zero-measured support) and chooses a rational strategy α(t, x) accounting for the agent's own position m(t, x) : [0, T ] × Ω → R relative to other agents.

In [10] , it was shown when σ in (2),(3) characterized by the stochastic nature of agents interaction equilibrium is constant, the distribution of agents m(t, x)

obeys the Kolmogorov (Fokker-Planck) equation

with initial condition

and the Neumann boundary condition ∂m/∂x = 0 ∀t and x ∈ Γ Ω .

Here, the boundary condition (7) prevents the loss of density m(t, x) with time.

So, on the one hand, we have a differential SIRC model (1) that describes a population at a macroscopic level, divide it into several parts, but provides no individual description. On the other hand, the individuality can be introduced 85 and controlled using the MFG model (4)- (7) . Following the ideas proposed in [12] , the spatial SIRC model was combined with the MFG one and a man- to the following PDE system:

Here σ i , i ∈ {S, I, R, C} are non-negative parameters to, as before, characterize the stochastic processes within the population. If the initial value of m i is given:

and, similar to (7), it meets zero Neumann boundary conditions that no mass can flow in or out of Ω, one obtains:

Note that

It is noteworthy that summing up the equations in system (8) , one obtains the equation (5) with respect to the total population mass m(t, x), which leads to the following requirement:

i.e., the total mass of the four groups will be conserved ∀t. 90 Now, following the MFG principles, an assumption that agents are rational and tend to choose a strategy to maximize their own benefit in accordance with the cost functions (4) produces the management problem that minimizes the

where the index SIRC denotes the set of all possible values from {S, I, R, C}.

Here F i function is the running cost to implement strategy α i ; g i is payment

for current position or state and 1 0 m I (T, x)dx is the terminal cost of epidemic spread.

As the function of strategy implementation cost F i (ᾱ i , t, x) (an agent's payment for state changing) a piecewise continuous function was considered:

for all admissible values ofᾱ i ∈ R and ∀i ∈ {S, I, R, C}. Note that function has the same form for any part of the population. This can help to account for infected, but asymptotic population masses. If an individual is an infection carrier but does not know about it, in fact, they continue their everyday life for it is more beneficial for them. The piecewise continuous form of (12) means that the transition to self-isolation for an agent (whenᾱ i ≤ 0) is more expensive than maintaining close contact with other people. It is noteworthy that for (12) the following properties satisfy

These properties ensure the unique solvability of the equation

For the g i (t, x, m i ) function:

Here c 1i , c 2i ∈ [0, 1] are positive constants. Thus, the function g i together with Thus, the MFG approach leads to the optimization problem that requires 105 you to find the optimal strategy α i by minimizing the function (11) with the restriction in the form of PDE system (8) with initial (9) and boundary (10) conditions.

To ensure the strategy chosen by a person is optimal the Lagrange multiplier method [10] was used. For that the first equation in (8) 

Do the same with other equations in (8) 

where i ∈ {I, R, C}:

Note that expressions (15)-(18) are valid when the following boundary condi-

and

Now, the Lagrange function corresponding to the optimization problem under consideration can be written down as:

As the result, the minimization of (11) together with the (8)-(10) can be represented [10] as the problem of finding a saddle point:

The variation (15)- (18) with respect to the components

to find a stationary point produces a system analogous to the Hamilton-Jacobi-

where δ(T − t) is the delta function and zero initial condition:

The variation (15)- (18) with respect to the components

gives the following optimality conditions forᾱ i ∈ R in addition to system (23), Thus, two coupled PDE systems with initial and boundary conditions (8) - (10) and (19), (23), (24) together with (25) give the necessary conditions for the minimization of (11).

In our case, agents are considered as "rational" but "selfish" i.e., inclined to choose a strategy to maximize their own benefits based on other agents' strategies, which can lead to a faster spread of the infection. In reality, however, the situation is controlled by some external influence such as introduction of strict quarantine and other restrictions or allocation of social benefits, etc. These quarantine regulations not only determine agents' behavior but, as a rule, are strict requirements and therefore cannot be accounted for in the cost functional where their non-fulfillment becomes possible. In this way, an assumption can be made that the total control of a system's behavior can be determined by two processes:

Hereα i , as before, is a person's strategy and ρ i (α i , t, x) is the external (corrective) control which adjusts the player's strategy depending on the fulfillment of a certain condition. For example, we can put the following expression:

where c 3i ∈ [1, 2] are positive constants and (27) 

A series of numerical experiments was performed to demonstrate the differ-120 ence between the standard differential SIRC model (1) and its MFG extension.

To search for a numerical implementation of the differential problem mentioned above, discrete uniform grids were introduced in in time 

To satisfy the boundary condition (10), it was put that

To construct the computational scheme for equations (8) , the idea firstly proposed by Prof. V.V.Shaydurov [14] was modified for the numerical solution of FPK equation (5)and applied to various MFG statements in [15, 16, 17] written in collaboration with V. Petrakova (formerly Kornienko), one of the authors of this article. The idea was to split the approximation of the FPK equation (5) into two parts: advection and diffusion ones. Note, that the advection part can be written in (8) as:

For diffusion part the same property can be used in (8):

where 

The following notation for the approximation of the differential expression σ 2 i Δm i /2 ∀i ∈ {S, I, R, C} was intoduced:

For approximation of (33), the following finite-difference scheme was applied:

Note that 

where

and 

and initial conditions corresponding to (9)

Thus, instead of differential equations (8) ∀i ∈ { §, I, R, C} the following systems of algebraic equations were obtained:

Hereinafter, the subscript points indicate acceptable values in a corresponding position of the mesh function. A and B i k indicate the following matrices: A and B i k are the matrices of the form is provided by the monotonously property of the M-matrix [18] . Put that ∃ i 0 , k 0 , j 0 for which f i0 k0−1,j0 is negative and m i0 k0,j0 change its value to neg-130 ative one. From the non-negativity of m i0 (40) it follows that f i0 k0−1,j0 is non-negative for non-negative parameters β, μ, ε, δ. The obtained contradiction proves the proposition. 

Proof. Use the key property (45) of the approximation scheme (38). Put

is the solution of (38) for all non-negative f i,h k−1,j+1/2 . The componentsm i,h k,j+1/2 will be non-negative due to M-matrices properties [18] . Multiply (38) by τ and h and sum up over j = 0, ..., N −1 for non-negative f i,h k−1,j+1/2 .

Taking into account (45), the following equality is obtained: ·) into the solution of system (38), where f i,h k−1,j+1/2 can be of any sign. Note that

which is due to the monotone property of the M-matrix [18] that occurs after 

Having taken the maximum of both parts of (46) leads to the required estimate.

Replace the integral cost function in (11) by the discrete one:

Here r h k,i+1/2 is carried out for F (α, t, x) in the following way:

with α i,h k,j := α i,h (t k , x j ) , so the differential minimization problem is replaced by the grid one:

140 ence problem (49) approximates the initial differential formulation by the order

Also, introduce the following:

for any vectors a = a j+1/2 j=0,...,N −1 and b = b j+1/2 j=0,. ..,N −1 . To formulate a discrete optimal control problem, introduce the grid-function set: 

Now write down the Lagrangian for the grid optimization problem (49)

Here (B i k ) * = (B i k ) T means the matrix conjugates to B i k . Then the problem of finding a saddle point (22) in the grid case is written as:

Differentiate the Lagrangian (52) with respect to the individual components to obtain the following system of algebraic equations:

Here

and 

for ∈ {S, R, C} and for i = I:

Proof. Put |ṽ i,h (t k , x j+1/2 )| as the component reaching its maximum absolute value on the layer t k so that |ṽ i,h (t k , x j+1/2 )| = v h,i (t k , ·) ∞,h . Use again the key property of the coefficients (45) to obtain the inequality

For i ∈ {S, R, C} with zero "initial" conditions for (54) the required estimate can be obtained after taking a maximum over k in (58). For i = I from (57)

It is known that for an M-matrix with strict diagonal dominance the following statement holds (Theorem 2 from [19] ):

where R is the amount of diagonal dominance, which is the same for each row of the matrix A and equaled to 1/τ . Then for i = I (58) can be rewritten as

Taking a maximum over k in (59) one obtains the required estimate for i = I.

Also (53) gives the following grid analogue for optimality conditions (29) at the point (t k , x j ) :

Thus, the solution to the discrete optimization problem (49) can be found 145 by the following iterative algorithm:

1. Put the initial value of grid functions α i,h ·,· equal to zero.

2. For the zero functions α i,h ·,· obtain the initial value of functions m i,h ·,· using (38)-(42) and J h (m i,h ·,· , α i,h ·,· ) using (47), (48).

3. Calculate the value of v i,h ·,· functions solving the system (54)-(57). 6. Obtain a new value for the functions by solving J h (m i,h ·,· , α i,h ·,· ) by (47), (48).

·,· , α i,h ·,· ) reaches its minima with the given accu-155 racy then choose functions m i,h ·,· , α i,h ·,· obtained on the previous steps as a solution of (49). Otherwise go to step 3 for new iteration.

For the discrete statement of the optimal control problem with external influence on agent's strategy, the schemes (38)-(43) and (54)-(57) introducing into consideration the following representation can be used: 

wherer h k,i+1/2 is carried out for F (α, t, x) in the following way:

Then, the grid optimization problem can be rewritten as:

The solution of (64) can be found by the iterative algorithm described in the previous section taking into account the representation (61). 

For numerical implementation of the algorithm the parameters presented below were applied. Table 1 describes the epidemiological constants obtained by solving the inverse problem for a statistical data set describing COVID-19

incidence in the Novosibirsk region [20] . 165 For the MFG model the parameters presented in Table 2 were applied. Note that the combination of constants c 1i and c 2i is quite close in its physical meaning, the so-called "index of self-isolation" introduced by the Yandex company The initial distributions of agents mass were written as:

where A i is proportion of current fraction in relation to the total population at the initial time moment; B i is normalization factor equaled to integral over [0, 1] from expression in brackets;

ary conditions (7) for m 0i . Note that A i , ∀i ∈ {S, I, R, C} coincides with the corresponding initial value S(0), I(0), R(0), C(0) for the differential SIRC x c S = 0.8; σ c S = 0.1;

x c I = 0.2; σ c I = 0.1;

x c R = 0.7; σ c R = 0.2;

The physical meaning of constants (66) is that non-infected (S) people do not seek to comply with restrictions and self-isolation for economic reasons while For the A i constants, the following parameters presented in Table 3 were 180 chosen for different time periods.

Here parameters A I are also a solution for the inverse problem. The figure shows the statistical data has two areas describing incidence rise 190 and fall. Since SIR-type models are bad for long-term prediction, the two curves with different epidemiological parameters were put together to approximate the epidemiological situation. The model's epidemiological parameters were restored based on the statistical data and the optimization method presented in the work [21] for two intervals: from May 1 to June 30 and from June 30 to 195 August 8, 2020. The obtained parameters are presented in Table 1 . These parameters were used for numerical implementation of the two models in question.

The spatial and stochastic constants for the MFG implementation are given in Section 4.4. As can be seen from the figure, MFG gives a more accurate approximation to statistical data, but like the parent SIRC model, it does not work 200 well for long-term forecasts. Figure 2 also shows that introducing spatial distribution into consideration has a significant impact on population behavior. Figure 3 demonstrates the dependence of the obtained solution on the selected initial spatial distribution.

The comparison was made on based on a coronavirus report for the city of 205 Novosibirsk from December 1, 2020 when a decline in the incidence was observed. For this time period, the restored epidemiological parameters are given in Table 1 while the stochastic parameter of the system σ 2 was chosen equal to 0.5, and the following curves were used as initial distributions for different population groups:

(1) -initial distributions m 0i determined by (65) with parameters (66);

(2) -initial distributions m 0i determined by (65) The figure 4 shows that taking into account external restrictive measures that are not the agent's choice has a significant impact on population dynamics.

The article describes a way to apply the well-known economic mean field models for forecasting a pandemic spread of epidemics, in particular that of COVID-19. This approach has been opted because traditional epidemiological SIR-based models do not take into account population heterogeneity and there-230 fore cannot be used for long-term forecasts. Another well-known approach to a virus spread, the so-called agent-based models, allow taking into account nonepidemiological factors, but lead to computationally complex systems [22] . In turn, MFG models' structural simplicity and the ability to account for spatial characteristics can solve both problems.

Numerical experiments into modeling the dynamics of COVID-19 spread in Novosibirsk, Russia have shown that despite the fact that the dynamics of the population is determined by the differential SIRC model, taking into account the spatial agent characteristics has a huge impact on the final result. This what gives this model an advantage over SIR-type models but also generates a 

This computational approach is applicable for any SIR-type model.

The numerical algorithm considered in this paper is based on the ideas proposed in the works [14, 15, 16, 17] and modified for application to epidemiology problems. It gives a direct and simple rule for minimizing the cost functional, ensures the fulfillment of the law of conservation of the entire agents mass, and 250 allows to take into account more complex control functions F i (α, t, x) instead of those that are quadratically dependent on α [10] .

This work was supported by the Russian Science Foundation (project No. 18-71-10044).

Mathematical epidemiology: Past, present, and future

A contribution to the mathematical the-260 ory of epidemics

The SIRC model and influenza A

Dynamical Behavior of a Stochastic SIRC Model for Influenza A

The Fractional SIRC Model and Influenza A, Mathematical Problems in Engineering

Global Dynamics of a Nonautonomous SIRC Model for Influenza A with Distributed Time Delay

Dynamics of a SIRC epidemiological model, Elec-275 tronic

Immunological memory to SARS-CoV-2 assessed for up to 8 months after infection

On the connection between symmetric N-player games and mean field games

Mean Field Games and Mean Field Type Control Theory

Mean field games

Controlling Propagation of epidemics via mean-field control

Economic. Principles, Problems, and Policies

Conservative Difference Schemes for the Computation of Mean-field Equilibria, AIP Conf. Proc. 1895

The Euler-Lagrange Approximation of the Mean Field Game for the Planning Problem

A finite-difference solution of Mean Field problem with a predefined control resource

Numerical Solution of Mean Field Problem with Limited Management Resource

M-matrix characterizations. I-nonsingular M-305 matrices

Norm estimates for the inverses of matrices of monotone type and totally positive matrices

Official information on the situation with coronavirus in the Novosibirsk region

Mathematical modeling and forecasting of COVID-19 in Moscow and Novosibirsk region

Agentbased modeling of COVID-19 outbreaks for New York state and UK: parameter identification algorithm