Replicator equations induced by microscopic processes in nonoverlapping population playing bimatrix games
RReplicator equations induced by microscopic processes in nonoverlapping populationplaying bimatrix games
Archan Mukhopadhyay ∗ and Sagar Chakraborty † Department of Physics, Indian Institute of Technology Kanpur, Uttar Pradesh 208016, India
This paper is concerned with exploring the microscopic basis for the discrete versions of the stan-dard replicator equation and the adjusted replicator equation. To this end, we introduce frequency-dependent selection—as a result of competition fashioned by game-theoretic consideration—intothe Wright–Fisher process, a stochastic birth-death process. The process is further considered to beactive in a generation-wise nonoverlapping finite population where individuals play a two-strategybimatrix population game. Subsequently, connections among the corresponding master equation,the Fokker–Planck equation, and the Langevin equation are exploited to arrive at the deterministicdiscrete replicator maps in the limit of infinite population size.
I. MOTIVATION
The substantial literature on the replicator equationsacross disciplines epitomizes the pithy aphorism that “Allmodels are wrong, but some are useful”: In spite of be-ing devoid of many realistic details of the respective sys-tems, the replicator equations—the stylized equations ofevolutionary dynamics—have brought forth an insightfulplayground for understanding the theoretical aspects ofthe evolutionary game theory as studied in physics [1–4],chemistry [5, 6], mathematics [7], social science [8–12],economics [13–17], biology and ecology [18–22], and re-inforcement learning [23–30]. Strictly speaking, the de-terministic replicator equations, whether continuous ordiscrete in time, are models for the case of infinite-sizedpopulation of agents or individuals; but for a reason-ably large enough population, the replicator dynamicsare quantitatively decent approximations of the corre-sponding stochastic models [31–38].Probably the most fundamental aspect of the replica-tor equations is the direct connection of their asymptoticdynamical states with the Nash equilibria of the normalform games’ payoff matrices that are used in the equa-tions. This fact is known as the folk theorem [39, 40] ofthe evolutionary game theory. The theorem essentiallyasserts that a stable fixed point or a ω -limit point of thereplicator dynamics is a Nash equilibrium; and a strictNash equilibrium corresponds to a asymptotically stablefixed point. The connection of the evolutionarily stablestrategy/state with the fixed points of the replicator dy-namics is less general; results differ depending on whetherthe equation is discrete or continuous. In this context, itshould be borne in mind that it is rather unexplored [41]how the non-fixed point equilibrium sets (like periodic or-bits and chaotic orbits) correspond to the game theoreticequilibria (like Nash equilibria or evolutionarily stablestrategies).It is natural to relate the discrete replicator equationswith the populations having nonoverlapping generations ∗ [email protected] † [email protected] and the continuous one with the overlapping generationscase. But it is not uncommon to see the discrete replica-tor equation for the overlapping generations case [32, 34];in fact, there exists a discrete overlapping-generationsreplicator equation [42] that spans the entire range ofcases from strictly nonoverlapping case to overlapping (todifferent degrees) case. Nevertheless, the discrete repli-cation equations that are not for generation-wise strictlynonoverlapping populations, in effect, means that thestate of the population is observed after every finite in-terval of time—a factor that explicitly enters the discreteequations and when it goes to zero, yields a continuousreplicator equation. A continuous differential equationcan always be discretized in time (e.g., when solving itnumerically) to yield a discrete difference equation thathas the exact the same interpretation: It quantifies theevolution of the state of the system that is observed aftera finite interval of time that can suitably be taken as thetime interval between two successive generations. How-ever, this exercise may mask the fundamental fact thata discrete equation is a model on its own; the profoundimportance of the logistic map [43]— vis-a-vis its contin-uous version [44]—in the development of the theory ofchaos is there for all to see. II. INTRODUCTION
In line with the opening sentence of this paper, a veryuseful discrete replicator equation [9, 23, 41, 45–49] is asfollows: x ( k +1) i = x ( k ) i + x ( k ) i [( Π x ( k ) ) i − x ( k ) T Π x ( k ) ] . (1)It models the evolution of a population with n differ-ent types and x i stands for the fraction of i th phenotypefrom the finite countable set of phenotypes. Π is the pay-off matrix of the corresponding population game and x is the state vector. The superscripts in brackets on x denote the time step or generation number, and super-script T represent the transpose operation. This is thesimplest known replicator equation that is in compliancewith the folk theorem, and simultaneously gives rise toperiodic and chaotic orbits for the simple case of × a r X i v : . [ q - b i o . P E ] F e b normal form games [47, 48]. Therefore, it is a very usefulstylized model that can be investigated to find the con-nections of the periodic and the chaotic orbits with thepurely game theoretic concepts. Henceforth, for the sakeof convenience, we call it type-I replicator map. Notethat Eq. (1) looks like the Euler forward discretization ofthe continuous (standard) replicator equation [18, 50–55], dx i dt = x i (cid:2) ( Π x ) i − x T Π x (cid:3) . (2)The type-I replicator map (Eq. (1)) should be contrastedwith the more popular discrete replicator equation, whichwe term type-II replicator map [18, 19, 50, 56, 57], givenbelow: x ( k +1) i = x ( k ) i ( Π x ( k ) ) i x ( k ) T Π x ( k ) . (3)The continuous dynamic, emanating from this discreteform, leads to the following equation: dx i dt = x i (cid:2) ( Π x ) i − x T Π x (cid:3) x T Π x , (4)which is known as the adjusted replicator equation [19,46].It is very well known [34, 39, 42, 56, 58, 59] thatparticular results, e.g., stability and convergence in thetemporal evolution, of a discrete replicator dynamic arenot necessarily implied by the corresponding continuousreplicator dynamic. The topologies of the solutions neednot be same even if we take two different discrete repli-cator dynamics, viz., the type-I and the type-II repli-cator dynamics. We must recall [19] that topologicallythe solutions of both the continuous replicator equationand the adjusted replicator equation are equivalent forsymmetric normal form games but it need not be so forasymmetric games with bimatrix normal form. Asym-metric games are frequently realized in natural scenarios,e.g., inter-species interactions, interactions between sub-populations, and intra-species interactions due to assign-ment of social roles [60]. A few specific real life examplesof inter-role bimatrix games are: competition betweenthe owner and an invader of a territory [61, 62], the costassociated with sex for a male and a female [63, 64], anddonation game between parent and offspring [65]. Thus,the comparative cases of the corresponding discrete dy-namics in the bimatrix games excites even more curiosity.The microscopic basis of the evolutionary game the-ory relies on the probabilistic behavioural model [66].Both the continuous and the adjusted replicator equa-tions have been shown to be the governing equations inthe large population limit in a generation-wise overlap-ping population where a Moran process with respectiveselection mechanisms is in play [55]. The question weare primarily asking in this paper is the following: Un-der what conditions can the discrete replicator equations(the type-I and the type-II replicator maps) be obtainedas a large population limit of a generation-wise nonover-lapping population where the evolutionary dynamics is described by a microscopic stochastic birth-death pro-cess. In view of the above mentioned subtlety of thebimatrix games, we seek the answer to this question inthe more general case of the bimatrix games.Since it is clearly essential to start with a stochasticbirth-death model that captures nonoverlapping gener-ations, our choice is to work with the seminal Wright–Fisher (WF) model [67, 68] that describes a biologi-cal population with strictly nonoverlapping generations.Here the reproduction process is synchronous, i.e., all re-produce at the same time. Next generation is chosen fromthe pool of these newly born offsprings. There are manyreal life examples of generation-wise non-overlappingpopulations of organisms such as periodical cicada [69],annual plants [70], pink salmon [71], and squid [71]. Theoriginally proposed WF model is devoid of any selectionmechanism. Mathematical structure of the genetic driftin WF model is known in detail [72]. Moreover, the sta-tistical aspects of the WF process with a selection mech-anism incorporated is also well-studied [73, 74].The WF model has been successful in finding some ofthe statistical features of real life biological systems, e.g.,the effective population size in molecular evolution [75],the population divergence time in chimpanzee [76], andthe evolution of drug resistance of the influenza virus [77].Moreover, a recent paper [78] has shown that the WFprocess can model the genomic data of recent past quiteaccurately. These are remarkable practical applicationsof the WF model that, in essence, is a rather simplemathematical model of the birth-death process in a non-overlapping population.In this paper, we show that starting with a selectiondriven WF process—appropriately modified to modeltwo-strategy bimatrix games—we can reach to the type-I and the type-II replicator maps as the governing de-terministic dynamics in the infinite population limit.The WF process, being a discrete time Markov pro-cess [79, 80], facilitates the use of the standard linksamong the master equation, the Fokker–Planck equation,and the Langevin equation to arrive at the maps. III. THE MODEL
We consider a haploid population of size N consistingof individuals with two types of roles denoted by α and β .The entire collection of individuals of α role constitutea subpopulation that we term α -population of size N α and, remaining N β = N − N α individuals constitute β -population. Now, being confined to the one-locus-two-allele theory, the state of the entire population is givenby the allele frequency at the locus under consideration;any individual (of either role) has one of the two alleles,either A or B , at the locus. Hence, for the purpose ofour calculations, an individual can be fully specified byan allele or (pheno-)type, and a role.We denote the state of θ -population (here and hence-forth, θ ∈ { α, β } ) at generation k as i θ which stands forthe number of alleles belonging to type A and role θ .Hence, the joint state of the system can be representedas ( i α , i β ) . For θ -population with a total of N θ num-ber of alleles, i θ can have any value between and N θ .Let the probability of choosing type A individual from θ -population as a parent is given by p θ ; consequently, theprobability of choosing type B is − p θ . As in the case inthe standard WF model, the transition probability T θi θ ,j θ for θ -population from the current state i θ to state j θ inthe next generation is taken as T θi θ ,j θ = (cid:18) N θ j θ (cid:19) ( p θ ) j θ (1 − p θ ) N θ − j θ ; θ ∈ { α, β } . (5)Thus, it is a microscopic stochastic evolutionary dynam-ics having synchronous reproduction where the offspringsin a generation of a subpopulation are chosen through abinomial random sampling with replacement from theimmediately preceding generation of the same subpopu-lation. In the absence of any selection, i.e., for the caseof pure genetic drift, p θ = i θ /N θ . A. Incorporation of selection
In the presence of selection, the probability p θ of choos-ing allele A with role θ as parent must be dependent onthe fitness of the individuals with role θ . In the paradigmof natural selection, the fitness of a type is measured bythe expected reproductive growth rate of that type [81].There may be different ways in which p θ may dependon the fitness profile of the whole population. Below wepresent two simplest logical possibilities:1. In the most frequently used form [73, 82], p θ is theratio of effective reproductive fitness of all individu-als belonging to allele A of θ -population to the effec-tive reproductive fitness of the whole θ -population: p θ = i θ f Aθ i θ f Aθ + ( N θ − i θ ) f Bθ ; θ ∈ { α, β } . (6)Here, f Aθ and f Bθ are respectively the fitnesses oftype A and type B individual in θ -population. Itis clear from Eq. (6) that states ( i α , i β ) , where i α ∈ { , , · · · , N α − } and i β ∈ { , , · · · , N β − } ,are transient states; and for i α ∈ { , N α } and i β ∈ { , N β } , the ( i α , i β ) are absorbing states. Ifthe initial state of the system is in any one of thetransient states then it ultimately reaches to one ofthe absorbing states in finite time and stays thereforever.2. In the presence of selection, it is of higher prob-ability that an offspring has a type A individ-ual as parent if the fitness of type A is morethan the fitness of type B . This comparison be-tween fitnesses should be done in a pair of ran-domly chosen individuals—drawn sequentially with replacement—when the pair consists of individualsof both types. The probability of choosing such apair is [ i θ /N θ ] [( N θ − i θ ) /N θ ] and the comparisonbetween fitnesses can be quantified by f Aθ − f Bθ .In conclusion, p θ (which is only i θ /N θ during ex-clusively random genetic drift) gets augmented bya contribution from the selection and mathemati-cally it may be written down as follows to provideus with an alternate choice: p θ = i θ N θ + (cid:18) i θ N θ (cid:19) (cid:18) N θ − i θ N θ (cid:19) (cid:18) f Aθ − f Bθ ∆ f max θ (cid:19) ; (7) θ ∈ { α, β } . Here the positive proportionality con-stant ∆ f max θ must be chosen in such a way that p θ always stays non-negative and not more than unity.One possible choice is that ∆ f max θ be the maximumpossible value of absolute difference between the fit-nesses of the two types in the θ -population. It canbe easily noted that the transient states and the ab-sorbing states are same as the ones correspondingto Eq. (6). One can easily verify that if f Aθ = f Bθ for θ = α and θ = β , then we get back the standardWF model with no selection and only drift. This istrue for the other case (Eq. (6)) as well.The next question is how exactly one should quan-tify the fitnesses used above? This reproductive fitnessis derived by comparing how one type with one specificrole performs against the subpopulation of the other rolewhen randomly matched. Independent of this interac-tion, there may exist a baseline fitness. Hence, an ef-fective reproductive fitness of the individuals should bedefined as, f Aθ = 1 − w + wπ Aθ , (8a) f Bθ = 1 − w + wπ Bθ , (8b)where θ ∈ { α, β } and w is the strength of selectionweighing the contributions, π Aθ and π Bθ (respectively forthe type A and the type B individuals in θ -population),owing to the inter-subpopulation interaction. Such anexpression for fitness assumes that w is the probabilitythat the fitness (i.e., the expected growth rate) is solelyderived from the inter-role interactions whereas − w isthe probability that only the baseline fitness (fitness inthe absence of any interaction) is contributing. For thepurpose of simplicity, the baseline fitness is assumed tobe independent of the type; by construction, it is samefor both the types and has no role in selection. Hence,the contribution towards the selective advantages (whichdrive the natural selection process) are solely derivedfrom the fitness term due to the inter-role interactionswhich is weighted by w . Hence, the choice of the phrase—strength of selection—to describe w is apt. B. Frequency-dependent selection
While the fitness of an individual of one subpopulationdeceptively seems to be independent of the state of theother subpopulation, it should definitely be not so. Thefitness must depend on the environment that an individ-ual is in. All other individuals in the entire populationeffectively act as an environment. The evolutionary gametheory is arguably an appropriate and efficient techniqueto grasp the essence of frequency dependent selection [81].For a matrix game—having a finite number of strategies(or types)—each payoff element stands for the contribu-tion towards the growth rate of a type resulting from ainteraction corresponding to a particular strategy pro-file. Thus, we resort to the standard path of invokinggame-theoretic consideration that naturally suits such ascenario.The inter-subpopulation interactions, which we haveconsidered up to now, are purely asymmetric in nature.The information about these interactions is contained inthe following two-strategy bimatrix normal form game:Role βA B
Role α A a α , a β b α , c β B c α , b β d α , d β where the first element and the second element in eachbox stand for the payoffs of an individual from α -population and β -population respectively. The respec-tive payoff matrices are more explicitly and compactlygiven below: Π θ = (cid:20) a θ b θ c θ d θ (cid:21) ; θ ∈ { α, β } . (9)If we assume that every individual interacts with a sam-ple of individuals in the opposite role, then the averagepayoff of the individuals must be functions of the frac-tions of the two alleles in the other sub-population. Inother words, if we consider the sample to be the entiresub-population in the opposite role, then the average pay-off per interaction can be written down as, (cid:20) π Aα π Bα (cid:21) = Π α (cid:20) i β /N β − i β /N β (cid:21) , (10a) (cid:20) π Aβ π Bβ (cid:21) = Π β (cid:20) i α /N α − i α /N α (cid:21) . (10b)In the light of these expressions, it is now explicit thatthe fitness of an individual (see Eqs. (8)) of one sub-population is completely dependent on the state of theother sub-population as it should be. It is further clearthat the strength of selection, w , decides the relative con-tribution between baseline fitness and the fitness fromthe game theoretic interactions [83]. One may note that,mathematically, the case w < is equivalent to the case w = 1 with a different payoff matrix; the payoff elementsof the resultant payoff matrix obviously depends both on the underlying game and the strength of selection w [84].This expression of game-theoretic fitness has been suc-cessfully used to model many real life events such as theevolution of human intestinal microbiota [85] and chim-panzee choice rates in competitive games [86].Equipped with the information about how the differ-ential reproduction modelled by the differential fitnessesleads to (natural) selection mediated evolution, we nowturn our attention towards writing the master equationfor the (modified) WF process that is a discrete timeMarkov process. We are driven by the expectation thatin the limit of infinite population, there must be a govern-ing dynamics in discrete time that captures the essenceof the evolutionary game under consideration. IV. TOWARDS DETERMINISTIC DYNAMICS
The master equation [79, 80] of the WF process we aredealing with is a stochastic process and can be writtenfollows: P ( k +1) i θ − P ( k ) i θ = j θ = N θ (cid:88) j θ =0 P ( k ) j θ T θj θ ,i θ − j θ = N θ (cid:88) j θ =0 P ( k ) i θ T θi θ ,j θ , (11)where θ ∈ { α, β } and P ( k ) i θ stands for the probabilitythat the θ -population is in the state i θ at generation k . Subsequently, we rescale the system variables of the θ -population to define the new variables: x θ = i θ /N θ , ˜ x θ = j θ /N θ , and t θ = k/N θ . The corresponding prob-ability density thus reads as ρ θ ( x θ , t θ ) = N θ P ( k ) i θ . Fur-thermore, for the convenience of mathematical manipu-lations, we write the rescaled transition matrix as a func-tion of the state immediately before jump and the jumpsize: (cid:101) T θ (˜ x θ , r θ ) = T θN θ ˜ x θ ,N θ x θ , where r θ = x θ − ˜ x θ . Thus,finally Eq. (11) can be expressed as follows: ρ θ ( x θ , t θ + N − θ ) − ρ θ ( x θ , t θ ) = x θ − (cid:88) r θ = x θ [ ρ θ ( x θ − r θ , t θ ) × (cid:101) T θ ( x θ − r θ , r θ ) (cid:105) − x θ − (cid:88) r θ = x θ ρ θ ( x θ , t θ ) (cid:101) T θ ( x θ , − r θ ) , (12)where θ ∈ { α, β } .Progressing further with Eq. (12) is quite challeng-ing [72, 74, 87, 88] because the transition matrix is nottridiagonal (e.g., the ones for the random walk and theMoran process) for the case of the binomial sampling.However, our primary aim in this paper being the deriva-tion of the corresponding deterministic dynamic, thescope of our paper lies in the limit, N θ → ∞ . This isvery fortunate because of the fact [89] that the binomialdistribution is rigorously approximated by a Gaussiandistribution for any p θ in this limit and the analyticaltractability rendered by the Gaussian distribution is veryconvenient for the required mathematical manipulations.Thus, to being with, we approximate the transition ma-trix when N θ → ∞ as follows: (cid:101) T θ (˜ x θ , r θ ) → G ( N θ x θ ) = 1 (cid:112) πσ θ exp (cid:20) − ( N θ x θ − µ θ ) σ θ (cid:21) , ⇒ (cid:101) T θ (˜ x θ , r θ ) → (cid:112) πσ θ exp (cid:20) − ( N θ ˜ x θ + N θ r θ − µ θ ) σ θ (cid:21) , ⇒ (cid:101) T θ ( x θ , r θ ) → (cid:112) πσ θ exp (cid:20) − ( N θ x θ + N θ r θ − µ θ ) σ θ (cid:21) = G ( N θ r θ ) , (13)to the leading order as other higher order terms arecomparatively negligible. Here, G ( N θ x θ ) is the proba-bility density function of the Gaussian distribution in N θ x θ with mean µ θ = N θ p θ and standard deviation σ θ = (cid:112) N θ p θ (1 − p θ ) . From Eq. (13) we note that fora fixed value of x θ , G ( N θ x θ ) can be equivalently seen asthe Gaussian distribution ( G ( N θ r θ ) being the probabil-ity density function) in N θ r θ with mean µ θ − N θ x θ andstandard deviation σ θ .With this in mind, in Eq. (12), we expand the probabil-ity densities in the Taylor series about t θ in the left handside and about x θ in the right hand side; and similarlyexpand the transition probabilities about x θ to finallyarrive at the Fokker–Planck equation, ∂ρ θ ∂t θ = − ∂∂x θ (cid:34) ρ θ (cid:40) lim N θ →∞ (cid:90) + N θ − N θ r θ G ( N θ r θ ) d ( N θ r θ ) (cid:41)(cid:35) + 12 ∂ ∂x θ (cid:34) ρ θ (cid:40) lim N θ →∞ (cid:90) + N θ − N θ r θ G ( N θ r θ ) d ( N θ r θ ) (cid:41)(cid:35) , (14)having ignored the higher order terms in N − θ . ρ θ isshorthand for ρ θ ( x θ , t θ ) . Here we have made use ofEq. (13) to replace (cid:101) T θ ( x θ , r θ ) with its correspondingGaussian approximation G ( N θ r θ ) and consequently, thesummations become integrals. The validity of such diffu-sion approximation to derive the Fokker–Planck equationis well investigated in population genetics [90]. Assumingthe stochastic noise to be uncorrelated in time, we canuse the It ˆo calculus to reach the following form of theLangevin equation for θ ∈ { α, β } [79, 80, 91, 92], x ( k +1) θ − x ( k ) θ = lim N θ →∞ (cid:90) + N θ − N θ r θ G ( N θ r θ ) d ( N θ r θ )+ ξ θ (cid:115) lim N θ →∞ (cid:90) + N θ − N θ r θ G ( N θ r θ ) d ( N θ r θ ) , (15)where ξ θ is the Gaussian white noise. While writingEq. (15), we have once again reminded ourselves thatthe WF model fashions the birth-death process in a pop-ulation with nonoverlapping generations, and hence thediscretized version of the Langevin equation is what weshould use. The stochastic term has / √ N θ dependance,implying that only the deterministic term contributes forthe infinite population. Therefore, the discrete dynam-ics for the evolution of the states of θ -population in the very large, generation-wise nonoverlapping population isgiven as x ( k +1) θ = p θ ( x ( k ) α , x ( k ) β ); θ ∈ { α, β } . (16)Here, considering Eqs. (6)–(10), we have explicitly high-lighted that p θ is a function of both x α and x β .Subsequently, in the light of Eq. (8), Eq. (10)and Eq. (16), the two different models of selection—viz., Eq. (6) and Eq. (7) as discussed in Sec. III A—respectively yield x ( k +1) α = x ( k ) α (cid:16) Π α x ( k ) β (cid:17) x ( k ) Tα Π α x ( k ) β , (17a) x ( k +1) β = x ( k ) β (cid:16) Π β x ( k ) α (cid:17) x ( k ) Tβ Π β x ( k ) α ; (17b)and x ( k +1) α = x ( k ) α + x ( k ) α ∆ f max α (cid:104) ( Π α x ( k ) β ) − x ( k ) Tα Π α x ( k ) β (cid:105) , (18a) x ( k +1) β = x ( k ) β + x ( k ) β ∆ f max β (cid:104) ( Π β x ( k ) α ) − x ( k ) Tβ Π β x ( k ) α (cid:105) . (18b)The subscript, , denotes the first row of the correspond-ing column vectors. Thus, we have successfully provideda microscopic basis for the discrete versions of the stan-dard and the adjusted continuous replicator equationsfor the case of a two-strategy bimatrix population gamein a generation-wise nonoverlapping population; we alsovalidate our findings numerically in FIG. 1.In passing, through FIG. 2, we emphasize once morethat the replicator maps are good approximations of thegoverning dynamics for the selection driven WF processesonly in the large population limit—note that the stochas-tic contribution is inversely proportional to the squareroot of the population size. The presence of stochastic-ity, in the finite population, leads the population to fixatein any one of its absorbing states—the sub-populationseventually fixate either into all A or into all B players—with a non-zero probability. However, it is a challeng-ing task to find these fixation probabilities analyticallyfor a selection driven WF process because of the non-tridiagonal form of the transition matrix. In future, itmay be insightful to investigate if one can make use ofthe available approaches [79, 93–95]—using either thediscrete master equation or the corresponding Fokker–Plank equation—adopted for finding the fixation proba-bility of a Moran process involving symmetric games (oreven asymmetric games [96]).Coming back to Eqs. (18a)–(18b), please note the ex-plicit presence of the positive definite quantities— ∆ f max α and ∆ f max β — that, at first glance, may be thought to beabsorbable into the respective payoff matrices. However, n . . . . . . . . x α (a) Type-I map N α = N β =1000 n . . . . . . x β (b) Type-I map N α = N β =1000 n . . . . . . . . x α (c) Type-II map N α = N β =1000 n . . . . . . x β (d) Type-II map N α = N β =1000 FIG. 1. Two dimensional replicator dynamics, type-I bi-matrix dynamics (Eq. (18)) and type-II bimatrix dynamics(Eq. (17)), are mean field dynamics of the WF process inthe very large population limit. Subplots (a) and (b) exhibitthe dynamics of α -population and β -population respectively;black solid lines are the evolution under Eq. (18) and blackdashed lines with red error bars are the corresponding simu-lated WF process averaged over 50 trials for N α = N β = 1000 .Subplots (c) and (d) are the analogous plots correspondingto the other case, viz., Eq. (17). For illustrative purposes,we consider that the α -population is playing Harmony game( a α = 1 , b α = 1 , c α = 0 . , and d α = 0 ) in both the caseswhereas the β -population is playing Prisoner’s Dilemma game( a β = 1 , b β = − , c β = 0 . , and d β = 0 ) and Leadergame ( a α = 1 , b α = 5 , c α = 6 and d α = 0 ) respectivelyin the cases corresponding to the type-I map and the type-II map. In subplots (a) and (b), the initial conditions cho-sen is ( x α = 0 . , x β = 0 . while in subplots (c) and (d) ( x α = 0 . , x β = 0 . is chosen as the initial condition. there is an interesting caveat that we discuss now. It isstraightforward to see that if we remove the concept ofthe roles in the population and allow for game-theoreticinteraction of any two individuals in the resultant pop-ulation of size N (= N α + N β ) in the form of a × normal form game with payoff matrix, Π = Π α = Π β ,we get back Eq. (1) and Eq. (3). (This leads to a one-dimensional dynamics in a simplex, Σ , in contrast witha two-dimensional dynamics in Σ × Σ accessed by thebimatrix game.) Specifically, we arrive at the type-I repli-cator map (cf. Eq. (1)), x ( k +1) = x ( k ) + 1∆ f max x ( k ) [( Π x ( k ) ) − x ( k ) T Π x ( k ) ] , (19)where ∆ f max is the maximum possible value of absolutedifference between the fitnesses of the two types in thepopulation with no roles. This equation, in contrast toEq. (1), is incapable of showing any non-fixed point out- n . . . . . . x α (a) N α = N β =10 n . . . . . . x β (b) N α = N β =10 n . . . . . . x α (c) N α = N β =10 n . . . . . . x β (d) N α = N β =10 FIG. 2. For a finite population size, the effect of stochasticityleads the sub-population to randomly fixate in one of the twotypes, all A or all B . Subplots (a) and (b) exhibit the dynam-ics of the α -population and the β -population respectively forthe same underlying games and initial conditions as used inFIG. 1(a) and FIG. 1(b); each coloured line is the simulatedWF process generated in a random trial for N α = N β = 10 .Similarly, subplots (c) and (d) exhibit the same scenario forthe α -population and the β -population respectively for theunderlying games and initial conditions used in FIG. 1(c) andFIG. 1(d). comes (like periodic orbits and chaotic orbits); further-more, unlike Eq. (1), any real payoff matrix can be usedin Eq. (19) without violating the constraint x ∈ [0 , atany time. Details are presented in Appendix A. V. DISCUSSION AND CONCLUSION
We emphasize that similar derivation [55] of the stan-dard and the adjusted continuous replicator equationsfrom the microscopic birth-death process is based onthe Moran model that models birth-death process ina generation-wise overlapping population; but for thegeneration-wise nonoverlapping population, fundamen-tally different birth-death process—e.g., the WF model—must be used. It is not obvious a priori that the deter-ministic equations in the large population limit in thetwo birth-death processes should be connected. Conse-quently, in this paper, we have undertaken to analyticallyprove that the mean field dynamics that approximatesthe selection-driven evolution in a microscopic stochasticbirth-death process in a population, where generationsare nonoverlapping, can indeed be modelled by replica-tor maps. We have found that the type-I replicator mapand the type-II replicator map, which are dynamicallydrastically different [48] even in the simplest case of × normal form games, are the manifestations of the differ-ent ways in which individual fitness affects the selectionin the birth-death process.It is interesting to note that the mean-field dynam-ics of a selection driven WF process is the discretizedversion of the mean-field dynamics corresponding to aselection driven Moran process [55] although these twoprocesses have two different microscopic justifications.The literature of these two microscopic processes has of-ten highlighted this similarity; in fact, it is known thatfor the limit of infinite population many of the statisti-cal and genealogical properties of these two microscopicprocesses are approximately same [74, 97]. The modelof selection—mediated through the choice of p θ —plays acrucial role in this aforementioned similarity of the mean-field dynamics. Note that the increment of the fraction ofa type in the offspring pool of the WF process is similarto the excess probability that this type is chosen as anoffspring in the corresponding Moran process [55]. Thus,the reason behind the similarity in the mean-field dynam-ics for these two processes is the particular incorporationof selection in the neutral model. However, this similar-ity is only for the mean-filed dynamics—the stochasticcomponents of the two processes behaves differently (seeAppendix B). We remark that a different form of p θ couldgive rise to an altogether different mean field dynamicsthat may be investigated in future; in this paper, our goalhas been to find a microscopic model that leads to thewidely used versions of the replicator maps.We remark that the WF model has the known advan-tage over the Moran model in allowing for diploidy in thepopulation; in fact, one can easily redo the entire exerciseof this paper by considering a monoecious diploid pop-ulation to arrive at the same replicator maps. It mustbe explicitly highlighted that instead of type-I equationgiven by Eq. (1), the above mentioned microscopic pro-cess yields a slightly different Eq. (19) that is very attrac-tive as it always gives rise to convergent fixed-point out-comes. From a different perspective, it also means thatthe problem of deriving (from some microscopic birth-death process) the type-I equation given by Eq. (1) whichcan yield chaotic outcomes [47, 48] remains an open ques-tion. From our present study, it appears that one needsto go beyond both the Moran and the WF models tosolve this problem.While relatively less in vogue in the biological systems,Eq. (1)—as itself or in related forms—also appears inmodelling intergenerational cultural transmission [9, 45],boundedly rational players’ imitational behaviour in bi-matrix cyclic games [46], and reinforcement learning [23].It is interesting to recall that a good behaviour rule doesnot require aggregate population behaviour implicit inthe natural selection to induce the replicator map; in-stead the map is arrived at based on the rational be-haviour of the players [54]. In such contexts, the ‘strate-gies’ A and B are not hardwired in the individuals andthey are free to choose strategies. In such a scenario,we can put w = 1 , i.e., the interactions are only game- theoretic and note that the selection process in the type-I replicator map (in contrast with the type-II replicatormap) intriguingly requires the underlying microscopic in-teractions to depend only on location information (seeEq. 7): At every time step a pair of randomly chosen in-dividuals with opposite strategies compare each other’sexpected payoffs; in other words, the selection dependson f Aθ − f Bθ = π Aθ − π Bθ . (The other selection model re-quires the knowledge of the global information for thecalculation of the average payoff; see the denominator ofthe right hand side of Eq. (6).)Some of the basic driving forces of evolution are se-lection, mutation, drift, migration, and recombination.All the replicator equations discussed herein are con-cerned with selection driven evolution. The determin-istic replicator equations, by construction, can not cap-ture the effect of drift. Migration is best modelled intothe equations after imparting some spatial structure inthe system. Extensions of the replicator equations toincorporate the effect of mutation has been extensivelyresearched [53, 98–104]; mathematically, one can mod-ify the fitness functions and derive different replicator-mutator maps in a generation-wise nonoverlapping popu-lation of haploid (or monoecious diploid) individuals (seeAppendix C). The mechanism of recombination, and re-lated effects, are beyond the scope of the present inves-tigation; but they could be a potential avenue of futureresearch in the context of this paper.It is interesting to note that while the two different se-lection mechanisms leads to two different replicator mapsin the infinite population limit, some of their effects onthe finite population are identical. For example, considerthe simpler case of the symmetric normal form × games with no different roles: If the selection favourstype A individual in invading a population with type B individuals and opposes the opposite scenario, then forevery possible transient initial state, the probability that A fixates is more than for the case where only drift (noselection) is in action. This follows [73] from the fact thatthe submatrix formed by deleting the first and the lastcolumns and rows of the transition matrix (correspond-ing to any of the two selection mechanisms) is positive-definite, and its every row-sum—consequently, the spec-tral radius—is less than unity. Similarly, one may alsobe interested in comparing between the two update rules,in terms of the stochastic gains due to the intrinsic noisearising as a result of finiteness of the population size [105].More exhaustive studies on the comparison between thetwo selection mechanisms for the case of finite populationmay be worth pursuing in future.Before we end, let us remind ourselves that which repli-cator equation is a better model of evolution is an is-sue open for debate. One of the earliest hint aboutthis non-unanimity may be found in the seminal bookby Maynard-Smith [19] where it is discussed whetherthe adjusted replicator dynamics is a better model thanthe standard replicator dynamics, given the fact thatthese two continuous dynamics cease to have topologi-cally identical solutions for the asymmetric games. How-ever, one must take this non-unanimity positively andconstructively since the replicator equations are merelyminimal models that easily incorporate the Darwiniantenet of natural selection to investigate mathematicalsystem mimicking exclusively the replication-selection as-pect of the real evolutionary systems. In this spirit, whatwe have achieved in this paper is to provide a more sound(microscopic) footing to the type I replicator map thathas immense potential to serve as a convenient testbed forinvestigating the connections [41, 48] between the com-plex dynamical outcomes and game-theoretic equilibria. ACKNOWLEDGMENTS
The authors are grateful to Jayanta K. Bhattachar-jee, Debashish Chowdhury and Samrat Sohel Mondal forhelpful discussions.
AIP PUBLISHING DATA SHARING POLICY
The data that support the findings of this study areavailable from the corresponding author upon reasonablerequest.
Appendix A: No Non-Fixed Point Outcomes inEq. (19)
For the purpose of this section and without any lossof generality, it is sufficient to work with the followingform of the payoff matrix [48]:
Π = ( ST ); S, T ∈ R .For further convenience, we explicitly define a function g ( x ( k ) ) as follows to rewrite Eq. (19) as x k +1 = g ( x ( k ) ) = x ( k ) + 12 H x ( k ) (cid:18) Cx ( k ) + S ∆ f max (cid:19) , (A1)where C = 1 − T − S and H x ( k ) = 2 x ( k ) (1 − x ( k ) ) . Onemay note that, by definition, ∆ f max = max x | Cx + S | = max {| S | , | − T |} . (A2)We immediately note that for all S, T ∈ R in Eq. (A1), x ( k ) ∈ [0 , for all k ∈ N ∪ { } because ≤ | ( Cx ( k ) + S ) / ∆ f max | ≤ and ≤ H x ( k ) / ≤ (1 − x ( k ) ) ≤ for all x ( k ) ∈ [0 , and for all S, T ∈ R .Now to prove that no periodic orbit is possible inEq. (A1), we only have to show that g ( x ) − x and g m ( x ) − x ( m ∈ N ), where g m ( x ) is the m th iterate of g ( x ) , have same sign for all x ∈ [0 , \ F ; here, F is the setof all fixed points or 1-period points, i.e., the solutionsof g ( x ) = x . (Note that, in the immediately precedingparagraph, we have already shown that g m ( x ) ∈ [0 , forall x ∈ [0 , .) This is so because it means that in theplot of g ( x ) versus x , the number of intersections made with the line g ( x ) = x is same as the number of inter-sections made between the line and the plot of g m ( x ) versus x ; in other words, no prime m -period point (nonew intersection) can appear.To this end, we first define δ m = g m ( x ) − x and useEq. (A1) to find δ m +1 = δ m (cid:18) H g m ( x ) C ∆ f max (cid:19) + δ H g m ( x ) H x , (A3)where g ( x ) = x . Note that for the case of m = 1 ,Eq. (A3) can be simplified as, δ = δ (cid:18) H g ( x ) C ∆ f max + H g ( x ) H x (cid:19) . (A4)We have already shown that the map, expressed byEq. (A1), is forward invariant: g m ( x ) ∈ [0 , for all m ∈ N . By definition, ≤ H g m ( x ) ≤ / . If C ≥ , it istrivial to see that the term in the parentheses of Eq. (A4)is always positive which implies that δ and δ have samesign. Whenever C < , it is straightforward to see that | C/ ∆ f max | < and hence, (1 / H g m ( x ) | C/ ∆ f max | < / for all m ∈ N . Therefore, δ and δ have same signs forall x ∈ [0 , \ F .Following the similar arguments, it is easy to see fromEq. (A3) that δ and δ m +1 have same sign if δ and δ m have same sign. In the preceding paragraph, wehave shown that δ and δ have same sign, and so fromEq. (A3)—with m = 2 —we conclude that δ and δ havesame sign. Thus, using the method of induction, oneconcludes that δ and δ m have same sign for all m ∈ N and no prime m -period points are possible. Furthermore,since there is no periodic orbit in the system, no chaoticattractor—that is supposed to have countably infinitenumber of unstable periodic orbits—can be present inthe system. Appendix B: A brief note on the stochastic dynamics
The master equation captures the time evolution ofthe probability distribution of a system and has the fullinformation of the stochastic evolution of the model un-der consideration. Thus, in principle, one should be ableto use it to derive the replicator map directly from themaster equation [11]. To this end, we multiply Eq. (11)by i θ and sum over its all possible values to write thefollowing: i θ = N θ (cid:88) i θ =0 i θ P ( k +1) i θ − i θ = N θ (cid:88) i θ =0 i θ P ( k ) i θ = j θ = N θ (cid:88) j θ =0 P ( k ) j θ i θ = N θ (cid:88) i θ =0 i θ T θj θ ,i θ − i θ = N θ (cid:88) i θ =0 i θ P ( k ) i θ j θ = N θ (cid:88) j θ =0 T θi θ ,j θ . (B1)On using the form of transition matrix elements givenby Eq. (5) and the properties of a binomial distribu-tion, we find that (cid:80) i θ = N θ i θ =0 i θ T θj θ ,i θ = N θ p θ ( x ( k ) α , x ( k ) β ) and (cid:80) j θ = N θ j θ =0 T θi θ ,j θ = 1 . Using these relations, we rewriteEq. (B1) as, (cid:104) x ( k +1) θ (cid:105) = (cid:104) p θ ( x ( k ) α , x ( k ) β ) (cid:105) , (B2)where the pair of angular brackets denotes the ensem-ble average (average over the probability distribution).For the case of infinite population size ( N θ → ∞ ) wecan safely approximate that the fluctuation around themean is negligible. Hence, we can remove the angularbrackets used in Eq. (B2) to reach Eq. (16) which furtherleads to Eq. (17) and Eq. (18) depending on the choiceof p θ ( x ( k ) α , x ( k ) β ) . We observe once more that Eq. (B2),being the mean field dynamics, can not explicitly revealthe implications of the stochastic effects that is presentin the system.However, in order to find the stochastic componentof the dynamics analytically, one must resort to furtherapproximations. We note that the master equation un-der consideration has time-independent transition ratesand hence the system it is describing is Markovian in na-ture. Consequently, its Kramer–Moyal expansion up tosecond order (an approximation leading to the Fokker–Planck equation) corresponds to the Langevin equationwith a Gaussian delta-correlated noise [106]. Of course, itwould be more realistic to work ab initio with a Langevindynamics having a coloured noise with a finite corre-lation time; however, not only the corresponding non-Markovian dynamics is generally analytically intractablebut also it means that the corresponding master equa-tion has time-dependent transition rates [107] at oddswith the basic premise of the ideal Wright–Fisher pro-cess considered in this paper.Within the approximation of the Markov process—modelled by the Langevin equation with white noise ξ θ —it is easy to find the stochastic component analytically:Calculations simplify Eq. (15) to x ( k +1) θ = p θ ( x ( k ) α , x ( k ) β ) + (cid:115) p θ ( x ( k ) α , x ( k ) β )[1 − p θ ( x ( k ) α , x ( k ) β )] N θ ξ θ . (B3)This clearly indicates that only the mean field dynamics(deterministic part) of the selection driven Moran processand the selection driven WF process are similar while thestochastic parts can be different, e.g., the amplitude ofthe stochastic contribution for a selection driven Moranprocess with local update rule (which leads to the repli-cator dynamics as the mean field dynamics) only dependson the state of the system [55] whereas the same for thecorresponding WF process (which has the type-I map as the mean-field dynamics) depends both on the state ofthe system and the underlying game (see Eq. (B3)). Appendix C: Replicator-Mutator Maps
Within the scope of the replicator equations, the muta-tion effectively means that a particular type of individualreproduces the other type of individual. Let the proba-bility that some A type offsprings are born from an A type individual is given by q θ ( ≤ q θ ≤ ) in the θ -population leaving the possibility that with probability − q θ , some B type offsprings are reproduced by an A type individual. In general, in the case of a populationwith n different types of individual, the probability thatsome of j th type individuals are reproduced by an i thtype individual is given by q ij ( i, j ∈ { , , · · · , n } )— n × n elements of the mutation matrix that is row stochastic.For the θ -population, further assuming that the mutationis symmetric between the A and the B type individuals,we note that the single parameter q θ captures the muta-tion process. The resulting replicator-mutator map canbe obtained by modifying the expressions of the fitnessesas follows: f Aθ = 1 − w + w (cid:2) q θ π Aθ + (1 − q θ ) π Bθ (cid:3) , (C1a) f Bθ = 1 − w + w (cid:2) q θ π Bθ + (1 − q θ ) π Aθ (cid:3) . (C1b)Following the method detailed in the text, the corre-sponding replicator-mutator maps, viz., type I and typeII, for the two-strategy bimatix games are obtained re-spectively as: x ( k +1) α = x ( k ) α + q α (cid:16) Π α x ( k ) β (cid:17) x ( k ) α + (1 − q α ) (cid:16) Π α x ( k ) β (cid:17) × (1 − x ( k ) α ) − x ( k ) α (cid:104) x ( k ) Tα Π α x ( k ) β (cid:105) , (C2a) x ( k +1) β = x ( k ) β + q β (cid:16) Π β x ( k ) α (cid:17) x ( k ) β + (1 − q β ) (cid:16) Π β x ( k ) α (cid:17) × (1 − x ( k ) β ) − x ( k ) β (cid:104) x ( k ) Tβ Π β x ( k ) α (cid:105) ; (C2b)and x ( k +1) α = q α (cid:16) Π α x ( k ) β (cid:17) x ( k ) α + (1 − q α ) (cid:16) Π α x ( k ) β (cid:17) (1 − x ( k ) α ) x ( k ) Tα Π α x ( k ) β , (C3a) x ( k +1) β = q β (cid:16) Π β x ( k ) α (cid:17) x ( k ) β + (1 − q β ) (cid:16) Π β x ( k ) α (cid:17) (1 − x ( k ) β ) x ( k ) Tβ Π β x ( k ) α . (C3b) [1] C. Hauert and G. Szabo, Am. J. Phys. , 405 (2005). [2] E. G. Hidalgo, Physica A: Statistical Mechanics and itsApplications , 393 (2006). [3] D. Helbing and A. Johansson, Phys. Rev. E , (2010).[4] S. Sadhukhan, R. Chattopadhyay, and S. Chakraborty,Phys. Rev. Research , 013009 (2021).[5] P. F. Stadler, BioSystems , 1 (1991).[6] B. M. R. Stadler and P. F. Stadler, Adv. Complex Syst. , 47 (2003).[7] S. Shahshahani, A New Mathematical Framework forthe Study of Linkage and Selection (Mem. Amer. Math.Soc., Volume 17, Number 211, AMS, Providence, RI,1979).[8] J. Wang, F. Fu, T. Wu, and L. Wang, Phys. Rev. E ,(2009).[9] J. D. Montgomery, Am. Econ. J.: Microeconomics ,115 (2010).[10] J. S. Weitz, C. Eksin, K. Paarporn, S. P. Brown, andW. C. Ratcliff, Proc. Natl. Acad. Sci. U.S.A. , E7518(2016).[11] Y.-H. Lin and J. S. Weitz, Phys. Rev. Lett. , 148102(2019).[12] A. R. Tilman, J. B. Plotkin, and E. Akcay, Nat. Com-mun. , 915 (2020).[13] E. V. Damme, Eur. Econ. Rev. , 847 (1994).[14] J. W. Weibull, Eur. Econ. Rev. , 868 (1994).[15] L. Samuelson, Evolutionary Games and Equilibrium Se-lection , 1st ed., MIT Press Books, Vol. 1 (The MITPress, 1998).[16] D. Friedman, J. Evol. Econ. , 15 (1998).[17] J. Newton, Games , 31 (2018).[18] J. Hofbauer and K. Sigmund, Evolutionary Gamesand Population Dynamics (Cambridge University Press,Cambridge, 1998).[19] J. M. Smith,
Evolution and the Theory of Games (Cam-bridge University Press, Cambridge, 1982).[20] M. A. Nowak and K. Sigmund, Science , 793 (2004).[21] K. Tokita, Phys. Rev. Lett. , 178102 (2004).[22] L. David and T. T. D, Interface Focus , (2014).[23] T. Börgers and R. Sarin, J. Econ. Theory , 1 (1997).[24] Y. Sato, E. Akiyama, and J. D. Farmer, Proc. Natl.Acad. Sci. U.S.A. , 4748 (2002).[25] Y. Sato and J. P. Crutchfield, Phys. Rev. E , 015206(2003).[26] A. Traulsen, T. Röhl, and H. G. Schuster, Phys. Rev.Lett. , 028701 (2004).[27] T. Galla, Phys. Rev. Lett. , 198702 (2009).[28] A. Kianercy and A. Galstyan, Phys. Rev. E , 012815(2013).[29] J. Jost and W. Li, Phys. Rev. E , 022113 (2014).[30] D. Bloembergen, K. Tuyls, D. Hennes, and M. Kaisers,J. Artif. Intell. Res. , 659 (2015).[31] D. Foster and P. Young, Theor. Popul. Biol. , 219(1990).[32] K. Binmore, Fun and Games , 1st ed. (D. C. Heath, Lex-ington, MA, 1992).[33] R. T. Boylan, J. Econ. Theory , 473 (1992).[34] A. Cabrales and J. Sobel, J. Econ. Theory , 407(1992).[35] D. Fudenberg and C. Harris, J. Econ. Theory , 420(1992).[36] R. T. Boylan, J. Econ. Theory , 615 (1995).[37] L. Samuelson, J. Econ. Perspect. , 47 (2002).[38] D. Fudenberg and T. Strzalecki, Econometrica , 651(2015).[39] J. H. Nachbar, Int. J. Game Theory , 59 (1990). [40] R. Cressman and Y. Tao, Proc. Natl. Acad. Sci. ,10810 (2014).[41] A. Mukhopadhyay and S. Chakraborty, J. Theor. Biol. , 110288 (2020).[42] J. W. Weibull, Evolutionary Game Theory , 1st ed.,Vol. 1 (MIT Press Books, The MIT Press, 1997).[43] R. May, Nature , 459 (1976).[44] H.-O. Peitgen, H. Jürgens, and D. Saupe,
Chaos andfractals - new frontiers of science (2. ed.). (Springer,2004).[45] A. Bisin and T. Verdier, Q. J. Econ. , 955 (2000).[46] J. Hofbauer and K. H. Schlag, J. Evol. Econ. , 523(2000).[47] D. Vilone, A. Robledo, and A. Sánchez, Phys. Rev.Lett. , (2011).[48] V. Pandit, A. Mukhopadhyay, and S. Chakraborty,Chaos , 033104 (2018).[49] A. Mukhopadhyay and S. Chakraborty, Chaos: An In-terdisciplinary Journal of Nonlinear Science , 121104(2020).[50] P. D. Taylor and L. B. Jonker, Math. Biosci. , 145(1978).[51] P. Schuster and K. Sigmund, J. Theor. Biol. , 533(1983).[52] P. Schuster and K. Sigmund, Berichte der Bunsenge-sellschaft für physikalische Chemie , 668 (1985).[53] K. M. Page and M. A. Nowak, J. Theor. Biol. , 93(2002).[54] R. Cressman, Evolutionary Dynamics and ExtensiveForm Games , 1st ed., MIT Press Books, Vol. 1 (TheMIT Press, Cambridge, MA, 2003).[55] A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev.Lett. , 238701 (2005).[56] E. Dekel and S. Scotchmer, J. Econ. Theory , 392(1992).[57] E. van Damme, Stability and Perfection of Nash Equi-libria (Berlin: Springer Berlin Heidelberg, 1991).[58] F. J. Weissing, “Evolutionary stability and dynamic sta-bility in a class of evolutionary normal form games,” in
Game Equilibrium Models I: Evolution and Game Dy-namics , edited by R. Selten (Springer Berlin Heidelberg,Berlin, Heidelberg, 1991) pp. 29–97.[59] J. Björnerstedt, M. Dufwenberg, P. Norman, andJ.W.Weibull, “Evolutionary selection dynamics and ir-rational survivors,” in
Understanding Strategic Interac-tion , edited by W. Albers, W. Güth, P. Hammerstein,B. Moldovanu, and E. V. Damme (Springer Berlin Hei-delberg, Berlin, Heidelberg, 1997) pp. 128–148.[60] A. McAvoy and C. Hauert, PLOS Comput. Biol. , 1(2015).[61] O. Leimar and M. Enquist, J. Theor. Biol. , 475(1984).[62] A. Grafen, Anim. Behav. , 462 (1987).[63] T. H. Clutton-Brock, “The evolution of parental care,”(Princeton University Press, Princeton, 1991).[64] M. D. Jennions and L. Fromhage, Philos. Trans. R. Soc.B , 20160312 (2017).[65] J. A. Marshall, J. Theor. Biol. , 386 (2009).[66] D. Helbing, Theor Decis , 149 (1996).[67] W. S., Genetics , 97 (1931).[68] R. A. Fisher, The genetical theory of natural selection (Oxford University Press: Clarendon, 1930).[69] L. D. English, J. J. English, R. N. Dukes, and K. G.Smith, Environ. Entomol , 245 (2006). [70] M. C. Albani and G. Coupland, in Plant Development ,Current Topics in Developmental Biology, Vol. 91,edited by M. C. Timmermans (Academic Press, 2010)pp. 323 – 348.[71] S. Ratner and R. Lande, Ecology , 3093 (2001).[72] T. Tran, J. Hofrichter, and J. Jost, Theor. Biosci. ,73 (2013).[73] L. Imhof and M. Nowak, J. Math. Biol. , 667 (2006).[74] P. Tataru, M. Simonsen, T. Bataillon, and A. Hobolth,Syst. Biol. , e30 (2016).[75] B. Charlesworth, Nat Rev Genet , 195 (2009).[76] P. Tataru, T. Bataillon, and A. Hobolth, Genetics ,1133 (2015).[77] A. Ferrer-Admetlla, C. Leuenberger, J. D. Jensen, andD. Wegmann, Genetics , 831 (2016).[78] D. Nelson, J. Kelleher, A. P. Ragsdale, C. Moreau,G. McVean, and S. Gravel, PLOS Genetics , 1 (2020).[79] C. Gardiner, Handbook of stochastic methods forphysics, chemistry and the natural sciences , 3rd ed.(Springer-Verlag, Berlin, 2004).[80] N. V. Kampen,
Stochastic processes in physics andchemistry , 3rd ed. (North Holland, 2007).[81] S. J. Brown, Proc. R. Soc. B. , 20160847 (2016).[82] A. Traulsen, J. M. Pacheco, and L. A. Imhof, Phys.Rev. E , 021905 (2006).[83] M. A. Nowak, A. Sasaki, C. Taylor, and D. Fudenberg,Nature , 646 (2004).[84] J. Claussen and A. Traulsen, Phys. Rev. E , (2005).[85] A. Wu and D. Ross, Games , 3:26 (2016).[86] C. F. Martin, R. Bhui, P. Bossaerts, T. Matsuzawa, andC. Camerer, Sci. Rep. , 5182 (2014).[87] N. C. Weber, J. Appl. Probab. , 504 (1986).[88] T. D. Tran, J. Hofrichter, and J. Jost, Theory Biosci. , 83 (2015).[89] A. Papoulis, Probability, Random Variables, andStochastic Processes , 3rd ed. (McGraw-Hill, New York,1991).[90] E. Aalto, J. Theor. Biol. , 317 (1989). [91] J. Grasman and O. A. van Herwaarden,
AsymptoticMethods for the Fokker–Planck Equation and the ExitProblem in Applications (Springer-Verlag, Berlin, 1999).[92] Z. Czechowski, Chaos: An Interdisciplinary Journal ofNonlinear Science , 053109 (2016).[93] W. Ewens, Mathematical population genetics (Springer,Berlin, 2004).[94] T. Antal and I. Scheuring, Bull. Math. Biol. , 1923(206).[95] M. Assaf and M. Mobilia, J. Stat. Mech.: Theory Exp. , P09009 (2010).[96] T. Sekiguchi and H. Ohtsuki, Dyn Games Appl , 93(2017).[97] A. Bhaskar and Y. S. Song, Bioinformatics , i187(2009).[98] M. Eigen and P. Schuster, Naturwissenschaften , 541(1977).[99] K. P. Hadeler, SIAM J. Appl. Math. , 1 (1981).[100] P. F. Stadler and P. Schuster, J Math Biol. , 597(1992).[101] I. Bomze and B. Reinhard, Games Econ. Behav. , 146(1995).[102] N. L. Komarova, J. Theor. Biol. , 227 (2004).[103] M. Eigen, J. McCaskill, and P. Schuster, “The molec-ular quasi-species,” in Advances in Chemical Physics (John Wiley & Sons, Ltd, 2007) pp. 149–263.[104] S. Mittal, A. Mukhopadhyay, and S. Chakraborty,Phys. Rev. E , 042410 (2020).[105] T. Röhl, A. Traulsen, J. C. Claussen, and H. G. Schus-ter, Phys. Rev. E , 026108 (2008).[106] M. S. Miguel and R. Toral, in Nonlinear Phenom-ena and Complex Systems , Instabilities and Nonequi-librium Structures VI., Vol. 5, edited by E. Tirapegui,J. Martínez, and R. Tiemann (Springer, Dordrecht,2000) pp. 35 – 127.[107] P. Häunggi and P. Jung, “Colored noise in dynamicalsystems,” in