
==== Front
Genetics
Genetics
genetics
Genetics
0016-6731
1943-2631
Oxford University Press US

38639307
10.1093/genetics/iyae055
iyae055
Investigation
Population and Evolutionary Genetics
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
Evolutionary graph theory beyond single mutation dynamics: on how network-structured populations cross fitness landscapes
Kuo Yang Ping Computational Biology Department, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA 15232, USA

https://orcid.org/0000-0002-4296-2326
Carja Oana Computational Biology Department, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA 15232, USA

Ralph P Editor
Corresponding author: Computational Biology Department, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA 15232, USA. Email: oana.carja@gmail.com, ocarja@andrew.cmu.edu
Conflicts of interest. The author(s) declare no conflicts of interest.

6 2024
18 4 2024
18 4 2024
227 2 iyae05507 2 2024
01 4 2024
06 5 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of The Genetics Society of America.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs licence (https://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial reproduction and distribution of the work, in any medium, provided the original work is not altered or transformed in any way, and that the work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

Spatially resolved datasets are revolutionizing knowledge in molecular biology, yet are under-utilized for questions in evolutionary biology. To gain insight from these large-scale datasets of spatial organization, we need mathematical representations and modeling techniques that can both capture their complexity, but also allow for mathematical tractability. Evolutionary graph theory utilizes the mathematical representation of networks as a proxy for heterogeneous population structure and has started to reshape our understanding of how spatial structure can direct evolutionary dynamics. However, previous results are derived for the case of a single new mutation appearing in the population and the role of network structure in shaping fitness landscape crossing is still poorly understood. Here we study how network-structured populations cross fitness landscapes and show that even a simple extension to a two-mutational landscape can exhibit complex evolutionary dynamics that cannot be predicted using previous single-mutation results. We show how our results can be intuitively understood through the lens of how the two main evolutionary properties of a network, the amplification and acceleration factors, change the expected fate of the intermediate mutant in the population and further discuss how to link these models to spatially resolved datasets of cellular organization.

evolutionary graph theory
spatial structure
rugged fitness landscapes
stochastic tunneling
cell graphs
National Institute of General Medical Sciences 10.13039/100000057 R35GM147445 United States-Israel Binational Science Foundation 10.13039/501100001742 2019266 NIH 10.13039/100000002 T32 EB009403
==== Body
pmcIntroduction

In recent years, evolutionary graph theory has started to reshape our understanding of how spatial structure can direct evolutionary dynamics (Lieberman et al. 2005; Ohtsuki and Nowak 2006; Ohtsuki et al. 2006; Santos et al. 2006; Ohtsuki et al. 2007; Poncela et al. 2007; Allen et al. 2013; Jiang et al. 2014; Maciejewski et al. 2014; Leventhal et al. 2015; Kuo et al. 2021; Su et al. 2022; Kuo and Carja 2024). Using tools from network theory, these models allow us to tune the heterogeneity of spatial structure, beyond what is possible with deme-based or lattice-based models (Wright 1943; Kimura and Weiss 1964; Maruyama 1970a; Slatkin 1981; Whitlock and Barton 1997; Whitlock 2003; Carja et al. 2014). The nodes of the graph represent individuals in the population and edges are proxies for the local pattern of replacement and substitution. Nodes can also be interpreted as genetically homogeneous subpopulations where genetic drift is ignored and a new beneficial mutation arriving in the node subpopulation is assumed to fix immediately (but see also studies where this assumption is relaxed Marrec et al. 2021; Yagoobi and Traulsen 2021; Yagoobi et al. 2023).

By incorporating heterogeneity in spatial structure, these network models have been shown to greatly extend the range of possible evolutionary outcome of a population and can boost the selective benefit of new mutations and, reversely, suppress the spread of deleterious mutants (Adlam et al. 2015; Hindersin and Traulsen 2015; Tkadlec et al. 2020; Allen et al. 2021; Kuo et al. 2021). Network structures have also been shown to affect times to fixation of new mutants in the population (Hindersin and Traulsen 2014; Tkadlec et al. 2019; Kuo and Carja 2024). One of the main obstacles to theoretical progress in this area has been the analytical difficulty of multidimensional stochastic models and previous results have been restricted to studying times and probabilities of fixation for one single mutation appearing in the population. These previous single-mutation frameworks can therefore only be used to predict sequential fixation dynamics, in the limit of small mutation rates.

However, complex traits often arise from interactions between multiple gene products and, in order to reach a fitness peak, populations often need to cross fitness valleys or plateaus (Wright et al. 1932; Burch and Chao 1999; Wood et al. 2007; Kvitek and Sherlock 2011; Huang 2013; Vogelstein et al. 2013; Rogers et al. 2018; Acar et al. 2020; Salehi et al. 2021). Well-mixed large populations have been shown to cross wide fitness valleys remarkably quickly, suggesting valley-crossing dynamics is common, even when mutations that directly increase fitness are available (Weinreich and Chao 2005; Jain and Krug 2007; Weissman et al. 2009). Furthermore, under elevated mutation rates, crossing fitness valleys occurs not through sequential fixation, but mostly through stochastic tunneling, especially when the population size is large and when the intermediate mutant is deleterious (Nowak et al. 2002; Komarova et al. 2003; Iwasa et al. 2004; Weissman et al. 2009). This makes previous evolutionary graph theory results in the limit of sequential fixation no longer predictive for populations where stochastic tunneling and fitness valley crossing is prevalent.

In addition to the work done assuming a well-mixed population, previous spatial models using deme-based or lattice-based structures have shown that these spatial structures can accelerate the crossing of fitness valleys and plateaus, while delaying the evolution of complex traits with advantageous intermediates (Bitbol and Schwab 2014; Komarova 2014; Durrett and Moseley 2015). However, unless local differences in the nodes or demes are assumed, these regular, lattice-based structures do not change probabilities of fixation for intermediate mutants, and only affect valley crossing through changes to the extinction or fixation time of the intermediate (Pollak 1966; Maruyama 1970b; Lande 1979). In contrast, more heterogenous graph spatial topologies can shape both probabilities and times to fixation for intermediate mutants, and we have yet to understand how the complex interplay between structural properties that amplify or suppress selection on the intermediate mutant, combined with the effects of motifs that accelerate or decelerate mutational spread shape evolutionary outcomes.

Here we study the role of network topologies in shaping multi-mutational dynamics and, in particular, probabilities of fitness valley crossing and stochastic tunneling. In agreement with previous lattice and deme-based models (Bitbol and Schwab 2014; Komarova 2014), our results show that heterogenous spatial structure can promote fitness landscape crossing by allowing intermediate mutants to persist for longer, until the final beneficial mutation can appear in the population. However, in contrast with previous models, we show that more complex population structures rarely purely promote or suppress fitness landscape crossing and instead, switch across regimes as a function of the fitness of the intermediate mutants. We study the role of network properties in shaping rates of fitness valley crossing and compare these rates across well-studied families of graphs.

We also discuss how to apply these network-based approaches to large-scale datasets of spatial organization. We use previously studied datasets of the cellular networks of the stem cell niches in the bone marrow and apply our evolutionary model to study rates of mutation accumulation and leukemia initiation. Our results show that these cellular spatial architectures reduce the probability of neoplasm initiation across biologically relevant mutation rates and fitness distributions, compared to well-mixed populations. However, we show that there exists a threshold mutation rate (under exposure to carcinogens, for example) above which the bone marrow structure shifts from a suppressor into a promoter of neoplasm initiation. Our results provide the groundwork for further exploration of the role of heterogeneous spatial topology in shaping rates of mutation accumulation in both engineering and biological settings.

Model

We consider an asexual population of N haploid individuals and study the process by which this population acquires a beneficial trait that requires mutations at multiple loci. This process is usually referred to as the “K-hit” mutational process. Here we study the K=2 case (illustrated in Fig. 1). We assume that, at time t=0, a single first mutant, with reproductive fitness (1+s), appears in a population consisting of only wild-type individuals of fitness equal to one. Depending on the sign of the selection coefficient s, this first mutation can be beneficial (i.e. s>0), neutral (i.e. s=0), or deleterious (i.e. s<0). With the probability of mutation upon reproduction μ, individuals carrying the first mutation can acquire a second mutation assumed to be extremely beneficial, such that fixation of double-mutant individuals is guaranteed with probability one. Therefore, the probability that all the individuals acquire the beneficial trait is equal to the probability of acquiring the second mutation. The second mutant can spread through the population either through a step-wise fixation process or through stochastic tunneling, without ever visiting a population fixed on the intermediate mutant (Fig. 1c).

Fig. 1. Illustration of the 2-hit process on network-structured populations. Panel a) A two-step mutational process. The first mutation can be beneficial, neutral, or deleterious. The second mutation is assumed to be extremely beneficial, such that fixation is guaranteed. Panel b) The Bd (Birth–death) update rules. Panel c) There are two ways for the second mutation to fix in the population. The first is through step-wise fixation, where the intermediate mutant reaches fixation before the second mutation occurs (I). The second is through stochastic tunneling, where the second mutation is acquired before fixation of the intermediate (II).

To represent heterogeneous population spatial structure, we use unweighted and undirected graphs, where each node in the graph represents one individual in the population and edges are proxies for the local pattern of replacement and substitution. A node can also represent 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. We study the role of population structure in shaping rates of fitness valley crossing by analyzing how graph properties determine the probability that the population reaches the fitness peak.

To study graph properties independently of one another, we use well-known graph generation algorithms implemented in NetworkX (Hagberg et al. 2008), as well as algorithms that allow us to systematically tune network parameters (Kuo et al. 2021; Kuo and Carja 2024). We analyze a wide variety of graph families including preferential attachment graphs (Barabási and Albert 1999), bipartite graphs (Asratian et al. 1998), detour graphs (Möller et al. 2019), star-like graph (Tkadlec et al. 2019), small-world networks (Watts and Strogatz 1998) used to model properties of social networks, and random geometric networks (Waxman 1988; Penrose et al. 2003), mathematical representations of populations embedded in Euclidean space. We also introduce regular island graphs which are created by adding random edges to two separate k-regular graphs and shortcut graphs which are created by adding edges connecting the loop of a detour graph. We present a detailed list of the graphs used in this study in the Materials and methods section.

We write analytic predictions for the probability of fitness valley crossing as a function of the population size N, the mutation rate μ, the intermediate mutant selection coefficient s, and the statistical properties of the spatial structure and compare these analytic predictions with results from Monte Carlo simulations. These analytic results allow us to generalize our insights beyond the graph families discussed here. To simulate the evolutionary trajectories of the population, at time t=0, we introduce one intermediate mutant on a random node of the network and use a Moran Birth–death (Bd) process to track mutant frequency changes (Fig. 1b). The birth–death process assumes that at each time step, an individual from the population is selected to reproduce proportional to fitness and a random neighbor node is selected to die, leaving an unoccupied node for the offspring of the reproducing node. This is the essential difference that allows us to study the role of local population structure, in comparison to a well-mixed population: in a well-mixed population, a node from the entire network would be randomly selected for death. Each time an intermediate mutant is selected to reproduce, with probability μ, it can acquire the second beneficial mutation. The simulation ends when either the intermediate mutant dies out or the second mutation is acquired. The fitness landscape crossing probability is then determined by the fraction of simulations that end with acquisition of the second mutation. Landscape crossing time is determined by averaging over the time it takes for the beneficial trait to fix, given fixation. Landscape crossing probabilities and times are computed using at least 107 simulations, for any given network structure.

Results

Landscape crossing is fundamentally shaped by both the probability and the time it takes for the first mutant to spread in the population and create opportunities for the appearance of the second mutant. Previous work has shown that, for single mutations, the evolutionary role of complex spatial structure can be quantified by analyzing two essential network properties: the network amplification factor, which shapes probabilities of fixation compared to well-mixed populations (Lieberman et al. 2005; Kuo et al. 2021) and the network acceleration factor, which shapes the time to fixation for new mutations in the population (Kuo and Carja 2024). The network amplification factor α quantifies how to rescale the selection coefficient s for the well-mixed model to obtain the same probability of fixation as an allele with selection coefficient s on a network population. In other words, a mutation with a selective benefit s in a graph population, would have a probability of fixation corresponding to an equivalent selection coefficient αs in a well-mixed population. If α>1, the graph amplifies selection, if α<1 the graph suppresses selection and if α=1 the graph does not change probabilities of fixation compared to the well mixed. Similarly, we previously introduced the network acceleration factor, λ (Kuo and Carja 2024), and showed that this network property shapes times to fixation, making the evolutionary dynamics λ times faster than in a well-mixed population, for λ>1.

While the amplification and acceleration have closed-form approximations using network descriptors, such as the mean and variance in degree (Kuo et al. 2021; Kuo and Carja 2024), they can also be computed using simulations and the definition from Lieberman et al. (2005) and solving

(1) Φ=1−(1+s)−α1−(1+s)−αN

for α, where Φ is the single mutation fixation probability. The amplification factor is constant only for weak selection (McAvoy and Allen 2021), and varies with selection when selection is strong for many graphs (Voorhees 2013; Voorhees and Murray 2013; Allen et al. 2020). In this paper, we assume amplification factor is in the constant regime and we estimate it using Ns=1. See more detailed description of how to compute network amplification and acceleration in Materials and methods.

The acceleration factor is defined as the ratio of conditional mean time to fixation for the equivalent well-mixed model (with the mutant having selection coefficient αs) and the network model with selective coefficient s (Kuo and Carja 2024), λ=Twm(αs)Tgraph(s).

We will show that the interplay between how the network structure changes the fixation probability, the fixation time of the intermediate mutant and the mutation rate μ leads to rich landscape crossing dynamics that cannot be predicted or interpreted through the simpler lens of single mutation results. We start by presenting heuristic descriptions of the roles of these two evolutionary network parameters, α and λ, in shaping fitness landscape crossing and, using this intuition, we then present our full analytic results, with a detailed description of the analytic approach in the Supplementary Material.

The role of the network acceleration factor

Let us assume a network that only affects time to fixation and not the probability of fixation (α=1), for a single mutant. Any k-regular graph, for example, satisfies this property. There are two independent evolutionary scenarios that can lead to the acquisition of the second mutant.

The first scenario occurs when the intermediate mutant lineage is poised for extinction and the second mutation appears in the population before the lineage goes extinct. In a well-mixed population, the lineage reaches a size T, smaller than the establishment threshold of 1/s, in T generations and then drifts to extinction in another T generations, with probability 1/T. The overall shape of such a trajectory is shown in Fig. 2a and the total number of mutational opportunities to acquire the second mutation depends on the area under the trajectory, W(T)=T2 (Weissman et al. 2009). The probability that a second mutation arises on the background of this mutant lineage is the expectation of the probability that a second mutation occurs given a trajectory that reaches maximum size of T, Φ(W(T)), over all possible such trajectories,

(2) P=∫Φ(W(T))p(T)dT.

When the magnitude of s is strong enough (|s|>μ), the rate at which double-mutants are produced is dominated by rare lucky intermediate mutant lineages that survived for T∼1/s (Weissman et al. 2009). Therefore, the probability of landscape crossing under this scenario is

(3) P∼μW(T=1s)p(T=1s)=μs.

A network-structured population with acceleration factor λ effectively stretches the time the lineage reaches a size of T to T/λ generations, and also the time to extinction to 2T/λ generations. The area under the trajectory thus becomes T2/λ and, using equation (3), the rate of landscape crossing is changed by a factor of λ−1, to μ/λs.

Fig. 2. The role of the network amplification and acceleration factors in shaping the frequency of the intermediate mutant and fitness valley crossing probability. Top panels show dynamics for networks of constant amplification equal to 1. Bottom panels show dynamics for networks with acceleration equal to 1. Panel a) Typical trajectory of the intermediate mutant frequency, given that the intermediate mutant goes extinct. Population structures that increase extinction time lead to an increased probability of acquiring the second mutation. Panel b) Typical trajectory of the intermediate mutant frequency, given that the intermediate mutant fixes. Panel c) Ratio of crossing probability between k-regular graphs and the well-mixed model. Here, μ=10−4 and the lines represent our numerical approximation to equation (4) in the limit of infinite time, using Supplementary equation (54), while the dots are results from 107 simulations per parameter and network. Network acceleration factors calculated using equation (10). Panel d) Typical trajectory of the intermediate mutant frequency, given the intermediate mutant goes extinct. A suppressor increases the line of transition between the stochastic and the deterministic dynamics, allowing the mutant to drift for longer and to a higher frequency. The structure would, therefore, increase the probability of acquiring the second mutation. Panel e) Typical trajectory of the intermediate mutant frequency, given the intermediate fixes. Panel f) Ratio of crossing probability between shortcut graphs with acceleration factor equal to one and the well-mixed model. Here, μ=10−4 and the lines represent our numerical approximations to equation (4) (Supplementary equation (54)), while the dots are results from 107 simulations per parameter and network. Network amplification factors empirically determined using equation (9).

The second scenario occurs when the intermediate mutant lineage would reach fixation on its own. In this scenario, the trajectory behaves stochastically below the establishment threshold, 1/s. As soon as the mutant frequency crosses this threshold, the population dynamics becomes essentially deterministic (Fig. 2b). At any point along the trajectory, the second mutation can occur and carry the mutant lineage to fixation. Since this scenario is conditioned on the first mutant fixing, the population is thus guaranteed to eventually acquire the second mutation and therefore, in this regime, the network structure does not change the probability of crossing the fitness valley. The total probability of crossing the fitness landscape is the sum of the probabilities of acquiring the second mutation under the two independent evolutionary processes (Fig. 2c). When the first mutant is strongly deleterious, the population depends on the second mutation appearing in time to cross the fitness landscape and the acceleration factor of the network changes the rate of fitness valley crossing by a factor of λ−1. Otherwise, the probability of stochastic tunneling remains unchanged by the network structure.

The role of the network amplification factor

Let us now assume a network structure that changes the fixation probability of a new mutant, but not the time to fixation (λ=1) (Fig. 2). An example of this type of structure are the shortcut graphs we present in Fig. 2f.

For a network with amplification factor α, the mutant experiences an effective selection strength of αs in a well-mixed population. Therefore, if the intermediate mutant lineage is on its way to extinction, with probability αT, the mutant reaches an intermediate mutant size of Tα≪1αs in Tα generations and goes to extinction in another Tα generations. The network amplification factor does not change the expected trajectory that reaches size Tα, but instead changes the limit on the size that an intermediate mutant trajectory can reach, from a maximum of 1s to 1αs. Using equation (2), the rate of landscape crossing under this scenario becomes μ/αs.

In contrast, if the intermediate mutant lineage is destined to reach fixation, the crossing probability is the same as the probability of fixation of the first mutant with effective fitness αs, since the population is guaranteed to acquire the second mutation. In this regime, the network structure does not change the probability of crossing the fitness valley, beyond modifying the effective selective coefficient of the intermediate mutant. Considering the two scenarios together, when the first mutant is strongly deleterious, the population depends on the first scenario to cross the fitness landscape, and the rate of landscape crossing is changed by a factor of 1/α. In the limit where the first mutant is strongly beneficial, fixation of the intermediate becomes the dominant process, and the crossing probability is approximately αs (Fig. 2f).

The general analytic approximation

The previous sections provide heuristic intuition into the independent roles of the network amplification and acceleration factors in shaping probabilities of valley crossing. We detail the general analytic derivation in the Supplementary Material and show how to derive approximations for the two quantities of interest: the probability that a second mutation arises in the population, and the expected time for that to occur. To do so, we use the diffusion approximation. Direct application of the diffusion approximation is difficult due to the added complexity of the network pattern of replication and replacement. Analytic approaches that make use of the adjacency matrix of the network (which uniquely identifies the graph) and its associated transition probabilities become intractable for large networks, since they track a Moran process with 2N states. The approach we take here and the core idea behind our approximation is that we can greatly reduce the complexity of the problem by using the node degree distribution, and only keep track of the mutant frequencies for all groups of nodes of the same degree. While the degree distribution might not uniquely represent the network and some of the graph information is lost, this approach nonetheless greatly reduces the number of possible states in the Moran model. For that reason, the analysis is split into two parts.

In Supplementary Section 1.1, we first show that if we assume nodes of the same degree in the network to be ’evolutionarily’ identical, we can derive the expected change in the intermediate mutant frequency for each group of nodes of equal degree. To this end, we first observe that, since wildtype to mutant replacements depend on the number of edges between mutants and wildtype individuals, we can write the expected change in mutant frequencies at nodes of the same degree as a function of the frequencies of the different edge types possible in the graph. By first deriving the expected change in edge type frequencies, we show that, if we assume weak selection and weak mutation, the intermediate mutant frequency for each group of nodes with distinct degrees is the same. We define this quasi-equilibrium node frequency as qa, which depends on the initial mutant placements on the network. The quasi-equilibrium edge frequencies also depend on qa. Selection and mutation, however, will slowly alter this quasi-equilibrium.

In the second part of the analysis, Supplementary Section 1.2, we write the expected change in qa due to selection and mutation. We write the variance in the change of qa and, in the limit of weak selection s≪1, weak mutation Nμ≪1, and large population size 1/N≪1, write the Kolmogorov backward equation

(4) ∂∂tΦ(qa,t)=λqa(1−qa)(1N∂2∂qa2+αs∂∂qa)Φ(qa,t)+Nμqa[1−Φ(qa,t)],

where Φ(qa,t) is the probability of landscape crossing after time t, when starting with quasi-equilibrium node frequency qa (see also Supplementary equation (39)). The probability of landscape crossing when starting with one intermediate mutant in a random node on the graph is then defined as Φ=limt→∞Φ(qa=1N,t).

We next show that equation (4) is similar to the one governing the Moran process for a well-mixed population, if the appropriate mapping changes are done that account for the roles of the two network evolutionary parameters, α and λ: 1. an effective mutant selection coefficient αs and 2. an effective extended fixation time Tfix/λ (see Supplementary Section 1.2.2 and, specifically, equation (48)).

Thus, even though seemingly the network structure can introduce a lot of complexity and heterogeneity, we show that it is sufficient to know the amplification and acceleration parameters of the population structure in order to use the arsenal of tools developed for well-mixed evolutionary theory to understand landscape crossing in a network population. The crossing probability and time can be found by simply utilizing the transition matrix for the rewritten, analogous well-mixed system of equations. Lastly, we show how to use finite difference numerical approximations to solve for the crossing probability (Supplementary equation (54)) and crossing time (Supplementary equation (64)).

We also show that the dynamics of fitness landscape crossing reduce to three simple cases when the population size becomes large. As N becomes large, we can approximate the tunneling probability using continuous time branching processes (Weissman et al. 2009) (see Supplementary Section 1.2.3) and write

(5) Φ=αs−λ−1μ+(αs+λ−1μ)2+4λ−1μ2(1+αs).

Under weak mutation, series expansion in μ (see Supplementary equations (52)–(56)) reduces equation (5) to

(6) Φ={μαλ|s|fors≪−2μλ−1μfor−2μ≪s≪2μαsfors≫2μ.

In the case of a deleterious intermediate mutant, s≪−2μ, equation (6) recaptures the crossing probability of μ/λ|s|, when α=1, and μ/α|s| when λ=1. A network structure with amplification parameter α and acceleration parameter λ changes the crossing probability of the well-mixed population, μ/|s|, by 1/αλ. Therefore, an amplifier, as defined in the single mutant case, can promote the fixation of a mutant lineage with deleterious intermediates, as long as the product of amplification and acceleration parameters is less than one. This is in stark contrast to the case of sequential fixation, where, in the case of a deleterious intermediate mutant, an amplifier always reduces the acquisition probability of the second mutation (Fig. 3a). With neutral intermediate mutants, −2μ≪s≪2μ, selection on these mutants has a negligible effect on the crossing probability: equation (6) only depends on λand not α. In the case of a beneficial intermediate mutation, 2μ≫s, the population crosses the fitness landscape through intermediate mutants that are destined to fix. Since λ only affects the probability of acquiring the second mutation before fixation of the first, the crossing probability only depends on α and the crossing dynamics reduce to the sequential mutation case.

Fig. 3. Network structures can switch from promoting to suppressing fitness valley crossing as a function of their amplification and acceleration properties. Panel a) We compare fixation probabilities of the second mutation in a structured population (a bipartite graph, with α=1.7, λ=0.1) and a well-mixed population, under stepwise fixation dynamics (μ=10−12). Here, dots represent simulation results from 107 replicates, and lines are calculated using Supplementary equation (54). The population size is N=100. Panel b) Fixation probabilities of the second mutation in a bipartite graph, α=1.7, λ=0.1, compared to those for a well-mixed population under a stochastic tunneling scenario (μ=10−7). Here, dots represent simulation results from 107 replicate simulations per parameter and network, and lines represent analytic approximations from Supplementary equation (54). In this regime, the population structure crosses the well-mixed line twice, for negative selection strengths, and switches from promoting to suppressing fitness landscape crossing. Zoomed in version around Ns=0 in the insert. The population size is N=100. Panel c) The general amplifier suppressor classification for the single-mutation case. Panel d) The seven different scenarios for the ratio of crossing probabilities between a network over the well-mixed population, as a function of selection strength Ns and the two main evolutionary properties of the networks, the amplification and the acceleration factor. Here, the lines represent analytic approximations using equation (6) and the red and blue colors denote suppression and amplification of selection, respectively. The population size is N=1,000.

These results highlight that the intersection between the crossing probability of the structured population and that of the well-mixed case can shift from neutrality. For example, for the amplifier of selection in Fig. 3b, there exists a region between the new intersection and neutrality, where the population structure behaves as a suppressor (see insert). In addition, for very deleterious intermediates, there exists an additional intersection, where the population structure acts as a suppressor on the left and an amplifier on the right. Single mutation fixation probabilities can therefore lead to erroneous inferences of multi-mutational landscape crossing dynamics.

The network structure can still act as an overall suppressor or amplifier, but only in a small region of the parameter space where the acceleration factor is exactly one. Beyond this knife-edge case, the general dynamics become more complicated, with the existence of piece-wise amplification or suppression, and we observe seven additional categories of graphs as the acceleration factor shifts away from one (Fig. 3d).

These additional categories of graphs can be understood using the intuition developed in Fig. 2 and equation (6). We start from the largest values on the x-axis and go clockwise (light blue regime). In the regime where λ<1, α>1 and αλ>1, the population structure is an amplifier except for weakly deleterious intermediate mutation. As λ decreases and αλ starts decreasing below one, a population structure can act as an amplifier for moderately deleterious mutations, but behaves like a suppressor for weakly and strongly deleterious intermediate mutants. As λ decreases further, αλ≪1 and the population structure will increase crossing probability compared to a well-mixed population, regardless of the fitness of the intermediate (darker blue regime).

A network structure with α<1 and αλ≪1 acts as a suppressor, increasing the crossing rate for a deleterious intermediate mutant. For a neutral mutation, the crossing probability depends on 1/λ and, for λ<1, the population structure acts as an amplifier for weakly beneficial mutations. As the strength of selection keeps increasing however, the crossing probability depends on 1/α and the population structure starts acting as a suppressor, decreasing the rate of valley crossing (red regime). As λ increases above 1, but αλ<1, the population structure becomes an amplifier for weakly deleterious mutations (orange regime). There is a further sub-category of behavior where the ratio of crossing probabilities is higher than one for a moderately deleterious intermediate mutant. As λ increases even further, αλ is no longer smaller than one (light orange regime). The population structure starts behaving like an amplifier for deleterious mutations, crosses neutrality with a probability below the well mixed, and behaves as a suppressor for beneficial mutations. No graphs exist in the upper right corner of the parameter space, since networks that are amplifiers have longer fixation time than well-mixed populations (Tkadlec et al. 2019).

Fitness valley crossing for different network families

Our results provide a unifying way of predicting multi-mutational dynamics by distilling properties of complex network organization through two single parameters: the acceleration and amplification factors. We can use this framework to compute the rate of fitness landscape crossing across known network families and we can design networks with the right ratio of amplification and acceleration to optimize for the rate of valley crossing. For example, for the case in equation (6) where the dynamics is controlled by 1/αλ, the graph with the highest crossing probability is the star graph (Fig. 4a). Since the intermediate mutant lineage must wait for ∼1/αλ for the second mutation to occur, the crossing time also scales with αλ (Fig. 4b).

Fig. 4. The rate of fitness valley crossing across network families. Panel a) Comparison of our analytic approximation (using Supplementary equation (54)) and simulations (using 107 replicate Monte Carlo simulations per parameter and network) for known graph families, colors as in Panel c) Panel b) Comparison of our analytic approximation for crossing time (Supplementary equation (64)) and simulations (107 replicate Monte Carlo simulations per parameter and network) for known graph families. Panel c) The rate of fitness valley crossing, as a function of the amplification and acceleration factor for graphs of size 100. The mutation rate is set to 10−4. Here, Ns=−1. The black lines are contour lines of the ratio of crossing rate between the network and a well-mixed population and the dots represent individual networks of the families depicted. Network amplification factors empirically determined using equation (9). Network acceleration factors empirically determined using equation (10).

Using the crossing probability Φ and crossing time τ, we can approximate the effective rate of crossing the fitness landscape by (1NμΦ+τ)−1 (Frean et al. 2013). Using this approximation for the rate of evolution, we compare the main families of networks in Fig. 4c (see also Supplementary Section 2). Lattice graphs and detour graphs have the highest rate of crossing the valley, out of all graphs we explored, overtaking strong amplifiers like the star graph. Another graph family of interest is the bipartite graph. These graphs have the lowest fixation time, for a given fixation probability, in the case of a single mutation. However, we show that they can span a wide range of crossing rates, showcasing that a network that is Pareto optimal for a single mutation in terms of fixation time and probability does not guarantee the best performance when crossing a fitness valley. The best-performing graph is inherently dependent on the underlying mutational landscape.

The roles of other network parameters

While the network amplification and acceleration factors play a unifying role in shaping evolutionary dynamics on the network, we also analyze how other network properties, in particular network connectivity and node heterogeneity, shape the fitness landscape crossing probability.

We use networks with preferential attachment, which allow for easy independent tuning of the number of edges and the shape of the degree distribution. We find that strong amplifiers of selection (preferential attachment graphs with low mean and high variance in degree) can increase the probability of fitness valley crossing compared to well-mixed populations (Fig. 5a). This is contrary to expectation, since these graphs suppress deleterious intermediates in the case of sequential fixation. In evolutionary scenarios where the mutation rate is high or the population size is large, the rate of evolutionary dynamics in the population depends not only on the probability of crossing, but also on the expected time to cross. For the same mutational landscape as Figure 5a, the crossing time decreases with increasing mean degree and increases with increasing network node heterogeneity (Fig. 5b, see Supplementary Material for the analytic derivation).

Fig. 5. Amplifiers of selection can promote crossing of fitness valleys, even if the intermediate mutation is deleterious. The dots represent ensemble averages across 107 replicate Monte Carlo simulations per parameter set and network. Here we use preferential attachment graphs and approximations are obtained by fitting the degree standard deviation to the network amplification and acceleration factors (see Supplementary Fig. S3). Lines show the results of substituting the fitting functions into Supplementary equation (54). Network amplification factors empirically determined using equation (9). Network acceleration factors empirically determined using equation (10). Panel a) Here, for dots of the same color, the mean degree of the network is held constant (as in the legend) and we vary the variance of the degree distribution. Population size N=100, μ=10−5, and s=−0.01, such that Ns=−1. Panel b) The crossing time for networks with different variance of the degree distribution. Here, N=100, μ=10−5, and s=−0.01, such that Ns=−1. Lines show the results of substituting the fitting functions into Supplementary equation (64). Panel c) The ratio of the crossing probability for the graph-structured population over the crossing probability of a well-mixed population, as a function of mutation rate. Here, N=100, mean degree is 19 and s=−0.01. Panel d) The ratio of the crossing probability for the graph-structured population over the crossing probability of a well-mixed population, as a function of the selective coefficient of the intermediate mutant. Here, N=100, μ=10−5, and mean degree is 19.

When the mutation rate is low, sequential fixation dominates and the relative crossing probability is equal to the relative fixation probability of the intermediate mutant. However, as the mutation rate increases, stochastic tunneling starts to become a dominant contributor to the crossing probability and the amplifier networks begin to act as suppressors of a deleterious intermediate mutation, increasing the crossing probability compared to the well-mixed model (Fig. 5c). Our approximation holds for a range of mutational landscapes, ranging from fitness valleys to fitness plateaus and fitness hills (Fig. 5d). In the strongly deleterious valley crossing regime, we see the largest deviation from the known theory of single mutation dynamics, since strongly deleterious mutations are unlikely to fix and successful crossing of the landscape depends solely on stochastic tunneling.

Application to bone marrow stem cell population architectures: rates of neoplasm initiation

We apply our model to study the propagation of somatic mutations in the stem cell architectures of the bone marrow to understand spatial factors shaping rates of leukemia initiation. We study a simple two-hit process. While this is a highly simplified model, it nonetheless highlights how the properties of spatial structure, in interplay with details of the mutational landscape, offer predictive power that is impossible to obtain with single mutation models. We assume there exists a set of cancer driver mutations with different effects on cell fitness and that the acquisition of any pair of these genes in the set can lead to neoplasm initiation (Fig. 6a). The second cancer-initiating mutant is assumed to have an extremely high somatic fitness and to fix upon introduction and we assume a negative gamma distribution for the distribution of intermediate mutational effects (Fig. 6b), shifted to the right to account for beneficial mutations (Piganeau and Eyre-Walker 2003):

(7) P(s∣α,β,s0)=βαΓ(α)(s0−s)α−1eβ(s−s0),

where s0 is the maximum fitness effect of the mutation (s<s0), and α and β are parameters of the Gamma distribution.

Fig. 6. A 2-hit model of tumor initiation for bone marrow cellular architectures. Panel a) Illustration of the hematopoietic stem cells niche architecture (for a full description, see Kuo et al. 2021) and a 2-hit model. In this example, there are three genes where any combination of two mutations leads to neoplasm initiation. The letters represent a genotype and the arrows represents mutation paths. Panel b) We assume the intermediate mutations to have different fitness effects, as represented by a shifted negative gamma distribution. Here, the gamma distribution has parameters s0=4×10−4, α∈{2,2.5,4,8}, and β∈{4×10−4,3.2×10−4,2×10−4,1×10−4}, such that the mean of all four distributions presented is equal to (s0−αβ)=−4×10−4. The fraction of beneficial mutations is proportional to parameter β, such that β=4×10−4 corresponds to the darkest purple line and β=1×10−4 corresponds to the lightest purple line. The population size is N=50,000. Panel c) Ratio of cancer initiation rate between the bone marrow tissue cellular architecture and the well-mixed population, as a function of the shape of the fitness distribution and the mutation rate. Ratios are determined using Supplementary equation (54). Here, as above, the gamma distribution has a fixed mean selection strength of s=−4×10−4. The population size is N=50,000. Panel d) The ratio of cancer initiation rate between a 4-regular graph and the well-mixed population, as a function of the shape of the fitness distribution and the mutation rate. Here, the 4-regular graph has 0 triangles. The acceleration factor is 0.66, calculated using Supplementary equation (73). As above, the gamma distribution has a fixed mean selection strength of s=−4×10−4. The population size is N=50,000. Here, the well-mixed line overlaps the stepwise fixation line, since a regular graph does not affect single mutation probability of fixation.

We use published datasets that provide the spatial location of hematopoietic stem and progenitor cells in four samples of mouse tibia (Coutu et al. 2018) and the spatial locations of eight bone marrow samples of CXCL12-abundant reticular cells (which critically modulate hematopoiesis at various levels, including hematopoietic stem cell maintenance), each with two images of two anatomically distinct regions (the diaphysis and the metaphysis), in total 16 cellular populations (Gomariz et al. 2018).

We use previously built cellular networks of the stem cell niches of the bone marrow (Kuo et al. 2021), where every Hematopoietic Stem Cell (HSC) niche constitutes a node in the graph and an edge is added between two nodes if the distance between them is less than a cut-off radius, similar to the generation of a random geometric graph. Kuo et al. (2021) showed that, across a wide variety of parameters and regardless of the birth–death process used, these networks are strong suppressors of selection, potentially delaying mutation accumulation in this tissue.

For a range of biologically relevant mutation rates, and using these previous estimates of amplification and acceleration factor for these networks, we use our analytic results for the rate of fitness landscape crossing to compute the total rate of neoplasm initiation,

(8) Nμ∫p(s)Φ(1N|s,μ,N)ds.

We find that the change in the rate of cancer initiation, compared to a well-mixed population, depends heavily on the assumed fraction of incoming beneficial mutations into the population. At smaller, native mutation rates, the bone marrow mostly suppresses the rate of leukemia initiation, but at higher mutation rates, under environmental carcinogens, for example, the cellular structure can increase the rate of cancer initiation compared to the well-mixed population, for a much wider range in the fraction of incoming beneficial intermediates (Fig. 6c). This is in contrast to lattice-based populations, where the population topology promotes mutation accumulation compared to well-mixed, regardless of the mutation rate and the fitness distribution of the mutation (Fig. 6d). These results show that spatial heterogeneity in tissue architectures can reduce the rates of mutation accumulation and tumor initiation in ways that can be overlooked by previous lattice-based spatial models.

Discussion

Here we study how network structure shapes fitness landscape crossing, compared to well-mixed populations. Our results extend previous well-mixed, lattice, or deme-based models (Bitbol and Schwab 2014; Komarova 2014) and analyze non-symmetric spatial structures of controllable complexity. These previous models show that, by changing the time intermediate mutants can persist until the beneficial final mutation occurs, these spatial structures can help populations cross fitness valleys, compared with well-mixed populations. Here, we recover previous results for lattice-based populations by studying regular graphs with amplification factor α=1 (Bitbol and Schwab 2014; Komarova 2014). However, more complex network structures can change not just the time to fixation, but also the probability of fixation for intermediate variants and the interplay between how the structure affects probabilities and times to fixation and the mutation rate towards the second mutant can give rise to complex evolutionary dynamics, not observed in previous models.

Using simulation and analytic approximations, we show that these complex evolutionary dynamics can be intuitively explained by analyzing how two main network properties, the amplification and acceleration factors, change the expected fate of the intermediate mutant in the population. This intuition can be used to inform on the design of population topologies that optimize rates of crossing fitness landscapes: 1) network-structured populations with low acceleration and amplification factors cross fitness valleys and plateaus more effectively, while 2) network populations with high amplification factors cross fitness hills more effectively. The first observation can be understood and interpreted using Wright’s shifting balance theory, since through either decreased effective selection or increasing hindrance to gene flow, an increase in the relative force of drift leads to an increased window in which the final beneficial mutation can rescue the mutant lineage (Wright et al. 1932).

We find that most network structures rarely purely promote or purely suppress fitness valley crossing, instead transitioning across regimes depending on the selective coefficients of the intermediate mutants. For example, our results show that single-mutation amplifiers, like bipartite graphs and preferential attachment graphs, can increase valley crossing probability for strongly and weakly deleterious intermediate mutations, yet decrease valley crossing for nearly neutral intermediates. Overall, we find seven different network categories that classify evolutionary dynamics of fitness landscape crossing for network-based populations.

For mathematical tractability, obtaining analytical approximations on large, heterogeneous graph structures comes at the cost of simplifying the evolutionary process assumed on the graph. While in this study we compare our results with deme-based or lattice-based models of the Moran process, where all nodes are assumed to be equivalent and probabilities of fixation for the mutant intermediate are not changed, there is also a large literature studying more complex dynamics on populations of demes or on regular graph structures (Barton 1993; Marrec et al. 2021; Yagoobi and Traulsen 2021; Marrec 2023). For example, in Barton (1993) random deme extinction is introduced as an extra component of sampling drift, which can shape intermediate mutant fixation probabilities. Other models, such as Bitbol and Schwab (2014), allow for the uncoupling of migration, birth, and death and, thus, exploration over a broader parameter space. Further work is needed to understand how relaxing the modeling assumptions we make here interplays with heterogenous network structure to shape landscape crossing dynamics.

One key limitation of our analysis is the focus on low selection strength for the intermediate mutant. In this selection regime, the amplifier/suppressor framework, with constant strength of amplification with respect to selection is a valid first-order approximation (McAvoy and Allen 2021). Under a broader range of selection however, there is a wide range of graphs that are piece-wise or transient amplifiers/suppressors, with graphs having selection-dependent amplification factors (Voorhees 2013; Voorhees and Murray 2013; Allen et al. 2020). A notable example is a reducer, a network that consistently decreases the fixation probability compared to a well-mixed model, irrespective of the selection strength of the mutant (Hindersin et al. 2016; Allen et al. 2020). Finer categorizations of amplifiers also exist when considering different mutant initializations, where the mutant does not appear uniformly in the population (Allen et al. 2021). Extending the analytic treatment of crossing dynamics to these selection regimes is challenging and further study is required. Previous work studying single mutation probabilities of fixation has also shown that the type of birth–death process assumed on the network can significantly shape the observed dynamics and further work is needed to understand landscape crossing for other processes, such as death–Birth (Hindersin and Traulsen 2015; Kuo et al. 2021).

Our model can be applied to understand the evolutionary trajectory of cellular systems with complex spatial architectures, and we use an existing dataset of cellular structure to study how the spatial organization of the stem cell niches in the bone marrow shapes rates of leukemia initiation. We find that the relative rate towards neoplasm initiation in this spatial cellular architecture depends critically on the mutation rate, as well as the proportion of beneficial mutations coming into the population. We show that a structure that suppresses cancer initiation can instead transition into one promoting cancer initiation, depending on the distribution of mutational effects. A previous analysis of nearly 50,000 healthy individuals showed that most somatic mutations confer a fitness advantage (Watson et al. 2020). Under a regime of prevalent beneficial driver mutations, we show that the bone marrow architecture suppresses accumulation of these mutations compared to well-mixed populations. However, under exposure to carcinogens (increased rate of mutation Tomlinson et al. 1996; Nowak et al. 2002; Pino and Chung 2010), the architecture of the bone marrow can instead increase rates of neoplasm initiation. As spatial datasets of cellular and molecular organization become increasingly available, linking these datasets with theoretical models of evolutionary spread will be crucial for a rigorous understanding of how complex spatial structure shapes evolutionary dynamics at these levels of organization.

Materials and methods

List of network families used in the study

Linking network topology to evolutionary dynamics and understanding which network properties shape rates of evolution is complicated by the fact that these properties are often correlated, hard to tune independently and differ across many network families (Kuo et al. 2021). In this study, we explore both well-known network families using built-in generators from NetworkX (Hagberg et al. 2008) and also design graphs that allow us to tune properties independently, as detailed below.

k-regular graphs

A k-regular graph is a graph where each node has the same number of neighbors, k. We generate k-regular graphs using built-in generators from NetworkX (Hagberg et al. 2008), with k∈{3,5,10,20}. For some of these graphs, to generate Fig. 4, we further tune the number of triangles and change their acceleration factor (while keeping mean degree constant at k), using methods previously outlined in Kuo and Carja (2024). The algorithm ensures uniform sampling across all connected regular graphs with the target fraction of triangles.

Erdős Rényi random networks

The Erdős Rényi model starts with a set of N isolated nodes, and connects each pair of two nodes with probability p, independently. We generate Erdős Rényi random networks using built-in generators from NetworkX (Hagberg et al. 2008).

Preferential attachment graphs

For graphs with preferential attachment, nodes are added sequentially starting from one initial node until the population reaches size N. Each new node is added to the network and connected to other individuals with a probability proportional to the individual’s current degree to the power of a given parameter β. Using k and β, this family of graphs allows for straightforward independent tuning of mean and variance in degree. We generate these networks using NetworkX (Hagberg et al. 2008) with k∈{3,5,20}, and β∈[−3,3].

Small world networks

We generate small world networks using built-in generators from NetworkX (Hagberg et al. 2008) and the Watts–Strogatz model with parameters: number of nodes N, the mean degree k, and rewire probability p. We use N=100,  k∈{8,12,16}, and p∈[0,1].

Bipartite graphs

Bipartite graphs have two node sets n1 and n2. Edges in the graph only connect nodes from opposite sets. We generate bipartite graphs using built-in generators from NetworkX (Hagberg et al. 2008). We vary n1 from 1 to 50, and n2=100−n1. For each parameter set, there exists one single complete bipartite graph.

Random geometric graphs

The random geometric graph model places N nodes at random following a probability distribution. Two nodes are joined by an edge if the distance between the nodes is below a predefined cut-off radius. For the graphs in Fig. 4, we place nodes in 2-dimensional space. For the x-position, 50 nodes are drawn from a normal distribution N(0,0.25), and 50 nodes are drawn from N(3,0.25). The y-position is drawn from a normal distribution N(0,0.25). We start with the lowest cut-off radius, resulting in a connected graph, and increase the cut-off radius until the complete graph is formed.

Detour graphs

A detour graph is formed by starting with a complete graph of size n1 and replacing one of the edges with a path of length n2+1≥2 (Möller et al. 2019). To generate the detours in Fig. 4, we vary n1 from 3 to 99, and n2=100−n1.

Shortcut graphs

We design the shortcut graphs that allow us to vary amplification, while keeping acceleration equal to one used in Fig. 2f and Fig. 4, by starting with a detour graph and adding edges (shortcuts) between two nodes. We add 1–255 edges to the detour graph with n1=90 and n2=10.

Star and star-like graphs

The star graph consists of one center node connected to N−1 outer nodes. We generate star graphs using built-in generators from NetworkX (Hagberg et al. 2008). A star-like graph is constructed by adding random edges to the star graph (Tkadlec et al. 2019). We generate 20 graphs each for 50, 100, 200, and 300 random edges added to a star-graph of size 100.

Regular islands

A regular island is the result of connecting two k-regular graphs. A random edge is picked from each k-regular graph, and edge swap is performed to connect the two k-regular graphs.

Computing the amplification and acceleration parameters for a given network

Here we show that the amplification and acceleration factors of a network, computed from single mutation analyses, if used together, can be predictive of the landscape crossing behavior of network-structured populations.

In previous work, we showed how to analytically compute the amplification and acceleration factors (Kuo et al. 2021; Kuo and Carja 2024). The approximation for the amplification factor from Kuo et al. (2021) is very accurate for amplifiers and can slightly deviate from the true amplification factor for suppressors (we also write that approximation as equation (37) in the Supplementary Material). Similarly, the approximations described in Kuo and Carja (2024) for computing the acceleration factor are very accurate if triangles in the network are distributed equally amongst the triplet types (we write this approximation as equation (38) in the Supplementary Material).

While these approximations provide intuition into how other network parameters affect the amplification and acceleration of a network, since they can slightly deviate from the true number under certain scenarios, here we determine them empirically, for increased accuracy.

Using the Birth–death process, and starting with one single mutant with fitness (1+s) invading a population of wild-type individuals with fitness 1, the amplification factor is computed using the definition from (Lieberman et al. 2005) and solving

(9) Φ=1−(1+s)−α1−(1+s)−αN

for α, where Φ is the mutant’s fixation probability.

For graphs in this paper, we estimate α using Ns=1. When graph size is N=100, this equates to s=0.01. This value is large enough such that the difference in fixation probability compared to the well-mixed model is big enough, while also small enough to be in the constant regime. The amplification factor for a given network is estimated once and used to study landscape crossing with varying selection on the intermediate mutant.

For the acceleration factor, we use the following definition

(10) λ=Tfix,wm(αs)Tfix,graph(s),

which is the ratio of conditional mean fixation time for an equivalent mutant on the well-mixed population (with selection coefficient αs) and the one on the network population with mutant selection coefficient s (Kuo and Carja 2024). For the graphs presented here, we evaluate λ at s=0, since there is a closed-form expression for Tfix,wm(s=0), (N−1) generations (Ewens 2004). For weak s,  λ varies no more than O(s).

Building the cellular networks of the stem cell niches of the bone marrow

The samples vary in dimensions, number of cells, and segmentation techniques. We normalize the data by expressing the distance in units of the average shortest distance between pairs of cells (62.72 μm for hematopoietic stem cells). We use networks generated using cut-off distance of 15, which is the closest to our estimated biological interaction range. We interpret the cut-off distance as the maximum distance an HSC could travel in its entire lifespan. Live-animal tracking of individual hematopoietic stem cells in their niche showed MFG cells, a largely quiescent population with long-term self-renewal capability, displacing an average distance of 8.69μm in a 2.5 hour period (Christodoulou et al. 2020). HSCs have median replication time (the time when 50% of HSCs have divided) of 1.7 weeks (Abkowitz et al. 2000). During homeostasis, the rate of replication should balance the rate of depletion. This leads to the estimated interaction range of 1,028μm which corresponds to 16.4× the distance between the shortest pairs. More details on network generation and robustness results are presented in Kuo et al. (2021).

The amplification and acceleration factors of the bone marrow network structures are then calculated using Monte Carlo simulation of single mutation fixation of 107 replicate simulations at Ns=0.01. The two factors are used to derive the analytical approximation for the crossing probability in these cellular populations, shown together with simulations in Supplementary Fig. S6a. Although the image data provide precise spatial description of the hematopoietic stem cell niches, these only represent a small portion of the HSC population. Previous models estimate the number of hematopoietic stem cells that are actively making white blood cells at any one time to be in the range of 50,000–2,00,000 (Lee-Six et al. 2018). Here, we assume the overall spatial distribution of the HSC population across the bone marrow is identical to our sample. We use linear regression to extrapolate the amplification factor to the lower bound population of 50,000 (Supplementary Fig. S6b). The acceleration factors inferred from data are close to 1, so we assume the larger population also has acceleration factor of 1.

Supplementary Material

iyae055_Supplementary_Data

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.

Data availability

The graphs used in this study, information about the C++ and Python code used to simulate and analyze the model presented here are available at https://github.com/yangpingkuo/Evolutionary-graph-theory-beyond-single-mutation-dynamics. Supplementary material available at GENETICS online.

Funding

We gratefully acknowledge support from the NIH National Institute of General Medical Sciences (award no. R35GM147445), the United States-Israel Binational Science Foundation (award no. 2019266) and from the NIH T32 training grant (no. T32 EB009403).
==== Refs
Literature cited

Abkowitz  JL, Golinelli  D, Harrison  DE, Guttorp  P. 2000. In vivo kinetics of murine hemopoietic stem cells. Blood J Am Soc Hematol. 96 (10 ):3399–3405.
Acar  A, Nichol  D, Fernandez-Mateos  J, Cresswell  GD, Barozzi  I, Hong  SP, Trahearn  N, Spiteri  I, Stubbs  M, Burke  R, et al. 2020. Exploiting evolutionary steering to induce collateral drug sensitivity in cancer. Nat Commun. 11 (1 ):1923. doi:10.1038/s41467-020-15596-z 32317663
Adlam  B, Chatterjee  K, Nowak  MA. 2015. Amplifiers of selection. Proc R Soc A: Math Phys Eng Sci. 471 (2181 ):20150114. doi:10.1098/rspa.2015.0114
Allen  B, Gore  J, Nowak  MA. 2013. Spatial dilemmas of diffusible public goods. Elife. 2 :e01169. doi:10.7554/eLife.01169 24347543
Allen  B, Sample  C, Jencks  R, Withers  J, Steinhagen  P, Brizuela  L, Kolodny  J, Parke  D, Lippner  G, Dementieva  YA. 2020. Transient amplifiers of selection and reducers of fixation for death-birth updating on graphs. PLoS Comput Biol. 16 (1 ):e1007529. doi:10.1371/journal.pcbi.1007529 31951612
Allen  B, Sample  C, Steinhagen  P, Shapiro  J, King  M, Hedspeth  T, Goncalves  M. 2021. Fixation probabilities in graph-structured populations under weak selection. PLoS Comput Biol. 17 (2 ):e1008695. doi:10.1371/journal.pcbi.1008695 33529219
Asratian  AS, Denley  TM, Häggkvist  R. 1998. Bipartite Graphs and Their Applications. Vol. 131 . Cambridge: Cambridge University Press.
Barabási  A-L, Albert  R. 1999. Emergence of scaling in random networks. Science. 286 (5439 ):509–512. doi:10.1126/science.286.5439.509 10521342
Barton  NH . 1993. The probability of fixation of a favoured allele in a subdivided population. Genet Res. 62 (2 ):149–157. doi:10.1017/S0016672300031748
Bitbol  AF, Schwab  DJ. 2014. Quantifying the role of population subdivision in evolution on rugged fitness landscapes. PLoS Comput Biol. 10 (8 ):e1003778. doi:10.1371/journal.pcbi.1003778 25122220
Burch  CL, Chao  L. 1999. Evolution by small steps and rugged landscapes in the RNA virus ϕ6. Genetics. 151 (3 ):921–927. doi:10.1093/genetics/151.3.921 10049911
Carja  O, Liberman  U, Feldman  MW. 2014. Evolution in changing environments: modifiers of mutation, recombination, and migration. Proc Natl Acad Sci USA. 111 (50 ):17935–17940. doi:10.1073/pnas.1417664111 25427794
Christodoulou  C, Spencer  JA, Yeh  S-CA, Turcotte  R, Kokkaliaris  KD, Panero  R, Ramos  A, Guo  G, Seyedhassantehrani  N, Esipova  TV, et al. 2020. Live-animal imaging of native haematopoietic stem and progenitor cells. Nature. 578 (7794 ):278–283. doi:10.1038/s41586-020-1971-z 32025033
Coutu  DL, Kokkaliaris  KD, Kunz  L, Schroeder  T. 2018. Multicolor quantitative confocal imaging cytometry. Nat Methods. 15 (1 ):39–46. doi:10.1038/nmeth.4503 29320487
Durrett  R, Moseley  S. 2015. Spatial Moran models I. Stochastic tunneling in the neutral case. Ann Appl Probab: off J Inst Math Stat. 25 (1 ):104. doi:10.1214/13-AAP989
Ewens  WJ . 2004. Mathematical Population Genetics: theoretical Introduction. Vol. 27 . New York: Springer.
Frean  M, Rainey  PB, Traulsen  A. 2013. The effect of population structure on the rate of evolution. Proc R Soc B: Biol Sci. 280 (1762 ):20130211. doi:10.1098/rspb.2013.0211
Gomariz  A, Helbling  PM, Isringhausen  S, Suessbier  U, Becker  A, Boss  A, Nagasawa  T, Paul  G, Goksel  O, Székely  G, et al. 2018. Quantitative spatial analysis of haematopoiesis-regulating stromal cells in the bone marrow microenvironment by 3D microscopy. Nat Commun. 9 (1 ):2532. doi:10.1038/s41467-018-04770-z 29955044
Hagberg  AA, Schult  DA, Swart  PJ. 2008. Exploring network structure, dynamics, and function using networkx. In: Varoquaux G, Vaught T, Millman J, editors, Proceedings of the 7th Python in Science Conference (SCIPY 08), Pasadena (CA) USA. p. 11–15.
Hindersin  L, Traulsen  A. 2014. Counterintuitive properties of the fixation time in network-structured populations. J R Soc Interface. 11 (99 ):20140606. doi:10.1098/rsif.2014.0606 25142521
Hindersin  L, Traulsen  A. 2015. Most undirected random graphs are amplifiers of selection for birth-death dynamics, but suppressors of selection for death-birth dynamics. PLoS Comput Biol. 11 (11 ):e1004437. doi:10.1371/journal.pcbi.1004437 26544962
Hindersin  L, Werner  B, Dingli  D, Traulsen  A. 2016. Should tissue structure suppress or amplify selection to minimize cancer risk?  Biol Direct. 11 (1 ):1–11. doi:10.1186/s13062-016-0140-7 26738889
Huang  S . 2013. Genetic and non-genetic instability in tumor progression: link between the fitness landscape and the epigenetic landscape of cancer cells. Cancer Metastasis Rev. 32 (3-4 ):423–448. doi:10.1007/s10555-013-9435-7 23640024
Iwasa  Y, Michor  F, Nowak  MA. 2004. Stochastic tunnels in evolutionary dynamics. Genetics. 166 (3 ):1571–1579. doi:10.1534/genetics.166.3.1571 15082570
Jain  K, Krug  J. 2007. Deterministic and stochastic regimes of asexual evolution on rugged fitness landscapes. Genetics. 175 (3 ):1275–1288. doi:10.1534/genetics.106.067165 17179085
Jiang  C, Chen  Y, Liu  KJR. 2014. Evolutionary dynamics of information diffusion over social networks. IEEE Trans Signal Process. 62 (17 ):4573–4586. doi:10.1109/TSP.2014.2339799
Kimura  M, Weiss  GH. 1964. The stepping stone model of population structure and the decrease of genetic correlation with distance. Genetics. 49 (4 ):561–576. doi:10.1093/genetics/49.4.561 17248204
Komarova  NL . 2014. Spatial interactions and cooperation can change the speed of evolution of complex phenotypes. Proc Natl Acad Sci USA. 111 (Supplement 3):10789–10795. doi:10.1073/pnas.1400828111 25024187
Komarova  NL, Sengupta  A, Nowak  MA. 2003. Mutation–selection networks of cancer initiation: tumor suppressor genes and chromosomal instability. J Theor Biol. 223 (4 ):433–450. doi:10.1016/S0022-5193(03)00120-6 12875822
Kuo  YP, Arrieta  CN, Carja  O. 2021. A theory of evolutionary dynamics on any complex spatial structure. bioRxiv. 10.1101/2021.02.07.430151, preprint: not peer reviewed.
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. doi:10.1371/journal.pcbi.1011905 38489353
Kvitek  DJ, Sherlock  G. 2011. Reciprocal sign epistasis between frequently experimentally evolved adaptive mutations causes a rugged fitness landscape. PLoS Genet. 7 (4 ):e1002056. doi:10.1371/journal.pgen.1002056 21552329
Lande  R . 1979. Effective deme sizes during long-term evolution estimated from rates of chromosomal rearrangement. Evolution. 33 (1 ):234–251. doi:10.2307/2407380 28568063
Lee-Six  H, Øbro  NF, Shepherd  MS, Grossmann  S, Dawson  K, Belmonte  M, Osborne  RJ, Huntly  BJ, Martincorena  I, Anderson  E, et al. 2018. Population dynamics of normal human blood inferred from somatic mutations. Nature. 561 (7724 ):473–478. doi:10.1038/s41586-018-0497-0 30185910
Leventhal  GE, Hill  AL, Nowak  MA, Bonhoeffer  S. 2015. Evolution and emergence of infectious diseases in theoretical and real-world networks. Nat Commun. 6 (1 ):6101. doi:10.1038/ncomms7101 25592476
Lieberman  E, Hauert  C, Nowak  MA. 2005. Evolutionary dynamics on graphs. Nature. 433 (7023 ):312–316. doi:10.1038/nature03204 15662424
Maciejewski  W, Fu  F, Hauert  C. 2014. Evolutionary game dynamics in populations with heterogenous structures. PLoS Comput Biol. 10 (4 ):e1003567. doi:10.1371/journal.pcbi.1003567 24762474
Marrec  L . 2023. Quantifying the impact of genotype-dependent gene flow on mutation fixation in subdivided populations. bioRxiv. 10.1101/2023.11.29.569213, preprint: not peer reviewed
Marrec  L, Lamberti  I, Bitbol  AF. 2021. Toward a universal model for spatially structured populations. Phys Rev Lett. 127 (21 ):218102. doi:10.1103/PhysRevLett.127.218102 34860074
Maruyama  T . 1970a. Effective number of alleles in a subdivided population. Theor Popul Biol. 1 (3 ):273–306. doi:10.1016/0040-5809(70)90047-X 5527634
Maruyama  T . 1970b. On the fixation probability of mutant genes in a subdivided population. Genet Res (Camb). 15 (2 ):221–225. doi:10.1017/S0016672300001543
McAvoy  A, Allen  B. 2021. Fixation probabilities in evolutionary dynamics under weak selection. J Math Biol. 82 (3 ):1–41. doi:10.1007/s00285-021-01568-4 33475794
Möller  M, Hindersin  L, Traulsen  A. 2019. Exploring and mapping the universe of evolutionary graphs identifies structural properties affecting fixation probability and time. Commun Biol. 2 (1 ):1–9. doi:10.1038/s42003-019-0374-x 30740537
Nowak  MA, Komarova  NL, Sengupta  A, Jallepalli  PV, Shih  I-M, Vogelstein  B, Lengauer  C. 2002. The role of chromosomal instability in tumor initiation. Proc Natl Acad Sci USA. 99 (25 ):16226–16231. doi:10.1073/pnas.202617399 12446840
Ohtsuki  H, Hauert  C, Lieberman  E, Nowak  MA. 2006. A simple rule for the evolution of cooperation on graphs and social networks. Nature. 441 (7092 ):502–505. doi:10.1038/nature04605 16724065
Ohtsuki  H, Nowak  MA. 2006. The replicator equation on graphs. J Theor Biol. 243 (1 ):86–97. doi:10.1016/j.jtbi.2006.06.004 16860343
Ohtsuki  H, Pacheco  JM, Nowak  MA. 2007. Evolutionary graph theory: breaking the symmetry between interaction and replacement. J Theor Biol. 246 (4 ):681–694. doi:10.1016/j.jtbi.2007.01.024 17350049
Penrose  M . 2003. Random Geometric Graphs. Vol. 5 . Oxford, UK: Oxford University Press.
Piganeau  G, Eyre-Walker  A. 2003. Estimating the distribution of fitness effects from DNA sequence data: implications for the molecular clock. Proc Natl Acad Sci USA. 100 (18 ):10335–10340. doi:10.1073/pnas.1833064100 12925735
Pino  MS, Chung  DC. 2010. The chromosomal instability pathway in colon cancer. Gastroenterology. 138 (6 ):2059–2072. doi:10.1053/j.gastro.2009.12.065 20420946
Pollak  E . 1966. On the survival of a gene in a subdivided population. J Appl Probab. 3 (1 ):142–155. doi:10.2307/3212043
Poncela  J, Gómez-Gardenes  J, Floría  LM, Moreno  Y. 2007. Robustness of cooperation in the evolutionary prisoner’s dilemma on complex networks. New J Phys. 9 (6 ):184. doi:10.1088/1367-2630/9/6/184
Rogers  ZN, McFarland  CD, Winters  IP, Seoane  JA, Brady  JJ, Yoon  S, Curtis  C, Petrov  DA, Winslow  MM. 2018. Mapping the in vivo fitness landscape of lung adenocarcinoma tumor suppression in mice. Nat Genet. 50 (4 ):483–486. doi:10.1038/s41588-018-0083-2 29610476
Salehi  S, Kabeer  F, Ceglia  N, Andronescu  M, Williams  MJ, Campbell  KR, Masud  T, Wang  B, Biele  J, Brimhall  J, et al. 2021. Clonal fitness inferred from time-series modelling of single-cell cancer genomes. Nature. 595 (7868 ):585–590. doi:10.1038/s41586-021-03648-3 34163070
Santos  FC, Pacheco  JM, Lenaerts  T. 2006. Evolutionary dynamics of social dilemmas in structured heterogeneous populations. Proc Natl Acad Sci USA. 103 (9 ):3490–3494. doi:10.1073/pnas.0508201103 16484371
Slatkin  M . 1981. Fixation probabilities and fixation times in a subdivided population. Evolution. 35 (3 ): 477–488. doi:10.2307/2408196 28563585
Su  Q, Allen  B, Plotkin  JB. 2022. Evolution of cooperation with asymmetric social interactions. Proc Natl Acad Sci USA. 119 (1 ):e2113468118. doi:10.1073/pnas.2113468118 34983850
Tkadlec  J, Pavlogiannis  A, Chatterjee  K, Nowak  MA. 2019. Population structure determines the tradeoff between fixation probability and fixation time. Communications Biology. 2 (1 ):1–8. doi:10.1038/s42003-019-0373-y 30740537
Tkadlec  J, Pavlogiannis  A, Chatterjee  K, Nowak  MA. 2020. Limits on amplifiers of natural selection under death-birth updating. PLoS Comput Biol. 16 (1 ):e1007494. doi:10.1371/journal.pcbi.1007494 31951609
Tomlinson  IP, Novelli  M, Bodmer  W. 1996. The mutation rate and cancer. Proc Natl Acad Sci USA. 93 (25 ):14800–14803. doi:10.1073/pnas.93.25.14800 8962135
Vogelstein  B, Papadopoulos  N, Velculescu  VE, Zhou  S, Diaz  Jr  LA, Kinzler  KW. 2013. Cancer genome landscapes. Science. 339 (6127 ):1546–1558. doi:10.1126/science.1235122 23539594
Voorhees  B . 2013. Birth–death fixation probabilities for structured populations. Proc R Soc A: Math Phys Eng Sci. 469 (2153 ):20120248. doi:10.1098/rspa.2012.0248
Voorhees  B, Murray  A. 2013. Fixation probabilities for simple digraphs. Proc R Soc A: Math Phys Eng Sci. 469 (2154 ):20120676. doi:10.1098/rspa.2012.0676
Watson  CJ, Papula  A, Poon  GY, Wong  WH, Young  AL, Druley  TE, Fisher  DS, Blundell  JR. 2020. The evolutionary dynamics and fitness landscape of clonal hematopoiesis. Science. 367 (6485 ):1449–1454. doi:10.1126/science.aay9333 32217721
Watts  DJ, Strogatz  SH. 1998. Collective dynamics of small world networks. Nature. 393 (6684 ):440–442. doi:10.1038/30918 9623998
Waxman  BM . 1988. Routing of multipoint connections. IEEE J Sel Areas Commun. 6 (9 ):1617–1622. doi:10.1109/49.12889
Weinreich  DM, Chao  L. 2005. Rapid evolutionary escape by large populations from local fitness peaks is likely in nature. Evolution. 59 (6 ):1175–1182.16050095
Weissman  DB, Desai  MM, Fisher  DS, Feldman  MW. 2009. The rate at which asexual populations cross fitness valleys. Theor Popul Biol. 75 (4 ):286–300. doi:10.1016/j.tpb.2009.02.006 19285994
Whitlock  MC . 2003. Fixation probability and time in subdivided populations. Genetics. 164 (2 ):767–779. doi:10.1093/genetics/164.2.767 12807795
Whitlock  MC, Barton  N. 1997. The effective size of a subdivided population. Genetics. 146 (1 ):427–441. doi:10.1093/genetics/146.1.427 9136031
Wood  LD, Parsons  DW, Jones  S, Lin  J, Sjoblom  T, Leary  RJ, Shen  D, Boca  SM, Barber  T, Ptak  J, et al. 2007. The genomic landscapes of human breast and colorectal cancers. Science. 318 (5853 ):1108–1113. doi:10.1126/science.1145720 17932254
Wright  S . 1932. The roles of mutation, inbreeding, crossbreeding and selection in evolution. In: Proceedings of the XI International Congress of Genetics. Vol. 8 . p. 209–222.
Wright  S . 1943. Isolation by distance. Genetics. 28 (2 ):114–138. doi:10.1093/genetics/28.2.114 17247074
Yagoobi  S, Sharma  N, Traulsen  A. 2023. Categorizing update mechanisms for graph-structured metapopulations. J R Soc Interface. 20 (200 ):20220769. doi:10.1098/rsif.2022.0769 36919418
Yagoobi  S, Traulsen  A. 2021. Fixation probabilities in network structured meta-populations. Sci Rep. 11 (1 ):17979. doi:10.1038/s41598-021-97187-6 34504152
