Counterintuitive properties of the fixation time in network-structured populations

Evolutionary dynamics on graphs can lead to many interesting and counterintuitive findings. We study the Moran process, a discrete time birth–death process, that describes the invasion of a mutant type into a population of wild-type individuals. Remarkably, the fixation probability of a single mutant is the same on all regular networks. But non-regular networks can increase or decrease the fixation probability. While the time until fixation formally depends on the same transition probabilities as the fixation probabilities, there is no obvious relation between them. For example, an amplifier of selection, which increases the fixation probability and thus decreases the number of mutations needed until one of them is successful, can at the same time slow down the process of fixation. Based on small networks, we show analytically that (i) the time to fixation can decrease when links are removed from the network and (ii) the node providing the best starting conditions in terms of the shortest fixation time depends on the fitness of the mutant. Our results are obtained analytically on small networks, but numerical simulations show that they are qualitatively valid even in much larger populations.


Introduction
Most analytical results for evolutionary dynamics have been obtained for well-mixed populations, regular networks or deme-structures. However, the consideration of non-regular networks has shown that there is a wealth of evolutionary phenomena that is not captured by these approaches. For example, while the fixation probability of a single mutant is the same on all regular networks [1], some non-regular networks can increase this probability and serve as amplifiers of selection, or decrease it and serve as suppressors of selection. It seems to be tempting to use amplifiers of selection to speed up experimental evolution, as the probability of the fixation of beneficial mutants could be increased. However, it turns out that amplifiers of selection can at the same time slow down the time until fixation. Thus, such an approach would have the drawback that while the time until a mutant appears is decreased, the time until it takes over is increased [2,3].
Evolutionary dynamics on network structures can also serve as a powerful abstraction when studying the somatic evolution of cancer [4][5][6]. In this context, the idea is that directed networks can substantially decrease the number of cells at risk and thus inhibit the accumulation of cancer mutations.
While most previous work on network structured populations has focused on the fixation probability [1,4,[7][8][9][10][11], the time to fixation has received considerably less attention so far [2,3,12]. Several questions that appear somewhat obvious are still open: does the time to fixation always increase when a link is removed from a network? Do amplifiers of selection always change the time to fixation?
We study constant selection, where the fitness does not depend on the frequencies of the types. There exists a huge body of closely related research on evolutionary game theory on graphs (e.g. [12][13][14][15][16][17][18][19][20][21]). It has been shown there that the graph structure can substantially affect evolutionary game dynamics. But as the fixation time on networks already leads to interesting and counterintuitive results for constant selection, we focus on this case.

& 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution
As the number of possible states in non-regular networks rapidly increases with population size N, we focus on the smallest population size that allows us to obtain any interesting results, N ¼ 4. We consider all six possible undirected networks in detail and calculate (i) the probability of fixation, (ii) the time to fixation, and (iii) the sojourn times depending on the fitness of the mutant. This approach illustrates that the time to fixation can actually increase when a link is added to the network. It particularly allows us to infer in which states the system spends this additional time. Moreover, our approach shows that the optimal starting point for a novel mutation in terms of its fixation time also depends on the fitness of this mutation.

The Moran process in well-mixed populations
The Moran process is a birth-death process in a well-mixed population [22]. We start from N 2 1 wild-type individuals with fitness 1 and one mutant with fitness r . 0. If r . 1, the mutant is advantageous compared to the wild-type, whereas r , 1 implies a disadvantage. One can also study neutral evolution, where r ¼ 1. At each time step, one of the N individuals is selected for birth with probability proportional to its fitness. This individual gives birth to an identical offspring which replaces another randomly chosen individual.
The probabilities of increasing and decreasing the number of mutants are given by a tridiagonal transition matrix T (Nþ1)Â(Nþ1) . The probability to increase the number of mutants by one is given by for 0 i N. The probability to decrease the number of mutants from i to i 2 1 in one time step is We assume a mutation rate of zero, which is reflected by T 0,1 ¼ T N,N21 ¼ 0. As T 0,21 ¼ T N,Nþ1 ¼ 0, the boundaries i ¼ 0 and i ¼ N are absorbing.
If a mutant replaces another mutant or a wild-type replaces another wild-type, the population does not change. This happens with probability T As the ratio of the transition probabilities is for 1 i N 2 1, the fixation probability f i,N for i mutants in a well-mixed population is given by [23][24][25] The time which one mutant needs to take over the whole population, given that it does succeed, is called conditional fixation time because it is conditioned on the success of the mutants. By contrast, the unconditional fixation time measures how long it takes for the mutants to either go extinct or fixate, starting from one mutant. For birth-death processes (and the Moran process in particular), the expected conditional fixation time t 1,N is given by [23][24][25] T m,mÀ1 T m,mþ1 : (1:5)

The Moran process in structured populations
Population structure can be introduced into the Moran model by assuming that individuals are represented by the nodes in a network and assuming that reproduction and replacement is only possible between connected nodes (see [1] and [24, ch. 8]. The fixation probability can then be assessed in the following way: individuals of the wild-type with fitness 1 are placed on the nodes of a network. A single mutant with fitness r is placed on one of the nodes at random. In each time step, one individual is chosen for birth with probability proportional to its fitness. Then one of its neighbours is chosen at random to be replaced by the new offspring. Thus, the links of a node determine into which of the neighbouring sites the individual on that node can reproduce. Throughout our work, we only consider connected undirected networks with equal weights. The standard Moran process corresponds to the complete network, where every node is adjacent to all other nodes, which implies that the probability of being replaced is equal for all individuals. On isothermal networks, where every node has the same number of neighbours in the case of undirected networks, it has been shown in [1] that the fixation probability is the same as in the well-mixed population. Examples for isothermal networks are the ring and the two-dimensional lattice with periodic boundary conditions. But in general, the probability of replacement can vary crucially between different nodes, resulting in fixation probabilities that are distinct from the well-mixed population. Non-isothermal networks that increase (decrease) the fixation probability for advantageous mutants and decrease (increase) it for disadvantageous mutants are called amplifiers (suppressors) of selection.

A general approach to calculate probabilities and times of fixation
While the analytical results for birth-death processes can be tailored to some network structures [3,26], this does not always work. For assessing absorption probabilities and times of Markov chains, a more general approach, which is typically used in numerical considerations, described in ch. 11 of [27] can be exploited (see also [28] for an application in another context). Let T sÂs be the transition matrix of a discrete-time Markov chain with at least one absorbing state. We renumber the states such that the t transient states are first and the a absorbing states are last, where t þ a ¼ s is the total number of states (for our problem, we will always consider the case of two types, which implies a ¼ 2). The transition matrix now has the following canonical form: where Q tÂt consists of the transition probabilities between transient states and R tÂa describes the transitions from the transient into the absorbing states. As transitions are not possible from an absorbing state to a transient state, the lower left block is zero. Once absorbed, the process will stay there forever, therefore the lower right block of T is an identity matrix I aÂa . For a starting distribution x 1Âs , the product xT m gives the distribution after exactly m time steps. For m!1, we can recover the fixation probabilities from this. Let us call F ¼ S 1 n¼0 Q n ¼ (I À Q) À1 the fundamental matrix of the Markov chain. The i,jth entry of F, given by rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20140606 is the expected sojourn time in the transient state j, given that the process started in the transient state i [27, ch. 11].
Multiplying F with the transition probabilities to the absorbing states provides the absorption probabilities. The absorption probability in state j after starting in state i, To assess the total time the process spends before absorption, we look at the so-called unconditional fixation time. The unconditional fixation time, which is the expected number of time steps before the process is absorbed in either of the absorbing states, can be calculated from the expected sojourn time. Summing over the ith row of F yields the unconditional fixation time t i , after starting in state i. For N 2 1 transient states, this is given by To calculate the conditional fixation time from the unconditional sojourn times, we use an approach described in [29,30]. If f i,N is the fixation probability to reach state N starting in state i and F i,j is the unconditional sojourn time in state j after starting in state i, then the conditional fixation time of the process on any network, starting in state i, is given by (1:10) This allows us to go beyond unconditional times even on heterogeneous networks.

Analytical results for small networks
First, we analyse small networks that allow a fully analytical approach. There are six different connected undirected graphs with four nodes, illustrated in figure 1.
Out of these six networks in figure 1, only the complete graph and the ring are isothermal. On the other four networks, the nodes vary in degree, i.e. number of neighbours. The diamond has two nodes with degree two and two nodes with degree three. On the shovel, there are three types of nodes: having either one, two or three neighbours. The least possible number of links in a connected network of size four is three links. The two associated networks are called the line and the star. The two outer nodes of the line have one neighbour and the two inner nodes have two neighbours. On the star, the nodes vary even more in degree: the centre node has three neighbours, whereas the three leaf nodes are only connected to the centre.

States and transition matrices
To calculate the fixation probability, we first look at the different possible states and the transitions between states. Then we rearrange the states in the transition matrix as discussed in §1.3, such that the transient states are first.
Complete graph and ring of size four. Let I, II, III and IV be the states with 1, 2, 3 and 4 mutants, respectively, and V the state with only wild-type individuals. The states of this Markov chain are shown in figure 2. Transient state numbers are highlighted in blue, whereas absorbing states are shaded in grey.
In figure 2, the five different states of the Moran process on a well-mixed population and a ring of size four are displayed. State I is the initial state. If the absorbing state IV is reached, this means that the mutant reached fixation in the population. Owing to the special structure of the ring, mutants can only invade in a cluster. Therefore, a change in the number of mutants can only happen if nodes at the boundary of the mutant cluster are selected for birth.
The canonical form of the transition matrix for the process on the complete graph is given in equation (2.1). Transient states are highlighted in blue and absorbing states in light grey. (2.1) The diagonal of the transition matrix T mix is positive. The Moran process stays in the same state, meaning that the number of mutants does not change whenever a mutant replaces a mutant or a wild-type individual replaces another wild-type individual.
With the approach given in §1.3, we reproduce the wellknown fixation probability of a mutant in the well-mixed population [23,24] For the ring, we obtain the canonical transition matrix Recall that the ring is isothermal. Hence, the fixation probability must be the same as for the well-mixed population and given by equation (2.2) as well. Indeed, the However, the second line of T ring in equation (2.3) is different from the second line of T mix in equation (2.1). For example, whenever there are two mutants on the ring network, the probability to stay with two mutants in the next time step is T II,II ¼ 1/2, whereas the corresponding probability to stay at the state with two mutants on the complete graph is only T II,II ¼ 1/3. This indicates that fixation takes longer on the ring, because the probability of no change in state II is higher [3].
Diamond. Next, let us consider the diamond, which has nine states I, II,. . ., IX shown in figure 3. There are several possible mutant -node configurations in this network, such that the states are no longer just determined by the number of mutants. Instead, they are also determined by the degree of the respective node.
Owing to the larger state space, the process on the diamond is not a simple birth-death process. Thus, the transition matrix does not have a tridiagonal shape and the conditional fixation time is not given by equation (1.5), but has to be computed from equation (1.10). The canonical form of the transition matrix for the diamond is given by The state space and transition matrix of the other three networks of size four (shovel, line and star, see figure 1) can be analysed in a similar fashion, but a numerical approach is often more efficient to implement. In §2.1.2, we use the transition matrices to compute the fixation probabilities on these networks.

Fixation probabilities
The fixation probability on the six different networks of size four is shown in figure 4. It can be seen that, compared with the well-mixed population, the star increases fixation probability for advantageous mutants and decreases it for disadvantageous mutants. Thus, the star is called an amplifier of selection. As the ring is isothermal, i.e. every node has the same number of neighbours, a mutant has the same fixation probability on the ring as in the well-mixed population [1]. From figure 4, we can conclude that the diamond, the shovel, the line and the star are amplifiers of selection. Thus for size four, all non-isothermal networks are amplifiers, which means that there are no suppressors of selection. Calling a network an amplifier of selection seems to imply that evolution proceeds faster on this network than on the complete graph. This arises from the idea that the waiting time for a successful mutation to occur is much longer than the fixation time which the mutation needs to spread through the population. In the following, we focus on the fixation time in order to see how it is affected by the removal and addition of links.

Fixation times
With the approach from §1.3, we can calculate the expected conditional fixation time for the well-mixed population depending on the fitness r of the mutants, The fixation time t diamond 1,N can be calculated in the same fashion, but it is a rational function with both numerator and denominator of degree 15 with up to 13-digit coefficients. Therefore, only the Taylor approximation for weak selection up to second order is included here, t diamond 1,N % 10:7 þ 0:4(r À 1) À 2:5(r À 1) 2 : (2:6) The structures shovel, line and star have a substantially higher conditional fixation time than the complete graph, the ring and the diamond. For neutral evolution, these are t shovel 1,N (1) % 16, t line 1,N (1) % 21:3 and t star 1,N (1) % 23:2. This shows that on every network of size four, neutral fixation is slower than on the complete network.
Let us now focus on the complete network, the ring and the diamond to analyse the effect of the removal and addition of one link. We first compare the analytical results to simulations. Figure 5 shows the expected conditional fixation time of one advantageous mutant in a population of size four. The analytical result is plotted together with the averages over 10 6 independent realizations. Figure 5 shows that on the diamond and on the ring, fixation time is higher than on the complete graph. On the ring, the increase is approximately one time step, whereas the diamond increases fixation time even more. Thus, figure 5 reveals a surprising aspect of the fixation times: intuitively, one may speculate that the removal of a link should always prolong the process of fixation, because fewer possible paths for the mutant to spread should slow down the process. This intuition can be illustrated by road networks, where one would think that more connections speed up overall traffic. But this is not true for our case and even for road networks, there are paradoxical situations, where the traffic flow can be increased by closing a road [31]. From equations (2.5) and (2.6), we see that although the ring is constructed by dropping one link from the diamond, fixation is faster on the ring than on the diamond. This is also seen in figure 5, where the fixation time is plotted against the mutant's fitness r. This result seems counterintuitive and needs a closer investigation, which can be provided by analysing the sojourn times.

Sojourn times
Intuitively, the more links a network has, the faster the fixation of mutants should happen. But this is not true in general, as shown in figure 5. An additional link may not only provide more possibilities for the mutants to place their offspring; it may also delay fixation when the mutants replace each other more often. We look at the sojourn times to understand this process in more detail.
The sojourn times measure how much time is spent in the transient states on average before one of the absorbing states  Figure 3. Owing to the different degrees of the nodes, the diamond has more possible states than the complete graph. There are two nodes with degree three and two nodes with degree two. Therefore, one must distinguish between those two types of nodes. For example, this leads to the distinction between states I and II, which would be the same if all nodes had the same degree.  of the system is reached [32]. Omitting the cases in which the mutants become extinct, the conditional sojourn times measure the expected time spent in the transient states before mutant fixation. Summing over the conditional sojourn times before absorption into the all-mutant state yields the total time it takes to go from one to N mutants, the conditional fixation time. By comparing the conditional sojourn times in the transient states of different networks, we can infer which states cause the delayed fixation.
In figure 6, we plot the expected time the system spends in states with one, two and three mutants for the wellmixed population, the ring and the diamond. For neutral evolution, r ¼ 1, the well-mixed population and the ring have exactly the same sojourn time in the states with one and three mutants. The ring prolongs the sojourn in the state with two mutants for one time step on average. The diamond has a shorter sojourn time in the two-mutant states than the ring; however, it increases the sojourn times in all states compared with the complete graph.
In figure 6b, we consider the limit of strong selection, r ! 1. The process on the ring stays in the states with one and three mutants for the same amount of time as the process on the complete graph. So the only difference is the sojourn in the two-mutant state, which causes the conditional fixation time to be higher on the ring than in the well-mixed population. We saw in figure 5 that a mutant on the diamond has a higher conditional fixation time than on the ring. As we can see in figure 6b, this increase in fixation time is mainly due to the increased sojourn time in the two different states with three mutants. The three states with two mutants on the diamond have a lower sojourn time than the two-mutant state on the ring. However, the average sojourn in states with three mutants is almost one time step longer. Interestingly, under strong selection more time is spent in these two states than under neutral selection. This is because the probability to leave states VI and VII decreases with higher intensity of selection.
So far, we only know that the total sojourn time in all of the three-mutant states is increased on the diamond. Let us now address the question as to which of the two different states with three mutants causes the prolonged sojourn time, VI or VII. For this purpose, we visualize the sojourn time as a function of the mutant's fitness for all seven transient states.
Looking at the sojourn time in all transient states, see figure 7, we see that for neutral evolution and slight fitness difference, the system spends most time in state IV. This changes at r % 1.65. For higher selective advantage, most time is spent in state VII (where the one wild-type individual is situated at a node with two neighbours, see figure 3). This can be explained in the following way: in state VII, the one wild-type individual has lower chances of being replaced, because it is only connected to two of the three mutants. The three mutants keep replacing each other for a longer time than in state VI, therefore the process stays longer in this state before going to fixation.     Instead of randomly choosing a node to place the first mutant on, we now assess the effect of the initial node on the fixation time. On the complete graph and the ring, all nodes have the same number of links and therefore all initial conditions are identical. However, the diamond consists of two nodes with three neighbours and two nodes with two neighbours. The conditional fixation time for those two initial conditions is plotted in figure 8. The average of the two curves is identical to the conditional fixation time shown in figure 5, because the probability to place the first mutant at either one of the initial states is 1/2.
In figure 8, the initial position of the mutant has a non-trivial impact on the fixation time. For small values of r, starting from a node with degree 2 takes longer to fixate than starting from a node with degree 3. This changes at r % 5.8. For larger values of r, starting at a node with degree 3 leads to a slightly higher fixation time than starting at a node with degree 2.

Numerical simulations for larger networks
So far, we have only considered very small networks in order to allow for an analytical consideration of the fixation times. Thus, we still have to show that this is not only an effect of these extremely small systems. To explore whether the removal of links can also decrease the fixation time, we use our small networks as motifs of a larger network that is constructed from these motifs (see figure 9). As several links are removed simultaneously, the effect size is not expected to decrease rapidly with the system size. For these larger networks, the analytical approach is not feasible and we therefore perform simulations instead. We start the birthdeath Moran process by putting one mutant on a random node of a network consisting solely of wild-type individuals. At each time step, one individual is chosen for birth with probability proportional to its fitness. The offspring then replaces one of its neighbours at random. The fixation time of the mutants is averaged over 10 5 independent realizations. Figure 9 shows the mean conditional fixation time on quadratic lattices with and without periodic boundary conditions for mutant fitness r ¼ 2. Fixation occurs much faster with periodic boundary conditions than on the corresponding networks without periodic boundary conditions. Note that by connecting the boundaries, we make the lattices isothermal, i.e. all nodes have the same number of neighbours. Hence, the fixation probability on these networks is the same as on the complete network of the same size. However, the fixation time is increased, compared with the well-mixed population on a complete network.
Intuitively, the fewer links a graph has, the longer fixation should take. And as figure 9 shows, the lattices with periodic boundary conditions confirm that intuition. Meaning that starting from network F, by dropping links to obtain network E, the fixation time increases. By dropping even more links (network D), the fixation time increases even more. A similar result was found by Whigham & Dick [33]. They simulated the process on different isothermal ring structures of size N ¼ 100 and showed that the fixation time decreases with increasing node degree. But this ordering property by the number of links seems to only hold for isothermal graphs.
Interestingly, the lattices without periodic boundary conditions do not behave according to that intuition. Instead, fixation takes longer on network A than on network B, even though A has more links. This shows that the results obtained for size N ¼ 4 are qualitatively still valid for larger lattices without periodic boundary conditions.

Discussion
Evolutionary dynamics on networks is a way to consider population structures that can lead to many non-trivial findings. For example, one can construct amplifiers or suppressors of selection which lead to fixation probabilities that deviate from results in a systematic way. The isothermal  theorem states that for the birth -death Moran process, all regular networks have the same fixation probability, which is a remarkable finding [1]. Also temporal aspects of this dynamics can be highly interesting: a population structure that leads to higher probabilities of fixation can at the same time increase the time of the fixation process itself. Based on the results of Frean et al. [3], we originally set out to show that the fixation time increases on undirected networks if links are removed. However, as we have shown here, it is possible to construct counter-examples in which the removal of links decreases the fixation time.
A similar phenomenon has been reported in transportation systems, where the flow through a system can be increased if connections are removed [34]. The Braess Paradox describes the situation where adding a seemingly helpful link can have a negative effect on the flow through the network [31].
Our analytical approach for small networks allows us to infer how much time the system spends in which states, meaning that the fixation time can be dissected in great detail.
By analysing the sojourn times, we have shown why our intuition of decreasing fixation time by adding links is not true in general. As we have seen, adding links to a network not only increases the possibilities for invasion of the mutants, but also increases the likelihood of mutants replacing each other. Therefore, a general statement about the behaviour of the fixation time when removing or adding links is not possible.
Furthermore, we have shown that the starting position of the first mutant has a crucial impact on the fixation time. Depending on the fitness of the mutant, it can be faster if the mutant is placed at a highly connected node or at a node that has only very limited connectivity. On the diamond, a disadvantageous or slightly advantageous mutant does better in terms of fixation time when starting at the more highly connected nodes. However, if the mutant is very advantageous, the more isolated nodes provide a shorter fixation time.
To investigate the counterintuitive result obtained for size N ¼ 4, the small networks were used to create larger quadratic lattices. We have shown that this decrease of the fixation time with the removal of links is not just a finite size effect, but can also be found for large population sizes.
In particular, heterogeneous networks seem to be most relevant for the real world. They are found among humans and different animal species [35][36][37][38][39][40][41]. However, to directly transfer our results to such systems seems premature, as it is unclear whether the dynamics in our model is a good approximation for the processes in these networks. Moreover, our analysis reveals that there are still open challenges for a full theoretical understanding of evolutionary dynamics on networks. Unfortunately, the size of the state space increases rapidly with the system size and requires a tedious construction of transition matrices. While analytical approaches are still being developed [42][43][44], exploring these networks will largely rely on numerical simulations.