{'fraction_non_alphanumeric': 0.03716837404780667, 'fraction_numerical': 0.0112949829261886, 'mean_word_length': 4.208162334701322, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation of Mutation Rates from Fluctuation Experiments via Probability Generating Functions

This paper calculates probability distributions modeling the Luria-Delbr\"uck experiment. We show that by thinking purely in terms of generating functions, and using a 'backwards in time' paradigm, that formulas describing various situations can be easily obtained. This includes a generating function for Haldane's probability distribution due to Ycart. We apply our formulas to both simulated and real data created by looking at yeast cells acquiring an immunization to the antibiotic canavanine. This paper is somewhat incomplete, having been last significantly modified in March 29, 2014. However the first author feels that this paper has some worthwhile ideas, and so is going to make this paper publicly available.


Introduction
The famous experiment of Luria and Delbrück [5], [12] determined whether mutations in bacteria arise via Darwinian evolution or some kind of Lamarckian process. The experiment consisted of taking small samples of bacteria, and then allowing them to replicate up to a known large number of cells on many plates, and then looking at the distribution of number of cells on each plate that had acquired an immunity to an antibiotic. The mutation was carefully chosen to be a forward mutation, that is, the probability of the mutation taking place was very much larger than the probability of the mutation reversing itself.
If the immunity was acquired by some kind of Lamarckian process, then one would naturally expect the distribution of surviving cells in the plates to follow a Poisson distribution. In particular, the variance would have a value very close to the mean, and the probability of a plate having a so called 'jackpot,' that is, a much larger number of surviving cells than the mean, would be vanishingly small.
However, Luria and Delbrück argued that if the immunity was acquired by Darwinian evolution, then the variance of the number of surviving cells would be much larger than the mean, and furthermore, the probability of any plate having a 'jackpot' would be large enough that it would be observed reasonably often.
Their experiment conclusively validated the latter assumption. Soon after this, scientists sought after a formula for the distribution of surviving cells in each plate. According to [11], the biologist J.B.S. Haldane produced a formula for the distribution in 1946. But his work was never widely published at that time, and in 1949 Lea and Coulson [4] produced a different distribution, which seems to be the basis of most research these days.
A great advantage of their distribution is that it can be given by a simple generating function, whereas the calculations described by Haldane seem to be very time consuming, and requires enumerating combinatorial structures. For example, the paper [6] gives an algorithm whereby the Lea-Coulson distribution can be calculated quickly. The survey paper [9] gives many ways in which one can calculate the fluctuation rate, that is, the probability of a single offspring cell acquiring the mutation. And the paper [1] gives a Bayesian approach to estimating the fluctuation rate.
The primary goal of this paper is to promote an approach which uses generating functions from the very beginning. The other change we propose to the approaches taken in other papers is to think 'backwards in time' rather than 'forwards in time,' that is, instead of considering how mutants might have developed from initial conditions, look at the final situation and reason out where the mutations must have come from. We will illustrate this approach under several conditions. In particular we will be able to obtain a generating function for the Haldane distribution which allows for rapid calculations.

The experimental setup to be modeled
We assume that in each plate we start from a small number of cells, none of which have the mutation. Then the cells replicate repeatedly. When each cell replicates, one of the offspring can acquire the mutation with a probability µ, which is extremely small, effectively infinitesimal. We continue the replication process until we have obtained a population of n cells, where n is effectively infinite. We assume that m = µn is of order 1.
We shall use the following terminology (which is common in the literature). We shall call a cell that acquires the mutation, but not from its parent, a mutation. A cell that has the mutant gene will be called a mutant. We will assume that the probability of a mutant have offspring that are not mutants is zero. This can either come from supposing that the genetic change is much more likely to happen than to be reversed (that is, it is a forward mutation), or by realizing that the population of mutants is so much smaller than the number of nonmutants, that even if the probability of a mutant having non-mutant offspring were of the order µ, we wouldn't observe this anyway.) Thus all mutations are mutants, but not vice versa.
The quantity m is the fundamental parameter that describes the probability distribution of the random variable R which is the number of cells that have acquired the mutation. All the formulas described in this paper will be of the form where m = µn is the fundamental parameter that we try to estimate, α is a real number, and q r (m) is a polynomial of degree r.
We assume that all cells take approximately the same amount of time between cell divisions. We use the letter g to denote the number of generations backwards in time, setting g = 0 to denote when the experiment finishes. Under this assumption, the cell divisions that took place recently (that is, with g small), are assumed to take place synchronously. However the bulk population may or may not be assumed to be dividing asynchronously. We give each mutant a generation number g ≥ 1, which says this mutant was created g generations ago. In future papers we will drop the assumption of the same amount of time between cell divisions.
Another assumption we make is that the population of non-mutant cells is effectively infinite. Thus the number of non-mutant cells is effectively the same as the number of all cells g generations ago, and does not depend upon the number of mutations created.
Finally, we will allow for the possibility that cells may die or become non-functional. If cell death does not occur, then the total number of cells g generations ago is 2 −g n. Figure 1 illustrates the process looking only four generations back, but in our formulas we assume that the total number of generations is effectively infinite.
Haldane produced a distribution for the Luria-Delbrück experiment based on the assumption that all the cells divide with perfect synchronicity, and no cells die. The formula can be described thus: given an integer r ≥ 0, let P r denote the set of sequences (a 0 , a 1 , . . . ) of non-negative integers, with only finitely many non-zero terms, such that ∞ s=0 a s 2 s = r. Then (2.2) (Note that [11] has a typographical error in stating this formula. Note also that his g is our m/2.)

Generating functions for probability distributions
If R is a random variable that takes values in the non-negative integers, then the probability generating function of R is defined by the formula Note that G R is an analytic function, with radius of convergence at least as big as 1. We will make much use of the following well known results.
be a sequence of non-negative integer valued independent random variables.
1. If the sum Z = ∞ k=1 X k converges almost surely, then 2. If X k are identically distributed, and N is a non-negative integer valued random variable that is independent of the X k , and if Z = N n=1 X k , then 3. If N is a Poisson random variable with mean λ, then All the distributions we shall consider will have a generating function of the form where α > 0, and H R (x) is an analytic function satisfying H R (0) = 0 with Taylor series For example if the mutation acquisition were Lamarckian, then R would have the Poisson distribution with mean m, that is and it is well known that the generating function is Another example is the Lea-Coulson distribution [4]. This assumes that the mutation acquisition is Darwinian, but it also approximates the discrete replication process by a continuous process. They obtain a distribution that satisfies (3.5) with α = 1 and There are two ways to compute the polynomials q r from the power series for H R (x). One way is to place H R (x) into the Taylor's series for e x . This involves multiplication of polynomials. We discuss this method in Section 10.
The other algorithm is described in Ma, Sandri, and Sarkar [6]: In practice, this algorithm worked very well. Note that to compute q r (m), one only needs to know h s for s ≤ r.

The generating function approach
The approach adopted throughout this paper is to compute the generating function for R, the number of mutations. Let us first illustrate this method to derive the Haldane distribution. We assume the cell division is completely synchronous, and that no cells die or malfunction. While this paper was being prepared, we discovered that the following result is a special case of a result that appeared in 2013 [15].
Hence from the Ma-Sandri-Sarkar algorithm equation (3.10) we obtain the rapid formula From equation (4.2) we obtain the polynomials given in Table 1.
Proof. We know that g generations ago, the population count is n/2 g . All of these were created from cells dividing g + 1 generations ago, and therefore half of these have copied genotypes, the rest having the original genotype. Thus the number of mutants created at  that time is N g , a Poisson random variable with parameter 1 2 µn/2 g = 2 −g−1 m. The number of mutations that arise from the gth generation mutants is 2 g N g . Hence Furthermore, because we suppose the total population of cells to be very large, we can assume that the random variables N g are independent. Thus we obtain the generating function This approach can be generalized as follows. First, we define times t k , (k = 0, 1, 2, . . . ), with 0 = t 0 < t 1 < t 2 < . . . . The units of time will be generations, and as with generation number, a positive number denotes the number of generations before the experiment is completed. Let us suppose that a time t k , (k ≥ 1), and only at those times, some or all of the cells divide. For each k, we let n k denote the total population of cells at time t k , and let π k = n k /n. Thus 1 = π 0 ≥ π 1 ≥ π 2 ≥ · · · > 0. Because we have assumed that the total population of cells is effectively infinite, we can assume t k → ∞ as k → ∞.
We consider two random variables: 1. for each possible k ≥ 0, the number of mutants N k created at negative time t k , which will always be a Poisson random variable with parameter π k − π k+1 ; 2. for each possible k ≥ 0, the number of mutations M k that arise from any one mutant which was created at time t k .
Then the total number of mutations will always be given be the formula where M (l) g denotes independent copies of M l . Hence we obtain the generating function Let's explain first how the Haldane distribution fits into this paradigm. Since cell division is synchronous, we set t g = g where the time is measured in multiples the time between cell divisions. Since no cells die, we have that π g = 2 −g , and M g = 2 g , so that G Mg (x) = x 2 g . Substituting into equation (4.6), we obtain equation (4.1). We can also obtain the Lea-Coulson distribution using this paradigm. Note the argument we present is in their paper as a "second proof" (and indeed the whole "generating function approach" is inspired by their paper). We still suppose that no cells die, but we assume that the formulas π k = 2 −t k and M k = 2 t k still hold when t k is not an integer. To get the most even spacing of the generating of cells, we set t k = log 2 (k + 1). Then π k − π k+1 = 1/(k + 1) − 1/(k + 2) = 1/(k + 1)(k + 2), and M k = k + 1. Thus equation (3.8) follows after substituting the summation variable k by k − 1.

Haldane's distribution versus the Lea-Coulson distribution
The Lea-Coulson assumes that the generations are somehow being created continuously. From this point of view, it would seem that their distribution is unrealistic. However there also are problems with the Haldane distribution. It is obvious that the probability of obtaining 128 mutants is far higher than the probability of obtaining 127 mutants. If one had a single mutation, say, seven generations ago, this will result in 128 mutants. But it is possible that a few of these mutants will either die (and perhaps any of their mutant ancestors died), or fail to be transferred to the plate properly. Thus while Haldane's distribution puts a very low probability of getting 127 mutants compared to 128 mutants, nevertheless from an experimental point of view, it is very likely that not all of the 128 mutants will be observed, and the chances of seeing 127 or less mutants is comparable with the probability of seeing 128 mutants.
To show this effect, in Figure 2 we give log/log plots of p r (m) to r, for the value m = 2, for both distributions. On the x-axis we show log 2 (r), and on the y-axis we show log 10 (p r (2)). Notice that for the Lea-Coulson distribution that this graph approaches a straight line as r → ∞. In fact, in [6], it is shown that p r (m) ≈ c m r −2 as r → ∞, where c m depends only on m.
We will look at several ways to modify the Haldane distribution.
1. Allow that some cells may die.
2. Examine the situation when only a fixed proportion of cells are plated with the antibiotic.
4. Allow that the mutants replicate at a different rate than the non-mutants [3].
We will explore the third case only superficially, as we have not yet obtained formulas that fully handle this issue.

Accounting for cell death
Let us now assume that a cell dies with probability θ before it replicates. Thus from one cell, the expected number of cells that come from this cell one generation later is 2(1 − θ). Thus π g = (2(1 − θ)) −g (6.1) The computation of G Mg is given by a recursion relationship. From a mutant created g generations ago, with probability θ this mutant will die, and with probability 1 − θ there will by X + Y cells, where X and Y are independent and have the same distribution as M g−1 .
Thus we obtain These formulas are then substituted into equation (4.6) to obtain equation equation (2.1) with It is unlikely we will find an explicit formula for solving the recurrence relation (6.3), because this is in essence the same recurrence relation that is used to define the famous Mandelbrot set [7].

Plating only a fixed proportion of the cells
This is discussed for the Lea-Coulson distribution in [13,14].
Let us now assume that we perform an experiment such that the probability of obtaining R = r mutants is given by p r , with corresponding generating function F R (x) = e −αm exp(mH R (x)). Now suppose that we only plate a proportion θ ∈ (0, 1] of the cells. How many mutants will we observe? Let us call the number of mutants observed R θ .
Theorem 7.1. The random variable R θ has generating function Proof. We have that R θ has the same distribution as R r=1 I r , where I r is a sequence of independent random variables taking the value 1 with probability θ, and the value 0 with probability 1 − θ. The result follows from applying Proposition 3.1.
Another way to obtain this formula is to use the binomial distribution and then rearrange the resulting double sum. Usually, the most timely and accurate numerical method for computing α θ , and coefficients of the power series for H θ , is by direct evaluation from the power series of the derivatives of H(x) at x = 1 − θ. For the Lea-Coulson probability measure, it makes sense to compute α directly by substituting equation (3.8). One might also be tempted to compute the coefficients of H θ by symbolically differentiating H R given in equation (3.8). However the symbolic derivatives become very unwieldy.

Asynchronous cell division
There are two notions of asynchronous cell division that we shall consider.
1. slight asynchronous cell division We suppose that cell division is asynchronous enough so that at time g generations before the experiment ends, the cell population is given by 2 −g n, even if g is not an integer, but not so asynchronous so that over a few generations the cell division is well approximated by synchronous cell division.
2. full asynchronous cell division We suppose that cell division is sufficiently asynchronous that we should base our model on a probability distribution that describes the time for a cell to divide.
Accounting for full asynchronous cell division seems to be very difficult, and we have not yet solved this problem 1 . In this chapter we shall consider slight asynchronous cell division, and show that it makes no difference to the Haldane distribution. We assume that cell splitting takes place every g/n time units, where at the end we let n → ∞. Therefore in equation (4.6), we set t k = k/n, π k = 2 −k/n , and G M k = 2 [k/n] , where [x] denotes the integer part of x. Then writing k/n = g + j/n where g = [k/n], we obtain where the last step comes from summing the geometric series. Thus the formula is independent of n.

Mutants replicating at a different rate than non-mutants
The analogous formula to the Lea-Coulson was calculated by [3]. We will do the analog for the Haldane distribution. While this paper was being prepared, we found out that the this result appeared in 2013 [15].
Let g denote the number of generations back in time as counted in the time for a mutant to replicate. Suppose that the time for a non-mutant to replicate is τ times the time for a mutant to replicate. Then we can apply equation (4.6) with M g = 2 g , and π g = 2 τ g , to obtain equation equation (2.1) with 10 Comparison of the Ma-Sandri-Sarkar algorithm with using the Fast Fourier Transform Another way to compute the probability generating function given by equation (3.5) is to expand it using Taylor's series for e x , and obtain In order to calculate p r (m), it is only necessary to sum the series to the rth term, and to calculated [H R (x)] n up to x r . It is well known that [H R (x)] n can be rapidly calculated using the Fast Fourier Transform. If one performs a algorithm complexity analysis, one finds that the time to compute p r (m) for m ≤ n is O(n 3 ) using the Ma-Sandri-Sarkar algorithm, but is O(n 2 log n) when using the Fast Fourier Transform. (The exception is when using the Ma-Sandri-Sarkar for Haldane's distribution, when the time taken is also O(n 2 log n).) A difficulty with the Fast Fourier Transform is getting sufficient accuracy in the coefficients of G R (x). If one calculates [H R (x)] n using long multiplication (and this is, in effect, what the Ma-Sandri-Sarkar algorithm is doing), then because all the coefficients of H R (x) are positive, there will be no canceling large numbers to produce small numbers in computing the coefficients. However, the Fast Fourier Transform combines the various coefficients using complex numbers, and so the accuracy of the final answer is limited by the size of the largest coefficient.
However there is also the fast polynomial package FLINT [2], which comes standard with the Sage software package [10]. While all our computations were performed with Sage, we found that in our situations that the Ma-Sandri-Sarkar was much faster than using the built in polynomial multiplication.

Analyzing data using likelihood functions/Bayesian statistics
In this section, we will show how to compute the likelihood function to estimate the value of m from experimentally obtained data. Suppose we have performed n separate experiments, using identical but independent protocols, and obtain counts of plates r 1 , r 2 , . . . , r n . Then the likelihood function is where R has the distribution that we believe most correctly represents our situation. The Bayesian approach is to suppose that L(m) represents a probability distribution for a random variable M, where Pr(M ∈ [m 1 , m 2 ]) represents the 'belief' that we have that the parameter m lies between m 1 and m 2 . This is related to L(m) via the formula where here we think of dm is being infinitesimal. Here f (m)dm is usually called the prior distribution, and the constant C is chosen so that the integral of the right hand side is zero. The paper [1]  However, since we are interested in graphing these using a log scale on the x axis, to get an idea of which range of m we should use for plotting, we want to compute E(log(M)) and var(log(M)). These are not so easily computed, but numerical integration does an outstanding job.
In any case, we would hope that the amount of data collected should be enough so that the formulas E(M) or exp(E(log(M))) or the maximum of L(m), should give similar answers, even if different prior distributions are used. And we would hope to see a bellshape distribution, so we should get similar answers for the estimate of the middle 68.3%-ile of belief, that is, E(M) ± var(M) or exp(E(log(M)) ± var(log(M))).
Finally, in performing the computations, the numbers involved tend to be much smaller than much computer software can represent. For example, the commonly used double precision IEEE floating point representation cannot represent numbers much smaller than 10 −308 , and values far smaller than this will often arise in calculating the likelihood function. However we used the software package sage [10], and this was well capable of handling very small numbers.

Simulations 13 Analysis of real data
The following experiments were performed by students who were part of the 'Mathematics in Life Science' program at the University of Missouri.
The mutations in question confer resistance to an antibiotic (canavanine) in yeast cells. There were two types of resistant mutants, one giving red colonies, the other giving white colonies. Our experiment focused on the red mutants, because we were able to sequence their DNA to find the exact mutation responsible for resistance.
In Experiment A, the students plated 70 60μl cultures, each arising from a 25ml culture. In Experiment B was to plate 50 60μl cultures from a single 25ml culture. Each culture started at a cell density of 10 5 cells/ml and ended at a cell density of 10 8 cells/ml. See Figure 4.
The counts for the 50 60μl cultures suggest we hit a jackpot. However in looking at the distribution created by the Bayesian method, the data seems to rather confidently give a large figure for m, which we strongly suspect is too large. One way to interpret this data is   to use a log-log plot for the likelihood function. In this case it can be seen ...... issues with how to interpret probabilities obtained using the Bayesian method.