# Fixation Probabilities on Complete Star and Bipartite Digraphs

Authors: Chunming Zhang; Yingjiang Wu; Wanping Liu; Xiaofan Yang
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

This paper exactly formulates thekth-order fixation probabilities on complete star digraphs (CSDs), which extend the results from Broom and Rychtář (2008). By applying these probability formulae, some asymptotic properties on CBDs are analyzed, and certain CSDs are determined to be amplifiers of selection for arbitrary relative fitness larger than 1, while all the CSDs are proved to be amplifiers of selection for fixed relative fitness slightly larger than 1. A numerical method for fixed population structure (by solving a linear system) is developed to calculate the fixation probabilities on complete bipartite digraphs (CBDs), and some conjectures are finally given through simulations.


---

## Body

## 1. Introduction

As a newly emerging branch of population dynamics, evolutionary graph theory studies the evolution of structured populations and explores the effect of population structure on evolutionary dynamics [1–6]. Primarily, evolutionary dynamics were investigated on homogeneous infinite populations however, populations in the real world are neither infinite nor homogeneously mixed [2]. In recent years, there is a growing interest in the investigation of evolutionary dynamics on spatial structures or populations with certain nonhomogeneous structure [5, 7–10]. The study of evolutionary dynamics on graphs was popularized by Lieberman et al. [4]. In this framework, the structure of a population is modeled by a weighted digraph on vertices 1,2,…,N, which are occupied by individuals (residents and mutants). The graph can describe the architecture of cells in a multicellular organism and represent spatial structure among animals or plants in an ecosystem [8]. Moreover, graphs can also represent relationships in a social network of humans, which means that the dynamics describes cultural selection and the spread of new inventions and ideas. It is supposed that every individual of the population occupies a unique vertex of the graph. In each iteration, a random individual i is chosen for reproduction with a probability proportional to its fitness, and the resulting offspring will occupy an adjacent vertex j with the probability wij, which represents the weight on edge 〈i,j〉, if any. The intrinsic weights of a digraph are defined this way: for a vertex i with l outgoing edges, let wij=1/l if edge {i,j} exists. In this paper, we focus our attention on the structures of this kind of digraphs with intrinsic weights.Consider a homogeneous population on a weighted graph whose individuals all have baseline fitness 1. Each individual is chosen as the reproducing one with probability proportional to its fitness. Supposek new mutants with relative fitness r (the new mutant called advantageous if r>1, while disadvantageous if 0<r<1) are introduced by placing them on k randomly chosen vertices of the graph. These mutants have a certain chance of fixation, that is, to generate a lineage that takes over the population. It is an issue in population dynamics to find the kth-order fixation probability, that is, the fixation probability of k mutants, in a population [7, 11, 12].An unstructured population can be modeled by a complete digraph where all edges have the same weight. The evolution of an unstructured population is often modeled by theMoran process, whose kth-order fixation probability is
(1.1)ρk(M)(N;r)=(1-1/rk)(1-1/rN),
where resident individuals have fitness 1, while mutants have relative fitness r. The Moran process defines a balance between natural selection and random drift. Lieberman et al. [4] gave the approximation of the first-order fixation probability for stars with large size N by
(1.2)P=(1-1/r2)(1-1/r2N).
Broom and Rychtář [7] obtained the exact average fixation probability for a mutant, given by ϱ=[n2r/(nr+1)+r/(n+r)]/[(n+1)[1+n/(n+r)∑j=1n-1‍((n+r)/(r(nr+1)))j]], which recovers (1.2) because for large N we derive ϱ≈1/(1+∑j=1n-1‍1/r2j)=(1-1/r2)/(1-1/r2n), where n represents the total size of a star except the central individual. In this paper, we further extend these results by deriving the explicit kth-order fixation probability on complete star digraphs (there exists a center vertex which connects each other vertex by two directed edges) as follows:
(1.3)ϱ=n2r/(nr+1)+r/(n+r)(n+1)[1+n/(n+r)∑j=1n-1((n+r)/r(nr+1))j].Thetemperature of a vertex is defined as the sum of all weights that lead into that vertex, that is, the temperature of vertex j is given by Tj=∑j=1N‍wij[8]. If all the vertices have the same temperature, then a graph is isothermal. It is known that a structured population has the same first-order fixation probability as the corresponding Moran process if and only if the structure is an isothermal digraph [8]. A structure is referred to as an amplifier of selection (resp., suppressor of selection) if the first-order fixation probability of one advantageous mutant on this structure is greater than (resp., less than) that for the corresponding Moran process. An important issue in evolutionary graph theory is to answer whether a given structure is an amplifier of selection or a suppressor of selection [4, 13, 14].This paper proceeds as follows. In Section2, we formulate the exact kth-order fixation probability on complete star graphs, then apply it to answer whether a given complete star graph is an amplifier of selection or not. In Section 3, we obtain the first-order fixation probability on complete bipartite digraphs (all vertices are divided into two sets, and each vertex is connected to every vertex in the other set by a directed edge) by solving a linear system through numerical methods, and pose several conjectures based on the simulations results.In this paper, we further develop the theory of evolutionary processes on graphs first developed by Lieberman et al. [4], by following the approach of proving analytical results for simple systems in a similar way to Broom and Rychtář [7]. Some useful results advance the formal underpinnings of the modeling of evolution on graphs, since nearly no attempt is made to prove general theoretical results in studies of evolution on graphs.


## 2. On Complete Star Digraphs

In this section, we obtain explicit formulae for thekth-order fixation probabilities on complete star digraphs and explore their properties and applications in nature and society.


### 2.1. Basic Concepts and Notations

Acomplete star digraph (CSD) is a digraph with a single central vertex such that (I) there exists an edge from the center vertex to each peripheral vertex; (II) there exists an edge from each peripheral vertex to the center vertex; (III) there exist no other edges. Therefore, a CSD of size N, denoted by 𝒮N, is a digraph with vertex set 𝒱={1,2,…,N} and directed edge set ℰ={〈1,j〉,〈j,1〉∣j=2,…,N}. Figure 1 illustrates one CSD with intrinsic weights. Note that the term “CSD” means “CSD with intrinsic weights” in the sequel.Figure 1
One complete star digraph𝒮6 with intrinsic weights, the central vertex is indexed 1.Thekth-order fixation probability on 𝒮N is denoted by ρk(N;r), which represents the probability of the event that these k mutants generate a lineage that takes over the whole population.For technical reasons, we need the following notations in the sequel.The configuration of a population on𝒮N, at time t, is depicted by a vector Θ(t)=(θ1(t),θ2(t)), where θ1(t)=1 or 0 if a mutant occupies the central vertex or not, respectively, and θ2(t) denotes the number of mutants occupying vertices 2,…,N. The total number of mutants at time t is denoted by θ(t), that is, θ(t)=θ1(t)+θ2(t). The probability of the event that the mutants finally fixate by starting with Θ(0)=(k1,k2), is denoted by ρk1,k2(N;r) abbreviated to ρk1,k2 (without ambiguity).


### 2.2. Explicit Fixation Probability Formulae

Before exactly formulating thekth-order fixation probability, the following two necessary lemmas are presented.Lemma 2.1.
Consider a CSD of sizeN≥3 and denote by r the relative fitness of a mutant. Then for k=1,…,N-1, the difference equation system holds:
(2.1)pk=Apk-1
where pk=(ρ1,k,ρ0,k)T∈ℝ2, and
(2.2)A=(α1-ααβ1-αβ),
with α:=1+(N-1)/r,β:=r(N-1)/[r(N-1)+1].Proof.
This proof proceeds by calculating some conditional probabilities and by the total probability formula. LetΘ(0), at time t=0, represent the initial configuration of a population on 𝒮N, Θ(1) represent the new configuration after one step-time, and k=θ(0)=θ1(0)+θ2(0) represent the initial total number of mutants. Then two cases are discussed with respect to Θ(0).
Case1.  (The initial configuration Θ(0)=(1,k-1)). In order to reduce the k mutants by one, a resident individual must be selected to reproduce and, meanwhile the center mutant (staying at vertex 1) has to be chosen for death, whose conditional on the resident’s selection, happens with probability 1. Therefore, the conditional probability of going from Θ(0)=(1,k-1) to Θ(1)=(0,k-1) is given by
(2.3)Pr{Θ(1)=(0,k-1)∣Θ(0)=(1,k-1)}=(N-k)(rk+N-k).
On the other hand, to increaseω by one, the center mutant has to be chosen for reproduction, and one resident individual connected to the center vertex must be selected to die. Thus the conditional probability of going from Θ(0)=(1,k-1) to Θ(1)=(1,k) is given by
(2.4)Pr{Θ(1)=(1,k)∣Θ(0)=(1,k-1)}=r(N-k)[(rk+N-k)(N-1)].
By employing (2.3) and (2.4), the probability of the event that the configuration stays unchanged is given by
(2.5)Pr{Θ(1)=(1,k-1)∣Θ(0)=(1,k-1)}=1-(N-k)(rk+N-k)-r(N-k)[(rk+N-k)(N-1)].
Case2. (The initial configuration Θ(0)=(0,k)). On the one hand, in order to increase k by one, the center individual has to be replaced by a new mutant produced by another. Therefore, the probability of going from Θ(0)=(0,k) to Θ(1)=(1,k) is given by
(2.6)Pr{Θ(1)=(1,k)∣Θ(0)=(0,k)}=rk(rk+N-k).
On the other hand, to reducek by one, the center resident individual has to be chosen for reproduction and its offspring must replace one mutant. Therefore, the probability of going from Θ(0)=(0,k) to Θ(1)=(0,k-1) is given by
(2.7)Pr{Θ(1)=(0,k-1)∣Θ(0)=(0,k)}=k[(rk+N-k)(N-1)].
It follows from (2.6) and (2.7) that the probability of the configuration remaining unchanged is
(2.8)Pr{Θ(1)=(0,k)∣Θ(0)=(0,k)}=1-rk(rk+N-k)-k[(rk+N-k)(N-1)].
By the total probability formula, we derive that for k=1,…,N-1,
(2.9)ρ1,k-1=Pr{Θ(1)=(0,k-1)∣Θ(0)=(1,k-1)}ρ0,k-1+Pr{Θ(1)=(1,k)∣Θ(0)=(1,k-1)}ρ1,k+Pr{Θ(1)=(1,k-1)∣Θ(0)=(1,k-1)}ρ1,k-1,(2.10)ρ0,k=Pr{Θ(1)=(1,k)∣Θ(0)=(0,k)}ρ1,k+Pr{Θ(1)=(0,k-1)∣Θ(0)=(0,k)}ρ0,k-1+Pr{Θ(1)=(0,k)∣Θ(0)=(0,k)}ρ0,k.
Plugging (2.3)–(2.8) into (2.9)-(2.10) and simplifying the results, we deduce
(2.11)ρ1,k=αρ1,k-1+(1-α)ρ0,k-1,(2.12)ρ0,k=βρ1,k+(1-β)ρ0,k-1,
where α:=1+(N-1)/r,β:=r(N-1)/[r(N-1)+1].
By substituting (2.11) into (2.12) and rearranging the terms, we get
(2.13)ρ0,k=αβρ1,k-1+(1-αβ)ρ0,k-1.
Equation on (2.1) follows by combining (2.11) and (2.13).Lemma 2.2.
LetN and r be given. Then for 1≤k≤N-1, the following hold
(2.14)ρ1,k-1=[β-(α-1)αk-2(1-β)k-1][β-(α-1)αN-2(1-β)N-1],ρ0,k=[β(1-αk(1-β)k)][β-(α-1)αN-2(1-β)N-1].Proof.
We solve the linear difference equation system (2.1) in view of the standard technique given in Elaydi [15]. The matrix A in (2.2) has two real eigenvalues λ1=1,λ2=α(1-β).
Furthermore,η1=(1,1)T and η2=(α-1,αβ)T are eigenvectors of A corresponding to λ1,λ2, respectively. Denote
(2.15)D=(1α-11αβ),
and its inverse matrix
(2.16)D-1=[1(1-α+αβ)](αβ1-α-11).
Through the method of diagonalizing matrix A, we make the variable change:
(2.17)(ρ1,kρ0,k)=D(xkyk),
and thus
(2.18)(xkyk)=D-1AD(xk-1yk-1)=(100α(1-β))(xk-1yk-1)=(100αk(1-β)k)(x0y0).
Therefore,
(2.19)(ρ1,kρ0,k)=D(100αk(1-β)k)D-1(ρ1,0ρ0,0)=ρ1,0(1-α+αβ)(αβ-(α-1)αk(1-β)kαβ-αk+1β(1-β)k).
Particularly, we know
(2.20)1=ρ1,N-1=[ρ1,0(1-α+αβ)]×[αβ-(α-1)αN-1(1-β)N-1],
which implies
(2.21)ρ1,0=(1-α+αβ)(αβ-(α-1)αN-1(1-β)N-1).
Equation on (2.14) follows from substituting (2.21) into (2.19). The proof is complete.Next, we present the main theorem of this section.Theorem 2.3.
Letr and N>2 be given. Then for each k=1,…,N-1, the explicit kth-order fixation probability on 𝒮N is given by
(2.22)ρk(N;r)=rN[r(N-1)+1]N-k-2[rN(r(N-1)+1)N-2-(N-1+r)N-2]×[((N-1+r)r)k(r(N-1)+1)k-(krk)(N-2)(r2-1)(N-1+r)k-2-((N-1+r)r)k].Proof.
In light of (2.21), we derive
(2.23)ρk(N;r)=kNρ1,k-1+(1-kN)ρ0,k=k/Nαk-2(1-β)k-1[α2β(1-β)-(α-1)]+β[1-αk(1-β)k]β-(α-1)αN-2(1-β)N-1.
By substitutingα and β (see Lemma 2.1) into this equation, we have
(2.24)ρk(N;r)={k/N[1+(N-1)/r]k-2[1/(r(N-1)+1)]k-1[(1+(N-1)/r)2(r(N-1))/((r(N-1)+1)2)-(N-1)/r]+r(N-1)/(r(N-1)+1)[1-(1+(N-1)/r)k/(r(N-1)+1)k]}/[(r(N-1)+1)N-1(r(N-1))/(r(N-1)+1)-(N-1)/r[1+(N-1)/r]N-2/(r(N-1)+1)N-1]=[rN-1(r(N-1)+1)N-1]{((r(N-1)+1)k+1)k(N-1)(r2-1)(2-N)(N+r-1)k-2/[rk-1(r(N-1)+1)k+1]+r(N-1)/((r(N-1)+1)k+1)[(r(N-1)+1)k-((r+N-1)/r)k]}/{(N-1)[rN(r(N-1)+1)N-2-(r+N-1)N-2]}=rN[r(N-1)+1]N-k-2/[rN(r(N-1)+1)N-2-(N-1+r)N-2]×[(r(N-1)+1)k-(k/rk)(N-2)(r2-1)(N-1+r)k-2-((N-1+r)/r)k].
The proof is complete.Particularly, for the casek=1, (2.22) becomes
(2.25)ρ1(N;r)=rN-1(r2-1)(r(N-1)+1)N-3[rN(r(N-1)+1)N-2-(N-1+r)N-2]×[(N2+(r-3)(N-1))(N-1+r)].Here we point out that (2.25) is just the same as the formula ϱ=[n2r/(nr+1)+r/(n+r)]/[(n+1)[1+n/(n+r)∑j=1n-1‍((n+r)/(r(nr+1)))j]], which is given in Broom and Rychtář [7], because
(2.26)ϱ=n2r/(nr+1)+r/(n+r)(n+1)[1+n/(n+r)∑j=1n-1((n+r)/(nr2+r))j]=n2r/(nr+1)+r/(n+r)(n+1)[1+[rn-1(nr+1)n-1-(n+r)n-1]/[(r2-1)(r(nr+1))n-1]]=[n(n+1)r2+(n+1)(n2-n+1)r](r2-1)[r(nr+1)]n-1(nr+1)(n+r)(n+1)[rn+1(nr+1)n-1-(n+r)n-1]=rN-1(r2-1)[r(N-1)+1]N-3rN[r(N-1)+1]N-2-[N-1+r]N-2×(N2+(r-3)(N-1))(N-1+r)=ρ1(N;r),
where n represents the number of all individuals except the central one, that is,  n=N-1.


### 2.3. Applications and Properties of the Fixation Probability

Theorem2.3 allows us to calculate the fixation probability for a given CSD, so here we use it to answer whether a given CSD is an amplifier of selection or not. We have the following result.Theorem 2.4.
𝒮
3
,
𝒮
4
,
𝒮
5
,
𝒮
6
,
𝒮
7
,
𝒮
8
,
𝒮
9, and 𝒮10 are all amplifiers of selection.Proof.
Combining (1.1) and (2.25) and simplifying, we find that it suffices to confirm that for r>1 and N=3,4,…,10, the inequality ρ1(N;r)>ρ1(M)(N;r) holds, which is equivalent to J(N,r)>0, where
(2.27)J(N,r):=(r+1)(rN-1)[N2+(r-3)(N-1),r(N-1)+1]N-3-(N-1+r)[rN(r(N-1)+1)N-2-(N-1+r)N-2].
By certain algebraic calculations, it is easy to derive
(2.28)J(3,r)=r3-r2-r+1=(r-1)2(r+1)>0,J(4,r)=12r5+4r4-8r3-24r2-4r+20=(r-1)2(12r3+28r2+36r+20)>0,J(5,r)=144r7+72r6+9r5-63r4-288r3-252r2+135r+243=(r-1)2(144r5+360r4+585r3+747r2+621r+243)>0.
By using Maple 10, we derive the polynomials J(N,r)/(r-1)2 for N=6,7,…,10, as shown in Table 1.
It follows directly from Table1 that J(N,r)>0 holds for N=6,7,8,9,10, which along with (2.28) leads to this result.
By (2.28) and Table 1, we guess that for N≥3 and (2.27) have the following form:
(2.29)J(N,r)=(r-1)2∑j=02N-5‍ajrj,N≥3,
where each aj(j=0,1,…,2N-5) is a nonnegative integer. However, it seems very difficult to prove this inequality in its general form even by the mathematical softwares, since we confront complicated polynomials as N increases, as shown in Table 1.Table 1
N
J
(
N
,
r
)
/
(
r
-
1
)
2
6
2000
r
7
+
5200
r
6
+
8640
r
5
+
12096
r
4
+
14928
r
3
+
14160
r
2
+
8992
r
+
3104
7
32400
r
9
+
86400
r
8
+
145800
r
7
+
205800
r
6
+
265825
r
5
+
318075
r
4
+
317225
r
3
+
243475
r
2
+
139125
r
+
46625
8
605052
r
11
+
1642284
r
10
+
2802996
r
9
+
3981348
r
8
+
5160960
r
7
+
6340608
r
6
+
7402608
r
5
+
7540272
r
4
+
6332004
r
3
+
4444596
r
2
+
2468988
r
+
823500
9
12845056
r
13
+
35323904
r
12
+
60813312
r
11
+
86804480
r
10
+
112842688
r
9
+
138883248
r
8
+
164923857
r
7
+
188867315
r
6
+
194198613
r
5
+
171318455
r
4
+
133184793
r
3
+
91162491
r
2
+
50328733
r
+
16777159
10
306110016
r
15
+
850305600
r
14
+
1473863040
r
13
+
2112117120
r
12
+
2752004160
r
11
+
3392000064
r
10
+
4032000000
r
9
+
4672000000
r
8
+
5268953280
r
7
+
5440222400
r
6
+
4946130304
r
5
+
4076784000
r
4
+
3098682560
r
3
+
2109150400
r
2
+
1162256640
r
+
387420416
⋮
⋮Another utility of Theorem2.3 is to study the asymptotic properties of the fixation probabilities on complete star digraphs. For that purpose, we rewrite (2.22) in the following form:
(2.30)ρk(N;r)=[1+1(r(N-1))]-k×{k(N-2)(r2-1)((N-1)2r2k)(1+rN-1)k-2[1+1(r(N-1))]k-1r2k(1+r(N-1))k-k(N-2)(r2-1)((N-1)2r2k)(1+rN-1)k-2}/[1-1r2N-2(1+(r-1r))/((N-1+1r))N-2].
By letting r>1 and taking limits on both sides of (2.30) in N, then ρk(N;r)→1-1/r2k as N→∞. Lieberman et al. [4] declared this result without rigorous argument. This statement tells us that for sufficiently large N, 𝒮N is an amplifier of selection because limN→∞ρ1(N;r)=1-1/r2>1-1/r=limN→∞ρ1(M)(N;r).Figure2 shows how the first-order fixation probability goes to the limit. For comparative purposes, the first-order fixation probability for the Moran process is also given in this figure. By Theorem 2.4 and a close look at Figure 2, we pose the conjecture: let r>1. Then, 𝒮N(N≥3) is an amplifier of selection.First-order fixation probabilities of CSDs and Moran process withN by fixing relative fitness r.
(a)
r
=
1.05
(b)
r
=
1.3
(c)
r
=
1.9
(d)
r
=
2.7In the proof of Theorem2.3, we argue that the conjecture holds for 3≤N≤10 by computing some polynomials; however, we confront complex polynomials when N becomes large. Thus, it seems very difficult to generally prove the conjecture. Fortunately, here we derive the following weak but useful result through rigorous reasoning.Theorem 2.5.
Letτ(N)>0, depending on N, be sufficiently small. Then for any r∈(1,1+τ(N)), 𝒮N(N≥3) is an amplifier of selection.Proof.
It suffices to prove that forr∈(1,1+τ(N)) the following inequality holds:
(2.31)ρ1(N;r)>ρ1(M)(N;r),
which is equivalent to ρ1(N;r)>ρ1(M)(N;r) by applying (1.1) and (2.25). We prove this as follows.
Letf(r)=(r+1)[N2+(N-1)(r-3)][r(N-1)+1]N-3, then
(2.32)f′(r)=f(r)+(N-1)(r+1)[r(N-1)+1]N-3+(N-3)(N-1)(r+1)[N2+(N-1)(r-3)][r(N-1)+1]N-4.
It follows from (2.27) that
(2.33)J(N,r)=(rN-1)f(r)-(N-1+r)[rN(r(N-1)+1)N-2-(N-1+r)N-2]
whose first- and second-order derivatives in r are, respectively, given by
(2.34)∂J(N,r)∂r=NrN-1f(r)+(rN-1)f′(r)-rN[r(N-1)+1]N-2-N(N-1+r)rN-1[r(N-1)+1]N-2+(N-1)(N-1+r)N-2-(N-2)(N-1)(N-1+r)rN[r(N-1)+1]N-3,∂2J(N,r)∂r2=N(N-1)rN-2f(r)+2NrN-1f′(r)+(rN-1)f′′(r)-2NrN-1[r(N-1)+1]N-2-N(N-1)(N-1+r)[r2(N-1)+r]N-2-(N-2)(N-1)r2[r2(N-1)+r]N-3[2N(N-1)+(2N+1)r]+(N-2)(N-1)[(N-1+r)N-3-rN(r(N-1)+1)N-3]-(N-3)(N-2)(N-1)2(N-1+r)rN[r(N-1)+1]N-4.
Through certain calculations, we derive
(2.35)f(1)=2[N2-2N+2]NN-3,f′(1)=NN-1+2(N-3)(N-1)(N2-2N+2)NN-4,
which imply
(2.36)∂J(N,r)∂r∣r=1=0,∂2J(N,r)∂r2∣r=1=2(N-1)(N-2)3NN-3>0.
It follows from the Taylor’s theorem and the fact J(N,1)=0 that
(2.37)J(N,r)=J(N,1)+∂J(N,r)∂r∣r=1(r-1)+12!∂2J(N,r)∂r2∣r=1(r-1)2+O((r-1)2)=(N-1)(N-2)3NN-3(r-1)2+O((r-1)2)=(r-1)2[(N-1)(N-2)3NN-3+O(r-1)].
Therefore, there exists a small enough τ(N)>0 such that the inequality (N-1)(N-2)3NN-3+O(r-1)>0 holds for arbitrary r∈(1,1+τ(N)). The proof is complete.Theorem2.5 shows that no matter the size of a complete star digraph, it is always an amplifier of selection if the mutant’s relative fitness is slightly larger than 1.In the following, let us examine the asymptotic behavior ofρk(N;r) provided that k=kN depends on N. If r>1 and SN(N≥3) is an amplifier of selection, then by taking limits on (2.30), we have (I)if kN→+∞, then ρkN(N;r)→1; (II) if kN→k*<+∞, then ρkN(N;r)→1-1/r2k*. The first assertion is interesting because it demonstrates that no matter how slowly the initial number of mutants increases (e.g., proportional to [lnlnN]) the corresponding fixation probability will always approach one.Whenr=rN depends on N, one may imagine that the asymptotic behavior of ρk(N;r) would become much more complex. For example, given k and suppose r=rN. Then we have the assertions: (I) if rN→r*>1, then ρk(N;rN)→1-1/(r*)2k; (II) if rN→1, then ρk(N;rN)→0. Note that the first assertion follows directly from (2.30). Considering the second assertion, assume that rNN→δ>1, then it follows from (2.30) that limN→∞ρk(N;rN)=limN→∞(1-1/rN2k)/(1-1/rN2N-2)=0/(1-1/δ2)=0. Finally, in the case rNN→1, the second assertion follows naturally.Figure3 shows that how the kth-order fixation probability changes as the increasing number mutants k.Figure 3
The curves ofkth-order fixation probability on 𝒮50 with respect to k by fixing r. These curves show the values of ρk dramatically converge to 1, and they increase as k increases, which is consistent with the fact that the more initial mutants there are, the larger the fixation probability is.


## 2.1. Basic Concepts and Notations

Acomplete star digraph (CSD) is a digraph with a single central vertex such that (I) there exists an edge from the center vertex to each peripheral vertex; (II) there exists an edge from each peripheral vertex to the center vertex; (III) there exist no other edges. Therefore, a CSD of size N, denoted by 𝒮N, is a digraph with vertex set 𝒱={1,2,…,N} and directed edge set ℰ={〈1,j〉,〈j,1〉∣j=2,…,N}. Figure 1 illustrates one CSD with intrinsic weights. Note that the term “CSD” means “CSD with intrinsic weights” in the sequel.Figure 1
One complete star digraph𝒮6 with intrinsic weights, the central vertex is indexed 1.Thekth-order fixation probability on 𝒮N is denoted by ρk(N;r), which represents the probability of the event that these k mutants generate a lineage that takes over the whole population.For technical reasons, we need the following notations in the sequel.The configuration of a population on𝒮N, at time t, is depicted by a vector Θ(t)=(θ1(t),θ2(t)), where θ1(t)=1 or 0 if a mutant occupies the central vertex or not, respectively, and θ2(t) denotes the number of mutants occupying vertices 2,…,N. The total number of mutants at time t is denoted by θ(t), that is, θ(t)=θ1(t)+θ2(t). The probability of the event that the mutants finally fixate by starting with Θ(0)=(k1,k2), is denoted by ρk1,k2(N;r) abbreviated to ρk1,k2 (without ambiguity).


## 2.2. Explicit Fixation Probability Formulae

Before exactly formulating thekth-order fixation probability, the following two necessary lemmas are presented.Lemma 2.1.
Consider a CSD of sizeN≥3 and denote by r the relative fitness of a mutant. Then for k=1,…,N-1, the difference equation system holds:
(2.1)pk=Apk-1
where pk=(ρ1,k,ρ0,k)T∈ℝ2, and
(2.2)A=(α1-ααβ1-αβ),
with α:=1+(N-1)/r,β:=r(N-1)/[r(N-1)+1].Proof.
This proof proceeds by calculating some conditional probabilities and by the total probability formula. LetΘ(0), at time t=0, represent the initial configuration of a population on 𝒮N, Θ(1) represent the new configuration after one step-time, and k=θ(0)=θ1(0)+θ2(0) represent the initial total number of mutants. Then two cases are discussed with respect to Θ(0).
Case1.  (The initial configuration Θ(0)=(1,k-1)). In order to reduce the k mutants by one, a resident individual must be selected to reproduce and, meanwhile the center mutant (staying at vertex 1) has to be chosen for death, whose conditional on the resident’s selection, happens with probability 1. Therefore, the conditional probability of going from Θ(0)=(1,k-1) to Θ(1)=(0,k-1) is given by
(2.3)Pr{Θ(1)=(0,k-1)∣Θ(0)=(1,k-1)}=(N-k)(rk+N-k).
On the other hand, to increaseω by one, the center mutant has to be chosen for reproduction, and one resident individual connected to the center vertex must be selected to die. Thus the conditional probability of going from Θ(0)=(1,k-1) to Θ(1)=(1,k) is given by
(2.4)Pr{Θ(1)=(1,k)∣Θ(0)=(1,k-1)}=r(N-k)[(rk+N-k)(N-1)].
By employing (2.3) and (2.4), the probability of the event that the configuration stays unchanged is given by
(2.5)Pr{Θ(1)=(1,k-1)∣Θ(0)=(1,k-1)}=1-(N-k)(rk+N-k)-r(N-k)[(rk+N-k)(N-1)].
Case2. (The initial configuration Θ(0)=(0,k)). On the one hand, in order to increase k by one, the center individual has to be replaced by a new mutant produced by another. Therefore, the probability of going from Θ(0)=(0,k) to Θ(1)=(1,k) is given by
(2.6)Pr{Θ(1)=(1,k)∣Θ(0)=(0,k)}=rk(rk+N-k).
On the other hand, to reducek by one, the center resident individual has to be chosen for reproduction and its offspring must replace one mutant. Therefore, the probability of going from Θ(0)=(0,k) to Θ(1)=(0,k-1) is given by
(2.7)Pr{Θ(1)=(0,k-1)∣Θ(0)=(0,k)}=k[(rk+N-k)(N-1)].
It follows from (2.6) and (2.7) that the probability of the configuration remaining unchanged is
(2.8)Pr{Θ(1)=(0,k)∣Θ(0)=(0,k)}=1-rk(rk+N-k)-k[(rk+N-k)(N-1)].
By the total probability formula, we derive that for k=1,…,N-1,
(2.9)ρ1,k-1=Pr{Θ(1)=(0,k-1)∣Θ(0)=(1,k-1)}ρ0,k-1+Pr{Θ(1)=(1,k)∣Θ(0)=(1,k-1)}ρ1,k+Pr{Θ(1)=(1,k-1)∣Θ(0)=(1,k-1)}ρ1,k-1,(2.10)ρ0,k=Pr{Θ(1)=(1,k)∣Θ(0)=(0,k)}ρ1,k+Pr{Θ(1)=(0,k-1)∣Θ(0)=(0,k)}ρ0,k-1+Pr{Θ(1)=(0,k)∣Θ(0)=(0,k)}ρ0,k.
Plugging (2.3)–(2.8) into (2.9)-(2.10) and simplifying the results, we deduce
(2.11)ρ1,k=αρ1,k-1+(1-α)ρ0,k-1,(2.12)ρ0,k=βρ1,k+(1-β)ρ0,k-1,
where α:=1+(N-1)/r,β:=r(N-1)/[r(N-1)+1].
By substituting (2.11) into (2.12) and rearranging the terms, we get
(2.13)ρ0,k=αβρ1,k-1+(1-αβ)ρ0,k-1.
Equation on (2.1) follows by combining (2.11) and (2.13).Lemma 2.2.
LetN and r be given. Then for 1≤k≤N-1, the following hold
(2.14)ρ1,k-1=[β-(α-1)αk-2(1-β)k-1][β-(α-1)αN-2(1-β)N-1],ρ0,k=[β(1-αk(1-β)k)][β-(α-1)αN-2(1-β)N-1].Proof.
We solve the linear difference equation system (2.1) in view of the standard technique given in Elaydi [15]. The matrix A in (2.2) has two real eigenvalues λ1=1,λ2=α(1-β).
Furthermore,η1=(1,1)T and η2=(α-1,αβ)T are eigenvectors of A corresponding to λ1,λ2, respectively. Denote
(2.15)D=(1α-11αβ),
and its inverse matrix
(2.16)D-1=[1(1-α+αβ)](αβ1-α-11).
Through the method of diagonalizing matrix A, we make the variable change:
(2.17)(ρ1,kρ0,k)=D(xkyk),
and thus
(2.18)(xkyk)=D-1AD(xk-1yk-1)=(100α(1-β))(xk-1yk-1)=(100αk(1-β)k)(x0y0).
Therefore,
(2.19)(ρ1,kρ0,k)=D(100αk(1-β)k)D-1(ρ1,0ρ0,0)=ρ1,0(1-α+αβ)(αβ-(α-1)αk(1-β)kαβ-αk+1β(1-β)k).
Particularly, we know
(2.20)1=ρ1,N-1=[ρ1,0(1-α+αβ)]×[αβ-(α-1)αN-1(1-β)N-1],
which implies
(2.21)ρ1,0=(1-α+αβ)(αβ-(α-1)αN-1(1-β)N-1).
Equation on (2.14) follows from substituting (2.21) into (2.19). The proof is complete.Next, we present the main theorem of this section.Theorem 2.3.
Letr and N>2 be given. Then for each k=1,…,N-1, the explicit kth-order fixation probability on 𝒮N is given by
(2.22)ρk(N;r)=rN[r(N-1)+1]N-k-2[rN(r(N-1)+1)N-2-(N-1+r)N-2]×[((N-1+r)r)k(r(N-1)+1)k-(krk)(N-2)(r2-1)(N-1+r)k-2-((N-1+r)r)k].Proof.
In light of (2.21), we derive
(2.23)ρk(N;r)=kNρ1,k-1+(1-kN)ρ0,k=k/Nαk-2(1-β)k-1[α2β(1-β)-(α-1)]+β[1-αk(1-β)k]β-(α-1)αN-2(1-β)N-1.
By substitutingα and β (see Lemma 2.1) into this equation, we have
(2.24)ρk(N;r)={k/N[1+(N-1)/r]k-2[1/(r(N-1)+1)]k-1[(1+(N-1)/r)2(r(N-1))/((r(N-1)+1)2)-(N-1)/r]+r(N-1)/(r(N-1)+1)[1-(1+(N-1)/r)k/(r(N-1)+1)k]}/[(r(N-1)+1)N-1(r(N-1))/(r(N-1)+1)-(N-1)/r[1+(N-1)/r]N-2/(r(N-1)+1)N-1]=[rN-1(r(N-1)+1)N-1]{((r(N-1)+1)k+1)k(N-1)(r2-1)(2-N)(N+r-1)k-2/[rk-1(r(N-1)+1)k+1]+r(N-1)/((r(N-1)+1)k+1)[(r(N-1)+1)k-((r+N-1)/r)k]}/{(N-1)[rN(r(N-1)+1)N-2-(r+N-1)N-2]}=rN[r(N-1)+1]N-k-2/[rN(r(N-1)+1)N-2-(N-1+r)N-2]×[(r(N-1)+1)k-(k/rk)(N-2)(r2-1)(N-1+r)k-2-((N-1+r)/r)k].
The proof is complete.Particularly, for the casek=1, (2.22) becomes
(2.25)ρ1(N;r)=rN-1(r2-1)(r(N-1)+1)N-3[rN(r(N-1)+1)N-2-(N-1+r)N-2]×[(N2+(r-3)(N-1))(N-1+r)].Here we point out that (2.25) is just the same as the formula ϱ=[n2r/(nr+1)+r/(n+r)]/[(n+1)[1+n/(n+r)∑j=1n-1‍((n+r)/(r(nr+1)))j]], which is given in Broom and Rychtář [7], because
(2.26)ϱ=n2r/(nr+1)+r/(n+r)(n+1)[1+n/(n+r)∑j=1n-1((n+r)/(nr2+r))j]=n2r/(nr+1)+r/(n+r)(n+1)[1+[rn-1(nr+1)n-1-(n+r)n-1]/[(r2-1)(r(nr+1))n-1]]=[n(n+1)r2+(n+1)(n2-n+1)r](r2-1)[r(nr+1)]n-1(nr+1)(n+r)(n+1)[rn+1(nr+1)n-1-(n+r)n-1]=rN-1(r2-1)[r(N-1)+1]N-3rN[r(N-1)+1]N-2-[N-1+r]N-2×(N2+(r-3)(N-1))(N-1+r)=ρ1(N;r),
where n represents the number of all individuals except the central one, that is,  n=N-1.


## 2.3. Applications and Properties of the Fixation Probability

Theorem2.3 allows us to calculate the fixation probability for a given CSD, so here we use it to answer whether a given CSD is an amplifier of selection or not. We have the following result.Theorem 2.4.
𝒮
3
,
𝒮
4
,
𝒮
5
,
𝒮
6
,
𝒮
7
,
𝒮
8
,
𝒮
9, and 𝒮10 are all amplifiers of selection.Proof.
Combining (1.1) and (2.25) and simplifying, we find that it suffices to confirm that for r>1 and N=3,4,…,10, the inequality ρ1(N;r)>ρ1(M)(N;r) holds, which is equivalent to J(N,r)>0, where
(2.27)J(N,r):=(r+1)(rN-1)[N2+(r-3)(N-1),r(N-1)+1]N-3-(N-1+r)[rN(r(N-1)+1)N-2-(N-1+r)N-2].
By certain algebraic calculations, it is easy to derive
(2.28)J(3,r)=r3-r2-r+1=(r-1)2(r+1)>0,J(4,r)=12r5+4r4-8r3-24r2-4r+20=(r-1)2(12r3+28r2+36r+20)>0,J(5,r)=144r7+72r6+9r5-63r4-288r3-252r2+135r+243=(r-1)2(144r5+360r4+585r3+747r2+621r+243)>0.
By using Maple 10, we derive the polynomials J(N,r)/(r-1)2 for N=6,7,…,10, as shown in Table 1.
It follows directly from Table1 that J(N,r)>0 holds for N=6,7,8,9,10, which along with (2.28) leads to this result.
By (2.28) and Table 1, we guess that for N≥3 and (2.27) have the following form:
(2.29)J(N,r)=(r-1)2∑j=02N-5‍ajrj,N≥3,
where each aj(j=0,1,…,2N-5) is a nonnegative integer. However, it seems very difficult to prove this inequality in its general form even by the mathematical softwares, since we confront complicated polynomials as N increases, as shown in Table 1.Table 1
N
J
(
N
,
r
)
/
(
r
-
1
)
2
6
2000
r
7
+
5200
r
6
+
8640
r
5
+
12096
r
4
+
14928
r
3
+
14160
r
2
+
8992
r
+
3104
7
32400
r
9
+
86400
r
8
+
145800
r
7
+
205800
r
6
+
265825
r
5
+
318075
r
4
+
317225
r
3
+
243475
r
2
+
139125
r
+
46625
8
605052
r
11
+
1642284
r
10
+
2802996
r
9
+
3981348
r
8
+
5160960
r
7
+
6340608
r
6
+
7402608
r
5
+
7540272
r
4
+
6332004
r
3
+
4444596
r
2
+
2468988
r
+
823500
9
12845056
r
13
+
35323904
r
12
+
60813312
r
11
+
86804480
r
10
+
112842688
r
9
+
138883248
r
8
+
164923857
r
7
+
188867315
r
6
+
194198613
r
5
+
171318455
r
4
+
133184793
r
3
+
91162491
r
2
+
50328733
r
+
16777159
10
306110016
r
15
+
850305600
r
14
+
1473863040
r
13
+
2112117120
r
12
+
2752004160
r
11
+
3392000064
r
10
+
4032000000
r
9
+
4672000000
r
8
+
5268953280
r
7
+
5440222400
r
6
+
4946130304
r
5
+
4076784000
r
4
+
3098682560
r
3
+
2109150400
r
2
+
1162256640
r
+
387420416
⋮
⋮Another utility of Theorem2.3 is to study the asymptotic properties of the fixation probabilities on complete star digraphs. For that purpose, we rewrite (2.22) in the following form:
(2.30)ρk(N;r)=[1+1(r(N-1))]-k×{k(N-2)(r2-1)((N-1)2r2k)(1+rN-1)k-2[1+1(r(N-1))]k-1r2k(1+r(N-1))k-k(N-2)(r2-1)((N-1)2r2k)(1+rN-1)k-2}/[1-1r2N-2(1+(r-1r))/((N-1+1r))N-2].
By letting r>1 and taking limits on both sides of (2.30) in N, then ρk(N;r)→1-1/r2k as N→∞. Lieberman et al. [4] declared this result without rigorous argument. This statement tells us that for sufficiently large N, 𝒮N is an amplifier of selection because limN→∞ρ1(N;r)=1-1/r2>1-1/r=limN→∞ρ1(M)(N;r).Figure2 shows how the first-order fixation probability goes to the limit. For comparative purposes, the first-order fixation probability for the Moran process is also given in this figure. By Theorem 2.4 and a close look at Figure 2, we pose the conjecture: let r>1. Then, 𝒮N(N≥3) is an amplifier of selection.First-order fixation probabilities of CSDs and Moran process withN by fixing relative fitness r.
(a)
r
=
1.05
(b)
r
=
1.3
(c)
r
=
1.9
(d)
r
=
2.7In the proof of Theorem2.3, we argue that the conjecture holds for 3≤N≤10 by computing some polynomials; however, we confront complex polynomials when N becomes large. Thus, it seems very difficult to generally prove the conjecture. Fortunately, here we derive the following weak but useful result through rigorous reasoning.Theorem 2.5.
Letτ(N)>0, depending on N, be sufficiently small. Then for any r∈(1,1+τ(N)), 𝒮N(N≥3) is an amplifier of selection.Proof.
It suffices to prove that forr∈(1,1+τ(N)) the following inequality holds:
(2.31)ρ1(N;r)>ρ1(M)(N;r),
which is equivalent to ρ1(N;r)>ρ1(M)(N;r) by applying (1.1) and (2.25). We prove this as follows.
Letf(r)=(r+1)[N2+(N-1)(r-3)][r(N-1)+1]N-3, then
(2.32)f′(r)=f(r)+(N-1)(r+1)[r(N-1)+1]N-3+(N-3)(N-1)(r+1)[N2+(N-1)(r-3)][r(N-1)+1]N-4.
It follows from (2.27) that
(2.33)J(N,r)=(rN-1)f(r)-(N-1+r)[rN(r(N-1)+1)N-2-(N-1+r)N-2]
whose first- and second-order derivatives in r are, respectively, given by
(2.34)∂J(N,r)∂r=NrN-1f(r)+(rN-1)f′(r)-rN[r(N-1)+1]N-2-N(N-1+r)rN-1[r(N-1)+1]N-2+(N-1)(N-1+r)N-2-(N-2)(N-1)(N-1+r)rN[r(N-1)+1]N-3,∂2J(N,r)∂r2=N(N-1)rN-2f(r)+2NrN-1f′(r)+(rN-1)f′′(r)-2NrN-1[r(N-1)+1]N-2-N(N-1)(N-1+r)[r2(N-1)+r]N-2-(N-2)(N-1)r2[r2(N-1)+r]N-3[2N(N-1)+(2N+1)r]+(N-2)(N-1)[(N-1+r)N-3-rN(r(N-1)+1)N-3]-(N-3)(N-2)(N-1)2(N-1+r)rN[r(N-1)+1]N-4.
Through certain calculations, we derive
(2.35)f(1)=2[N2-2N+2]NN-3,f′(1)=NN-1+2(N-3)(N-1)(N2-2N+2)NN-4,
which imply
(2.36)∂J(N,r)∂r∣r=1=0,∂2J(N,r)∂r2∣r=1=2(N-1)(N-2)3NN-3>0.
It follows from the Taylor’s theorem and the fact J(N,1)=0 that
(2.37)J(N,r)=J(N,1)+∂J(N,r)∂r∣r=1(r-1)+12!∂2J(N,r)∂r2∣r=1(r-1)2+O((r-1)2)=(N-1)(N-2)3NN-3(r-1)2+O((r-1)2)=(r-1)2[(N-1)(N-2)3NN-3+O(r-1)].
Therefore, there exists a small enough τ(N)>0 such that the inequality (N-1)(N-2)3NN-3+O(r-1)>0 holds for arbitrary r∈(1,1+τ(N)). The proof is complete.Theorem2.5 shows that no matter the size of a complete star digraph, it is always an amplifier of selection if the mutant’s relative fitness is slightly larger than 1.In the following, let us examine the asymptotic behavior ofρk(N;r) provided that k=kN depends on N. If r>1 and SN(N≥3) is an amplifier of selection, then by taking limits on (2.30), we have (I)if kN→+∞, then ρkN(N;r)→1; (II) if kN→k*<+∞, then ρkN(N;r)→1-1/r2k*. The first assertion is interesting because it demonstrates that no matter how slowly the initial number of mutants increases (e.g., proportional to [lnlnN]) the corresponding fixation probability will always approach one.Whenr=rN depends on N, one may imagine that the asymptotic behavior of ρk(N;r) would become much more complex. For example, given k and suppose r=rN. Then we have the assertions: (I) if rN→r*>1, then ρk(N;rN)→1-1/(r*)2k; (II) if rN→1, then ρk(N;rN)→0. Note that the first assertion follows directly from (2.30). Considering the second assertion, assume that rNN→δ>1, then it follows from (2.30) that limN→∞ρk(N;rN)=limN→∞(1-1/rN2k)/(1-1/rN2N-2)=0/(1-1/δ2)=0. Finally, in the case rNN→1, the second assertion follows naturally.Figure3 shows that how the kth-order fixation probability changes as the increasing number mutants k.Figure 3
The curves ofkth-order fixation probability on 𝒮50 with respect to k by fixing r. These curves show the values of ρk dramatically converge to 1, and they increase as k increases, which is consistent with the fact that the more initial mutants there are, the larger the fixation probability is.


## 3. Fixation Probabilities on Complete Bipartite Digraphs

This section gives a recursive equation regarding the fixation probabilities on complete bipartite digraphs and obtains some results through a numerical method by simulation.


### 3.1. Basic Concepts and Notations

Acomplete bipartite digraph (CBD) is a digraph whose vertices are partitioned into two partite sets so that there is no edge connecting any two vertices in the same partite set, and there is an directed edge connecting each vertex in one partite set to each vertex in the other partite set [16]. A CBD denoted by ℬN1,N2(here N1  and N2 represent the sizes of two partite sets, respectively, and due to symmetry, let N1≤N2) is a digraph with vertex set 𝒱={1,2,…,N} (N=N1+N2 is the total size) and edge set ℰ={〈i,j〉,〈j,i〉∣i=1,2,…,N1,j=N1+1,…,N1+N2}. A CBD is balanced if its two partite sets have the same size, that is, N1=N2, otherwise it is unbalanced. The unbalance degree of ℬN1,N2 can be measured by |N2-N1|. Figure 4 depicts two CBDs (the second one’s unbalance degree is 1). The term “CBD” implies “CBD with intrinsic weights” in the sequel.Two complete bipartite digraphs. (a) Shows a general CBD, while (b) shows an unbalanced one with intrinsic weights.
(a)
ℬ
N
1
,
N
2
(b)
ℬ
2,3 with intrinsic weightsWithout loss of generality, we suppose here thatN1≤N2. Obviously, when N1=1, the complete bipartite digraph is just the complete star digraph 𝒮N2. A CBD is an isothermal digraph if and only if it is balanced, thus the first-order fixation probability on a balanced CBD is already known by the isothermal theorem.For technical reasons, the following notations are needed in the sequel.Consider a homogeneous population on a CBD where all individuals have fitness 1. Suppose thatk new mutants with relative fitness r are introduced by placing them on k randomly chosen vertices. Let ρk(N1,N2;r) represent the kth-order fixation probability, that is, the probability of the event that these k mutants generate a lineage that takes over the population.At timet, the configuration of a population on ℬN1,N2 is described by a vector Ω(t)=(m1(t),m2(t)), where m1(t) and m2(t) represent the number of mutants staying at vertices 1,2,…,N1 and vertices N1+1,N1+2,…,N1+N2, respectively.Letm(t) represent the total number of mutants at time t, that is, m(t)=m1(t)+m2(t), and ρk1,k2(N1,N2;r) (without ambiguity, ρk1,k2) represent the probability of the event that, starting from Ω(0)=(k1,k2), the mutants finally fixate.


### 3.2.  A Numerical Method

Here we will give an approach for calculating the fixation probability on a CBD and thus three theorems are established as follows.Theorem 3.1.
Consider a CBDℬN1,N2 with N2≥N1≥2. Then, we have ρk1,k2(N1,N2;r)(0≤k1≤N1,0≤k2≤N2) satisfying the recursive equation:
(3.1)χ1ρk1-1,k2+χ2ρk1,k2-1+χ3ρk1+1,k2+χ4ρk1,k2+1=(χ1+χ2+χ3+χ4)ρk1,k2,
where χ1:=k1(N2-k2)/N1,χ2:=k2(N1-k1)/N2,χ3:=rk2(N1-k1)/N1,χ4:=rk1(N2-k2)/N2, and ρ0,0=0,ρN1,N2=1,ρ-1,i=ρj,-1=ρN1+1,i=ρj,N2+1=0,i=0,1,…,N2,j=0,1,…,N1, are boundary conditions.Proof.
LetΩ(0), at time t=0, represent the initial configuration of a population on ℬN1,N2 and Ω(1) represent the new configuration after one step-time. Thus we have m1(0)=k1,m2(0)=k2,m(0)=k1+k2. It is simple to calculate the following conditional probabilities:
(3.2)Pr{Ω(1)=(k1+1,k2)∣Ω(0)=(k1,k2)}=rk2(N1-k1)(N1(rk+N-k)),Pr{Ω(1)=(k1,k2+1)∣Ω(0)=(k1,k2)}=rk1(N2-k2)(N2(rk+N-k)),Pr{Ω(1)=(k1-1,k2)∣Ω(0)=(k1,k2)}=k1(N2-k2)(N1(rk+N-k)),Pr{Ω(1)=(k1,k2-1)∣Ω(0)=(k1,k2)}=k2(N1-k1)(N2(rk+N-k)),
where k=k1+k2,N=N1+N2.
Therefore, the probability that the configuration does not change is given by(3.3)Pr{Ω(1)=(k1,k2)∣Ω(0)=(k1,k2)}=1-Pr{Ω(1)=(k1+1,k2)∣Ω(0)=(k1,k2)}-Pr{Ω(1)=(k1,k2+1)∣Ω(0)=(k1,k2)}-Pr{Ω(1)=(k1-1,k2)∣Ω(0)=(k1,k2)}-Pr{Ω(1)=(k1,k2-1)∣Ω(0)=(k1,k2)}.
By the total probability formula, we derive
(3.4)ρk1,k2=Pr{Ω(1)=(k1+1,k2)∣Ω(0)=(k1,k2)}ρk1+1,k2+Pr{Ω(1)=(k1,k2+1)∣Ω(0)=(k1,k2)}ρk1,k2+1+Pr{Ω(1)=(k1-1,k2)∣Ω(0)=(k1,k2)}ρk1-1,k2+Pr{Ω(1)=(k1,k2-1)∣Ω(0)=(k1,k2)}ρk1,k2-1+Pr{Ω(1)=(k1,k2)∣Ω(0)=(k1,k2)}ρk1,k2.
Note that ρ0,0=0,ρN1,N2=1 and ρ-1,i,ρj,-1,ρN1+1,i,ρj,N2+1,i=0,1,…,N2,j=0,1,…,N-1 are meaningless, thus we assume ρ-1,i=ρj,-1=ρN1+1,i=ρj,N2+1=0, i=0,1,…,N2, j=0,1,…,N-1. Then (3.1) follows immediately by plugging (3.2) and (3.3) into (3.4) and simplifying. The proof is complete.Theorem 3.2.
For1≤k≤N-1, the following holds
(3.5)ρk(N1,N2;r)=∑k1=max⁡⁡{k,N2}-N2min⁡⁡{k,N1}(N1k1)(N2k-k1)(N1+N2k)ρk1,k-k1(N1,N2;r).
Particularly,
(3.6)ρ1(N1,N2;r)=N2(N1+N2)ρ0,1(N1,N2;r)+N1(N1+N2)ρ0,1(N1,N2;r).
Sketch. The randomness of the initial configuration of mutants gives
(3.7)Pr{Ω(0)=(k1,k-k1)}=(N1k1)(N2k-k1)(N1+N2k),
where k1 represents the number of mutants staying at vertices 1,2,…,N1.
Equation on (3.5) follows from this equation and the following total probability formula:
(3.8)ρk(N1,N2;r)=∑k1=max⁡⁡{k,N2}-N2min⁡⁡{k,N1}Pr{Ω(0)=(k1,k-k1)}ρk1,k-k1(N1,N2;r).From the above theorems, calculatingρ1(N1,N2;r) reduces to solving the difference equation system (3.1) with (N1+1)(N2+1)-2=N1N2+N1+N2-1 variables, because ρ0,0,ρN1,N2 are known. It is known [17] that solving a linear system in n variables with the Gaussian elimination method would take O(n2) time. So, we have the result.Theorem 3.3.
ρ
1
(
N
1
,
N
2
;
r
) is numerically obtainable in O(N12N22) time.
In the best caseρ1(1,N2;r) is computed in O(N22) time, whereas in the worst case ρ1(N1,N1;r) consumes up to O(N14) calculation time.


### 3.3. Numerical Experiments

Based on the theory established above, this section aims to experimentally reveal how a CBD, affects the first-order fixation probability (3.6) of an advantageous mutant. The experiments are done using Matlab.At first sight,ρ1(N1,N2;r) depends on three parameters. However, we will focus our attention on how ρ1 is affected by N1 and N2 because these two parameters characterize the CBD. Therefore, we will fix the fitness r in the subsequent experiments.Experiment 1.
To verify the theory’s correctness and guarantee the effectiveness of the following experiments, this experiment illustrates (see Figure5) that the fixation probability of balanced CBDs is equal to that of the corresponding Moran process. Here, the isothermal theorem of Lieberman et al. [4] is referenced.Figure 5
Comparisons between the first-order fixation probability (FOFP) of CBDs and that of the Moran process with different fixedr. The lines and little rectangles represent the FOFPs of the Moran process and CBDs over N, respectively. This figure shows that FOFPs of CBDs are always equal to those of the corresponding Moran process with the same size N.Experiment 2.
Givenr>1 and 𝒩k={2,3,…,k}, we wish to know how ρ1 depends on N1 and N2. For the purpose, we evaluate ρ1 for each pair (N1,N2) with N1,N2∈𝒩50. The resulting ρ1-(N1,N2) lattice graphs are plotted in Figure 6.The lattice graphs ofρ1 with respect to N1,N2 by fixing r. Note that the figures are symmetric due to the symmetry of CBDs.
(a)
r
=
1.1
(b)
r
=
1.3
(c)
r
=
2
(d)
r
=
2.5Experiment 3.
GivenN and r>1, it is interesting to know how ρ1 depends on N1 (without loss of generality, again let N1≤N2). Thus for each pair (r,N) with r∈{1.1,1.2,1.4,1.7,2} and N∈{52,102}, we evaluate ρ1 for all N1=2,…,[N/2]. The resulting ρ1-N1 curves are plotted in Figure 7, which suggests that ρ1 declines as N1 increases, that is,
(3.9)ρ1(1,N-1;r)>ρ1(2,N-2;r)>⋯>ρ1([N2],N-[N2];r).
In particular,ρ1 attains its maximum and minimum when N1=1 and N1=[N/2], respectively. This implies that all unbalanced CBDs under consideration are amplifiers of selection. Moreover, the more unbalanced a CBD, the more significantly it will amplify selection. On this basis, we pose the following conjecture.First-order fixation probability over the size of smaller partite set by fixing the total sizeN.
(a)
N
=
52
(b)
N
=
102Conjecture 3.4.
For anyr>1 and N≥4, the inequality chain (3.9) holds.Experiment 4.
Letr>1 and N1 be given, we expect to know how ρ1 depends on N. For that purpose, for each pair (r,N1) with r∈{1.2,1.6} and N1∈{2,3,4,5} we compute ρ1 for 2N1≤N≤250. The resulting ρ1-N curves are given in Figure 8. For comparison, then ρ1-N curves for the corresponding Moran processes are also depicted in this figure. The following conjecture is posed here.First-order fixation probability over the total population size. This figure shows the phenomenon that for each pair(r,N1) under consideration, ρ1(N1,N-N1;r) converges. Furthermore, the limit of ρ1 declines as N1 increases with fixed r.
(a)
r
=
1.2
(b)
r
=
1.6Conjecture 3.5.
For anyr>1 and any positive integer N1≥2, ρ1(N1,N-N1;r) converges. Moreover, the following inequality chain holds:
(3.10)lim⁡N→∞ρ1(2,N-2;r)<lim⁡N→∞ρ1(3,N-3;r)<lim⁡N→∞ρ1(4,N-4;r)<⋯.


## 3.1. Basic Concepts and Notations

Acomplete bipartite digraph (CBD) is a digraph whose vertices are partitioned into two partite sets so that there is no edge connecting any two vertices in the same partite set, and there is an directed edge connecting each vertex in one partite set to each vertex in the other partite set [16]. A CBD denoted by ℬN1,N2(here N1  and N2 represent the sizes of two partite sets, respectively, and due to symmetry, let N1≤N2) is a digraph with vertex set 𝒱={1,2,…,N} (N=N1+N2 is the total size) and edge set ℰ={〈i,j〉,〈j,i〉∣i=1,2,…,N1,j=N1+1,…,N1+N2}. A CBD is balanced if its two partite sets have the same size, that is, N1=N2, otherwise it is unbalanced. The unbalance degree of ℬN1,N2 can be measured by |N2-N1|. Figure 4 depicts two CBDs (the second one’s unbalance degree is 1). The term “CBD” implies “CBD with intrinsic weights” in the sequel.Two complete bipartite digraphs. (a) Shows a general CBD, while (b) shows an unbalanced one with intrinsic weights.
(a)
ℬ
N
1
,
N
2
(b)
ℬ
2,3 with intrinsic weightsWithout loss of generality, we suppose here thatN1≤N2. Obviously, when N1=1, the complete bipartite digraph is just the complete star digraph 𝒮N2. A CBD is an isothermal digraph if and only if it is balanced, thus the first-order fixation probability on a balanced CBD is already known by the isothermal theorem.For technical reasons, the following notations are needed in the sequel.Consider a homogeneous population on a CBD where all individuals have fitness 1. Suppose thatk new mutants with relative fitness r are introduced by placing them on k randomly chosen vertices. Let ρk(N1,N2;r) represent the kth-order fixation probability, that is, the probability of the event that these k mutants generate a lineage that takes over the population.At timet, the configuration of a population on ℬN1,N2 is described by a vector Ω(t)=(m1(t),m2(t)), where m1(t) and m2(t) represent the number of mutants staying at vertices 1,2,…,N1 and vertices N1+1,N1+2,…,N1+N2, respectively.Letm(t) represent the total number of mutants at time t, that is, m(t)=m1(t)+m2(t), and ρk1,k2(N1,N2;r) (without ambiguity, ρk1,k2) represent the probability of the event that, starting from Ω(0)=(k1,k2), the mutants finally fixate.


## 3.2.  A Numerical Method

Here we will give an approach for calculating the fixation probability on a CBD and thus three theorems are established as follows.Theorem 3.1.
Consider a CBDℬN1,N2 with N2≥N1≥2. Then, we have ρk1,k2(N1,N2;r)(0≤k1≤N1,0≤k2≤N2) satisfying the recursive equation:
(3.1)χ1ρk1-1,k2+χ2ρk1,k2-1+χ3ρk1+1,k2+χ4ρk1,k2+1=(χ1+χ2+χ3+χ4)ρk1,k2,
where χ1:=k1(N2-k2)/N1,χ2:=k2(N1-k1)/N2,χ3:=rk2(N1-k1)/N1,χ4:=rk1(N2-k2)/N2, and ρ0,0=0,ρN1,N2=1,ρ-1,i=ρj,-1=ρN1+1,i=ρj,N2+1=0,i=0,1,…,N2,j=0,1,…,N1, are boundary conditions.Proof.
LetΩ(0), at time t=0, represent the initial configuration of a population on ℬN1,N2 and Ω(1) represent the new configuration after one step-time. Thus we have m1(0)=k1,m2(0)=k2,m(0)=k1+k2. It is simple to calculate the following conditional probabilities:
(3.2)Pr{Ω(1)=(k1+1,k2)∣Ω(0)=(k1,k2)}=rk2(N1-k1)(N1(rk+N-k)),Pr{Ω(1)=(k1,k2+1)∣Ω(0)=(k1,k2)}=rk1(N2-k2)(N2(rk+N-k)),Pr{Ω(1)=(k1-1,k2)∣Ω(0)=(k1,k2)}=k1(N2-k2)(N1(rk+N-k)),Pr{Ω(1)=(k1,k2-1)∣Ω(0)=(k1,k2)}=k2(N1-k1)(N2(rk+N-k)),
where k=k1+k2,N=N1+N2.
Therefore, the probability that the configuration does not change is given by(3.3)Pr{Ω(1)=(k1,k2)∣Ω(0)=(k1,k2)}=1-Pr{Ω(1)=(k1+1,k2)∣Ω(0)=(k1,k2)}-Pr{Ω(1)=(k1,k2+1)∣Ω(0)=(k1,k2)}-Pr{Ω(1)=(k1-1,k2)∣Ω(0)=(k1,k2)}-Pr{Ω(1)=(k1,k2-1)∣Ω(0)=(k1,k2)}.
By the total probability formula, we derive
(3.4)ρk1,k2=Pr{Ω(1)=(k1+1,k2)∣Ω(0)=(k1,k2)}ρk1+1,k2+Pr{Ω(1)=(k1,k2+1)∣Ω(0)=(k1,k2)}ρk1,k2+1+Pr{Ω(1)=(k1-1,k2)∣Ω(0)=(k1,k2)}ρk1-1,k2+Pr{Ω(1)=(k1,k2-1)∣Ω(0)=(k1,k2)}ρk1,k2-1+Pr{Ω(1)=(k1,k2)∣Ω(0)=(k1,k2)}ρk1,k2.
Note that ρ0,0=0,ρN1,N2=1 and ρ-1,i,ρj,-1,ρN1+1,i,ρj,N2+1,i=0,1,…,N2,j=0,1,…,N-1 are meaningless, thus we assume ρ-1,i=ρj,-1=ρN1+1,i=ρj,N2+1=0, i=0,1,…,N2, j=0,1,…,N-1. Then (3.1) follows immediately by plugging (3.2) and (3.3) into (3.4) and simplifying. The proof is complete.Theorem 3.2.
For1≤k≤N-1, the following holds
(3.5)ρk(N1,N2;r)=∑k1=max⁡⁡{k,N2}-N2min⁡⁡{k,N1}(N1k1)(N2k-k1)(N1+N2k)ρk1,k-k1(N1,N2;r).
Particularly,
(3.6)ρ1(N1,N2;r)=N2(N1+N2)ρ0,1(N1,N2;r)+N1(N1+N2)ρ0,1(N1,N2;r).
Sketch. The randomness of the initial configuration of mutants gives
(3.7)Pr{Ω(0)=(k1,k-k1)}=(N1k1)(N2k-k1)(N1+N2k),
where k1 represents the number of mutants staying at vertices 1,2,…,N1.
Equation on (3.5) follows from this equation and the following total probability formula:
(3.8)ρk(N1,N2;r)=∑k1=max⁡⁡{k,N2}-N2min⁡⁡{k,N1}Pr{Ω(0)=(k1,k-k1)}ρk1,k-k1(N1,N2;r).From the above theorems, calculatingρ1(N1,N2;r) reduces to solving the difference equation system (3.1) with (N1+1)(N2+1)-2=N1N2+N1+N2-1 variables, because ρ0,0,ρN1,N2 are known. It is known [17] that solving a linear system in n variables with the Gaussian elimination method would take O(n2) time. So, we have the result.Theorem 3.3.
ρ
1
(
N
1
,
N
2
;
r
) is numerically obtainable in O(N12N22) time.
In the best caseρ1(1,N2;r) is computed in O(N22) time, whereas in the worst case ρ1(N1,N1;r) consumes up to O(N14) calculation time.


## 3.3. Numerical Experiments

Based on the theory established above, this section aims to experimentally reveal how a CBD, affects the first-order fixation probability (3.6) of an advantageous mutant. The experiments are done using Matlab.At first sight,ρ1(N1,N2;r) depends on three parameters. However, we will focus our attention on how ρ1 is affected by N1 and N2 because these two parameters characterize the CBD. Therefore, we will fix the fitness r in the subsequent experiments.Experiment 1.
To verify the theory’s correctness and guarantee the effectiveness of the following experiments, this experiment illustrates (see Figure5) that the fixation probability of balanced CBDs is equal to that of the corresponding Moran process. Here, the isothermal theorem of Lieberman et al. [4] is referenced.Figure 5
Comparisons between the first-order fixation probability (FOFP) of CBDs and that of the Moran process with different fixedr. The lines and little rectangles represent the FOFPs of the Moran process and CBDs over N, respectively. This figure shows that FOFPs of CBDs are always equal to those of the corresponding Moran process with the same size N.Experiment 2.
Givenr>1 and 𝒩k={2,3,…,k}, we wish to know how ρ1 depends on N1 and N2. For the purpose, we evaluate ρ1 for each pair (N1,N2) with N1,N2∈𝒩50. The resulting ρ1-(N1,N2) lattice graphs are plotted in Figure 6.The lattice graphs ofρ1 with respect to N1,N2 by fixing r. Note that the figures are symmetric due to the symmetry of CBDs.
(a)
r
=
1.1
(b)
r
=
1.3
(c)
r
=
2
(d)
r
=
2.5Experiment 3.
GivenN and r>1, it is interesting to know how ρ1 depends on N1 (without loss of generality, again let N1≤N2). Thus for each pair (r,N) with r∈{1.1,1.2,1.4,1.7,2} and N∈{52,102}, we evaluate ρ1 for all N1=2,…,[N/2]. The resulting ρ1-N1 curves are plotted in Figure 7, which suggests that ρ1 declines as N1 increases, that is,
(3.9)ρ1(1,N-1;r)>ρ1(2,N-2;r)>⋯>ρ1([N2],N-[N2];r).
In particular,ρ1 attains its maximum and minimum when N1=1 and N1=[N/2], respectively. This implies that all unbalanced CBDs under consideration are amplifiers of selection. Moreover, the more unbalanced a CBD, the more significantly it will amplify selection. On this basis, we pose the following conjecture.First-order fixation probability over the size of smaller partite set by fixing the total sizeN.
(a)
N
=
52
(b)
N
=
102Conjecture 3.4.
For anyr>1 and N≥4, the inequality chain (3.9) holds.Experiment 4.
Letr>1 and N1 be given, we expect to know how ρ1 depends on N. For that purpose, for each pair (r,N1) with r∈{1.2,1.6} and N1∈{2,3,4,5} we compute ρ1 for 2N1≤N≤250. The resulting ρ1-N curves are given in Figure 8. For comparison, then ρ1-N curves for the corresponding Moran processes are also depicted in this figure. The following conjecture is posed here.First-order fixation probability over the total population size. This figure shows the phenomenon that for each pair(r,N1) under consideration, ρ1(N1,N-N1;r) converges. Furthermore, the limit of ρ1 declines as N1 increases with fixed r.
(a)
r
=
1.2
(b)
r
=
1.6Conjecture 3.5.
For anyr>1 and any positive integer N1≥2, ρ1(N1,N-N1;r) converges. Moreover, the following inequality chain holds:
(3.10)lim⁡N→∞ρ1(2,N-2;r)<lim⁡N→∞ρ1(3,N-3;r)<lim⁡N→∞ρ1(4,N-4;r)<⋯.


## 4. Conclusions

This paper derives the explicit formulae for thekth-order fixation probabilities on complete star digraphs and applies them to study some asymptotic properties of populations on CSDs. In order to better understand how the structure of a population affects its evolution, we also chose to study how complete bipartite digraphs affect the first-order fixation probability of an advantageous mutant. For fixed population structure, a method (by solving a linear system) has been proposed to calculate this fixation probability. Numerical experiments have revealed some interesting phenomena, which need a reasonable explanation. After careful consideration, we find that, except isothermal digraphs, rooted digraphs and multiplerooted digraphs, complete star digraphs may be the only digraphs whose kth-order fixation probabilities can be expressed analytically. Therefore, a wise strategy for studying the fixation probabilities on other kinds of structures is to numerically determine them [18]. We are currently working toward this direction and also attempting to modify the model under study to adapt the evolution of the trustworthiness of large-scale distributed systems. It is also worth extending this work to, say, the cooperation on one graph, on a pair of graphs, or on hypernetworks [5, 19–23].



---

*Source: 940465-2012-08-13.xml*