Replicator dynamics with diffusion on multiplex networks

In this study we present an extension of the replicator equation with diffusion to multiplex graphs. We derive an exact formula for the diffusion term, which shows that, while diffusion is linear for numbers of agents, it is necessary to account for non-linear terms when working with fractions of individuals. We also derive the transition probabilities that give rise to such macroscopic behavior, completing the bottom-up description. Finally, it shown that the usual assumption of constant population sizes induces a hidden selective pressure in the system.


I. INTRODUCTION
During the last decade agent based modeling has increased its importance as a powerful tool to model situations in which the complexity of the interactions of many-agent systems makes it impossible, or at least very difficult, to make analytic predictions of the dynamical behavior of the system [1,2]. Furthermore, the agent based models complement the analytic approach allowing for an exploration of the coarse-grained dynamics and the connection between micro-scale and macro-scale behavior [3][4][5][6]. Such agent based modeling is of special importance in evolutionary game theoretical studies [7,8] which aim to capture the intricacies of biological, social and economical systems, where the non-linearity and feed-backs of the systems cannot be easily foreseen [9][10][11][12]. In the middle of such framework [13], and connected with the micro-evolutionary dynamics [3][4][5][6], stands the replicator equation [14,15] The replicator equation, which was introduced shortly after the foundation of the evolutionary game theoretical framework [7,16,17], has been extensively used during the last decades to model the evolution of the fractions x α = n α /N of traits of type α, α = 1, ..., L in large well-mixed populations with frequency dependent selection, i.e. when the fitness f α (reproductive potential or capacity) of the agents traits and the mean population fitnessf = α f α x α depend on the population composition. Such traits may represent different phenotypes and genotypes in biological settings, or different behavioral strategies in a cultural evolutionary framework.
In addition to classical biological applications related to the evolution of gene frequencies and phenotypic traits, the replicator equation and agent based simulations using microscopic updating rules that lead to it have been * Electronic address: rrequejo@ffn.ub.edu used to study how evolution, with a stress on the evolution of cooperation, is affected by physical entities. These studies include the interplay between different temporal scales of interaction and selection [18,19], the effect of network structures, both spatial lattices [19][20][21] and scale-free networks [22,23], and the effect of linking fitness and resource availability [9,10,12,24], showing a self-organized feedback similar to homeostatic regulation [25]. Furthermore, the richness of dynamical portraits [26,27] has also been shown in minimalist scenarios, as exemplified by the introduction of loner [28,29] and joker [10,30] strategies, which allow for several phase transitions and complex dynamical behavior [10,31,32], or the use of the discrete time version of the replicator equation for two strategies, which has been proven to show periodic and chaotic behavior [11].
The replicator equation is on the core of the framework of evolution [13], linked to the quasispecies equation [13], the game dynamical equation [3,4], adaptive dynamics [13], Lotka-Volterra dynamics [33,34] and the Price equation [35][36][37][38]. The extension of the replicator dynamics to regular networks has been developed in a situation where each node is an agent [39], finding that for weak selection it corresponds to a transformation of the payoff matrix. Furthermore, an extension of the replicator dynamics has been developed to represent a two dimensional world where the agents diffuse [32], which shows the appearance of Turing patterns and other complex behavior. However, many real world situations cannot be modeled as a simple network or a two dimensional space, and need the introduction of several kinds of links to represent different properties of the system [40], as in the air transportation networks, where each airline represents a different network [41].
Let us focus on a cultural evolutionary framework along the rest of the paper. In such a framework the individual traits being selected are behavioral traits, called strategies. The main aim of this paper is to develop the mathematical tools that allow to model diffusion of strategies in multiplex networks in a compatible way with the selection dynamics described by the replicator equation (with or without mutations), irrespective of the microscopic dynamics that give rise to the replicator equa-tion (see appendix A). Let us provide a simple example to illustrate the situation: imagine several cities which are connected by bus, train and plane, each kind of connection with a different network structure. Individuals in each city are interacting between them, both directly or indirectly: they have different jobs and different incomes, which determines, at least partially, their choices on how to travel, and they also have information about the transportation ways. Based on such interactions and information, they may decide to travel in one or another transport; such choice determines the individual strategy. As the different transportation ways determine different networks structures, as well as the travel speeds are different, a multi-layer network is necessary to represent the full transportation system; hence, as a first approach to such simplified situation, we need to extend the replicator dynamics (describing strategy changes) to diffusive individuals in this kind of multi-layer.
We may think that introducing diffusion on multiplex models of evolution is a straightforward task, as the diffusive process has extensively been studied on such networks [42,43]. However, it poses a challenge: as the replicator dynamics describe the evolution of fractions of individuals, the diffusive models need to be rewritten in a compatible way, accounting for the constraint on the addition of the fractions to one, as well as for the conservation of the number of agents. We develop such extension in this paper, showing that working with fractions introduces dependencies which are not present for the diffusion of the numbers of agents. These new dependencies can only be taken into account introducing a non-linear term, which has to be added to the linear one, in order to represent the general dynamics of diffusion of fractions of individuals in the multiplex. Furthermore, we discuss some situations in which the linear scenario can be recovered, as when population sizes or population size ratios between sites are constant, and show that in such situations hidden selective pressures act on the system, even when they do not appear explicitly on the equations.
The paper is outlined as follows: We start in section II by defining multiplex networks and showing the problem to overcome when working with fractions; after that, in section III we derive the diffusion term compatible with the replicator dynamics and infer the transition probabilities that give rise to it; then, in section IV we discuss some situations in which the linear dynamics are recovered (and the analytic calculations simplified), and show that some of this situations carry attached the appearance of hidden selective pressures; finally, in section V we discuss the results. In a general multi-layer network [40] each node in one layer can be connected to any node in any other layer. Example of multiplex network formed by three layers (blue, gray and red). The set of nodes is the same in all layers and the inter-layer links connect them in a one to one basis. For clarity, the inter-layer connection between the first and third layer has been omitted.
Hence, the connectivity may be given by a tensor with four indices, M αβ ij , whose entries are 1 if nodes i α -node i in layer α-and j β are connected and zero otherwise. In multiplex networks [42][43][44][45][46], however, the set of nodes i = 1, . . . , N is the same in all layers α = 1, . . . , L, and the connectivity of each layer is given by a matrix A α = {a α ij } (see Fig.1). Furthermore, only connections among one node i α and its counterpart in another layer i β are allowed. This inter-layer connectivity is the same for all nodes and is given by the inter-layer connectivity matrix Λ αβ [43]. Hence, multiplex networks have a connectivity defined by with i = j and α = β. Each site i of a multiplex may be regarded as one structural entity, and the different layers α represent different interconnection structures between these entities. Structural entities (sites) may represent specific locations (cities), and the different connectivity of each transportation system or mobility pattern would define the different layers through which the agents move, while the interlayer connections refer to the possibility to reach another transport (layer) from a specific location (for simplicity we will assume that Λ αβ = 1). Individuals changing from one layer to another can be represented as changing their strategy due to selection (due to prices, availability, or other competitive reasons) or mutation (random trial) in an evolutionary framework.
In order to describe the state of the entire multiplex system, we may define a set of vectors { n i }, one per site i, each one with L (number of layers) components. From an evolutionary game theoretical point of view the components of such vectors represent the number of agents n α i in a given position i α of the multiplex, and hence n i is the population composition at i. Then, the evolution of the state of node i α may be assumed as given by a functional F α i [{ n i (t)}, a α ij , Λ αβ ], which depends, respectively, on the state of the system, and on the intra-and inter-layer connectivities.
As the state of the system is instantaneously defined by a set of quantities n i for each site i, we may write the dynamics as a set of coupled differential equations, one per component (layer), where the terms on which the functional depends are, respectively, the state of i α , the state of connected nodes in the same layer (intra-layer neighborhood of i α ), and the state of equivalent nodes connected through interlayer connections (inter-layer neighborhood of i α ).
In order to approach a replicator-equation like functional (see Eq.(1)), we have to normalize n α i with respect to N i = α n α i , the number of agents (population size) on site i, which is only site dependent. The normalized quantity is the fraction x α i = n α i /N i , which must fulfill the restriction α x α i = 1, and hence its derivative satisfies αẋ α i = 0. Note that, whenever the agents diffuse from j β to i β (or in the opposite direction), they are modifying the value of the fraction in i α through the modification of the population size N i . Hence, in order to account for the dynamics in a position i α it is no longer enough to take into account the direct neighborhood of i α , given by its connectivity, but it is necessary to take into account the entire neighborhood of i, introducing an extra dependence. The functional in terms of fractions is thus where the extra dependencies are given by the first term between brackets, which now accounts for the entire state of i, and the last term, which accounts for the neighborhood of i β . As we show in the following, such extra dependencies require a modification of the diffusion term, which will no longer be linear, but include a non-linear term.

III. REPLICATOR DYNAMICS WITH DIFFUSION IN THE MULTIPLEX
In this section we derive a diffusion term in the multiplex compatible with the replicator dynamics, i.e. describing the evolution of the fractions of agents x α i at site i with diffusive pattern given by layer α, and keeping the conservation of the total number of agents in the multiplex, i N i = N where N is constant, as well as the constraint α x α i = 1. We will assume that diffusion and evolution are uncoupled, and hence the diffusion term can simply be added to the dynamics.
As the diffusion of agents has been already studied, let us start trying to derive the equations for diffusion of the fractions from those for the number of agents, and discuss its use in game theoretical studies. If the agents are diffusing across a network structure with adjacency matrix of layer α given by the elements a α ij = 1 for connected nodes and 0 otherwise, the transition probabilities determining the microscopic dynamics of diffusion of (numbers of) individuals i given j in layer α are where D α is the diffusion coefficient of layer α and is the transition rate of increase (decrease, by changing all signs) of the number of agents in i α in one unit triggered by the movement of one agent in j α . Note that, as the process of diffusion impilies a redistribution of agents, the aggregated number of agents N i +N j is preserved in each diffusive event among i and j, and hence the constraint (T +α i|j ) diff = (T −α j|i ) diff holds, which also ensures the global conservation of agents in the entire multiplex system by linking the processes n α i → n α i + 1 and n α j → n α j − 1. In the following, for simplicity, we will write the transition rates in a simplified form, explicitly stating the process to which it refers in the focal variables and the terms involved, but not the linked process; in this way, Eq.(6) would be T [n α i → n α i + 1 | n α j ] The microscopic dynamics can be connected with the macroscopic behavior of the system by expanding a Fokker-Planck equation and truncating high order terms (this happens naturally for N 1 when working with fractions, see appendix A for the one dimensional derivation), which results in the Langevin equation, where ξ is uncorrelated Gaussian noise, η is either the number of individuals or the fraction of individuals, and the drift and diffusion terms (do not confuse the latter with the diffusive dynamics studied in this section) are respectively where the transition probabilities have to be written in terms of numbers or fractions depending on the choice in Eq.(7). The first term in Eq.(8) (drift) accounts for the deterministic behavior of the system and the latter (diffusion term) for stochastic effects. Note that the stochastic effects disappear in the thermodynamic limit N i → ∞, or whenever N i j (T +α i|j +T −α i|j ), i.e. when the transition rates are small compared to the population size.
Whenever we introduce the transition probabilities Eq.(5) into the drift term in Eq.(8) expressed for numbers of individuals, the deterministic diffusive dynamics for the number of agents are well known equation for the diffusion of particles on a network, where the tensor is the graph Laplacian of the corresponding layer α (with δ αβ the Kronecker's delta) and k α i is the degree of site i in layer α.
In order to derive the equation for the diffusion of fractions instead of numbers, it is important to note that, if we differentiate the fraction x α i = n α i /N , we geṫ and hence a non-linear term that depends on x α i appears. Since we are trying to find the exact description of the diffusive process which is compatible with the replicator equation for mobile agents in a multiplex, we can follow this approach: first, construct a compatible macroscopic equation and then, infer the transition probabilities that give rise to it. In order to do this we may introduce Eq.(9) into Eq.(10) (note that the latter is a replicator equation of the form of Eq.(A1)), and use β n β i = N i , obtaining a term of the form, where ρ ij = N j /N i is a population size dependence between neighboring sites. As it can be observed, in addition to a linear term, the first one, a second term appears. This term ensures that the normalization of the fractions is the proper one, i.e. fractions always add up to one, as shown in Fig.2(a),(b) for a system of two layers and two nodes in each layer. Note that, as previously foreseen in section II, this term implies a dependence of the dynamics at i α on positions j β to which it is not directly connected, but represents the neighborhood of i β nodes.
The diffusion term (11) may be rewritten in a more compact form Note that, if one extracts the node degree from the second parenthesis -the network Laplacian term-, then it becomes k β i (δ ij − a β ij /k β i ), which is similar to the first parenthesis term, ρ ij (δ αβ − x α i ), but with constant values instead of variables. In this way the latter term may be interpreted as a population composition dependent Laplacian, which accounts for the instantaneous heterogeneity in population sizes and fractions of individuals in the network. Now, it is possible to infer the transition probabilities of the microscopic process from the emergent macroscopic dynamics (Eq. (12)), which may be used in Markovian analysis [12]. The macroscopic deterministic dynamics emerge from the drift term in Eq. (8), which by extrapolation to two dimensions results iṅ where is the transition rate of increase (decrease, by changing all signs) of the fraction of agents in i α given the fraction of agents in j β (see appendix C). Note that x β j may also vary in the same process, but such variation is specified in its related transition rate, and not explicitly written here for simplicity.
Then, we can compare Eq. (12) and Eq.(13) (see appendix C for the derivation), and infer the bi-dimensional transition probabilities describing the variation of the fractions of individuals x α i which give rise to Eq. (12), resulting in where we have omitted the subscript 'diff' for simplicity.
Note that the transition probabilities depend on both the layers (α) and the sites (i), and that the symmetry which is present for mutation transition probabilities, which makes them keep local population sizes constant (see appendix B), is broken in general settings, T −α|β = T +β|α . Such asymmetry is however expected, as diffusion implicitly needs variable local population sizes (although the total number of particles in the multiplex is conserved), and assuming constant population sizes locally may have unexpected effects, as shown in the next section.
The transition probabilities in Eq. (15), together with those that give rise to the replicator equation with mutations (appendices A,B) complete the microscopic description of the evolutionary process for agents in a multiplex. Let us finally write down the deterministic replicator dynamics (see appendix A) with additive mutations (see appendix B) and diffusion in a multiplex: where f α i is the fitness of individuals in position i α ,f i is the mean fitness of individuals in i across layers and q αβ i is Eq.(11) (Full eq.). Fraction of individuals, x γ i . FIG. 2: Intra-layer diffusive dynamics of a system consisting of two layers with diffusion coefficients D α = 10 −1 and D β = 10 −2 , and two sites, i = 1 and j = 2 . The initial configuration is with N individuals in nodes i α and j β which are allowed to diffuse within its own layer (note that N may take any value in the situation depicted, provided that it is the same in both layers, as the fractional character of ρ makes it vanish); no inter-layer process is acting on the system. In order to fulfil the constraint  Remarkably, the non-linear effects due to diffusion are all contained in the x α i in the first parenthesis of the latter term in Eq. (16). There are however some situations in which the non-linearity disappears, and in which the analytic calculations can be simplified. Let us discuss them, as well as their implications, with special emphasis on the appearance of hidden selective pressures.

IV. RECOVERING LINEAR DIFFUSION AND SIMPLIFYING THE ANALYTICS.
Some situations allow us to recover linear diffusion, as well as to simplify the analytic calculations. Here, we analyze three of such scenarios. The first one refers to a situation in which population sizes are all forced to be constant across sites, and hence the second term in Eq.(10) disappears. Remarkably, in this situation a hidden selective pressure appears favoring the increase of fast diffusing strategies. The second scenario refers to slow population change, situation which approximates the constant population sizes scenario. The last subsection explores the situation in which the population ratio can be expressed as some analytic function. In that case, time scales separation allows us to write some formulas which simplify the calculations.
A. Constant population ratios (or sizes) induce hidden selective pressure.
Let us first study the case in which some mechanism makes ρ ij = N j /N i → c, where c is a constant; for simplicity, we will assume c = 1, which includes the usual assumption in evolutionary game theoretical studies of constant population sizes [3][4][5][48][49][50]). This case is shown in Fig.2(b), where the restriction on the addition of fractions of individuals to one during the dynamical evolution of the system is fulfilled. However, the final state is not a symmetric one in which there are N/2 individuals in each node, as it happens without the restriction on the quotient of population sizes. Why does this happen?
The case in which ρ ij is forced to be equal to a certain value in Eq.(11) may happen only if the conservation of particles due to the diffusive process does no longer hold. In such case, which may be due to a fast (compared to diffusion) and neutral evolutionary process acting on the system, the difference in diffusive velocities transforms into different "diffusive pressures", which induce different selective pressures while the system decays to the equilibrium.
The different selective pressures are induced by the fact that fast diffusing individuals are increasing their frequency in a site compared to slower diffusing ones, and then the increased fraction is fixated by the neutral selective process, which only renormalizes the population size without altering the proportions. This is similar to a Wright-Fisher process, in which the population reproduces during the reproduction period according to their fitness and then, keeping the proportions of individuals, the population size is renormalized to its initial value. In our case, however, the variation in the fractions is due to the diffusion of the agents and the continuous renormalization of the population size is due to the neutral evolutionary process, thus favoring the fast diffusive strategy, as shown in Fig.2(b).
B. Slow population size change and quasi-neutral selection.
The equation describing the evolution of the fraction of particles present at each point may in principle be simplified whenever, starting from a situation near the equilibrium, the variation of the number of agents is so slow that it can be neglected,Ṅ i /N i → 0. This is equivalent to assuming thatẋ α i ≈ṅ α i /N i , which results in 17) However, Eq.(17) does not generally keep the proper normalization for the fractions β x β i = 1, as it can be easily proven with a simple example. Let us assume a multiplex network formed by identical networks with identical diffusion coefficients (for two-dimensional spatial networks, this equals the Fisher-Kolmogorov reactiondiffusion scenario for gene wave-front propagation [52][53][54]). In this case, which is analogous to a mono-layer network, there is no dependence of the adjacency matrix terms and diffusion coefficients on the layer index α, which now serves only to identify the different strategies present in each node. Hence, by summing Eq.(17) over layers (or strategies) and noting that β x β i = 1, we obtain the condition βẋ which is only equal to zero, i.e. satisfies the dynamical constraint, whenever These two restrictions are equivalent to requiring that j L ij N j = 0, which could be implemented or engineered in the system, but it is not a priori expected to happen as a self-organizing feature. Furthermore, simulations using a system with two nodes and two layers confirm that the normalization is not fulfilled using Eq.(17), as shown in Figs.(2)(c),(d). Hence, the construction of a microscopic model that describes the macroscopic diffusive dynamics of fractions of individuals cannot be done by simply rewriting the transition probabilities in Eq. (5) in terms of fractions.
The case of slow population change is equivalent to assuming that the second term in Eq.(11) vanishes and is hence only slightly influencing the dynamics described by the first term. In such case, the diffusion term approaches Eq. (17). Note however that, if we assume that the term above is strictly zero, we recover the case of strictly constant population sizes, and hence hidden selective pressures may appear, as explained in the previous subsection. From an evolutionary perspective, small perturbations introduced in the neutral selection limit satisfy the conditions leading to the slow population size change approximation whenever diffusion is slow. If we define the quasi-neutral selection limit as represented by f α i → 0, then the population size at each site varies slowly, given thatṄ i = N ifi → 0 and D γ → 0. The quasi-neutral selection limit is important for two reasons: first, it allows for analytic calculations to be carried out in a similar way to the weak selection limit [55,56] and, second, this limit approaches the neutral theory of molecular evolution proposed by Kimura [57]. However, care should be taken when approaching this limit, as explained above.
Let us finish exploring the situations in which ρ ij can be written as a function of the fractions of individuals (or their fitness, which are determined by such fractions once the payoff matrix is known) and time due to time-scales separation.
Whenever we are able to write the population size dependence of the diffusion term as the entire replicator dynamics in the multiplex, including replication and deaths, mutations and diffusion, can be written as a function of the state of the system, given by the fractions of individuals in each node, and of the structural terms (diffusion coefficients and Laplacian of the multiplex).
In the cases in which it is not possible to write an explicit dependence for ρ ij as above, there are at least two situations in which time scales separation allows for approximations that take such form. The first one is whenever diffusion is slow and most of the population size change is due to replication and death. In this case it is easy to prove that the differential equatioṅ governs the evolution of population sizes. Note that the solution of this equation is an exponential integral, which implies that the system has memory. More precisely, the entire history of the difference of mean population fitness differences, which depends on the population compositions, is influencing the present state. Hence, Eq. (22) is a memory kernel of the past states of the system. As the memory kernel has an exponential form, although the entire history is contained in it, the influence of past states decays very fast with time compared to present states. This can be easily proved noting that, given two time lapses beginning at t = 0 and ending at t 1 and t 2 > t 1 , the memory ker- (fj −fi)dt , which is a memory kernel of the time lapse between t 1 and t 2 and does not take into account the time lapse between 0 and t 1 . Hence, we may always make t 2 = t 1 + dt to express the quotient as an instantaneous integral of the fitness differences, which allows for computation of the dynamics without keeping track of the entire history of the system.
The opposite limit to slow diffusion is the fast diffusion limit: In this case we may assume that, after a short transient, the population size variations fulfillṄ i /Ṅ j = κ ij (t) = 0, ∞; then the approximation can be used. Furthermore, in the fast diffusion limit we may assume that the population size variations due to replication and death may be negligible compared to diffusion, as in the previous subsection; in such case the conditions in Eq.(18) -equal populations across all sites or a local mean field describing the population sizes-may be prone to happen, although again, this may introduce hidden selective pressures.

V. DISCUSSION AND CONCLUSIONS
We have presented a complete description of the diffusive process in terms of fractions of individuals which is consistent with the replicator dynamics. Such term can be added to the replicator equation in order to make a full description of a system of diffusing and evolving agents in a multiplex. We have also found the transition probabilities that describe the microscopic dynamics from which the macroscopic behavior emerges. As we have shown, due to the multiplex structure it is necessary to include in the diffusion term a dependence on the quotient of population sizes of neighboring sites, as well as to add an extra diffusive term which is non-linear; this happens in order to keep the restriction imposed by the addition of fractions to one, as well as to account for the conservation of particles due to its diffusion.
The extra non-linear diffusive term in Eq.(11) takes into account the state of the vicinity of the focal node and its equivalent nodes in all layers (extended neighborhood), and not only its direct neighbors (directly connected nodes). This is relevant for the calculation of the evolution of the fractions of individuals at each site, and contrasts with the usual diffusion term, which only takes into account the state of directly connected nodes within each layer.
We have finally explored the recovery of the linear scenario, finding that when population sizes are constant across sites (due to environmental saturation, for instance), the non-linearity disappears. However, in this case the diffusion process induces an extra evolutionary pressure acting while the system reaches the equilibrium. This has been argued to happen because, even if the evolutionary process does not alter the fractions of individuals at each site, and only acts so as to maintain the population size constant, the faster diffusing strategy increases its proportion in the neighboring sites faster than the slower strategy in the focal site, and hence the renormalization process favors its increase. This suggests that, depending on the network architecture, the induced evolutionary selective pressures may work so as to create extra gradients of selection acting while the system is out of the equilibrium, which may induce, depending on the multiplex architecture, an unexpectedly complex phenomenology.
The replicator equation in a well-mixed population (no network structure, Eq.(1)) can be derived -for large population sizes accepting smooth derivatives-by differentiating the fraction of individuals x α , and assuming that the fitness of the individuals corresponds to the instantaneous per-capita growth rate of each strategy due to replication and deaths [48], where x = {x α } is the state vector of the population, and the mean fitness of the aggregated population fulfills f = α x α f α =Ṅ /N . If the population size is constant, N = 0, then the per-capita growth rate is proportional to the difference between trait and mean population fitness, being hence the replicator equation applicable for constant and variable populations with slightly different fitness definitions. As it will be proven in the following, the replicator dynamics also emerge as the macroscopic description of some microscopic dynamics. Let us start the derivation of the replicator equation by assuming that there is a microscopic process, which can be described by some transition probabilities between states, and that such states are well defined. As we will assume that individuals of different types diffuse through different network architectures, each of such networks being part of a multiplex structure, let us introduce now the related notation: As before, each agent type will be labeled by the superscript α of the layer to which it belongs (related with its strategy), and the subscript i will refer to a site in such layer.
The probability for node i α to be at time t in a state with n α i individuals of α type will be denoted as P (n α i , t), and the probability of increasing or decreasing such number of individuals by one individual will be T + (n α i ) and T − (n α i ) respectively, where (similarly for T − with a sign change). With this notation, it is possible to write the master equation which describes the evolution of n α i . Now, let us assume that the total number of agents in site i is N i = α n α i 1. We do not require it to be infinite, but just large enough, so that we can make a continuous approach without neglecting finite size fluctuations. In this case it is possible to define the re-scaled variables For simplicity, let us assume that there are only two strategies present in the population, and the constraint α x α i = 1; we can then expand the master equation in a one dimensional Taylor expansion for N i 1 (the derivation is similar for the three- [49] and n-strategies cases [50] by using a multivariate Taylor expansion), giving rise to As the previous equation has the form of a Fokker-Planck equation, it is possible to transform it into the Langevin equationẋ where ξ is uncorrelated Gaussian noise and the drift and diffusion terms are respectively accounting for the deterministic behavior and the stochastic effects.
Note that, whenever the transition probabilities can be written as the drift term e(x α i ) looks like a replicator equation. This is indeed the only term acting in the thermodynamic limit N → ∞, and whenever the size of large populations does not increase or decrease too fast, N i x α i (R(x α i ) + D(x α i )). In this cases the Langevin equation simplifies toẋ and the terms R(x α i ) and D(x α i ) act as the replication and death components of the fitness difference f α i −f i (compare Eqs. (A10) and (1)). Hence, any process in which the transition probabilities can be factorized as shown in Eq.(A9), can be written as a replicator like system with fitness obtained as the solution of the equation This property has a particularly useful value: by using the latter equation and decomposing fitness in payoff components, it is possible to combine several processes, as virus spread (linear models) and cultural reproduction (frequency dependence, usually non-linear), into a unified evolutionary framework, as shown in [47].

Appendix B: Microscopic local updating rules and mutations
Whenever there is an arbitrary number of strategies in the population, the transition rates in Eq.(A9) can be decomposed in additive terms in a way in which each term relates to the contribution of each of the strategies. In this way Eq.(A7) can be generalized tȯ where T +α|β i is the transition rate of increase of the fractions of individuals in i, α corresponding to the increase in one unit due to the action of agents in β, and decreases them by one unit for the minus sign (note that this may imply a variation in the local population size N i ). The exact choice of the transition rates depends on the microscopic dynamics. These transition rates may be written as [49], whenever the transitions depend on a big number of random interactions between agents in the same site (local well-mixing assumption) in which one individual replaces another (defined by the second condition) and there are mutations between γ and α individuals at a rate q γα i . Factor g ±γ|β contains the information about how the microscopic updating rule acts: it states the exact mechanism by which one strategy increases or decreases due to the action of another.
Two probabilistic microscopic dynamics are usually investigated which give rise to the replicator dynamics (assuming q γα i = δ γα ). The first one is the modified Moran process, where one randomly chosen individual is assumed to die, and another individual, chosen according to a probability proportional to its fitness, reproduces. This process is defined by g The second process is proportional imitation, where one individual compares its strategy to another one, both chosen at random in the mean field limit, and changes its strategy with a probability increasing linearly with the difference of the payoffs between them. This process is defined by g where ∆f s,max is the maximum fitness difference and keeps the proper normalization. In both cases the symmetry condition g −α|β i = g +β|α i ensures that the evolutionary process maintains a constant population, and the dynamics result in the replicator equation, up to a multiplicative factor which relates to the temporal scale of the dynamics.
Whenever mutations happen as a strategy change at any point during the lifetime of the individuals, they can be introduced as additive terms [4] of the form to the transition rates in Eq.(B3), where the coupled mutations may be eliminated by setting q αβ i = δ αβ (with δ αβ the Kronecker's delta). The introduction of the additive term in the transition rates gives rise to the extra additive term in the replicator dynamics This term describes the effect of random mutations or equivalently of random exploration of strategies [51] in the evolutionary process. For the case of equal symmetric mutations, i.e. q αβ i = q βα i = µ, between all strategies this term simplifies to where L is the number of strategies [30]. Let us finally recall that, when mutations are coupled to the reproductive dynamics as in Eq.(B3), then Eq.(B1) gives rise to the replicator-mutator equation, The diffusion term in Eq.(12) could also be added to this equation to represent situations where only newborns mutate.
fractions of individuals (Eq. (12)) is defined by the transition probabilities determining the increase or decrease of the number of individuals, as shown in Eq. (13). In order to infer such transition probabilities, which complete the bottom up description, we can start expanding Eq.(12) as Then, we can split the double sum on the four contributions corresponding to the terms arising from i α , j α , i β and j β (assuming that α = β and i = j), obtaining By performing the same kind of split in Eq. (13) we finḋ Then, since the terms relate to the influence of different nodes k γ on the focal node i α , we can compare term by term both equations above, finding that each term in the first equation corresponds to the substraction of a pair of terms in the second.
In order to finally find the right transition terms, we need to take into account the physical constraints in the system. Let us analyze them one by one.
First, we can focus on the first term on the above equations. This term relates to the influence of i α on its own dynamics. Since the agents in i α are diffusing away from that position, the influence is necessarily negative, and hence T +α|α i|i = 0, and i denotes a dependence, not implying that such term is constant (indeed, it varies in the process). Note also that, whenever x α i = 0, the transition rate is zero, as expected due to the absence of agents diffusing away. Furthermore, if x α i = 1, the transition rate is again zero, as expected due to the fact that agents diffusing away decrease the number of agents, but leave the fraction unchanged.
Now, let us focus on the influence of j α on the dynamics of i α (j = i). Since agents diffusing away from j α are incresing the number of agents in i α , the influence is necessarily possitive, and the transition rates are T −α|α i|j = 0 and The limits, as before, can be easily proven to behave in the correct way.
Then, for the influence of i β on i α (β = α), we have to note that agents diffusing away from i β decrease the denominator of the fraction x α i = n α i / γ n γ i , since they are decreasing n β i , and hence the fraction x α i increases. Therefor, the transition rates are T −α|β i|i = 0 and Finally, individuals diffusing away from j β (j = i and β = α) are incresing the number of individuals in i β , and hence increasing the denominator in the fraction of x α i (opposite as in the previous case). The influence is then so as to decrease the fraction of individuals at i α and hence T +α|β i|j = 0 and From the transition probabilities above we can derive the diffusion term corresponding to the stochastic effects in the Langevin equation (extrapolating s to two dimensions in Eq.(8)), which is This term can be added to Eq.(12) -multiplied by white Gaussian noise-in order to account for the stochastic effects introduced by the diffusive process in situations in which the thermodynamic limit cannot be assumed, or population sizes vary fast.