Abstract.
We demonstrate a vast expansion of the theory of evolutionary stability to finite populations with mutation, connecting the theory of the stationary distribution of the Moran process with the Lyapunov theory of evolutionary stability. We define the notion of stationary stability for the Moran process with mutation and generalizations, as well as a generalized notion of evolutionary stability that includes mutation called an incentive stable state (ISS) candidate. For sufficiently large populations, extrema of the stationary distribution are ISS candidates and we give a family of Lyapunov quantities that are locally minimized at the stationary extrema and at ISS candidates. In various examples, including for the Moran and WrightFisher processes, we show that the local maxima of the stationary distribution capture the traditionallydefined evolutionarily stable states. The classical stability theory of the replicator dynamic is recovered in the large population limit. Finally we include descriptions of possible extensions to populations of variable size and populations evolving on graphs.
Stationary Stability for Evolutionary Dynamics in Finite Populations
Marc Harper^{*}^{*}*corresponding author, email: marc.harper@gmail.com
AMS keywords: 91A22, 92D25, 37B25^{1}, Dashiell Fryer^{2}
^{1}Department of Genomics and Proteomics, UCLA,
^{2}Department of Mathematics, Pomona College
Graphical Abstract

1. Introduction
Finite population evolutionary dynamics are of broad interest for both application and theory [60] [46] [13] [14] [42] [21] [45] [44]; in particular, concepts of evolutionary stability have been widelystudied [62] [59] [69] [66], and recently applied to finite populations [70] with mutation (and weak selection) [5]. Many evolutionary dynamics, such as the replicator equation, effectively assume infinite population sizes which allows a variety of tools from dynamical systems to be applied. Lyapunov stability theorems are crucial tools in the study of dynamical systems and evolutionary dynamics and there is a long history of Lyapunov stability results in evolutionary game theory and mathematical biology [69] [7] [1] [53] [55]. Evolutionary dynamics in finite populations are commonly modeled as Markov processes, such as the Moran and the WrightFisher processes [41] [40] [35] [61] [12], which are not deterministic; in the case of the Moran process, the replicator equation can be recovered in a large population limit [63] [64]. Such stochastic models lack Lyapunov stability theorems. The theory of Markov processes, however, provides powerful analytic tools, notably the stationary distribution, capturing the longrun behavior of a process. We will show that there is an intimate connection between the stationary distributions of these processes and the celebrated theory of evolutionary stability.
In particular, for populations of finitely many replicating types, we define a Lyapunovlike quantity, namely a quantity that is positivedefinite, decreasing toward an “equilibrium” locally, and minimal at the equilibrium, for a large class of birthdeath processes we call the incentive process [27], including the Moran process and the Fermi process [65]. We then show local maxima of the stationary distribution of the Markov processes are local minima of the Lyapunov quantity for sufficiently large populations, which in turn are generalized evolutionarily stable states. In this manner, we have effectively extended the folk theorem of evolutionary dynamics [11] [32] to these stochastic processes, including those with mutation (in fact we generally require some mutation for the stationary distribution to exist). Similar yet more nuanced results hold for the WrightFisher process on two types. For dynamics on populations of three types, such as the famous rockscissorspaper dynamics, we show that similar results hold for the Moran process, and that the stationary distribution and Lyapunov quantity recapitulate the phase portraits of the replicator equation. Our main results apply to populations with any finite number of types. Finally we demonstrate some extensions to the dynamics of birthdeath processes in structured populations, i.e. in the context of evolutionary graph theory, and of populations of variable size.
The use of informationtheoretic quantities such as cross entropy and relative entropy [37] date back at least to [1] and [7]. Recently author Harper used generalizations from statistical thermodynamics [43] to extend the wellknown Lyapunov stability result for the replicator dynamic to a large class of geometries [24], simultaneously capturing the facts that the relative entropy is a Lyapunov function for the replicator dynamic and the Euclidean distance is a Lyapunov function for the projection dynamic [53]. (Local versions of such results were given first in [30] using Riemannian geometry.) Author Fryer extended the local Lyapunov result for the relative entropy from the replicator dynamic to a class of dynamics known as the incentive dynamics [17], which includes the logit, Fermi, and bestreply dynamics, and introduced the concept of an incentive stable state (ISS) as a generalization of evolutionary stability [69]. The authors results were combined and further extended to various timescales (including the standard discrete and continuous timescales) in [28], yielding a vastly general form of the classical Lyapunov stability theorem of [1]. Now we extend some of these results to finite populations.
We proceed by applying a local/discrete variant of relative entropy to the study of Markov processes and drawing inspiration from methods from statistical inference. Bayesian inference is analogous to the discrete replicator dynamic [58] and relative entropy is a commonly used measure of convergence of an inference process. Simply stated, given a particular population state, we look at the expected next state of the system, formed by weighting the adjacent states by the transition probabilities of moving to those states, and compare this expected state to the current state via the standard relative entropy function (and others) from information theory. Intuitively, a population state is stable if the expected population state is close to current state. Indeed, we show that extrema of the stationary distribution minimize relative entropy between the expected state and the current state for a variety of processes. In particular we focus on the incentive process with mutation, a mapping of Fryer’s incentive dynamic to finite populations, first introduced in [26], a generational version called the fold incentive process, and the WrightFisher process.
We can say more, in fact. For the distance between the expected next state and the current state, we can use any of the informationtheoretic divergences [28], which range from the Euclidean distance for to the relative entropy for . The Euclidean distance will give the best results for finite populations since it is welldefined on the boundary of the simplex as well as the interior, and will detect stable points everywhere. The relative entropy, however, will yield the best connection to deterministic dynamics in the limit that , where we recapture the classical results for the replicator equation.
The authors explored the incentive process without mutation in [27], where we investigated the fixation probabilities of the process in the absence of mutation. Together these two approaches give a complete characterization of the behavior of the incentive process, whether the “equilibria” of the process are the boundary fixation states or if there are interior maxima of the stationary distribution. Stationary distributions of the Moran process and variants with mutation have been studied by Claussen and Traulsen [10] and others e.g. [60]. We note that the incentive processes captures a similar set of processes as Sandholm’s microfoundations approach using revisions protocols [49] [51] [52].
2. Preliminaries
2.1. Incentive Proportionate Selection with Mutation
The incentive process was briefly introduced in the appendix of [26], as a generalization of the Moran process incorporating the concept of an incentive. An incentive is a function that mediates the interaction of the population with the fitness landscape, describing the mechanisms of many dynamics via a functional parameter, including replicator and projection dynamics, and logit and Fermi processes. A Fermi incentive is frequently used to avoid the general issue of dividing by zero when computing fitness proportionate selection [67] [65]. The authors described the fixation probabilities of the incentive process without mutation in [27]. We now describe the incentive process with mutation, which captures a variety of existing processes, such as those used in [10] and [60].
Let a population be composed of types of size with individuals of type so that . We will denote a population state by the tuple and the population distribution by . Define a matrix of mutations where may be a function of the population state for our most general results, but we will typically assume in examples that for some constant value , the mutation matrix takes the form for and . Finally we assume that we have a function which takes the place of in the Moran process. See Table 2 for a variety of example incentives. We will denote the normalized distribution derived from the incentive function as .
The incentive process is a Markov process on the population states defined by the following transition probabilities, corresponding to a birthdeath process where birth is incentiveproportionate with mutation and death is uniformly random. To define the adjacent population states, let be the vector that is 1 at index , 1 at index , and zero otherwise, with the convention that is the zero vector of length . Every adjacent state of state for the incentive process is of the form for some . At a population state we choose an individual of type to reproduce proportionally to its incentive, allowing for mutation of the new individual as given by the mutation probabilities. The distribution of incentive proportionate selection probabilities is given by ; explicitly, the th component is
(1) 
We also randomly choose an individual to be replaced, just as in the Moran process. This yields the transition probabilities
(2) 
We will mainly consider incentives that are based on fitness landscapes of the form for a game matrix ; some authors do not allow selfinteraction and use fitness landscapes of the form (e.g.):
for a game matrix defined by
For our purposes the difference will have little impact. Though we primarily investigate two oneparameter families of incentives defined in terms of a fitness landscape, we note that the incentive need not depend on a fitness landscape or the population state at all. The Fermi process of Traulsen et al is the Fermi for [65], and is called the logit incentive, which is used in e.g. [3]. The replicator incentive has previously been studied in the context of evolutionary game theory [24] [28] and derives from statisticalthermodynamic and informationtheoretic quantities [68]. Recently human population growth in Spain has been shown to follow patterns of exponential growth with scalefactors [29]. Other than the bestreply incentive, which we will consider primarily as a limiting case of the Fermi, we will assume that all incentives are positive definite if to avoid the restatement of trivial hypotheses in the results that follow. For the Fermi incentives, positivity is of course guaranteed for all landscapes. This is particularly convenient for threetype dynamics since the meanfitness can frequently be zero (for zerosum games such as the rockpaperscissors game), which would cause the transition probabilities to be illdefined. We will also typically assume that incentives are continuous in the population distribution, that is, that they are discretized versions of continuous functions on the probability simplex.
At first glance the introduction of incentives might seems like just an alteration of the fitness landscape. While it is true that on the interior of the probability simply many incentives can be rewritten as nonlinear fitness functions, there are significant advantages to the change in perspective from fitnessproportionate selection to incentive proportionatereproduction. In particular, the authors showed in [28] that incentives ultimately lead to a deeper understanding of evolutionary stability and a substantial improvement in the ability to find Lyapunov functions (via formal and geometric considerations) for a widerange of evolutionary dynamics. Moreover, when dynamics interact with the boundary, such as innovative and nonforward invariant dynamics, incentives yield a better description of evolutionary stability. Finally, we can study how a population interacts with the same fitness landscape by varying the incentive function and its parameters rather than the fitness landscape itself.
Dynamics  Incentive 

Projection  
Replicator  
Replicator  
Fermi  
Best Reply 
2.2. The WrightFisher Process
In contrast to the Moran process, which models a population in terms of individual birthdeath events, the WrightFisher process is a generational model of evolution [35] [12]. Each successive generation is formed by sampling, proportionally to incentive, the current generation. Define the WrightFisher Process with mutation for evolutionary games by the following multinomial transition probabilities:
(3) 
where is defined as before in Equation 1. This is a slight generalization of the basic process as given by Imhof and Nowak to include mutation, though we will not consider parameters for intensity of selection via the , linear combination approach as in [35]. In contrast to the Moran process, the WrightFisher process is not tridiagonal, rather every state is accessible from every other state, so long as the incentiveproportionate probabilities are nonzero.
2.3. fold Incentive process
Since the WrightFisher process is a generational process, replacing the entire population in each iteration, and the incentive process is atomic process, they can exhibit very different behaviors. Consider the following process, which will be referred to as the fold incentive process, as an intermediate between the two types of processes. Define each step of the process as steps of the Moran process, so that is the Moran process, and , where is the population size, yields a generational processes derived from the incentive process. The crucial difference is the simultaneity of the replication events, as all occur instantaneously for the WrightFisher process, so the transition probabilities for the fold process differ in general.
The transition probabilities of the fold process can be computed directly from the transition matrix of the incentive process (equation 2) by simply computing the th power of the transition matrix. The entries of the transition matrix, correspond to the probability of moving from population state to population state in exactly steps of the Moran process.
2.4. Stationary Distributions
For a deterministic dynamic there are notions of longterm behavior, such as an limit set[31], and powerful tools such as the Lyapunov stability theorems. For Markov processes we lack any notion of deterministic behavior since the population trajectories can only be determined probabilistically from previous states. We do, however, have a tool that does not exist in the deterministic setting: the stationary distribution of the Markov process. The stationary distribution is a probability distribution on the states of the process, indicating the likelihood of finding the population in any particular state in the longrun. It is characterized in multiple ways, of which we will use that facts that (1) the stationary distribution is the left eigenvector of the transition matrix with eigenvalue one, that is , and (2) the rows of the limit convergence to the stationary distribution for any starting distribution as .
Stationary distributions for the Moran process in two dimensions and some recentlystudied generalizations are given in several recent works [10] [4] [60] [36] [23] [2] [71]; stationary distributions for both the Moran process and the WrightFisher process have been studied in [33]; a number of formulas for various finite population processes are given in [50] and [54]. The solution only relies on the fact that the transition matrix is tridiagonal, and so applies to the incentive process with mutation without modification. The components of the stationary distribution satisfy the detailedbalance condition (also known as the local balance condition) and such processes have a stationary distribution given by the following formula. For any sequence of states ,
(4) 
where can be obtained from the normalization . This particular formulation relies on the fact that the transitions between neighboring states are nonzero, which we assume is valid throughout. For games with more than two types, none of the processes described so far satisfy the detailedbalance condition in general [61], though the neutral landscape is reversible for some choices of mutation matrix. It is, however, always the case that the global balance conditions are satisfied:
(5) 
which follows easily from the fact that the stationary distribution is the right eigenvector of the transition probability matrix . Though we will not be able to give explicit closed forms in the higher dimensional cases (in general they are not known), the global balance equations will be sufficient for our analytic purposes. Exact analytic forms for the stationary distribution of the WrightFisher process are similarly not known. Nevertheless, from a computational perspective, for any concrete values of the various parameters of these processes, the stationary distribution can be computed fairly efficiently for relatively large populations.
3. Notions of Stability
Let us pursue notions of stability for these Markov processes. Intuitively, a population state is stable if there is a reluctance of the population to move from that state. We characterize this in multiples ways. Evolutionary stability has been studied extensively and is closely related to the concept of a Nash equilibrium, which can be loosely described as a state in which no strategy has incentive to deviate from. In [16], Fryer defines an incentive stable state (ISS) as a generalization of evolutionary stability, and shows that there is a generalization of the wellknown Lyapunov theorem for the replicator equation. The authors ported the notion of an incentive stable state to the incentive process in [27]. In some cases an ISS for a particular incentive is again an ESS, such as for the bestreply incentive, and the Fermi incentive. We will give a generalization of those definitions in this work for completeness.
From the Moran process, an analogous form for the standard equation for an interior evolutionary stable state (ESS) for a 2x2 game can be obtained by equating the fitness functions of the two types in the population, matching the fact that for the replicator equation any interior rest point has equal fitness across all types [45]. This can similarly be obtained by equating the transition probabilities for the Moran process, which also reproduces the correct form of Fryer’s ISS in equality for the incentive process, and adds and to set of potential stable points (for the Moran process). Since we do not (yet) have inequality in a neighborhood (as in the definition of a standard ESS) from this derivation we will refer to such states as candidate ISS. Explicitly, we define a candidate ISS as a population state such that . For and the replicator incentive, this definition implies that for all and , which is typical for an ESS. Note that solutions of the ISS candidate equation in a finite population will often not be integral.
Although we lack a Lyapunov theory for Markov processes, it is nevertheless desirable to have a quantity that is locally minimal at stable states and decreasing locally otherwise as an indicator of stability. Such a function yields another method of computation of candidate ISS without having to solve the ISS candidate equation, which might be very difficult. We combine this desire with the idea that a stable state for a stochastic process should remain close to itself by defining the expected distribution resulting from a single iteration of the incentive process
We then consider distance functions such as where is the relative entropy or KullbackLeibler divergence, a commonlyused measure of distance between probability distributions.
For this purpose, we define the following oneparameter distance function
(6) 
This family has the property that if , which is an simple consequence of a more general definition (see [24]), and all have the property that iff (positive definite). We will also make use of the fact that they are all bounded above (on our state spaces) by the squared distance [56]
In what follows, we will write for . This brings us to our first propositions. The first gives two forms for the expected distribution of the incentive, both of which are simple algebraic consequences of the definitions given so far. For twotype populations, the second form is, explicitly,
Proposition 1.
For the incentive process, we have two forms for the expected distribution:
Proposition 2.
For both the incentive process and the WrightFisher process, the following are equivalent:

(ISS Candidate definition)

, that is, is a fixed point of

for all
Proof.
(2) For WrightFisher follows from the fact that the expected value of a multinomial distribution is just . The others follow easily from the previous proposition and Equation 1. ∎
For the fold incentive process the natural notion of expected distribution is the th iterate of as defined for the incentive process. We will not give a closed form for ; the following proposition will suffice for our purposes, and is also an easy consequence of our definitions.
Proposition 3.
Let be defined as for the incentive process. Then for the fold incentive process, if and only if .
Hence for all three processes we have that ISS candidates are equivalent to fixed points of the expected distribution function. This definition captures the usual notion of ISS/ESS as follows. We can write (loosely) that , where is the number of types in the population, so that if e.g. or is otherwise relatively small, we have asymptotically that , i.e. a fixed point of (and so a candidate for incentive stability). For the replicator incentive, this implies that for all , which is the case for an ESS, as above. Our results are true regardless of the form of the mutation matrix, but the stable points need not coincide with the ISS of the incentive for large or imbalanced mutations (which is a feature). We note that Garcia et al have found that the form of the mutations can significantly alter to locations of stable points [22].
The stationary distribution yields our final notion of stability. We define a candidate ISS to be an ISS if the stationary distribution at the candidate is locally maximal, and we call such local maxima stationary stable in general. Intuitively, a stable state of the dynamic should occur at some notion of a maximum of the stationary distribution since it indicates the likelihood of finding the process in any particular state, in the long run, and a tendency to remain in the state. Our main results are that in many cases, for sufficiently large populations, local maxima and minima of the stationary distribution are ISS candidates.
We note that several authors (e.g. [4] [23] [2] [71]) have used the stationary distribution as a solution criterion, looking at average abundance of each population type computed by weighting the population states by the stationary distribution. While the authors of this paper find this to be an interesting approach, we use the stationary distribution to measure the stability of particular population states rather than as a measure on individual types. To see how these differ, consider that for any symmetric process on two types, the average abundance is equal for the two types [4]. For instance, the neutral landscape and the landscape given by , for the replicator incentive would have both have equal abundance for the two types and a stationary distribution with maximum at , but so would the neutral landscape as , in which the stationary distribution would be locally maximal at the two fixation states and [18] and zero at the central point. In all these cases strategy abundance gives an intuitive aspect of “equallylikely to be represented” in the long run, yet taking the average loses information about the manner in which this representation manifests (coexistence of all types or domination by a single type). In all cases the given states are stationary stable.
4. Main Theorem for Twotype Incentive Processes
First we show that when the stationary distribution is locally maximal or minimal, the distance between and is minimal. We require a technical lemma to prove the main theorem that concerns continuity of incentiveproportionate replication probabilities. The crux of the proof is essentially that from the stationary distribution we have that if is local maximum or minimum of the stationary distribution then
Note there are two directions for the path through the states in Equation 4, ending with or , so there is a second solution to the ISS max/min criteria as well (see Table ). To satisfy we need
So we need “equality” of the righthand sides to ensure that these criterion are both satisfied, and we can argue that this is the case for sufficiently large by appealing to continuity. This should not be surprising – in both cases, the population is shifting by one in the same direction, and as gets larger, these points are closer together since for large .
Lemma 1.
Proof.
Since is a continuous function of , so is . Let and . For we simply need to choose large enough so that the set of population states is sufficiently finegrained enough so that . Continuous functions on compact spaces (the simplex) are uniformly continuous, so continuous implies there is such that if then for all and . Since , is sufficient. ∎
If the incentive function is continuous as a function of , and the mutations are constant or otherwise continuous as functions of , then we have that is also continuous, satisfying the hypothesis of the lemma.
Theorem 1.
For a continuous incentive and mutation matrix, there is a sufficiently large so that local maxima and minima of the stationary distribution of the incentive process are local minima of the distance for ; interior extrema are local minima of the relative entropy .
Proof.
If is a maximum or minimum of the stationary distribution then by the formula for the stationary distribution (Equation 4) we must have that the last term in the product passes through 1 so that the stationary distribution switches from increasing to decreasing (or vice versa). By continuity there is where the ratio is exactly equal to one by the intermediate value theorem, and so the state satisfies .
Let . Then by the lemma and the triangle inequality, for sufficiently large (scaling if necessary),
For the distance, we have that, by Proposition 2, (at interior states)
and similarly at boundary states for , . This shows that vanishes for large ; at neighboring states, the difference in transitions is larger, and is locally minimal. ∎
A maximum of the stationary distribution need not occur when is maximal, hence it is not a good definition of stability. The converse of the theorem is not true; the Lyapunov quantity may be minimized but fail to have a local maximum or minimum of the stationary distribution. While we have only proven that the theorem for large for arbitrary continuous incentives, it may be possible to make sharper estimates based on the form of the incentive. As we will see from the examples for dimensions greater than 2, for linear landscapes the theorem seems to hold for almost any with few exceptions.
For the replicator equation, it is known that there is at most one asymptotically stable interior state for linear fitness landscapes. We have a similar result in this context as a corollary.
Corollary 4.
For twotype populations, linear fitness landscapes, the Moran process, if as then there is at most one interior stationary stable state as .
Note that the hypothesis that is necessary, otherwise one can find counterexamples with multiple interior stationary maximum for specific parameters, e.g. , and game matrix given by , produces a process with interior local maxima at . In fact for this game, any produces two interior solutions, , so long as is large enough for the solutions to be more than from the boundaries. (Note that Theorem 1 still holds.)
4.1. fold Incentive Process
These two propositions show that the main theorem for locallybalanced birthdeath processes will apply to the fold incentive process for all .
Proposition 5.
The stationary distribution of a Markov process is the same as the stationary distribution of the process , defined by the th power of the transition matrix.
Proof.
The stationary distribution of a Markov chain can be obtained by the rows of the matrix defined by
and so the stationary distributions of the fold incentive process are the same for all , given a fixed transition matrix . ∎
The next proposition follows easily from the fact that the fold incentive process has the same stationary distribution and the definition of matrix multiplication.
Proposition 6.
If a Markov process satisfies the detailedbalance condition, so does the process defined by the th power of the transition matrix.
Corollary 7.
Theorem 1 holds for the fold incentive process for twotype populations.
Examples
As an illustration of how large may need to be, we consider the replicator incentive and game matrix defined by , . Direct substitution shows that is a solution to the ISS candidate equation for all and . (Note that for , and are also solutions.) For the two stationary extrema equations, we have that . For large , these solutions all converge to . Even for small , these solutions differ only by a small amount (). Similarly, for the landscape given by , , we can give an explicit formula for the stationary distribution for arbitrary (see [26]):
(7) 
Since only the binomial factor depends on , it is clear that is the local maximum of the stationary distribution for any and any . So we see explicitly that the property of being stationary stable is dependent on the population size .
Let the fitness landscape be defined by the replicator incentive and the game matrix , (the Moran landscape with relative fitness ) and let for . As shown in [27], for there is an internal candidate ISS. We consider three cases. For , there is no ISS (the transition probabilities never intersect), and no internal local maximum of the stationary distribution. The latter is true for as well, but there is a state where the transitions intersect, and is locally minimal at this state (), which shows that a candidate ISS need not be stable in the sense of the stationary distribution. An imbalanced mutation rate can force an internal equilibrium for , for instance with and . For , the projection incentive, there is a local maximum of the stationary distribution (at ). It is unique and is minimal.
It is possible to have multiple interior stationary stable states for the replicator incentive with . For example, the following parameters give a process with two local maxima of the stationary distribution: , , , and the neutral fitness defined by . It is also possible to have an internal local stationary maximum as well as a boundary maximum, e.g. for , , , , fitness landscape defined by , , , .
5. Main theorem for Globally Balanced Incentive Processes
In general, incentive processes for more than twotypes need not be locallybalanced, so we cannot use Equation 4 to establish the main theorem. In fact, it appears to rarely be the case (the neutral landscape is an exception in some cases), though computations indicate that in some cases, the detailed balance condition is nearly satisfied, and Equation 4 produces a good approximation of the stationary distribution. Regardless, we now generalize Theorem 1 to globallybalanced incentive processes. The proof strategy is similar to the locallybalanced case, using the global balance equation instead.
Any probability inflowoutflow balanced state, i.e. any state such that , will minimize the distance between the expected next state and the current state. We will call such states probability flow neutral. This could happen in principle at e.g. a saddle point, but since we are mostly concerned with stationary maxima for stability purposes, we will focus on maxima.
Proposition 8.
For a continuous incentive and mutation matrix and sufficiently large , local maxima and minima of the stationary distribution are probability flow neutral. Precisely, for , there is an such that for ,
Proof.
We prove just the maximal case. Suppose we have local maxima of the stationary distribution and let . Since all the transition probabilities are continuous in the population state, so is the stationary distribution, which is therefore uniformly continuous and bounded by . For sufficiently large we have that for states adjacent to . By the local balance equation,
so dividing through by and using the above inequality for ,
∎
We take it as given that stationary extrema occur at probability flow neutral states, knowing that the population states may differ by a small amount. Now we show that stationary extrema minimize the distance functions . While the relative entropy is only welldefined on the interior of the simplex, the same result holds when the state space is restricted to a boundary simplex of lower dimension. This means that the relative entropy can still detect boundary stable points of the same process on the restricted subspace, but the other do so on the full simplex.
Theorem 2.
For a continuous incentive and mutation matrix on types, there is a sufficiently large so that probability flow neutral states of the incentive process are local minima of the distance functions for ; interior extrema are local minima of the relative entropy .
Proof.
Since the proof is similar to the twotype case, we just sketch the proof. Let and suppose that is a probability flow neutral state of the process:
We have that , and we use continuity to choose a sufficiently large so to replace righthand terms with (as in the lemma before Theorem 1). Then we have that
∎
Since the fold incentive process has the same stationary distribution as the case, and since the expected value functions agree at equilibrium via Proposition 1, we have the next theorem as an immediate corollary. This shows in particular that the connection between the stationary distribution and evolutionarily stable states can hold for generational processes.
Theorem 3.
Theorem 2 holds for the fold incentive process.
5.1. Three Type Examples
Here we give several interesting examples using I.M. Bomze’s classification of threetype phase portraits for the replicator dynamic with linear fitness landscape (see [6] and the additions and corrections [8]). We use the Fermi incentive for all examples. A full list corresponding to each of the 48 phase portraits is available online at http://people.mbi.ucla.edu/marcharper/stationary_stable/3x3/. Figure 5 gives stationary distributions and expected distances for the following game matrices, to illustrate Theorem 2.
5.2. RockScissorsPaper Games
The rockscissorspaper game is given by a matrix of the form:
In [3] (online supplement), Andrae et al give steadystate distributions for the rockscissorspaper game in the context of entropy production. Figure 6 contains stationary distributions for the rockscissorspaper game; these are qualitatively similar to those in Figure 4 of [3]. Traulsen et al consider the average drift of a relative entropy equivalent Lyapunov quantity in [65] for the Moran process in a finite population, and find that convergence depends on the value of . The rockscissorpaper game for the Moran process does not yield a detailbalanced Markov process [3], and because of the cyclicnature of the process, the transitions are particularly removed from detailedbalance. This is reflected in the fact that the stationary distribution can take many iterations of the transition matrix to converge. These plots also illustrate that the stability of particular states depends significantly on the mutation rate.
5.3. Four Types
For completeness, we include a higher dimensional example. Consider the incentive process defined by the Fermi incentive, , and the game matrix
Although we cannot easily plot the full stationary distribution, we can plot the twodimensional boundary simplices. Figure 7 gives the expectation distance and the stationary distribution for two of the four faces. Three are the same; the fourth, where , is distinct, and similar to the stationary distribution in Figure 1. From inspection of the game matrix, we would predict that is an ISS for the distinct face. The three similar faces have boundary ISS such as . Note that the stationary distribution and expected distances are computed for the full process on the threedimensional state space; only the twodimensional faces are plotted.
6. The WrightFisher Process
The WrightFisher process behaves differently than the Moran process for many combinations of parameters; see, e.g. [35]. In general, even for twotypes, the WrightFisher process is not locally balanced [61]. Because of the combinatorial complexity of the process, the WrightFisher process is often studied via diffusiontype approximations. In some cases, these differential equations produce a qualitatively similar result to our main theorem for the incentive process. For instance in [9], the authors show that for the replicator incentive, the limiting distribution of the diffusion approximation concentrates on ESS. Combined with the fact that is the same as for the incentive process (Prop 1), one can reasonably expect some subset of the parameters to produce a similar result to Theorem 1 for the WrightFisher process. See Figure 5 for examples.
We give an partial analog of Theorem 1 for the WrightFisher process. Crucially, note that a local maximum of the stationary distribution of the WrightFisher process is in fact a global maximum since every state is connected directly to every other state (assuming that that transition probabilities are never zero), and therefore is necessarily unique.
Theorem 4.
Suppose a given incentive and mutation matrix are continuous. For the WrightFisher process, suppose that the stationary distribution has a global maximum at , is symmetric about the maximum, and is otherwise vanishingly small (if the maximum is not central). Then for sufficiently large , the state is an ISS candidate and a minimum of the distance .
Proof.
We only give the case for in detail. Since the stationary distribution satisfies , we have that
where denotes an adjacent state and by continuity for some we must have that (at the maximum since the distribution is symmetric). A bit of algebra shows that this occurs when
For sufficiently large , solutions to this equation are the same as those of the ISS candidate equation. For , looking at all the immediately adjacent states leads to for large , which implies by summation that . ∎
Computationally, we can compute the stationary distribution of the WrightFisher process for any combination of parameters for smaller easily. Figures 8 and 8 give a variety of computational examples, some of which suggest that a more general version of Theorem 4 holds. Note the sensitivity to the mutation rate . Qualitatively, Figure 1 is similar for the WrightFisher process on the interior. Figures for all 48 of the landscape in Bomze’s classification are also available online at http://people.mbi.ucla.edu/marcharper/stationary_stable/3x3/.
7. Extensions
7.1. Nonconstant Populations
Like the Moran process, the incentive and WrightFisher processes keep the population at a fixed size. Here we give an example that suggests that the methods used in this manuscript could be useful for processes with variably sized populations. We use the process defined in [25], and note that there are other models of variable population size process and dynamics in the literature, such as [39].
We modify the Moran process by inserting a coinflip in each step to determine if a birth event or death event is to take place. The probability of this intermediate step can depend on the population state, and so can define a carrying capacity for the population. For instance, we could take an curve that has probability one of a birth event if the population has size , is for , and is 0 for . Clearly there are many variations on this probability, and the Moran process is a special case, if we separate the birth and death events.
Figure 9 shows the stationary distribution (left plot) of such a process with the replicator incentive for the fitness landscape defined by and and the birth probabilities described in the previous paragraph. Now the ternary plot shows the population states where ; the rightward coordinate is the population size . The right plot is of the expected next state from the current state, where the neighboring states now may have a total population size that is greater or less by one. On the right hand plot, the Lyapunov quantity is minimal when , as we might expect from the fitness landscape and the examples earlier in the paper. These minima occur at the maxima of the stationary distribution as expected, since for , we would expect the stationary max for fixed to be . We leave a precise extension of the main theorem to such processes as an open question.
7.2. Evolutionary Graph Theory
Many authors now study Moranlike processes for populations distributed on graphs [38] [47]. So far we have been working in a fullyconnected population, that is, a population on a complete graph in the sense that every replicator interacts with every other replicator (not to be confused with the graph defined by the transitions of the Markov processes). We outline how a result analogous to Theorem 1 could hold for populations on graphs. A birthdeath process on a graph is again a Markov process, and we consider an incentive process on twotypes. Let the graph be a cycle. We argue for the stationary stability of some states over others. Suppose that the replicators are distributed about the cycle (i.e. a proper coloring), as in Figure 10 (left). Then any reproductive event or will necessarily change the state, since every vertex has adjacent neighbors of the other type. Such a state is inherently unstable, and will have a small presence in the stationary distribution.
Similarly, we can follow this reasoning to conclude that the most stable configurations will have the least possible interactions between the two types. In the case of a cycle, this means that the population is segregated into two connected subpopulations, as depicted in Figure 10 (right). All the remains is to determine the relative sizes of the two subpopulations. Proceeding in the spirit of Theorem 1, we look for states such that the population state’s probability of moving up is the same as the probability of moving down . For very small , only replication events at the interface of the two subpopulations will affect the graph state. Suppose that subpopulations are of distribution . Then we have that the probability of increasing the size of the subpopulation of type (red) is simply since we must choose an individual to reproduce, it has to be one of the two on the interface between the subpopulations, and we have to randomly chose the individual of type (blue) to be replaced. Similarly for the subpopulation of individuals, and these transitions are equal when , i.e. when the subpopulations are in the ISS proportions. We have verified computationally that this is in fact the case for some two player games, such as , and the neutral landscape for . Specifically, there is a particular value of such that these configurations are stationary stable for each tested (when consolidated via cyclical symmetry). For , this value is analytically computable to be , and appears sufficient computationally in other cases. For larger the number of states in the Markov chain is too large for a direct approach.
A similar line of reasoning applies to the random graph, and so we conjecture that analogous results should hold for the random graph and kregular graphs (and likely others).
8. Discussion
8.1. Vanishing Mutation
Throughout, we have explicitly assumed that the mutation rates are (at least) nonzero so as to guarantee the existence of the stationary distribution. In [19] [21] [20], it is shown that as , the stationary distribution concentrates on the boundary states for the Moran process (and some others), and is given in terms of the fixation probabilities of the process. These results complement those of this manuscript, connecting stationary stability with fixation probabilities. In other words, it is trivially the case that for , the fixation states are stationary stable for e.g. the Moran process on linear landscapes, since they are locally maximal for the stationary distribution. That they are ISS/ESS also depends on the incentive – this would be true for the the Moran process (e.g. in the case) but not necessarily for e.g. the projection incentive on the neutral landscape.
8.2. Existence of ISS Candidates
We have left the question of whether stationary stable states or ISS candidates exist as a prerequisite for the results presented thus far. For fixed and other parameters fixed, there may not be a value of such that for an arbitrary incentive, but the converse is true: , where is the number of types, gives for all , so the central point of the simplex is an ISS candidate, and it is stationary stable (the process is equivalent to the neutral landscape, which has a known stationary distribution for all for this mutation matrix [36]). Computationally we have observed for fixed population size, there is often a critical value of that determines if there is an internal stationary maximum, which is not surprising given that the processes often fixate when . Hence there are may be analogs to Theorems 1 2 such that for fixed , there exists a mutation matrix that yields a stationary maximum that is also evolutionarily stable, or some other relationship between and that achieves this.
8.3. Large Population Limit
We have shown above that the relative entropy is locally minimal for evolutionary stable states, and so we can view Theorem 2 as a finite population analog of the fact that is a Lyapunov function for the replicator dynamic when is evolutionarily stable. Explicitly, Traulsen et al show [63] [64] [66] [48] [72] essentially that the Moran process maps to the stochastic replicator equation with mean dynamic (in the case). The stable rest points of this dynamic are those that satisfy the ISS candidate condition when we replace with , and we have that the relative entropy is locally minimal at these points. Moreover, this limiting process shows that any stable rest point of the replicator dynamic must be a limit of stationary stable states for some sufficiently large (the stationary states are by definition rational and the rest point could be irrational).
What happens to the stationary distribution in the large population limit? Consider the explicit stationary distribution in Equation 7. For large , we have that the stationary distribution converges to a binomial distribution with . Scaling to the simplex , we have that the stationary distribution is, by the normal approximation to the binomial, a normal distribution centered at with variance . As , the stationary distribution becomes a delta function centered on the ESS. Furthermore, since the replicator equation is deterministic, the expected next state is also a deterministic quantity, and naturally replaced by the ESS , and so the Lyapunov quantity becomes the function . This is essentially why the stationary distribution and the relative entropy capture the properties of the vector field of the replicator equation, as in Figure 1. We highly encourage the reader to compare Figure 3 to Figure 1.4 of [65] for a striking visualization of the relationship between a Lyapunov quantity (essentially equivalent to the relative entropy) for the replicator dynamic and the functions that we have defined in this work. Sandholm et al also derive a number of large population limits of the stationary distribution for similar processes [50] and [54].
The key difference between the distance functions is the natural geometry of the simplex, the Shahshahani geometry [57], is closely associated with the relative entropy. Forwardinvariance of the replicator equation is a property acquired because of the limit , and not true for the Markov processes considered above. This substantially changes the character of the limits as since for integral population states, the process gets no closer to zero than .
8.4. Conclusion
We have introduced the concept of stationary stability for finite populations and have shown that it captures the traditional notion of evolutionary stability. In the large population limit, our results limit to the known Lyapunov theory of evolutionary stability via Traulsen et al’s mapping of the Moran process to the replicator equation. We have demonstrated that stationary stability works for all mutation rates and for a large class of incentives, connects with the fixation theory of the Moran process, and likely extends to finite populations on graphs.
Crucially, we have given Lyapunov quantities that are very easy to compute as a means of understanding the stationary distribution, which can be very difficult to study analytically, and for which the stationary distribution can not be characterized in closed form (currently). These quantities derive from informationtheoretic distance functions ranging from the Euclidean distance to the relative entropy (or KullbackLeibler divergence).
Methods
All computations were performed with python code available at https://github.com/marcharper/stationary. This package can compute exact stationary distributions for detailbalanced processes and approximate solutions for all other cases mentioned in this manuscript. All plots created with matplotlib [34], except the vector field in the graphical abstract, created with Dynamo [15].
Acknowledgments
This work was partially supported by the Office of Science (BER), U. S. Department of Energy, Cooperative Agreement No. DEFC0202ER63421.
References
 [1] Ethan Akin and Viktor Losert. Evolutionary dynamics of zerosum games. Journal of mathematical biology, 20(3):231–258, 1984.
 [2] Benjamin Allen and Corina E Tarnita. Measures of success in a class of evolutionary models with fixed population size and structure. Journal of mathematical biology, pages 1–35, 2012.
 [3] Benjamin Andrae, Jonas Cremer, Tobias Reichenbach, and Erwin Frey. Entropy production of cyclic population dynamics. Physical review letters, 104(21):218102, 2010.
 [4] Tibor Antal, Martin A Nowak, and Arne Traulsen. Strategy abundance in 2 2 games for arbitrary mutation rates. Journal of theoretical biology, 257(2):340–344, 2009.
 [5] Tibor Antal, Arne Traulsen, Hisashi Ohtsuki, Corina E Tarnita, and Martin A Nowak. Mutationselection equilibrium in games with multiple strategies. Journal of theoretical biology, 258(4):614–622, 2009.
 [6] Immanuel M Bomze. Lotkavolterra equation and replicator dynamics: a twodimensional classification. Biological cybernetics, 48(3):201–211, 1983.
 [7] Immanuel M Bomze. Cross entropy minimization in uninvadable states of complex populations. Journal of Mathematical Biology, 30(1):73–87, 1991.
 [8] Immanuel M Bomze. Lotkavolterra equation and replicator dynamics: new issues in classification. Biological cybernetics, 72(5):447–453, 1995.
 [9] Fabio ACC Chalub and Max O Souza. The frequencydependent wright–fisher model: diffusive and nondiffusive approximations. Journal of mathematical biology, pages 1–45, 2011.
 [10] Jens Christian Claussen and Arne Traulsen. Nongaussian fluctuations arising from finite populations: Exact results for the evolutionary moran process. Physical Review E, 71(2):025101, 2005.
 [11] Ross Cressman. Evolutionary dynamics and extensive form games, volume 5. the MIT Press, 2003.
 [12] Warren J Ewens. Mathematical population genetics: I. Theoretical introduction, volume 27. Springer, 2004.
 [13] S Ficici, J Pollack, et al. Effects of finite populations on evolutionary stable strategies. In Proceedings of the 2000 genetic and evolutionary computation conference, pages 927–934. MorganKaufmann, 2000.
 [14] Gary B Fogel, Peter C Andrews, and David B Fogel. On the instability of evolutionary stable strategies in small populations. Ecological Modelling, 109(3):283–294, 1998.
 [15] Francisco Franchetti and William H Sandholm. An introduction to dynamo: Diagrams for evolutionary game dynamics. Biological Theory, pages 1–12, 2012.
 [16] Dashiell EA Fryer. On the existence of general equilibrium in finite games and general game dynamics. arXiv preprint arXiv:1201.2384, 2012.
 [17] D.E.A. Fryer. The kullbackliebler divergence as a lyapunov function for incentive based game dynamics. arXiv preprint arXiv:1207.0036, 2012.
 [18] Drew Fudenberg, Lorens Imhof, Martin A Nowak, and Christine Taylor. Stochastic evolution as a generalized moran process.
 [19] Drew Fudenberg and Lorens A Imhof. Imitation processes with small mutations. Journal of Economic Theory, 131(1):251–262, 2006.
 [20] Drew Fudenberg and Lorens A Imhof. Monotone imitation dynamics in large populations. Journal of Economic Theory, 140(1):229–245, 2008.
 [21] Drew Fudenberg, Martin A Nowak, Christine Taylor, and Lorens A Imhof. Evolutionary game dynamics in finite populations with strong selection and weak mutation. Theoretical Population Biology, 70(3):352–363, 2006.
 [22] Julián García and Arne Traulsen. The structure of mutations and the evolution of cooperation. PloS one, 7(4):e35287, 2012.
 [23] Chaitanya S Gokhale and Arne Traulsen. Strategy abundance in evolutionary manyplayer games with multiple strategies. Journal of Theoretical Biology, 283(1):180–191, 2011.
 [24] Marc Harper. Escort evolutionary game theory. Physica D: Nonlinear Phenomena, 240(18):1411–1415, 2011.
 [25] Marc Harper. Inferring fitness in finite, variablysized, and dynamicallystructured populations. arXiv preprint arXiv:1303.4566, 2013.
 [26] Marc Harper. Inherent randomness of evolving populations. Physical Review E, 89(3):032709, 2014.
 [27] Marc Harper and Dashiell Fryer. Incentive processes in finite populations. arXiv preprint arXiv:1306.2389, 2013.
 [28] Marc Harper and Dashiell Fryer. Lyapunov functions for timescale dynamics on Riemannian geometries of the simplex. Dynamic Games and Applications, pages 1–16, 2014.
 [29] A Hernando, R Hernando, A Plastino, and AR Plastino. The workings of the maximum entropy principle in collective human behaviour. Journal of The Royal Society Interface, 10(78), 2013.
 [30] Josef Hofbauer and Karl Sigmund. Adaptive dynamics and evolutionary stability. Applied Mathematics Letters, 3(4):75–79, 1990.
 [31] Josef Hofbauer and Karl Sigmund. Evolutionary games and population dynamics. Cambridge University Press, 1998.
 [32] Josef Hofbauer and Karl Sigmund. Evolutionary game dynamics. Bulletin of the American Mathematical Society, 40(4):479–519, 2003.
 [33] Thierry E Huillet. On discretetime multiallelic evolutionary dynamics driven by selection. Journal of Probability and Statistics, 2010, 2010.
 [34] J. D. Hunter. Matplotlib: A 2d graphics environment. Computing In Science & Engineering, 9(3):90–95, 2007.
 [35] L.A. Imhof and M.A. Nowak. Evolutionary game dynamics in a wrightfisher process. Journal of mathematical biology, 52(5):667–681, 2006.
 [36] Kshitij Khare and Hua Zhou. Rates of convergence of some multivariate markov chains with polynomial eigenfunctions. The Annals of Applied Probability, pages 737–777, 2009.
 [37] Solomon Kullback and Richard A Leibler. On information and sufficiency. The Annals of Mathematical Statistics, 22(1):79–86, 1951.
 [38] Erez Lieberman, Christoph Hauert, and Martin A Nowak. Evolutionary dynamics on graphs. Nature, 433(7023):312–316, 2005.
 [39] Anna Melbinger, Jonas Cremer, and Erwin Frey. Evolutionary game theory in growing populations. Physical review letters, 105(17):178101, 2010.
 [40] Patrick Alfred Pierce Moran. Random processes in genetics. In Mathematical Proceedings of the Cambridge Philosophical Society, volume 54, pages 60–71. Cambridge Univ Press, 1958.
 [41] Patrick Alfred Pierce Moran et al. The statistical processes of evolutionary theory. The statistical processes of evolutionary theory., 1962.
 [42] Christina A Muirhead and John Wakeley. Modeling multiallelic selection using a moran model. Genetics, 182(4):1141–1157, 2009.
 [43] Jan Naudts. Generalised thermostatistics. Springer, 2011.
 [44] Andrew E Noble, Alan Hastings, and William F Fagan. Multivariate moran process with lotkavolterra phenomenology. Physical Review Letters, 107(22):228101, 2011.
 [45] Martin A Nowak. Evolutionary dynamics: exploring the equations of life. Belknap Press, 2006.
 [46] Martin A Nowak, Akira Sasaki, Christine Taylor, and Drew Fudenberg. Emergence of cooperation and evolutionary stability in finite populations. Nature, 428(6983):646–650, 2004.
 [47] Hisashi Ohtsuki, Jorge M Pacheco, and Martin A Nowak. Evolutionary graph theory: breaking the symmetry between interaction and replacement. Journal of Theoretical Biology, 246(4):681–694, 2007.
 [48] Jorge M Pacheco, Arne Traulsen, and Martin A Nowak. Active linking in evolutionary games. Journal of theoretical biology, 243(3):437–443, 2006.
 [49] William H Sandholm. Excess payoff dynamics and other wellbehaved evolutionary dynamics. Journal of Economic Theory, 124(2):149–170, 2005.
 [50] William H Sandholm. Simple formulas for stationary distributions and stochastically stable states. Games and Economic Behavior, 59(1):154–162, 2007.
 [51] William H Sandholm. Stochastic evolutionary game dynamics: Foundations, deterministic approximation, and equilibrium selection. American Mathematical Society, 201(1):1–1, 2010.
 [52] William H Sandholm. Stochastic imitative game dynamics with committed agents. Journal of Economic Theory, 147(5):2056–2071, 2012.
 [53] William H Sandholm, Emin Dokumacı, and Ratul Lahkar. The projection dynamic and the replicator dynamic. Games and Economic Behavior, 64(2):666–683, 2008.
 [54] William H Sandholm and Mathias Staudigl. Large deviations and stochastic stability in the small noise double limit, i: Theory. Institute of Mathematical Economics Working Paper, (505), 2014.
 [55] Yuzuru Sato, Eizo Akiyama, and James P Crutchfield. Stability and diversity in collective adaptation. Physica D: Nonlinear Phenomena, 210(1):21–57, 2005.
 [56] Abdolreza Sayyareh. A new upper bound for kullbackleibler divergence. Appl. Math. Sci, 67:3303–3317, 2011.
 [57] Siavash Shahshahani. A new mathematical framework for the study of linkage and selection. AMS Bookstore, 1979.
 [58] Cosma Rohilla Shalizi. Dynamics of bayesian updating with dependent data and misspecified models. Electronic Journal of Statistics, 3:1039–1074, 2009.
 [59] John Maynard Smith. Evolution and the Theory of Games. Springer, 1993.
 [60] Christine Taylor, Drew Fudenberg, Akira Sasaki, and Martin A Nowak. Evolutionary game dynamics in finite populations. Bulletin of mathematical biology, 66(6):1621–1644, 2004.
 [61] Christine Taylor, Yoh Iwasa, and Martin A Nowak. A symmetry of fixation times in evoultionary dynamics. Journal of theoretical biology, 243(2):245–251, 2006.
 [62] Peter D Taylor and Leo B Jonker. Evolutionary stable strategies and game dynamics. Mathematical biosciences, 40(1):145–156, 1978.
 [63] Arne Traulsen, Jens Christian Claussen, and Christoph Hauert. Coevolutionary dynamics: from finite to infinite populations. Physical review letters, 95(23):238701, 2005.
 [64] Arne Traulsen, Jens Christian Claussen, and Christoph Hauert. Coevolutionary dynamics in large, but finite populations. Physical Review E, 74(1):011901, 2006.
 [65] Arne Traulsen and Christoph Hauert. Stochastic evolutionary game dynamics. Reviews of nonlinear dynamics and complexity, 2:25–61, 2009.
 [66] Arne Traulsen, Jorge M Pacheco, and Lorens A Imhof. Stochasticity and evolutionary stability. Physical Review E, 74(2):021905, 2006.
 [67] Arne Traulsen, Jorge M Pacheco, and Martin A Nowak. Pairwise comparison and selection temperature in evolutionary game dynamics. Journal of theoretical biology, 246(3):522–529, 2007.
 [68] Constantino Tsallis. Possible generalization of boltzmanngibbs statistics. Journal of statistical physics, 52(1):479–487, 1988.
 [69] Jörgen W Weibull. Evolutionary game theory. MIT press, 1997.
 [70] Geoff Wild and Peter D Taylor. Fitness and evolutionary stability in game theoretic models of finite populations. Proceedings of the Royal Society of London. Series B: Biological Sciences, 271(1555):2345–2349, 2004.
 [71] Bin Wu, Arne Traulsen, and Chaitanya S Gokhale. Dynamic properties of evolutionary multiplayer games in finite populations. Games, 4(2):182–199, 2013.
 [72] Da Zhou and Hong Qian. Fixation, transient landscape, and diffusion dilemma in stochastic evolutionary game dynamics. Physical Review E, 84(3):031907, 2011.