Introduction

Genetic drift is defined as the random changes in frequency of genetic variants (Crow and Kimura 1970; Hartl and Clark 1997; Lynch, et al. 2016). Much, or even most, of DNA sequence evolution is governed by the random drift of neutral mutations. Nevertheless, even non-neutral variants are subjected to genetic drift as advantageous mutations are not deterministically fixed. If the strength of genetic drift is under-estimated, random changes that are missed in the analysis would result in the over-estimation of other evolutionary forces of greater biological interest, including selection, mutation, migration, meiotic drive and so on. In particular, the conclusion of pervasive selection at the molecular level has been suspected to be due to the failure to fully account for genetic drift (Crow and Kimura 1970; Fu 1997; Li 1997; Charlesworth 2009; Lynch, et al. 2016; Chen, Yang, et al. 2022).

Hagedoorn and Hagedoorn (1921) may be the first to suggest that random forces (accidents, for example) could impact long-term evolution. In the same decade, models of genetic drift were developed, obviously with some disagreements, along two lines. The one approach that became widely adopted is the Wright-Fisher (WF) model (Crow and Kimura 1970; Hartl and Clark 1997), formalized by Fisher (1930) and Wright (1931). They define genetic drift to be solely due to the random sampling of genes. Thus, in generation transition, the variance in the frequency of a variant allele, x, would be

where N is the population size. The effective population size, N_e, is a function of varying N’s in space, time and demography. For example, when N fluctuates, N_e would be the harmonic mean of N’s across time (Crow and Kimura 1970; Caballero 1994; Charlesworth 2009). Weakly or strongly, N_e should still reflect the actual population size, N. In the WF model, genetic drift is governed by 1/ N_e in haploid populations, or 1/(2N_e) in diploid populations.

At the same time, JBS Haldane introduced the branching process to model genetic drift (Haldane 1927; Haldane 1932). In the Haldane model, each gene copy is independently transmitted to K descendants with the mean and variance of E(K) and V(K). In its essence, genetic drift is V(K), i.e., the variance in transmission among neutral variants. There would be no drift if V(K) = 0. To scale drift at the population level, the variance becomes V(K)/N (Kimura and Crow (1963); Chen, et al. (2017); see also the Supplementary Information). In short, V(K) is the reproductive component of genetic drift while 1/N is the sampling component or the scaling factor.

It should be noted that the WF model has been modified in many directions (Wright 1969; Gillespie 1975; Eldon and Wakeley 2006; Chen, et al. 2017; Sackman, et al. 2019) and, in particular, with the inclusion of V(K) (Kimura and Crow 1963; Sjodin, et al. 2005; Der, et al. 2011; Cannings 2016). These modifications share the basic features as well as the difficulties in resolving the paradoxes encountered by the WF model (see Results). Such difficulties stem from the sampling step that is tied to a pre-determined population size. These difficulties are shared by other models with population sampling, such as the Moran model (Moran 2008). In Discussion, we will present an overview of the modified WF models where these paradoxes will be discussed in relation to both the standard and modified WF models.

Among the paradoxes, a most curious one is genetic drift in relation to changing population size. The WF model dictates stronger drift when N is smaller. However, when N is very small (say, N = 10 in a bacteria or yeast population) and increases to 20, 40, 80 and so on, there is in fact little drift in this exponential phase. Drift will intensify when N grows to near the carrying capacity. This trend is the exact opposite of the tenet of the WF model. The second paradox concerns the genetic drift of sex-linked genes in relation to sex-dependent breeding successes (Charlesworth and Charlesworth 2000; Bachtrog 2013; Cortez, et al. 2014; Wilson Sayres, et al. 2014; Makova, et al. 2024). A third paradox is about drift strength when selection is also at work. A new mutation with a selective advantage would always be fixed if there is no random drift. With drift, fixation is probabilistic but the probability does not depend on N or N_e. This curious property echoes the view that N or N_e is a scaling factor of drift, which is determined mainly by V(K) (see Results).

The power of the Haldane model in resolving these paradoxes is the focus of this study. Furthermore, the companion study (Wang et al. 2024) addresses a fourth paradox – the evolution of multi-copy gene systems, whereby evolution proceeds in two stages - between as well as within individuals. Multi-copy gene systems including viruses, mitochondria, transposons and ribosomal genes, are hence beyond the power of the WF models. Strictly speaking, even diploidy is also multi-copy systems (Silver 1985; Wu, et al. 1988; Wu, et al. 1989; Lindholm, et al. 2016; Courret, et al. 2019).

In this study, we present the integration of the Haldane and WF models into the WFH model. It is essentially the extension of the Haldane model with the capacity for N regulation. In this integration, the wealth of results obtained by the WF model are often adequate approximations for the more exact results. Most important of all, the WFH model can resolve the many paradoxes of molecular evolution.

Results

In the WFH model, genetic drift is defined as in the Haldane model of branching process (see Eq. (4) of Chen, et al. (2017)). Let K be the number of progenies receiving the gene copy of interest from a parent; K would follow the same distribution for all neutral variants. In haploids, K is equivalent to the progeny number of each parent. In diploids, K is the transmission success of each gene copy, which is separately tracked. By using 2N for the gene copy number in the population, the WF as well as Haldane model treats each diploid individual as two-haploids. By the Taylor approximations of the ratio of two variables (Kendall, et al. 2006), we obtain the approximation as

where x is the variant frequency (see Supplementary Information). This is the same approximation, when V(K) is not very large, as obtained by the WF model (Kimura and Crow 1964; see Chen et al. 2017 for details). For Eq. (2), we provide further extensive simulations on the accuracy in both the fixation probability and fixation time of neutral variants (see Supplementary Information). The equivalence makes a simple point that the WF model represent special cases of the WFH model, and the results are often, but not always, shared.

In the WFH model, there would be no genetic drift if V(K) = 0, regardless of the value of N. When E(K) = 1 and N stays constant, we obtain

Eq. (2) or Eq. (3) as applied to the WF model (Kimura and Crow 1963; Crow and Denniston 1988; Chen, et al. 2017) is biologically dubious. This is because the sampling scheme of the WF model dictates Poisson distribution for K, hence, V(K) = E(K). To permit V(K) ≠ E(K), the WF model would have to be modified to resemble the Haldane model (Kimura and Crow 1963; Caballero 1994). In Discussion, we shall analyze this “fix” of the WF model that still fails to resolve the paradoxes raised in this study.

In Table 1, the data of progeny number from the literature indeed show over-dispersion in all taxa surveyed, i.e., V(K) > E(K). In fact, V(K) > 10 E(K) can often be found, e.g., the ratio of V(K)/E(K) among males of the mandrill (Mandrillus sphinx), the great reed warbler (Acrocephalus arundinaceus), and the rhesus macaque (Macaca mulatta) at 19.5, 12.6 and 11.3, respectively (Hasselquist 1995; Setchell, et al. 2005; Dubuc, et al. 2014). We have also found that V(K) may be far larger than E(K) in other biological systems such as viruses (Ruan, Luo, et al. 2021; Ruan, Wen, et al. 2021; Hou, et al. 2023).

Table 1.

I. The paradox of genetic drift in relation to N changes (Paradox of changing N)

1. Empirical demonstration and simulation

This “paradox of changing N” is that “genetic drift increases when N increases”. Fig. 1a shows the results in a cell culture in the exponential growth phase where each cell doubles every 13 hours. Hence, V(K) ∼ 0 with almost no drift when N is as small as < 50. Genetic drift is expected to increase as N approaches the carrying capacity. The paradox is exemplified through computer simulations in Fig. 1b-c. These panels show the drift pattern in the Haldane model and the WF model, respectively, when the populations are growing as the logistic growth model (see Methods). In the WF model, the drift is strong initially but weakened as N increases. The pattern is reversed in the Haldane model. The contrast between the two panels is the paradox.

Fig. 1.

To verify the simulations of Fig. 1b-c, we study the growth of yeast populations under laboratory conditions where cells appear to increase following the logistic growth curve. As shown in Fig. 1d, the lower portion of the S-curve portrays the exponential growth in a low-density population, while the upper curve indicates a slowdown in growth near the population’s carrying capacity. We directly recorded the individual output (K) of each yeast cell in the early (n = 25) and late (n = 65) stage of population growth under real-time high-content fluorescence imaging (Methods).

It’s evident that V(K) is decoupled from E(K). In the early stage, E(K) is 10.04 with V(K) = 0.60 over five one-hour intervals. In the high-density phase, E(K) decreases to 2.83 but V(K) increases to 1.09 in four intervals (Fig. 1e-f). Most interestingly, during the late stage of population growth, V(ΔK) (adjusted to the four-hour time span for comparison; see Methods) indeed increases as N increases (Supplementary Table 3). At the last time point, the population is closest to the carrying capacity and V(ΔK) experiences a substantial leap from around 1.12 to 1.82. Overall, the variance in progeny number V(K) in a high-density population is always greater that in a low-density population (Fig. 1g). Indeed, V(K) may exhibit an increase with the increase in N, sometimes outpaces the latter.

Measurements of these two growth phases demonstrate the “paradox of changing N” by showing i) there is almost no drift when N is very small; ii) as N increases, V(K), hence, genetic drift would also increase. Furthermore, V(K) may increase sharply when N is close to the carrying capacity. This trend appears to echo the field observations by Coltman, et al. (1999) who show that V(K) increases more than 5-fold when N increases only 2.5-fold in the United Kingdom (UK) population of Soay rams (Ovis aries).

2. The density-dependent Haldane (DDH) model

The “paradox of changing N” can be explained by N_e = N/V(K) if the numerator and denominator move in the same direction. Most interestingly, V(K) may increase more than N itself. This has been shown when V(K) is close to zero while N is growing exponentially. In addition, when N is very near the carrying capacity, even a small increase in N would intensify the competition and inflate V(K). This is shown in the yeast experiment of Fig. 1d and supported by Coltman, et al. (1999) field study of ram populations in UK. Thus, the relationship between N and V(K) may realistically lead to this paradox.

We now extend the Haldane model by incorporating density-dependent regulation of N in order to define the conditions of this paradox. The model developed here is on an ecological, rather than an evolutionary time scale. The simplest model of N regulation is the logistic equation:

with C_k being the carrying capacity. In the ecological time scale, a changing N would mean that the population is departing from C_k or moving toward it (Fig. 2a). (As will be addressed in Discussion, changing N in the WF model is depicted in Fig. 2b whereby N is at C_k. C_k may evolve too, albeit much more slowly in an evolutionary time scale.) Here, we consider changes in N in the ecological time scale. Examples include the exponential growth when N ≪ C_k, seasonal fluctuation in N (Frankham 1995; Charlesworth 2009) and competition modeled by the Lotka-Volterra equations (Bomze 1983).

Fig. 2.

Details of the DDH model are given in the Supplementary Information. A synopsis is given here: We consider a non-overlapping haploid population with two neutral alleles. The population size at time t is N_t. We assume that expected growth rate E(K) is greater than 1 when N < C_k and less than 1 when N > C_k, as defined by Eq. (5) below:

The slope of E(K) vs. N (i.e., the sensitive of growth rate to changes in population size), as shown in Fig 3a, depends on z. To determine the variance V(K), we assume that K follows the negative binomial distribution whereby parents would suffer reproduction-arresting injury with a probability of p_t at each birthing (Supplementary Information). Accordingly, V(K) can then be expressed as

By Eq. (6), the ratio of V(K)/E(K) could be constant, decrease or increase with the increase of population size. With E(K) and V(K) defined, we could obtain the effective population size by substituting Eq. (5) and Eq. (6) into Eq. (3).

Eq. (7) presents the relationship between effective population size (N_e) and the population size (N) as shown in Fig. 3. The density-dependent E(K) could regulate N with different strength (Fig. 3a). The steeper the slope in Fig. 3a, the stronger the regulation.

Fig. 3.

The main results of this study are depicted in Fig. 3b. First, with no or little regulation of N, N_e and N are strongly correlated. The green dashed lines portray the conventional view of decreasing drift with increasing N. Second, under a particular strength of N regulation, the red line shows no dependence of N_e on N, meaning that genetic drift is independent of N. Finally, as N becomes strongly regulated, N_e and N would be negatively correlated as the blue dashed line shows. This trend is the paradox of changing N.

In summary, genetic drift effect can indeed decrease, increase or stay nearly constant as the population size increases, depending on the strength of N regulation near C_k. We further note that the V(K)/E(K) ratio of Eq. (7) can be independent (Fig. 3c), positively dependent (Fig. 3d), or negatively dependent (Fig. 3e) on N. Interestingly, the results of Fig. 3b are nearly identical when the V(K)/E(K) ratio changes (Supplementary Figs 1-3). These results show that the strength of genetic drift depends on the ecology that governs E(K), V(K) and N.

In the WFH framework, the paradox of changing N may not be an exception but a common rule in natural populations. Note that the WF approximation yielding Eq. (2) assumes nearly constant N. The assumption would imply very strong N regulation near C_k and, according to the DDH model, this is the condition for realizing the paradox of changing N. Coltman et al. (1999)’s observations of the reproduction in rams may be such a realization.

II. The paradox of genetic drift in sex chromosomes

Since the relative numbers of Y, X and each autosome (A) in the human population are 1:3:4, the WF model would predict the genetic diversity (θ_Y, θ_X and θ_A, respectively) to be proportional to 1:3:4. In a survey of human and primate genetic diversity, θ_Y is almost always less than expected as has been commonly reported (Hammer, et al. 2001; Wang, et al. 2014; Wilson Sayres, et al. 2014; Makova, et al. 2024). The low θY value has been used to suggest that human Y chromosomes are under either positive or purifying selection (Wang, et al. 2014; Wilson Sayres, et al. 2014; Makova, et al. 2024). Similarly, the reduced diversity X-linked genes was interpreted as a signature of positive selection (Pan, Liu, et al. 2022).

As pointed out above, under-estimation of genetic drift is a major cause of over-estimation of selective strength. Our goal is hence to see if genetic drift of different strength between sexes can account for the observed genetic diversities. Details will be presented in Wang et al. (2024). Below is a synopsis.

Let V_m and V_f be the V(K) for males and females respectively and α’ = V_m/V_f is our focus. (We use α’ for the male-to-female ratio in V(K) since α is commonly used for the ratio of mutation rate (Miyata, et al. 1987; Makova and Li 2002).) The three ratios, θ_Y /θ_X (denoted as R_YX), θ_Y /θ_A(R_YA) and θ_X /θ_A (R_XA), could be expressed as the functions of α’, which incorporate the relative mutation rates of autosomes, X and Y (Supplementary Information). These relative mutation rates can be obtained by interspecific comparisons as done by Makova and Li (2002).

We note that there are many measures for the within-species genetic diversity, θ = 4N_eμ. Under strict neutrality and demographic equilibrium, these measures should all converge. In the neutral equilibrium, the infinite site model dictates the frequency spectrum to be ξ_i = θ/i, where ξ_i is the number of sites with the variant occurring i times in n samples. Since every frequency bin is a measure of θ, different measures put different weights on the i-th bin (Fay and Wu 2000; Fu 2022). While Π, the mean pairwise differences between sequences, is most commonly used in the literature, we use several statistics to minimize the possible influences of selection and demography (Wang et al. 2024). In this synopsis, we used the Watterson estimator (Watterson 1975) as the measure for θ_Y, θ_X and θ_A by counting the number of segregating sites.

Using any of the three ratios (θ_Y /θ_X, θ_Y /θ_A and θ_X /θ_A), we can obtain α’; for example, R_YA = θ_Y /θ_A and, hence,

where y is the mutation rate of Y-linked sequences relative to autosomal sequences. With rearrangement,

Similar formulae of α’ can be obtained for R_YX and R_XA but the accuracy for estimating α’ is highest by R_YA whereas R_XA is the least accurate for this purpose (Supplementary Information).

Table 2 presents the α’ estimates in chimpanzees and bonobos. This is part of a general survey in mammals with a strong emphasis on primates and humans (Wang et al. 2024). It is almost always true that α’ > 5 in primates. Sources of data used are also given in Table 2. As shown for chimpanzees, α’ is often far larger than 5, above which the resolution is very low as can be seen in Eq. (8). Note that, when R_YA is under-estimated and approaching y/8, α’ would increase rapidly when the denominator is close to 0. That is why α’ often becomes infinity (Supplementary Fig. 5). The estimated α’ (MSE) in Table 2 alleviates the problem somewhat by using all three ratios to calculate the mean square error (MSE) (Supplementary Information).

Table 2.

Among primates surveyed, bonobo is the only exception with α’ < 1. While chimpanzees and bonobos are each’s closest relatives, their sexual behaviors are very divergent (de Waal 1995; De Waal and Lanting 2023). With unusually strong matriarchs, bonobo society seems to stand out among primates in sexual dominance. The α’ value is important in behavioral studies and, particularly, in primatology and hominoid research. If 0.1% of males sire 100 children and the rest have the same K distribution as females, V_m/V_f would be ∼ 10. Such outlier contribution is the equivalent of super-spreaders in viral evolution which may easily be missed in field studies. In this brief exposition, we highlight again that V(K) or, more specifically, V_m and V_f are key to genetic drift.

III. The paradox of genetic drift under selection

Genetic drift operates on neutral as well as non-neutral mutations. Let us assume a new mutation, M, with a frequency of 1/N is fixed in the population with the probability of P_f. Fisher (1930) first suggested that the fixation probability of an advantageous mutation should increase in growing populations, while decrease in shrinking populations. If there is no genetic drift (i.e., infinite population size in WF model), a beneficial mutation will always be fixed and P_f = 1. In the WF model, it is well known that P_f ∼ 2s for a new advantageous mutation, with fitness gain of s, while population size is large (Haldane 1932; Kimura 1962). This seems paradoxical that the determinant of genetic drift, 1/N_e, does not influence P_f (Otto and Whitlock 1997; Lanfear, et al. 2014).

In the Method section, we show that P_f under the Haldane model is given by

where M₀ and W₀ are the initial number of mutant allele (M allele) and wildtype allele (W allele), with N = M₀ + W₀. We note that u_M(t) and u_W(t), respectively, are the extinction probabilities of M allele and W allele by generation t with the initial number of alleles of 1. While obtaining the direct analytical solution of Eq. (9) may not be feasible, we could obtain its numerical solution due to its convergence as t increases (see Methods). The accuracy of the numerical solution from Eq. (9) is confirmed through simulation (Supplementary Fig. 4). Moreover, with the aid of the WF model (see Supplementary Information), we could obtain the approximation of Eq. (9) as follows.

We verify Eq. (10) by both the numerical solution from Eq. (9) and simulations based on the branching process (Supplementary Fig. 4). The fixation probabilities obtained by numerical solution vs. those inferred from Eq. (10) are shown in Fig. 4. The salient feature is that the fixation probability of 2s, as in the classical formula, would be a substantial over-estimate when V(K) is larger than E(K). Eq. (10) is sufficiently accurate as long as N≥50. When N is as small as 10, the theoretical result is biased. Indeed, at such a low N value, the population is prone to extinction. The DDH model presented above should rectify this deficiency.

Fig. 4.

The main message of Eq. (10) is that genetic drift under positive selection is influenced by s and V(K) but not by N. This independence from N seems intuitively obvious: when an advantageous mutation increases to say, 100 copies in a population (depending mainly on s), its fixation would be almost certain, regardless of N. This may be a most direct argument against equating genetic drift with N, or Ne (which is supposed to be a function of N’s).

Discussion

Genetic drift should broadly include all random forces affecting evolutionary trajectories. Many of these forces are excluded from the WF model, thus resulting in the many paradoxes addressed in the two companion papers, as well as the over-estimation of selection in the literature. We shall now dissect these paradoxes in detail in relation to the standard and modified WF models.

Compared with the WF models of genetic drift, the Haldane model is more intuitive as well as inclusive in depicting real populations. However, analytical solutions of the branching process can be more challenging to obtain. The WF model is often more tractable mathematically. For example, P_f ∼ 2s/V(K) is a good WF approximation (Kimura and Crow 1963) to the Haldane model when N is above 100 (Fig. 4). The integrated WFH model is hence both biologically realistic and analytically tractable (or at least computationally feasible).

Nevertheless, in many circumstances where the WF model is not operative, genetic drift defined by the Haldane model would demand new solutions. The DDH model is such an example, developed to resolve the paradox of changing N’s. The DDH model presents the conditions under which the WF model arrives in the approximate solution (i.e., Eq. (2)). The substitution of E(N) for N makes the approximation possible but effectively requires strong regulation of N. Interestingly, if N is strongly regulated, genetic drift is positively correlated with N, opposite of the WF model prediction (see Fig. 3). In short, the Haldane model can incorporate population size regulation, which is beyond the power of the WF model.

We shall further clarify the “paradox of changing N’s” on two fronts. First, the DDH model is on an ecological time scale whereby a population is either approaching or departing the carrying capacity (Fig. 2a). However, in many studies, the population size changes in the evolutionary time scale (Fig. 2b). For example, by the PSMC model (Li and Durbin 2011), both the European and Chinese populations experience a severe bottleneck 10-60 kyr ago. Presumably, various environmental forces may have reduced C_k drastically. Averaged over the long-time span, N should be at or very near C_k. Second, Gillespie has proposed the concept of “genetic draft” that may also reverse the effect of population size on the rate of substitution due to selection on linked variations (Gillespie 2000, 2001). Genetic draft of Gillespie is thus distinct from the strictly neutral process discussed in this study.

The standard WF model has been extended in several directions (overlapping generations, multiple alleles, ploidy, etc.). The modification most relevant to our studies here is the introduction of V(K) into the model, thus permitting V(K) ≠ E(K). While the modifications are mathematically valid, they are often biologically untenable. Kimura and Crow (1963) may be the first to offer a biological mechanism for V(K) ≠ E(K), effectively imposing the Haldane model on the WF model. Other models (Kimura and Crow 1963; Lynch, et al. 1995; Sjodin, et al. 2005; Der, et al. 2011; Cannings 2016) indeed model mathematically the imposition of the branching process on the population, followed by the WF sampling. The constructions of such models are biologically dubious but, more importantly, still unable to resolve the paradoxes. It would seem more logical to use the Haldane model in the first place by having two parameters, E(K) and V(K).

Even if we permit V(K) ≠ E(K) under the WF sampling, the models would face other difficulties. For example, a field biologist needs to delineate a Mendelian population and determine its size, N or N_e. In all WF models, one cannot know what the actual population being studied is. Is it the fly population in an orchard being sampled, in the geographical region, or in the entire species range? It is unsatisfactory when a population biologist cannot identify the population being studied. The Haldane model is an individual-output model (Chen, et al. 2017), which does not require the delineation of a Mendelian population.

We shall now review the paradoxes specifically in relation to the modified WF models, starting with the multi-copy gene systems such as viruses and rRNA genes covered in the companion study (Wang, et al. 2024). These systems evolve both within and between hosts. Given the small number of virions transmitted between hosts, drift is strong in both stages as shown by the Haldane model (Ruan, Luo, et al. 2021; Ruan, Wen, et al. 2021; Hou, et al. 2023). Therefore, it does not seem possible to have a single effective population size in the WF models to account for the genetic drift in two stages. The inability to deal with multi-copy gene systems may explain the difficulties in accounting for the SARS-CoV-2 evolution (Deng, et al. 2022; Pan, Liu, et al. 2022; Ruan, Wen, et al. 2022; Hou, et al. 2023; Ruan, et al. 2023).

We now discuss the first paradox of this study, which is about the regulation of N. In the general WF models, N is imposed from outside of the model, rather than self-generating within the model. When N is increasing exponentially as in bacterial or yeast cultures, there is almost no drift when N is very low and drift becomes intense as N grows to near the carrying capacity. As far as we know, no modifications of the WF model can account for this phenomenon that is opposite of its central tenet. In the general WF models, N is really the carrying capacity, not population size.

The second paradox of sex chromosomes is rooted in V(K) ≠ E(K). As E(K) is the same between sexes but V(K) is different, clearly V(K) = E(K) would not be feasible. The mathematical solution of defining separate N_e’s for males and females (Kimura and Crow 1963; Lynch, et al. 1995; Sjodin, et al. 2005; Der, et al. 2011; Cannings 2016) unfortunately obscures the interesting biology. As shown in Wang et al. (2024; MBE), the kurtosis of the distribution of K indicates the presence of super-breeder males. While the Haldane model can incorporate the kurtosis, the modified WF models are able to absorb only up to the variance term, i.e., the second moment of the distribution. The third paradox of genetic drift is manifested in the fixation probability of an advantageous mutation, 2s/V(K). As explained above, the fixation probability is determined by the probability of reaching a low threshold that is independent of N itself. Hence, the key parameter of drift in the WF model, N (or N_e), is missing. This paradox supports the assertion that genetic drift is fundamentally about V(K) with N being a scaling factor.

As the domain of evolutionary biology expands, many new systems do not fit into the WF models, resulting in the lack of a genetic drift component in their evolutionary trajectories. Multi-copy gene systems are obvious examples. Others include domestications of animals and plants that are processes of rapid evolution (Diamond 2002; Larson and Fuller 2014; Purugganan 2019; Chen, Yang, et al. 2022; Pan, Zhang, et al. 2022; Wang, et al. 2022). Due to the very large V(K) in domestication, drift must have played a large role. Somatic cell evolution is another example with “undefinable” genetic drift (Wu, et al. 2016; Chen, et al. 2017; Chen, et al. 2019; Ruan, et al. 2020; Chen, Wu, et al. 2022). The Haldane (or WFH) model, as an “individual output” model, can handle these general cases of genetic drift.

The Haldane model and the WF model are fundamentally different approaches to random forces of evolution. While the WF models encounter many biological contradictions, they have provided approximate mathematical solutions to more realistic scenarios. In systems such as in viral evolution (Ruan, Hou, et al. 2022; Hou, et al. 2023) or somatic cell evolution (Chen, Wu, et al. 2022; Zhai, et al. 2022) whereby the WF solution is absent, further development of the WFH model will be necessary.

Methods

Cell culture and image analysis

NIH3T3 cells, a fibroblast cell line that was isolated from a mouse NIH/Swiss embryo, were stably transfected with the fluorescent, ubiquitination-based cell cycle indicator (Fucci) (Sakaue-Sawano, et al. 2008) plasmid using Lipofectamine 3000 Transfection Reagent (Invitrogen) following the manufacturer’s specified instructions. The Fucci-labeled cells exhibited distinct fluorescent signals indicative of G1, S, and G2/M phases, represented by red, yellow, and green, respectively. NIH3T3-Fucci cell was derived from single-cell colony and cultured in DMEM supplemented with 10% Calf Bovine Serum and penicillin/streptomycin. All cells were maintained at 37◻ with 5% CO₂. Subsequently, the cells underwent extended time-lapse imaging using high-content fluorescence microscopy (PerkinElmer Operetta CLS) equipped with a 10x objective lens. Images were captured hourly over a 100-hour period, and the analysis was conducted using ImageJ (Fiji) (Schneider, et al. 2012) to count the number of cells (Supplementary Table 1).

Yeast strain construction

Strains were constructed on the genetic background of Saccharomyces cerevisiae strain BY4741. A GFP or BFP fluorescent protein expression cassette, under the control of the TDH3 promoter, was inserted into the pseudogene locus YLL017W. Transformations followed a published protocol (Gietz and Schiestl 2007). Transformants were plated on synthetic complete medium without uracil (SC-Ura), and from these, single colonies were selected. Confirmation of replacements was achieved through PCR, and cassette verification was performed using fluorescence microscopy. Subsequently, the constructed strains were cultivated non-selectively in YPD medium (1% Yeast extract, 2% Peptone, 2% D-glucose) at 30 ◻ on a rotary shaker.

Estimation of V(K) and E(K) in yeast cells

To discern division events, even under high concentration, we conducted co-cultures of the GFP-yeast and BFP-yeast at ratios of 1:1 and 1:25, with the initial cell concentration of 0.1% and 12.5%, respectively. Then the yeast cells were then continuously imaged under high-content fluorescence microscopy (PerkinElmer Operetta CLS) for 10 hours with 1-hour intervals to observe the individual offspring of GFP-yeast. Yeast cells with a distinct offspring number (K) within the initial 5 hours (with the high-density group limited to the first 4 hours) were documented (Supplementary Table 2). Subsequently, the mean and variance of K for each cell were calculated across various time intervals as follows.

According to the law of total variance,

where I_t is the offspring number at time t for an initial single cell (i.e., the total number of progeny cells for a single cell after t hours), as documented in Supplementary Table 2. And E_{t-1, t} represents the average of offspring number after single time interval (1 hour here) for each single cell from time t-1 to time t, while V_{t-1, t} is the corresponding variance. Utilizing the documented total cell count at time t (denoted by N_t), we could calculate the average of offspring number for a single yeast cell in a specific time interval, e.g., E_{t-1, t} = N_t / N_t-1 and E[I_t-1] = N_t /N₀. With some rearrangement,

Applying the aforementioned equations, we observed that both V_{t-1, t} (i.e., V(K) within a one-hour interval) and V(I_t) (i.e., the V(K) from 0 to time t) for yeast cells under high density (12.5%) are generally greater than the values under low density (Supplementary Table 3). This suggests that the variance of offspring number in high-density conditions could be larger than that in a low-density context. It’s noteworthy that the estimated value of V_{t-1, t} could be less than zero (Supplementary Table 3), implying a very small variance in the number of offspring. Furthermore, the occurrence of negative values for V_{t-1, t} is more frequent under low density than high density, indicating the higher variance of offspring number in high density as previously suggested.

To show the variance of offspring number overtime, we also calculated the change of offspring number within successive one-hour intervals, denoted as V(I_t - I_t-1). This value is intricately linked to the variance of offspring overtime. To facilitate the comparison with the total variance from 0h to 4h, we set V(ΔK) = 4V(I_t - I_t-1) to account for the four-hour time span.

Simulation of genetic drift in the Haldane model and the Wright-Fisher (WF) model

In both models, interactions between individuals are implicitly included through the dependency of the average number of offspring on population size, as defined by Eq. (5). This dependency leads to the logistic population growth, reflecting the density-dependent interactions. The initial population size is set to 4, with initial gene frequency (i.e., the frequency of mutant allele) of 0.5. And the carrying capacity is established at 100, 000, with r = 1 and z = 1. In WF model, the population size at next generation is N_t+1 = N_t ×E(K_t). Note if the calculated value of N_t+1 is not an integer, we round it to the nearest whole number. The number of mutant alleles follows a binomial distribution, allowing the simulation of gene frequency overtime in WF model. For Haldane model, the number of offspring number is assumed to follow a beta-binomial distribution. The mean of offspring number E(K_t) is obtained from Eq. (5). The variance V(K_t) is obtained from Eq. (6) with parameters a = 1/300 and b = 1.2. Given the mean and variance of the beta-binomial distribution, we simulated the number of mutant alleles and population size overtime and then traced the gene frequency overtime.

Fixation probability of a new mutation in Haldane model

Here we obtain the fixation probability of advantageous mutation in Haldane model governed by a branching process. The Haldane model considers a well-mixed population of N_t haploid parents with only two types of alleles (W for wildtype, M for mutant). There is no migration and new mutations in this model. Each individual is assumed to independently reproduce in discrete and non-overlapping generation, t = 0, 1, 2, …. Thus, the number of offspring per allele (also an individual in this haploid population) at any generation is represented by a set of identically and individually distributed random variables. In particular, the numbers of offspring of M and W alleles (denoted as K_M and K_W respectively) can be represented by following distributions.

For a particular distribution, we can obtain the average number of offspring as follows.

With the selection coefficient of s (0 for neutral mutation, positive for advantageous mutation, and negative for deleterious mutation) of M allele, average offspring number of M allele and W allele will follow the relationship.

Now, the evolution of the number of M allele (denoted as M_t) and the number of W allele (denoted as W_t) as time process is a branching process.

To compare Haldane model with WF model (the offspring number follows Poisson distribution with variance and mean equal to 1), we will set E(K_W) = 1, E(K_M) = 1 + s. And then let the offspring number follows a more general and realistic distribution (including negative binomial distribution (negBin), beta-binomial distribution (betaBin)), which can let the ratio of variance to mean range from 1 to a large number. (For simplicity, we let V(K_M)/E(K_M) = V(K_W)/E(K_W).) Based on the branching process, we obtained the fixation probability M alleles (see Supplementary Information for more details on the derivation).

Where

Note u_M(t = 1) = P(K_M = 0) = i₀, u_W(t = 1) = P(K_W = 0) = j₀. And both {u_M(t), t = 0, 1, 2, …} and {u_W(t), t = 0, 1, 2, …} are a bounded monotonic sequence. Thus, although we cannot obtain the direct analytical solution for P_f, we could easily obtain their numerical solution by iteration to the case when both u_M(t) and u_W(t) converge.

Acknowledgements

We thank Weiwei Zhai, Xionglei He for helpful comments. The work was supported by the National Natural Science Foundation of China (32150006, 32200493, 32293193, 81972691), the National Key Research and Development Projects of the Ministry of Science and Technology of China (2021YFC2301300, 2021YFC0863400).

Competing interest statement

Authors declare that they have no competing interests.