Relevant Assumptions

The fluctuation experiment is conceptually simple. Although neither mutations nor mutants are directly observable, a mutant can be made visible by allowing it to form a colony on a solid culture. As a result, analysis of fluctuation assay data revolves around inferring the number of mutations from the number of mutants in a test tube. This challenging task relies critically on mathematical models that bridge the gap between mutations and mutants. Like mathematical models for any other purposes, a mathematical model for fluctuation assay data depends on simplifying assumptions. Acquiring a sound understanding of the key assumptions not only helps experimentalists to better design experiments, but also enables them to analyze their data more confidently. Major assumptions underlying common mathematical models for fluctuation assay data include the following.

• A1. Cells (nonmutants and mutant alike) undergo unimpeded growth in a liquid culture (contained in a test tube). In other words, cells are in logarithmic phase (also known as exponential phase).

• A2. When a cell divides, the probability is a constant p that one of the two daughter cells is a mutant. The probability is 1 − p that both daughter cells are nonmutants. Therefore, there is no possibility of a nonmutant splitting into two daughter mutants.

• A3. Mutants grow at the same rate as nonmutants. As a consequence, the relative fitness is unity.

• A4. Back mutation is negligible, and hence is not considered.

• A5. Cell death is negligible, and hence is not considered. As Kendal and Frost (1988) rightly asserted, this assumption is not a problem when dealing with bacterial cells (in contrast to somatic cells). But there is an important exception. When the Luria-Delbrück protocol is used to study antibiotic-induced mutations, as in the study of Cairns et al. (2017), cells must grow in the presence of an antibiotic. The dose administered should be sufficiently low so as not to cause considerable death of nonmutants. Otherwise, owing to death of large numbers of nonmutants, unimpeded growth of antibiotic-resistant mutants will grossly inflate the mutant frequency. This caution should not be confused with the use of a lethal dose of an antibiotic after the cells are plated, as a selective agent.

• A6. Every mutant exiting prior to plating is capable of growing into a colony after the plating process. In other words, the plating efficiency is perfect.

• A7. Cell life span obeys the negative exponential distribution.

• A8. The total numbers of nonmutants immediately before plating, denoted by N_t, are approximately the same across tubes within an experiment.

• A9. Nonmutants die immediately after plating. In other words, selection after plating is so lethal that postplating mutation is negligible. Note that postplating mutation was of interest to Lang and Murray (2008), and to Ford et al. (2013).

• A10. When a nonmutant undergoes a mutation, the resulting daughter mutant is capable of immediately manifesting the mutation’s phenotype. This assumption thus excludes the possibility of phenotypic delay, a concept that can be traced to Newcombe (1948).

• A11. N₀, the number of cells used to seed a test tube, is small enough to guarantee that the N₀ seeding cells contain no mutants.

A useful mutation model yields a so-called mutant distribution, which describes the probability distribution of the number of mutants, and which usually involves a parameter customarily denoted by m, the expected number of mutations per tube. Once an accurate estimate of m is inferred from fluctuation assay data using an appropriate mutant distribution, the mutation rate is determined by(1)For a detailed discussion of the mutation rate p, the reader is referred to Zheng (2017c).

Almost by definition, a simplifying assumption is a deviation from biological reality. One is doomed to disappointment when checking any of the above assumptions outside the context of what the fluctuation experiment aims to achieve. As the point of the fluctuation experiment is to estimate the fundamental parameter m, the validity of the assumptions (A1)–(A11) should be assessed by whether they allow acceptable estimates of m to be extracted from fluctuation assay data. Although limited research has been conducted in this regard, the assumptions seem to be less worrisome than originally thought. For example, the assumption (A7), the most contrived at first glance, was recently found to be an acceptable assumption.

A peculiar feature of all mutant distributions is a lack of explicit analytic expressions. A mutant distribution must therefore be identified by its probability generating function (PGF). If Y denotes the random number of mutants, the PGF of Y is the function defined by . rSalvador encompasses five mutation distributions. The first mutant distribution is due to Lea and Coulson (1949). This distribution depends on all the above 11 assumptions, and its PGF is(2)Here, is strictly less than unity (Zheng 2002). But φ is assumed to be unity by convention. Strictly speaking, the Lea-Coulson mutant distribution refers to the special case where the PGF is of the form (3)Nádas et al. (1996) were the first to suspect that setting ∅ = 1 could have adverse effects on the estimation of m. However, in practice, ∅ > 0.999 almost always holds true. Under such a condition, the effects of ∅ are negligible. rSalvador is the only software tool that permits the user to specify a value smaller than unity for ∅, mainly for educational purposes.

The second mutant distribution resulted from an attempt to relax the assumption (A3) by Mandelbrot (1974). Because Koch (1982) made a similar attempt to allow for differential growth rates between nonmutants and mutants, the mutant distribution is usually called the Mandelbrot-Koch distribution. The following PGF is a generalization of the PGF

.(4)Here B denotes the usual beta function. The parameter w is relative fitness, that is, the ratio of the growth rate of mutants to that of nonmutants.

The third mutant distribution can be traced to Armitage (1952), who was the first to attempt to relax the assumption (A6). The resulting mutant distribution was further studied by Stewart et al. (1990). The PGF of that distribution is

(5)Here, ξ = / (1− ). The parameter ϵ is called the plating efficiency, which denotes the portion of culture plated or the probability that a mutant forms a visible colony after the plating process.

The fourth mutant distribution is a special case of the relatively new distribution (Zheng 2010). The distribution has two parameters, A and k, and is denoted by The PGF of a variable is

(6)The distribution’s usefulness lies in its ability to relax the assumption (A8).

The last mutant distribution derives from Haldane’s mutation model (Zheng 2005). Because cell growth is synchronous in the Haldane model, cell life spans are nonrandom. This is in violation of the assumption (A7). The Haldane mutant distribution is indexed by the mutation rate p and the number of elapsed cell generations g. rSalvador includes the Haldane model mainly for educational purposes, and it relies on algorithms developed by Zheng (2007).

Classical Analysis

In most published applications of the fluctuation experiment, inference about m relied on the Lea-Coulson mutant distribution (Lea and Coulson 1949) defined by the PGF given in Equation (3). The Lea-Coulson distribution owes its popularity partly to the work by Ma et al. (1992). SALVADOR 1.0 (Zheng 2002) was the first publicly available software package allowing the user to calculate maximum likelihood (ML) estimates of m under the Lea-Coulson model. The web tool FALCOR (Hall et al. 2009) rendered this capability more easily accessible to bench scientists. But FALCOR does not automate the process of computing confidence intervals (C.I.s) for m. As a result, the user is forced to compute C.I.s outside FALCOR by following instructions given on FALCOR’s website.

rSalvador has a unique advantage. Since SALVADOR 2.0 (Zheng 2005), the expected Fisher information has been replaced by the observed Fisher information in all algorithms for computing point and interval estimates of m. As a consequence, all C.I.s produced by SALVADOR (or rSalvador) are based on the likelihood ratio method. Efron and Hinkley (1978) were among the first to observe that likelihood ratio C.I.s are preferable to the commonly used Wald-type C.I.s. Raue et al. (2009) referred to likelihood ratio C.I.s as finite sample C.I.s, to underscore the superiority of likelihood ratio C.I.s over Wald-type C.I.s. rSalvador (including SALVADOR) is so far the only fluctuation assay tool that advocates likelihood ratio C.I.s. One distinctive feature of the likelihood ratio C.I. for m is that it is usually asymmetric about the ML estimate of m. As an example, consider the well-known Demerec data (Demerec 1945). The ML point estimate of m is

> newton.LD(demerec.data)

[1] 10.84383

But a 95% likelihood ratio C.I. for m is

> confint.LD(demerec.data)

[1] 8.650538 13.194765

The center of this C.I. is at 10.92, not 10.84—the ML estimate of m. Unlike the Wald-type C.I.s, a likelihood ratio C.I. requires iterative computing. The user can view the iterative process by setting the option show.iter = T in the above command.

As an educational feature, rSalvador allows the user to specify an arbitrary value between 0 and 1 for the parameter φ given in Equation (2). In the Demerec experiment, N₀ = 90 and N_t = 1.9 × 10⁸ (Zheng 2002). However, if the experiment were terminated prematurely at N_t = 900, one would have φ = 1 − 90/900 = 0.9. A noticeably different estimate of m would result, as predicted by Nádas et al. (1996).

> newton.LD(demerec.data,phi = 0.9)

[1] 18.92394

The simple Lea-Coulson model relies on each of the assumptions (A1)–(A11), which may belie the model’s practical usefulness. One source of skepticism is the assumption (A7), as cell life spans do not obey the negative exponential distribution. It was this assumption that triggered the rhetorical question “Can estimates be trusted?” In a study of the Haldane model, Zheng (2007) compared the Lea-Coulson model with the Haldane model and found that the two models produced comparable estimates of mutation rates. This finding suggests that the assumption (A7) is acceptable, because the two models assume drastically different cell life distributions. A simulation study by Ycart (2013) reached essentially the same conclusion. Still, Ycart (2013) suggested using a cell life distribution derived from a particular yeast study as a universal cell life distribution. In response to Ycart’s suggestion, Gillet-Markowska et al. (2015) emphasize that “there is no universal cellular division time model as it depends on experimental conditions like the strain or the media.”

Accounting for Relative Fitness

Classical analysis of fluctuation assay data excludes the possibility of differential growth rates between mutants and nonmutants. Recognition of the importance of allowing for this possibility can be traced to Lieb (1951), who highlighted the P₀ method of Luria and Delbrück (1943) for its capability to accommodate differential growth rates. Using the mutant distribution defined by the PGF in Equation (4), Zheng (2002) developed algorithms for computing ML estimates of m and w. These algorithms rely on the expected Fisher information that must be computed by truncating infinite series. Because computing infinite series could be prohibitively expensive in practice, Zheng (2005) improved these algorithms by replacing the expected Fisher information with the observed Fisher information, which not only obviated the problem of infinite series but also made possible the computation of likelihood ratio C.I.s. A comparison between Equations (32) and (46) in Zheng (2005) clarifies the distinction between the expected and the observed Fisher information. The elegant work of Efron and Hinkley (1978) provides helpful information on this topic.

However, Hamon and Ycart (2012) drew renewed attention to the Mandelbrot-Koch model by making the following observation.

Indeed, the likelihood and its derivatives can be computed by iterative algorithms: theoretically at least, the problem could be considered solved. This is not so in practice, mainly because the multiple sums that must be computed by the optimization algorithm make it quite unstable. According to the numerous tests that we have made, the ML estimates cannot be reliably computed for samples whose maximum exceeds 1000.

Such a disastrous outcome is unreproducible, as it was an artifact of the computer code of Hamon and Ycart. The improved algorithms for computing ML estimates given by Zheng (2005) had been thoroughly tested using simulated fluctuation experiments in which the maximum mutant counts far exceeded 5000. A diagnostic clue to the enigma is on p. 1258 of Hamon and Ycart (2012) below an expression for the expected Fisher information, denoted by I (α, ρ).

Fortunately the partial sums increase, so that when computing the inverse I⁻¹(α, ρ) the sum of the first m terms yields conservative confidence intervals; yet we do not consider it satisfactory.

Unaware of the advantages of using the observed Fisher information, Hamon and Ycart (2012) still used algorithms reliant on the expected Fisher information to compute ML estimates. This then-outdated practice led to an exaggerated sense of disaster, so Hamon and Ycart recommended the generating function (GF) method as a remedy. They envisioned future experiments in which cultures would contain enormous numbers of mutants to allow the experimentalist to harness the GF method’s potential. This would imply the abandonment of the Luria-Delbrück protocol, as large numbers of mutants must be directly counted, without the aid of plating. Such views caused confusion. For example, a simulation study of the GF method led Gillet-Markowska et al. (2015) to emphasize the importance of avoiding exceedingly large mutant counts.

[T]he precision on the estimation of m (and by consequence the estimation of μ) is higher for the smallest values of m. Therefore, users should not outgrow the cultures in order to limit the number of mutants that grow on selective plates.

The exploration of the GF method by Hamon and Ycart (2012) offers valuable lessons, which I here illustrate using the experimental data given in Table 3 of Rosche and Foster (2000). Eight of the 60 mutant counts were excluded in the analysis by Hamon and Ycart (2012), but I shall use the intact data. One reason for choosing this data set is that it contains the largest mutant count ever reported—3000. At least three versions of the GF method are available: that given by Hamon and Ycart (2012) via the R code attendant to their paper, that given by the web tool bz-rates (Gillet-Markowska et al. 2015), and that given by the R package flan (Mazoyer et al. 2017). ML estimates of m and w, along with estimates produced by the three versions of the GF method, are displayed in Figure 1. Judging by the contours of the log-likelihood function, the ML method outperforms all three versions of the GF method. The large mutant count did not derail the ML method, contradicting the observation of Hamon and Ycart (2012).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Contours of the log-likelihood function for the data given in Table 3 of Rosche and Foster (2000). The four colored dots represent estimates of m and w. The green dot indicates the ML estimates, the blue dot indicates the GF estimates via flan, the purple dot indicates the GF estimates via bz-rates, and the red dot indicates the GF estimates by an earlier implementation of the GF method.

> newton.joint.MK(cairns.foster.data)

[1] 1.3027909 0.7281044

Likelihood ratio C.I.s for m and w can be computed separately as follows.

> confint.profile.m(cairns.foster.data)

[1] 0.9855115 1.6749828

> confint.profile.w(cairns.foster.data)

[1] 0.5209636 1.0298620

The exploration of the GF method drew attention to the issue of parameter identifiability. Only recently did investigators in systems biology and related fields recognize practical identifiability as a problem distinct from the better-known problem of structural identifiability (Raue et al. 2009). Practical identifiability considers the effects of the amount and quality of data. Computational difficulties encountered with the Mandelbrot-Koch mutant distribution may also be viewed as practical nonidentifiability. The inability to estimate w in the examples given by Mazoyer et al. (2017) is indicative of practical nonidentifiability. In addition, the C.I.s that extend to zero in these examples are a symptom of the same problem. No multi-parameter mutant distribution seems to be immune to practical nonidentifiability. Angerer (2001) was perhaps the first to document computational difficulties that were due to practical nonidentifiability. In my own simulation exploration of the mutant distribution used by Lang and Murray (2008) and by Ford et al. (2013), I also found practical nonidentifiability a stubborn problem.

A motivation for studying the GF method was to improve the estimation of the fitness parameter w. Simulation suggests that it is uneconomical (in terms of sample size) to use the Luria-Delbrück protocol as a tool for measuring bacterial fitness. The traditional fitness assay (also known as the competition assay) is a more efficient tool, which has been increasingly conducted in tandem with fluctuation assays in recent studies (Scanlan et al. 2015). The Mandelbrot-Koch mutant distribution enables the investigator to incorporate information on w in the analysis of fluctuation assay data. For instance, if a fitness assay yields an estimated w of 0.73 in the above example, one obtains an estimate of m by

> newton.MK(cairns.foster.data,w = 0.73)

[1] 1.30215

In the absence of attendant fitness assays, the investigator may use joint ML estimation cautiously as a tool to obtain an estimate of m adjusted for the nuisance parameter w. In doing so, the issue of parameter identifiability should be kept in mind. While the ML approach is unsuitable for unrealistically large mutant counts, it remains the method of choice for real-world applications. An experimentalist is rarely able to count 500 mutant colonies per culture (Boe et al. 1994), let alone 5000 mutant colonies per culture. The experimentalist resorts to partial plating if a poorly designed experiment yields unwieldy mutant counts.

Accounting for Plating Efficiency

When an experimentalist plates only a portion of each culture, we refer to the actually plated proportion as the experiment’s plating efficiency, often denoted by ϵ. A culture before plating is densely populated by cells, but mutants constitute only a minute fraction of the cell population. As a result, mutants are sparsely and randomly dispersed in the culture. If 10% (i.e., ) of a culture is plated, in general, the experimentalist does not transfer exactly 10% of the mutants onto a selective culture—the actual number of mutants transferred is a random number. Although a liquid culture resides in a three-dimensional tube, the randomness induced by partial plating is more easily understood by analogy with a slice of raisin bread, as illustrated by Figure 2.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

This diagram draws an analogy between a three-dimensional liquid culture and a slice of raisin bread, to help explain the randomness induced by partial plating. The green points, each symbolizing a mutant, are randomly dispersed in the square area. The two circles indicate two possible ways of sampling (plating) an equal portion of the square area that represents a three-dimensional liquid culture.

Partial plating is common, as it is a practical measure to circumvent overwhelmingly large mutant counts. For example, except for one experiment, all experiments reported by Luria and Delbrück (1943) had a plating efficiency smaller than unity. Their experiment 16 is perhaps the best known by biology students. In that experiment a portion of 0.08 ml was plated from each of the twenty 0.2-ml cultures. Hence, the plating efficiency was . The popular genetics textbook of Griffiths et al. (2000, p. 481) missed the partial plating feature of that experiment. A recent tutorial (Meneely 2016, p. 373) also mistakenly omitted the partial plating feature of the same experiment. Because 11 of the 20 cultures had no mutants, the tutorial argues by the P₀ method that the mean number of mutations per culture is m ∼ 0.6.

Why is the above claim incorrect? There is a crucial distinction between the number of zero mutant counts and the number of actual null cultures—cultures having no mutants prior to plating. Suppose a culture contains one mutant prior to plating. When 40% of that culture is plated, the mutant has a 60% chance of not being counted. Thus, a zero mutant count may result from a nonnull culture. In other words, the experimentalist is uncertain about the number of actual null cultures. To account for this additional source of uncertainty, one can apply a modified P₀ method, which was proposed by Stewart et al. (1990, p. 184).

> p0.LD.plating(luria.16.data,e = 0.4)

[1] 0.9786801

Hence, the mean number of mutations per culture is not 0.60. Luria and Delbrück reported that each culture contained 5.6 × 10⁸ bacteria. Is this the number of cells in the plated portion or in the whole culture? To remove the ambiguity, note that Luria and Delbrück also reported an average mutant count of 28.4 mutants per culture. Because the mean number of the observed mutants in the plated portion was 11.35, the meaning of “per culture” is now clear. There must be 5.6 × 10⁸ cells in the whole culture. That is, N_t = 5.6 × 10⁸. Therefore, the mutation rate is p = 0.98/N_t = 1.75 × 10⁻⁹ mutations per cell division, according to Equation (1).

rSalvador relies on the mutant distribution defined by the PGF in Equation (5) to compute ML estimates and likelihood ratio C.I.s for the fundamental parameter m. It employs the algorithms developed by Zheng (2008). For example, an ML estimate of m for Luria and Delbrück’s experiment 16 is

> newton.LD.plating(luria.16.data,e = 0.4)

[1] 1.18636

and an LR-based 95% C.I. for m is

> confint.LD.plating(luria.16.data,e = 0.4)

[1] 0.5803079 2.0908012

Many recent investigators were unappreciative toward the effects of partial plating. For example, in the reanalysis of the experiments of Luria and Delbrück (1943), Hamon and Ycart (2012, Table 1) and Ycart (2013, Table 3) ignored the effects of partial plating. In the reanalysis of data from David (1970) and from Werngren and Hoffner (2003), Ycart and Veziris (2014) also disregarded the effects of partial plating. Moreover, unusual care may be required to determine the plating efficiency of an experiment that involves sophisticated procedures such as culture condensation. For instance, the plating efficiency in the study by Werngren and Hoffner (2003) was 0.4, as explained in Zheng (2016a, p. 357). But at first glance the plating efficiency may seem to be 0.2, as was assumed by Mazoyer et al. (2017).

Investigators who were aware of the effects of partial plating often relied on an adjustment due to Stewart et al. (1990), as this was the only approach recommended on the FALCOR web site. To apply the Stewart approach, the experimentalist first computes an initial estimate m* of m by assuming perfect plating efficiency.

> newton.LD(recent.expt)

[1] 15.36162

Because the adjustment factor is , an estimate of m is × A = 15.36 × 21.5 = 330.24. Note that this estimate of m is ∼56% of the corresponding ML estimate (593.6). The reason for this inadequacy of the Stewart approach was given by Zheng (2015b). The above discussion is limited to cases where w = 1. Gillet-Markowska et al. (2015) proposed applying the Stewart method also to cases where w ≠ 1, which is inappropriate in view of the above example. rSalvador currently is unable to deal with such cases.

Comparing Mutation Rates

Rarely did a researcher perform fluctuation experiments solely to estimate mutation rates. The researcher’s task is often to compare mutation rates across strains or under distinct experimental conditions. Early investigators applied a modified t-test to accomplish this task. This practice can lead to misleading p-values (Zheng 2015a). Recent investigators turned to the Mann-Whitney test (Ford et al. 2013), which also is inappropriate. When applying the Mann-Whitney test to fluctuation assay data, the investigator actually compares the mean numbers of mutations (m) between two experiments, not the two mutation rates that need comparing. To compare mutation rates, the investigator must take into consideration the final cell numbers N_t. It is a common feature that N_t values are consistent within each group (e.g., a strain or a particular experimental condition), but they differ noticeably between groups.

rSalvador provides two approaches to mutation rate comparison. The recommended approach is the likelihood ratio test (LRT) for fluctuation assay data, developed by Zheng (2016a). The second is purely empirical; it tests for equality of two mutation rates by checking whether two 84% C.I.s for the mutation rates overlap. If the two 84% C.I.s overlap, the test is not significant at the 5% level. Empirical evidence supporting the latter approach can be found in Zheng (2015a). Consider now the following two experiments from the study of Krašovec et al. (2014) (data were kindly provided by R. Krašovec).

expt1 = c(0, 2, 0, 3, 0, 15, 21, 0, 1, 0, 0,

4, 4, 1, 3, 2, 1, 0, 0, 0, 1, 0, 9, 0)

expt2 = c(8, 2, 4, 3, 6, 11, 2, 2, 0, 13, 6, 8)

The investigators measured the relative fitness in both experiments and obtained w₁ = 1.47 and w₂ = 1.45. They also measured the two values: N_t,1 = 2.27 × 10⁸ and N_t,2 = 5.15 × 10⁸. A Mann-Whitney test would yield a p-value of 0.01. But this p-value is misleading, as cells in the second population underwent about twice as many cell divisions as cells in the first population (the precise ratio is R = N_t,2 / N_t,1 = 2.29). The LRT approach gives a likelihood ratio statistic of 0.19, corresponding to an approximate p-value of 0.66.

> LRT.MK(expt1,expt2,w1 = 1.47,w2 = 1.45,R = 2.29)

[1] 0.1924875 0.6608543

Not surprisingly, the 84% C.I.s for the two mutation rates overlap.

> confint.MK(expt1,w = 1.47)/2.27e8

[1] 1.798407e–09 5.244954e–09

> confint.MK(expt2,w = 1.45)/5.15e8

[1] 2.155918e–09 6.105526e–09

Both approaches are applicable to cases where partial plating is adopted. Consider mutant count data from the first (strains H37Rv) and third (strain E 729/94) experiments in the study of Werngren and Hoffner (2003). The plating efficiencies in both experiments were 0.4. The final cell density for the H37Rv strain was 2.3 × 10⁸ cells per ml and that for the E 729/94 strain was 1.3 × 10⁸ cells per ml. Because cultures in both experiments were of the same size (5 ml), the ratio of the two N_t values is therefore R = 1.3/2.3. An LRT gives a p-value of 0.62.

> LRT.LD.plating(wh.data[[1]], wh.data[[3]], e1 = 0.4,

e2 = 0.4, R = 1.3/2.3)

[1] 0.2435538 0.6216511

Mazoyer et al. (2017) proposed a new method based on asymptotic normality of estimated mutation rates. The following simulation study compares statistical power of the three methods. Seven groups of 20-culture fluctuation experiments were simulated using the Lea-Coulson mutant distribution. Each group comprises 10,000 experiments. In the baseline group, N_t = 2 × 10⁸ and the mutation rate is 1.0 × 10⁻⁸ . In the other groups, N_t = 1 × 10⁸ and the mutation rates are k × 10⁻⁸ for and 3.0. In all groups, N₀ = 50. Each experiment in the baseline group is compared with the corresponding experiments (i.e., experiments having the same serial number) in all other six groups. Table 1 shows statistical power of the three methods at the 5% significance level. The LRT method and the method of checking overlapping of 84% C.I.s performed almost equally well, but the asymptotic normality method was less powerful.

Accounting for Variability in N_t

There was deep concern that variability in N_t among cultures could interfere with the analysis of fluctuation assay data (Rosche and Foster 2000). However, findings in a recent study (Zheng 2016b) have somewhat lessened worries about that issue. In practice, the experimentalist can use a few cultures to gauge the variability in N_t by calculating the coefficient of variation (CV) for N_t. The new study found that a CV of 0.2 or smaller has a negligible effect on the estimation of mutation rates. More importantly, larger CVs for N_t were rarely encountered in practice. Therefore, if the experimentalist is facing a more critical issue, e.g., partial plating, consideration of variability in N_t should give way to addressing the more critical issue.

Several methods to account for variability in N_t were reviewed in Zheng (2016b). For example, FALCOR refers to a method that requires the experimentalist to measure N_t for all cultures, which often is too laborious to be practical. rSalvador offers a method based on the distribution defined by the PGF in Equation (6). The method assumes that the expected number of mutations in a culture is proportional to the N_t of that culture. It further assumes that, conditional on a culture’s N_t, the number of mutants in that culture obeys a Lea-Coulson distribution. If the CV for N_t in an experiment is C, then, under additional mild assumptions, the unconditional distribution of the number of mutants is a distribution (Zheng 2016b). Here, m₀ is regarded as the overall mean number of mutations per culture, and one can use rSalvador to make inferences about the parameter m₀. If is an estimate of m₀, then an estimate of the mutation rate is , where is the mean of N_t of all cultures in the experiment.

As an example, suppose that the CV for N_t in the Demerec experiment (Demerec 1945) is 0.15. Point and interval estimates of m₀ are obtained as follows.

> newton.B0(demerec.data,cv = 0.15)

[1] 11.09696

> confint.B0(demerec.data,cv = 0.15)

[1] 8.765365 13.665749

Determining Sample Size

The question of “How many cultures are needed” must have been asked by investigators countless times, but it has received little attention in the literature. Investigators often choose a sample size on the basis of a published example in which the sample size was determined intuitively. There were no guidelines for determining sample size until recently (Zheng 2017b).

In planning an experiment, the investigator may form a preliminary idea about the magnitude of m. A large m may reduce the sample size, but it may result in cultures containing more mutants than the investigator is able to count. Zheng (2017b) proposed the following η₅₀₀ index to guide experimental design:

(7)As a rule of thumb, the experimentalist may choose η₅₀₀ < 0.01. Thus, the experimentalist allows no >1% of the cultures to have >500 mutants per culture. As an illustration, assume that the mutation rate to be determined is in the neighborhood of 2 × 10⁻⁷ mutations per cell division. Assume further that the mutants have a relative fitness of 0.75. If the investigator plans to allow the final number of cells per culture to reach 2 × 10⁷, then m would be ∼4.0 according to Equation (1). Therefore, the η₅₀₀ index is

> 1 – sum(prob.MK(m = 4,w = 0.75, n = 500))

[1] 0.001243678

But if the mutants have a relative fitness of 1.2, the η₅₀₀ index will be larger.

> 1 – sum(prob.MK(m = 4,w = 1.2, n = 500))

[1] 0.02221585

If the anticipated value of m meets the η₅₀₀ criterion, the experimentalist determines a required sample size by choosing an appropriate ψ score. The ψ score was introduced by Zheng (2017b) as a convenient yardstick for judging the quality of C.I.s. The ψ score is defined by

(8)Consider a case where has been determined via a fitness assay. The experimentalist further believes that Choosing a ψ score of 0.25 would lead to an estimated sample size of 31.

> samp.size.MK(m = 4,w = 0.75,psi = 0.25)

[1] 31

If an anticipated value of m is large enough to lead to an unacceptable η₅₀₀ index, the experimentalist may resort to partial plating. For example, in the case where and the η₅₀₀ index is

> 1 – sum(prob.MK(m = 50,w = 1,n = 500))

[1] 0.1894214

However, by plating 10% of each culture, the experimentalist reduces η₅₀₀ to

> 1 – sum(prob.LD.plating(m = 50,e = 0.1,n = 500))

[1] 0.01086645

And the corresponding sample size is

> samp.size.LD.plating(m = 50,e = 0.1)

[1] 15