Inferring Fitness in Finite Populations with Moran-like dynamics

03/19/2013
by   Marc Harper, et al.
0

Biological fitness is not an observable quantity and must be inferred from population dynamics. Bayesian inference applied to the Moran process and variants yields a robust inference method that can infer fitness in populations evolving via a Moran dynamic and generalizations. Information about fitness is derived solely from birth-events in birth-death and death-birth processes in which selection acts proportionally to fitness, which allows the method to be applied to populations on a network where the network itself may be changing in time. Populations may also be allowed to change size while still allowing estimates for fitness to be inferred.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 12

page 13

page 16

page 18

01/03/2001

Adaptive evolution on neutral networks

We study the evolution of large but finite asexual populations evolving ...
11/01/2015

Limiting fitness distributions in evolutionary dynamics

Darwinian evolution can be modeled in general terms as a flow in the spa...
01/28/2019

Galton-Watson process and bayesian inference: A turnkey method for the viability study of small populations

1 Sharp prediction of extinction times is needed in biodiversity monitor...
03/22/2021

The dynamical regime and its importance for evolvability, task performance and generalization

It has long been hypothesized that operating close to the critical state...
04/17/2020

"Perchance to dream?": Assessing effect of dispersal strategies on the fitness of expanding populations

Unraveling patterns of animals' movements is important for understanding...
07/12/2021

The Fundamental Theorem of Natural Selection

Suppose we have n different types of self-replicating entity, with the p...
10/01/2013

EVOC: A Computer Model of the Evolution of Culture

EVOC is a computer model of the EVOlution of Culture. It consists of neu...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

The theory of natural selection dominates biological explanation and thinking. Despite the elegance of Darwin’s approach to descent with modification and later refinements and abstractions, accurate models and measurements of real populations of replicating organisms can be very difficult to obtain. Analogies of fitness and physical potentials notwithstanding, fitness is more like a hidden variable than an observable quantity. There is no instrument that measures fitness in a manner analogous to a voltmeter. Determining a hidden quantity can be achieved with statistical inference. In this paper we discuss such a method in finite populations with geometric structure, and address critical theoretical questions such as: How much information does a single birth or death event yield about the fitness of a replicating organism?

Statistical inference has been successfully applied in many areas of biology and bioinformatics, e.g. in phylogenetics [3]. The approach in this paper uses models that differ from some other recent approaches to inferring fitness, such as [21] and [18], yet the goal is much the same: practically quantify evolutionary parameters such as fitness that are very useful intuitively. In particular, we will be concerned with traditional models from population biology, like the Moran process[14] [14] and recent extensions that include graphical distribution [12]

. This is to bring together the knowledge gained from evolutionary game theory, population dynamics, and statistical inference rather than a statement on which models are more appropriate. Throughout this investigation, it will become clear that separating birth and death events yields greater information gain per population state transition, so we will consider some new variants of the Moran model. We study the effectiveness of Bayesian inference applied to these models by simulating large numbers of population trajectories.

Let us first consider a basic version of the general problem of inferring fitness. Consider a population of two types of replicating entities and , with individuals of type and individuals of type , with a total finite population size of . Type can refer to genotype or phenotype as is appropriate. Suppose that replicators of type have fitness and replicators of type have , where is the relative fitness between the two types. Given that the population is in the state where both and (so that both types of replicators are represented in the population), let an individual be selected for replication proportionally due to fitness. More precisely, an individual of type

is selected with probability

and an individual of type is selected with probability . The state of the population is not generally informative about the relative fitness because it is not possible to know without the history of the population what the relative differences imply about . Concretely, if , if could be that type is dominating the population and is declining, or that type has a relatively new phenotype that is currently in the process of fixating and has yet to take over the majority of the population. If it were known that the state was stable under some suitable definition of stability then could be estimated by (assuming , which occurs when ). This, however, assumes many things that are likely to be unrealistic in a variety of scenarios – for example (1) that the state of the population can be known precisely (in the laboratory perhaps, but less likely for natural populations), and (2) that the population has converged to a stable state (which implies that the population has been under observation for some time, or was in a stable configuration upon observation). Ideally we would prefer a method that allows one to come upon a population, observe it for a short time, capturing possibly only a random sample of selection events and population states, and make an accurate estimate of the relative fitness . This is achievable.

One way to approach the problem of inferring would be to use the fixation probabilities of type for the population in configuration , which is given by [15]. Then by preparing many populations and recording the proportion of populations that fixation with all members of type or type , the ratio of such events could be used to solve numerically for , by inferring as a binomial parameter. The approach has several flaws: (1) it is likely costly and very time consuming, if even practically possible; (2) it is not applicable to populations on graphical structures for which fixation formulas are not known and may be practically impossible to calculate; and (3) it ignores all the information revealed by the population trajectory – that is, all the information obtainable from replication events – using only the information obtainable from the final state.

What does an observation of a replication of either type or tell us about the relative fitness ? Let us suppose first that . A replication of an individual type then suggests that . How much information such an event gives about about depends on the relative sizes of the quantities and . If is much larger than , a replication of an individual of type is evidence that is larger than one (since based on and alone, a replication of a type replicator should have been more likely), whereas if then the selection landscape was approximately uniform, and a replication of type was just as likely. In this case, the replication event yields little information about the relative fitness. Of course, we cannot put too much stock into a single event (in case of a fluke). Repeated observations strengthen the information yield.

Now let us make the claim that replication events yield information more precise. Assume that we can keep the population in the state , where , and observe arbitrarily many replication events. Let and count the number of observed replications of type and respectively. We can then estimate the value of with a maximum likelihood estimate , and using the values of and determine . This is simply a version of Bayes’ original method of inference for determining the value of a binomial parameter

using Beta distributions with parameters

and , with , and distribution function:

where is the generalized factorial function with for integers . The estimate is the mean of the Beta distribution. To quantify how much information is gained about from each individual replication event, fix and and observe one birth event. This yields and where either or . Then the information gain associated to the birth event is given by the Kullback-Leibler information divergence [11]. To be clear, this is the information gained experimentally by the observer about , not the information gained by the population about its environment. Information about the relative fitness is inferred from the population trajectory (in this case, the sequence of states of the population had we not artificially fixed the state).

Unlike in this example, a population will vary among a sequence of population states as birth and death events occur, possibly changing in total size and changing its spatial distribution. At each state, the probabilities of each type replicating will differ since they depend on the population state . We no longer have a single (which depends on

) to infer and it will be more direct to work with probability distributions for

on or for some maximum value . It is also no longer the case that a sample of identically distributed observations is possible – observations of birth and death events will depend on the state and change the state as a population evolves (altering the fitness landscape). Nevertheless, the conclusions of some typical results for similar inference problems hold in simulations despite the lack of these common assumptions.

The Moran process is a birth-death process that models natural selection in finite populations [13] [14]. Study of the Moran process has been extensive, including fixation probabilities for various landscapes and starting states [6] [23] [1] [15], evolutionary stability [20] [5] [4], the evolution of cooperation [16], and in the context of multi-level selection [24]. Unlike deterministic models of selection such as the replicator equation [7]

, the Moran process is a stochastic process, specifically a Markov chain. A Markov process does not have a fixed trajectory for a given starting state, so when we refer to a trajectory we mean one particular sequence of states of a population modeled by the process. Markov chains are typically analyzed in terms of transition probabilities and quantities derivable from the transitions, such as absorption probabilities and mean convergence times. Such analysis for the Moran process can be found in many sources, such as

[15]. For our purposes, we will need actual population sequences so as to infer values of from population trajectories. We will consider the Moran process and variations that separate birth and death events, including generalizations of the Moran process to populations distributed on graphs, which are well known in evolutionary graph theory [12] [19]. In particular, we will consider birth-death processes on graphs in which a member of the population is selection for reproduction from the entire population, replacing an outbound neighbor, and death-birth processes in which a member is randomly selected to be replaced by an inbound neighbor (which is chosen proportionally to fitness).

2. Preliminaries

2.1. Markov Processes for Birth and Death

The Moran process models selection in a well-mixed finite population of replicating entities. As before, consider a population of replicators, of type A (A-individuals) and of type B (B-individuals) where is fixed as before. Individuals of type A and B have fitness and respectively, and may depend on the population parameters and . Although we could dispense with one of the parameters or since is fixed, we will need both for a modification that does not maintain a fixed value for . We will only consider populations of integral size at least 3.

The population is updated by selecting an individual at random proportionally to fitness and selecting an individual at random uniformly to be replaced. For the population to change state, individuals of different types must be chosen. The transition probabilities between states are given by [23] [16]

with . The fitness landscape is given by

for a game matrix defined by

The Moran process is given by . The process has two absorbing states and , corresponding to the fixation of one of the two types. We will first consider the Moran process for populations with and where is the relative fitness of type B versus type A. Given , much can be said about the Moran process in this case, including fixation probabilities starting from any population state. We consider a different problem. Given a trajectory of the Markov process as a series of states can we accurately determine the value of ? We could directly measure the rates of reproduction and compare those values, which we will refer to as counting. From the Moran process we can only detect birth of either type if the population changes state; otherwise we could only guess based on the population distribution which type both had a birth and a death. The combined birth-death transition probabilities wash out the information obtained in the case that the population state stays the same (e.g. individuals of the same type replicates and then dies). For this reason we also consider a modification of the process that breaks the transitions into separate birth and death events. The results of the process are the same as the Moran process but we will know in each step which type replicated and which type died. In this case, the probabilities for fitness proportionate reproduction are given by:

We can also view this process as conditional on the size of the population. If , we select a replicator to reproduce. If the population is of size , then we randomly select a replicator to be removed. Inference benefits from this modification because the transitions which leave the state unchanged now yield information. Since the death event is random, it carries no information about , only about the population size, which is assumed known in this case. We could alternatively consider a death-birth process. Results in this case are similar but slightly more susceptible to stochastic noise for very small population states. Because the death event occurs first, the reproductive event is at a state of smaller population size , but otherwise basically the same as the birth-death case. For populations on graphs, birth-death and death-birth processes can be quite different.

2.2. Bayesian Inference

To use Bayesian inference to infer the value of

from a sequence of states we first choose a prior probability distribution for the possible values of

. In theory , but for computational purposes we will restrict to an interval , where

is chosen to be sufficiently large. Choosing a prior probability distribution is a subjective process, though in this case there are a few reasonable choices. We could use a uniform distribution on

in attempt to chose an uninformative distribution. The uniform prior, however, puts too much weight on the interval which should have approximately the same weight as on since is a ratio. By the neutral theory of evolution [10], most mutations will have little effect on the relative fitness , so priors which are somehow centered at are reasonable choices in many scenarios. It is also reasonable to assume that since both types and

exist in the population (and presumably can reproduce). A gamma distribution with well-chosen parameters fits both assumptions and can have mean, median, or mode at

. We will not dwell on the problem of prior choice other than to say that a different prior could further improve the accuracy for some combinations of parameters and states. In simulations, gamma distributions based on the assumptions of the neutral theory perform well. Later in the text we will see that there is another choice for the prior distribution that is computationally convenient. The gamma distribution takes the form

where is unrelated to the previous usage. Note that there are many distributions that could take the place of the gamma distribution and there is no particular reason for it over other distributions that also have the desired properties, though it is noteworthy that the tail of the distribution drops off exponentially. Simulations reported in this paper use parameters and .

Bayesian inference takes a distribution on the possible values of the parameter and an observation to produce an updated distribution that takes into account the information gained by the observation. Starting with the prior distribution, inference proceeds until the data is exhausted or the parameter is known to a sufficient accuracy. To update the distribution using the sequence of states we use Bayes’ Theorem, where

indicates a transition from :

and we have that from the definition of the Markov process. can be calculated by integrating over , but it is only necessary to normalize at the final step unless we wish to have intermediate estimates. To form an estimate over the entire sequence of states, pair the states into transitions (each state is in two adjacent pairs), and form the posterior distribution . Finally, to extract an estimate for , we normalize the posterior and compute the mean value or mode numerically.

2.3. Conjugate distributions

Just as in the introductory example, the observations of birth events (or population state changes) can be used to determine probability distributions for the purposes of inference. In the case of a binomial distribution, each observation of a success or failure multiplies the inferred distribution by

or (increasing the value of or ). In the case of a birth-death process as described above, we multiply by the corresponding transition factor. Now and

are vector parameters indexed by

(with ) and the distribution takes the form

Rewriting, and introducing a normalization constant yields the form

The normalization constant can be computed symbolically using partial fractions but the formulas are, even for simple values of and , very complex for even relatively small values of ; similarly so for the integration to determine the mean of the distribution. Moreover, analytic formulas for the maximum likelihood estimator and the Fisher information are not easily obtainable. For this reason we proceed numerically, noting that using this conjugate distribution still prevents integration at every step in the inference process. We will refer to this distribution as the FPS distribution (fitness proportionate selection distribution).

The FPS distribution provides more options for the prior distribution . If any then . Similarly, if any then as . For the distribution to have a finite integral in the case that the distribution is supported on , we need that . So, one possible choice of parameters for a prior is for all . For these parameters, the FPS distribution has mode equal to 1 (the maximum likelihood estimate) and mean approximately . The mode is equal to one if for all , and more generally if for all . To see this, notice that the equation for derivative of the logarithm of to be zero can be rewritten as

For close to one, we have the approximation . Merely counting the replication events for each type would lead to an estimate of the form , so it is clear how the position in the population affects the estimate (as well as the true value of ). For this reason we will modify the counting procedure in a later section to weight the value of a replication event by the population state, which significantly improves the results. In general, the maximum likelihood method is a given by the real positive root of a polynomial and is unique in , which can be shown with basic techniques from calculus. Given a prior, the estimate for

given by the mode of the posterior is known as the maximum a posteriori probability estimate (equal to the maximum likelihood estimate if the prior is uniform). This estimate can be easily computed numerically. Note that the full distribution can be used to give a credible interval in addition to a point estimate.

For the Moran process, the analogous distribution includes a third vector parameter corresponding to transitions in which the population state is unchanged:

Similarly to the previous case, we can rewrite this as

Figure 1. Examples of Priors. Blue: Gamma with . Green: FPS distribution with N=30, for all . Red: FPS with and all others 0.

2.4. Conjugate Distributions for Variably-Sized Populations

We will also consider the case that the population size may change. To do so, we merely need to include terms for different values of , where and are now triangular matrix parameters and . The FPS distribution is given by

where is the possible set of total population sizes for the particular process. Call this distribution the variable-population FPS distribution. The previously described FPS distribution is the special case when . This distribution would be used, for instance, with a Markov process in which, rather than have a birth and death event in each cycle, has a birth or death event with some probability that potentially depends on the population state (but not the parameter to be inferred). In this case multiple birth or death events could occur in sequence without the other; still only the birth events will be used to infer the value of . In particular, suppose a population has carrying capacity , so that the probability of a death event is (so that the population cannot exceed its carrying capacity), and where for all population states in which

, perhaps given by a sigmoid function with inflection point at

.

From these examples it should be clear that any Markov chain on population states with parameter dependent transition probabilities could be treated in a similar manner. This paper will only cover the distributions described so far, but the reader could consider more complex examples where the fitness functions and transition probabilities depend on more variables (such as mutation probabilities or fitness functions deriving from unknown game matrices).

2.5. Conjugate Distributions for Populations on Graphs

Suppose now that rather than a well-mixed population, we have populations in which the replicators occupy vertices on a directed graph, such as a directed cycle or a graph with an undirected star topology. For the birth-death process on a graph, a replicator is chosen proportionally to fitness and replaces a randomly selected outgoing neighbor. Since the replicator that reproduces is selected from the entire population, the appropriate probability distribution is the fixed-population FPS distribution. If the graph can change size (i.e. lose or gain vertices), then the variable-population FPS distribution is needed.

For death-birth processes, the variable-population FPS distribution is needed for both static graphs and those in which vertices and edges may be added or removed. In this case, the set consists of the sizes of the sets of inbound neighbors for each vertex (since these are the subpopulations that in which replicators are selected for reproduction).

3. Computations and Simulations

Trajectories for the Markov process were computed with mpsim (Markov process simulator), a flexible and parallelized simulator for any Markov process for which the transition graph can be specified and stored in computer memory, created specifically for the purpose of generating large numbers of trajectories. Source code for mpsim is available on Github at https://github.com/marcharper/mpsim

. Additional code containing more simulations and code to process trajectories, compute posterior probability distributions and relative fitness estimates, is also available on Github at

https://github.com/marcharper/fitness_inference. Parameter estimates are ultimately determined by numerical integration.

Figure 2. The End Result of Inference. In this diagram, the population began at . The top plot shows the posterior distribution after 500 transitions (iterations). The third plot shows the probabilities of selecting type and to reproduce, and respectively, as well as the probability of choosing a type individual to be replaced (). The lowest plot shows the mean, median, and mode of the posterior distribution over all 500 iterations.
Figure 3. PDFs for various stages of one example of inference. Blue: Prior (FPS with ); Green: After 15 steps; Red: After 60 steps; Cyan: Final posterior (920 steps). Population size , true value of is 4; the inferred value is .
Figure 4. Histogram of inferred final values for 10,000 simulations with , , starting state for the Moran process. The prior distribution is a Gamma with . A small portion of trajectories estimate .

3.1. Simulation Results, Moran Process

For the purposes of comparison, for each trajectory the total counts of observations of birth events of the two types is recorded to yield estimates of the form

mimicking the maximum likelihood method as discussed in previous sections. These estimates are compared to numerically computed means and modes of posterior distributions.

Estimates for the value of by counting and by inference via simulation reveal several commonalities. 1000 sample runs for ranging from 3 to 50, for ranging from 0.1 to 2 in steps of 0.1, and for each possible starting state were generated and the methods of determining compared. It is not possible to present results from all possible starting states for all combinations of variables, so we will focus on the most interesting initial starting points: a single mutant, and an initial state with fitness proportionate selection probabilities equal (or as nearly so as possible). The full data set is available on request (and can be regenerated with the previously referred to software). See Figures 5 and 6.

(a) Counting
(b) Bayes Mean
(c) Bayes Mode
Figure 5. Top Row: Color is given by the inferred values of r minus the true value for Counting, Bayesian Mean, and Bayesian Mode (left to right) for 1000 simulations at each pair

. Bottom: Standard deviations of inferred values for same methods. Initial point:

where . Although the deviation pattern looks similar, the magnitude is very different from counting to either inference estimate.
(a) Counting
(b) Bayes Mean
(c) Bayes Mode
Figure 6. Moran Process. Top Row: Color is given by the inferred values of r minus the true value for Counting, Bayesian Mean, and Bayesian Mode (left to right) for 1000 simulations at each pair . Bottom: Standard deviations of inferred values for same methods. Initial point: (i.e. a single mutant). Although the deviation pattern looks similar, the magnitude is very different from counting to either inference estimate.

(1) Starting states adjacent to absorbing states: For all values of and , starting states close to either absorbing state have a tendency to produce short trajectories, yielding little information for either method due to a shortage of data, especially if near a state that favors the true value of

. Moreover, for the counting method, estimates can be completely nonsensical since there may be zero births of one of the two types, leading to estimates of zero or infinity. Inference performs much better in this situation, both on such trajectories and generally, being more accurate and having a smaller variance, but depending heavily on the prior distribution.

(2) Small : for small values of , trajectory length can be short since the process can converge quickly. Nevertheless, both methods perform well if the starting state if is not adjacent to an absorbing state and

, with inference performing better generally. The transition probabilities for small populations are more skewed (as functions of

) than for large populations (e.g. compare with ). This means that observations for small populations can have a larger impact on the inferred value than for larger populations, and that more observations are required to lower variance. Figure 6 shows that smaller population sizes have greater average deviation from the true value of and larger variances in inferred value.

(3) close to 1: For values of close to 1, the lower variance of inference produces significantly better estimates. Because trajectories can ultimately converge to either absorbing state with relatively similar probabilities (depending on the starting state), the additional information obtained by inference in each transition produces far better estimates. For , the variance from counting can be ten times greater than the variance for inference in addition to inference producing a more accurate mean value. Note also that for values of , a prior that has more weight away from can improve estimates, and the choice of prior may favor this case.

Although it appears that counting performs similarly to inference in the extreme case of a single mutant, it is critical to note that these values are only computed in cases that counting could give an estimate. That is, in many cases, especially for extreme relative fitnesses and starting states, only one type registers any replication events, which can lead to estimates of infinite relative fitness. This affect is slightly less extreme for processes that separate death and birth in accounting.

3.2. Separating Birth and Death

Counting suffers from the fact that birth and death events are only clear in the case that the population changes state. Since the most likely transition is very often to stay in the current state, even if the information gained from such a transition is small, inference has a big advantage due to the sheer number of such transitions, especially when close to an absorbing state. Trajectories of the Moran process tend to oscillate near the absorbing state because while the higher fit type typically dominates the population, it is also most likely to be randomly selected for death. For instance, suppose and . Then the population will have a tendency to cycle in the states . The reason is that while and individual of type is far more likely to be chosen to reproduce because and nearly all individuals are of type , it is also unlikely that the lone individual of type will be chosen for death: . This can lead to a large number of replication events near fixation which can lead the estimate to drift upward near the end of the trajectory.

To understand the strength of the effect of separating birth and death on the ability to infer fitness, compare the results from the Moran process to the modified process in which birth and death events are recorded separately. The method of counting also benefits from this process, but care must be taken. Naively proceeding as before is unwise because of the often quite small transition probabilities near absorbing states. To manage this over-abundance of births of the more fit type, we count weighted by population size. If an individual reproduces in a population in state , this is counted as ; similarly births are counted as . This is justified by the estimate given previously and produces much improved estimates (since the transitions in which the population stays in the same state now count as replication events and the large numbers of oscillator transitions near the absorbing states would lead to substantial overestimations). For the modified process, both methods perform well, significantly improving versus the Moran process. Inference still maintains a smaller variance, and overall has lower variance and more accurate estimates for this process than the Moran process. In general, counting performs poorly (see Figures 8 and 7).

(a) Inv. Counting
(b) Counting
(c) Bayes Mode
Figure 7. Separated Moran Process. Top Row: Color is given by the inferred values of r minus the true value for Inverse Counting, Counting, and Bayesian Mode (left to right) for 1000 simulations at each pair . Bottom: Standard deviations of inferred values for same methods. Initial point: where . Although the deviation pattern looks similar, the magnitude is very different from counting to either inference estimate.
Figure 8. Comparison of Inverse counting (blue) and Bayesian mean (green). Curves indicate the deviation from the true value of and the standard deviation over 1000 simulations for and population size along the x-axis.

3.3. Sampling

As a practical matter, the full trajectory of may not be available or practically obtainable, so we also compare the two methods by sampling the trajectories for a subset of transitions in each trajectory. For sample sizes of 10 and 20 (with the full trajectory used if shorter than the sample size), the inference method suffers relatively little and mostly in variance in comparison to counting, which performs poorly even in cases where is large and the starting state is near the central state. So as a practical means of estimating fitness, it is not necessary to know the population trajectory to a high degree of precision, nor is the method dependent on the fact that there is dependence between states and observations. In other words, the fact that an observation of a transition from to is dependent on having had observations that got the population to the state in the first place would lead nonzero indices of and to have nonzero neighbors (but this is not the case for a sampled trajectory).

(a) Inv. Counting
(b) Bayes Mean
(c) Bayes Mode
Figure 9. Separated Moran Process, Sample Size 20. Top Row: Color is given by the inferred values of r minus the true value for Inverse Counting, Counting, and Bayesian Mode (left to right) for 1000 simulations at each pair . Bottom: Standard deviations of inferred values for same methods. Initial point: where . Although the deviation pattern looks similar, the magnitude is very different from counting to either inference estimate.

3.4. Variably-sized Populations and Death-Birth updating

For well-mixed populations, death-birth updating differs little from birth-death updating other than to effective reduce the population size by one. It also makes fixation slightly more likely when the population is in a state adjacent to an absorbing state. For variable size populations, there are multiple effects in play. Replication events when the population is smaller can carry more information but smaller populations mean that fixation is more likely.

Fixation probabilities can be difficult to derive analytically for variably-sized population processes. The population can now fixate in many ways (the process may have many more than two absorbing states, depending on how variable the population is allowed to be). A discussion of the many ways in which to assign probabilities to birth and death in each round of the process would be lengthy; simulations indicate performance vary similar to the fixed-population size case and will be omitted.

4. Results: Populations on Graphs

It will not be possible to cover a comprehensive set of graphs, so we will focus on several interesting cases. Note that to make inferences about , it is not necessary for the graph to be connected or for a particular type to be able to fixate in the population. A sample of a non-absorbing trajectory is enough to make an estimate. The only requirements for directed graphs for each vertex to have an outgoing neighbor (if birth-death is the updating process) or for each vertex to have an incoming vertex (if death-birth is the updating process). For an undirected graph we simply require that each vertex have at least one neighbor. Even these requirements can be relaxed if desired.

4.1. Cycles and -regular graphs

Population trajectories on a graph depend on both the number of replicators of each type and the manner in which they are initially distributed. Consider a population on a cycle. One initial case would be for all the replicators of type to be on a semicircle and all the replicators of type to be on the other semicircle. In this case, only the replication events at the boundaries of either semicircle will alter the population from its initial state. Similarly, suppose the replicators are initially distributed as around the cycle. Every initial replication event will change the population state in this case. Whether or not this favors one type over the other depends on the true value of . Simulations indicate that inferences from populations on cycles are better than those in well-mixed populations. One reason for this is simply that more replication events may occur on average (depending on the initial distribution) than in the well-mixed case because replicators may have a tendency to replace their own types due to the cycle structure (e.g. in the semicircle initial state). This produces longer trajectories, which can yield better estimates.

The directed cycle has a degenerate special case. For death-birth updating on a directed cycle, the estimate of will be completely dependent on the prior. This is because each vertex has a single incoming neighbor so no fitness proportionate reproduction occurs during birth events, and so yield no information. In this case, the FPS distribution is improper, with everywhere.

A -regular graph is a graph in which each vertex has exactly edges. A cycle is a 2-regular graph. For connected -regular graphs, death-birth processes have FPS distribution identical to that of populations of size ( if undirected). For birth death processes, the distribution is that of a population of size , where is the number of vertices. As noted earlier, smaller populations can have more average deviation from the true value of the parameter and more variance in estimates, so the processes of birth-death and death-birth can have substantially different behaviors on a -regular graph if and are significantly different. This is in contrast to a complete graph, in which the two processes are nearly identical (effective size vs. ).

4.2. Star Topology

Birth-death processes on star topologies have previously been shown to enhance the strength of selection [12] versus a well-mixed population. Simulations indicate the star topologies yield more stable inferred values of in both mean and standard deviation versus well-mixed populations (complete graphs), again likely due to longer trajectories. Again birth-death processes differ significantly from death-birth processes, since in the former case every birth event is from a population of size , and replaces the central vertex with high probability. In the death-birth case, death events occurring on the non-central states are always replaced by the central vertex (which carries no information), and death events at the central vertex act as if the population size were . Since death events are equiprobable at every vertex for the death-birth case, the type occupying the central vertex will replicate with probability regardless of its fitness and increase its proportion in the population! See Figure 10 for a comparison of the star, cycle, and complete graph.

Figure 10. Left: Fitness vs. Mean fixation time (steps), Population size for complete graph (blue), undirected cycle (green), and star (red). Right: Fitness vs. Inferred fitness. 200 birth-death simulations per . Initial State: and each occupy half of the graph. Longer trajectories produce better estimates with less variance, with the star graph performing best.

4.3. Dynamic Graphs

It is also possible to infer fitness of birth-death and death-birth processes on graphs with dynamic structure, such as those with active linking [19]. In this case, one would again use the variable-population FPS distribution. As vertices and edges are added or removed, the subpopulations in which fitness proportionate selection occurs may change, and fundamentally little is changed from the case in which an well-mixed population of variable size evolves, from the point of view of computing and estimate from the data. If a vertex has no outgoing neighbors, it replaces itself in the case of a birth-death process. For a death-birth process, there must be at least one incoming neighbor.

The preceding examples of the cycle and the star topology indicate that the effectiveness of inferring the fitness is affected the population structure, and in particular can improve the estimates obtained. Consider a random graph in which the probability that any two vertices are connected depends on a fixed probability . Means and standard deviations for random graphs for a range of probabilities are given in Figure 11. Notice that the estimates are better for small values of , and quickly tend toward a similar distribution as increases (as the graph becomes “more complete”). This case is somewhat similar to the case of a -regular graph. The population gets small-population selection events without the small population tendency to fixate quickly.

Figure 11. Means and standard deviations of parameter estimates (Bayesian mean) for 1000 birth-death trajectories on a random graph with on the horizontal axis with . Initially the population is in the state . Outgoing vertices are selected for each iteration (not fixed) with probability . If no outgoing neighbors are randomly selected, the vertex remains occupied in a birth-death event. This can lead to longer trajectories for smaller and hence better estimates.

5. Quantifying Information Gain

Let us now make good on the promise to quantify how much information is gained by a replication event. Consider the case of a well-mixed fixed size population evolving via the FPS transition probabilities, which in this case form a Bernoulli random variable with regard to the choice of type to replicate. At each step of the trajectory we know the population state

with , the estimate of the value of the fitness , and the true value of (since it was used to create the trajectory). Hence at each step in the process, we can compare the Bernoulli random variables formed by the true value of and the estimate using the formula and . The information divergence between these random variables is

This essentially measures the prediction power of the estimate – divergences close to zero indicate there is little information left to gain. See Figures 12 and 13 for examples with estimates from Bayesian inference. In both cases, the divergence is near zero after just 20% of the trajectory data. This partially explains the surprisingly good estimates arising from small samples of trajectories discussed previously.

Figure 12. Information gain for a population of size , , and starting state (33, 17) for a single trajectory.
Figure 13. Information gain for a population of size , , and starting state (54, 46) for a single trajectory.

6. Discussion

We have seen that it is possible to fairly accurately infer the unknown fitness of a replicator in Moran-like processes in finite populations on graphs. In particular, even if the trajectory of the population is reduced to a small sample of population state transitions, estimates for the relative fitness can be accurate. Inference is more accurate and has substantially less variation than counting methods, and is much more efficient than using fixation probabilities. Moreover, the method of inference allows the incorporation of prior information, such as a prior belief in the neutral theory of evolution, that is not possible (or not obviously so) with counting methods. In many cases counting methods fail to give a useful estimate because one of the types yields no replication events. This could be worked around by using pseudocounts, adding a nonzero constant to the numerator and denominator so as to avoid zero/infinite estimates. It is not clear a priori what values would be appropriate choices for , and in any case, a choice of amounts to an attempt to incorporate prior information into the estimate. Finally, note that Bayesian inference typically yields estimates with far lower variance than the counting methods considered in this work.

At the end of a particular trajectory, the posterior distribution takes on the shape of a normal distribution. This fits the conclusion of the Bernstein von-Mises theorem for Markov processes

[2] which says that the posterior should be a normal distribution with mean at the estimate and variance given by the inverse of the Fisher information of the distribution. This has a natural interpretation here. Simply put, from a single trajectory there are many values of the unknown fitness parameter that are likely to produce the trajectory. The posterior distribution indicates the likelihood that the trajectory was generated by any particular value of the fitness (given the prior distribution). Hence the posterior distribution could be used to estimate the probability that the true fitness is greater than one, that it lies in a particular interval, and other similar statistical calculations.

It is possible to infer the fitness of replicators evolving on a structured spacial distribution, such as for evolutionary games on graphs [22]. We saw that the characteristics of the graph alter the effectiveness of inference, with some graphs improving the accuracy and precision of estimates. This may be related to the fact that some graphs can amplify selection [12]. Dynamic graphs, such as the random graph, behave similarly.

Analogously, it is likely possible to infer fitness of replicating entities in similar dynamical systems. Inferring multiple parameters simultaneously would be of use for Moran processes for more than two types (requiring relative fitness parameters for a -type process), birth-death processes including mutation parameters in addition to (or instead of) fitness parameters, fitness landscapes dependent on game matrices or other additional parameters, Moran-like processes with multiple levels of selection [24], reproductive processes with mechanisms other than fitness proportionate reproduction. It should also be possible to completely determine fitness landscapes, for instance by observing all the entries of a game matrix on which a fitness landscape is based. All of these examples are straightforward variations of the method described in this manuscript.

Methods

All computations were performed with python code available at https://github.com/marcharper/mpsim and https://github.com/marcharper/fitness_inference. All plots created with matplotlib [8]. The software used throughout the results sections depends substantially on the python libraries SciPy [9] and Numpy [17].

Acknowledgments

This research was supported by the Office of Science (BER), U. S. Department of Energy, Cooperative Agreement No. DE-FC02-02ER63421.

References

  • [1] Tibor Antal and Istvan Scheuring. Fixation of strategies for an evolutionary game in finite populations. Bulletin of mathematical biology, 68(8):1923–1944, 2006.
  • [2] J Borwanker, G Kallianpur, and BLS Prakasa Rao. The bernstein-von mises theorem for markov processes. The Annals of Mathematical Statistics, pages 1241–1253, 1971.
  • [3] Joseph Felsenstein and Joseph Felenstein. Inferring phylogenies, volume 2. Sinauer Associates Sunderland, 2004.
  • [4] 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. Morgan-Kaufmann, 2000.
  • [5] 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.
  • [6] Drew Fudenberg, Lorens Imhof, Martin A Nowak, and Christine Taylor. Stochastic evolution as a generalized moran process. Unpublished manuscript, 2004.
  • [7] Josef Hofbauer and Karl Sigmund. Evolutionary game dynamics. Bulletin of the American Mathematical Society, 40(4):479, 2003.
  • [8] J. D. Hunter. Matplotlib: A 2d graphics environment. Computing In Science & Engineering, 9(3):90–95, 2007.
  • [9] Eric Jones, Travis Oliphant, Pearu Peterson, et al. SciPy: Open source scientific tools for Python, 2001–.
  • [10] Motoo Kimura. The neutral theory of molecular evolution. Cambridge University Press, 1985.
  • [11] Solomon Kullback and Richard A Leibler. On information and sufficiency. The Annals of Mathematical Statistics, 22(1):79–86, 1951.
  • [12] Erez Lieberman, Christoph Hauert, and Martin A Nowak. Evolutionary dynamics on graphs. Nature, 433(7023):312–316, 2005.
  • [13] 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.
  • [14] Patrick Alfred Pierce Moran et al. The statistical processes of evolutionary theory. The statistical processes of evolutionary theory., 1962.
  • [15] Martin A Nowak. Evolutionary dynamics: exploring the equations of life. Belknap Press, 2006.
  • [16] 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.
  • [17] Travis E Oliphant. Python for scientific computing. Computing in Science & Engineering, 9(3):10–20, 2007.
  • [18] Jakub Otwinowski and Ilya Nemenman. Genotype to phenotype mapping and the fitness landscape of the e. coli lac promoter. PloS one, 8(5):e61570, 2013.
  • [19] Jorge M Pacheco, Arne Traulsen, and Martin A Nowak. Active linking in evolutionary games. Journal of theoretical biology, 243(3):437–443, 2006.
  • [20] Mark E Schaffer. Evolutionarily stable strategies for a finite population and a variable contest size. Journal of theoretical biology, 132(4):469–478, 1988.
  • [21] Ruth G Shaw and Charles J Geyer. Inferring fitness landscapes. Evolution, 64(9):2510–2520, 2010.
  • [22] György Szabó and Gábor Fáth. Evolutionary games on graphs. Physics Reports, 446(4):97–216, 2007.
  • [23] 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.
  • [24] Arne Traulsen, Anirvan M Sengupta, and Martin A Nowak. Stochastic evolutionary dynamics on two levels. Journal of theoretical biology, 235(3):393–401, 2005.