# On some dynamical features of
the complete Moran model for neutral evolution
in the presence of mutations

G. Gaeta

Dipartimento di Matematica, Università degli Studi di Milano
v. Saldini 50, 20133 Milano (Italy)

giuseppe.gaeta@unimi.it
and

INFN, Sezione di Milano

12/2/2024

## Abstract

We present a version of the classical Moran model, in which mutations are taken into account; the possibility of mutations was introduced by Moran in his seminal paper, but it is more often overlooked in discussing the Moran model. For this model, fixation is prevented by mutation, and we have an ergodic Markov process; the equilibrium distribution for such a process was determined by Moran. The problems we consider in this paper are those of first hitting either one of the “pure” (uniform population) states, depending on the initial state; and that of first hitting times. The presence of mutations leads to a nonlinear dependence of the hitting probabilities on the initial state, and to a larger mean hitting time compared to the mutation-free process (in which case hitting corresponds to fixation of one of the alleles).

# 1 Introduction

The Moran model [1–3] is a classical tool in studying *neutral evolution* [4–12], i.e. the evolution of genome when several alleles are present but none of these presents any advantage over the others; thus in this framework evolution – and in particular fixation of a single allele, or at intermediate times elimination of some of the different alleles – is due uniquely to random fluctuations in the reproduction rates of individuals carrying the different alleles. Needless to say, the relevance of these fluctuations is intimately related to the fact any population is of *finite size*.

The Moran model deals with the simple case in which only two alleles are present at a locus, and the population size is fixed; this captures the essen-
tial features of neutral evolution yet is simple enough to be amenable to exact solution.

On the other hand, the standard Moran model only considers reproduction fluctuations, and does not take into account the presence of mutations.

It should be mentioned in this regard that the original paper by Moran [1] also took mutations into account and provided results for this case, see below for detail. Most of the literature dealing with the Moran mode (including that considering generalizations of the Moran model, see e.g. [13–30]), however, disregards the possibility of mutations; we will hence refer to the Moran model without mutations as the standard Moran model.

Taking into account mutations does of course raise immediately a substantial problem. That is, if mutations are possible, fixation becomes impossible. The question to be asked concerns then on the one hand the equilibrium distribution of the process (this was determined by Moran [1]) and on the other the dynamics of the model, in particular, given an initial state, the probability first hitting either one of the extreme states, i.e. those corresponding to a uniform population and which will therefore also called pure states (all other states will be called mixed states); and the hitting times for these. That is, we are led to study a problem of stopping times, or of exit from a domain; in this case, the domain is that of “mixed states”, and exit from this domain corresponds to hitting of the “pure states”, representing the border of the set of allowed states.

In the present paper, we consider a version of the Moran model, in which mutations are taken into account, and evolution – in the same sense as in the standard Moran model – is the result of fluctuations in the mutation rates as well as in reproductive rates for the two alleles; this will hence be dubbed as the “complete” Moran model (as opposed to the “standard” Moran model, see above). As already mentioned, this is not a new model, but corresponds to the full model considered by Moran in 1957. We will then consider the two problems mentioned above (first hitting either one of the pure states and hitting times), see below for more detail.<fnref n="1" />

**Remark 1.** When mutation is present, the Moran process is an ergodic Markov chain, hence it has a unique equilibrium distribution; this is attributed by Moran to Wright (see eq.(9) in Moran’s paper [1]), and is given by

$$P(x) = \frac{\Gamma(\beta_1 + \beta_2)}{\Gamma(\beta_1) \Gamma(\beta_2)} (1 - x)^{\beta_1 - 1} x^{\beta_2 - 1}, \quad (1)$$

where

$$x = k/n, \quad \beta_i = n \alpha_i,$$

with n the population size and αᵢ the mutation rates for the allele i = 1, 2, i.e. the probability that a descendant of an individual with allele i has the different allele. ⊙

<footnote n="1">¹These problems were already considered in the literature, but the available results are mostly of numerical nature (we will not discuss them); we will instead adopt an analytic approach, and provide exact formulas – albeit in recursive form. See however the Appendix for closed-form expressions in terms of special functions.</footnote>
When discussing neutral evolution one is of course interested in the case where mutations do not favor any of the alleles, hence – in the framework where only two alleles are present, as in the Moran model – they are symmetric in the two alleles. We will thus consider, in the notation of Remark 1, the case $α_1 = α_2 = α$, which also entails $β_1 = β_2 = β$.

The complete Moran model we consider is amenable to a limited random walk, like many biological models [31,32]; at difference with the standard Moran model, however, it does not correspond to a symmetric random walk (even for symmetric mutation rates).

We will state more precisely the questions we will investigate:

1. How the introduction of mutations affects the *first hitting probabilities* (that is, the probability that one “pure” state, in which all the population shares the same allele, is reached before the other one) for the two alleles, and more precisely their dependence on the initial state $k$, $0 < k < n$;

2. How the introduction of mutations affects the *mean hitting time* (this is the time at which either one of the “pure” states is first reached) for a process starting from the initial state $k$, $0 < k < n$.

It turns out that for both questions one could obtain analytical formulas in terms of hypergeometric functions; these explicit formulas are however rather involved and of little practical use (they are given in the Appendix). On the other hand, the two questions mentioned above may also be discussed in terms of exact (analytical, not numerical) computations for arbitrary but given population number $n$; these will allow to get to the point with minimal mathematical complications.

**Remark 2.** Needless to say, the “complete” model considered in this note (like the standard Moran model usually considered in the literature) is far too simple to be directly applicable to concrete current biological problem; it should instead be seen as a way to get a qualitative – and partially quantitative – understanding of the balance among reproductive versus mutational fluctuations, and their role in neutral evolution. $\circled{0}$

**Remark 3.** Finally, we note that the fixation problem can of course also be considered when one of the two alleles is advantageous, in particular in the case of beneficial mutations; this framework is outside of the scope of the present work, and is discussed e.g. in the review paper by Patwa and Wahl [33]. $\circled{0}$

# 2 The standard Moran model

In the standard Moran model [1–3], we have a haploid population of $n$ individuals, each of them carrying an allele *A* or *B* of the only gene under study. At each time step, two individuals from the population are chosen (they may coincide). One is allowed to reproduce, and the other is eliminated; the one which is reproducing transmits the allele *A* or *B* to its heir. Denoting by $n_A = k$,
$n_B = n - k$ the number of individuals carrying the alleles A and B respectively this dynamics means that at each step, depending on which alleles are carried by the individuals chosen for reproduction and for elimination, the number $k$ can either remain unchanged, $\hat{k} = k$, or change to $\hat{k} = k \pm 1$. Denoting by $p_k$ the probability for $k \to \hat{k} = k + 1$ and by $q_k$ the probability for $k \to \hat{k} = k - 1$, it is easily seen that

$$p_k = q_k = \frac{k}{n} \left(1 - \frac{k}{n}\right). \qquad (2)$$

Note that for $k = 0$ or $k = n$, there will be no further change in the composition of the population: all its members carry the same allele (B for $k = 0$, A for $k = n$), and this is definitely fixed.

The fixation time $T(k, n)$, i.e. the mean time to fixation for a process starting from state $k$ in a population of $n$ individuals, is computed to be

$$T(k, n) = n \left[ \sum_{m=1}^{k} \frac{n-k}{n-m} + \sum_{m=k+1}^{n-1} \frac{k}{m} \right]. \qquad (3)$$

For large $n$ we have $T(k, n) \approx \hat{T}(k, n)$, where

$$\hat{T}(k, n) := -n^2 \left[ \left(1 - \frac{k}{n}\right) \log \left(1 - \frac{k}{n}\right) + \frac{k}{n} \log \left(\frac{k}{n}\right) \right]. \qquad (4)$$

## 3 The Moran model in the presence of mutations

The standard Moran model aims at illustrating how fixation emerges as a result of fluctuations in the reproductive rate of the two alleles; on the other hand, in general there will also be mutations making that descendants of individuals with one allele will actually have the other allele. Even in the case where the mutation rates are the same for the $A \to B$ and the $B \to A$ mutations, there will be fluctuations in mutations rates – exactly as in reproduction rates – and these can also lead to reaching transiently a “pure” state (fixation is not possible, since mutations can halter a population sharing the same allele).

More precisely, in the neutral case of interest here and in the presence of mutations, reaching a pure state will result from a combination of reproduction and mutation fluctuations. Thus we would like to understand how the combination of the reproductive and mutation fluctuations affects the dynamics, and in particular the hitting of the extreme, pure states, in neutral evolution; in particular, how the hitting probabilities for the two pure states and the mean hitting time respond to the presence of mutations.

We will thus consider a simple Moran-like model in which mutations are also present; this will be dubbed the “complete” Moran model (as mentioned in the Introduction, this was actually considered by Moran in his seminal 1957 paper).

The process undergone by this model will be described as follows. We have a homozygous population of $n$ individuals, and the gene looked upon can exist
in two alleles, A and B; the state of the system is labeled by the number k of
individuals with the A allele (correspondingly, there will be n − k with the B
allele).

At each step, several operations are performed, i.e.:

(i) A single individual is chosen, with equal probability for each of the n
individuals;

(ii) The chosen individuals reproduces, giving birth to a single descendant; this has the same allele of the parent with probability (1 − μ), and the other allele with probability μ ;

(iii) Then again a single individual is chosen from the pre-existing population, again with uniform probability, and this is eliminated (the newborn is excluded from this choice, but the chosen parent can coincide with the eliminated individual).

We stress that the process is quite obviously a Markov chain, and actually
an ergodic one. Note that for $\mu$ = 0 we recover exactly the standard Moran
model (which is not ergodic).

It is clear that the population size remain constant, by construction; more-
over at each step the state *k* can only be changed to *k* = *k* ± 1, or remain equal
(unless *k* = 0, in which case it cannot pass to *k* − 1; or *k* = *n*, in which case
it cannot pass to *k* + 1). We thus have to evaluate the probabilities for these
three events. Actually we are only interested in evaluating such probabilities
for *k* ≠ 0 and *k* ≠ *n*: in fact once either one of the pure states is reached we
stop following the process.

We note that (for $k \neq 0, n$) the process can jump to $\hat{k} = k + 1$ in two ways: either the reproducing individual is of type A, mutation does not take place in the reproduction, and the eliminated individual is of type B; or the reproducing individual is of type B, mutation does occur in reproduction, and the eliminated individual is again of type B. The probability for these to occur are easily evaluated, and the transition probability $P(k \rightarrow k + 1) := p_k$ is just their sum.

Similarly, the process can jump (for $k \neq 0, n$) from $k$ to $\hat{k} = k - 1$ in two ways: either the reproducing individual is of type B, mutation does not take place in the reproduction, and the eliminated individual is of type A; or the reproducing individual is of type A, mutation does occur in reproduction, and the eliminated individual is again of type A. The probability for these to occur are easily evaluated, and the transition probability $P(k \rightarrow k - 1) := q_k$ is just their sum.
In explicit terms,

$$ 
\begin{aligned}
p_k &= \frac{k}{n} (1 - \mu) \left(\frac{n-k}{n}\right) + \left(\frac{n-k}{n}\right) \mu \left(\frac{n-k}{n}\right) \\
&= \left(1 - \frac{k}{n}\right) \left[\mu + \frac{k}{n} - 2\mu \frac{k}{n}\right];
\end{aligned} 
\qquad (5) $$

$$ 
\begin{aligned}
q_k &= \left(\frac{k}{n}\right) \mu \left(\frac{k}{n}\right) + \left(\frac{n-k}{n}\right) (1-\mu) \left(\frac{k}{n}\right) \\
&= \left(\frac{k}{n}\right) \left[1 - \left(\mu + \frac{k}{n} - 2\mu \frac{k}{n}\right)\right].
\end{aligned} 
\qquad (6) $$

Note again that for $\mu \to 0$ we recover the probability transitions of the standard Moran model.

At difference with the standard Moran model, we have in general (the exception being for $k = n/2$, if $n$ is an integer) $p_k \neq q_k$, i.e. a non-symmetric random walk.

In this case we have for the drift $\delta k = \hat{k} - k$

$$ \langle \delta k \rangle = p_k - q_k = \mu \left(1 - 2 \frac{k}{n}\right); \qquad (7) $$

thus for $\mu \neq 0$ this dynamics drives the process towards equilibrium, and we expect that hitting the borders requires larger times than in the case of the standard Moran model.

In the following, the ratio $\xi_k = q_k/p_k$ will be of interest; it follows from (5) and (6) and simple algebra that this is given by

$$ \xi_k = \frac{q_k}{p_k} = \frac{k(2\mu - 1) - \mu n + n}{(k - n)(-2\mu k + k + \mu n)}. \qquad (8) $$

Writing $k/n = x \in [0, 1]$, this also reads

$$ F(x) = \left(\frac{x}{x-1}\right) \left(1 + \frac{1}{2\mu x - x - \mu}\right); \qquad (9) $$

this is plotted as a continuous function in Figure 1 (this kind of representation is justified for large $n$, and will be used again in the following).

It should be noted that we always have $F(1/2) = 1$ and moreover, unless $\mu$ is large, $F(x)$ is substantially different from the unit function only for $x \simeq 0$ or $x \simeq 1$; in fact, expanding $F$ in a $\mu$ series, we get

$$ F(x) = 1 - \frac{(1-2x)}{x(1-x)} \mu + O(\mu^2). \qquad (10) $$

Thus for a small mutation rate $\mu \ll 1$ the properties of the model – or at least those controlled by the parameter $\xi_k$, which we will see is a central one – are expected to be quite similar to those of the standard Moran model, except near the pure states $k = 0$ and $k = n$ (corresponding to $x = 0$ and $x = 1$).
<figure><img src="image_1.png" /><figcaption>Figure 1: The function $F(x)$ given by (9), plotted here for $\mu = 0$ (black, constant and equal to unity), $\mu = 10^{-3}$ (blue), $\mu = 10^{-2}$ (green), and $\mu = 10^{-1}$ (red). For small but nonzero $\mu$, the curve flattens on $F = 1$ (the behavior of the standard Moran model) except near the borders $\alpha = 0$ and $\alpha = 1$; that is – going back to $k$ and $n$ – except for $k \approx 0$ and $k \approx n$.</figcaption></figure>

The behavior in this region is of course essential for actually reaching a pure state, so this deviation in a limited region will suffice to provide rather different hitting properties, in particular concerning the mean hitting time.

We also note that the $F(x)$ with smaller $\mu$ will be more flat in the central region, and steeper at the borders, see Figure 1. They become of course singular in the limit $\mu \to 0$.

**Remark 4.** Note that one could in principles also consider a model in which only mutational fluctuations – and not reproductive ones – are at play; this would be a “purely mutational” Moran model. This has some mathematical interest, but is less interesting biologically, in that in the limit of vanishing mutation rate it yields a trivial model in which whatever initial composition of the population remains fixed. $\circ$

# 4 Hitting probability: a general discussion

We want to evaluate the hitting probability $\pi(k_0)$ and $\chi(k_0)$, depending on the initial state $k_0$, for the two alleles; these are respectively the probability $\pi(m)$ that if the initial state is $k_0 = m$, the system hits the state $k = n$ before passing through the state $k = 0$; and the probability $\chi(m)$ that it hits the state $k = 0$ before passing through the state $k = n$. We note preliminarily that by their
definition, these quantities satisfy

$$ \pi(0) = 0, \quad \pi(n) = 1; \quad \chi(0) = 1, \quad \chi(n) = 0. \tag{11} $$

We should also note that, in principles, the system could also oscillate indefinitely along the (finite) lattice; thus we are not guaranteed *apriori* that $\pi(m) + \chi(m) = 1$.

It is well known that actually random walks on a finite one-dimensional lattice always hit the border (in the case of absorbing boundaries, as for the standard Moran model, this means reaching fixation), as the measure of forever oscillating paths is null. It will be useful to give a simple scheme to prove this, and to actually compute $\pi(k)$ and $\chi(k)$.

We consider a general random walk on the lattice $\mathbb{Z}_{n+1} = \{0, 1, \dots, n\}$, with transition probabilities $p_k$ for $k \to k+1$, $q_k$ for $k \to k-1$, and $r_k$ for $k \to k$, with $p_k + q_k + r_k = 1$, with $p_k \neq 0 \neq q_k$ for all $k$ but $k=0$ and $k=n$; $\pi(k)$ and $\chi(k)$ are as above. The situation is trivial for $k_0 = 0$ (by definition, $\pi(0) = 0$, $\chi(0) = 1$) and for $k_0 = n$ (by definition, $\pi(n) = 1$, $\chi(1) = 0$), so we will from now on assume $k_0 \neq 0, n$.

In order to compute $\pi(k)$ and $\chi(k)$, $1 \le k \le n-1$, one filters the process according to its first step. We will first focus on $\pi$. Recalling that $p_k+q_k+r_k = 1$ we write

$$ (p_k + r_k + q_k) \pi(k) = p_k \pi(k+1) + r_k \pi(k) + q_k \pi(k-1). \tag{12} $$

Introducing the function (defined for $0 \le k \le n-1$)

$$ \eta(k) := \pi(k+1) - \pi(k), \tag{13} $$

and writing for ease of notation

$$ \xi_k := q_k / p_k, \tag{14} $$

eq. (12) is promptly rearranged as the recursion relation

$$ \eta(k) = \xi_k \eta(k-1). \tag{15} $$

Thus, starting from $\eta_0 = \eta(0) = \pi(1) - \pi(0)$, we have

$$ \eta(k) = \left( \prod_{h=1}^{k} \xi_h \right) \eta_0. \tag{16} $$

Now we note that, due to (13),

$$ \pi(k) = \pi(0) + \sum_{h=0}^{k-1} \eta(h); \tag{17} $$

recalling also that $\pi(0) = 0$, specializing to $k=n$, and using (16), we get

$$ \pi(n) = \sum_{k=0}^{n-1} \eta(k) = \sum_{k=0}^{n-1} \left( \prod_{h=1}^{k} \xi_h \right) \eta_0. \tag{18} $$
Using now $\pi(n) = 1$, we get in the end

$$ \eta_0 = \left[ \sum_{k=0}^{n-1} \left( \prod_{h=1}^{k} \xi_h \right) \right]^{-1} . \qquad (19) $$

This yields immediately

$$ \eta(k) = \left( \frac{\prod_{h=1}^{k} \xi_h}{\sum_{k=0}^{n-1} \left( \prod_{h=1}^{k} \xi_h \right)} \right), \qquad (20) $$

and hence

$$ \pi(k) = \left( \frac{\sum_{h=0}^{k-1} \left( \prod_{m=1}^{h} \xi_m \right)}{\sum_{h=0}^{n-1} \left( \prod_{m=1}^{h} \xi_m \right)} \right). \qquad (21) $$

In order to evaluate $\chi(k)$, we could proceed along the same lines; but we can as well – and more economically – use a symmetry argument. Mapping $(k; p_k, q_k) \to (n-k, q_k, p_k)$ we have just reversed the role of the A and B alleles. This shows at once that we will get

$$ \chi(k) = \left( \frac{\sum_{h=k}^{n-1} \left( \prod_{m=1}^{h} \xi_m \right)}{\sum_{h=0}^{n-1} \left( \prod_{m=1}^{h} \xi_m \right)} \right). \qquad (22) $$

Note this implies immediately that $\pi(k) + \chi(k) = 1$ for all $k$.

**Remark 5.** Albeit the probabilities $\pi(k)$ and $\chi(k)$ can sometimes – and in particular for the model discussed here – be expressed in terms of hypergeometric functions, in concrete computations it may be convenient to fix $n$ and $\mu$, and use this scheme to provide exact values of the quantities of interest through a symbolic manipulation program rather than going through complex and rather obscure analytic expressions. $\odot$

**Example 1.** For the standard Moran model, see (2), and more generally for any model in which $p_k = q_k$, we have $\xi_k = 1$ for all $k \neq 0, n$; then (20) and (21) grant that $\pi(k)$ grows linearly from $\pi(0) = 0$ to $\pi(n) = 1$, each increment $\eta(k)$ being of course given by $\eta(k) = \eta_0 = 1/n$.

# 5 Hitting time: a general discussion

We now come to consider the mean hitting time $T(k)$, i.e. the time (measured in time unit such that there is one step per unit) needed to reach either one of the end states $k = 0$ or $k = n$. Note that, obviously, $T(0) = 0 = T(n)$. In our discussion we consider $n$ as fixed.

We denote the average time spent in state $j$ before hitting the border by a random walk starting in site $k$ as $\tau(j|k)$ (this is defined for $j \neq 0, n$); by definition,

$$ \tau(j|0) = 0 = \tau(j,n). \qquad (23) $$
Obviously we have

$$T(k) = \tau(1|k) + \dots + \tau(n-1|k) = \sum_{j=1}^{n-1} \tau(j|k). \quad (24)$$

Now we have to evaluate $\tau(j|k)$; we will again (as in the discussion of hitting probability) filter the random walk w.r.t. its first step, which takes the system from state $k$ to state $k+1$ with probability $p_k$, to state $k-1$ with probability $q_k$, and leaves it in state $k$ with probability $r_k = 1 - p_k - q_k$. Thus – noticing that if the walk starts from $k=j$ by the time we make a first step it will have spent a unit time in site $j$ – we have

$$\tau(j|k) \left( \begin{aligned} & p_k + r_k + q_k = \delta_{jk} + p_k \tau(j|k+1) \\ & + r_k \tau(j|k) + q_k \tau(j|k-1) \end{aligned} \right). \quad (25)$$

We write as above $\xi_j = q_j/p_j$ and also, for ease of notation,

$$\rho_j := 1/p_j. \quad (26)$$

Moreover we introduce

$$\theta(j, k) := \tau(j|k+1) - \tau(j|k). \quad (27)$$

Then the relation (25) reads

$$\theta(j, k) = \xi_k \theta(j, k-1) - \rho_j \delta_{jk}. \quad (28)$$

Thus the $\theta(j, k)$ can be determined recursively. We will write

$$\beta_j := \tau(j|1) - \tau(j|0) = \theta(j, 0). \quad (29)$$

With this, we easily get

$$\theta(j, k) = \begin{cases} \beta_j \prod_{i=1}^k \xi_i - \rho_j \prod_{i=j+1}^k \xi_i & \text{for } j < k, \\ \beta_j \prod_{i=1}^k \xi_i - \rho_j & \text{for } j = k, \\ \beta_j \prod_{i=1}^k \xi_i & \text{for } j > k. \end{cases} \quad (30)$$

Note that here the $\rho_j = 1/p_j$ are given, while $\beta_j$ are constants yet to be determined.

We next observe that, by (27) and recalling (23),

$$\tau(j|k) = \theta(j, 0) + \dots + \theta(j, k-1). \quad (31)$$

Recall that, again by (23), we have $\tau(j|n) = 0$ as well; thus specializing (31) to $k=n$ we get

$$\sum_{\ell=0}^{n-1} \theta(j, \ell) = 0. \quad (32)$$
This allows to determine $\beta_j$, as we now discuss.

Using (31), the relation (28) is promptly converted into an equivalent one for the $\tau$. Using $\tau(j|0) = 0$ and (31), we get

$$ \tau(j|k) = \begin{cases} \left[1 + \sum_{m=1}^{k-1} \left(\prod_{i=1}^{m} \xi_i\right)\right] \beta_j & (j < k), \\ \left[1 + \sum_{m=j+1}^{k-1} \left(\prod_{i=j+1}^{m} \xi_i\right)\right] \rho_j & (j = k), \\ \left[1 + \sum_{m=1}^{k-1} \left(\prod_{i=1}^{m} \xi_i\right)\right] \beta_j & (j > k). \end{cases} \quad (33) $$

In order to determine $\beta_j$, we should once again recall (23). Specializing (33) to the case $k=n$ (which is always in the class $j < k$, as $j=1, \dots, n-1$) and requiring this to vanish finally determines $\beta_j$ as

$$ \beta_j = \left[ \frac{1 + \sum_{m=j+1}^{n-1} \left( \prod_{i=j+1}^{m} \xi_i \right)}{1 + \sum_{m=1}^{n-1} \left( \prod_{i=1}^{m} \xi_i \right)} \right] \rho_j . \quad (34) $$

Now, inserting this in (33), one is able to determine explicitly – modulo the problem of evaluating the sums and products appearing in our formulas – the $\tau(j|k)$; and hence, using (24), also the mean hitting time $T(k)$.

Needless to say, the evaluation of sums and products can be trivial or extremely difficult (or impossible), depending on the expression of $\xi_k$ for the model at hand.

**Remark 6.** These formulas are substantially simplified for symmetric random walks, i.e. for $\xi_k = 1$. In this case, in fact, we get

$$ \theta(j, k) = \begin{cases} \beta_j - \rho_j & \text{for } j \le k, \\ \beta_j & \text{for } j > k. \end{cases} \quad (35) $$

This in turn yields

$$ \tau(j|n) = n \beta_j - (n-j) \rho_j; \quad (36) $$

hence requiring this to vanish we just get

$$ \beta_j = \left( \frac{n-j}{n} \right) \rho_j . \quad (37) $$

In the end, it results

$$ T(k) = k \sum_{j=1}^{n-1} \beta_j - \sum_{j=1}^{k-1} (k-j) \rho_j . \quad (38) $$
Combining this with (37), we get

$$ T(k) = \sum_{j=1}^{k-1} \left[ k \left( \frac{n-j}{j} \right) - (k-j) \right] \rho_j + \sum_{j=k}^{n-1} \left[ k \left( \frac{n-j}{j} \right) \right] \rho_j . \quad (39) $$

⊙

**Example 2.** For the simple symmetric random walk in $[0, n]$, $p = p_k = q_k = q = p$ ($0 < p \le 1/2$), and hence all the $\xi_i$ appearing in the formula reduce to unity, while $\rho_j = 1/p_j = 1/p$. Then

$$ \beta_j = \left( \frac{n-j}{n} \right) \left( \frac{1}{p} \right), \quad \tau(j|k) = \left( \frac{j(n-k)}{np} \right), $$

and in conclusion

$$ T(k) = \frac{(n-k)(n-1)}{2p}. $$

**Example 3.** For the simple non-symmetric random walk in $[0, n]$, $p = p_k$, $q_k = q$ ($0 < p + q \le 1$), all the $\xi_i$ are equal to $\xi = q/p$. In this case we get

$$ \begin{aligned} \beta_j &= \frac{1}{p} \left( \frac{1 - \xi^{n-j}}{1 - \xi^n} \right), \\ \tau(j, k) &= \frac{1}{p} \left( \frac{1 - \xi^j}{1 - \xi^n} \right) \left( \frac{1 - \xi^{n-k}}{1 - \xi} \right) \xi^{k-j}. \end{aligned} $$

In conclusion

$$ T(k) = \frac{1}{p} \left( \frac{[n\xi^n(1-\xi) - \xi(1-\xi^n)](1-\xi^{k-n})}{(1-\xi)^2(1-\xi^n)} \right). $$

# 6 Hitting probability: the complete Moran model

We will now apply the general discussion of Section 4 to the complete Moran model discussed in Section 3, i.e. to the case $\xi_k$ is given by (8). We will proceed as mentioned in Remark 5, i.e. analytically with given values for $n$ and $\mu$, and consider the dependence of the emerging results both on the total population number $n$ and on the mutation rate $\mu$.

## 6.1 Dependence on the total population number

The dependence of the form of $\pi(k)$ on the total population number is quite substantial. This is shown in Figure 2 and Figure 3. In particular, it appears from these that for small $n$ the hitting probability grows basically as in the standard Moran model (i.e. linearly) from $\pi(0) = 0$ to $\pi(n) = 1$, while for larger $n$ the behavior becomes more and more different; in particular, in the limit of large $n$ the probability of first hitting each of the two pure states is substantially 1/2 unless the initial state is very near to $k = 0$ or to $k = n$.
<figure><img src="image_2.png" /><figcaption>Figure 2: Plots of $\pi(k)$ for $\mu = 0.01$ and different values of $n$. For small $n$ the values of $\pi(k)$ grow nearly linearly (as for $\mu = 0$), while for higher values of $n$ the function $\pi(k)$ gets flatter in the central region and steeper near the borders.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: Comparison of the behavior of $\pi(k)$ for different values of $n$ and the same value of $\mu$ (here $\mu = 0.01$ as in Figure 2). The interpolating functions $f_n(x)$ are defined for the values of $\pi(k)$ reported in Figure 2; we plot $f_n(x*n)$ in order to compare the shape obtained for different values of $n$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4: The values of $\pi(k)$ for $n = 100$ and for different values of the mutation rate $\mu$: (a) $\mu = 0.001$ (black points); (b) $\mu = 0.01$ (blue points); (c) $\mu = 0.1$ (red points). For small $\mu$ the behavior is very similar to the linear one holding for the standard ($\mu = 0$) Moran model, while as $\mu$ is increased the behavior is more and more flattened towards a constant $1/2$ probability unless the process starts very near to the borders. [Note that the red line is perfectly horizontal in its core, but may appear tilted due to an optical effect related to the nearby diagonal lines.]</figcaption></figure>

## 6.2 Dependence on the mutation rate

From the biological point of view, the most interesting question is the dependence on the mutation rate $\mu$. This is shown in Figure 4. This shows a transition between the behavior of the standard Moran model, with the first hitting probability of an allele growing linearly with its initial weight in the population, to a situation in which the probability is essentially the same for the two alleles unless the initial state is extremely near to the borders. In other words, now – unless the initial state is extremely strongly skewed towards one of the alleles – first hitting one or the other of the alleles is a random process with equal probabilities.

This is easily understood in terms of our previous discussion: in this case the drift will lead the process towards equilibrium ($k = n/2$), so hitting the pure states will occur as the result of a large fluctuations from equilibrium, which can of course take place in both directions.

This observation also suggest that *the mean hitting time will substantially increase in the presence of mutations* (unless the mutation rate is negligible); this point will be discussed in the next Section.
<figure><img src="image_5.png" /><figcaption>Figure 5: The value of $f(x) = \pi(xn)$ for $\alpha \in [0, 1]$ for $\sigma = n\mu = 1$ obtained in two different ways: (a) $n = 1000$ and $\mu = 0.001$ (blue points); (b) $n = 100$ and $\mu = 0.01$ (red points). The functions obtained are similar but not equal. See the discussion in the text.</figcaption></figure>

## 6.3 The leading control parameter

We have considered the dependence of the hitting probability $\pi(k)$ (and intrinsically also of $\chi(k) = 1 - \pi(k)$) on the two parameters appearing in the complete Moran model; that is, the population number $n$ and the mutation rate $\mu$.

Actually, one can consider a combination of these, i.e.

$$ \sigma := n \mu, \qquad (40) $$

which plays the role of leading control parameter. By this we mean that the behavior of the model depends to a large extent on $\sigma$, independently on how this is obtained in terms of $n$ and $\mu$. For small $\sigma$ the behavior of the model is similar to that of the standard Moran model, while for large $\sigma$ the behavior is different, showing the flattening of the $\pi(k)$ function in the core and its steepening at the borders.

Equivalently, we can say that the mutation rate is small or large by comparing it with $\mu_0 = 1/n$. We stress that a mutation rate $\mu_0$ means that there is one mutation once each individual has reproduced once (on the average); thus it is a very natural measure unit for the mutation rate.

Note that there is no simple scaling relation in the complete Moran model, as appears from the functional form of $F(x)$ defined in (9). This is also shown in Figure 5, where we plot $f(x) := \pi(xn)$ in order to compare the situation for different values of $n$ and $\mu$ yielding the same value for $\sigma$.
<figure><img src="image_6.png" /><figcaption>Figure 6: The mean hitting time $T(k)$ for the complete Moran model with $\mu = 0.01$ and different values of $n$. Together with an increase of the hitting time, there is a variation in the shape of $T(k)$, see also Figure 7.</figcaption></figure>

# 7 Hitting time: the complete Moran model

We can now apply our general formalism developed in Section 5 to the complete Moran model, i.e. to the case where (5) and (6) hold. We proceed in the same way as in our general discussion above, and will work with given population number $n$ and mutation rate $\mu$, which allow to determine exactly all the quantities of interest.

## 7.1 Dependence on the total population

We will first consider the situation where the mutation rate is fixed, and look at how the mean hitting time varies with the population number. The result of computations for $\mu = 0.01$ and different population numbers are shown in Figure 6.

Needless to say, an increase in $n$ leads to larger hitting times. Together with this, we observe that the shape of the function $T(k)$ also changes. In order to compare the results for different $n$, we consider interpolating functions $\Theta_n(k)$ for $T(k)$ with total population $n$, and normalize these so that their maximum is one; in other words, we consider

$$ \vartheta_n(x) = \Theta_n(x * n)/\Theta_n(n/2), \quad (41) $$

defined for $x \in [0, 1]$. The different functions $\vartheta_n$ corresponding to the data plotted in Figure 6 are compared in Figure 7. This confirms that when $n$ is
<figure><img src="image_7.png" /><figcaption>Figure 7: The functions $\vartheta_n(x)$ defined in eq.(41) for the different values of $n$ considered in Figure 6: $n = 10$ (black curve), $n = 25$ (blue), $n = 50$ (orange), $n = 100$ (red). For higher $n$, the function shape becomes flatter in the core and steeper near the borders.</figcaption></figure>

increased, the function $T(k)$ gets a flatter form in the core, and a steeper one near the borders.

## 7.2 Dependence on the mutation rate

The main motivation for this work was to consider the effect of introducing mutations in the standard Moran model; this is considered as a paradigmatic model for neutral evolution, and also used to get a rough estimate of the timescales for allele fixation.

Albeit in the complete Moran model fixation is not possible, and is substituted by reaching (temporarily) a “pure” state, the dependence of the mean hitting time on the mutation rate is specially interesting from the point of view of Evolutionary Dynamics (see the discussion in Sect. ?? below).

Our discussion of the hitting probabilities for the two alleles suggest that on the one hand the presence of mutations will lead to an increase of the mean hitting time; on the other hand we have seen that when $\mu$ increases the hitting probabilities are essentially independent of the initial state (and equal for the two alleles), unless the initial state is very near to the borders. This suggests that the process spends a long time wandering back and forth in the core of the state space, until a large fluctuation leads to fixation. This picture suggests that the shape of the mean hitting time function $T(k)$ will also be changed as $\mu$ is increased, and in particular that its dependence on the initial state $k$ will be not so relevant except when $k$ is near to the borders.
<figure><img src="image_8.png" /><figcaption>Figure 8: Dependence of mean hitting time $T(k)$ on the mutation rate for constant population size. Here we consider a population of $n = 100$ individuals, and plot the function $T(k)$ for different values of the mutation rate $\mu$; in detail: for $\mu = 0$ (black points), corresponding to the standard Moran model; for $\mu = 0.001$ (red points); for $\mu = 0.005$ (green points); and for $\mu = 0.01$ (blue points). The increase in $\mu$ leads to a change is the shape of the curve, which becomes wider, and above all to a raising in the curve, i.e. to an increase in the hitting time.</figcaption></figure>

The exact evaluation of $T(k)$ for a population of fixed size ($n = 100$) and different values of $\mu$, whose results are reported in Figure 8, confirms these expectations.

## 7.3 Reversal time in the presence of mutations

In the absence of mutations, once the process reaches either one of the pure states the situation is fixed forever. If mutations are possible (but staying within the simple frame in which only two alleles can exist) the allele which at some point is absent from the population will appear again due to a mutation; and the process will at some point, after a time $T$ from hitting the pure state, reach the opposite pure state.

The typical time scale for this can be estimated, at first order, as

$$ \mathcal{T} = [\mu \pi(1)]^{-1} T(1), \qquad (42) $$

where, as before, $\pi(k)$ is the probability to fix $A$ starting from the state $k_0 = k$ and $T(k)$ is the mean time needed to fixation starting from state $k_0 = k$.

A comparison between $\mathcal{T}$ and the maximal mean hitting time $T_M$ (which
<table><caption>Table 1: Comparison of reversal time $\mathcal{T}$ and of the maximal mean hitting time $T_M = T(n/2)$ in the complete Moran model for $\mu = 0.01$ and different values of the total population $n$. In all cases, $\mathcal{T}$ is two orders of magnitude larger than $T_M$.</caption><thead><tr><th>n</th><th>π(1)</th><th>T(1)</th><th>T</th><th>T<sub>M</sub></th></tr></thead><tbody><tr><td>10</td><td>0.11</td><td>30.63</td><td>27845</td><td>68</td></tr><tr><td>25</td><td>0.06</td><td>131.85</td><td>219750</td><td>496</td></tr><tr><td>50</td><td>0.06</td><td>548.51</td><td>914183</td><td>2558</td></tr><tr><td>100</td><td>0.10</td><td>5290.58</td><td>5290580</td><td>19977</td></tr></tbody></table>

corresponds to $T(n/2))$ is provided in Table 1. It appears from this that $\mathcal{T}/T_M \approx 10^2$ in all considered cases.

**Remark 7.** As stressed in the Introduction, the complete Moran model is an ergodic Markov chain, admitting a single equilibrium state; in particular, no fixation is possible. One can however consider, from the mathematical point of view, a related model with absorbing boundary condition. This corresponds to stopping the process once a pure state is reached. In biological terms, this would correspond to the – extremely unrealistic – hypothesis that mutations are not possible if the system is in a pure state, and in this sense such a model should be seen as a mathematical one, possibly with real-world applications but surely not relevant in the context we have been discussing and which provides motivation for the present work. It is obvious that all the computations presented above go unaltered in such a case, albeit with a different meaning: hitting probabilities become fixation probabilities, and hitting time becomes fixation time. $\circ$

# 8 Discussion and conclusions

We have considered a modification of the standard Moran model for neutral evolution, actually already considered in the original paper by Moran [1]. This model takes into account mutational as well as reproductive fluctuations, and is hence called the “complete” Moran model.

The presence of mutations forbids that an allele is fixed forever, so the questions to be investigated concern (a) the probability $\pi(k)$ to first hit one of the pure states or the other if initially there are $k$ individuals with the first allele in a population of $n$ individuals, and (b) the time $T(k)$, under these initial conditions, for a pure state to be first hit. Note that the hit condition corresponds, in the standard (mutation-less) Moran model, to fixation of an allele.

As for point (a), we have observed that the main difference is that while in the standard Moran model $\pi(k)$ grows linearly from $\pi(0) = 0$ to $\pi(1) = 1$, in the complete model the behavior is nonlinear. More precisely, there is a sharp increase for $k \simeq 0$ followed by a substantially flat behavior and by a new sharp increase for $k \simeq n$. In particular this means that unless the initial conditions are strongly skewed towards one of the pure states, each of them will have
a probability 1/2 to be selected for first hit by fluctuations. This effect gets quickly more important with the growth of the mutation rate $\mu$ for constant $n$, and also increases with the growth of the population size $n$ for constant $\mu$.

As for point (b), we have observed two effects. On the one hand, the presence of mutations slows down the hitting of pure states, i.e. for a given population size $n$ and initial condition $k$, the corresponding mean hitting time $T(k)$ increases (as obvious) with increasing $\mu$. On the other hand, the shape of the function $T(k)$ is also modified as $\mu$ increases; here too, for larger $\mu$ the function gets steeper near the borders $k=0$ and $k=n$ and flatter in the core.

As for the methods employed, we have used exact recursive relations to evaluate $\pi(k)$ and $T(k)$ at given $n$. General formulas providing these quantities in terms of special functions (in particular, hypergeometric ones) can be obtained, and are given in the Appendix, but the recursion formulas appear to be more manageable for our purposes.

# A Explicit expressions for the complete Moran model

As mentioned above, the hitting probability and mean hitting time for the complete Moran model can be explicitly computed; the resulting formulas are a bit complex and involve special functions, in particular hypergeometric functions and the digamma function [34].

## A.1 Hitting probability

Recalling that the general expression for $\pi(k)$ is given by eq.(21), where $\xi_m := q_m/p_m$, see eq. (14), and that for the complete Moran model $p_m$ and $q_m$ are given by eqs. (5) and (6) respectively, we obtain – by a direct computation with Mathematica – that in this case the function $\pi(k)$ can be expressed in terms of the Hypergeometric function

$$ {}_2F_1(a, b; c; z) := \sum_{k=0}^{\infty} \left( \frac{P(a, k) P(b, k)}{P(c, k)} \right) \frac{z^k}{k!}, $$

where $P(a, k)$ is the function$^2$

$$ P(a, k) := \frac{\Gamma(a + k)}{\Gamma(a)} $$

with $\Gamma$ the Euler Gamma function.

In order to have a lighter notation, we will set

$$ \mathcal{F}(b; c; z) := {}_2F_1(1, b; c; z) $$

<footnote>$^2$This is also known in the literature as Pochhammer symbol; a more convenient name for it is the *rising factorial*, and it is associated to Stirling numbers of the first kind [35].</footnote>
and moreover write

$$ \nu = 1 - n; \quad \eta = \frac{\mu - 1}{\mu}. $$

We then have

$$ \pi(k) = \frac{\Gamma(k + \nu) \mathcal{F}(1; \nu; \eta) - \eta^k \Gamma(k + 1) \Gamma(\nu) \mathcal{F}(k + 1; k + \nu; \eta)}{\Gamma(k + \nu) [\mathcal{F}(1; \nu; \eta) - \mu(\eta\mu)^n n\pi \csc(n\pi)]} . \quad (43) $$

## A.2 Mean hitting time

Recalling that the general expression for $T(k)$ is given by eq.(39), where $\rho_j$ is given by eq.(26), and that for the complete Moran model $p_j$ is given by eq.(5), we obtain – again by a direct computation with Mathematica – an explicit formula. In order to simplify notation, we will denote by $\Psi(z)$ the digamma function – which is often denoted as $\psi^{(0)}(z)$ – that is,

$$ \Psi(z) = \psi^{(0)}(z) = \frac{\Gamma'(z)}{\Gamma(z)} $$

with $\Gamma(z)$ the usual Euler Gamma function.

Moreover, we denote by $\gamma_e$ the Euler-Mascheroni gamma constant,

$$ \gamma_e \approx 0.577216, $$

and write for short

$$ \lambda = \frac{\mu}{1 - 2\mu}, \quad \sigma = 1 + \lambda = \frac{1 - \mu}{1 - 2\mu}. $$

With these notations, it results

$$ T(k) = \left( \frac{n}{\mu(1-\mu)(1-2\mu)} \right) \times \left[ \begin{aligned} & \gamma_e (1-\mu)(1-2\mu)k \\ & + \mu(1-2\mu)(n-k)[\Psi(1-n) - \Psi(k-n)] \\ & + k(1-\mu)(1-2\mu)[\Psi(n) - \Psi(n\sigma)] \\ & + (n\mu^2 - 2k\mu)[\Psi(1+n\lambda) - \Psi(k+n\lambda)] \\ & + k[\Psi(1+n\lambda) - \mu\Psi(k+n\lambda)] \end{aligned} \right] . \quad (44) $$

This might be somewhat simplified using the relations

$$ \Psi(1+z) = \frac{1}{z} + \Psi(z), $$

$$ \Psi(1-z) = \Psi(z) + \pi \cot(\pi x); $$

and its asymptotic behavior for large $n$ can be inferred from

$$ \Psi(z) \approx \log(z) - \frac{1}{2z} $$

which holds for $|z| \to \infty$ for $|\arg(z)| < \pi - \varepsilon$.
### A.3 Mean hitting time for small mutation rate

Finally, let us consider the series expansion for small mutation rate $\mu$. Reverting to the usual polygamma notation [34], i.e. with

$$ \psi^{(k)}(z) := \frac{d^k \Psi(z)}{dz^k} = \frac{d^{k+1}}{dz^{k+1}} \log[\Gamma(z)], $$

the series expansion for small $\mu$ yields

$$ T(k) \approx \left[ \frac{\pi^2}{6} kn^2 - \gamma_e kn + n^2 (\psi^{(0)}(1-n) - \psi^{(0)}(k-n) - k\psi^{(1)}(n)) \right. \\ \left. - kn(\psi^{(0)}(k) + \psi^{(0)}(1-n) - \psi^{(0)}(k-n)) \right] \\ + \left[ \frac{\pi^2}{2} kn^2 - \gamma_e n(n+k) \right. \\ + \frac{1}{2} kn^3 (\psi^{(2)}(1) - \psi^{(2)}(n)) - kn^2 (\psi^{(1)}(k) - 2\psi^{(1)}(n)) \\ + n^2 (\psi^{(0)}(1-n) - \psi^{(0)}(k) - \psi^{(0)}(k-n)) \\ \left. + kn (\psi^{(0)}(k-n) - \psi^{(0)}(k) - \psi^{(0)}(1-n)) \right] \mu \\ + O(\mu^2). \qquad (45) $$

## Acknowledgements

This work was started while I was on sabbatical leave from Università degli Studi di Milano, enjoying the hospitality and the relaxed work atmosphere of SMRI; in particular I thank L. Peliti for sharing with me his notes on Knuth’s “Two notes on notation” [35].

This work was supported by “iniziativa specifica MMNLP” of INFN; I am also a member of IAMP and of ISNMP. My work is supported by GNFM-INdAM.

## References

[1] P.A.P. Moran, “Random processes in genetics”, *Math. Proc. Cambridge Phil. Soc.* **54** (1958), 60-71

[2] S. Karlin and J. McGregor, “On a genetics model of Moran”, *Math. Proc. Cambridge Phil. Soc.* **58** (1962), 299-311

[3] P. Donnelly, “The transient behaviour of the Moran model in population genetics”, *Math. Proc. Cambridge Phil. Soc.* **95** (1984), 349-358

[4] M. Kimura, “Diffusion Models in Population Genetics”, *J. Appl. Prob.* **1** (1964), 177-232
[5] M. Kimura, "Evolutionary rate at the molecular level", *Nature* **217** (1968), 624-626

[6] M. Kimura, "DNA and the neutral theory", *Phil. Trans. Royal Soc. London B* **312** (1986), 343-354

[7] M. Kimura, "Molecular evolutionary clock and the neutral theory", *J. Molec. Evol.* **26** (1987), 24-33

[8] M. Kimura, "The neutral theory of molecular evolution: A review of recent evidence", *Jap. J. Genetics* **66** (1991), 367-386

[9] M. Kimura, *Population genetics, molecular evolution, and the neutral theory: selected papers*, University of Chicago Press 1994

[10] T. Ohta, "The nearly neutral theory of molecular evolution", *Ann. Rev. Ecology and Systematics* **23** (1992), 263-286

[11] T. Ohta and J.H. Gillespie, "Development of neutral and nearly neutral theories", *Theor. Popul. Biol.* **49** (1996), 128-142

[12] E.G. Leigh, "Neutral theory: a historical perspective", *J. Evol. Biol.* **20** (2007), 2075-2091

[13] E. Baake and R. Bialowons, "Ancestral processes with selection: Branching and Moran models", *Banach Center Publications* **80** (2008), 33-52 (arXiv:q-bio/0702002)

[14] F.A. Chalub and M.O. Souza, "On the stochastic evolution of finite populations", *J. Math. Biol.* **75** (2017), 1735-1774

[15] F.A. Chalub and M.O. Souza, "Fitness potentials and qualitative properties of the Wright-Fisher dynamics", *J. Theor. Biol.* **457** (2018), 57-65

[16] C.A. Muirhead and J. Wakeley, "Modeling Multiallelic Selection Using a Moran Model" *Genetics* **182** (2009), 1141-1157

[17] M.A. Nowak, C.E. Tarnita and T. Antal, "Evolutionary dynamics in structured populations", *Phil. Trans. R. Soc. B* **365** (2010), 19-30

[18] A. Traulsen, J.C. Claussen, and C. Hauert, "Coevolutionary Dynamics: From Finite to Infinite Populations", *Phys. Rev. Lett.* **95** (2005), 238701

[19] A. Traulsen, A.M. Sengupta and M.A. Nowak, "Stochastic evolutionary dynamics on two levels", *J. Theor. Biol.* **235** (2005), 393-401

[20] A. Traulsen, M.A. Nowak and J.M. Pacheco, "Stochastic dynamics of invasion and fixation" *Phys. Rev. E* **74** (2006), 011909

[21] C. Vogl and L. Caitlin Mikula, "A nearly-neutral biallelic Moran model with biased mutation and linear and quadratic selection", *Theor. Popul. Biol.* **139** (2021), 1-17
[22] J. Dalmau, “Convergence of a Moran model to Eigen’s quasispecies model”, *J Theor. Biol.* **420** (2017), 36-40

[23] R. Cerf and J. Dalmau, “The distribution of the quasispecies for a Moran model on the sharp peak landscape”, *Stoch. Proc. Appl.* **126** (2016), 1681-1709

[24] R. Cerf, “Critical population and error threshold on the sharp peak landscape for a Moran model”, *Ann. Appl. Prob.* **25** (2015), 1936-1992

[25] R. Cerf, “Critical population and error threshold on the sharp peak landscape for the Wright-Fisher model”, *Ann. Appl. Probab.* **25** (2015), 1936-1992

[26] J. Dalmau, “The Wright-Fisher model for class-dependent fitness landscapes”, *El. J. Prob.* **26** (2021), 1-44

[27] R. Cerf and J. Dalmau, “Quasispecies on class-dependent fitness landscapes”, *Bull. Math. Biol.* **78** (2016), 1238-1258

[28] R. Cerf and J. Dalmau, *The quasispecies equation and classical population models*, Springer 2022

[29] D.B. Saakian, M.W. Deem and C.K. Hu, “Finite population size effects in quasispecies models with single-peak fitness landscape”, *EuroPhys. Lett.* **98** (2012), 18001

[30] D.B. Saakian, T. Yakushkina and E.V. Koonin, “Allele fixation probability in a Moran model with fluctuating fitness landscapes”, *Phys. Rev. E* **99** (2019), 022407

[31] H.C. Berg, *Random Walks in Biology*, Princeton University Press 1993

[32] E.A. Codling, M.J. Plank and S. Benhamou, “Random walk models in Biology”, *J. Royal Soc. Interface* **5** (2008), 813-834

[33] Z. Patwa and L. M. Wahl, “The fixation probability of beneficial mutations”, *J. R. Soc. Interface* **5** (2008), 1279-1289

[34] M. Abramowitz and I.A. Stegun eds., *Handbook of Mathematical Functions*, National Bureau of Standards, 1972

[35] D.E. Knuth, “Two notes on notation”, *Amer. Math. Monthly* **99** (1992), 403-422