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

37226886
10.1093/genetics/iyad092
iyad092
Investigation
Population and Evolutionary Genetics
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
Featured
Self-contained Beta-with-Spikes approximation for inference under a Wright–Fisher model
Guerrero Montero Juan SUPA, School of Physics and Astronomy, University of Edinburgh, Edinburgh, EH9 3FD, UK

Blythe Richard A
Ralph P Editor
Corresponding author: SUPA, School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD, UK. Email: j.a.guererro-montero@sms.ed.ac.uk
Juan Guerrero Montero and Richard A Blythe contributed equally to this work.

Conflicts of interest The author(s) declare no conflicts of interest.

10 2023
25 5 2023
25 5 2023
225 2 iyad09210 3 2023
10 5 2023
19 8 2023
© The Author(s) 2023. Published by Oxford University Press on behalf of The Genetics Society of America.
2023
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

We construct a reliable estimation method for evolutionary parameters within the Wright–Fisher model, which describes changes in allele frequencies due to selection and genetic drift, from time-series data. Such data exist for biological populations, for example via artificial evolution experiments, and for the cultural evolution of behavior, such as linguistic corpora that document historical usage of different words with similar meanings. Our method of analysis builds on a Beta-with-Spikes approximation to the distribution of allele frequencies predicted by the Wright–Fisher model. We introduce a self-contained scheme for estimating parameters in the approximation, and demonstrate its robustness with synthetic data, especially in the strong-selection and near-extinction regimes where previous approaches fail. We further apply the method to allele frequency data for baker’s yeast (Saccharomyces cerevisiae), finding a significant signal of selection in cases where independent evidence supports such a conclusion. We further demonstrate the possibility of detecting time points at which evolutionary parameters change in the context of a historical spelling reform in the Spanish language.

Guerrero Montero and Blythe introduce a method for estimating effective population size and selection coefficients from time-series of allele frequency data; the method is based on the Beta-with-Spikes approximation to the Wright-Fisher process. The authors demonstrate its robustness with synthetic data and apply it to detecting selection in previously published allele frequency data for Saccharomyces cerevisiae. They extend the method to account for abrupt changes in model parameters, applying this extension to detecting language reforms in word-frequency data in Spanish.

Wright–Fisher model
time-series analysis
effective population size
selection strength
University of Edinburgh 10.13039/501100000848
==== Body
pmcStochastic models have a long history in population genetics as a tool to understand the fate of populations undergoing evolution and to draw inferences about the demographic and evolutionary forces that have shaped a present-day population. The starting point of many analyses is the Wright–Fisher model (Fisher 1930; Wright 1931; Crow and Kimura 1970) which characterizes the fluctuations that naturally occur as genetic material is replicated (genetic drift) and how these interact with mutation and selection. One key question is the extent to which variation in a population can be explained by a neutral model (Kimura 1983), that is one that appeals to genetic drift operating in the absence of selection. To this end, a variety of classical statistical tests were developed (see e.g. Kreitmann 2000, for a review) to detect departure from predictions of the neutral theory. Traditionally, these are based on quantities that can be ascertained from a sample of genetic material taken from a population at a single time, such as the number of nucleotide differences.

Latterly, interest has turned to the analysis of time-series data, in particular measurements of the frequency of alleles in a population at multiple points in time. Historically, such information has been obtained from evolution experiments conducted in the laboratory (Lenski et al. 1991). Meanwhile, advances in high-throughput sequencing technologies admit the collection of large data sets (Reuter et al. 2015), and there is particular interest in sampling microbial populations at multiple points in time (Dominguez-Bello et al. 2011). In addition to the expectation that additional time points will boost inferential power, the ability to analyze time series further opens the door to applications to other evolutionary paradigms, for example cultural evolution (Cavalli-Sforza and Feldman 1981; Boyd and Richerson 1988), in which analogs of nucleotide differences might not exist, but counterparts to allele frequencies do. The evolutionary process of language change is a case where one can identify processes that parallel mutation and selection (Croft 2000), which can furthermore be represented by a Wright–Fisher model (Baxter et al. 2006; Reali and Griffiths 2010; Blythe and Croft 2021) and where large-scale historical records of variation are available (Davies 2010; Michel et al. 2011). This broader context motivates further questions, such as whether one can detect changes in the sociocultural environment—for example, shifting attitudes towards specific behaviors—through changes in the intensity of evolutionary forces over time.

In this work, we develop a method to estimate the effective population size and selection strength in the Wright–Fisher model, to include the possibility that they each change over time, from a set of allele (or variant) frequencies obtained at different points in time. As in other works with a similar aim (see Tataru et al. 2017; Paris et al. 2019, for reviews), the basic strategy is to construct the likelihood of the observed parameter values and maximize with respect to the model parameters. This likelihood function involves the distribution of allele frequencies (DAF), conditioned on some initial value. Its construction has proved challenging: where exact solutions for the Wright–Fisher model exist, they are restricted to specific regimes and are cumbersome to work with, comprising for example an infinite series over special functions in the case of the diffusion approximation to the neutral Wright–Fisher model (Crow and Kimura 1970).

These difficulties have motivated a variety of schemes for approximating the DAF (Tataru et al. 2017; Paris et al. 2019). First, one can numerically integrate the diffusion equation corresponding to model of interest (Bollback et al. 2008). Although still computationally demanding, this has the advantage of being able to incorporate arbitrary evolutionary forces, such as frequency-dependent selection. It is possible to reduce these demands by expressing the DAF in terms of a sequence of orthogonal polynomials and truncating at a low order (Lukić and Hey 2012), although this requires some care in dealing with artifacts arising from the truncation. Another way to proceed is to approximate the DAF with an appropriate distribution function, such as a Gaussian (Lacerda and Seoighe 2014) or a Beta distribution (Hui and Burt 2015), and fix parameters by matching its moments to those obtained from the Wright–Fisher model. Such distribution functions are well behaved by construction, but may fail to adequately approximate the true DAF (Paris et al. 2019). A separate line of attack is offered by the coalescent process (Kingman 1982), which is dual to the Wright–Fisher model and particularly well adapted to reconstructing genealogies (Sirén et al. 2011). However, such methods lend themselves most naturally to neutral evolution, and become somewhat more complex in the presence of selection.

Here, we pursue the approach of approximating the DAF with a distribution function that is sufficiently rich to capture the properties of the underlying Wright–Fisher model, but has a small number of parameters that can be estimated efficiently. Specifically, we adopt the Beta-with-Spikes (BwS) distribution introduced by Tataru et al. (2015) as the functional form, and introduce a self-contained scheme that iteratively generates parameter values over multiple generations. The Beta distribution has been found to better describe changes in allele frequencies than a Gaussian distribution (Paris et al. 2019), primarily because the requirement that allele frequencies lie between 0 and 1 means that frequency differences are necessarily non-Gaussian as these boundary points are approached. Replacing the Gaussian with a Beta distribution rectifies this problem, but fails to account adequately for the accumulation of probability at the boundaries as individual realizations of the evolutionary dynamics cause an allele to reach fixation. It is precisely this shortcoming that augmenting the Beta distribution with delta functions (spikes) at the boundary points seeks to address (Tataru et al. 2015). By estimating fixation probabilities and moments of the Wright–Fisher DAF, the parameters in the BwS distribution can then be chosen to match: this is then found to improve the approximation to the exact DAF relative to the Gaussian or standard Beta approximations (Paris et al. 2019).

In Tataru et al. (2017), the moment-matching procedure is based on recursion relations for the mean and variance of the Wright–Fisher DAF that are not exact: they are based on a Taylor-series expansion of the fitness function around the mean allele frequency in the previous generation. This approximation breaks down when selection is large (Lacerda and Seoighe 2014), and errors accumulate over multiple generations, sometimes to the point of exiting the parameter regime for which the BwS distribution is well defined. The self-contained approach to estimating the BwS parameters that we introduce here avoids these problems. The basic idea is to start with a BwS distribution at the beginning of one generation, and to determine the parameter values that best approximate the distribution that results after one generation of evolution within the Wright–Fisher model. This leads to a self-contained recursion, in the sense that we map the BwS parameters directly from one generation to the next via averages with respect to the BwS distribution, rather than via the moments of the target Wright–Fisher distribution.

Methods

In this section, we set out the method for obtaining the self-contained approximation for parameters in the BwS distribution, and use synthetic data to compare the quality of the approximation with that obtained with previous moment-based approaches. We also outline how to use this approximation to obtain maximum-likelihood estimates for the effective population size and selection coefficient in the Wright–Fisher model, including the case where these change over time. In the following section, we will apply these methods to both genetic and cultural data sets.

Wright–Fisher model

The Wright–Fisher model describes the evolution of a population of randomly replicating individuals belonging to two or more types. In genetics, these types correspond to different alleles; in cultural evolution to variant forms of some socially learned behavior. We consider here the case of two variants, denoted by A and a, with xt representing the proportion of A individuals at generation t. The total population size, N, is assumed fixed, and each of the individuals in generation t+1 has probability g(xt) of being assigned type A (otherwise they are assigned type a). We refer to g(x) as the fitness function, and it can be interpreted as the mean number of offspring that each A leaves in the next generation. Given the above, the probability that the proportion of type A individuals in generation t+1 is equal to xt+1 is given by the binomial distribution

(1) PrWF(xt+1|xt)=(NNxt+1)g(xt)Nxt+1(1−g(xt))N(1−xt+1).

In population genetics, it is well understood that structured populations, where individuals are divided into classes by age, sex, location, or other characteristics, can be approximated by a Wright–Fisher model by setting N equal to an appropriate effective population size (Charlesworth 2009). The interpretation of N is less obvious in the cultural case; however concrete models of language use have found N to depend on factors like the size and structure of the speech community, the memory lifetime of individual speakers and the number of tokens produced in an utterance (Baxter et al. 2006; Reali and Griffiths 2010; Blythe and Croft 2021). In any case, N quantifies the effects of drift in the transmission: the lower N, the higher the uncertainty in the transmission to the next generation.

The form of g(x) can be chosen to incorporate a variety of evolutionary forces, such as mutation, migration, and selection of different types. Here, we focus on the case of frequency-independent selection, which is traditionally modeled by assigning a weight (relative reproductive success) of 1+s0 to A and 1 to a. This yields a fitness function (1+s0)x(1+s0x). Here, we have found it helpful to adopt a different parametrization, in which the weights are es/2 and e−s/2, respectively. This formulation has two particularly appealing features.

First, changing the sign of s is equivalent to exchanging A and a: that is, there is a symmetry between positive and negative values of s. This further implies that while we need s0≥−1 for the model to be well defined, s can take any positive or negative value. The second feature is that the fitness function

(2) g(x;s)=xx+(1−x)e−s

satisfies the functional relation g(g(x;s1),s2)=g(x;s1+s2). Thus, in the absence of fluctuations, the evolution over k generations with selection coefficients s1,s2,…,sk is the same as a single generation of evolution with selection coefficient s1+s2+⋯+sk. Below, we will explain how this property can be exploited to speed up the inference of evolutionary parameters by aggregating multiple generations into one. Meanwhile, we note that the two different specifications of the fitness function g(x;s) can be mapped onto each other via the relation s=ln(1+s0). For small values of s (or s0), we further have that s≈s0.

BwS approximation

Allele frequency data sampled at different time points may not be one generation apart. In this case, it is necessary to sum (1) over multiple intermediate generations to obtain the appropriate DAF. This is not viable in practice for large population size N and large numbers of intermediate generations k, as the memory requirements for this procedure scale as O(N2), and its computational complexity as O(kN3) (Paris et al. 2019). These considerations motivate approximating the DAF with some distribution that has a small number of parameters, and using (1) to determine how those parameters change over multiple generations.

Tataru et al. (2015) introduced the BwS distribution for this purpose. Starting at generation t with a fixed frequency xt, the distribution after k generations is assumed to be well described by the form

(3) PrBwS(xt+k|xt,N,s)=P0,kδ(xt+k)+P1,kδ(1−xt+k)+(1−P1,k−P0,k)xt+kαk−1(1−xt+k)βk−1B(αk,βk),

which has four parameters αk, βk, P0,k, and P1,k. The central part of this distribution is a Beta distribution, whose shape is controlled by αk and βk and whose normalization is given by the Beta function B(αk,βk). This distribution has been employed by itself as an approximation to the Wright–Fisher DAF (Hui and Burt 2015), and is well adapted for that purpose by virtue of being defined on the interval 0≤x≤1 (unlike the Gaussian distribution which assigns a finite probability to unattainable frequencies). Furthermore, a variety of shapes can be accessed by tuning αk and βk, including a uniform distribution, distributions that are strongly peaked around the mean, and those that have an integrable divergence at the boundaries.

This flexibility is however insufficient to capture the accumulation of probability at the boundary points which occurs when an allele goes to fixation. This possibility is incorporated into Equation (3) through the two Dirac delta function contributions at the extremes of the interval. The quantities P0,k and P1,k then correspond to the probability that alleles a and A have fixed after k generations, respectively. The Beta distribution contribution then describes the DAF conditioned on fixation not yet having occurred.

The crucial step in applying the BwS approximation to data is to estimate the parameters αk, βk, P0,k, and P1,k. The general approach is to estimate moments of the DAF and fixation probabilities in the Wright–Fisher model, and choose the parameters in the BwS distribution so that they match up. To this end, we note first that by defining ak and bk as

(4) αk=(Ek*(1−Ek*)Vk*−1)Ek*

(5) βk=(Ek*(1−Ek*)Vk*−1)(1−Ek*),

the mean and variance of the Beta part of the BwS approximation will match the mean Ek* and variance Vk* of the Wright–Fisher DAF after k generations, conditioned on fixation not having occurred. These can be obtained from the mean Ek and variance Vk of the full distribution, as well as the fixation probabilities P0,k and P1,k, via

(6) Ek*=Ek−P1,k1−P0,k−P1,k

(7) Vk*=Vk+Ek2−P1,k1−P0,k−P1,k−(Ek*)2.

It remains to estimate Ek, Vk, and the two fixation probabilities.

Truncated Taylor-series estimation scheme

Previous approaches to estimation (Tataru et al. 2017; Paris et al. 2019) have been based around a Taylor-series expansion of the fitness function. Since we will use results from this method in a comparison below, we briefly summarize the procedure.

By recursively applying the law of total probability to the Wright–Fisher transition probability (1), the DAF after k+1 generations with starting frequency xt takes the form:

(8) PrWF(xt+k+1|xt)=∑xt+kPrWF(xt+k+1|xt+k)PrWF(xt+k|xt).

As set out in Supplementary Section S1, one can use this expression to derive exact recursions for the mean, variance, loss probability and fixation probability. These are

(9) Ek+1=E[g(xk)]

(10) Vk+1=(1−1N)Var[g(xk)]+1NE[g(xk)](1−E[g(xk)])

(11) P0,k+1=E[(1−g(xk))N

(12) P1,k+1=E[g(xk)N],

where the mean and variance on the right-hand side are with respect to the Wright–Fisher DAF after generation k. These recursions are closed in the first two moments only for linear fitness functions g(x), and therefore cannot be computed exactly for nonzero selection coefficients. Closure can be obtained by Taylor expanding g(xk) about Ek to second order and dropping any higher order moments that appear (see Supplementary Section S1.1). This yields

(13) Ek+1≈g(Ek)+12Vkg″(Ek)

(14) Vk+1≈1NEk+1(1−Ek+1)+(1−1N)Vkg′(Ek)2.

This, however, is not exact, with the error at each step increasing with the selection coefficient s. Consequently, this recursion is expected only to be valid for small s. The fixation probabilities can be estimated by considering the probability at each boundary within the Beta part of the distribution (Tataru et al. 2017)

(15) P0,k+1≈P0,k+(1−P0,k−P1,k)B(αk,βk+N)B(αk,βk)

(16) P1,k+1≈P1,k+(1−P0,k−P1,k)B(αk+N,βk)B(αk,βk).

Despite the approximations made, these recursion relations benefit from being simple and quick to apply.

Self-contained estimation scheme

Here, we take a different approach to estimating the BwS parameters, which is motivated by the expectation that it will keep the accumulation of error under control. The basic idea is to take the BwS distribution obtained after k generations, and generate the intermediate distribution

(17) Print(xt+k+1|xt)=∫PrWF(xt+k+1|xt+k)PrBwS(xt+k|xt)dxt+k

by applying just one step of the Wright–Fisher process (1). We then examine the moments and fixation probabilities of this intermediate distribution, and use the values obtained to set the BwS parameters for generation k+1. We view this as a self-contained estimate, as it maps directly from one set of BwS parameters to the next.

Figure 1 compares the Wright–Fisher transition probability with the intermediate and BwS distributions for the case of high selection (s=0.2) and high drift (N=50). Even in this challenging regime, we find the intermediate and BwS distributions remain similar to the exact Wright–Fisher transition probability after k=8 generations. The continuous BwS distribution is generated by mapping its mean, variance, loss and fixation probabilities to those of the intermediate distribution. These are all derived in Supplementary Section S1.2.

Fig. 1. Top panel: Comparison of the intermediate and Wright–Fisher distributions. Bottom panel: Comparison of the BwS and Wright–Fisher distributions. Distributions generated with N=50, s=0.2, x0=0.5 after k=8 generations.

For the mean and variance, we find

(18) Ek+1=P1,k+(1−P0,k−P1,k)Ek[g(x)]

(19) Vk+1=(1−1N)[P1,k+(1−P0,k−P1,k)Ek[g(x)2]]+1NEk+1−Ek+12,

where Ek[⋅] represents the expectation value under the Beta distribution with parameters αk and βk:

(20) Ek[f(x)]=∫01f(x)xαk−1(1−x)βk−1dxB(αk,βk).

Meanwhile, the loss and fixation probabilities are obtained as

(21) P0,k+1=P0,k+(1−P0,k−P1,k)Ek[(1−g(x))N]

(22) P1,k+1=P1,k+(1−P0,k−P1,k)Ek[g(x)N].

This method does involve more computation than the Taylor-series approach, in that four integrals of the type (19) have to be performed in each iteration, which can increase the computation time by a up to a factor of order 10, compared to an iteration of the Taylor-series approach. However, we have found this effort is manageable in practice, and if necessary can be reduced by appealing to scaling properties identified below. In any case, some care is needed to evaluate these numerically when αk<1 or βk<1 (or both), as well as when the integrand in Equation (19) is sharply peaked around its mode. We set out the numerical integration algorithms in Supplementary Section S2, and furthermore provide a link to the code used to obtain our results in Data availability.

Comparison of estimation schemes

To assess the relative quality of the two estimation schemes, we construct a baseline BwS distribution which is obtained by computing the moments and fixation probabilities within the Wright–Fisher model using numerical methods that are exact to machine precision. Then, we can measure the distance between this baseline and each of the BwS distributions obtained through the estimation schemes set out above. For this purpose, we use the Wasserstein distance, bearing in mind that a smaller distance indicates a better approximation to the baseline. (Note that comparing to the exact Wright–Fisher distribution is difficult, because this is a discrete distribution while the BwS distribution is defined over continuous frequencies.) The results are shown in Fig. 2.

Fig. 2. Wasserstein distance between the BwS distribution with numerically exact moments and with approximated moments for two approximation schemes and three values of the selection strength s, as a function of the initial frequency x0 and the generation k. Left: results for the self-contained approximation. Right: results for the approximation based on the truncated Taylor expansion. Top figures: weak selection (s=0.01). Middle figures: intermediate selection (s=0.1). Low figures: strong selection (s=0.6). Sudden increase of the Wasserstein distance to a value of 1.0 in the approximation based on the truncated Taylor expansion for intermediate and strong selection is due to accumulation of error that leads to an undefined DAF. Results for N=100.

In all cases the population size N=100, and we compare performance as a function of the number of generations k, the initial frequency x0 and for three different values of the selection coefficient. The estimation based on the truncated Taylor expansion quickly accumulates error in its estimation of the moments when the selection coefficient is sufficiently large, which may lead to negative values of α and β and an undefined DAF. This happens at generation k=23 for the DAF with intermediate selection and as early as k=7 for the DAF with strong selection. This is reflected in the figure as points where the Wasserstein distance reaches a maximal value of 1. The self-contained estimation never leads to an undefined DAF, provided sufficiently accurate integration schemes are used, as α and β as defined in Equation (4) are always nonnegative given that P0, P1, E, and V are all obtained from the same well-defined distribution.

Figure 3 provides a closer look at the accuracy in the estimation of the mean E, variance V, loss probability P0 and fixation probability P1 for both the self-contained and truncated Taylor estimation schemes, with N=100 and s=0.1. For each one of the four parameters, the absolute difference between their estimated and exact values is plotted as a function of the initial probability x0 and generation k. Again, the plots demonstrate the robustness of the self-contained estimation, keeping an absolute error below 0.1 for all data points, whereas the truncated Taylor expansion surpasses this value for all parameters before generation k=30.

Fig. 3. Absolute value of the difference between exact and approximated values of the mean, variance, loss probability and fixation probability, as a function of the initial frequency x0 and the generation k. Left subpanels: self-contained approximation. Right subpanels: truncated Taylor expansion approximation. Results for N=100, s=0.1.

Maximum-likelihood inference from time-series data

We now discuss how to infer the effective population size, N, and selection coefficient, s, through likelihood maximization. The situation we have in mind is when one has samples drawn from the population at a sequence of times t1,t2,…,tn. One of the problems we will have to contend with is how many generations of the Wright–Fisher model a particular time interval ti+1−ti corresponds to, which may or may not be known a priori. We will set out below a procedure that allows us to deal with this uncertainty, and relate a specific time interval Δt to a number of generations k.

Another source of uncertainty is in the allele frequencies themselves, as these will be subject to a sampling error that decreases in magnitude as the sample size is increased. In the first instance, we assume that samples are sufficiently large that such error can be neglected. Later, under Regulated and unregulated language change, we set out a scheme that can mitigate against changes in the sample size over time, and will outline extensions to our approach that account for finite sample sizes more rigorously in the Discussion.

Given the above, we can identify the likelihood of the observed data within the BwS approximation as

(23) L(X|N,s)=∏iPrBwS(xti+1|xti,N,s),

where X is the sequence of frequency measurements X=(xt1,xt2,…,xtn). Typically, the likelihood function has a single maximum, meaning that the optimal values of N and s can be located straightforwardly with a standard optimization algorithm (Press et al. 2007).

It is often desirable to distinguish neutral from nonneutral evolution, that is, whether the maximum-likelihood value of s is significantly different from zero. This can be achieved by examining the likelihood ratio

(24) λ=2ln(L(X|N,s)L(X|N0,0))

in which N and s are the optimal parameters when s is unconstrained, and N0 is the optimal effective population size when the constraint s=0 is imposed. The set of all trajectories generated by genetic drift with an effective population size of N0 provide a null distribution for λ. If, for any empirical trajectory, λ is found to lie in the tail of that distribution, we can regard it as significantly different from drift. For long time series, Wilk’s theorem (Wilks 1938; Casella and Berger 2001) may hold, and λ can be assumed to follow a χ2 distribution. In the general case, this distribution can be constructed by generating a set of artificial pure drift time series following a Wright–Fisher process with N=N0 and s=0. The P-value can then be computed as the fraction of these time series whose likelihood ratio is higher than that of the data of interest, as outlined by Feder et al. (2014).

Figure 4 compares the performance of both the self-contained and the truncated Taylor approximations of the moments in the estimation of the selection parameter s from artificially generated time series, using maximum-likelihood inference. While both approximations perform similarly at low ks (the product of the true selection parameter and the number of generations between data points) the numerical stability of the self-contained scheme keeps the error in the estimation decidedly lower starting at ks=0.8.

Fig. 4. Relative error in the estimation of the selection parameter s in artificially generated time series, using a BwS-based maximum-likelihood inference with both the self-contained and the tructated Taylor approximations of the moments, as a function of ks, the product of the true selection parameter s and the number of generations k between data points.

We now turn to the issue of matching a real time interval, Δt, to a number of generations, k, in the Wright–Fisher model. In situations where this is not known, we can appeal to scaling properties of N and s with k to obtain values whose scales are set primarily by Δt and only weakly by k. As previously mentioned, the fitness function defined by Equation (2) satisfies in the deterministic limit (N→∞) an exact scaling relation whereby k generations, each lasting Δtk, with selection coefficient sk are equivalent to a single generation, lasting Δt, with coefficient s1=ksk. Meanwhile, in the diffusion limit (N≫1) and pure drift (s=0), k generations with drift coefficient Nk are equivalent to a single generation with drift coefficient N1=Nk/k (e.g. Crow and Kimura 1970). In the general case, we propose that for any number of generations k and k~ between two data points separated by a time interval Δt, we have the scaling behaviors

(25) ksk=k~sk~,Nk/k=Nk~/k~.

These relations allow us to divide a time interval into k steps for the purpose of performing the analysis, and quote effective population sizes and selection coefficients appropriately for a standardized time interval determined by k~. For example, data could be presented at intervals of ten years, divided into k=2 steps of five years for the purposes of analysis, and quoted for a standardized interval of one year (k~=10) to facilitate comparison of analyses performed for different time series.

The success of this approach depends on Equation (24) holding with reasonable accuracy for general N and s, beyond the special limits described above. Figure 5 confirms this for the case of synthetic data generated by iterating Equation (1) 10 times between updates. In this case, the true number of generations between sample points is k~=10, but we choose to analyze with a different number, k. The error in the maximum likelihood estimates of Nk/k and ksk, relative to the true values, is shown as a function of k in Fig. 5. We find that the error to be modest (around 10% or less) even when the parameters are far from the values that make the scaling behavior exact (low N for the scaling behavior of s, high selection strength Ns for the scaling behavior of N). What this means in practice is that one can reduce the number of iterations of the Self-contained estimation scheme to a small number k, by dividing the time between data points into k generations, while retaining reasonable parameter estimates for some chosen standardized generation time.

Fig. 5. Relative error in the scaling behavior of N and s as time series with 10 generations between data points are reanalyzed as having k<10 generations between data points. Statistics generated as the average over 2,000 artificially generated time series with s=0.05 and N=100, N=1,000, N=10,000.

Time-dependent evolutionary parameters

Many previous works on analyzing evolutionary time series have assumed constant effective population size and selection coefficient. We may however be interested in situations where these parameters change. For example, in cultural evolution, the selection coefficient may change due to one form of behavior gaining social prestige or being stigmatized. Similarly, in genetic evolution, the appearance of a new predator or pathogen could affect an organism’s fitness.

It is straightforward to extend the maximum-likelihood approach outlined above to the case where different parameters apply over different time intervals. In particular, a single abrupt change can be modeled as two sets of parameters (N,s) that apply before and after a transition time T. For t<T, the parameters take values (N1,s1), and for t>T, they take values (N2,s2). The optimal parameters N1, s1, N2, s2 and T for a frequency time series X can be found by maximizing the likelihood function L(X|N1,s1,N2,s2,T) with respect to all five parameters.

Again, one can use likelihood ratios to determine whether the time division provides a significantly better explanation of the data. Specifically, we consider the ratio

(26) λ=2ln(L(X|N1,s1,N2,s2,T)L(X|N,s))

which compares the optimal likelihood of a model where both N and s change at a time T with one where N and s are fixed for the entire trajectory. Since these models are not nested (on account of the simpler model being found by setting the division point T to its maximum or minimum possible value), λ cannot be assumed to be χ2-distributed. Instead, a P-value must be obtained by constructing the empirical distribution of λ from trajectories in which there is no time division. This approach can be extended to multiple abrupt changes by further subdividing the trajectories.

Sensitivity to changes in selection strength

To test the ability of this algorithm to detect changes in selection strength, we generate artificial time series by iteratively sampling allele frequencies from Equation (1) for T generations. At generation T/2, the selection strength goes from s=0 to s=Δs. For each set of values of 500≤N≤10,000, 6≤T≤50 and 0.001≤Δs≤0.3, we generate 2,000 time series, compute their likelihood ratios using Equation (25) and find the associated P-values. If a P-value is under the standard significance threshold, P<0.05, the change in selection is considered detected.

Figure 6a shows the fraction of time series for which changes in s are detected as a function of the change in selection strength Δs at fixed N=2,000 and T=20. As expected, this fraction grows monotonically with Δs. Empirically, we find that it is well fit by the logistic function

(27) f(Δs)=11+exp(−a−bΔs)

which allows us to identify a characteristic Δs through the value for which the fraction of detected changes equals one half. That is above this value, we are more likely to detect the change than not.

Fig. 6. a) Fraction of significant selection as a function of Δs for N=2,000, T=20, together with fitted logistic function and estimated characteristic value of Δs. The fitted logistic function (Equation (26)) has parameters a=−3.1±0.3, b=51±5, r2=0.992. b) Characteristic Δs as a function of N for fixed T=20. Power law has proportionality constant c=0.284±0.015, exponent d=−0.52±0.02, and r2=0.994. c) Characteristic Δs as a function of T for fixed N=2,000. Power law has parameters c=2.33±0.12, d=−0.481±0.008, r2=0.999. Dots in subfigure (a) represent empirical points obtained as the average of the detection of change in 2,000 artificially generated time series. The crossed dot in (a), as well as all dots in (b) and (c) represent characteristic values of Δs, obtained from interpolated logistic functions.

Figure 6b and c shows how this characteristic value varies with effective population size N and the number of generations T, respectively. We find that the larger the effective population size or the longer the time series, the smaller the change that can be detected. In both cases this is expected: fluctuations diminish as N increases, thereby increasing the signal-to-noise ratio. Similarly, longer time series provide more information and allow stronger inferences to be drawn. Figure 6b indicates that relatively small changes are detectable within trajectories of a modest length (e.g. around T=10 generations).

Applications to empirical data

Having validated the methods of the previous section with synthetic data, we now apply them to empirical data from previously published studies. In doing so, we demonstrate their applicability to both genetic and cultural evolution.

Beneficial mutations in yeast populations

We first analyze data from an experiment carried out by Lang et al. (2011), in which 592 populations of baker’s yeast (Saccharomyces cerevisiae) were evolved over 1,000 generations in a rich environment. Several of these populations were deep sequenced, revealing the presence of many adaptive mutations (Lang et al. 2013). Feder et al. (2014) identified three mutations, affecting genes STE11, IRA1, and IRA2, which were likely to be beneficial, this based on their appearance and spreading in several populations (Lang et al. 2013). The mutant allele frequency trajectories are shown in Fig. 7. To test for selection, Feder et al. (2014) used two methods, each based on a Gaussian approximation to the DAF. The first of these uses the distribution of the likelihood ratio, as described in the section Maximum-likelihood inference from time-series data. The second is a simpler test, based around the idea that rescaled differences between allele frequencies at subsequent time points should all be drawn from a standard Normal distribution. Despite the independent evidence that all three mutations were being selected for, and in spite of applying the methods to arbitrarily selected subsets of the time series, only one of the six analyses (the likelihood ratio test applied to the last four data points in the IRA1 time series) showed a significant P-value for selection (Feder et al. 2014, Supplementary Table S2).

Fig. 7. Left: trajectory of allele frequencies of mutation D579Y in gene STE11 in population RMB2-F01. Centre: trajectory of allele frequencies of mutation Y822* in gene IRA1 in population RMS1-D12. Right: trajectory of allele frequencies of mutation A2698T in gene IRA2 in population BYS2-D06.

Here, we repeat the analysis using the likelihood ratio test combined with the BwS approximation and both the self-contained and the truncated Taylor estimation schemes. One advantage of methods based on the BwS approximation is that there is no need to truncate the time series to exclude fixation events, a step that was required in the analysis of Feder et al. (2014). There is also no incentive to exclude problematic points close to boundary values because the BwS approximation stays robust. We use a time unit of k~=5 generations (the greatest common divisor of the number of generations between data points for all three time series), and rescale s and N as described under Maximum-likelihood inference from time-series data.

Our results using the self-contained approximation of the moments on the untruncated time series are shown in Table 1, indicating that we find a significant P-value in two cases, consistent with the independent evidence of selection. The analysis of gene STE11 produces a low, albeit not below the threshold of significance, P-value of 0.084. This additional sensitivity to selection compared to previous analyses most likely derives from improved handling of the DAF when allele frequencies approach the boundary values. Note how models with selection have a greater optimal population size N than those that rely solely on drift to explain the behavior of the data, as greater effective population size corresponds to lower drift and more deterministic trajectories.

Table 1. Results of the analysis of the allele-frequency trajectories of genes STE11, IRA1, and IRA2 using the self-contained scheme for the estimation of the moments.

Gene	N 0	N	s	P-value	
STE11	1,350	1,860	0.010	0.084	
IRA1	1,240	1,970	0.0091	0.006	
IRA2	987	1,980	0.019	0.0	
N0 , optimal population parameter under the null-model of pure drift. N, optimal population parameter under model with selection. s, optimal selection strength.

Our results using the truncated Taylor-series scheme, shown in Table 2, show qualitatively similar results, but higher P-values for all three time series. This is likely due to the numerical instability of the truncated Taylor method, which produces diverging likelihood ratios in the empirical statistics used to compute the P-value, artificially increasing its value. These observations, combined with the results for synthetic data, suggest that the self-contained estimation scheme allows selection to be more reliably detected. It is also of note that the computation time for these frequency time series under the self-contained scheme was only 7 times slower than its truncated Taylor counterpart, perfectly manageable for this type of analysis.

Table 2. Results of the analysis of the allele-frequency trajectories of genes STE11, IRA1, and IRA2 using the truncated Taylor-series scheme for the estimation of the moments.

Gene	N 0	N	s	P-value	
STE11	627	1,590	0.0097	0.19	
IRA1	526	1,440	0.015	0.014	
IRA2	377	794	0.018	0.022	
N0 , optimal population parameter under the null-model of pure drift. N, optimal population parameter under model with selection. s, optimal selection strength.

Regulated and unregulated language change

We now turn to cultural evolution and examine the dynamics of historical changes occurring in the 2019 update to the Spanish Google Books corpus (Michel et al. 2011). We look at both a regulated and unregulated change that occurred between the 19th and 20th centuries (Amato et al. 2018). The regulated change was a spelling reform introduced by the Real Academia Española, the central regulatory institution of the Spanish language, in their Gramática de la lengua castellana (Real Academia Española 1911), and entailed the accentuated á (meaning to), ó (meaning or), é (alternative form of y, meaning and), and ú (alternative form of ó) being replaced by their unaccented forms a, o, e and u. The unregulated change occurred in the absence of such an intervention, and involved the competition dynamics between two completely equivalent forms of the past subjunctive tense, with verbal affixes -ra- and -se-. Thus, the third person singular of the past subjunctive of the verb evolucionar (to evolve) could be either evolucionara or evolucionase. Both forms of the past subjunctive are considered completely equivalent in all contexts. In spite of this, in the last 150 years, there has been a steady transition in the corpus, from a clear preference of the -se- form to a clear preference of the -ra- form (see lower panel of Fig. 8). As of yet, there is no agreed upon explanation of this phenomenon, from either corpus-based or sociolinguistic perspectives (Kempas 2011; Guzmán Naranjo 2017)

Fig. 8. Top: regulated change. Trajectory of frequency of usage in the Spanish Google Books corpus of old, accentuated spellings of words a, e, o and u, with detected years of change in the selection parameter marked with vertical lines. The first detected year in 1910 is only a year before the true year of introduction the orthographic reform by the RAE which declared the old spellings nonstandard. Bottom: unregulated change. Trajectory of frequency of usage in the Spanish Google Books corpus of the -ra- form of the past subjunctive, as opposed to the -se- form, for which the time-divided model is not significant.

The processes of regulated and unregulated change have previously been modeled by the cultural analogs of mutation and migration in a large population (Amato et al. 2018). In the notation of the present work, this corresponds to a linear fitness function g(x)=ax+b and a large fixed value of the effective population size N, with residuals modeled by a standard Normal distribution rather than the Wright–Fisher model. Given that we are dealing with a competition between pre-existing variants, we consider it more natural to view the evolution as being driven by selection, albeit where the selection coefficient may change over time, for example, due to the imposition of the reform, or because social preferences and norms can change over time. To this end, we turn to the method described under Time-dependent evolutionary parameters to detect changes in the selection coefficient with fluctuations in variant frequencies accounted for through the cultural analog of genetic drift (whose amplitude may also change over time).

An issue that we have to contend with when dealing with historical language data is that the sample sizes change over time. Specifically, the general trend is for data to become scarcer as earlier time periods are examined. The increased sampling fluctuations at early times could then be misattributed to drift, that is, the intrinstic fluctuations in the cultural transmission process, rather than the sampling of linguistic data from the population. One way to address this is to create subsamples of the larger data sets in the time series, with the subsample size chosen in such a way that the contribution from sampling is of equal magnitude across the time series. Then, any detected change in the effective population size must be due to changes in the intrinsic fluctuations, rather than sampling. In practice, we achieve this by generating for each time point a binomial random variable with a success probability equal to that of the original sample, but a sample size m given by

(28) m=m01−m0n,

where n is the original sample size and m0 is the smallest original sample size across the entire time series. This formula is derived in Supplementary Section S3.

After constructing the resampled time series, we apply the method of Time-dependent evolutionary parameters to estimate parameter and P-values were found for models with and without a single time division. When this time-divided model has a P-value below 0.05, we repeat the process for each of the subseries, accepting subsequent time divisions whenever P<0.05. The result of this analysis is shown in Fig. 8, with estimated parameter and P-values given in Tables 3 and 4.

Table 3. Results for the analysis of unregulated change in the affixes of past subjunctive verbal forms in Spanish between the years 1850 and 2000, using time-divided models.

T	N 1	N 2	s 1	s 2	P	
1961	46.9	124	0.016	0.026	0.15	
The time division is found not to be significant.

Table 4. Results for the analysis of regulated change in the accentuation of single-letter words in Spanish between the years 1850 and 2000, using time-divided models.

T 1	T 2	N 1	N 2	N 3	s 1	s 2	s 3	p 1	p 2	
1910	1920	63	72	197	− 0.002	− 0.45	− 0.013	0.0	0.0	
Two time divisions are found to be significant, in 1910 and 1920, delimiting the transition process between the old and new spelling rules introduced by the RAE in 1911.

A single time division model is not found to be significant for the process of an unregulated change, thus suggests that it has not been driven by any abrupt change in the social perception of either the -ra- or -se- forms of the past subjunctive. By contrast, we find that a first division at T=1910 followed by a second division at T=1920 are both significant in the case of a regulated change. The first time division delimits an early period where the accented spellings (á, ó, é, ú) were widely used, and one where they rapidly fell out of use. The division point falls at the start of the decline, and is in fact only one year before the introduction of the reform (Real Academia Española 1911). Thus, it seems likely that the reform caused individual language users to change their attitude towards the accented forms. We also note that the estimated effective population size does not change across this first time division. The second time division falls at the end of the period of decline, and we note from Table 4 that the selection coefficient is estimated to be much smaller than during the transition period. It is perhaps the case that modern Spanish speakers encounter the accented forms sufficiently rarely that they do not hold any particular disposition towards it. The significance or otherwise of the change in effective population size is somewhat less clear, and we do not speculate further. We conclude this section by noting that if one does not account for the possibility of the selection coefficient changing in time, one does not find a significant effect of selection (according to the likelihood ratio test discussed under Maximum-likelihood inference from time-series data).

Discussion

In this work, we have introduced a method for obtaining reliable maximum-likelihood estimates of parameters within the Wright–Fisher model from time-series data for allele frequencies. Our approach is underpinned by the BwS distribution (Tataru et al. 2015, 2017) which, despite not exactly matching the DAF within the Wright–Fisher or related models, captures its essential features. These are the possibility of extinction or fixation of an allele, accounted for by the spikes, and that unfixed alleles are governed by a continuous distribution that is well characterized by its mean and variance.

The challenge in utilizing the BwS approximation is accurately determining appropriate parameter values. Earlier works (Tataru et al. 2017; Paris et al. 2019) used Taylor-series expansions to estimate how parameter values should change from one generation to the next. These have the benefit of being simple to evaluate, but the truncation of the Taylor series results in the approximation being unreliable when the selection coefficient is large. In particular, it can generate parameter values that cause the BwS distribution to be ill-defined.

Here, we have turned to a self-contained approximation, where the BwS distribution is used as the initial condition for one step of Wright–Fisher evolution, and a fresh BwS distribution is fit to the intermediate distribution that results. In this approach, the approximating distribution remains well defined, and provides an adequate approximation to the Wright–Fisher model even when selection is strong. We have demonstrated the reliability of the method for the Wright–Fisher model with frequency-independent selection by comparing distributions directly, and by determining the error on the maximum-likelihood estimate of the selection coefficient for artificial time series where the true value is known.

The method is however more computationally intensive than the Taylor-series approach, since it is necessary to compute four integrals over a BwS distribution at each generation. Nevertheless, we find that the method can be applied to data for both genetic and cultural evolution without undue computational effort. In particular, we can reduce the number of calculations that need to be performed by aggregating multiple generations into a single effective generation, and appeal to the scaling properties of the effective population size and the selection strength when this is done, as described under Maximum-likelihood inference from time-series data.

In the section Beneficial mutations in yeast populations, we found that we were able to obtain a significant signal of selection for two out of three genes for which there is independent evidence of selection (Lang et al. 2013; Feder et al. 2014), whereas other methods either break down or do not yield a uniformly significant result, even when time series are truncated. Although we cannot be certain that the gene frequencies were driven by selection in all cases, our results suggest that the inability to reject the null hypothesis of drift in Feder et al. (2014) may lie in the sensitivity of the test that was applied.

As noted in the introduction, cultural evolutionary processes, such as language change, can also be couched in evolutionary terms (Cavalli-Sforza and Feldman 1981; Boyd and Richerson 1988; Croft 2000) and furthermore represented mathematically by the Wright–Fisher model (Baxter et al. 2006; Reali and Griffiths 2010; Blythe and Croft 2021). Until recently, the analysis of historical corpus data within this framework has been hampered by the limited availability of methods that can be applied to variant frequency time-series data. In a pioneering work, Newberry et al. (2017) applied the method of Feder et al. (2014) to assess the relative contributions from drift and selection in historical changes, and found that drift was a likely explanation in many cases. However, this analysis suffers from the same potential lack of sensitivity as was seen in the application to genetic data (Feder et al. 2014).

Furthermore, our approach lends itself to extensions that allow for the possibility that evolutionary parameters may change over time. In the section Regulated and unregulated language change, we demonstrated this in the context of regulated and unregulated change, showing that changes in selection strength that might reasonably be expected in regulated change are detected by our method, whereas no such changes were found in the the case of unregulated change.

A limitation of the method we have employed here is the assumption that sample sizes are large enough that the uncertainty on allele frequency estimates drawn from them can be neglected. In reality, this may not be the case, under which circumstances one would normally turn to a hidden Markov framework, as proposed by Bollback et al. (2008) in the context of genetic time series data. A naïve numerical implementation of this scheme would involve integrating over each of the (now hidden) frequencies xti in Equation (22), which dramatically increases the computational demands. One way to circumvent the additional integrals is to employ an expectation-maximization algorithm. However, we have found that if the effective population size is considered a free parameter, expectation-maximization tends to push this towards infinity due to piecewise deterministic trajectories being favored by the algorithm. Therefore, some further work is needed to develop tractable methods for jointly estimating effective population size and the selection coefficient when working with data subject to sampling uncertainty.

In the meantime, we have shown how one can account for known variation in sample sizes over the course of the time series, which is important when trying to determine if the effective population size (which governs fluctuations intrinsic to the evolutionary process) changes over time. The basic idea is to reduce the size of the larger samples so that the uncertainty due to sampling is then uniform across the time series. Although the resulting estimates of the effective population size then contain a contribution from both intrinsic fluctuations and sampling, any detected changes in the effective population size are most likely to arise from a change in the amplitude of the intrinsic fluctuations.

In summary, despite certain limitations, the method introduced here allows evolutionary parameters to be reliably estimated, and when combined with empirical likelihood ratio tests, can be used to test departure from a variety of null hypotheses. Although we have focused on the Wright–Fisher model with frequency-independent selection, it could be extended to models that involve other evolutionary processes. Extensions to processes involving more than two alleles are likely also possible in principle, although may involve higher-dimensional integrals that become difficult to perform numerically.

Supplementary Material

iyad092_Supplementary_Data

Data availability

The code, as well as data used in section Regulated and unregulated language change are available at https://datashare.ed.ac.uk/handle/10283/4811. Data used in section Beneficial mutations in yeast populations, originally obtained by Lang et al. (2011), is available in the supplementary materials of Feder et al. (2014). Supplemental material is available at GENETICS online.

Funding

JGM holds Principal’s Career Development Scholarship awarded by the University of Edinburgh. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.
==== Refs
Literature cited

Amato  R, Lacasa  L, Díaz-Guilera  A, Baronchelli  A. The dynamics of norm change in the cultural evolution of language. Proc Natl Acad Sci USA. 2018;115 :8260–8265. doi:10.1073/pnas.1721059115 30072428
Baxter  GJ, Blythe  RA, Croft  W, McKane  AJ. Utterance selection model of language change. Phys Rev E. 2006;73 :046118. doi:10.1103/PhysRevE.73.046118
Blythe  R, Croft  W. How individuals change language. PLoS ONE. 2021;16 :1–23. doi:10.1371/journal.pone.0252582
Bollback  JP, York  TL, Nielsen  R. Estimation of 2Nes from temporal allele frequency data. Genetics. 2008;179 :497–502. doi:10.1534/genetics.107.085019 18493066
Boyd  R, Richerson  PJ. Culture and the Evolutionary Process. Chicago: University of Chicago Press; 1988.
Casella  G, Berger  RL. Statistical Inference. 2nd ed. Boston: Cengage Learning; 2001.
Cavalli-Sforza  LL, Feldman  MW. Cultural Transmission and Evolution: A Quantitative Approach. Princeton: Princeton University Press; 1981.
Charlesworth  B . Effective population size and patterns of molecular evolution and variation. Nat Rev Genet. 2009;10 :195–205. doi:10.1038/nrg2526 19204717
Croft  W . Explaining Language Change: An Evolutionary Approach. London: Pearson Education; 2000.
Crow  J, Kimura  M. An Introduction in Population Genetics Theory. New York: Harper and Row; 1970.
Davies  M . The Corpus of Historical American English (COHA).  2010. https://www.englishcorpora.org/coha
Dominguez-Bello  MG, Blaser  MJ, Ley  R, Knight  R. Development of the human gastrointestinal microbiota and insights from high-throughput sequencing. Gastroenterology. 2011;140 :1713–1719. doi:10.1053/j.gastro.2011.02.011 21530737
Feder  A, Kryazhimskiy  S, Plotkin  J. Identifying signatures of selection in genetic time series. Genetics. 2014;196 :509–522. doi:10.1534/genetics.113.158220 24318534
Fisher  RA . The Genetical Theory of Natural Selection. Oxford: Clarendon Press; 1930.
Guzmán Naranjo  M . The se-ra alternation in Spanish subjunctive. Corpus Linguist Linguist Theory. 2017;13 :97–134.
Hui  TYJ, Burt  A. Estimating effective population size from temporally spaced samples with a novel, efficient maximum-likelihood algorithm. Genetics. 2015;200 :285–293. doi:10.1534/genetics.115.174904 25747459
Kempas  I . Sobre la variación en el marco de la libre elección entre cantara y cantase en el español peninsular. Moenia. 2011;17 :243–264.
Kimura  M . The Neutral Theory of Molecular Evolution. Cambridge: Cambridge University Press; 1983.
Kingman  JFC . The coalescent. Stoch Process their Appl. 1982;13 :235–248. doi:10.1016/0304-4149(82)90011-4
Kreitmann  M . Methods to detect selection in populations with applications to the human. Annu Rev Genomics Hum Genet. 2000;1 :539–559.11701640
Lacerda  M, Seoighe  C. Population genetics inference for longitudinally-sampled mutants under strong selection. Genetics. 2014;198 :1237–1250. doi:10.1534/genetics.114.167957 25213172
Lang  GI, Botstein  D, Desai  MM. Genetic variation and the fate of beneficial mutations in asexual populations. Genetics. 2011;188 :647–661. doi:10.1534/genetics.111.128942 21546542
Lang  G, Rice  D, Hickman  M, Sodergren  E, Weinstock  G, Botstein  D, Desai  M. Pervasive genetic hitchhiking and clonal interference in forty evolving yeast populations. Nature. 2013;500 :571–574. doi:10.1038/nature12344 23873039
Lenski  RE, Rose  MR, Simpson  SC, Tadler  SC. Long-term experimental evolution in Escherichia Coli. I. Adaptation and divergence during 2,000 generations. Am Nat. 1991;138 :1315–1341. doi:10.1086/285289
Lukić  S, Hey  J. Demographic inference using spectral methods on SNP data, with an analysis of the human Out-of-Africa expansion. Genetics. 2012. 192 :619–639.22865734
Michel  JB, Shen  YK, Aiden  AP, Veres  A, Gray  MK, Pickett  JP, Hoiberg  D, Clancy  D, Norvig  P, Google Books Team, et al  Quantitative analysis of culture using millions of digitized books. Science. 2011;331 :176–182. doi:10.1126/science.1199644 21163965
Newberry  M, Ahern  C, Clark  R, Plotkin  J. Detecting evolutionary forces in language change. Nature. 2017;551 :223–226. doi:10.1038/nature24455 29088703
Paris  C, Servin  B, Boitard  S. Inference of selection from genetic time series using various parametric approximations to the Wright–Risher model. G3 (Bethesda). 2019;9 :4073–4086. doi:10.1534/g3.119.400778 31597676
Press  WH, Teukolsky  SA, Vetterling  WT, Flannery  BP. Numerical Recipes: The Art of Scientific Computing. 3rd ed. Cambridge: Cambridge University Press; 2007.
Real Academia Española . Gramática de la lengua castellana. 27th ed. Madrid: Perlado, Páez y Cía;  1911.
Reali  F, Griffiths  TL. Words as alleles: connecting language evolution with bayesian learners to models of genetic drift. Proc R Soc B. 2010;277 :429–436. doi:10.1098/rspb.2009.1513
Reuter  JA, Spacek  DV, Snyder  MP. High-throughput sequencing technologies. Mol Cell. 2015;58 :586–597. doi:10.1016/j.molcel.2015.05.004 26000844
Sirén  J, Marttinen  P, Corander  J. Reconstructing population histories from single nucleotide polymorphism data. Mol Biol Evol. 2011;28 :673–683. doi:10.1093/molbev/msq236 20819907
Tataru  P, Bataillon  T, Hobolth  A. Inference under a Wright–Fisher model using an accurate beta approximation. Genetics. 2015;201 :1133–1141. doi:10.1534/genetics.115.179606 26311474
Tataru  P, Simonsen  M, Bataillon  T, Hobolth  A. Statistical inference in the wright–fisher model using allele frequency data. Syst Biol. 2017;66 :e30–e46.28173553
Wilks  SS . The large-sample distribution of the likelihood ratio for testing composite hypotheses. Ann Math Stat. 1938;9 :60–62. doi:10.1214/aoms/1177732360
Wright  S . Evolution in mendelian populations. Genetics. 1931;16 :97–159. doi:10.1093/genetics/16.2.97 17246615
