(C) PLOS One This story was originally published by PLOS One and is unaltered. . . . . . . . . . . Evolutionary graph theory beyond pairwise interactions: Higher-order network motifs shape times to fixation in structured populations [1] ['Yang Ping Kuo', 'Computational Biology Department', 'School Of Computer Science', 'Carnegie Mellon University', 'Pittsburgh', 'Pennsylvania', 'United States Of America', 'Joint Carnegie Mellon University-University Of Pittsburgh Ph.D. Program In Computational Biology', 'Oana Carja'] Date: 2024-04 Abstract To design population topologies that can accelerate rates of solution discovery in directed evolution problems or for evolutionary optimization applications, we must first systematically understand how population structure shapes evolutionary outcome. Using the mathematical formalism of evolutionary graph theory, recent studies have shown how to topologically build networks of population interaction that increase probabilities of fixation of beneficial mutations, at the expense, however, of longer fixation times, which can slow down rates of evolution, under elevated mutation rate. Here we find that moving beyond dyadic interactions in population graphs is fundamental to explain the trade-offs between probabilities and times to fixation of new mutants in the population. We show that higher-order motifs, and in particular three-node structures, allow the tuning of times to fixation, without changes in probabilities of fixation. This gives a near-continuous control over achieving solutions that allow for a wide range of times to fixation. We apply our algorithms and analytic results to two evolutionary optimization problems and show that the rate of solution discovery can be tuned near continuously by adjusting the higher-order topology of the population. We show that the effects of population structure on the rate of evolution critically depend on the optimization landscape and find that decelerators, with longer times to fixation of new mutants, are able to reach the optimal solutions faster than accelerators in complex solution spaces. Our results highlight that no one population topology fits all optimization applications, and we provide analytic and computational tools that allow for the design of networks suitable for each specific task. Author summary Accelerating the rate of solution discovery has been a long-standing goal in directed evolution and evolutionary optimization problems. One way to accelerate evolutionary search is to design population structures that control the mode and tempo of evolutionary dynamics. Population structure can amplify the action of selection and accelerate rates of evolution, or reversely, suppress selection and slow evolution down. Here we use evolutionary graph theory to systematically study how higher-order motifs in a population’s network of interaction and replacement can shape the velocity of evolution. In particular, we show that by changing the number of triangles (motifs of degree three) in a network, while fixing lower-order structures, we can continuously tune times to fixation of new variants in the population, independently of their probabilities of fixation. This allows us to design population structures that find optimal solutions across a broad spectrum of evolutionary optimization problems. We show that the best solutions need not always be the ones with the smallest times to fixation and discuss what motifs found in real biological networks can teach us about designing highly evolvable artificial ones. Citation: Kuo YP, Carja O (2024) Evolutionary graph theory beyond pairwise interactions: Higher-order network motifs shape times to fixation in structured populations. PLoS Comput Biol 20(3): e1011905. https://doi.org/10.1371/journal.pcbi.1011905 Editor: Christian Hilbe, Max Planck Institute for Evolutionary Biology: Max-Planck-Institut fur Evolutionsbiologie, GERMANY Received: July 10, 2023; Accepted: February 12, 2024; Published: March 15, 2024 Copyright: © 2024 Kuo, Carja. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All relevant data are within the manuscript and its Supporting information files. In addition, simulation code is available to be published with the paper at the following link: https://github.com/yangpingkuo/Evolutionary-graph-theory-beyond-pairwise-interactions. Funding: We gratefully acknowledge support from the NIH National Institute of General Medical Sciences (award no. R35GM147445 to OC), the United States-Israel Binational Science Foundation (award no. 2019266 to OC) and from the NIH T32 training grant (no. T32 EB009403 to YPK). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist. Introduction The spatial structure of a population is a powerful determinant of a population’s evolutionary outcome. Some structures have topological properties that can speed up evolution and amplify the spread of a mutant with even the slightest selective advantage [1], while others work against the force of selection, increasing the role of evolutionary stochasticity and chance [2]. Evolutionary graph theory is the mathematical framework that formalizes the representation of complex population structure and its evolutionary effects [2–5]. In this framework, each individual occupies a node in a graph, and the edges represent the spatial or replacement patterns of interaction between neighboring nodes. The mode and tempo of evolution are studied through two main quantities: the probability of fixation of a new mutant in the population, which measures the likelihood that the mutant lineage takes over the entire population, and the mean time to fixation, the expected time until the population consists of only mutant descendants. Most prior theoretical work has focused on studying probabilities of fixation, with theoretical explorations of times to fixation restricted to small networks [6], symmetric topologies such as lattices, rings, and stars [7–10] or sparse networks, where mutants grow in clusters [11, 12]. This imbalance of focus can be partially attributed to the assumption that the waiting time until a successful mutant appears in the population is much larger than the time it takes for this mutant to sweep through it. This would make the rate of evolution of the population mostly depend on the mutant’s probability of fixation. However, this assumption does not always hold, especially for engineering applications in directed evolution or evolutionary optimization problems with the goal of evolutionarily increasing rates of solution discovery [13–19]. These applications often utilize an elevated rate of mutation in order to speed up the generation of new candidate solutions and this accelerates search time by several orders of magnitude [20, 21]. In these contexts, fixation times of new mutants start to shape rates of evolution even more than probabilities of fixation, and selecting population structures purely for amplification or suppresion of selection could lead to a substantial evolutionary slowdown [11]. This has lead to recent interest in studying times to fixation, especially in the context of trade-offs with the probability of fixation [22, 23]. These studies have explored how to topologically build population graphs that increase probabilities of fixation of beneficial mutations, at the expense, however, of longer fixation times, which can slow down rates of evolution under elevated mutation rate. It remains an open question how to change a network’s topology in order to optimize times to fixation of new variants in the population, with negligible change to probabilities of fixation, i.e. to the amplification of selection. While previous work has focused on the evolutionary role of pairwise node connections exclusively, many spatial patterns of replacement do not take place between pairs of nodes, but rather as collectives, at the level of groups of nodes [24, 25]. Networks exhibit higher-dimensional patterns of node interconnections, which can be organized by the number of nodes participating in forming the patterns (Fig 1A). Lower-order patterns of connectivity, that can be captured at the level of individual nodes (dimension d = 1) and edges (dimension d = 2), have been shown to significantly shape probabilities of fixation [5, 26–29] and are therefore unsuitable for shaping times to fixation, while keeping probabilities constant. Here we show that moving beyond dyadic structures is necessary to explain the trade-offs between probabilities and times to fixation of new mutations and to be able to understand the networks for which we can minimize or maximize times to fixation, without changing the probabilities of fixation. PPT PowerPoint slide PNG larger image TIFF original image Download: Fig 1. Illustration of the model. Panel A illustrates the levels of structural organization in the network. The list of motifs for d ≥ 4 is not exhaustive. Panel B illustrates the Bd (Birth-death) and the dB (death-Birth) update rules. Panel C shows the degree heterogeneous graphs we design, consisting of two groups of nodes with distinct degrees. Panel D illustrates the edge swap operation used to tune triangle fractions in the graph, without changing the degree distribution and mixing pattern of the network. Initially, there is no triangle consisting of mixed node degrees. We randomly select two edges of the same type (denoted by color) to be disconnected and nodes that were “parallel” with respect to the two disconnected edges are then connected, thus preserving the number of edges. After the rewiring step, there exists a triangle that connects two yellow nodes and a blue node. The degrees of the nodes and frequencies of edge type, however, are preserved. https://doi.org/10.1371/journal.pcbi.1011905.g001 We systematically explore how higher-dimensional network motifs (d ≥ 3), and, in particular, three-dimensional wedges and triangles, shape probabilities and times to fixation of a new mutant in the population and identify motifs that allow for continuous tuning of times to fixation, independent of fixation probabilities. We show that increasing the triangle count of a graph increases times to fixation, without influencing the probability of fixation. We also show that increasing the mean degree of the network not only decreases the time to fixation, but also diminishes the triangles’ ability to shape times to fixation. We find a weak increase with triangle count for the probability of fixation only for highly assortative graphs, but, since they are known suppressors of selection, this effect does not undermine the utility of tuning triangles counts for these networks [5]. This is because in applications where we want to reduce the rate of evolution, increasing time to fixation can achieve the same goal as decreasing the probability of fixation. Collectively, our results suggest that the fraction of triangle motifs in a graph is one of the most versatile parameters in controlling fixation acceleration. We further show how to apply our algorithms and analytic results to two evolutionary optimization problems and show that the rate at which a population discovers the optimum can be tuned near continuously by adjusting the higher-order topology of the agent population. We highlight that the effects of population structure on the rate of solution discovery are more subtle than previously recognized, and find that decelerators with longer times to fixation are able to reach the optimal solution faster than accelerators when the space of possible solutions is rugged and complex. No one population structure is perfect for all optimization tasks, and careful consideration of the problem landscape is necessary for selecting the correct approach. Our work reinforces the importance of understanding and exploring common design principles of biological networks for engineering highly evolvable artificial systems. Model We compute probabilities and times to fixation of a new mutation a that appears in an initially homogeneous population of A-type individuals and use a Moran birth-death process to track variant frequency changes. An individual with the A allele is assumed to have fitness one, while an individual with allele a has assigned fitness (1 + s). The population size is kept fixed at N. We use unweighted, undirected graphs to represent the population structure. A node in the graph represents an individual in the population. The network edges are proxies for the local pattern of replacement and substitution: an edge can thus represent spatial proximity or the social architecture of interaction between individuals in the population. A node can also be interpreted as a homogeneous subgroup of individuals and the edges as migration corridors between them, with the assumption that the timescale of a mutation traveling between nodes is much larger than the time scale of fixation within a node [30–32]. We study both the Birth-death and the death-Birth process on these networks. These two update rules are equivalent in well-mixed populations, however they have been shown to lead to drastically different long-term evolutionary dynamics on graph structures [2, 5, 33]. In the Birth-death process, at every time step, a random individual is chosen for reproduction, proportional to fitness. An individual occupying a neighboring node is then randomly chosen for death and replaced by the new offspring. In the death-Birth update, an individual is first randomly chosen for death and is replaced by the offspring of a neighboring node. The neighbors of the vacant node compete for this vacancy and one neighbor is chosen for reproduction with probability proportional to fitness (Fig 1B). Higher-order motifs create graph geometries produced by different patterns of interconnecting nodes. To describe and constrain random graphs by levels of organization in successively finer detail, [34] previously introduced the mathematical formalism of dK-distributions, where d denotes the dimension of an interaction involving d nodes. Thus, the 0K-distribution represents the average node degree. The 1K-distribution is the degree distribution of the graph, which specifies the fraction of nodes having degree i in a graph. 2K refers to the pattern of interactions, or mixing pattern, between two nodes. It is defined as the fraction of edges that connect a node of degree i to a node of degree j. Network assortativity, r, is also often used to describe mixing behavior in a network and is defined as the Pearson correlation coefficient of the degrees at either ends of an edge [35]. The 3K distribution describes interactions between groups of three nodes. For the 3K distribution there exist two possible connection topologies: wedges, chains of three nodes connected by two edges and their dimension three complement topology, triangles, cliques of three nodes (Fig 1A). Here we consider interactions of d ≥ 3 and particularly focus on the role of 3K interactions in shaping probabilities and times to fixation. To systematically study the role of higher order motifs, the challenge lies in tuning these higher levels of organization while keeping lower levels constant. This is complicated by the fact that, for example, varying the number of triangles in the network also changes the degree distribution and the network mixing pattern, which we have previously shown to significantly affect probabilities of fixation [5]. We start by studying the role of higher order interactions in random regular graphs, where all the nodes have the same degree [36–38]. This simplifies the problem, since all k-regular graphs of fixed degree k share the same mixing pattern. We use degree-preserving edge swap operations to computationally tune the number of higher-order network motifs, while maintaining the lower level of organization, the degree distribution, constant [39]. We build algorithms that combine these approaches with simulated annealing to produce a sampling network generation algorithm. For example, for tuning the transitivity or fraction of closed triangles ϕ over all triples in the network to a target fraction ϕ target , we accept an edge swap based on the Metropolis–Hastings criterion [40] (1) Here, ϕ′ is the fraction of triangles after the proposed edge-swap, and P is the probability density function of a Gaussian with mean ϕ target and variance σ2. The algorithm always accepts edge swaps that bring the network closer to the target and accepts other edge swaps randomly. We increase σ by a factor of 1.001 every Nk times, where N is the size of the graph. We run the algorithm for k−million steps. Since edge swapping can cause the graph to become disconnected, we use the heuristic outlined in [41] to ensure graph connectivity. For regular graphs, the above algorithm will sample uniformly all connected regular graphs with the target fraction of triangles. This is because without the criterion in Eq (1), the Markov chain converges to a unique stationary distribution, which is the uniform distribution over the state space of all connected graphs with the same degree distribution [41]. The Markov chain with the condition in Eq (1) has a stationary distribution of where |G(k, ϕ′)| is the cardinality of connected graphs G with degree k and fraction of triangles ϕ′. In Section 3 in the S1 File, we also present an analysis of the d = 4 topology. We expand this approach to degree-heterogeneous graphs by designing graphs with two distinct degrees (Fig 1C). This allows us to seamlessly tune the number of triangles without changing degree distributions and mixing patterns of the graphs. These graphs are designed starting with two random regular networks of size n 1 and n 2 and fixed degrees k 1 and k 2 , and we use edge swaps to connect them. Here, the edge swap operation we previously used for the analysis of k-regular graphs is unsuitable, since that algorithm also alters the graph mixing pattern. This is not an issue for k-regular graphs, but not controlling for the mixing pattern of a graph with variance in degree can potentially lead to erroneous interpretations (see Fig A in the S1 File). For studying degree-heterogeneous networks, we instead use the dK-preserving rewiring, which is a generalization of the degree-preserving edge swap for higher-order motifs [34]. Fig 1D illustrates a 2K-preserving rewiring. Here, we replace the condition in Eq (1) with the following and alway accept the edge swap if the change in the fraction of triangles is towards the target, as well as according to (2) if the change is in the other direction. Here, 1/γ is the annealing temperature which controls how stringent the criterion must be and is decreased as the graph generation algorithm proceeds. For degree-heterogeneous graphs, we are not guaranteed uniform sampling over the fraction of triangles. Once the network structure is set, we use ensembles of at least 10, 000 Monte Carlo simulations, as well as analytic approaches, as described in the next section, to compute the probabilities and times to fixation of the new mutation in the population. Discussion To design spatial structures with accelerated (or decelerated) rates of evolution, we must first build a systematic understanding of how the spatial arrangement of a population shapes probabilities and times to fixation of new mutants entering a population. Of general interest are graph structures that can greatly amplify selection, but pay minimal cost in the time it takes for the mutant to sweep through and fix in the population. This type of structure is theorized to be desirable for two reasons: 1) fixation probability is important when the rate-limiting step is waiting for an advantageous mutant to occur; 2) fixation time is important in high mutation regimes, where the abundance of beneficial mutations negates loss due to stochastic drift. There are various approaches to finding such structures. For example, [22] developed a genetic algorithm that explores a space of connected graphs to optimize for fixation probability and fixation time. Separately, “α-bipartite” graphs [23] and “selection reactors” [52] use weighted edges to push the boundary of tradeoff between fixation probability and time. Here, we take a different approach by focusing on undirected and unweighted graphs and build a systematic understanding of how the different levels of structural organization of a network shape evolutionary dynamics. Our results show that 3-dimensional structural motifs are the most effective network property for tuning temporal dynamics and times to fixation of new mutants in the population. While our most novel finding involves triangles and higher-order network motifs, our analysis also reveals connections between the evolutionary role of lower-levels of network organization (the node and edges) and previous population structure models. This is because the heterogeneous graph families we explore can, in some cases, be topologically mapped to classical deme-based population structures. For example, groups of nodes with distinct degrees can be thought of as population subdivisions, and the edges connecting these groups can be viewed as migration paths between demes. When the node degrees of the groups are identical (as in the case of k-regular graphs), we find that higher order networks motifs do not change probabilities of fixation, only times to fixation and this behavior qualitatively mirrors what has been observed in previous classical models [53, 54]. More generally, we find that the number of closed triad interactions influences fixation time with minimal impact on fixation probabilities when the graph is not assortatively mixed. Triad interactions alter evolutionary outcome similar to the migration rate in a deme model, although through a completely different mechanism. The key distinction between the network-based models and previous deme-based models lies in the fact that deme models do not consider the structural properties within the subdivision, beyond the deme size. In contrast, the Moran process on graphs inherently models properties such as the number of connections, triad interactions, and higher-order connections within node groups. When there are differences in these structural properties between node groups, in addition to the quantitative differences observed in fixation time, we also observe qualitative deviations in fixation probabilities. Structural properties like degree distribution and mixing patterns tend to simultaneously affect fixation probability and fixation time, as these two quantities are closely intertwined [5, 23]. Here we find that the number of triangle motifs can affect both probabilities and times to fixation for graphs with non-zero assortativity. While the field of evolutionary graph theory is usually interested in accelerators that minimize time to fixation, we also use tuned network structures for two very different optimization problems and show that accelerators don’t necessarily lead to faster rates of adaptation. Optimal structures depend heavily on the fitness landscape considered, as predicted by Sewell Wright’s shifting balance theory [55]. One case where faster time to fixation can hurt the long term adaptation of a population is when the fitness landscape consists of multiple fitness peaks separated by valleys. Here, an accelerator can prevent individuals from adequately exploring these valleys to eventually discover higher fitness peaks. We find that graphs with reduced fraction of triangles do well when evolving populations to solve optimization problems on smooth landscape, however, the same graphs can actually lead to slower solution discovery when the optimization landscape is produced from a rugged Rastrigin function. The observed fitness trajectory shows similar tortoise-and-hare pattern seen in structured bacterial population [56], indicating that our observation is the consequence of the shifting balance process. Our results highlight the need to search for structures outside of fast amplifiers, as the diversity of evolution and optimization landscapes dictates that certain structural properties are preferred in specific contexts, but no single set of properties is ideal for all scenarios. Moreover, we’ve demonstrated that degree distribution, mixing patterns, and triangle counts are effective network properties for systematically exploring the space of potential evolutionary outcomes within structured populations. Our theoretical treatment is limited to studying probabilities and times to fixation, but more factors contribute to both natural and artificial evolving systems. For example, evolutionary algorithms often use high mutation rates and multiple mutants are available simultaneously in a population. Different clones compete and can interfere with the fixation of one another [57]. Effective traversal of the fitness landscape also relies on the existence and maintenance of mutation islands. Fixation time might not be the most direct quantity that describes the clonal dynamics in such a population and tools that study clonal interference could potentially also provide additional insights into the problem. Our work nonetheless provides the first exploration of the role of higher-order network motifs in shaping evolutionary dynamics and highlights the importance of using evolutionary design principles to engineer highly evolvable, adaptable artificial systems. Supporting information S1 File. The Supplementary material contains the complete analytical derivations and the supplementary figures associated with the manuscript. https://doi.org/10.1371/journal.pcbi.1011905.s001 (PDF) Acknowledgments This research was done using resources provided by the Open Science Grid, which is supported by the National Science Foundation award 1148698, and the U.S. Department of Energy’s Office of Science. [END] --- [1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011905 Published and (C) by PLOS One Content appears here under this condition or license: Creative Commons - Attribution BY 4.0. via Magical.Fish Gopher News Feeds: gopher://magical.fish/1/feeds/news/plosone/