Modeling COVID-19 with Mean Field Evolutionary Dynamics: Social Distancing and Seasonality

I. INTRODUCTION

THE third zoonotic human coronavirus of the century, which is known as SARS-CoV-2, emerged at the end of 2019, having a wide-ranging and severe impact upon many aspects of our society, especially health, the economy, and social life [1]. According to the World Bank [2], the average annual decline in GDP of the world in 2020 was -3:595%, compared with the increase of 2:334% in 2019. Moreover, COVID- 19, the human disease caused by SARS-CoV-2, poses a very serious threat to human lives worldwide. As of July 23rd, 2021, 192,489,618 people had been infected and 4,134,432 people had died due to the COVID-19 outbreak [3], and these numbers have increased significantly since then as the Delta variant of the virus has taken hold in many parts of the world. To mitigate the impact on the economy and human lives, the development of an accurate dynamical model for the prediction of the evolution of COVID-19 within the population is of vital importance.

A. Epidemic Models

A classical model for epidemic spread is the SIS (acronym: susceptible (S), infected (I), and susceptible (S)) epidemic model [4], [5]. In this model, infected individuals are assumed to lose immunity and become susceptible immediately after recovering. Impulsive differential equations are used to analyze the complex dynamics of SIS epidemic models. Another classical model is the Kermack-McKendrick epidemic model, which is a special case of a more general model introduced in [6], which divides the population into three main categories: susceptible individuals, infected individuals, and recovered individuals (SIR). The evolution of the population distribution through these three states is described by replicator dynamics (RDs) with simple assumptions on the rates of flow between different classes of individuals. The basic SIR epidemic model has been extended into many distinct epidemic models by creating new individual states. For example, [7] has extended the basic SIR model to the SIVR model, where “V" denotes a “variation" state. The variation state can characterize the mutations of influenza viruses during their spreading process. Reference [8] considered a continuous-time epidemic model with sub-populations of susceptible-exposed-infectedrecovered (SEIR) under a general feedback vaccination control rule. Reference [9] further extended the SEIR model into the discrete SEIADR (acronym: susceptible (S), exposed (E), symptomatic infectious (I), asymptomatic infectious (A), dead infectious (D), and recovered (R)) epidemic model, by incorporating the asymptomatic infectious and the dying infective bodies as infectious extra populations within the standard populations of SEIR type models. However, people who are asymptomatic infectious are usually not disclosed at all or just recorded as infectious in statistical summaries. In the paper, to model the evolution of COVID-19 with data collected in Wuhan, China, we will use the SIDR model, extended from the SIR model by adding the “D” state to represent the deaths.

B. Epidemic Dynamics

The replicator dynamics, introduced in the mathematical biology literature by [10], are the primary dynamics used to describe epidemic evolution under distinct epidemic models. Developed based on imitation, RDs require a bare minimum of information: each agent only needs to know its payoff for the current state/strategy. However, RDs fail Nash stationarity (which refers to the condition in which the rest points always coincide with the Nash Equilibrium) [11], which means that the equilibrium cannot be obtained by simply finding the rest points.

Other evolutionary dynamics also have the potential to describe epidemic evolution, although they have not been explored in this context. The best response dynamics [12] satisfy Nash stationarity. However, the protocol that generates these dynamics are discontinuous, requiring knowledge of the payoffs to all available strategies to obtain the current best response. The Brown-von Neumann-Nash (BNN) dynamics [13] satisfy Nash stationarity, and are established on continuous revision protocols. But these protocols also require that agents know the average payoff to the population, which should be provided by a central source. Finally, pairwise comparison dynamics, which first appeared in the transportation science literature [14], and were later developed by [11], satisfy Nash stationarity while only making limited informational demand: payoffs to the current strategies and randomly chosen candidate strategies. The novel epidemic dynamics that we propose in this paper fall into the category of pairwise comparison dynamics.

Aside from the prediction of the spread of the pandemic, epidemic dynamics can also be useful in many other fields, such as computer science and the study of social networks. In computer science, the propagation of computer viruses can be studied with stochastic epidemic dynamics [15], and based on such analyses, practical measures can be taken to contain the spread of viruses [16]. Similarly, the popularity assessment of web content on social media with epidemic dynamics is useful in applications such as caching, recommendation, advertising, and prediction of economic trends.

C. Mean Field Games in Epidemic Modeling

Mean field games, introduced in [17], are designed to deal with control problems involving a large number of rational agents with limited information, such as the spread of epidemics. Reference [18] introduced a mean field game model for controlling the propagation of epidemics on a spatial domain. The spread of epidemics within a spatial domain is formulated as an optimal control problem with the target of the number of infectious agents and the amount of movement of the population. Reference [19] proposed a mean field game model where each agent selects a dynamic strategy of making contacts, considering the trade-off of gaining utility but also risking infection from additional contacts. However, the above optimal-control-based epidemic models need to design models that require specification of an appropriate utility function to govern the evolution of epidemics. This cost function is typically hard to specify because we do not know what is regarded as “utility” by the agents. For example, in terms of COVID-19, some people think freedom should be more important while others think health should be more important. With mean field evolutionary dynamics (MFEDs), we bridge mean field games and evolutionary games. Moreover, we derive payoff functions for different individual states from the typically used RDs to govern the evolution of epidemics. In this way, we circumvent the definition of a global utility function.

D. Contribution

The main contributions of this paper are summarized as follows:

We propose the MFEDs as a model for epidemic dynamics. Moreover, we compare it with traditional RDs numerically with simulation results in the SIR, SIVR, and SIDR models.

We design the SIDR model and achieve a good fit with the statistics of COVID-19 in Wuhan, China. MFED modeling outperforms RD modeling significantly in terms of coefficient of determination and bias of prediction.

We show the effect of one-time social distancing and the seasonality of COVID -19 through the post-pandemic period.

The remainder of this paper is organized as follows. In Section II, we propose RDs in SIR, SIVR, and SIDR epidemic models, and derive the corresponding payoff functions. In Section III, we propose the MFEDs and compare them with RD theoretically. In Section IV, we compare the convergence behaviors of MFEDs and RDs with numerical results. In Section V, we model the evolution of public COVID-19 data in Wuhan, China with MFEDs. Section VI draws our conclusion.

II. REPLICATOR DYNAMICS

In Section II-A, preliminaries of population games on graphs are given at first. In Section II-B, we analyze the general form of RDs. In Section II-C, we propose the RDs in the SIR, SIVR, and SIDR models. Finally, in Section II-D, we derive the payoff functions in the SIR, SIVR, and SIDR models.

A. Preliminaries of Population Game on Graph

In order to give the general form of RDs (also the general form of MFEDs), we need to clarify the following fundamental concepts. We consider a population game on graph [TeX:] $$\mathcal{G}=(\mathcal{S}, \mathcal{E}):$$

Nodes and Edges: Nodes of this graph are pure strategies from the discrete strategy set [TeX:] $$\mathcal{S}=\{1,2, \cdots, s\}.$$ Edges are connections between nodes. Node [TeX:] $$i \in \mathcal{S} \text { and node } j \in \mathcal{S}$$ are able to form an edge [TeX:] $$(i, j) \in \mathcal{E}$$ if players can directly switch from strategy i to strategy j.

Neighborhood: The neighborhood of node i is the set of all nodes which have a direct connection to node i. It is defined as follows:

Population State Space: The population state space consists of all available population distribution on the discrete strategy set S and it is defined in the following way:

where [TeX:] $$\rho_{i}$$ represents the fraction of population selecting strategy i. The interior of [TeX:] $$\mathcal{P}(\mathcal{S})$$ is denoted as [TeX:] $$\mathcal{P}_{o}(\mathcal{S}).$$

Payoff Function: The payoff function to strategy i, [TeX:] $$F_{i}: \mathcal{P}(\mathcal{S}) \rightarrow \mathbb{R}$$ is a mapping from the current population distribution to the reward of selecting strategy i. Each agent is only interested in maximizing its own reward by selecting different strategies.

Remark 1. In the epidemic model, each individual state (susceptible, infected, and recovered) is a pure strategy for the player, i.e., a node on the strategy graph. The transition between different individual states is the edge of the strategy graph. The evolution of the epidemic is equivalent to the evolution of the population in the population state space.

B. Replicator Dynamics in General Form

The RDs, introduced in the mathematical biology literature by [10], are the most thoroughly studied evolutionary dynamics and have been widely regarded as the epidemic dynamics in many epidemic models. Their general expression is given by

where [TeX:] $$\rho_{i}$$ represents the fraction of the population selecting strategy [TeX:] $$i, \rho=(\rho)_{i=1}^{s}$$ is the population distribution on all strategies, [TeX:] $$F_{i}: \mathcal{P}(\mathcal{S}) \rightarrow \mathbb{R}$$ is the payoff to strategy i, and [TeX:] $$\bar{F}(\rho)=\sum_{i=1}^{s} \rho_{i} F_{i}(\rho)$$ is the average payoff of the whole population. Under these dynamics, the growth rate of the population on each strategy i is equivalent to its excessive payoff, i.e., to the difference between its payoff and the average payoff of the population. Intuitively, if the payoff of a given strategy is higher than the average, it will be selected by more agents as every individual in the population is trying to maximize his/her payoff.

C. Replicator Dynamics in SIR, SIVR, and SIDR

Denoting the fraction of susceptible population as [TeX:] $$\rho_{s},$$ the fraction of the population infected by the original pathogen as [TeX:] $$\rho_{I},$$ and the fraction of the population recovered as [TeX:] $$\rho_{r},$$ RDs [6] which describe the evolution of the epidemic in the SIR model are

These equations are developed based on the following assumptions: (i) The total population size is a constant N and there is no entry into or departure from the population. (ii) One infected individual can make enough contact with [TeX:] $$\beta N$$ others to transmit infection per unit time, where [TeX:] $$\beta$$ is the infection rate. (iii) Infected individuals recover from the disease and get permanent immunity at a rate [TeX:] $$\alpha \rho_{I}$$ per unit time. (iv) There are no deaths due to the disease. And (v) pathogens do not evolve in response to changing environments and medical interventions, i.e., evolutionary adaptations [20] are not considered. The strategy graph for the SIR model is shown in Fig. 1(a).

In the SIVR epidemic model, the mutations of the pathogen are considered during the spreading process, and thus the strategy graph contains the extra “variation (V)" state as shown in Fig. 1(b). We denote the fraction of population that is infected by the mutated pathogen as [TeX:] $$\rho_{v}.$$ Then the RDs [7] that describe the evolution of the epidemic under the SIVR epidemic model are

where [TeX:] $$\sigma_{1} \text { and } \sigma_{2}$$ are the infection rates of the two pathogens respectively, [TeX:] $$\mu_{1} \text { and } \mu_{2}$$ are the recovery rates of the two diseases, respectively. The population size is also assumed to be a constant and infected individuals will get permanent immunity once recovered.

In order to depict the transmission of COVID-19, we need an individual state to represent the deaths. Therefore, we extend the SIR model into the SIDR model, where the state “D” refers to the fraction of the population who die from COVID-19. The strategy graph for the SIDR epidemic model has been shown in Fig. 1(c). The RDs for the SIDR model are

where [TeX:] $$\rho_{d}$$ represents the fraction of the population who die from COVID-19, [TeX:] $$\eta, \epsilon_{1}, \epsilon_{2}$$ are the infection rate, recovery rate, and death rate of COVID-19, respectively. The population size is assumed to be a constant. Infected individuals will either die of recover with permanent immunity.

D. Payoff Functions in SIR, SIVR, and SIDR

We now derive the payoff functions for each state/strategy under the SIR and SIVR epidemic models. Comparing (4) with the general form in (3), we obtain the following equations:

where [TeX:] $$F_{s}, F_{I}, \text { and } F_{r}$$ are the payoffs when an individual becomes susceptible, infected, and recovered, respectively. [TeX:] $$\rho= \left[\rho_{s}, \rho_{I}, \rho_{r}\right]$$ is a vector recording the population distribution. [TeX:] $$\bar{F}(\rho)$$ is the average payoff of the whole population. Solving the equation systems in (7), we obtain the payoff functions for the SIR model as follows:

Likewise, by fitting (5) into (3) and solving the corresponding equation system, we can obtain the payoff functions for the SIVR model as follows:

Fig. 1.

Strategy graphs for distinct epidemic models. Acronym: Susceptible (S), infectious (I), recovered (R), variation (V), death (D)

where [TeX:] $$\hat{F}_{s}, \hat{F}_{I}, \hat{F}_{v}, \hat{F}_{r},$$ are the payoffs when an individual becomes susceptible, infected by the original virus, infected by the mutated virus, and recovered, respectively.

The payoff functions for the SIDR model can be found in a similar way. By transforming the RDs in (6) into the general form in (3), we can obtain the payoff functions for the SIDR model as follows:

where [TeX:] $$\tilde{F}_{s}, \tilde{F}_{I}, \tilde{F}_{d}, \tilde{F}_{r},$$ are the payoffs when an individual becomes susceptible, infected, dead, and recovered, respectively.

In this section, we first give the general form of RDs and show the intuitive ideas behind it. Then we show three specific RDs in the SIR, SIVR, and SIDR epidemic models. Finally, we derive the payoff functions for different individual states in the SIR, SIVR, and SIDR epidemic models. In the next section, these derived payoff functions will be used to govern the evolution of epidemics based on MFEDs.

III. MEAN FIELD EVOLUTIONARY DYNAMICS

In this Section, we establish the MFED-based epidemic models. Inserting the derived payoff functions in (8), (9), and (10), MFEDs can describe how the agents on different states interact with each other and how the population evolves between different states of COVID-19. The remainder of this section is organized as follows: In Section III-A, we introduce the general form of MFEDs; in Section III-B, we insert corresponding payoff functions and derive the MFEDs in the SIR, SIVR, and SIDR models.

A. Mean Field Evolutionary Dynamics in General Form

MFEDs are evolutionary dynamics for population games with discrete strategy sets. They were originally inspired by optimal transport theory [21]–[23] and mean field games [24]–[26]. A general form of MFEDs is given in the following theorem.

Theorem 1. Suppose that a population game has strategy graph [TeX:] $$\mathcal{G}=(\mathcal{S}, \mathcal{E}),$$ the payoff functions [TeX:] $$F_{i}: \mathcal{P}(\mathcal{S}) \rightarrow \mathbb{R}$$ are continuous. Then for any initial condition [TeX:] $$\rho^{0} \in \mathcal{P}_{o}(\mathcal{S}),$$ the Fokker-Planck equation

describes evolutionary dynamics in [TeX:] $$\mathcal{P}_{o}(\mathcal{S}) \cdot[\cdot]^{+}=\max \{\cdot, 0\},$$ N(i) is the neighborhood of node i and [TeX:] $$\omega_{i j}$$ is the weight on [TeX:] $$\text { edge }(i, j) . \text { Here, } \zeta \geq 0$$ represents the strength of Shannon- Boltzman entropy, which is a regularized term for a good population diversity when we want to optimize a global cost functional [27]. However, for the MFED-based epidemic models in this paper, we do not need a global cost functional and thus [TeX:] $$\zeta$$ is set to 0.

Remark 2. The change of population in a generic state i is computed by the difference between the inflow of population from state i’s neighbors and outflow of population from state i. The direction of population flow between state i and state j is determined by the payoff of choosing state i and the payoff of choosing state j, where state j must be in the neighborhood of state i. The MFEDs are highly related to the structure of the strategy graph [TeX:] $$\mathcal{G} .$$

A proof of Theorem 1 is given in the appendix. The MFEDs, proposed in Theorem 1, are pairwise comparison evolutionary dynamics [11]. Unlike the RDs, MFEDs set the probability of switching from the current strategy i to strategy j proportional to the differences between these strategies’ payoffs. One agent will compare its current payoff with the payoffs of the strategies in its neighborhood to determine whether it will switch.

B. Mean Field Evolutionary Dynamics in SIR, SIVR, and SIDR

With the derived payoff functions in (8), (9), and (10) and the general form of MFEDs in the Theorem 1, we will construct SIR, SIVR, and SIDR epidemic models based on corresponding strategy graphs shown in Fig. 1.

Following the strategy graph in Fig. 1(a) and putting payoff functions (8) into (11), we obtain the MFEDs for the SIR as follows:

where [TeX:] $$\zeta=0$$ and all weights are set to be unity. The first partial differential equation (PDE) describes the change in the fraction of populatoin that is susceptible. [TeX:] $$\rho_{I}\left[F_{s}-F_{I}\right]^{+}$$ is the fraction of the population, who recover from the infectious disease without immunity and thus become susceptible in a unit of time. [TeX:] $$\rho_{s}\left[F_{I}-F_{s}\right]^{+}$$ is the fraction of the susceptible population who get infected in unit time. The second and third PDEs can be interpreted in a similar way.

Likewise, following the strategy graph in Fig. 1(b) and putting payoff functions (9) into (11), we obtain the MFEDs for the SIVR as follows:

where [TeX:] $$\zeta=0$$ and all the weights are set to be unity. Unlike the SIR model, a susceptible individual can not only become infected by the original virus but also by the mutated virus in the SIVR model. Therefore, as shown in the first PDE in (13), the change in the fraction of the population that is susceptible in a unit of time is determined by the fractions on the population who become infected with the original virus and the mutation. The fraction of the population that becomes infected with the mutated virus has a similar effect on the fraction of the population that has recovered.

Following the strategy graph in Fig. 1(c) and putting the payoffs functions (10) into (11), we obtain the MFEDs for the SIDR model as follows:

where [TeX:] $$\zeta=0$$ and all the weights are set to be unity. Compared with the SIR model, the SIDR model considers the people who die from the infectious disease. As shown in Fig. 1(c), the common neighbor of state S and state R is state I, which is the same as that in the SIR model. Therefore, the change in the fraction of the population that is susceptible and the change in the fraction of the population that has recovered are computed in the same way as those in (12). The change in the faction of the population that is infectious in a unit of time is also influenced by the number of deaths as shown in the second PDE in (14).

IV. NUMERICAL ANALYSIS

MFEDs are constructed with the payoff functions derived from the classical RDs. Both MFEDs and RDs can describe the transmission of COVID-19 among people. In this section, we compare MFEDs with traditional RDs under the classic SIR epidemic model and the extended SIVR epidemic model via simulations.

A. MFEDs and RDs in SIR Model

In this simulation, the infection rate [TeX:] $$\beta$$ is 0.0005 and the recovery rate is 0.2. The total population size N is 1,000, with 960 susceptible individuals, 30 infected individuals, and 10 recovered individuals in the beginning. As shown in Figs. 2(a) and 2(b), there are several differences between the results. The peak size of the infected population is 312 using the MFEDbased model, higher than 239 with the RD-based model. The convergence size of the recovered population is 719 with the MFED-based model, lower than 895 with the RDbased model. Moreover, the susceptible population decreases monotonically with RDs, while it will increase slightly with MFEDs when the infected population reaches its peak value.

In Fig. 2(c), we compare the MFED-based SIR model with the RD-based SIR model in terms of the infection rate , which varies from 0.0003 to 0:001. A higher infection rate will lead to a larger peak number of infections in both models. The increasing rate of the number of infections in the MFEDbased SIR model is higher than that in the RD-based SIR model. Moreover, the increase in the infection rate has a more obvious impact on the peak number of infections in the RDbased SIR model. In Fig. 2(d), we compare the MFED-based SIR model with the RD-based model in terms of the population size N, which varies from 800 to 1200. The peak number of infections in the MFED-based SIR model is higher than that in the RD-based SIR model.

B. MFEDs and RDs in SIVR Model

In this simulation, the infection rates [TeX:] $$\sigma_{1} \text { and } \sigma_{2}$$ are 0.0005 and 0.0016, respectively. The recover rates [TeX:] $$\mu_{1} \text { and } \mu_{2}$$ are 0.0666 (1\15) and 0.1428 (1\7), respectively. The total population size N is 1000 with 949 susceptible individuals, 30 individuals infected by the original virus (virus 1), 20 individuals infected by the variation of the original virus (virus 2), and 1 recovered individual in the beginning. Similar convergence behaviors have been depicted in Figs. 3(a) and 3(b). There are some major differences between RDs and MFEDs. First, the peak numbers of the persons infected by virus 1 and virus 2

Fig. 2.

Comparison between MFEDs and RDs in the SIR model.

with MFEDs are smaller than those with RDs. Second, with MFEDs, the number of susceptible persons will increase when the number of infected persons reaches its peak value, while the number of susceptible persons will remain 0 with RDs.

In Figs. 3(c) and 3(d), we analyze the effect of the population size N on the original virus (virus 1) and the mutated virus (virus 2), respectively. The infection rate of virus 1 is assumed to be lower than that of virus 2. In both models, a larger population size will lead to a larger peak number of infections. However, the variation of infection rates of virus 1 and virus 2 has a more obvious impact on the peak number of infections in the RD-based SIVR model.

V. MODEL THE TRANSMISSION OF COVID-19

In this section, we model the transmission of COVID-19 with the proposed MFEDs in (14). In Section V-A, our data set and evaluation metrics are given. In the following Section V-B, we analyze the overall evolution of COVID-19 in the SIDR model. In Section V-B, we fit our model into the real statistics from Wuhan, China. In Sections V-C and V-D, we analyze the effect of one-time social distancing and seasonality of COVID- 19, respectively.

A. Data Set and Evaluation Metrics

1) Data set: We obtained the COVID-19 statistics from the National Health Commission of the People’s Republic of China, the Chinese Center for Disease Control and Prevention, as well as the coronavirus research center at Johns Hopkins University. These statistics depict the evolution of COVID-19 in Wuhan, China from Feb. 12th, 2020 to Apr. 25th, 2020. We determine to use these statistics for the reasons as follows:

The number of infections reached its peak value on Feb. 18th and decreased to 0, for the first time, on Apr. 25th. The statistics remained unchanged between Apr. 25th and May. 8th. Therefore, the statistics between Feb. 18th and Apr. 25th can depict the complete transmission process of COVID-19 in Wuhan, China.

Between Feb. 18th and Apr. 25th, Wuhan was completely locked down so that there was no significant population flow, making it relatively acceptable for us to assume a constant population size. This condition is needed for both MFED-based and RD-based models.

The whole process of the transmission of COVID-19 in Wuhan, China, ended in around two months. In this relatively short period, we can reasonably assume that the people were immune to COVID-19 after recovery. This is another condition needed for both MFED-based and RD-based models.

Fig. 3.

Comparison between MFEDs and RDs in the SIVR model. Virus 1 refers to the original virus and virus 2 refers to the mutated virus.

2) Evaluation metrics: The coefficient of determination [TeX:] $$\left(R^{2}\right)$$ and the average relative bias are employed to evaluate the goodness of fit of the results. [TeX:] $$R^{2}$$ measures how successful the fit is in explaining the variation of the data, which is computed by

where [TeX:] $$y_{i}, \hat{y}_{i}, \bar{y}$$ are the observed value, estimated value, and mean of observed values, respectively. The average relative bias is computed by

In the simulation, the statistics of COVID-19 on Feb. 12th serve as the initial values for the SIDR model. The total population size is assumed to be N = 50.333, which is the final total number of confirmed cases on Apr. 20th, 2020. The effective reproduction number [TeX:] $$R_{0}$$ determines the potential of the pandemic, which is defined as the average number of secondary infections caused by a single infected individual. As suggested in [28], [TeX:] $$R_{0}=\eta N$$ of COVID-19 should be between 2 (during the summertime) and 2.5 (during the wintertime). Therefore, the infection rate [TeX:] $$\eta$$ is set to 0.00004. According to the best-fit model parameters, the initial recovery rate and the initial death rate are set to 0.021 and 0.03, respectively.

Fig. 4.

Transmission of COVID-19 in the SIDR model.

In Fig. 4, we depict the evolution of COVID-19 in the SIDR model. The transmission dynamics in (14) have been applied to update the population distribution on each state. The peak number of infected cases, the final number of recovered cases, and the final number of deaths are 38; 170, 46.086, and 3.609, respectively. These critical numbers are close to the corresponding observed numbers.

B. Data Fitting

The data fitting results of MFEDs and RDs are shown in Fig. 5. MFEDs outperform RDs in terms of the predictions on

Fig. 5.

Data fitting results of MFEDs and RDs.

the number of infections and the number of recoveries. The infection rate and recovery rate can vary significantly during the spread of COVID-19, mainly due to the implementation of social distancing and the improvement in the treatment of COVID-19. Therefore, dynamic adjustments to the recovery rate and the death rate, shown in Table I, are implemented to improve the accuracy of prediction for both MFEDs and RDs.

Both MFEDs and RDs can capture the overall trend of the evolution of COVID-19 in Wuhan, China. However, RDs fail to make accurate predictions on the number of infections after 40 days as shown in Fig. 5(c), which leads to a high relative bias as shown in Table II. The MFED-based SIDR model’s estimation of the final number of deaths is 3,609, while the RD-based SIDR model’s estimation of the final number of deaths is 2,697. Given that the real number of deaths is 3,869 on Apr. 25th, we can conclude that MFEDs also outperform RDs on the final number of deaths. In Table II, the performances of MFEDs and RDs are compared based on the two evaluation metrics. MFEDs outperform RDs significantly in terms of the prediction of infections and recoveries.

C. Effect of One-time Social Distancing

We analyze the effect of one-time social distancing (SD) in the absence of seasonality in this subsection. Social distancing measures have been taken at the beginning of the transmission process of the COVID-19. From Figs. 6(a) to 6(e), the duration of social distancing varies from one month to an indefinite period. The length of the period is denoted in the blue shaded region. The reduction of infection rate varies from 0 to 80% representing the distinct effectiveness of SD. The total population size is assumed to be N = 50; 000. The basic infection rate is [TeX:] $$\eta=0.00004$$ without considering any seasonality, i.e., [TeX:] $$\eta$$ is held as a constant. With [TeX:] $$\eta=0.00004,$$ the effective reproduction number is [TeX:] $$R_{0}=N * \eta=2.$$

As shown in Fig. 6, depending on the strength of the impact of the infection rate, the peak number of infections varies from 10.000 to 27.500. To restrain the spreading of COVID-19, we need the peak number of infections to be less than the total capacity of patients in the hospital. These significant reductions in the peak number, due to the SD, are thus of vital importance. For a slight degree of SD, which

Fig. 6.

One-time SD without seasonality. The blue shaded region refers to the quarantine period, which varies from one month to indefinite length. Due to the SD during the quarantine period, the infection rate of COVID-19 is assumed to be reduced by 40% to 80% depending on the level of SD. When the blue shaded region ends, the infection rate of COVID-19 is set to the normal value.

leads to a 40% reduction in the infection rate, there will be no significant increase in the number of infections when the quarantine period ends. For a severe degree of SD, which leads to a more than 60% reduction in the infection rates, there will be a significant increase in the number of infections when the quarantine period ends. The critical reason behind this phenomenon is that only when the number of immune individuals reaches a certain level can the spreading of the COVID-19 be stopped. Therefore, after a severe degree of SD, returning to the normal SD step by step and developing a vaccine should be of vital importance.

D. Seasonality of COVID-19 through the post-pandemic period

We analyze the seasonality behavior of COVID-19 in the post-pandemic period in this subsection. Specifically, we consider the effect of the duration of COVID-19 immunity [TeX:] $$d_{0},$$ the proportion of immunity loss [TeX:] $$L_{0},$$ and the seasonal variation [TeX:] $$v_{0}.$$ The duration of COVID-19 immunity [TeX:] $$d_{0}$$ refers to the length of the period when the recovered individuals can keep their immunity to COVID-19. After the immunity duration [TeX:] $$d_{0},$$ a certain proportion [TeX:] $$L_{0}$$ of the recovered population will lose their immunity. The seasonal variation [TeX:] $$v_{0}$$ refers to the variation of the infection rate of COVID-19 from wintertime to summertime.

In Fig. 7(a), periodic outbreaks of COVID-19 are depicted for [TeX:] $$d_{0}=6$$ months and [TeX:] $$d_{0}=12$$ months. The immunity loss is [TeX:] $$L_{0}=0.8, \text { i. e., } 80 \%$$ of the recovered population will lose their immunity after the duration [TeX:] $$d_{0}.$$ There is no seasonal variation [TeX:] $$\left(v_{0}=1\right).$$ The peak number of infections will gradually decrease through the outbreaks. The period of the outbreaks of COVID-19 is closely related to the immunity duration.

In Fig. 7(b), we show the impact of the proportion of immunity loss on the seasonality of COVID-19. The immunity duration is set as [TeX:] $$d_{0}=6$$ months and there is no seasonal variation [TeX:] $$\left(v_{0}=1\right).$$ The immunity loss [TeX:] $$L_{0}$$ varies from 0.4 to 0.8. Higher immunity loss will yield more severe outbreaks of COVID-19 every six months. However, the peak number of infections will decrease gradually.

In Fig. 7(c), the effect of seasonal variation is shown. The immunity loss is [TeX:] $$L_{0}=0.8$$ and the duration is [TeX:] $$d_{0}=6$$ months. Without loss of generality, we assume that the infection rate is higher in the winter. (One can also set a higher infection rate in the winter. In that case, the roles of summer and winter will switch but the results will be similar.)

In particular, the infection rate is [TeX:] $$\eta$$ during the wintertime and [TeX:] $$v_{0} \eta$$ during the summertime, where [TeX:] $$v_{0}$$ could be 0:2 or 0.8. The winter of 2019 is regarded as the starting point. A higher seasonal variation would significantly reduce the peak number of infections in the summertime but yield a more severe outbreak of COVID-19 in the following winter. The overall trend of the peak number of infections is also decreasing.

VI. CONCLUSION

Mean field evolutionary dynamics, inspired by optimal transport theory and mean field games on graphs, have been proposed to model the evolution of COVID-19. This approach has been compared with the commonly used replicator dynamics with numerical simulation results. Applied into the new SIDR model, the mean field evolutionary dynamics have been seen to outperform the replicator dynamics in fitting the COVID-19 statistics of Wuhan, China. In our simulations, we have observed that one-time social distancing can reduce

Fig. 7.

Seasonality Analysis of COVID-19: (a) The length of immunity period varies from 6 months to 12 months, (b) the proportion of population who will lose their immunity to COVID-19 varies from 40% to 80%, and (c) the ratio of the infection rate of COVID-19 in winter to the infection rate of COVID-19 in the next summer varies from 0.2 to 0.8.

the peak number of infections significantly while a second outbreak of COVID-19 can arise after a high level of social distancing. Finally, as per our model, people’s limited length of immunity will yield a periodic outbreak of COVID-19, and higher seasonal variation of infection rates will reduce the peak number of infections in summer but lead to a more severe outbreak in the following winter.

APPENDIX A PROOF OF THEOREM 1

In this appendix, we derive the MFEDs in Theorem 1. First, we construct the Riemannian manifold on which the MFEDs are defined. Then we prove that the dynamics shown in Theorem 1 are well-defined evolutionary dynamics on the Riemannian manifold.

To measure distance in the state space [TeX:] $$\mathcal{P}(\mathcal{S}),$$ we need to define the following Wasserstein metric [29].

Definition 1. Given two discrete probability functions [TeX:] $$\rho^{0}, \rho^{1} \in \mathcal{P}_{o}(\mathcal{S}),$$ the Wasserstein metric W is defined by

where [TeX:] $$\nabla \Phi: \mathcal{S} \times \mathcal{S} \rightarrow \mathbb{R}$$ is given by

where [TeX:] $$\Phi$$ is a function and [TeX:] $$\Phi: \mathcal{S} \rightarrow \mathbb{R} \text { and } \omega_{i j}$$ is the weight on edge (i,j).

Besides the Wasserstein metric, we also need the following inner product [TeX:] $$g^{W}$$ to construct the Remannian manifold [TeX:] $$\left(\mathcal{P}_{o}(\mathcal{S}), g^{\tilde{W}}\right).$$

Definition 2. For any two tangent vectors [TeX:] $$\sigma^{1}, \sigma^{2} \in T_{\rho} \mathcal{P}_{o}(\mathcal{S}),$$ define the inner product [TeX:] $$g^{W}: T_{\rho} \mathcal{P}_{o}(\mathcal{S}) \times T_{\rho} \mathcal{P}_{o}(\mathcal{S}) \rightarrow \mathbb{R}$$ by

where [TeX:] $$\sigma^{i}=-\operatorname{div}\left(\rho \nabla \Phi^{i}\right) \text { for } i=1,2 . T_{\rho} \mathcal{P}_{o}(\mathcal{S})=\left\{\left(\sigma_{i}\right)_{i=1}^{n} \in\right. \left.\mathbb{R}^{n}: \sum_{i=1}^{n} \sigma_{i}=0\right\}$$ is the tangent space at a point [TeX:] $$\rho \in \mathcal{P}_{o}(\mathcal{S}).$$ [TeX:] $$\theta_{i j}$$ is the discrete probability on edge (i,j) which is defined by

where [TeX:] $$F_{i}: \mathcal{P}(\mathcal{S}) \rightarrow \mathbb{R}$$ is the payoff function and [TeX:] $$d_{i}$$ is the degree of node i (i.e. the total number of nodes in [TeX:] $$N(i)).$$

With the state space in (2) and the inner product in Definition 2, we can construct the Riemannian manifold [TeX:] $$\left(\mathcal{P}_{o}(\mathcal{S}), g^{W}\right)$$ [30], [27]. In this regard, we can give the proof of Theorem 1 as follows:

Proof. Given the tangent space [TeX:] $$T_{\rho} \mathcal{P}_{o}(\mathcal{S})=\left\{\left(\sigma_{i}\right)_{i=1}^{n} \in \mathbb{R}^{n}:\right. \left.\sum_{i=1}^{n} \sigma_{i}=0\right\},$$ there exists [TeX:] $$\Phi$$ such that [TeX:] $$\sigma=-\operatorname{div}(\rho \nabla \Phi)$$ for any [TeX:] $$\sigma \in T_{\rho} \mathcal{P}_{o}(\mathcal{S}) . \mathrm{As} \frac{d \rho}{d t}=\left(\frac{d \rho_{i}}{d t}\right)_{i}^{n}$$ is in [TeX:] $$T_{\rho} \mathcal{P}_{o}(\mathcal{S}),$$ we have

The noisy potential is given by

which is the summation of the potential and the Shannon- Bolztman entropy. Then we have

With (17) and (19), and the definition of gradient flow of [TeX:] $$-\bar{F}(\rho)$$ on the Remannian manifold [TeX:] $$\left(\mathcal{P}_{o}(\mathcal{S}), g^{W}\right),$$ we derive

As the above is true for all [TeX:] $$\left(\Phi_{i}\right)_{i=1}^{n} \in \mathbb{R}^{n},$$ we finally obtain

Replacing [TeX:] $$\theta_{i j}$$ with (16), the mean field evolutionary dynamics in Theorem 1 are proved.