Identification and Control of Game-Based Epidemic Models

The effectiveness of control measures against the diffusion of the COVID-19 pandemic is grounded on the assumption that people are prepared and disposed to cooperate. From a strategic decision point of view, cooperation is the unreachable strategy of the Prisoner’s Dilemma game, where the temptation to exploit the others and the fear of being betrayed by them drives the people’s behavior, which eventually results in a fully defective outcome. In this work, we integrate a standard epidemic model with the replicator equation of evolutionary games in order to study the interplay between the infection spreading and the propensity of people to be cooperative under the pressure of the epidemic. The developed model shows high performance in fitting real measurements of infected, recovered and dead people during the whole period of COVID-19 epidemic spread, from March 2020 to September 2021 in Italy. The estimated parameters related to cooperation result to be significantly correlated with vaccination and screening data, thus validating the model. The stability analysis of the multiple steady states present in the proposed model highlights the possibility to tune fundamental control parameters to dramatically reduce the number of potential dead people with respect to the non-controlled case.


Introduction
The recent coronavirus pandemic enforced the application of control measures for restraining the virus diffusion. The most successful measures for weakening and reducing the infection aim at preventing people from getting infected or at least to get the disease in a weak form, and consist, on the one hand, of the adoption of medical treatment protocols and of vaccination, and, on the other hand, of avoiding behaviors which favour virus spreading, such as social distancing, movement reduction, mask wearing, and so on [1]. The application of any control measures requires the organization of suitable information campaigns aimed at inducing people to adopt correct behaviors against the pandemic. To make these campaigns effective, people must behave cooperatively with respect to the limitations and requirements imposed by the deputed institutions and governments. Unfortunately, more often, when taking decisions under strong pressure, such as, for example, in the initial phases of the pandemic, the requirement of strong efforts may activate in the population selfish and conservative mechanisms, which often lead to the temptation to exploit the cooperative behavior of the others and to the fear of being betrayed by their incorrect conduct. An example of these mechanisms is represented by the panic buying that arose at the beginning of the COVID-19 pandemic (see, for example, [2]). Similarly, correct campaigns should be able to convince citizens that the costs of vaccination are less than the rewards, resulting in the possibility to recover almost normal life styles and well-being. Avoiding selfish mechanisms and promoting good ones requires fostering altruistic and cooperative feelings.
Emerging of cooperation is recognized as challenging in many different fields [3,4], and Evolutionary Game Theory (EGT) represents a natural mathematical framework to model this problem. Indeed, it provides a rigorous methodology for studying strategic interactions among people evolving over time [5,6]. The evident drawbacks of selfish behavior are highlighted in the defective Prisoner's Dilemma game, in both the mean field formulation [7,8], and for interconnected populations [9,10]. The influence of social networks on the dynamics of cooperation in evolutionary games has been deeply investigated [11][12][13][14][15].
With specific reference to the people's behavior in the management of an epidemic, an extensive review of the literature on behavioral epidemiology is presented in [16], where the possibility to include its effect in mathematical models is discussed from a methodological point of view, pointing out that concrete models should account for many factors, thus significantly increasing the number of parameters and equations. Models grounded on game theory have been proposed for analysing and predicting vaccination behavior [17,18]. In [19], a suitable game is introduced to account for the reward of social distancing in vulnerable and non-vulnerable agents, while [20] studies the optimal strategies of people to contain the epidemic spreading by appropriately planning their movements to avoid strongly infected places.
Moreover, the emergence of the COVID-19 epidemic fostered the development of mathematical models for the integration of the standard [21] or adapted [22,23] SIR and SIRS models with suitable control measures for contrasting the infection in different ways. For example, in [24] the role of a network on disease spreading is analyzed, in [25] optimal control measures are identified by assuming networked populations, while in [26] an analysis of strategic behavior during the COVID-19 pandemic was developed from an economic perspective. Additionally, in [27] the strong influence of collective behavioral patterns of the population in epidemic spreading was highlighted by stochastic models.
In this paper, we propose the integration of a SIRS epidemic model with a replicator equation, describing the evolution over time of cooperation in large populations. Despite the above large literature taking reasonably into account the presence of a network and its role in the disease spreading, in this paper we focus on the interplay between epidemiological and social dynamics from a mean field point of view, since our main objective is to find simple mechanisms regulating people's behavior and cooperation during a pandemic event in a large and complex nation like Italy. Similar approaches have been taken into account in [28], and deeply studied in [29], where the relationship between control measures, vaccine hesitation, information retrieval, and policy making has been investigated to enforce cooperative behavior aimed to reduce the virus' spread.
A convincing motivation for developing models like the one proposed in this paper is provided in [30], where the authors highlighted that the effects of shield immunity and economic shutdowns are complementary, such that governments should pursue them in tandem. In a more theoretical setting [31], the SIR model has been coupled with the replicator equation of evolutionary games. Similarly to our study, people's behavior is influenced by the infection strength, and, in particular, by evaluating the infection rates of cooperators and defectors, including the corresponding risks and costs. Since evaluating the infection rates is a difficult task for people, the model proposed in this paper assumes that cooperative/defective decisions are driven by easily retrievable information, such as the number of dead or infected individuals. A similar reasoning for spontaneous adoption of quarantine states has been done in [32]. Developing this new model may help to formalize and understand better the twofold feedback mechanisms between epidemic spreading and people's behavior from a theoretical point of view. On one hand, the infection rate is assumed to depend on the propensity to cooperate by respecting the control and recovery measures taken by the governments. In our case, it decreases when cooperation increases. On the other hand, the parameters of the game payoff matrix depend on the strength of the epidemic at any time; thus, the higher is the gravity of the epidemic, the higher is the propensity to cooperate, despite the limitations and costs imposed by the government's directives. Overall, the switches from the Prisoner's Dilemma game to the fully cooperative Harmony game and vice versa are observed.
The remaining of this paper is organized as follows. In the Material and Methods section, the assumptions on which the mathematical model is grounded, together with its theoretical properties, are reported, while the Results section presents the main findings of the research. Finally, some concluding remarks are provided in the Conclusions section.

Materials and Methods
A population of M individuals, split into 6 classes, is considered: According to the literature on epidemic models [21], and to significant recent results [25,33], we start from the following epidemic model of COVID-19: where β is the infection rate, α is the drop immunity rate, γ is the spontaneous healing rate, τ is the detection rate, ω is the healing rate of detected infected persons, and ζ is the death rate. According to [25], we assume that, for the COVID-19 case, some model parameters and the reproduction number are time varying. For the sake of simplicity, in the rest of the paper, we will refer to (1) as Infection Model (IM).

How Epidemics Influence People's Behavior
In an epidemic context, it is reasonable that individuals behave by choosing one of the following two strategies: cooperating (C) or not (N). In our study, the first option consists of respecting the restrictions imposed by the deputed institutions for constraining the spread of the disease, while the second indicates people adopting habits not fully aligned with the sanitary guidelines. We denote by x ∈ [0, 1] the share of cooperators, and by z = 1 − x the share of non cooperators.
In Game Theory, the decision mechanism is ruled by a payoff matrix, describing the outcome of the interaction between couples of players randomly chosen in the population [34]. In the investigated context, the payoff matrix has been chosen following [24]: where, the first row contains the payoffs for cooperation (C) and the second one the reward for defection/non-cooperation (N) for player 1. Similarly, the columns of the payoff matrix are related to the choices of player 2. More specifically, T corresponds to the temptation to defect and S is the "sucker's payoff", embodying the fear of being betrayed by others. Player 1 earns 1 if both cooperate, T if he defects and the opponent cooperates, S if he cooperates and the opponent defects, or 0 if both players defect.
It is useful to introduce the net payoffs for cooperation, σ C = 1 − T , and for defection, σ N = −S. Then, for σ C < 0 (i.e., the temptation to defect is higher than the reward for mutual cooperation) and σ N > 0 (i.e., the fear of being betrayed is significantly high), the payoff matrix (3) embodies a Prisoner's Dilemma game, where defection is dominant, while for σ C > 0 and σ N < 0, the game switches to a Harmony game, and cooperation is the dominant strategy [34].
In the epidemic context, the Prisoner's dilemma game can be restated as follows: if two individuals meet, and adopt good practices (i.e., social distancing, mask wearing, getting vaccinated), then they both get a normalized reward equal to 1, corresponding to the preservation of their own health. In contrast, if two individuals meet, and, for instance, neither of them wear a mask, then there is a high probability that, if one is infected, the second will get infected too. In this case, both players earn a payoff equal to 0, i.e., their health has not been preserved. Now, suppose that one of them knows that the other is respecting the containment measure; then (s)he will be tempted not to respect the good practices. In this case, the "free rider" will earn T > 1, which corresponds to her/his own safety (the other will not infect her/him, and hence her/his payoff is equal to 1) plus, for instance, the feeling of freedom of staying without mask (payoff equal to T − 1 = −σ C > 0). On the other hand, the "good citizen" is not protected-the other can be infected and can infect her/him too. Thus, s(he) will earn 0 minus the feeling of being the only one respecting the sanitary guidelines (payoff equal to S = −σ N < 0). When we move to the Harmony scheme, the outcomes for the hybrid case (a cooperator meets a defector) change, since T < 1 and S > 0.
Given a generic game, the dynamical evolution of the strategies within a well mixed population can be described by the Replicator Equation (RE) [5,6,34], which reads as follows: where represent the average payoffs π C and π N collected by the share of cooperative and noncooperative individuals, respectively, while is the average payoff of the whole population. When π C is bigger than π, the time derivative of x is positive, and hence x increases over time. At the same time, by the definition of π it follows that if π C > π, then π N < π, and hence y decreases over time. Symmetric arguments hold for the opposite situation (π C < π). Coherently with the payoff matrix (3) and Equation (4), cooperation evolves towards full defection (x = 0 and z = 1) for the Prisoner's Dilemma case or towards full cooperation (x = 1 and z = 0) for the Harmony game. It is natural to assume that the ruling parameters T and S vary with respect to people's perception on the infection strength, reasonably measured by the state variables I and D of the IM model (U is not observable, since it corresponds to the number of undetected infected people). In this regard, we can assume that the payoff matrix (3) depends on I and D, i.e., B = B(I, D), and that it changes from a Prisoner's Dilemma game to a Harmony one when the value of I + D exceeds a given activation threshold a: and where works as a "game-switching" term, σ 0 C > 0, σ 0 N < 0 and 0 < a < M. The term G(I, D) embodies the concept of coevolutionary rules in the evolutionary games framework, thus accounting at the same time for the evolution of strategies, and other related variables, such as the pandemic status [13]. Additionally, the role of G(I, D) in the proposed model recalls the findings reported in [15,35] on the ability of finite and networked populations to self-regulate the behavior when a higher common good (in this case, the public health) should be preserved, thanks to the awareness of individuals.
For I + D < a, B(I, D) represents a Prisoner's Dilemma game, since σ C (I, D) < 0 and σ N (I, D) > 0. That is, when the epidemic is enough low, defection is the "natural" choice, since a perceived weakness of the disease does not justify the adoption of any containment measure. This is coherent with the fact that some governments have been laxer than others in imposing restrictions. In this sense, the threshold a corresponds to an important tuning parameter for the policy-makers, as it can be reduced by means of effective informative campaigns on the real risks of the current sanitary crisis, thus making I + D > a. In this case, B(I, D) turns to a Harmony game since σ C (I, D) > 0 and σ N (I, D) < 0.
For the sake of simplicity, we observe that, by exploiting z = 1 − x, and using Equations (5) and (6), we get: Then, the RE (4) is reduced to a monodimensional ODE [34] which reads as:

How People Behavior Influences the Epidemic
As mentioned before, cooperation corresponds to the adoption of good practices, aimed at reducing the disease spreading, such as, limitation of social interactions, wearing masks, getting vaccinated, and so on. As a consequence, the infection rate β of IM varies according to people's behavior. More specifically, cooperation produces a reduction of the infection rate, while defection leads to its increment. In this work, we assume that the infection rate β linearly depends on the cooperation x as follows: where β 0 is the natural infection rate of the disease, while e ∈ [0, 1] weights the influence efficacy of the cooperation on the infection rate. Notice that β(x) ≥ 0 ∀x ∈ [0, 1]. In this case, reproduction number is defined as follows: where represents the basic reproduction number of the disease. Since x is a function of time, then R(x) is time varying, and it reaches the maximum R 0 for no cooperation (x = 0). Moreover, its minimum is equal to and it is attained in presence of maximum cooperation (x = 1), i.e., when countermeasures against the disease are effective at the highest level.

The Infection Game Model
By coupling the standard IM (1) with the RE (8), we get the Infection Game Model (see Figure 1), hereafter called IGM. The model is characterized by the extended state vector y = [S, U, I, R, D, x] , and it reads as: where β(x) and G(I, D) are defined in Equations (9) and (7), respectively. It can be easily proven that all the trajectories of system (13) are confined in the feasible set for all times t > 0, provided that the initial condition y(0) belongs to the same set (see [33] for S, U, I, R and D, and [34] for x). Model (13) presents infinite equilibria in F , specifically: It is worthwhile noticing that the introduction of a new equation into the epidemic model does not alter the biological mechanisms on which it is grounded. Indeed, at steady state, U * , I * and R * are all null, and the lines y * 1 and y * 2 correspond to the one of model (1). The line y * 3 is introduced by the replicator Equation (8), and, differently from y * 1 and y * 2 , where the D * component is equal to k 1 or k 2 , and it depends on the initial conditions, it provides the policy-makers with the exact information on the total amount of dead people at the equilibrium since D * = a.
For investigating the dynamics of the proposed system, we linearize (13) nearby each equilibrium and the eigenvalues of the corresponding Jacobian matrices ∂ẏ ∂y are evaluated [36]. For y * For y * 2 (k 2 ), the eigenvalues are Finally, the eigenvalues for y * Following the theory of dynamical systems [37], it is well known that the Jacobian matrix of a n-dimensional linearized dynamical system, evaluated in a steady state y * , has n c eigenvalues with null real part, n u eigenvalues with positive real part and n s eigenvalues with negative real part, where n c + n u + n s = n. The corresponding eigenvectors span the center, E c (y * ), unstable, E u (y * ), and stable, E s (y * ), subspaces, respectively.
Since the Jacobian matrices of all equilibria of system (13) have at least two null eigenvalues, the corresponding center subspaces are not empty. This means that no steady state is hyperbolic for any values of the parameters. In this case, it is interesting to investigate the conditions for which the unstable subspaces are empty. The following propositions hold. Proposition 1. E u (y * 1 (k 1 )) = ∅ if and only if a > k 1 and k 1 ∈ (k 1 , a).
is negative if and only if k 2 > a. Moreover, since γ > 0 and Then, E u (y * 2 (k 2 )) = ∅ if and only if a > k 2 and k 2 ∈ (a, M), or a < k 2 and k 2 ∈ (k 2 , M).
Remark 3. a < k 1 guarantees that h > 0. Indeed, using Equation (14), we get M R 0 = M − k 1 , and then: Summing up, for a > k 1 then h < 0, and hence E u (y * 3 (h)) = ∅ for all the equilibria y * By analyzing the case where R 0 > 1 (i.e., the disease is sufficiently strong to diffuse in the population), we observe that, in agreement to Propositions 1-3, the plane (e, a) is subdivided into three regions: R1 (a > k 1 ), R2 (k 2 < a < k 1 ) and R3 (a < k 1 , a > k 2 ). In 1] are reachable. In Figure 2, the three regions are

Available Data
In this work, we use data collected in Italy by Protezione Civile (the national Italian agency dealing with the prediction, prevention and management of emergency events), updated daily on GitHub [38]. In particular, the data used in this paper correspond to the infected, recovered and dead people per each day in the period of observation from 24 March 2020 to 14 September 2021 (540 days). The initial amount of infected, recovered and dead is taken equal to the initial value of the corresponding variable in the dataset. Since only the sum of the initial conditions of susceptible and undetected infected is known, U(0) or S(0) will be estimated and the other one will be calculated by difference from the condition S(0) + U(0) + I(0) + R(0) + D(0) = M.

Parameters Estimation
The proposed model is characterized by three different sets of parameters. The first set (θ 1 ) includes parameters that are constant for the model but are assumed unknown, and hence they are estimated from data. The unknown initial condition of the effective susceptible individuals S(0) is also estimated together with the initial condition of the behavioral variable x(0). The second set (θ 2 ) is composed by unknown estimated timevarying parameters, allowing us to take into account the epidemic evolution over time windows of size T W = 30 days. Since the dataset is composed of 540 days, there are in total W = 18 time windows. The last set (θ 3 ) consists of parameters fixed according to the literature values, since these are related to specific biological properties of COVID-19.
The estimation problem is formalized as follows: whereÎ,R andD represent the real measurements of infected, recovered and dead individuals, respectively. The estimation problem (17) has been solved for both IM and IGM. For IM, θ 1 = {S(0)}, θ 2 = {β, α, ω, ζ} and θ 3 = {γ, τ} = 1 14 , 0.166 , according to [25]. For IGM, , ω, ζ, e, a} and θ 3 is the same as in IM. Notice that, since the variable e and x are absent in the IM model, then Equation (9) reduces to: and β 0 is assumed time varying in the estimation procedure (β = β 0 ∈ θ 2 ). The larger number of parameters of IGM is due to presence of the game dynamics. Suitable constraints for each parameter in θ 1 and θ 2 have been set to guarantee their physical and/or mathematical meaning. The specific ranges of the parameters, the initial conditions and the estimated values are reported in Table 1 in the Section 3. The same constraints are used in both IM and IGM. Table 1. Values of the estimated parameters in the set θ 1 and in the θ 2 , values of parameters in the set θ 3 , ranges and initial conditions used for the estimation. Moreover, t.v. stands for "time varying", while n.a. means that the parameter is not estimated in the corresponding model.

Range
Param The Matlab source code for the estimation and simulation of both IM and IGM described in this work is available at https://github.com/dariomadeo/IGM (created on 20 December 2021).

Results
In this Section, the results of the parameter estimation procedure are reported. Figure 3 shows the recorded dataÎ,R andD (yellow dots), and the simulations of the three state variables I, R and D of IM (red lines) and IGM (blue lines) obtained using the estimated parameters.
It can be seen that the fitting performances between simulations of both IM and IGM model and real data, evaluated in terms of average percentage errors over times of the three measured state variables, defined as: where Ξ ∈ {I, R, D}, are satisfactory (see Table 2).   Figure 4 shows the dynamics of the undetected infected individuals U of IGM (blue line, left ordinate scale), as well as the percentage of undetected infected individuals with respect to the total infected population U + I (violet line, right ordinate scale) over the whole time period. The percentage of undetected cases over the total is about 20% in average, with a peak of 40%, indicating that the number of unrecorded cases is high, and hence the adoption of more organized containment measures, such as screenings and vaccination campaigns are still desirable to prevent future disease resurgences. The fixed parameters estimated by solving the optimization problem (17) over W · T w = 540 days of recordings of the variables I, R and D, are reported in Table 1, which includes also the list of all the model parameters, constant, fixed, and time varying, and their meaning. The evolution of α, ω, and ζ, is drawn in Figure 5, where the red lines are related to IM, while the blue ones to IGM. Subplot (a) contains the dynamics of the reinfection rate α. In subplot (b), the time course of the recovery rate ω is drawn. In subplot (c), the dynamic of the death rate ζ is depicted, which decreases over time, thanks to the improvement of the caregivers' techniques and to the vaccination campaign. In all cases, we observe a perfect agreement of the parameter values estimated for the two models. , ω (subplot (b)) and ζ (subplot (c)). The piecewise shape of the curves is due to the fact that these parameters are constant over each time window of T w = 30 days.
Furthermore, subplot (a) of Figure 6 shows the dynamics of the control parameter a (blue line), compared with the term G(I, D) (violet line), which rules the game switching. The fluctuations of a over time indicate that this parameter embodies the sensitivity of people to the observation of the pandemic status. High values of a indicate that people are less sensitive to the pandemic status, thus in order to be effectively cooperative they need to know that the number of infected and dead people is very large. The opposite behavior is observed for G(I, D), which embodies the two-sided game. Indeed, people play the Prisoner's Dilemma game when a is bigger than the pandemic strength, or a Harmony game in the opposite case. In this sense, G(I, D) acts as an "information filter", balancing the effective information represented by I + D with people perception a of the pandemic. In this regard, it would be interesting to further investigate the quality and effectiveness of the information campaign performed by media and politics. Our results show that, people played Prisoner's dilemma from March to October 2020, then switched to Harmony game from October 2020 to June 2021 (around day 450), and then again to PD. In this case, the combined effect of spring and successful vaccination campaign reduced the infection entity and allowed people to relax the attention level. Finally, during fall 2021, the game was slightly turning back to Harmony. Subplot (b) of Figure 6 depicts the time course of the control parameter e (blue line), of the state variable x(t), corresponding to the average people cooperation (violet) and the combined effect of e and x, which represents the overall efficacy of the containment measures (green line). It is interesting noticing that the overall cooperation is in general increasing from fall 2020 to the end of the observation. This fact embodies a natural "learning process", leading people to the adoption of good practices even though the government impositions are almost steady. Moreover, recalling the theoretical results depicted in Figure 2, one can see that all estimated values of parameter a fall below k 1 = 16.6 × 10 6 , thus the worst case concerning a potential dramatic number of dead people was avoided along the whole pandemic. For low values of a, the parameter e is most effective in terms of epidemic containment (a deep discussion on this aspect is reported below). Finally, subplot (c) displays the evolution of β for IM (red line) and IGM (blue line). The two curves exhibit very similar time courses, meaning that, differently form the usual SIRS models, such as the IM considered in this study, the proposed IGM explains the pandemic trends through new mathematical terms, accounting for the people ability to cooperate with the aim of reducing the infection.
Furthermore, the trend of β is coherent with the values of the time varying reproduction number R(t), as drawn in Figure 7, where the evolution of R(t), evaluated according to different methods, is shown. The reproduction number calculated on real data (yellow circles) according to the following formula [39]: where γ = 1 14 and τ = 0.166 are assumed constant, according to [25]. We notice that the method implemented for calculating R(t) is in agreement with the one that supported the decisions of Italian Government along the whole COVID-19 pandemic [40]. The parameters R(t), R 0 and R have been also evaluated using Equations (10)- (12), and depicted in Figure 7 with blue, violet and green lines, respectively. According to the theory of dynamical systems, we notice that the trend of R(t), evaluated on the basis of the estimated parameters of the IGM, anticipates the R(t) calculated using data, thus confirming the good performances of the proposed model. The violet area represents the distance between the current situation (blue line) and the worst case (violet line), i.e., R(t) = R 0 , corresponding to the equilibrium state x = 0. Similarly, the green area corresponds to the distance between the current situation and the best case (green line), i.e., R(t) = R, corresponding to the equilibrium state x = 1. Moreover, in order R to be minimum, the parameter e should be as large as possible. Then, at each time window, the best containment of the epidemic spread corresponds to have the blue curve coinciding to the green one, both lower than 1.
Upon these theoretical considerations, jointly with the very good performances of the proposed model, we can use parameter e to evaluate how successful have been the containment measures. To this aim, Figure 8 depicts the percentage variation over two consecutive months between the maximum values of infected (yellow) and dead (violet) individuals (y-axis), obtained with the value of e (x-axis) set during the first month, as indicated by the label. A visual inspection of the Figure reveals that the best results have been achieved for e sufficiently high (see also Figure 2 as reference), such as during the second lockdown (month 9, November and December 2021). Moreover, at the end of the first lockdown (month 2, May and June 2020) the value of e was very high, coherently to the heavy efforts asked to people in that period. On the other hand, we highlight the fact that the second, third and fourth waves, representing the periods in which the infection reached the highest peaks in both absolute and relative terms, resulted from the use of low values of the parameter e. These results are highlighted by circles in Figure 8, where, for example, the periods indicated by labels 2 and 9 clearly belong to region R3 of Figure 2. The impact of the control parameter e on the pandemic has been inspected in two different periods. The first period starts on day 210 (beginning of the second wave), and the parameter e is assumed varying in the next time window between 0.1 and 0.6, around the estimated value e = 0.29. Subplots (a) and (b) of Figure 9 display the resulting infected and dead individuals, respectively. The second period starts on day 510, where in the next two time windows, the parameter e varies in the range [0.3, 0.6], which includes the value e = 0.34 estimated in the first window (from day 510 to day 539). Subplots (c) and (d) of Figure 9 show the simulated infected and dead individuals, respectively. In each subplot, yellow dots indicate real data, the red curve is the simulation of the IGM using the estimated parameters, and the light blue area represents the range of variations of the dynamics for different values of e. It is clear that reducing the value of e produces an increase of both infected and dead individuals. On the other hand, when e is fostered, positive effects on the pandemic status are observed. Overall, this analysis suggests that, even though a strong increase of the parameter e can heavily drop the current pandemic status, its values should be chosen in a wise way, in order not to stress too much the population. Indeed, very good results are achievable even by setting e to value significantly lower than its maximum, i.e., e < 1. Figure 9. Two different scenarios obtained by varying the control parameter e at the beginning of the second wave (subplots (a,b)), reporting infected and dead people, respectively) and at the beginning of the fourth (subplots (c,d)), depicting infected and dead people, respectively). Yellow dots indicate real data, red curve is the simulation of the IGM using the estimated parameters, and the light blue area corresponds to the range of variations of the dynamics for different values of e.
In order to deeply investigate the meaningfulness of the developed model, we recall that the time varying parameters common to IM and IGM, i.e., α, ω, ζ and β, have similar shapes (see Figures 5 and 6). As a consequence, it is interesting to assess the significance of the specific IGM state variables and parameters, such as x, e and a. To this aim, we report the correlation analysis between these parameters and additional data, available in [38,41], not used for the modeling and identification phases. In particular, we have chosen data related to people's decisions on vaccination and testing. Figure 10 reports the scatter plots of the number of daily swabs (subplots (a-g)), first vaccine (subplots (h-n)) and second vaccine (subplots (o-u)) administrations (y-axis) versus α, ω, ζ, β, x, e and a (x-axis). Figure 10. Scatter plots of the number of daily swabs (subplots (a-g)), 1st vaccine (subplots (h-n)) and 2nd vaccine (subplots (o-u)) dose administrations (y-axis) versus α, ω, ζ, β, x, e and a (x-axis).
The colors indicate the time flow (blue for the first days, yellow for the last). The corresponding regression lines are reported (dashed black), as well as the correlation coefficients ρ. Most significant correlations (|ρ| > 0.5) are highlighted in red.
In the first four columns, the time varying parameters α, ω, ζ and β are considered. A correlation of −0.51 between the drop infection rate α and the second vaccine administration (subplot (o)), and a correlation of 0.5 between the recovery rate ω and the first vaccine administration (subplot (i)) are found. The death rate ζ has a negative correlation of −0.73 with the second vaccine administration (subplot (q)), coherently with the fact that vaccine reduces dramatically the mortality. Finally, the infection rate β exhibits, as expected, a negative correlation of −0.61 with the first vaccine administration (subplot (k)).
The peculiar structure of the β in the IGM (see Equation (9)) allows one to take into account the social dynamics as additional detail, and hence it can be used to better understand the dynamics of vaccination and testing, as well as to plan more appropriate countermeasures. x exhibits high correlations with swabs (subplot (e), ρ = 0.62), 1st vaccine (subplot (l), ρ = 0.75) and 2nd vaccine (subplot (s), ρ = 0.61). Indeed, a strong cooperation implies a diffuse willingness of people to check and protect their own health. Furthermore, significant positive correlations are found between the parameter e and 1st vaccine administration (subplot (m), ρ = 0.53), and between a and 2nd vaccine administration (subplot (u), ρ = 0.63). In the fist case, the diffusion of cooperative behavior stimulated an increase of the 1st vaccination in the population. The last case implies that an increased tolerance to the epidemic status may produce a higher propensity to vaccination. These results confirm that the extended IGM is able to shade light on additional aspects related to the social dynamics during the pandemic.

Conclusions
This paper presents an extended model, called IGM, which integrates a standard epidemic system and evolutionary games, used to describe the propensity of people and institutions to be cooperative, thus adopting the measures taken to control the epidemic. The resulting model assumes that the infection rate decreases when the cooperation rises, by means of a control parameter e, while the propensity to cooperate is favoured by strong infection levels, ruled by the threshold parameter a, which embodies the awareness of people concerning the gravity of the current sanitary crisis.
A stability analysis of the IGM equilibria was carried out, highlighting the possibility to tune the parameters e and a in order to drive the system dynamics towards regions where the number of potential dead people is dramatically reduced with respect to the non-controlled case. An identification procedure, consisting of the minimization of the relative square difference between real measurements and model simulations to estimate system parameters, has been applied. The results indicate high fitting performances of the developed model. Furthermore, the model, which assumes the presence of a twofold feedback between the infection strength and the propensity of people to comply with the containment measures, turns out to be effective in producing the observed alternating periods of rising and damping of COVID-19 epidemic in Italy. The satisfactory results are also confirmed by the strong correlations found between the behavioral variable x, indicating the cooperation level, and swabs, first and second vaccine administrations.
Overall, thanks to its theoretical properties and to its good performance in terms of fitting, the proposed model aspires to be considered as a reference framework to account simultaneously for the biological aspects of an epidemic and for the social and decision issues related to its containment. Future works will concern the extension of the present model to include the vaccine dynamics, making it even more effective to be adopted by the policy-makers. As a further development, the authors are investigating the impact of "imperfect information" on the described process. Indeed, it is well known that some of the infectious individuals are pauci-or even a-symptomatic. The information provided on the status of the pandemics can be also incorrect due to the difficulties inherent in data collection, or due to the diffusion of fake news. As a natural consequence, it would be interesting to analyze the system behavior by considering only the fraction of (supposedly) known infectious individuals.

Data Availability Statement:
The data used for this study are openly available in GitHub, reference numbers [38,41].

Conflicts of Interest:
The authors declare no conflict of interest.