Phase transitions in evolutionary dynamics

The evolutionary dynamics of a finite population where resident individuals are replaced by mutant ones depends on its spatial structure. The population adopts the form of an undirected graph where the place occupied by each individual is represented by a node and it is bidirectionally linked to the places that can be occupied by its clonal offspring. There are undirected graph structures that act as amplifiers of selection increasing the probability that the offspring of an advantageous mutant spreads through the graph reaching any node. But there also are undirected graph structures acting as suppressors of selection where, on the contrary, the fixation probability of an advantageous mutant is less than that of the same individual placed in a homogeneous population. Here, we show that some undirected graphs exhibit phase transitions between both evolutionary regimes when the mutant fitness varies. Firstly, as was already observed for small order random graphs, we show that most graphs of order 10 or less are amplifiers of selection or suppressors that become amplifiers from a unique transition phase. Secondly, we give examples of amplifiers of order 7 that become suppressors from some critical value. For graphs of order 6 and 7, we apply computer aided techniques to exactly determine their fixation probabilities and then their evolutionary regimes, as well as the exact critical values for which each graph changes its regime. A similar technique is used to explore a general method to suppress selection in bigger orders, namely from 8 to 15, up to some critical fitness value. The analysis of all graphs from order 8 to order 10 reveals a complex and rich evolutionary dynamics, which have not been examined in detail until now, and poses some new challenges in computing fixation probabilities and times of evolutionary graphs.


Introduction
In recent times the evolutionary theory on graphs has become a key field to understand biological systems. Although evolutionary dynamics has been classically studied for homogeneous populations, there is now a wide interest in the evolution of populations arranged on graphs after mutant spread. The process transforming nodes occupied by residents into nodes occupied by mutants is described by the Moran model. Introduced by Moran [1] as the Markov chain that counts the number of invading mutants in a homogeneous population, it was adapted to subdivided population by Maruyama [2,3] and rediscovered by Lieberman et al. [4] for general graphs. For undirected graphs where links have no orientation, mutants will either become extinct or take over the whole population, reaching one of the two absorbing states, extinction or fixation. The fixation probability is the fundamental quantity in the stochastic evolutionary analysis of a finite population.
If the population is homogeneous, at the beginning, one single node is chosen at random to be occupied by a mutant individual among a population of N resident individuals. Afterwards, an individual is randomly chosen for reproduction, with probability proportional to its reproductive advantage (1 for residents and r ≥ 1 for mutants), and its clonal offspring replaces another individual chosen at random. In this case, the fixation probability is given by Φ 0 (r) = 1 − r −1 1 − r −N = r N −1 r N −1 + r N −2 + · · · + r + 1 . (1) If the population is arranged on nodes of an undirected graph, the replacements are limited to the nodes that are connected by links. According to the Isothermal Theorem [2][3][4], the fixation probability Φ(r) = Φ 0 (r) if and only if the graph is isothermal (i.e. the temperature T i = j∼i 1/d j of any node i is constant, where j is a neighbor of i and d j is the number of neighbors of j), or equivalently regular (i.e. the degree d i of any node i is constant). But there are graph structures altering substantially the fixation chances of mutant individuals depending on their fitness. As showed in [4], there are graph structures that amplify this advantage. This means the fixation probability function Φ(r) > Φ 0 (r) for all r > 1 for the same order N . Notice that Φ(1) = 1/N and the inequality must be reversed for r < 1. Due to the exact analytical computation of Φ(r) given by Monk et al. [5], it is known that star and complete bipartite graphs are amplifiers of natural selection whose fixation functions are bounded from above by On the other hand, there are also graph structures that suppress the reproductive advantage of mutant individuals so that Φ(r) < Φ 0 (r) for all r > 1. Examples of this kind of graph structures were known for some fitness values (see [6]). In [7], we introduced examples of suppressors of natural selection of order 6, 8 and 10 (see Fig 1), denoted by 6 , 8 and 10 , whose fixation probabilities remain smaller than Φ 0 (r) for every r > 1. , we called -graph the undirected graph of even order N = 2n + 2 ≥ 6 obtained from the clique K 2n by dividing its vertex set into two halves with n ≥ 2 vertices and adding 2 extra vertices. Each of them is connected to one of the halves of K 2n and with the other extra vertex.
The analysis of disadvantageous mutants (with r < 1) is also interesting when comparing amplification and suppression of selection for graphs, but here, for simplicity, we focus on the case of advantageous mutant (with r > 1). Different initialization and updating types have been also considered in [8] and [9], see also [10] for a comparative analysis of update mechanisms. If the initial distribution is uniform (i.e. the probability that a node will be occupied by the initial mutant is equal for all the nodes) and the graph evolves under Birth-Death updating, Hindersin and Traulsen showed in [9] that almost all small undirected graphs are amplifiers of selection. Assuming both conditions and focusing on the advantageous case, we distinguish two different evolutionary regimes out of the isothermal one: given values 1 ≤ r 0 < r 1 ≤ +∞, a graph is an amplifier of selection for r ∈ (r 0 , r 1 ) if the fixation probability function Φ(r) > Φ 0 (r) and a suppressor of selection for r ∈ (r 0 , r 1 ) if Φ(r) < Φ 0 (r) for all r 0 < r < r 1 . Star and bipartite complete graphs are amplifiers for r ∈ (1, +∞) and graphs 6 , 8 and 10 are suppressors for r ∈ (1, +∞). There are suppressors that become amplifiers from some critical value r c > 1 (see Fig 2). Borrowing the terminology from percolation theory, in this case, we say r c is a transition phase from the suppressor regime to the amplifier regime. In general, we say that a number r c > 1 is a transition phase between two different evolutionary regimes if there are values 1 ≤ r 0 < r 1 ≤ +∞ such that the graph is a suppressor (resp. an amplifier) for r ∈ (r 0 , r c ) and an amplifier (resp. a suppressor) for r ∈ (r c , r 1 ).  [11].
In this paper, we show that most graphs of order 10 or less are amplifiers or suppressors that become amplifiers from a unique transition phase r c > 1. However, we also exhibit two other types of transitions which are much more interesting: • There are amplifiers of order greater or equal to 7 that become suppressors from a unique transition phase r c > 1.
• There are graphs of order greater or equal to 8 exhibiting more than one transition phase.
Since the fixation probabilities Φ 0 (r) and Φ(r) are monotonically increasing, these results were not exactly expected, while numerical simulations should be considered with caution. In fact, for finite populations, it has been suggested that results obtained for weak selection may remain valid when the selection is no longer weak. However, in [12], Wu et al. showed that this is no the case for homogeneous populations under frequency dependent selection. Here, we show that the graph structure can completely alter the evolutionary regime of a structured population even under constant selection. In addition, the fact that the survival chances of mutant individuals may decrease with respect to a homogeneous population when their fitness increases might have interesting biological consequences. Our results also demonstrate that numerical simulation is not enough to determine the evolutionary regime of a graph for the whole range of fitness values. Moreover, the very slight difference between different regimes for some graphs often forces to consider an extremely high number of trials when using Monte Carlo method, like recently in [13] to show a family of graphs of order ≥ 12 that changes from amplifier to suppressor as r increases. For example, 10 10 trials were needed in [7] to determine the evolutionary regime of the graphs N for 12 ≤ N ≤ 24 and 1 ≤ r ≤ 4.
In [11], we presented an extremely precise database of fixation probabilities for all undirected graphs of order 10 or less for fitness values r varying from 0.25 to 10 with step size of 0.25 (see details in Materials and Methods). Here, we use this database to detect new evolutionary dynamics: in this way, we identify two amplifiers of order 7 with a transition to the suppressor regime at some r c ≤ 10. Then we adapt the technique described in [7] to symbolically compute the exact fixation probability Φ(r) of these examples for all r > 1 proving that they really have a unique transition phase. A similar procedure is used to exactly determine the evolutionary regime of any graph of order 7 or less for any fitness value. Thus, we find a third example of order 7 exhibiting a transition from the amplifier to the suppressor regime for r c > 10. The first method is also applied to show the existence of graphs of order 8, 9 and 10 with more than one transition phase. This reveals that undirected graphs have a complex and rich evolutionary dynamics that are worth studying in detail.

Mathematical model
Let G be an undirected graph with node set V = {1, . . . , N }. In fact, all graphs considered here will be assumed connected. Denote by d i the degree of the node i. The Moran process on G is a Markov chain X n whose states are the sets of nodes S inhabited by mutant individuals at each time step n. The transition probabilities are obtained from the matrix W = (w ij ) given by w ij = 1/d i if i and j are neighbors and w ij = 0 otherwise. More precisely, for each fitness value r > 0, the transition probability between S and S is given by where is the total reproductive weight of resident and mutant individuals. The fixation probability of each subset is a solution of the system of 2 N linear equations with boundary conditions Φ ∅ (r) = 0 and Φ V (r) = 1. The (average) fixation probability is given by Denoting by P (r) = (P S,S (r)) the transition matrix, Eq 6 can be written as where Φ(r) = (0, Ψ(r), 1) is the vector with coordinates Φ S (r), (1, b(r), 0) is the vector with coordinates P S,∅ (r), and (0, c(r), 1) is the vector with coordinates P S,V (r). It can be also rewritten as where I is the identity matrix of size 2 N − 2. This equation has a unique solution Ψ(r) whose coordinates are rational functions on r with rational coefficients [7, S1 Text]. But considering Eq 3, we can multiply the equation associated to each state S by w S (r) reducing Eq 9 to where the entries of Q * (r) and c * (r) are now degree one polynomials with rational coefficients.

Database
In [11], we presented an accurate database of the fixation probabilities for all connected undirected graphs with 10 or less vertices, which means 11,989,763 graphs excluding the trivial one with one single vertex. The generation of the edge lists was done with SageMath, whereas the computation of Φ(r) was written in the C programming language. Firstly, we compute the matrix Q * (r) and the vector c * (r) whose entries are represented as couples a and b of pairs of 64 bits integers. Next, we solve Eq 10 for each fitness value r varying from 0.25 to 10 with step size of 0.25. Each graph is identified (up to isomorphism) with a unique 64 bits unsigned integer, which allows us to recover the edge list without previous knowledge of its order or size, see again [11] and references therein for details. The database is presented as a single HDF5 file available at [14].

Computation method
A method to compute the exact (average) fixation probability Φ(r) of graphs of small order with symmetries was described in [7]. As we already said, Φ(r) = Φ (r)/Φ (r) is a rational function where the numerator Symmetries are used to reduce the degree d and therefore the number 2(d + 1) of coefficients involved in the computation of Φ(r). Since Φ(r) converges to 1 as r → +∞, we have a d = b d = 1 and then Eq 6 can be replace with a system of 2d linear equations that arise from evaluating Φ(r) for fitness values r ∈ {1, . . . , d + 1, 1/2, . . . , 1/d}. Instead of the SageMath program given in [7], a new C++ program available from S1 File allows us now to symbolically • compute the fixation probability Φ(r) for these values, and • solve the reduced linear system Eq 11.
Once Φ(r) has been calculated solving this system, we can determine the states and the phase transitions of the graph by computing the sign and the zeros of the rational function Φ(r) − Φ 0 (r).

Results
Contrary to a widespread idea that it is easier to construct suppressors of selection than amplifiers of selection, which is true only if the edges are oriented, Hindersin and Traulsen [9] showed that most undirected (Erdös-Rényi) random graphs of small order are amplifiers of selection under Birth-Death updating. Here, by regarding all graphs of order 10 or less, we confirm that most undirected graphs of order ≤ 10 are amplifiers of selection for fitness values r ≤ 10, see Table 1. In order 6, there are exactly one suppressor of selection, namely the graph 6 described in [7] (see Fig 1), five isothermal graphs, six suppressors that become amplifiers from a unique critical value, and the remaining 100 are amplifiers of selection (see S1 Table). Two of suppressors changing into amplifiers was already described in [11] (see Fig 2.(A)). All the graphs portrayed in the paper are gathered in S1  However, additional exact computations for some higher values of the fitness r (varying from 10 to 2, 000 with step size of 1) seem exclude more than two transitions in order 8 and 9 and more than three transitions in order 10. Exact results are marked in bold: a complete description of the phase transitions for graphs of order 6 and 7 is given in S1  that later become amplifiers of selection, namely 52, including the suppressors presented in [6] and depicted in Fig 2.(B). Moreover, we find indeed three amplifiers for weak selection Id 1151592835082, Id 1151860745228, and Id 1151592837126 with transitions at r c ≈ 4.98, r c ≈ 6.37 and r c ≈ 24.79 respectively. A quantitative resume is given in Table 1, while some graphs are depicted in S1 Fig and S2   The two first suppressors Id 1134281908237 and Id 1134281902105 are shown in Fig 4, where we also added a third graph Id 1151998648333 with a unique transition from the suppressor to the amplifier regime at r c ≈ 5.17. Their construction is very similar to -graphs defined in [7]. Recall that N = n,n N is an undirected graph of even order N = 2n + 2 ≥ 6 obtained from the clique K 2n by dividing its vertex set into two halves with n ≥ 2 vertices and adding 2 extra vertices. Each of them is connected to one of the halves of K 2n and with the other extra vertex (see Fig 1). More generally, we denote by n,m N the undirected graph obtained adding two interconnected extra vertices to the clique K N −2 and connecting each one to disjoint families of vertices in the clique having n and m elements with n + m ≤ N . We say that n,m N is balanced if n = m and unbalanced otherwise. As explained in S1 Fig, using

Phase transitions for -graphs
For the three graphs in Fig 4, the suppression mechanism seems directly related to the suppressor nature of -graphs [7] and clique-wheels [15]. Indeed, on a star graph, it is more likely that a peripheral mutant survives and reproduces occupying the central node. On the contrary, on a complete graph, it is much more unlikely that the initial mutant will survive surrounded by residents. But in -graphs and clique-wheels, the survival chances of the mutants placed at the central clique are reduced by its connection with peripheral nodes occupied by residents, as well those of the peripheral mutants connected with central residents, although a subtle balance seems to be needed to suppress the reproductive advantage of mutant individuals according to our data. To confirm this idea, we explore evolutionary regimes and transitions of some balanced and unbalanced -graphs. Unbalanced -graphs, namely 1,4 7 and 2,3 7 , were studied in [7] using Monte Carlo simulation. Both were shown as suppressors; the first one was shown as amplifier from a relatively small fitness value, whereas the second one was shown as suppressor for any fitness value r ≤ 10. Now we know they exhibit a unique transition phase (from the suppressor to the amplifier regime) at r c ≈ 1.80 and r c ≈ 25.47 respectively.
In order 8, we consider the graphs 2,2 8 , 2,3 8 , and 1,4 8 generalizing 2,2 7 , 2,3 7 , and 1,4 7 , which are depicted in Fig 5. They still are suppressors that become amplifiers from critical values r c ≈ 4.15, r c ≈ 5.32 and r c ≈ 1.89 respectively.  Due to the symmetries of the graphs 2,2 N , our computation method is also applicable to these graphs when N varies from 6 to 15. The exact differences Φ 0 (r) − Φ(r) are depicted in Fig 6 when r varies from 1 to 10, although the monotonous behavior of transitions can be better seen in S2 Table. As before, all these graphs are suppressor with a unique transition phase to the amplifier regime. Transitions from amplifier to suppressor regime As we have already said, there are three amplifiers Id 1151592835082, Id 1151860745228, and Id 1151592837126 that become suppressors from critical values r c ≈ 4.98, r c ≈ 6.37 and r c ≈ 24.79. The first and the last of these graphs are constructed from the same building blocks (see Fig 7). Thus, we say that Id 1151592835082 is a friendship cycle F C 2 7 and Id 1151592837126 is a friendship star F S 4 7 . A friendship cycle F C m N is a cycle C n where m disjoint pairs of neighbors are connected to an extra central vertex so that two vertices have at most one neighbor in common, the total order N being equal to n + 1. In a friendship star F S m N , two neighbors can share two different neighbors, but they cannot be connected either to each other, nor to another new neighbor. Additionally, one single vertex can belong to more than two 3-cycles composing a friendship subgraph (to distinguish it from a friendship ribbon F R m N where no vertex can belong to more than two 3-cliques, see S2 Fig). Here N is the number of vertices and m is the number of 3-cliques.
We symbolically compute the fixation probability Φ(r) for the friendship cycles F C 2 N with N = 7, 8, 9 proving that they are amplifiers transformed into suppressors from fitness values r c ≈ 4.98, r c ≈ 3.12, and r c ≈ 2.45 respectively (see Fig 8-(A)). For bigger orders N varying from 10 to 15, we compute the fixation probabilities for selected fitness values r between 1 and 10. The differences Φ 0 (r) − Φ(r) are shown in Fig 8-(B). Only F C 2 10 shows a transition in that interval and it is not clear that the others friendship cycles F C 2 N evolve from the amplifier to the suppressor regime. However, the friendship cycle F C 3 10 (see S2 Fig) also is an amplifier that becomes a suppressor from r c ≈ 4.09.
On the other hand, we also symbolically computed the Φ(r) for the friendship star F S 6 10 showing that it becomes a suppressor at r c ≈ 9.96. Since there are not enough symmetries to reduce the system of equations, this approach could not be pushed any further. The same happens with F R 6 10 , and therefore we have to use the method from [11] seeing that this is an amplifier for r ≤ 10. Since we extended the interval until r = 2000 without finding any transition, this is probably a global amplifier. These graphs are also depicted in S2

Amplifiers and suppressors of order 8 and more
Other interesting phenomena are also revealed by the barcode diagrams in order 8 and higher as shown in S3 Fig, S4 Fig, and S5 Fig. Now, the number of amplifiers and suppressors is only approximate since the fitness is limited to the interval [1,10]. Thus, we distinguished between suppressors and apparent suppressors, that is, those being suppressors for 1 < r ≤ 10. The existence of double transitions of type Suppressor/Amplifier/Suppressor is clearly visible in S3 Fig As we said before, these remarks are only apparent because we cannot be sure that there will not exist other transitions for these orders. However, supplementary exact computations for higher values of the fitness r (varying from 10 to 2, 000 with step size of 1) corroborate the idea that there is not new transitions in these few graphs. Subject to this caveat, Table 1 gives a detailed account of isothermal and suppressor graphs, as well as graphs exhibiting more than one phase transition. This clarifies the assertion we made above that most graphs of order 10 or less are amplifiers of selection (cf. [9]) or suppressors that become amplifiers from a unique transition phase. The differences computed for 10 ≤ N ≤ 15 and 1 ≤ r ≤ 10 by combining the method described in [11] with the use of symmetries.

Discussion
Initially motivated by our interest in the robustness of biological and technological networks against invasion [16], we found in [7] some graph structures suppressing the advantage of mutant individuals occupying their nodes for any fitness value. This seems particularly appealing for biological networks like brain and protein-protein interaction networks, but also in the tumor initiation process within healthy tissues as proposed in [17]. Most graph structures reduce in a very slight amount the advantage of a invading mutant, but some suppression mechanisms could be amplified by repetitive rules (such as those described in [18] and [19] for neuronal networks) involved in the modular architecture of many biological networks.
As explained in Matherials and Methods, we had previously computed the fixation probabilities for all nontrivial connected undirected graphs with 10 or less vertices, totaling 11, 989, 763 graphs, for fitness values varying from 0.25 to 10 with step size of 0.25. Collected data have been published in [14]. Here, we used these data to show that most graphs of order 10 or less are amplifiers of selection. In spite of the complexity of some amplifier structures [4] and the effort required to prove that nature [5,6,8,20,21], this confirms the observation previously made by Hindersin and Traulsen [9] for random graphs of small order. It was also known that some suppressors become amplifiers from some critical fitness value, but here we exhibited new examples (gathered in S1 Fig) whose suppression mechanism is similar to that of clique-wheels studied in [15] and -graphs exhibited in [7]. However, for graphs of order N = 7, there is another type of transition: we found two amplifiers that become suppressors from critical values r c ≈ 4.98 and r c ≈ 6.37. But for larger orders, the phase transitions are even more amazing because some graphs exhibit more than one transition. As before, for N = 8, these graphs were detected from the barcodes diagram S3 Fig and then identified and represented (see S6 Fig) by using the database [14]. All double transitions are identical of type Suppressor/Amplifier/Suppressor.
Contrary to the idea that results obtained for weak selection may remain valid out of this case (see [12] and references therein for a discussion about this problem), these observations indicate that the evolutionary regime of a structured population can be completely altered by the graph structure even under constant selection. In our opinion, this has important biological and theoretical implications.
As an analytical computation of the fixation probability for all these graphs does not seem feasible, we adapted the method described in [7] (included in Materials and Methods and whose new code is available from S1 File) to symbolically compute the fixation probabilities of small order graphs with enough symmetries. This allowed us to determine the evolutionary regime of each graph of order 7 or less, providing the exact part of Table 1, which is one of the main results of the paper. The exact value of each transition is included in S1 Table. Thus we found a third amplifier of order 7 with a transition phase at r c ≈ 24.79, and consequently we were interested in the evolutionary regime of friendship cycles and friendship stars portrayed in S2  S5 Fig for order N = 10. In the first case, there are 6 graphs with double transitions of type Amplifier/Suppressor/Amplifier among a total of 49 graphs with more than one transition. In the second one, we found 19 graphs with triple transitions, all of the same type Suppressor/Amplifier/Suppressor/Amplifier. All these remarks are part of Table 1. Given the number of graphs of order 9 and 10 having some transition, it would be tedious (and hard without enough symmetries) to determine the exact number of transitions for each graph, so it cannot be excluded the (unlikely) existence of new transitions. However, additional computations for higher values of the fitness r (varying from 10 to 2, 000 with step size of 1) when the graph has a triple transition seem exclude this possibility.
As we already said, although transitions between different evolutionary regimes were known from [6], [9], [11] and [13], theses results reveal that undirected graphs have a complex and rich evolutionary dynamics, and pose new challenges in computing fixation probabilities and times. Now we know that numerical simulation cannot provide accurate answers to extremal problems on fixation probabilities, nor probably on absorption or fixation times. Our techniques can be easily adapted to determine these times for graphs of small order with enough symmetries. Analyzing the case of disadvantageous mutants (with r < 1) is feasible with the same techniques. Also, study the influence of the initial placement of the mutant on the evolutionary outcome for both uniform and temperature initialization (which perhaps is biologically more plausible [8]). In fact, this could clarify better why some suppression mechanisms work effectively. More difficult, however, is to change the updating method, Birth-Death by Death-Birth, converting amplifiers into suppressors and vice versa.
Finally, all barcode diagrams, but specially those of graphs of order 10, show very particular patterns in the distribution of phase transitions. The existence of multiple transitions has been established, but if we think that the fixation probability is a monotonically increasing rational function, this fact is somewhat surprising. What are the features of the graphs leading to multiple phase transitions and their particular distribution is an important issue for future work. Though the real challenge is to translate differential features of suppressors (as has been achieved with some isothermal graphs and amplifiers) into features of their state spaces allowing to analytically compute their fixation probabilities and then understanding their suppression mechanisms.
Supporting Information S1 File. C++ program. To compute the fixation probabilities of small order graphs for any fitness value r > 1, we ran a new C++ program available at https://bitbucket.org/geodynapp/phasetransition. S1 Table. Exact phase transitions from order 2 to order 7. Evolutionary state of all connected graphs of order 7 or less whose fixation probability Φ(r) = Φ 0 (r) for some fitness value r > 1.
S2 Table. Exact phase transitions for graphs 2,2 N from order 6 to order 15.