Using evolutionary game dynamics to reproduce phase portrait diversity in a pandemic scenario

Modeling the time evolution of a pandemic in distinct social and economic environments has proven to be a great challenge for scientists due to the strong coupling of infection dynamics to complex collective human behavior. In this work, we propose a relatively simple extension of the classical Susceptible-Infected-Recovered (SIR) model that encompasses evolutionary game imitation dynamics and is capable of reproducing basic features of real data such as recurrent waves with suppression and re-bounce of cases. In particular, we depict diﬀerent quantitative and qualitative aspects of infection time series by representing them as trajectories in a modiﬁed phase portrait.


Introduction
Mathematical modeling of epidemic processes has a longstanding history providing important clues on several relevant aspects related to the likely outcome of outbreaks. Mathematical modeling strategies have provided a better understanding of the influence of distinct underlying epidemic characteristics such as the contagious rate, population behavior, vaccination, average contamination time, among others. Recent progress has significantly contributed to help public health intervention decisions and to project future growth patterns, especially in the presently ongoing COVID-19 pandemic.
The epidemiological modeling can be subdivided into two main branches: individual/agent-based (network mod-els) [1][2][3] and mean-field/compartmental models (ODEs) [4][5][6][7]. The latter can be derived from the first. Burda [2] developed an agent-based model that simulates the disease transmission of a COVID-like epidemic, where only non-pharmaceutical interventions (NPIs) are applied. That model distinguishes local and non-local transmission. When both are indistinguishable, the mean-field dynamics quite reproduce the SIR dynamics. To study the speed of epidemics spreading in such models, see Ref. [3]. In the case of compartmental models, Maier and Brockmann [8] study  in the early stages of China. The generalized SIR-X model takes into consideration the social distancing efforts such as quarantine and general confinement. They fairly reproduce the data of various provinces and explain the observed subexponential growth. Dehning et al. [9] investigate the effective growth rate related to the interventions in Germany. Adding a reporting delay and a time-dependent rate, they built a non-autonomous system (see also Ref. [10] and [11]). With the help of Bayesian inference based on Markov chain Monte Carlo, the inference of parameters was made. Through this paper, we are interested in the compartmental models.
Due to a combination of personal and socioeconomic benefits, individuals during a pandemic constantly choose between maintaining a normal behavior -as before the outbreak -or playing protective behaviors, namely mask-wearing, social-distancing, and home-office. That collective decision-making process occurs dynamically concurrently with the pandemic's evolution. Poletti et al. [12] studied such processes by proposing a modified SIR model motivated by evolutionary game theory (imitation dynamics). Susceptible individuals can adopt normal or altered behaviors. The model is able to produce a wide range of dynamical regimes, see also Ref. [13]. Game theory was extensively employed in biological [14] and economic [15] models but appears recently in the branch of epidemiological models. The main classifications of game modeling in that context are self-learning [16][17][18] and imitation [12,[19][20][21][22].
Behavioral science is a way to study the pandemic response [23]. Along this line, various works appeared quite recently devoted to shed light in available COVID-19 data [24][25][26][27]. Amaral et al. [24] extended Poletti's model to split the susceptible infectious stock into normal and altered. In [25], a SEQIHR model (susceptible, exposed, quarantined, infected, hospitalized, recovered) with behavioral dynamics is studied. The infected individuals were distinguished by symptomatic and asymptomatic. The analysis illustrates that the natural shield immunity is worthless in the actual scenario.
In this work, we advance in the study of theoretical epidemiology by introducing some models based on evolutionary game theory, especially imitation dynamics. Our goal is to consider a simple extension of Poletti's model capable of reproducing relevant dynamical features of a real pandemic scenario, such as the current COVID-19 pandemics. The model combines compartmental mean-field SIR-like dynamics with imitation dynamics from game theory behavioral modeling. In particular, we aim to fill the gap in the current literature of how recurrent waves and pandemic persistence can be modeled in a phase portrait without the inclusion of vital dynamics, nor waning immunity in the recovered population. Grounded on the recently learning from the COVID-19 pandemic, we add new contributions to the payoffs related to normal and altered behavior. Some research groups are focused on the patterns of the risk diagrams [28]. We start by revisiting the phase diagrams for the basic SIR model. A comparison with the phase-diagram of more sophisticated models is provided. Qualitatively, we show that the novel terms introduced in this work can capture real decisionmaking dynamics. Besides the data-driven analysis to infer the model parameters, the R t -S plane directly captures the size and time spacing between recurrent waves. Although we just consider a few compartments to leave the mathematical model as simple as possible, we provide a qualitative tuning of the basic parameters that reproduce actual risk diagrams of the COVID-19 epidemic. Finally, we do a convenient change of variables to make the model directly comparable with the reported data (such as in [28]).

Modeling human behavior in epidemics
The starting point for modeling epidemics is the understanding of models that takes into consideration a lesser number of assumptions. That way, the SIR model [10,  29] assumes that the typical interaction time in an epidemic is much faster than demographic dynamics such as birth and death (vital dynamics). Here, the susceptible population S is exposed and possibly infected by the infectious population I at a per capita rate linearly proportional to the I-population. That proportionality constant is just the typical time between contacts. On average, the population fraction that is infectious will decrease at a constant per capita rate. Such a rate is the so-called characteristic time of the infection. That proportion is removed from the disease dynamics increasing R, i.e. some of the infectious turn into deceased or recovered. Here we neglect vital dynamics, soṠ +İ +Ṙ = 0 and the compartments are normalized such that S+I +R = 1. Across this paper, we work with population fractions (dimensionless) that are bounded by one. Even though simple, this model is idealized as a free infection epidemiological model. The assumptions are based on standard characteristics of the spreading disease such as mean infectious period γ −1 and typical time between contacts β −1 . On the other hand, in the absence of effective pharmaceutical interventions, human behavioral change plays a key role in the infection dynamics.
During a state of pandemic awareness, there is a fraction of susceptible individuals spontaneously that changes their conduct from a normal pre-pandemic social behavior to an altered social behavior by performing protective efforts such as social distancing. On account of this, the susceptible stock should split into population fractions playing each distinct strategy, namely S n and S a . A schematic representation of our SIRlike compartmental model is shown in Fig. 1. The exchange between the different strategies could be visualized as an imitation dynamics well-known from evolutionary game theory. In that way, the players choose to stay or change their strategy by comparing the payoffs. Hence, the term take the general two-strategy form [22], is the population fraction playing strategy N at a time t and 1 − x(t) is the population fraction playing strategy A at the same time t. The constant τ represents the characteristic time of the imitation process. Moreover, the term ∆p is a function related to the payoff difference between each strategy. In the imitation process, the players only compare the payoffs, i.e. ∆p = p n (I) − p a (I). the number of contacts is reduced. The payoffs are modeled as functions of the perceived populations. All players have a risk of infection (a cost), clearly, that cost is higher for individuals playing normal behavior, assuming a linear dependency with the fraction of infectious I. In addition, altered players pay an additional constant cost due to individual issues, say k 0 . Therefore, the whole population fraction playing altered behavior S a undergoes extra issues that increase as that fraction increases, namely food shortages and mass layoffs. Especially in emerging countries like Brazil, these two issues are very relevant in this ongoing COVID-19 pandemic. To consider the simplest model, we assume that cost depends linearly on the fraction of susceptible individuals playing altered behavior S a .
Moreover, the likelihood of infection is reduced when a sufficiently high percentage of the population has been previously exposed to the virus and becomes recovered. Since in a pandemic this herd immunity phenomenon might be widely perceived by the overall population, the knowledge of how many people were already recovered can affect individual and collective considerations about social distancing. An increasing number of recovered generate a false perception of herd immunity for which the susceptible population tends to loosen up the protective measures. Thereby generating a flux from altered to normal, manifested by a benefit to the altered players and a cost to the normal ones.
To recap, we assume that susceptible individuals may choose one of two strategies: normal behavior with no socioeconomic cost and high risk of contagion, with a payoff given by p n (I) = −k n 1 I + k n 2 R; and altered behavior with high socioeconomic cost and low risk of Based on the above assumptions, all constants are real, k 0 > 0 and k n 1 > k a 1 . Differences in behavior of infected individuals are considered to be less relevant to overall dynamics and can be evaluated in previous works [24]. The overall dynamics is given by Since there can only be exchange between compartments S n and S a if they are both populated, we may rewrite the exchange rate as Φ = S n S a ∆p, where the payoff function ∆p is written as As a matter of fact, by setting k 2 = 0 and k 3 = 0 in Eq. 5 we recover Poletti's model [12]. One of our goals is to determine the coefficients k 0 , k 1 , k 2 , k 3 that are suitable for describing real COVID-19 data of a particular population. This will be done on section 5.

Phase portrait trajectories and risk diagrams
Through this paper we are particularly interested in studying the phase portrait trajectories of our system and in comparing it with real COVID-19 data. As a matter of fact, we were motivated by the risk diagrams developed in Ref. [28], which is similar to a phase portrait, given by time-implicit trajectories in R t − −I plane, where R t is the effective reproductive number. In the context of a SIR-like model, the effective reproductive ratio R t can be defined as is increasing, otherwise R t < 1 the outbreak is decaying. Consequently, each epidemic 'wave' has a peak (or valley) at R t = 1.
For the simple SIR model, this results in R t = βS(t)/γ. Since the initial condition S(0) is close as possible to one, R 0 = β/γ. The risk diagram can be obtained analytically for the SIR model by direct integration of dI/dS [10], which lead us to where we consider S close to one and I close to zero at the beginning of the outbreak. Therefore In the Fig. 2 we show the risk diagram for the SIR model and R 0 = 3. Note that the condition R t > 1 (< 1) the outbreak in increasing (decreasing) is satisfied. When I = 0, the effective reproductive number take two possible values: , where W is the Lambert W-function. We are also interested in another plane of the phase portrait given by R t vs. S(t), which is just a straight line with slope R 0 for the SIR model.
A richer variety of phase portrait trajectories can be obtained with Poletti's SIR-like model [12], which is given by setting k 2 = k 3 = 0 in Eqs.(1)-(4) and Eq. (5). Now, the payoff ∆p depends on I(t) and this implicit time dependence may induce several flips in the susceptible population behavior, which alternates between normal and altered compartments. This alternation of social behavior driven by imitation dynamics produces recurrent waves of infection which decrease with time. This phenomenon is well depicted by another section of the complete phase portrait: the R t − S plane, where S(t) = S n (t)+S a (t), as shown in Fig. 3. Effectively, the system exhibits an alternating behavior switching between two SIR-like lines with β n or β a as S(t) decreases in time, indicating three epidemic waves for this given set o parameters. Although this model is capable of modeling recurrent waves driven by social behavioral changes, it has no capability of reproducing qualitatively well real data on the risk diagram R t − I [28], nor on the R t − S plane (as it can be seen in section 5). Thus, in order to extend the imitation dynamics model and to be able to better reproduce phase portrait trajectories of a real pandemic, our model explores nonzero values of the parameters k 2 and k 3 and gives rise to more complex and realistic diagrams. In the next section, we make a change of variables in order to further compare our extended model with real data for COVID-19 in the form of phase portrait trajectories.

Change of variables
In this section we show there is a convenient change of variables to compare a two-state susceptible population with usual collected data. Instead of dealing with S n (t), S a (t), I(t), we can switch to S(t), I(t), R t (t) by choosing R t (t) = βn γ S n (t) + βa γ S a (t), hence the mathematical epidemiological model becomes, where S, I and R t are considered independent variables and ∆p is rewritten as where k 0 > 0, k 1 > 0, k 2 > 0, but the sign of k 3 could be either positive or negative.

Some analytical and numerical analysis
Let us study the behavior of the system in a linear approximation via linear stability analysis. We are interested in the disease-free fixed points, i.e. I * = 0. From Eq. 13 we find three fixed points, namely with the constraints 0 < S * < 1 and β a /γ < R * t < β n /γ, the feasibility of the fixed points is guaranteed. This represents a plane of fixed points that can be stable or not. Their linear stability depend on the parameters and also with R * t . The value of R * t , i.e. the asymptotic value of R t depend on the initial condition. An analytical treatment of its dependency is quite intractable due to the complexity of the model. In Fig. 4 we show the asymptotic value of R t for a given set of parameters. The numerical analysis is made with the initial conditions is a good approximation used when R(0) = 0 and I(0) ≈ 0 (see Eqs. 9 and 10). For such set of parameters we find that S * 2 is not feasible. In the thin range of parameters for which S * 3 is feasible, it is always unstable. Whereas S * 1 is feasible and is stable for R * t < 1. From the recent learning of COVID-19 and other outbreaks, we know that is favorable to be in the neighborhood of R t = 1. By fixing R 0 = 2.97, roughly the largest R * t in Fig. 5, we show that k 2 and k 3 can be responsible for approaching R * t to one.

Quantitative comparison with COVID-19 data
In order to compare our model with real data from COVID-19 spreading dynamics, we employ a simple data-driven approach to infer the observable functionŝ I(t),R(t),Ŝ(t) andR t for a given population. Within this scope, as input data for a given location, we will use the time series of cumulative confirmed cases C(t), cumulative number of confirmed COVID-19 deaths D(t) [30], and an estimation of the actual infection fatality ratio of that population δ [30,31]. We assume that confirmed cases possess a timedependent sub-notification factor f (t) that depends on the capability of a population to detect new cases, and it can be estimated by f (t) = δ −1 D(t)/C(t). Moreover, due to the quite universal biological recovery time for COVID-19 which is given by T inf ectious ≈ 10 days, we take γ = T −1 inf ectious = 0.10 as a constant parameter regardless of the population location. Then, we proceed by assuming that a reasonable proxy for the number o infected individuals in a given time can be evaluated aŝ From Eq. (4), one can writê The time series for cumulative number of cases and deaths was obtained from Ref. [30], and the estimated IFR was considered as δ = 0.01 [31].
By proceeding this way, we are able to extract from real data the daily time series corresponding toŜ(t), I(t) andR t (t) that describe our epidemiological model for a given population. Further, we search for the optimal fitting parameters for Eqs. (11)-(13) that best reproduce the numerical time series. We minimize the error between our model fit {S(t), I(t), R t (t)} and the data-driven time series {Ŝ(t),Î(t),R t (t)} by performing a non-linear least-squares minimization using the free Python package lmfit.
In Fig. 6 we show how real COVID-19 data from Spain compares with our model for a given set of optimal parameters. The top panel presents three graphs that illustrate the time series forŜ(t),Î(t) andR t (t), where data points are depicted as blue circles and the best fits are denoted by black lines. We observe an excellent agreement between real data and our model for both S(t) and I(t) curves, whereas the average behav-ior of R t (t) is also very well captured by our fit if we neglect fluctuations in smaller time scales.
The bottom panel gives us further insight on the distinct nonlinear behavior of the dynamic system by displaying the phase portrait 3-dimensional trajectories in 2-D sections: the risk diagram given by the R t −I plane, the R t − S plane, and the usual I − S plane. The risk diagram fit captures the main aspects of the phase portrait trajectory: a large wave peak corresponding to a rapid increase in I(t) and the first crossing of the R t = 1 axis; then a good suppression of cases where strong nonpharmaceutical interventions were implemented which kept R t < 1 until cases were very low; followed by a smaller second wave, and finally a re-bounce of cases represented by a smaller loop in the risk diagram. From a physical point of view, it matches very well the historical hallmarks observed in COVID-19 cases and deaths in Spain until November 2020 [30,31].
In Fig. 7 we inspect the case of the State of Sao Paulo, which accounts for approximately 22% of Brazil's total population. As a matter of fact, the top graphs show a very good qualitative and quantitative reproduction of real data (given by circles) by our model fit (given by full curves) for all the time series representing S(t), I(t) and R t (t). The bottom panel displays Ref. [32], and the estimated IFR was considered as δ = 0.01 [31]. a very good qualitative agreement of the fits despite the quite noisy nature of the data, which must come from the lack of consistent testing in Brazil. The I − S plane displays clearly two peaks of infection as S(t) decreases with time, where the first wave is very long and the second peak is smaller than the first. The risk diagram R t − I shows a quite distinct behavior from the previous case: the first wave was considerably flatter than the Spanish one, presenting smaller values of R t but for a longer period; the suppression of cases was much weaker, lasting much less time in the period where R t < 1 and barely reaching values of R t smaller than 0.9; finally, the second wave comes in higher values of I(t) and the phase portrait trajectory resembles a spiral that rapidly displays a third peak. We may interpret this result as pointing out that Brazil has implemented NPIs quite early in the epidemic course, but they were not strong enough to produce a strong suppression of cases, and the socioeconomic burden produced long and low-spaced waves.
In Fig. 8 we utilize our approach to fit data from the United Kingdom (UK). Although the top graphs still present a good quality fit, the lower panel shows only a qualitative agreement with the data. However, we point out that the model was able to capture important fea-tures in the phase portrait trajectory: a second wave larger than the first one, a very small loop in the risk diagram indicating a re-bounce of cases due to the lack of stringent behavior, and trajectories which oscillate quite close to the R t = 1 axis. These features are quite often present in real data trajectories in risk diagrams [28], and to the best of our knowledge, the Poletti-like models of Refs. [12,24] are not able to reproduce.

Concluding remarks
In pandemics such as the current COVID-19 one and in the absence of pharmaceutical interventions, recurrent waves of infections can be caused by several distinct factors, namely: (i) vital dynamics and the birth of new susceptible individuals; (ii) seasonality in infections caused by respiratory virus; (iii) appearance of new pathogenic strains capable of partial reinfection; (iv) waning immunity; (v) governmental NPI measures and dynamic social adhesion to preventive NPI recommendations [33][34][35]. In this work, we explore the possibility of explaining important qualitative and quantitative features of the nonlinear dynamics of a multiplewave pandemic based solely on the latter ingredient, i. e., item (v). Ref. [30], and the estimated IFR was considered as δ = 0.01 [31].
We show that a SIR-like model based on evolutionary game theory, with imitation dynamics for susceptible individuals, can be useful to understand how social behavior plays its role in a pandemic scenario. Moreover, we take as independent variables the time series of S(t), I(t) and R t (t), and explore their phaseportrait trajectories. One of the phase-portrait sections are the risk diagrams in the R t −I plane [28], which has shown practical relevance in depicting a great variety of non-linear behavior in the pandemic dynamics for distinct countries. We show that our model can reproduce important qualitative and quantitative features of data such as recurrent waves of different sizes, behaviorindiced suppression of infections, and re-bounce of cases.
The imitation dynamics is based on the fact that susceptible individuals may choose between maintaining a normal behavior with no socioeconomic cost but high risk of contagion or adopting an altered behavior with high socioeconomic cost and low risk of contagion. Our extended model takes into account that in an information abundant pandemic scenario, individuals will evaluate the best strategy by evaluating a payoff ∆p based not only on the socioeconomic cost (k 0 parameter) and perceived infection risk (term proportional to k 1 I), but also in the growth of new cases (perceived as k 2 R t ) and how many people are still susceptible (perceived as k 3 S).
With one extra compartment, our simple SIR-like model can reproduce the size and time spacing between recurrent waves of the COVID-19 epidemic which are not captured by the standard SIR model. Moreover, we show that the novel terms introduced in this work can reproduce real decision-making dynamics by comparing our model to real data. Therefore our model could be extended to more than two different strategies by separating the susceptible groups in many levels with an increase of socioeconomic cost and a decrease of risk of contagion. In addition, greater accuracy in data fitting can be achieved by combining other factors that produce recurrent epidemic waves (items (i) to (iv) in the first paragraph of this section).