Cooperation dynamics under pandemic risks and heterogeneous economic interdependence

The spread of COVID-19 and ensuing containment measures have accentuated the profound interdependence among nations or regions. This has been particularly evident in tourism, one of the sectors most affected by uncoordinated mobility restrictions. The impact of this interdependence on the tendency to adopt less or more restrictive measures is hard to evaluate, more so if diversity in economic exposures to citizens’ mobility are considered. Here, we address this problem by developing an analytical and computational game-theoretical model encompassing the conflicts arising from the need to control the economic effects of global risks, such as in the COVID-19 pandemic. The model includes the individual costs derived from severe restrictions imposed by governments, including the resulting economic interdependence among all the parties involved in the game. By using tourism-based data, the model is enriched with actual heterogeneous income losses, such that every player has a different economic cost when applying restrictions. We show that economic interdependence enhances cooperation because of the decline in the expected payoffs by free-riding parties (i.e., those neglecting the application of mobility restrictions). Furthermore, we show (analytically and through numerical simulations) that these cross-exposures can transform the nature of the cooperation dilemma each region or country faces, modifying the position of the fixed points and the size of the basins of attraction that characterize this class of games. Finally, our results suggest that heterogeneity among regions may be used to leverage the impact of intervention policies by ensuring an agreement among the most relevant initial set of cooperators.


Introduction
The COVID-19 pandemic has caused the most severe global economic breakdown of the recent history, mainly due to the measures to control the virus spread (e.g., quarantines, stay-at-home and social distancing policies), which have produced a dramatic shutdown of the economic activity. In this context, the coordination among countries is an essential instrument to efficiently control pandemic and enhance economic recovery [1,2] .
The pandemic control measures result in two types of negative economic effects on countries and regions. First, the direct effects, derived from internal mobility restrictions, which seriously injure many sectors in the country. They do not just include hospital-ity and entertainment, but also banks, stock market and education sector [3] . Second, indirect effects, induced by trade and financial linkages among countries, which produce economic contagion of the consequences of restrictions in a country to a foreign linked country [3,4] . An example of an indirect effect is the travel restrictions in tourist dependent countries, which produce a serious economic loss to destination countries [5] . Other pandemic consequences are found in international supply chains. Thus, exportoriented countries are influenced by demand shortfall of the importing countries and, at the same time, supply disruption from some countries such as China produces input shortages and inflationary pressure in import-countries [6] . In addition, the contact pattern in some regions are influenced by the social distancing policies followed in others [2] . All this phenomena stem from the economic interdependence (EI) among countries.
Furthermore, the economic impact of the pandemic control measures is not homogeneous for all the countries and regions. Those service and export-oriented countries, such as those depending on tourism, entertainment and transport, have been more affected than others [7] . For example, Panama, which depends mainly on tourism and transport services, suffered a serious GDP decrease in 2020 (19%), and other many tourism-dependent smallislands states, such as Maldivas, Fiji and Bahamas, experimented similar economic consequences. In Europe, the COVID-19 pandemic economic loss in Spain, around 10% decrease of the GDP, duplicated the economic loss in Germany. These evidences show that the economic consequences of the implementation of control measures are markedly heterogeneous.
The coordination problem of pandemic control measures can be represented through a collective-risk dilemma (CRD) [8,9] . This is a multiplayer public good game (PGG) where every player can contribute with some amount to avoid a certain risk of failure. Normally, it is not necessary that all players contribute to achieve the common goal and there is some space for free-riders who benefit from the others' contribution. Traditionally, cooperation to mitigate climate change by cutting carbon emissions has been one of the most important applications of CRDs [10][11][12][13][14][15][16] . More recently, CRDs have been proposed for the coordination of restrictions and reopen policies in the current context of the COVID-19 pandemic [17] .
The aim of our study is to analyze the conditions for achieving coordination among countries to control pandemic situations, taking into account the specific economic consequences of these measures. Specifically, we propose a new CRD model where EIs among the players are considered by altering the expected profits by regions given the restrictions applied by cooperating regions. The proposed CRD model also incorporates the heterogeneity found for regions and countries with respect to their income loss when cooperating by restricting the economic activity. Heterogeneity is omnipresent in reality but only few studies included heterogeneous features in social dilemmas [18] . In fact, recent studies have shown that heterogeneity of the social network is less relevant than the levels of cooperation in the neighborhood of the players [19] and this is where the heterogeneous features of the players can make a difference. Additionally, public goods might be better off when managed from the bottom up [20] and including the heterogeneity of the players is clearly relevant for this venue.
We fit a log-normal distribution for feeding the model with heterogeneous economic loss values by analyzing real data from tourism contribution to GDP of the European Union NUTS2 regions. The experiments comprise the evaluation of the final cooperation levels of the population with and without EI and validate them under the presence of heterogeneity. Our results show that the existence of EI among countries can favor cooperation for all the tested conditions. To understand the reason behind these observations, we extract the stable and unstable points of the new dilemma with and without heterogeneity. To this end, we propose a new way of calculating the internal roots from the agent-based simulation results.
Finally, we take advantage of the reality coming from the income loss heterogeneity of the model to analyze those initial conditions facilitating cooperation. Three scenarios are evaluated. We first set the initial cooperators at random but we also condition those initial cooperators by a positive and negative bias by they income losses of the regions. The reader will see how significant differences in the final cooperation levels are achieved when choosing the best initial conditions alternative. These results will help to engineer more effective governmental policies to increase cooperation.

Model
We present a CRD model to represent the cooperation game among regions or countries (we call them countries from now on) when adopting measures to control pandemic spread having into account their negative economic consequences. We resort to evolutionary game theory and stochastic population dynamics [21][22][23][24][25][26][27] when required, combined with agent-based computer simulations [28,29] . The inspiration for this model was found in previous CRD models for cooperative actions against climate change [9,14] and pandemic spread [17] .

Game definition including economic interdependence (EI)
The model includes a finite number Z of countries or regions.
Every player i can choose two strategies s i (t) at every time step t: cooperation ( C), which means adopting measures to control the pandemic spread and suffering income loss; and defection ( D ), which means not adopting any measure, continuing with the normal economic activity and eventually free-riding the correct public health conditions derived from others' cooperation.
Players interact in groups of size N, representing international or regional agreements, alliances or work-groups. By assumption, these groups are randomly formed. Every group faces a risk r ∈ [0 , 1] for the health care system to collapse and the resulting economic breakdown, if the epidemic is not controlled enough inside the group. This happens when the number of cooperators in the group does not achieve a minimum M ≤ N. When the number of cooperators is equal or higher than M in a group of size N countries, the pandemic is under control and the economic breakdown is not produced.
When a country i cooperates, economic activity critically stops and an income loss c i occur in the region ought to these restrictions. This loss is not necessarily constant throughout the population, but every country has its loss value, which depends on the specific economic structure of the country. In addition, the economic breakdown in the cooperator i is spread to the rest of countries due to the international economic linkages among them. Thus, defector countries are not free of economic negative spillovers or contagion from other countries. Instead, they are influenced by the number of cooperators in the total population. The more the cooperator countries in the total population, the higher the negative economic consequences in the defector country.
The conditions above are represented in the expected payoff i D ( i C ) of a defector (cooperator) country i . They are: where j is the number of cooperators in the group and k is the number of cooperators in the total population. The variable M = mN , where N is the group size and m is the minimum fraction of cooperators to avoid collapse. The Heaviside step function (x ) is equal to 0 whenever x < 0 and equal to 1 otherwise. The initial endowment or maximum payoff obtained in absence of any pandemic is normalized to 1. The first term of Eq. (1) shows that when the number of cooperators in the group is above the threshold M, the expected payoff of a defector is not maximum, but lowers an amount due to the economic interdependence with cooperator countries. Specifically, the loss of the defector's income is a fraction k Z of its own maximum loss c i . This term c i k Z is added to the cooperator's payoff ( Eq. (2) ) since we assume that the single cost of cooperation c i already includes the economic contagion derived from international linkages. The total loss in the cooperators' income is c i , as expected. When the necessary number of cooperators in a group is not achieved ( j < M), both cooperators and defectors have a risk r of global economic collapse and null economic activity.

Individual update of the game strategies
After playing a game round t, players can update their strategies according to the received payoffs. Here, we have considered the Fermi function as the evolutionary update rule [30,31] . The Fermi rule is a stochastic pairwise comparison rule in which strategies that do well, are more likely to be imitated, and spread throughout the population. In detail, at each time-step, a player i with a payoff i is randomly selected from the population for strategy revision. Player i will then randomly select another player j from the population as a potential role model; i will imitate the strategy of j with a probability p that increases with their payoff difference -( j − i ) -and can be written as in [24] : The free parameter β is the intensity of selection, encoding the chance of mistakes during the imitation process. This means that a player i can copy another player's strategy j despite having a lower payoff. We set β = 1 in all the experiments of the study. Additionally, the players of the game can randomly explore other strategies, adopting a strategy at random with probability μ. This mutation (or exploration) probability ( μ) equals 1 Z in all experiments. This update rule can be used in both synchronous and asynchronous paradigm. We have confirmed via numerical simulations that our conclusions remain valid for both cases.

Stochastic population dynamics
Let us start by assuming that the income loss c i is homogeneous for all players. Then we have a unique income loss c for all countries and the payoff derived from any strategy depends on the general parameters and the number of cooperators j in the group. Let us also assume that all players are equally likely to interact, a configuration known as a well-mixed population [22] . In other words, we have a random selection of partners to form groups. In this limit, we can analytically compute the expected payoff (or fitness) of a cooperator/defector for a given number of cooperators k in the population by using an hyper-geometric sampling [32,33] : where we have removed the super-index referring to player i in the payoff functions.
The update rule described in Section 2.2 defines a Markov process where, at every time step, the probability to increase the number of cooperators k in the population is [24] , and the probability to decrease the number of cooperators is Several tools can be used to analyze the evolutionary dynamics emerging from this ergodic Markov chain. First, the gradient of selection [9,33] , G (k ) = T + (k ) − T − (k ) , indicates the direction of change for every cooperation level. Second, as the process includes probabilistic mutation ( μ), the population does not fixate in any stationary state. Thus, instead of computing the probability of fixation in each absorbing state, we can make use of the stationary distribution of this Markov chain to analyze the asymptotic state of the population, and assess the pervasiveness in time of a each fraction of cooperators. To this end, we build the transition matrix p k,k −1 and 0 for the rest. This is a tridiagonal matrix and the stationary distribution is the normalized eigenvector ( p k ) k =0 .Z corresponding to eigenvalue 1 of the transpose of matrix . Third, the expected final number of cooperators is calculated from the stationary distribution, as n C = Z k =0 p k k . We can also use the hypergeometric distributions above to calculate the probability a G (k ) of having groups with M cooperators or more for every cooperation level k [14,34] . This is given by ing the average group achievement of the game in the stationary solution.

Agent-based computer simulations
While convenient for the case with homogeneous income losses ( c), the analytical mean-field approach described above is no longer valid in the case of heterogeneous exposures or losses. The inclusion of regions or countries' diversity in exposure introduces a higher level of complexity may, nonetheless, be conveniently described through Monte-Carlo (MC) agent-based simulations [28,29] , performed in computer clusters and resorting to parallel computing architectures. Evolution proceeds in discrete steps involving imitation and mutation, in line with the stochastic dynamics described above.
We fix the size of the groups N to 10 and the minimum thresh-

Results and discussion
In this section we analyze the results of three different experiments by considering the effects of the economic interdependence (EI) in the dilemma. First, we present the analytical study of the game when having homogeneous costs c in Section 3.1 . We later apply agent-based simulations to the dilemma with costs heterogeneity ( Section 3.2 ) and study the roots of the dilemma in Section 3.3 . Finally, Section 3.4 provides guidance on how to engineer interventions for increasing cooperation by taking into account heterogeneity.

Analytical study of the cooperation increase when considering EI
First, we analyze the effect of the income loss on the final stationary state of the game. Fig. 1 shows the average final percentage of cooperators as a function of the homogeneous income loss c for three risk levels. To analyze the effect of economic interdependence (EI), we represent the trajectories with two models: one adopting the payoff Eqs. (1) and (2) (model with EI) and other one using the same equations but removing the term c i k Z (model without EI).
As it can be observed, the income loss negatively influences on the final number of cooperators. In general, its effect is not strong for income loss values near zero, but up to a certain threshold, the expected number of cooperators dramatically decreases until  achieving almost null values. The lower the risk associated with the game, the lower the levels of income loss where null cooperation is achieved.
Moreover, Fig. 1 shows the positive effect of EI on the cooperation in the three risk levels shown. Therefore, the existence of economic linkages among countries in the pandemic context favors cooperation in the dilemma. The gap in the number of cooperators between the two models is larger for larger risk values. This result is apparently counter-intuitive, since including new economic costs to the players may in general deter the disposition to cooperate. However, this is not the case in this context, since the income loss due to EI reduces the overall economic benefit at stake for all players and therefore shortens the benefit margin of being defector.
A deeper insight of the analytical solution of the model with and without EI is presented in Fig. 2 , which shows the gradient of selection G (k ) for the same three risk levels in Fig. 1 , with and without EI. If the gradient of selection, G (k ) , is positive (negative), the number of cooperators is likely to increase (decrease) whenever the population has k cooperators. The roots of G (k ) offer the finite population analogues of fixed points in infinite populations [9,33] . Here, we can identify configurations with two internal roots typical of this class of dilemmas: one unstable root on the lefthand side ( x L ) and one stable root ( x R ) for higher fractions of co-operators, associated with two well-defined basins of attraction. It also suggests that a critical mass of cooperators ( x L ) needs to be surpassed such that the system naturally self-organizes into a coexistence of cooperators and defectors.
Interestingly, the gradients in Fig. 2 show that stable equilibria ( x R ) occur for higher levels of cooperators when EI is included. At the same time, the coordination point ( x L ) tends to move towards lower values of the fraction of cooperators k/Z when EI is in place. This suggests that EI reduces the requirements to sway to the cooperative basin of attraction where cooperators and defectors may co-exist. Naturally, the position of both ( x L ) and x R depends on the value of risk ( r). The higher the perception of the risk, the lower the unstable equilibrium x L and the higher the stable equilibrium x R , being more likely to overcome the dilemma. The position of these roots -but also the amplitude of G (k ) -determine the final stationary distribution and the expected prevalence of cooperators. Finally, the exploration (or mutation) probability ( μ) drive the system to the center of the cooperation axis by making G (0) > 0 , further favouring the shift towards the right-hand side of the simplex.

Simulation-based results for heterogeneous income losses based on real data
As commented in the introduction, the income loss due to the implementation of pandemic control measures is not homogeneously distributed among countries. For example, tourism was one of the most affected sectors by the COVID-19 pandemic and tourism GDP contribution of the European regions is one example of the heterogeneity of the players involved, as discussed in [17] . The number of regions in the EU NUTS2 classification [35] is 312 and we use real data from this classification to compute the exposure of the regions to a lockdown and tourism activity halt. The considered indicator for this contribution are the nights spend at tourist accommodation establishments per inhabitant. Although the averaged contribution of the regions is 0.04, we observe a clear heterogeneity in the distribution. Fig. 3 shows the 312 data points and two fitted log-normal and power-law distributions.
The log-normal distribution with fitted parameters μ ln = −4 . 39 and σ ln = 1 . 63 is the one with the best fitting. Therefore, we use this distribution to feed the costs or income loss c i for the players. Because of the heterogeneity in the costs, it is not possible to fully assess the overall dynamics through the mean-field analysis of the previous section. Thus, in the following, we resort to agent-based simulations to estimate the evolutionary dynamics of the model [17,29] . We have tested both log-normal and power law distributions of potential income losses without significant changes and therefore, the same conclusions apply for both distributions. Fig. 4 shows the averaged final cooperation frequency for the game with and without EI using the log-normal distribution of income loss. To obtain each point of both curves at every risk value, we average the simulations results from a set of sufficient discrete values for the initial frequency of cooperators of the system, n 0 C . We can see in this figure how the increase in cooperation is clear for all the risk levels after averaging all the possible initial conditions of cooperation. The inclusion of EI boosts cooperation even for very low risk values and the positive gap in cooperation remains practically equal for the whole range of risk values. Therefore, we see that the main behavior of the system is the same when injecting heterogeneity through the real distribution of income loss c i . The results are in line with the homogeneous setting of the model. The global increase in final cooperation when incorporating EI in the dilemma is robust, independently from the heterogeneity of the income loss.

Study on the internal fixed points of the game
The adoption of a different exposure for each agent makes it hard to assess the overall dynamics from an analytical perspective. Particularly, the inference of the effective gradient of selection and associated equilibrium from large-scale computer simulations is not always trivial (see, e.g., [36] ). In this section we aim at explaining why there is an increase in cooperation when simulating the model with heterogeneous costs and estimate the roots of the dynamics as in the analytical model. Therefore, we first propose a methodology to estimate, from the sensitivity analysis on risk r and initial cooperators n 0 C , the stable and unstable equilibrium points. Fig. 5 shows this sensitivity analysis on r and n 0 C where each cell is the averaged final frequency of cooperators n C from the MC realizations of the agent-based model using the log-normal distribution of income loss c i . The plot on the left represents the abso-lute n C values of the game without EI while the plot on the right shows the relative increase in n C of the game including EI with respect to the traditional dilemma.
From the heatmaps in Fig. 5 , we see a clear increase in cooperation for the majority of r values. We may also take profit from the knowledge obtained in Fig. 2 regarding the existence of two internal fixed points to, analogously, try to find these roots in the heterogeneous case. As before, the heatmaps portray two basins of attraction: one in which cooperation cannot be sustained (yellowish areas) and other where cooperators and defectors co-exist (blue areas). In order to estimate the stable and unstable points from the heatmaps, we define the next calculation method: • Unstable fixed points ( x L ): For every value of risk r, we compute the highest increase in the final cooperation values n C in ascending order by initial cooperators n 0 C . For this calculated maximum increase, we set the unstable point as the pair of initial ( n 0 C ) and final cooperators ( n C ). To avoid small variations or numerical simulation errors in the calculation process, we discard small increases for the final cooperation values (e.g., the area where risk is below 0.1 of the left plot of Fig. 5 ). • Stable fixed points ( x R ): For every value of risk r, we select the lower initial cooperation value ( n 0 C ) where the majority of the players of the population are cooperators. This is done by defining a threshold for final cooperators n C to consider the population as fully cooperator. We set this threshold to 0.85 for extracting the stable points of our experiments. Fig. 6 shows fixed points for the population dynamics in both analytical and simulation-based approaches, with and without EI. Empty circles represent unstable fixed points ( x L ), and full circles represent stable fixed points ( x R ), coming from the simulationbased outputs of the heterogeneous setting, following the abovementioned method. Lines are calculated from the analytical approach when using homogeneous costs (i.e., c i = 0 . 04 , ∀ i ∈ Z).
If we compare both analytical and simulation-based points, we see results are in line even when having different heterogeneity settings. The incorporation of EI when there are pandemic risks that can affect the payoffs of the players by others' flow restrictions clearly shifts the fixed points. Stable points are shifted to the right and unstable points are shifted to the left when including this  EI effects. The increase in cooperation is given by this extension of the internal points in the dilemma. We also show that the method to obtain the internal points by exploiting the simulation results is robust and the points are equivalent to the ones obtained by the analytical approach, obtained the same conclusions. Therefore, heterogeneity and the use of simulation techniques are not generating significant differences and results are solid.

Using heterogeneity to boost cooperation by fixing the initial cooperation conditions
Given the reality is heterogeneous and the fact we were able to incorporate this heterogeneity in the model by using simulationbased techniques, our aim is to exploit this information to glimpse ways of boosting cooperation. These insights can serve as a kick-off for employing policies by institutions. In this experiment we have fixed the initial conditions of the fraction of cooperators n 0 C based on the heterogeneous values of cost or income loss of the countries or players c i .
Specifically, we have defined three main scenarios. In the first one, the initial cooperators are selected at random from the members of the population and there is not any initial bias. For the Fig. 7. Comparison on different initial conditions n 0 C for achieving the best final frequency of cooperation n C in the population. Initial coordination of players with the lowest income loss can drive the entire population to higher levels of cooperation. On the opposite, if the initial cooperators face a significant income loss, cooperation tends to decrease compared with the random case. This result shows that heterogeneity leads to different liabilities, depending on the level of exposure. Moreover, it suggests that one may profit from heterogeneity to trigger prosociality at the global level. Parameters: Z = 2 × 10 3 , N = 10 , M = 7, β = 1 , and μ = 1 /Z. second one, we start fixing the initial cooperators from high to low values c i . And finally, the third scenario considers the initial cooperators from low to high values of c i . Thus, the second scenario has a positive bias of initial cooperation with respect to cooperation costs c i while the third scenario has a negative bias. Fig. 7 shows the n C comparison for a sensitivity analysis on risk r using the three scenarios for the dilemma with EI using heterogeneous income loss (i.e., running the agent-based simulations). In order to get the lines of the plot we averaged a sufficient set of values for initial frequency of cooperators n 0 C by following each of the three scenarios. The set of values is from n 0 C = 0 to n 0 C = 0 . 5 as policies should be focus on a reduced number of players in the population and it is not suitable for the comparison to restrict initial cooperators to a high number.
The results of the plot are clear. When the system starts with those players having the lowest income loss as cooperators, the final cooperation of the population increases significantly. On the opposite, when we fix the initial cooperators to those having the highest income loss, the final cooperation even decreases with respect to the random initial cooperators setting. Final cooperation increase of negative bias with respect to positive bias is between 10% and 50% depending on the risk values. Just when risk values are below 0.1, all the scenarios have equal results as the dilemma does not facilitate any final cooperation.

Conclusions
COVID-19 pandemic and other global risks have changed how regions or countries interact concerning global failures such as economic knockdowns. In these cases, there are economic interdependences (EIs) or spillover effects among players when facing public goods games. For instance, economic or mobility restrictions set by some players can affect the expected payoff by defecting ones or free-riders . To cope with this new reality, we proposed a collective risk dilemma (CRD) that includes EI effects among players. Real data from the tourism contribution of the EU regions is employed to enrich the model by setting a genuine heterogeneous distribution of cooperation costs for the players.
We show that EI robustly increases cooperation for both homogeneous and heterogeneous cases. EI is able to modify the (finite population analogues of) internal fixed points when compared with the dynamics in the absence of EI. We depart from a classic CRD characterized by defector dominance and coexistence among cooperators and defectors, each outcome associated with two welldefined basins of attraction. In the absence of any additional community enforcement mechanism, EI drops the minimum number of cooperators required to reach the cooperative basin of attraction, and increase the prevalence of cooperators in coexistence point. We have computed these fixed points analytically for scenarios with homogeneous costs and through agent-based simulations in the case of heterogeneous costs. To this end, in the latter, we proposed a new method to infer these finite population analogues of stable and unstable fixed points in infinite populations from the simulation's outputs.
Finally, we have discussed how biased initial conditions based on the level of exposure may alter the final expected outcome. Results showed that the entire population benefits from having cooperators within the sub-group of players with lower income loss. This result suggests that one may profit from heterogeneity in designing effective interventions or governance policies to trigger prosociality at the global level, a result of particular importance if we consider that individual strategies may also depend on the perceived risk by each party [37] . Interventions should focus on those players showing a lower exposure to the economic risk.
The analysis of the effect of EI on cooperation could be extended in several ways. For example, the model assumes a wellmixed interdependence among players that does not necessarily fit the reality. In fact, the exposure degree for every player estimated with the real data of the EU regions represents an aggregated value emerging from a specific network of interdependences among regions, although this information is not available. Considering a network structure of interrelationships among players would enrich the possible scenarios and solutions for cooperation. For example, further model extensions can analyze the effects of different network topologies (e.g., scale-free networks and data-generated networks as in [38] ) as well as local and global ways of peer-influence among players.
Additionally, future work can assess if rewarding and sanctioning activities can be applied [15,34,[39][40][41][42]] to a specific target sub-population and the features of this subset of individuals to target, or how positive, and negative incentives can be optimally distributed among groups and actors. Moreover, reactions to the COVID-19 pandemic have shown a wide range of (often polarized) responses. Recent results have shown how uncertainty may influence how each individual perceives the dilemma [10,[43][44][45] , potentially leading to polarized reactions [16] , a development yet to be studied in the context of cooperation dynamics related to managing economic losses under pandemic conditions. Finally, how leaders act and influence others by their example and reputation can affect the whole population outcome [46] , and this phenomenon can be studied for this dilemma. All these open questions remain critical in the current quest of understanding and promoting human cooperation, given the difficulty in assessing the advantages and disadvantages of each possible type of intervention policies.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.