Entropic Equilibria Selection of Stationary Extrema in Finite Populations

Abstract

We propose the entropy of random Markov trajectories originating and terminating at the same state as a measure of the stability of a state of a Markov process. These entropies can be computed in terms of the entropy rates and stationary distributions of Markov processes. We apply this definition of stability to local maxima and minima of the stationary distribution of the Moran process with mutation and show that variations in population size, mutation rate, and strength of selection all affect the stability of the stationary extrema.

1. Introduction

This work is motivated by the stationary stability theorem [1], which characterizes local maxima and minima of the stationary distribution of the Moran process with mutation in terms of evolutionary stability. Specifically, the theorem says that for sufficiently large populations, the local maxima and minima of the stationary distribution satisfy a selective-mutative equilibria criterion that generalizes the celebrated notion of evolutionary stability [2]. This means that the stationary distribution encodes the usual information about evolutionary stability. Precisely which equilibria are favored (i.e., are maxima or minima) is a natural question and depends on the choice of various parameters, such as the mutation rate $\mu$, the strength of selection $\beta$, and the population size N.

We propose the random trajectory entropy (RTE) of paths originating and terminating at a state as a measure of stability of the state [3,4]. This is an information-theoretic quantity that is easily computable from the entropy rate and stationary distribution of a process, and varies continuously with the critical evolutionary parameters (as does the stationary distribution). We will see that RTE captures the behavior of the Moran process with mutation intuitively, leading to a simple method for equilibrium selection for finite populations—generally a significant problem in evolutionary game theory [5,6].

2. Results

2.1. Stationary Distributions, Entropy Rates, and Random Trajectory Entropies

Our first goal is to establish the random trajectory entropy (RTE) of a state as a measure of stability of the state. We are particularly concerned with the local and global extrema of the stationary distribution, shown in [1] to have a close connection with evolutionary stability.

The stationary distribution of a Markov process gives the probability that the process will be in each state in the long run [7]. As such, it is a fundamental convergence concept for Markov processes. We take the weighted graph viewpoint of Markov processes on a finite set of states V. Let the transition probabilities be given by a function $T:V\times V\to [0,1]$ (viewed as a matrix or a function), and the stationary distribution by a function $s:V\to [0,1]$ (appropriately normalized to a probability distribution). We assume throughout that all processes are irreducible (there is a path between any two states) and have unique stationary distributions. Let $V^{\prime }\subset V$ and define a stationary maximum of $V^{\prime }$ to be a state $v\in V^{\prime }$ such that $s\left(v^{\prime }\right)<s\left(v\right)$ for all $v^{\prime }\in V^{\prime }\setminus v$. Then, we have a local maximum v if the set $V^{\prime }$ is the set of neighboring states of v and a global maximum if $V^{\prime }=V$ (similarly for minima).

Although the stationary distribution of a process is often quite useful, it does not tell the full story of the process. While the stationary distribution gives the long-run occupancy of any particular state, it does not explain how much the process moves among states, and so gives an incomplete description of the dynamic stability of a state. Entropy rate is a generalization of Shannon entropy to Markov processes, and is commonly described as the inherent randomness or information content of a process [3]. The entropy rate of a process encodes both long-term and short-term information about the process, defined for a process X as follows:

$$H\left(X\right)=-\sum _{i,j}s\left(v_i\right)T(v_i,v_j)\log T(v_i,v_j).$$

(1)

The entropy rate is a value attached to a process rather than individual states. To measure the stability of a state, we need a quantity associated to both the process and the individual states that can discriminate between equilibria. Following [3], define the probability of a trajectory $V:v_0\to v_1\to \cdots \to v_k$ with no intermediate state being $v_k$ as the product of the transitions along the path

$$Pr\left(V\right)=T(v_0,v_1)T(v_1,v_2)\cdots T(v_{k-1},v_k).$$

(2)

Since the process is irreducible, we have that the sum over all possible such trajectories from $v_0$ to $v_k$ is one, so they form a probability distribution. Let $\mathcal{T}(v_0,v_k)$ be the set of all such paths and define the random trajectory entropy (RTE) from $v_0$ to $v_k$ to be the entropy of the probability distribution on $\mathcal{T}(v_0,v_k)$; that is,

$$H_{v_0{v}_k}=-\sum _{v\in \mathcal{T}(v_0,v_k)}Pr\left(v\right)\log Pr\left(v\right).$$

(3)

It was shown in [3] (Theorem 1, p. 1419) that when the starting and ending states are the same, the RTE is determined by the entropy rate and the stationary probability:

$$H_v:=H_{vv}=\frac{H\left(X\right)}{s\left(v\right)}.$$

(4)

From this we immediately have the following theorem characterizing local and global extrema of the stationary distribution.

Theorem 1.

For an irreducible Markov process with stationary distribution s, a state s is a local (resp. global) maximum (resp. minimum) if and only if the RTE $H_v$ is a local (resp. global) minimum (resp. maximum).

Furthermore, we now recognize the random trajectory entropy as a measure of the stationary instability of a state, which we can now use to compare and select equilibria for the same process and for closely related processes. Intuitively, a smaller RTE means that trajectories tend to stay near a local maxima (i.e., that random walks tend to be short), which is a way of saying that the state is stable (note that $1/s\left(v\right)$ is the expected number of steps it takes to return to v).

2.2. Applications

We now consider several explicit examples of finite population processes.

2.2.1. Moran Process with Mutation

For the Moran process with mutation, we use a special case of the formulation [1]; see also [8,9,10]. Let a population be composed of n types $A_1,\dots ,A_n$ of size N with $a_i$ individuals of type $A_i$ so that $N=a_1+\cdots +a_n$. We will denote a population state by the tuple $a=(a_1,\dots ,a_n)$ and the population distribution by $\overline{a}=a/N$. We assume the existence of a fitness landscape f where $f_i\left(\overline{a}\right)$ gives the fitness of type $A_i$. Typically, $f\left(\overline{a}\right)=G\overline{a}$ for some game matrix G (see [11,12,13] for general references on evolutionary games). Define a matrix of mutations M where $0\le M_{ij}\le 1$ 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 $\mu$, the mutation matrix takes the form $M_{ij}=\mu /(n-1)$ for $i\ne j$ and $M_{ii}=1-\mu$. A typical mutation rate is $\mu \approx 1/N$ [1,14].

The Moran process with mutation is a Markov process on the population states defined by the following transition probabilities, corresponding to a birth–death process where birth is fitness-proportionate with mutation and death is uniformly random. To define the adjacent population states, let $i_{u,w}$ be the vector that is 1 at index u, −1 at index w, and zero otherwise, with the convention that $i_{u,u}$ is the zero vector of length n. Every adjacent state of state a for the Moran process is of the form $a+i_{u,w}$ for some $1\le u,w\le n$. At a population state a we choose an individual of type $A_i$ to reproduce proportionally to its fitness, allowing for mutation of the new individual as given by the mutation probabilities. The distribution of fitness proportionate selection probabilities is given by $p\left(\overline{a}\right)=M\left(\overline{a}\right)\overline{\phi }\left(\overline{a}\right)$. Explicitly, the i-th component is

$$p_i\left(\overline{a}\right)=\frac{{\sum }_{k=1}^n{\phi }_k\left(\overline{a}\right)M_{ki}}{{\sum }_{k=1}^n{\phi }_k\left(\overline{a}\right)},$$

(5)

where the function $\phi \left(\overline{a}\right)={\overline{a}}_i{f}_i\left(\overline{a}\right)$. We also randomly choose an individual to be replaced, just as in the Moran process. This yields the transition probabilities

$$\begin{array}{cc}\hfill T_a^{a+i_{u,w}}& =p_u\left(\overline{a}\right){\overline{a}}_w\mathrm{for}u\ne w,\hfill \\ \hfill T_a^a& =1-\sum _{b\mathrm{adj}a,b\ne a}{T}_a^b.\hfill \end{array}$$

(6)

We will also utilize a variant incorporating a strength of selection term $\beta$ called Fermi selection [15]:

$$\phi \left(\overline{a}\right)={\overline{a}}_i{e}^{\beta f_i\left(\overline{a}\right)}.$$

(7)

For our examples, we will restrict our attention to processes defined by $X=X(N,n,\mu ,\phi )$. Several explicit examples of stationary distributions for Moran processes with mutation are given in [9,16]. The entropy rate of the Moran process with mutation was computed in [16] (for $n=2$) and [17] (for $n>2$) along with the development of a number of theoretical results. For our purposes, analytic values of the entropy rate are not needed. Generally, as $\mu \to 0$, the entropy rate also goes to zero, and attains its maximum as $N\to \infty$ for the neutral fitness landscape (e.g., with mutations $\mu =1/N$). The entropy rate of is bounded below by zero and above by $\frac{2n-1}{n}\log n$ [17]. The RTE is bounded below by the entropy rate, justifying the description of the entropy rate as the inherent randomness of a process.

2.2.2. Comparison of Equilibria of a Single Process

Since the entropy rate is associated to the entire process, for two different states i and j, we have that $H_i=H\left(X\right)/s\left(i\right)$ and $H_j=H\left(X\right)/s\left(j\right)$, so if $H\left(X\right)\ne 0$ we need only consider the values of the stationary process in this case to compare the equilibria as $H_j/H_i=s\left(i\right)/s\left(j\right)$.

2.2.3. Small Mutation Limit

Consider the special case in which the rate of mutation parameter $\mu \to 0$ in a population of two types. Then, we have ${\lim }_{\mu \to 0}H\left(X\right)=0$ [17]. In this case, the stationary distribution becomes a delta distribution on the corner states. For two population types A and B, we can express the limiting stationary distribution in terms of the fixation probabilities of the two types ${\rho }_A$ and ${\rho }_B$ [8]:

$$\underset{\mu \to 0}{\lim }s(0,N)=\frac{{\rho }_B}{{\rho }_A+{\rho }_B}\mathrm{and}\underset{\mu \to 0}{\lim }s(N,0)=\frac{{\rho }_A}{{\rho }_A+{\rho }_B}.$$

(8)

Hence, we have that

$$\underset{\mu \to 0}{\lim }\frac{H_{(0,N)}}{H_{(N,0)}}=\underset{\mu \to 0}{\lim }\frac{s(N,0)}{s(0,N)}=\frac{{\rho }_B}{{\rho }_A}.$$

(9)

In other words, the state with the type having greater fixation probability is more stable. For the classical Moran process with game matrix $G=\left(\begin{array}{cc}r& r\\ 1& 1\end{array}\right),$ we have that (assuming $r\ne 1$): ${\rho }_A=(1-r^{-1})/(1-r^{-N})$ and ${\rho }_B=(1-r^{1-N})/(1-r^{-N})$, which gives

$$\underset{\mu \to 0}{\lim }\frac{H_{(0,N)}}{H_{(N,0)}}=\frac{{\rho }_B}{{\rho }_A}=\frac{1-r^{1-N}}{1-r^{-1}}.$$

(10)

As expected, whether $r>1$ determines which equilibrium is favored. If $r=1$, ${\rho }_A=1/N={\rho }_B$ and the RTEs are equal.

2.2.4. Large Populations and Neutral Landscapes

For arbitrarily many types, the stationary distribution for the neutral fitness landscape (matrix of all ones) and any mutation rate $\mu$ can be analytically computed. For large N, the neutral landscape attains the maximum entropy rate, so for large populations a sufficient condition for a state for a non-neutral landscape to be more stable than the same state for the neutral landscape is simply to have a larger stationary probability [17]. For non-neutral landscapes, the large population limit need not maximize the entropy rate [17].

2.2.5. Comparison of Equilibria for Separate Processes on the Same States

Two instances of the Moran process with mutation can have the same stationary maximum state but different entropy rates. Consider the one-parameter family corresponding to a Hawk–Dove matrix

$$G=\left(\begin{array}{cc}1& 2\\ 2& 1\end{array}\right)$$

(11)

and transition probabilities defined by Fermi selection (the parameter is the strength of selection $\beta$). For convenience, fix a population size $N\ge 10$ and N even, and let the rate of mutation be $\mu =1/N$. Then, we have that $(N/2,N/2)$ is the unique stationary maximum [1]. As $\beta$ increases, the stationary probability at the maximum increases more quickly than the entropy rate (which is not monotonic in this case). The net result is that the random trajectory entropy is decreasing as a function of $\beta$ and that the stationary maximum is “more stable” (as is expected for greater strengths of selection; see Figure 1). For other two-player games, the situation is analogous. For example, for the coordination game, the interior equilibrium RTE is decreasing as a function of $\beta$. For both games, we have the intuitive result that the stability (measured by the RTE) of the extrema varies monotonically with the strength of selection.

We now consider multiple examples for the landscape derived from the three-type game matrix:

$$G=\left(\begin{array}{ccc}0& 1& 1\\ 1& 0& 1\\ 1& 1& 0\end{array}\right).$$

(12)

This landscape typically has several local extrema. Let the population size $N^{\prime }=6N$. Then, we have extrema at the simplex corners $(6N,0,0),(0,6N,0),(0,0,6N)$, center of the boundary simplices $(3N,3N,0),(3N,0,3N),(0,3N,3N)$, and the center $(2N,2N,2N)$. Varying either $\beta$ (Figure 2) or $\mu$ (Figure 3) can change which equilibria have the smallest RTE. As $\mu$ increases, stationary probability moves from the corner points of the simplex to the midpoints of the boundary simplices, and also toward the center (similarly for the strength of selection $\beta$). In both cases, the rate of change of the stationary extrema dominates because the entropy rate varies slowly.

Though we have focused on equilibria, the stationary distributions of finite population games can exhibit a variety of complex dynamical behaviors, as depicted in Figure 4. Consider the rock-paper-scissors (RPS) landscape given by the matrix

$$G=\left(\begin{array}{ccc}0& 1& -1\\ -1& 0& 1\\ 1& -1& 0\end{array}\right).$$

(13)

For some parameter choices, RPS landscapes produce an interesting stationary distribution with discretized cycles of constant trajectory entropy, analogous to the concentric cycles for the replicator equation and the fact that the relative entropy is a constant of motion of the replicator equation [12]. Assuming symmetry of the cycle (a large population of size divisible by 3 seems to suffice to yield approximate symmetry), no value on any cycle is a local maximum and the values on the maximal cycle are all global maxima. Note that in this case the stationary stability theorem as stated in [1] does not apply to the cycles (only to local extrema).

2.2.6. Comparison of Equilibria for Process with Varying Population Size

For the final example, we consider the effect of altering the population size N. In this case, the underlying state spaces are different even though the equilibria are generally same for large enough N. For the same number of types n, the entropy rate has the same upper bound (though the entropy rate typically increases with N), and so in order to enable a fair comparison we normalize by the number of states (since the stationary distribution is spread out over a variable number of states). In general, the number of states is $\left(\genfrac{}{}{0pt}{}{N+n-1}{n}\right)$. As for both $\beta$ and $\mu$, varying the population size N changes the favored equilibrium (see Figure 5). However, we note that the RTEs are increasing in N, and so the issue of normalization is critical to the comparison of equilibria for processes with different population sizes.

3. Discussion

We have proposed random trajectory entropy as a measure of the stability of states of finite Markov processes, and considered several examples from finite population biology. Variations of fundamental evolutionary parameters alter the stability of equilibria, agreeing with intuitive expectations. In particular, stability is closely tracked by stationary probability in several example population dynamics. We did not consider RTEs for paths that originate and terminate at different states, but it is reasonable to expect, for example, that a local stationary maxima will have smaller RTE in some neighborhood (and similarly for local minima).

Furthermore, as described by Blume [18], high selection intensity describes a population that is more likely to select best responses. As such, the examples demonstrated in this paper show a bias to certain equilibria based on the application. In population biology, common assumptions are low mutation and low selection intensity, which may favor single type populations. For human learning and decision making, higher selection intensity is expected. Our examples show that mixed (interior) equilibria are more likely in this scenario.

Finally, we note that since the entropy rate is an invariant for Markov processes, both the entropy rate and our measure of the stability of a state of a Markov process allow comparisons between processes on the same state space across parameter space and when the state space varies (due to changing population size). The various examples in this manuscript show that this is a useful way to capture the impact of varying parameters in isolation and when the parameters are related, such as when $\mu =1/N$ and $\mu \to 0$. Altogether, we have shown that random trajectory entropies are simple yet powerful invariants for understanding evolutionary processes and the relationships between mutation, fitness, and drift, and that each of the corresponding parameters can determine which equilibrium is most stable.

4. Methods and Materials

All computations were performed with open source code, archived at [19] in a Python package called stationary. This package can compute exact stationary distributions and entropy rates for reversible processes and approximate solutions for all other cases mentioned in this manuscript. All plots were created with matplotlib [20] and python-ternary [21].

The computation of the RTE for every state of a given process back to itself is as follows:

• Given values for all relevant parameters (e.g., a fitness landscape, $\beta$, N, $\mu$), generate the matrix of transition probabilities;
• Compute the stationary distribution of the process, which can be done efficiently (in memory and time) with a power iteration and sparse matrix implementation, or simply standard matrix multiplication for processes with small state spaces;
• With the stationary distribution and the matrix of transition probabilities, compute the entropy rate using Equation (1);
• For each state, compute the RTE using Equation (4) using the entropy rate and the stationary distribution.