Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

This article is part of the supplement: Selected articles from The 5th IEEE International Conference on Systems Biology (ISB 2011)

Open Access Research

Absorbing phenomena and escaping time for Muller's ratchet in adaptive landscape

Shuyun Jiao12 and Ping Ao13*

Author Affiliations

1 Shanghai Center for Systems Biomedicine, State Key Laboratory of Oncogenes and Related Genes, Shanghai Jiao Tong University, 200240, Shanghai, China

2 Department of Mathematics, Xinyang Normal University, 464000, Xinyang, Henan, China

3 Department of Physics, Shanghai Jiao Tong University, 200240, Shanghai, China

For all author emails, please log on.

BMC Systems Biology 2012, 6(Suppl 1):S10  doi:10.1186/1752-0509-6-S1-S10

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/6/S1/S10


Published:16 July 2012

© 2012 Jiao and Ao; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

The accumulation of deleterious mutations of a population directly contributes to the fate as to how long the population would exist, a process often described as Muller's ratchet with the absorbing phenomenon. The key to understand this absorbing phenomenon is to characterize the decaying time of the fittest class of the population. Adaptive landscape introduced by Wright, a re-emerging powerful concept in systems biology, is used as a tool to describe biological processes. To our knowledge, the dynamical behaviors for Muller's ratchet over the full parameter regimes are not studied from the point of the adaptive landscape. And the characterization of the absorbing phenomenon is not yet quantitatively obtained without extraneous assumptions as well.

Methods

We describe how Muller's ratchet can be mapped to the classical Wright-Fisher process in both discrete and continuous manners. Furthermore, we construct the adaptive landscape for the system analytically from the general diffusion equation. The constructed adaptive landscape is independent of the existence and normalization of the stationary distribution. We derive the formula of the single click time in finite and infinite potential barrier for all parameters regimes by mean first passage time.

Results

We describe the dynamical behavior of the population exposed to Muller's ratchet in all parameters regimes by adaptive landscape. The adaptive landscape has rich structures such as finite and infinite potential, real and imaginary fixed points. We give the formula about the single click time with finite and infinite potential. And we find the single click time increases with selection rates and population size increasing, decreases with mutation rates increasing. These results provide a new understanding of infinite potential. We analytically demonstrate the adaptive and unadaptive states for the whole parameters regimes. Interesting issues about the parameters regions with the imaginary fixed points is demonstrated. Most importantly, we find that the absorbing phenomenon is characterized by the adaptive landscape and the single click time without any extraneous assumptions. These results suggest a graphical and quantitative framework to study the absorbing phenomenon.

Background

Muller's ratchet proposed in 1964 is that the genome of an asexual population accumulates deleterious mutations in an irreversible manner. It is a mechanism that has been suggested as an explanation for the evolution of sex [1]. For asexually reproducing population, without recombination, chromosomes are directly passed down to offsprings. As a consequence, the deleterious mutations accumulate so that the fittest class loses. For sexually reproducing population, because of the existence of recombination between parental genomes, a parent carrying high mutational loads can have offspring with fewer deleterious mutations. The high cost of sexual reproduction is thus offset by the benefits of inhibiting the ratchet [2]. Muller's ratchet has received growing attention recently. Most studies of Muller's ratchet are related to two issues. One is that without recombination, the genetic uniformity of the offspring leads to much lower genetic diversity, which is likely to make it more difficult to adapt [3]. So its adaptiveness arouses concern. The other is that population lacking genetic repair should decay with time, due to successive loss of the fittest individuals [4,5]. So the fixation probability arouses concern. In addition, Muller's ratchet is relevant to some replicators [6,7], endosymbionts [8], and mitochondria [9]. In order to assess the relevance of Muller's ratchet, it is necessary to determine the rate (or the time) for the accumulation of deleterious mutations [10]. It is widely recognized that the rate of deleterious mutations being much higher than that of either reverse or beneficial mutations results in a serious threat to the survival of populations at the molecular level [4]. Because models proposed must rest on the biological reality, which must be analyzed on their own without any injection of extraneous assumptions during the analysis [11]. Overall, it has been a long interest to develop a suitable and quantitative theory for the ratchet mechanism and the incidental absorbing phenomenon.

Biologists have suggested [12,10] that a quantitative framework is needed. The potential evolutionary importance of Muller's ratchet makes it desirable to carry out careful quantitative studies [12]. And the incidental absorbing phenomenon is investigated quantitatively in broad literature. The simplest and earliest mathematical model is the pioneering work in Ref. [13]. It described the same evolutionary process on the condition of deterministic mutation-selection balance according to the Wright-Fisher dynamics. And it indicated numerical evidence of relation between the total number of individuals and the average time between clicks of the ratchet, but it did not focus on the absorbing phenomenon. It treated the pioneering model as a diffusion approximation [14], and produced more accurate predictions over the relatively slow regime. It noted that the increasing importance of selection coefficients for the rate of the ratchet for increasing values of the total number of individuals. But it is represented as stochastic differential equations and did not get the predictions over all parameters regions. It employed simulation approaches to Muller's ratchet [15] and estimated how different between the distribution of mutations within a population and a Poisson distribution. But it did not emphasize the absorbing phenomenon. In Ref. [2] it obtained diffusion approximations for three different parameter regimes, depending on the speed of the ratchet. The model shed new light on [14]. But it mainly focused on the property of the solution for these stochastic differential equations. In Ref. [10] it mapped Muller's ratchet to Wright-Fisher process, and got the prediction of the rate of accumulation of deleterious mutations when parameters lie in the fast and slow regimes of the operation of the ratchet. But it put the constraints of Dirac function on the boundary.

Previous works mainly focused on the parameter regimes with lower or higher mutation rates. And models are represented as stochastic differential equation. In Ref. [16] authors imagined the population evolved on an adaptive landscape, but they could not analytically construct it. It described discrete birth-death model and its corresponding diffusion manner by the adaptive landscape [17]. But it did not discuss the absorbing phenomenon. The concept of adaptive landscape is proposed by Sewall Wright to build intuition for the complex biological phenomena [5]. In the present article, inspired by [16,10], we model Muller's ratchet as a Wright-Fisher process, analytically construct the adaptive landscape, where the non-normalizable stationary distribution occurs. Here the adaptive landscape is analytically quantified as a potential function from the physical point of view [18]. We give the adaptive and unadaptive states for the whole parameters region by the adaptive landscape. We give the formula for the single click time of Muller's ratchet in the face of infinite and finite potential. In addition, we can handle the absorbing phenomenon without extraneous assumptions.

The key concept in constructing the adaptive landscape is of potential function as a scalar function. There is a long history of definition, interpretation, and generalization of the potential. Such potential has also been applied to biological systems in various ways. The usefulness of a potential reemerges in the current study of dynamics of gene regulatory networks [19], such as its application in genetic switch [20-23]. The role of potential is the same as that of adaptive landscape. In this article, we do not distinct them.

We now make the obvious advance to Muller's ratchet. We analytically construct the adaptive landscape. We demonstrate the position and adaptiveness of fixed points. This makes the dynamical behaviors of the population to be investigated. In addition, we give the area with imaginary fixed points. This makes the explaining for the imaginary fixed points biologically possible. Infinite potential barriers can be crossed over under some cases. We handle the absorbing phenomenon without any extraneous assumptions under the condition of diffusion approximation. Inversely, we demonstrate the power of the adaptive landscape.

Methods

Discrete model and absorbing boundary

We consider here in population genetics an important and widely applied mechanism- Muller's ratchet. It is the process by which the genomes of an asexual population accumulate deleterious mutations in an irreversible manner [24,25]. It corresponds to the repeated irreversible loss of the fittest class of individuals because of the accumulation of the deleterious mutations, the effective absence of beneficial mutations, without any recombination, but with the random drift [25,12]. Consider a population of haploid asexual individuals with discrete generations t = 0,1, 2,.... The common point in a generation is regarded as an adult stage, after all selection has occurred and immediately prior to reproduction. New mutations occur at reproduction and all mutations are assumed to deleteriously affect viability but have no effect on fertility. Supposed population size is fixed for each generation. The viability of a newly born individual is taken to be determined solely by the alleles they carry. This allows us to divide the population into different classes with different genotypes.

Here in one dimensional case, we consider one locus with two alleles (for example, A and a), that is, there are two classes in the haploid asexual population, one class with allele A while the other with allele a, supposed mutation from allele A to a is deleterious. We assume fixed population size of N, which means we have N alleles in all. We also assume that N >1. Generations are non-overlapping. The lifecycle of the individuals in the population is from adults to juveniles, during which we consider irreversible mutation, selection, and random drift. The frequency of the allele A for generation t is <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M1">View MathML</a> while that of allele a is <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M2">View MathML</a>. Let μ be the probability that an offspring of an adult with allele A is an individual with allele a, labeled by M1,0, that is M1,0 = μ. Analogously, M0,0 = 1 - μ, M0,1 = 0, M1,1 = 1. The relative viability of individuals with allele A is ν0 = 1 while that of individuals with allele a is ν1 = 1 - σ. Where σ can be treated as an effective selection coefficient associated with deleterious mutations. So the values of parameters for μ and σ are from 0 to 1. Then in generation t + 1, when selection and deleterious mutation are active, the probability that the offspring of a parent with allele A is chosen to be with allele a is <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M3">View MathML</a>, the probability that the offspring with allele A is <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M4">View MathML</a>, the probability that the offspring of a parent with allele a is still with allele a is <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M5">View MathML</a>. So the frequency of allele A in generation t + 1 is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M6">View MathML</a>

(1)

Eq.(1) describes the deterministic process that ignores random drift. Under the mutation-selection balance, the fixed points is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M7">View MathML</a>

(2)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M8">View MathML</a>

(3)

This means the population ultimately arrives at the state with allele frequency <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M9">View MathML</a> or <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M10">View MathML</a> and no transition between the two states occur. It is evident that selection rates σ is greater than mutation rates μ. But populations always evolve randomly. Since each individual is assigned a parent independently, if by generation t the average value of allele frequency is pn,t, we have that the transition probability from <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M79">View MathML</a> to <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M80">View MathML</a> by generation t + 1 is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M11">View MathML</a>

(4)

The matrix of transition probabilities is W which is composed of elements Wnm. The dynamical rule can be written as the matrix Eq.(5)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M12">View MathML</a>

(5)

Where P(t) represents the probability distribution of allele A. It is composed of N + 1 elements <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M13">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M14">View MathML</a> denotes the probability that allele A has the frequency n/N in generation t in the presence of deleterious mutation, selection and random drift. The matrix of transition probabilities is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M15">View MathML</a>

(6)

Here vT means the transpose of vector v. Then Eq.(5) can be expressed as the following

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M16">View MathML</a>

(7)

Where 0 is a column vector with N vanishing components, and v is a column vector with N non-zero elements given by υm = W0m, where m is from 1 to N. The column vector w is an N × N matrix with elements wnm = Wnm, where n and m run from 1 to N. By spectrum decomposition we can derive the maximum eigenvector for w, this corresponding vector is the quasi-stationary probability density of allele A. The following is its derivation process.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M17">View MathML</a>

(8)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M18">View MathML</a>

(9)

Eq.(8) has the solution

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M19">View MathML</a>

(10)

Because the leading large time behavior is determined by the eigenvalue of matrix w. From Perron-Frobenius theory this issue is transformed to solve the leading large eigenvalue of w and corresponding eigenvector. But

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M20">View MathML</a>

(11)

Where IT= (1,..., 1), there are N elements 1 of this vector. If we denote the leading large eigenvalue of w is λ1,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M21">View MathML</a>

(12)

We can describe Eq.(8) as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M22">View MathML</a>

(13)

where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M23">View MathML</a>

(14)

Then the probability density is derived by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M24">View MathML</a>

(15)

Where λ1 is determined by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M25">View MathML</a>. But <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M26">View MathML</a>. In the end we get λ1 = 1-vTq.

Populations evolution is a natural and random process. We model the process as the discrete form. The boundary of the discrete model is determined by its transition matrix. Generally, The transition probabilities from generation t to generation t + 1 for the allele frequency are expressed as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M27">View MathML</a>

(16)

From the expression of transition probability, it can be seen that the transition probabilities are zero for any frequency <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M81">View MathML</a> under the condition of parameter σ = 1. It means the population stays at its initial state. In addition, the transition probabilities from the boundary 0 to its next are (1,0,..., 0)T. This means boundary 0 can not output any probability flow to its next, it only absorbs probability from next. We call absorbing phenomenon occurring at the boundary 0.

Continuous model and adaptive landscape

Diffusion approximation

Here we briefly outline the diffusion approximation from the discrete to continuous models. At generation t the frequency of allele A is i/N, after evolutionary force, at generation t + 1 the allele frequency becomes j/N. Here δt = 1, the probability that allele frequency becomes j/N is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M28">View MathML</a>

(17)

Where Wij is the transition probability. Diffusion approximation is a description of the process, valid when N is large, where the allele frequency n/N are replaced by real number x, 0 ≤ x ≤ 1. Given that A starts out at gene frequency x0. Let <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M29">View MathML</a> be the probability for allele frequency A after t generations. And <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M30">View MathML</a> be the probability of allele A after t + 1 generations, then

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M31">View MathML</a>

(18)

Letting <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M32">View MathML</a>, among this ρ(x,t) is the probability density. Define M(x) as the probability that x increases by systematic force that include mutation and selection. And define V(x) as the probability that x changes because of random drift, either decreasing by amount δx with the probability V(x)/2 or increasing by the amount δx with the probability V(x)/2.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M33">View MathML</a>

(19)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M34">View MathML</a>

(20)

In any time interval δt, the probability that x remains at x is 1 - M(x) - V(x). The changes in states are only to δx or -δx. So δx is 0, positive or negative.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M35">View MathML</a>

(21)

and we can treat Eq.(21) as the following

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M36">View MathML</a>

(22)

that is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M37">View MathML</a>

(23)

Because the function about the change of allele frequency in one generation is continuous and smooth enough, under the condition that <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M38">View MathML</a> is of order of magnitude small than <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M39">View MathML</a> in one generation. Put it differently the change of allele frequency is not more than 1/N, and the probability density is smooth enough during the time scale of one generation. So we represent the process as the following approximated diffusion equation

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M40">View MathML</a>

(24)

and according to the definition of M(x) and V(x), the explicit expressions of them are

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M41">View MathML</a>

(25)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M42">View MathML</a>

(26)

Among this M(x) is the symbol for the change in allele frequency [26,11] that occurs in one generation due to systematic force. The function V(x) is the variance in allele frequency after one generation of binomial sampling of N alleles [27].

Adaptive landscape

Under the general diffusion approximation, frequency pn,t is treated as continuous quantities x, and this leads to the distribution of the frequency for the allele A being the probability density. Let ρ(x,t) be the probability density of the frequency for the allele A being x at time t. The diffusion process can be expressed by the following symmetric equation

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M43">View MathML</a>

(27)

with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M44">View MathML</a>

(28)

With a prime denoting differentiation of a function with respect to its argument such as D'(x) = ∂xD(x). Where M(x) and V(x) is from Eqs.(25) and (26) respectively. Adaptive landscape is directly given when we consider natural boundary as Feller's classification. It is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M45">View MathML</a>

(29)

The symmetric Eq.(27) has two advantages. On the one hand, the adaptive landscape is directly read out when the detailed balance is satisfied. On the other hand, the constructive method is dynamical, independent of existence and normalization of stationary distribution. We call f(x) directional transition rate, integrating the effects of M(x) and the derivative of V(x). Directional transition rate can give equilibrium states when it has the linear form.

When the process lies at stationary state, the probability flux of the system is zero, and probability flux flows in x ∊ [0,1]. In general, the stationary distribution for the diffusion approximation satisfying natural boundary condition is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M46">View MathML</a>

It has the form of Boltzmman-Gibbs distribution [28], so the scalar function Φ(x) naturally acquires the meaning of potential energy [19]. The value of Z determines the normalization of ρ(x, t = ∞) from the perspective of probability, and the finite value of Z manifests the normalization of ρ(x,t = ∞). The stationary distribution is not true in the face of infinite Z. It demonstrates absorbing phenomenon occurs at the boundary. Together with the flux at the boundary, the true stationary distribution could be got. The constant holds the same position as temperature of Boltzmman-Gibbs distribution in statistical mechanics. But it does not hold the nature of temperature in Boltzmman-Gibbs distribution.

We are interested in the dynamical property of adaptive landscape, so we treat Φ and Φ/∊ no difference in this respect, that is, for convenience we can take ∊ = 1 of ∊D(x). So according to Eq.(29), we have adaptive landscape as the following

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M47">View MathML</a>

(30)

From the expression of adaptive landscape Φ(x), we may find there are two singular points 0 and 1 of adaptive landscape, characterized by infinite value, infinity means adaptiveness or unadaptiveness of the system. Here the adaptive landscape is composed of three terms. The first term and the third term quantify the effect of the effect of irreversible mutations and selection respectively, the second term quantifies the effect of random drift.

The stationary distribution can be expressed as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M48">View MathML</a>

(31)

Results

Fixed points and their adaptiveness

To understand the mechanism of Muller's ratchet, a full characterization of dynamical process is a

prerequisite for obtaining more accurate decaying time. Here we study the dynamical behaviors by

investigating the position and adaptiveness of all fixed points. We further derive the parameter regions for all possible cases.

According to general analysis of a dynamical system, letting

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M49">View MathML</a>

(32)

we get

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M50">View MathML</a>

(33)

We solved the Eq.(32) and found two fixed points. If we denote

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M51">View MathML</a>

(34)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M52">View MathML</a>

(35)

They are

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M53">View MathML</a>

(36)

For two singular points x = 0,1, if x → 1, and σ ∊(μ, (2- 1)/(2Nμ - μ)), Φ(x) → -∞. So the population is unadaptive at x = 1. When x → 1, and σ ∊ ((2- 1/(2Nμ - μ), 1), Φ(x) → +∞. So the population is adaptive at x = 1. For x → 0, Φ(x) → +∞ in almost parameters regimes except σ = 1. So the population is always adaptive at x = 0. When σ = 1, the viability of the sub-fittest class is zero, so populations stay at the initial state, the corresponding minimum of adaptive landscape demonstrates the state with allele frequency x = 0.

Here we address dynamical behavior by the positions of two real inequivalent fixed points x1 < x2 first.

I) We find two different real fixed points in the regimes of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M54">View MathML</a> and σ ∊ (μ, 1); in the regimes of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M56">View MathML</a> except the regime of με[(2N-1)/4N(N-1)], and δ = (2-1)/(2-μ). We discuss the position between them and the boundary points 0, 1 and adaptiveness of them in the following.

i) 1< x1 < x2

In the regimes of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M57">View MathML</a> and σ ∊ (μ, (2N μ - 1)/(2Nμ - μ)); in the regimes of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M58">View MathML</a>, (2N - 1)/4N(N - 1)) and σ <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59">View MathML</a>, (2-1)/(2Nμ-μ)), the fixed points satisfy 1< x1 < x2. At the same time the singular point x = 1 is adaptive. There is one adaptive state with allele frequency x = 0 in the system. Populations tend to evolve to the adaptive state.

ii) 1 = x1 < x2

In the regions of μ ∊ (1/(2N - 1), (2N - 1)/4N(N - 1)) and σ = (2- 1)/(2Nμ - μ), the two fixed points satisfy x1 = 1, 1< x2. The state with allele frequency x = 1 is unadaptive. There is one adaptive state with allele frequency x = 0 in the system.

iii) 0 < x1 <1< x2

In the regimes of μ ∊ (0,1/(2N - 1)) and σ ∊ (μ, 1); in the regimes of μ ∊ (1/(2N - 1), 1) and σ ∊ ((2- 1)/(2Nμ - μ), 1), the fixed points satisfy 0 < x1 <1, 1< x2. The fixed point x1 is unadaptive. There is only one unadaptive state with allele frequency x = x1 in the system, and two unadaptive states with allele frequency x = 1 and x = 0 occur in the system. Populations tend to evolve to the adaptive states dependent on the position of the initial state. If the initial state with allele frequency is greater than x1, populations tend to evolve to the adaptive state with allele frequency x = 1.

iv) 0 < x1 < x2 <1

In the regimes of μ∊ ((2N - 1)/4N(N - 1), (2N - 1)/(4N - 3)) and σ <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59">View MathML</a>; in the regimes μ ∊ ((2N - 1)/(4N - 3), 1) and σ <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59">View MathML</a>, the fixed points satisfy 0 < x1 < x2 <1. The state with allele frequency x1 is unadaptive while that with allele frequency x2 is adaptive. There are two adaptive states with allele frequency x = 0 and x = x2 and two unadaptive states with allele frequency x = 1 and x = x1 in the system. Populations evolve to which adaptive states dependent on the initial position.

v) 0 = x1 <1< x2

In the regime of μ ∊ (0,1) and σ = 1, the fixed points satisfy x1 = 0, 1< x2. When selection rate σ = 1, the process lies at the initial state because for this case, the viability of the sub-fittest class is zero.

vi) x1 <0 or x2 <0

The case x1 <0 is impossible, and the case x2 <0 is impossible.

II) Then we discuss the case of two equivalent real fixed points x2 = x1.

In the regimes of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60">View MathML</a>, we find two same fixed points

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M61">View MathML</a>

(37)

i) 1 < x1,2

In the regimes of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M58">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60">View MathML</a>, there are two same fixed points satisfy 1 < x1,2, and they are unadaptive. There is one adaptive state with allele frequency x = 0 in the process.

ii) 1 = x1,2

At the two points of ((2N - 1)/4N(N - 1),2N/(2N - 1)2) and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M62">View MathML</a>, there are two same fixed points satisfy x1,2 = 1, and they are unadaptive. There is one adaptive state with allele frequency x = 0 in the process.

iii) 0 < x,1,2 <1

In the regime of μ ∊ ((2N - 1)/4N(N - 1), (2N - 1)/(4N - 3)), <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60">View MathML</a>; in the regimes of μ ∊ ((2N - 1)/(4N - 3), 1) and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M60">View MathML</a>, there are two same fixed points satisfy 0 < x1,2 <1, and they are unadaptive. There is one adaptive state with allele frequency x = 0 in the process.

III) Finally we consider two imaginary fixed points |x1| < |x2|. Where the |.| denotes the length for an imaginary points.

In the regime of <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M63">View MathML</a>, there are two imaginary fixed points in the system. There is only one adaptive state with allele frequency x = 0. Populations always evolve to the adaptive state.

Especially Φ'(x) = 0 is a linear equation, the fixed point is read out from the expression of f(x). We can measure the adaptiveness by the value of adaptive landscape. The bigger the value of adaptive landscape is, the more adaptive the corresponding state would be. The corresponding area of the fixed points in parameters plane is the Figure 1.

thumbnailFigure 1. Relation of fixed points and parameters for the system in all regimes. The regime represented by III i) with parameters regions σ ∊ (σ1, σ2) and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55">View MathML</a> has two imaginary fixed points. The red curve denoted by II) with parameters satisfying σ = σ2 and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M55">View MathML</a> has two equivalent fixed points. The three cases of II) occur in the intervals. The regimes denoted by I) with parameters satisfying σ ∊ (σ25) and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M58">View MathML</a>; satisfying σ ∊ (σ2, σ5) and μ ∊ ((2N - 1)/4N(N - 1), 1); satisfying σ ∊ (σ 15) and μ ∊ (0,1/(2N - 1)); satisfying σ ∊ (σ 15) and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M77">View MathML</a> have two real fixed points. The five cases of I) occur in the regions.

Irreversible mutation, selection and random drift balance

Concretely we divide mutation rates into three regimes. One is with mutation rates μ ∊ (0,1/(2N - 1)) and selection rates σ ∊ (μ, 1). Another is mutation rates μ ∊ (1/(2N - 1), (2N - 1)/(4N(N - 1))) and selection rates σ ∊ (1/(2N - 1), (2- 1)/(2Nμ - μ)). Another is with mutation rates μ ∊ ((2N - 1)/4N(N - 1), 1) and selection rates σ ∊ (2N/(2N - 1)2, (2- 1)/(2Nμ - μ)). The first parameter regimes, a part of regions of I iii), corresponds to the case that mutation rates lie in the lower regime. And the adaptive landscape is U-shape. The second parameters regime, same of I i), is the middle regime. Adaptive landscape with these parameters demonstrates the monotonic decreasing. The third parameters region, same of I iv), demonstrates the clicking process. The adaptive landscape of full parameter regimes is visualized as Figure 2. From the expression and visualization of adaptive landscape Φ(x), we may find there are two singular points 0 and 1 of adaptive landscape, characterized by infinity values in Figure 2. Singularity from the maximum of adaptive landscape indicates the population being adaptive while singularity from the minimum of adaptive landscape indicates the population being unadaptive. Figure 2 demonstrates the whole process of the population evolution including the forming and losing the fittest class except σ = 1. Because the viability of the class with allele a is not zero and deleterious mutation is arbitrary. The adaptive state under the condition of the parameter σ = 1 only means the initial states. With increasing selection rates the fittest class A forms quickly while with increasing mutation rates the fittest class A loses. In the lower mutation rates regime black dotline describes the population is likely to move to the fittest class with increasing selection rates, the process is dominated by selection. Dashed line and black line manifest the losing process of allele A. Because the mutation rates are lower, selection rates dependent of mutation rates are lower, these factors result in the change of fittest class is not easy. In the end there are two adaptive states in the process. In the higher mutation rates regime, black dotline and black line describe the population is likely to move to the fittest class so that the population exists in the form of coexistence of A and a. Because deleterious mutation rates are higher, as a consequence allele a occurs. Selection rates dependent of mutation rates tends to survive allele A. But the evolutionary process is dominated by the irreversible mutations, the fittest class A loses. There are two adaptive states in the process under the balance of irreversible mutation and selection. In the middle mutation rates and selection rates it demonstrates the losing process. So we can draw the conclusion that the click process occurs when there are two stable states in the process.

thumbnailFigure 2. Adaptive landscape against allele frequency x with mutation rates and selection rates in all regimes. x label represents allele frequency μ label represents mutation rates while vertical label is corresponding value of Φ. Assume population size is constant N=50. The lower regime corresponds to a part of I iii). In the region fixed points satisfy 0 < x1 <1< x2, dashed line represents μ=0.000005, σ=0.00005, black dotline stands for μ=0.000005, σ=0.01010, black line represents μ=0.01, σ=0.01012. The medium mutation rates regime corresponds to the region I i) with fixed points satisfy 1< x1 < x2, black dotline represents μ=0.0101015, σ=0.010102, dashed line stands for μ=0.0101015, σ=0.0102, black line investigates μ=0.0101015, σ=0.01015. The higher mutation rates regime corresponds to the region I iv) with fixed points satisfy 0 < x1 < x2 <1, black dotline represents μ=0.4, σ=0.9, dashed line stands for μ=0.02, σ=0.1, black line describes the case for μ=0.02, σ=0.05.

Characterization of the single click time

We visualize the adaptive landscape, then one may wonder about how the population moves from one peak to another and how long it might be to move from one maximum to another. The process was first visualized by Wright in 1932. In addition, the problem of transition from metastable states is ubiquitous in almost all scientific areas. Most of previous works encounter finite potential barriers from the physical point of view. Interesting issue here is that we touch upon infinite potential barriers under the circumstance of well defined two adaptive states. Then we manifest the derivation of a single click time. The time of a click of the ratchet is recognized as the random time of loss of the fittest class [10]. The single click time is well defined when there are two fittest classes in the process. It means the interval time between extinction of the two fittest classes. The corresponding processes are that there are two well-defined adaptive states in the system. Corresponding graphs of adaptive landscape is Figure 3. To evaluate the single click time and show the further power of adaptive landscape, in the following we will demonstrate how the single click time from one adaptive state to another is derived in this framework.

thumbnailFigure 3. Adaptive landscape with two adaptive states. N=50, blue line represents μ=0.02,σ=0.05 corresponding the case for I iv) while green line stands for μ = 0.000005, σ=0.00005 corresponding the case for I iii).

After straightforward calculation, backward Fokker- Planck equation corresponding to Eq.(27) can be expressed with the property of time homogeneous in the following form [29,30]

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M64">View MathML</a>

(38)

General single click time dependent on initial Dirac function satisfies

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M65">View MathML</a>

(39)

With

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M66">View MathML</a>

(40)

Above treatment is valid. Because populations evolution is according to Muller's ratchet, that is in the presence of deleterious mutation, without any recombination, but with selection and random drift. And the model in discrete manner demonstrates the transition probabilities are 0 from the boundary x = 0 to its next. So the boundary x = 0 only absorbs flux from its next, the boundary is absorbing. The probabilities from boundary x = 1 is not zero because μ ≠ 0 and μ, < σ ≠ 1. And the population can not be out of the boundary x = 1. So the boundary x = 1 is reflecting. The general solution corresponding to Eq.(39) is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M67">View MathML</a>

(41)

here Φ(x) = x(f(x')/D(x'))dx'(∊ = 1).

Here the evolutionary process occurs when x ∊ [0,1]. We are more interested in the transition time between the two adaptive states x = 0 and x = 1. In the process, there are two important states x*, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68">View MathML</a> Interval (0,1) contains a potential well at x* and a potential barrier at <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68">View MathML</a>. The single click time is composed of two elements, one denotes forming process of fittest class, the other describes losing process of fittest class. In general, the time spent on forming process is much smaller than that spent on losing process. So the transition time approximates to the time spent on losing process. Because we assume that near <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68">View MathML</a> we can write

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M69">View MathML</a>

(42)

and near x*

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M70">View MathML</a>

(43)

At the same time, if the central maximum of Φ(x) is large compared with 1/N, then exp(Φ(z)) is sharply peaked at <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68">View MathML</a>, while exp(-Φ(y))/D(y) is very small near y = x*. Eq.(41) is evaluated as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M71">View MathML</a>

(44)

From the expression of Eq.(44), the single click time is not sensitive to the assumption of Eq.(40). In the higher mutation rates regime, where <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68">View MathML</a> approximates to a adaptive state which is near enough to 1, x* corresponds to the unadaptive state that the population lies between the adaptive states 0 and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68">View MathML</a>. The potential barrier is finite. According to classical derivation corresponding to Eq.(44) the single click time approximates to

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M72">View MathML</a>

(45)

Here <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M68">View MathML</a> is the fixed point x2, and x* is the fixed point x1, parameters α, β is the same as Eqs.(34) (35) respectively. The difference of potential is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M73">View MathML</a>

(46)

Where α and β are the same as Eqs. (34) and (35). The approximated single click time varies with mutation rates in Figure 4. The single click time T1increases with population size N in certain regime, decreases with mutation rates μ and selection rates σ in the parameters regime μ ∊ (2N/4N(N - 1), 1) and σ <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M59">View MathML</a>. Because in the regime, with selection rates increasing, the difference of potential between two fixed points decreases, the viability of sub-fittest class decreases, populations evolve to the fittest class. These results in the the single click time shorter. In another hand, with deleterious mutation increasing, the population of sub-fittest class increases, the difference of potential between two fixed points decreases. These also results in the single click time shorter.

thumbnailFigure 4. The approximated single click time decreases with mutation rates increasing in the regime μ ∊ (2N/4N(N - 1), 1) and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M78">View MathML</a> denoted by I iv). Assume N=100, blue, black and dashed lines, corresponding high mutation rates regime, represent μ=210/(400*99), μ=500/(400*99) and μ=1000/(400*99) respectively. Red and green lines, corresponding low mutation rates regime, represent μ=210/(400*99) and μ=500/(400*99) respectively. The approximated single click time increases with mutation rates increasing in the regime μ (1/(2N-1),1) and σ ((2Nμ-1)/(2Nμ-μ),1) denoted by I iii). Green dotline represents this case.

For the lower mutation rates regime, where the potential barrier is infinite. The single click time can be estimated also, x* corresponds to the fixed point x1 that the population lies at the lowest potential.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M74">View MathML</a>

(47)

From expression of Eq.(47), the single click time goes to infinity with mutation rates tends to zero in the parameters regimes of μ ∊ (0,1/(2N - 1)) and σ ∊ (μ, 1). From Figure 4, when parameters regions lie μ ∊ (1/(2N - 1), 1) and σ ∊ ((2- 1)/(2Nμ - μ), 1), the results of the single click time is not sensitive to the population size. Biologically if deleterious mutation accumulates, the viability of sub-fittest class increases, these results in the single click time longer.

Analogous to the derivation of T10, we can calculate

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M75">View MathML</a>

(48)

Compared with the singular point x = 1, the difference between two singular points x = 0 and x = 1 is the mutation rates μ. This results in the power of (1-z) is not negative, so the single click time is finite from it. The infinity of the single click time from x = 0 comes from that the mutation from unfavored allele a to favored allele A is zero. This results in the second integral nonintegrable because of the negative power of z. Biologically because the absence of back mutation, the accumulation of deleterious mutation, once the population arrives at the state that almost individuals are with allele a, the population is absorbed the state and can not leave with high probability. Mathematically, because the second integral is singular for the singular point x = 0. And the integrated function is a fraction respect to argument x, but the highest power of denominator is smaller than that of numerator. That results in the power of x is 2N - 2. As a consequence the second integral is singular.

Discussion

We analytically construct adaptive landscape. The constructive method is independent on the existence and normalization of stationary distribution. We demonstrate the position and adaptiveness of all fixed points for the whole parameters regimes under the condition of the diffusion approximation. An interesting thing is the imaginary fixed points occurring. We give the parameters regions of their occurrence. However, we have not found any study of Muller ratchet for the fixed points to give a complete description. In addition, we give the description of escape from infinite potential. However, intuitively infinite potential means the population lies at adaptive state. The transition from the adaptive state can not occur. Here we find that the escape from infinite potential can not occur when the boundary is absorbing. So we define the absorbing boundary by adaptive landscape and the single click time without any extraneous assumptions.

The model with discrete manner describes the nature of populations evolution. Here we give two special cases. One is that the population lies at that state with allele frequency x = 0, the other is that the population lies at the state with allele frequency x = 1. We compare the model with discrete and continuous manners to conclude the definition of boundary conditions without any extraneous assumptions. The model with continuous manner is derived under the condition of enough small space change δx in one generation, in another word, when population size N is much bigger. The model with continuous manner can correspond to the model with discrete manner.

When allele frequency x = 0, the matrix for transition probabilities demonstrates that the population only absorbs the flux from the next, and it can not output any flux. In the model with continuous manner the value of adaptive landscape with allele frequency x = 0 arrives at the maximum except σ = 1. This demonstrates the state with allele frequency x = 0 is almost always adaptive. Then we calculate the single click time from the boundary with allele frequency x = 0. We find the single click time is infinite. So we draw the conclusion that the state with allele frequency x = 0 absorbs flux, and does not output flux. The boundary is absorbing. When the allele frequency x = 1 in the model with discrete manner, the transition probability is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/S1/S10/mathml/M76">View MathML</a>

that is, when μ >0, the transition probabilities are not zero, the state with allele frequency x = 1 can input any flux to its next state. The single click time from the state is finite in the model with diffusion manner, however the potential at the state could be infinite, then we draw the conclusion this boundary is not absorbing. This is consistent with the biological understanding.

This article presents an approach to estimate the single click time of Muller's ratchet. Furthermore, it define the absorbing phenomenon by the single click time without any extraneous assumptions. Inspired by [16,10], we connect Muller's ratchet to one locus Wright-Fisher model with asexual population including N haploid individuals. And our model is represented as a Fokker-Planck equation. We give a complete description for the position and adaptiveness of all fixed points in the whole parameters regimes. This is first done bases on diffusion approximation. The investigated elements is at the allele level. This is different from Ref. [10]. Our method does not need the existence and normalization of the stationary distribution. Our constructive method is independent of the stationary distribution. Compared with the method based on diffusion approximation [15,2], mathematically it is described as stochastic differential equations. Our method investigates the global dynamical property of the system, and reduces the complexity of calculating stochastic differential equations. In addition, the boundary condition of these stochastic differential equations is prescribed. Compared with Ref. [10], They added Dirac function to the boundary. But this is not appropriate for the adding non-differential Dirac function to stationary distribution, and stationary distribution should satisfy diffusion equation. However, the treatment is convenient for computing the stationary distribution. The stationary distribution of theirs is equivalent to our adaptive landscape. They had not given the shape of adaptive landscape when the mutation rate lies in the lower regime. We use the model defined in the interval (0, 1) to describe the absorbing boundary. We check the biological phenomenon by the model with both discrete and continuous manners. This is a new method to handle the boundary condition. We investigate the absorbing phenomenon by it without any extraneous assumptions.

To summarize, we have obtained two main sets of results in the present work. Most importantly, we find that the absorbing phenomenon is characterized by the adaptive landscape and the single click time without any extraneous assumptions. First, we demonstrate the adaptive landscape can be explicitly read out as a potential function from general diffusion equation. This not only allows computing the single click time of Muller's ratchet straightforward, but also characterizes the whole picture of the ratchet mechanism. The adaptive landscape has rich structures such as finite and infinite potential, real and imaginary fixed points. We analytically demonstrate the adaptive and unadaptive states for the whole parameters regimes. We find corresponding parameters regimes for different shapes of adaptive landscape. Second, we give the formula about the single click time with finite and infinite potential. And we find the single click time increases with selection rates and population size increasing, decreases with mutation rates increasing. These results give a new understanding of infinite potential and allow us a new way to handle the absorbing phenomenon. In this perspective our work may be a starting point for estimating the click time for Muller's ratchet in more general situations and for describing the boundary condition. Such demonstration suggests that adaptive landscape may be applicable to other levels of systems biology.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

Shuyun Jiao carried the research and writing. Ping Ao oversaw the whole project, and participated in research and writing.

Acknowledgements

We would like to thank Yanbo Wang for drawing the figures, also thank Quan Liu for discussions and technical help, thank Song Xu for technical help. We thank Bo Yuan for some advice on writing and correcting some expression on language. This work was supported in part by the National 973 Projects No. 2010CB529200 (P.A.), and in part by No. 91029738 (P.A.) and No.Z-XT-003 (S.J.) and by the project of Xinyang Normal university No. 20100073 (S.J.).

This article has been published as part of BMC Systems Biology Volume 6 Supplement 1, 2012: Selected articles from The 5th IEEE International Conference on Systems Biology (ISB 2011). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/6/S1.

References

  1. Maynard Smith J: The Evolution of Sex. Cambridge University Press; 1978. OpenURL

  2. Etheridge A, Pfaffelhuber P, Wakolbinger A: How often does the ratchet click? Facts, heuristics, asymptotics. In In Trends in Stochastic Analysis. London Math. Soc. Lecture Note Ser. 353. Cambridge Univ. Press; 2009:1091-1148. OpenURL

  3. Lampert KP, Schartl M: A little bit is better than nothing: the incomplete parthenogenesis of salamanders, frogs and fish.

    BMC Biol 2010, 8:78. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  4. Maia LP: Analytical results on Muller ratchet effect in growing populations.

    Phys Rev E Stat Nonlin Soft Matter Phys 2009, 79:032903. PubMed Abstract | Publisher Full Text OpenURL

  5. Barton NH: Mutation and the evolution of recombination.

    Phil Trans R Soc B 2010, 365:1281-1294. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  6. Gordo I, Charlesworth B: The degeneration of asexual haploid populations and the speed of Muller's ratchet.

    Genetics 2000, 154:1379-1387. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  7. Gordo I, Charlesworth B: On the speed of Muller's ratchet.

    Genetics 2000, 156:2137-2140. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. Moran NA: Accelerated evolution and Muller's ratchet in endosymbiont bacteria.

    Proc Natl Acad Sci USA 1996, 93:2873-2878. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  9. Bergstrom C, Pritchard J: Germline bottlenecks and the evolutionary maintenance of mitochondrial genome.

    Genetics 1998, 149:2135-2146. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  10. Waxman D, Loewe L: A stochastic model for a single click of Muller's ratchet.

    J Theor Biol 2010, 264:1120-1132. PubMed Abstract | Publisher Full Text OpenURL

  11. Ewens WJ: Mathematical Population Genetics. Cornell University Press; 2004. OpenURL

  12. Felsenstein J: The evolutionary advantage of recombination.

    Genetics 1974, 78:737-756. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. Haigh J: The accumulation of deleterious genes in a population-Muller's Ratchet.

    Theor Popul Biol 1978, 14:251-267. PubMed Abstract | Publisher Full Text OpenURL

  14. Stephan W, Chao L, Smale JG: The advance of Muller's ratchet in a haploid asexual population: approximate solutions based on diffusion theory.

    Genet Res 1993, 61:225-231. PubMed Abstract | Publisher Full Text OpenURL

  15. Gessler DDG: The constraints of finite size in asexual populations and the rate of the ratchet.

    Genet Res 1995, 66:241-253. PubMed Abstract | Publisher Full Text OpenURL

  16. Woodcock G, Higgs PG: Population evolution on a multiplicative single-peak fitness landscape.

    J Theor Biol 1996, 179:61-73. PubMed Abstract | Publisher Full Text OpenURL

  17. Zhou D, Qian H: Fixation, transient landscape, and diffusion dilemma in stochastic evolutionary game dynamics.

    Phys Rev E Stat Nonlin Soft Matter Phys 2011, 84:031907. PubMed Abstract | Publisher Full Text OpenURL

  18. Ao P: Laws in Darwinian evolutionary theory.

    Phys Life Rev 2005, 2:117-156. Publisher Full Text OpenURL

  19. Ao P: Potential in stochatic differential equation: novel construction.

    J Phys A: Math Gen 2004, 37:25-30. Publisher Full Text OpenURL

  20. Zhu XM, Yin L, Hood L, Ao P: Robustness, stability and efficiency of phage lambda genetic switch: dynamical structure analysis.

    J Bioinform Comput Biol 2004, 2:785-817. PubMed Abstract | Publisher Full Text OpenURL

  21. Zhu XM, Yin L, Hood L, Ao P: Calculating biological behaviors of epigenetic states in the phage lambda life cycle.

    Funct Integr Genomics 2004, 4:188-195. PubMed Abstract | Publisher Full Text OpenURL

  22. Liang J, Qian H: Computational cellular dynamics based on the chemical master equation: A challange for understanding complexity.

    J Comput Sci Technol 2010, 25:154-168. Publisher Full Text OpenURL

  23. Cao YF, Lu HM, Liang J: Probability landscape of heritable and robust epigenetic state of lysogeny in phage lambda.

    Proc Natl Acad Sci USA 2010, 107:18445-18450. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. Muller HJ: Some genetic aspects of sex.

    Am Nat 1932, 66:118-138. Publisher Full Text OpenURL

  25. Muller HJ: The relation of recombination to mutational advance.

    Mutat Res 1964, 106:2-9. PubMed Abstract | Publisher Full Text OpenURL

  26. Kimura M: Diffusion models in population genetics.

    J Appl Prob 1964, 1:177-232. Publisher Full Text OpenURL

  27. Gillespie JH: Population Genetics: a Concise Guide. The Johns Hopkins University Press; 2004. OpenURL

  28. Ao P: Emerging of Stochastic Dynamical Equalities and Steady State Thermodynamics from Darwinian Dynamics.

    Commun Theor Phys 2008, 49:1073-1090. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  29. van Kampen NG: Stochatic Processes in Physics and Chemistry. Amsterdam Press; 1992. OpenURL

  30. Øksendal B: Stochatic Differential Equations: an Introduction with Applications. Springer; 2003. OpenURL