SALVADOR 2.3 User Manual

Qi Zheng
Department of Epidemiology and Biostatistics
School of Rural Public Health
Texas A&M Health Science Center
College Station, Texas 77843
qzheng@srph.tamhsc.edu
September 10, 2008

Introduction

SALVADOR  is a package intended as a tool for estimating mutation rates using data generated under the Luria-Delbruck experimental protocol (Luria and Delbruck 1941).  The package was originally designed  for use in mutation research, but later it also appears to be a useful tool for educational purposes.  For the most part, this manual adheres to the notation and terminology used in the review article by Zheng (1999), where numerous important references can be found. The focus of SALVADOR is on maximum likelihood estimators derived under the Lea-Coulson formulation (Lea and Coulson 1949).  This formulation has been extensively studied; among major contributions are those made by Koch (1982), Mandelbrot (1974), Crump and Hoel (1974), Stewart et al. (1990) and Ma et al. (1992).  SALVADOR also contains estimators derived under the Haldane formulation and the Bartlett formulation. The Haldane formulation was explained in Sarkar (1991), and later was more thoroughly examined  by Zheng (2007). The Bartlett formulation was proposed in Bartlett's discussion of the Armitage (1952) paper and was further elaborated in Bartlett's textbook on stochastic processes (Bartlett 1978). The algorithms implemented in SALVADOR for the Bartlett formulation are those due to Zheng (2008). New to the current version (v2.3) is a maximum likelihood estimator that accounts for plating efficiency under the Lea-Coulson formulation (Zheng 2008a). Some of the methods contained in SALVADOR appear outmoded, but they are of historical interest and some are employed by SALVADOR to generate excellent initial estimates for the maximum likelihood estimators. Methods for simulating experimental data are provided mainly for educational purposes. For reasons of computational efficiency, SALVADOR was written in a mixture of the Mathematica and the C language. Consequently,  two programs must be install from within Mathematica .

SetDirectory["c:\\sal2-3"]

Install["salvadorV2-3.sal"]

<<salvadorV2-3.m

(To avoid distraction by a plethora of warning messages against use of similar variable names, e.g., var1 and var2, we suppress such warning messages.)

Off[General::spell1]

Simulation

To simulate under  the Lea-Coulson formulation one needs to specify five basic parameters: 1) a mutation rate (per cell per unit time) μ; 2) a cellular  birth rate for nonmutants (Subscript[β, 1]); 3) a cellular birth rate for mutants (Subscript[β, 2]); 4) an initial population size (Subscript[N, 0]); and 5) a final population size (Subscript[N, T]). It is strongly recommended that, before running a simulation,  you check the parameters of interest by computing the mean and variance of the number of mutants, for unrealistic parameters can lead to very time-consuming computational tasks and useless results. We now use the following five parameters to simulate a fluctuation experiment .

mu=10^-8;b1=1.2;b2=1.1;N0=100;Nt=2.4*10^8;

?MeanOfMutants

MeanOfMutants[mu,b1,b2,N0,Nt]

?VarOfMutants

VarOfMutants[mu,b1,b2,N0,Nt]

?simuMutants

Now we simulate the numbers of mutants in 20 sister cultures.

Table[simuMutants[mu,b1,b2,N0,Nt],{20}]

The above simulation is based on a  continuous-time formulation, under which the Luria-Delbruck distribution was derived. There is an obvious discrete-time counterpart, which was well explained by Angerer (2001).

?simuAngerer

simuAngerer[0.3,12,ShowGrowth->True]

Important Distributions

SALVADOR can compute the probability mass functions (p.m.f.) of commonly used distributions in fluctuation analysis: the Poisson distribution, the LD or Luria-Delbruck  distribution, the M or Mandelbrot (1974) distribution, and the convolution of a Poisson and an LD distribution. The last distribution was  popularized in fluctuation analysis by Cairns et al.  (1988). Using pmf2cdf, you can transform a list of p.m.f. to a list of c.d.f. (cumulative distribution function).  As a demonstration we  first compute (Subscript[p, 1]), ..., (Subscript[p, 10]) of a Poisson distribution having mean 10.5.

?pmfPoisson

pmfPoisson[10.5,10]

Now we plot the c.d.f. of a Luria-Delbruck distribution, so its upper tail behavior can be more easily examined.

?pmfLD

LDprob=pmfLD[5.7,150]; LDprob[[Range[10]]]

LDcdf=pmf2cdf[LDprob]; LDcdf[[Range[10]]]

ListPlot[Transpose[{Range[0,150],LDcdf}],Joined->True,Frame->True]

Now we compare an approximate form and the exact form of an M distribution.

?pmfM

pmfM[10.5,0.8,10]

If we assume ϕ=1, we  can also use an algorithm devised by Koch (1982)  (based on an algorithm of Lea and Coulson (1949)) to compute the p.m.f. of an M distribution.

?pmfMKoch

pmfMKoch[10.5,0.8,10]

To use the exact form of the M distribution, the parameter ϕ<1 must be specified.

pmfM[10.5,0.8,0.9999,10]

Finally let us take a look at the convolution of a Poisson distribution and an LD distribution .

?pmfCairns

pmfCairns[2.5,1,8.5,10]

Point Estimation of Mutation Rates

Under the Lea-Coulson formulation attention is focused on the mean number of mutations m, from which the mutation rate (Subscript[μ, β])  can be easily obtained. The recommended estimation method is the maximum likelihood method. We use the experimental data of Demerec (1945) to illustrate the capabilities of SALVADOR.  The following are the known experimental parameters of Demerec's experiment.

N0=90;Nt=1.9*10^8;phi=1-N0/Nt;

The experimental data of Demerec's (1945) experiment are given below.

demerec={33,18, 839,47,13,126,48,80,9,71,196,66,28,17,27,37,126,33,12,44,28,67,730,168,44,50,583,23,17,24};

We can fit an LD(m,1) distribution to the data by calling  newtonLD, which uses the Newton-Raphson methods to compute the maximum likelihood estimate (m.l.e.) of m.

?newtonLD

newtonLD[demerec]

Sometimes it is necessary to choose an initial guess, to set convergence criteria or to specify the maximum number of iterations. Thus it is helpful to familiarize yourself with the options of newtonLD.

Options[newtonLD]

For example, we can choose two distinct starting points and see how the two converge to the same m.l.e. Note that if InitialM is set to a negative number,  SALVADOR  uses its own rules to determine an initial guess (Zheng  2005).

newtonLD[demerec,InitialM->3,ShowIterations->True]

newtonLD[demerec,InitialM->15,ShowIterations->True]

If you believe that mutants grow at a rate different from that for nonmutants, you may try to fit a M(m,r) model to your data. Note that r in the current version of SALVADOR is the reciprocal of ρ in SALVADOR 1.0.

?newtonM

newtonM[demerec]

Options of newtonM are similar to those of newtonLD.

Options[newtonM]

Interval Estimation of Mutation Rates

SALVADOR relies on large-sample theory to construct confidence intervals (CIs).   For the equal growth model under the Lea-Coulson formulation,  the function for constructing a CI for m is LRIntervalLD.

?LRIntervalLD

A 95% asymptotic confidence interval for m is obtained by the command

LRIntervalLD[demerec]

And a 90% confidence interval would be constructed by using the Alpha option. Note that if a (1-α)100% CI is desired, then the options Alpha should be set to α, not α/2.

LRIntervalLD[demerec,Alpha->0.1]

In addition to Alpha, there are a few other options.

Options[LRIntervalLD]

InitialM sets an initial guess for computing the maximum likelihood estimate of m, which is needed for computing a CI for m. The two options, InitialLowerM and initialUpperM, are starting points for computing the lower and upper boundary points of the desired CI for m. When these  options are set to a negative number, SALVADOR adopts its own algorithms as described in Zheng (2004) to determine  starting values.  For example, if we choose 3 as an initial guess for the lower boundary point, and 16 as an initial guess for the upper boundary point, the iteration process will proceed as follows.

LRIntervalLD[demerec,InitialLowerM->3,InitialUpperM->16,ShowIterations->True]

You can also construct a Wald type CI for m by first computing expected Fisher information as suggested by Zheng (2002).

?fisherInfoLD

mhat=newtonLD[demerec]

info=fisherInfoLD[mhat,50000]

{mhat-1.96 Sqrt[1/(30 info)], mhat+1.96 Sqrt[1/(30 info)]}

This Wald type CI for m is very similar to the likelihood ratio based CI we just found.

When fitting a differential growth model to data, use the function profileIntervalM to construct a profile likelihood based CI for m, and use the function profileIntervalR to construct a CI for r.

?profileIntervalM

profileIntervalM[demerec]

To get a Wald type CI for m, we have to compute the expected Fisher information matrix approximately. Since we know little about the converging speed of computing this information matrix,  it is strongly recommended that a Wald type CI be  corroborated by a likelihood ratio based CI.

?fisherInfoMatrix

{mhatm,rhatm}=newtonM[demerec]

infoMat=fisherInfoMatrix[mhatm,rhatm,5000]

{{var11,var12},{var21,var22}}=Inverse[infoMat]

{mhatm-1.96 Sqrt[var11/30],mhatm+1.96 Sqrt[var11/30]}

Again the above Wald type CI is similar to the likelihood ratio based CI we just found.

To get a 99% CI for M, set Alpha to 0.01, not 0.005.

profileIntervalM[demerec,Alpha->0.01]

To see if the differential growth model is necessary, we can construct a 95% CI for r.

profileIntervalR[demerec]

Because unity is  contained in the above CI, the equal growth model appears sufficient.

Sometimes it is necessary for the user to supply initial guesses  when constructing a likelihood based CI.

Options[profileIntervalM]

profileIntervalR[demerec,InitialLowerR->0.5,InitialUpperR->1.5,ShowIterations->True]

Functions for extremely large observations

Occasionally,  the user may encounter data having extremely large observations. This type of data may force SALVADOR to process numbers of exceedingly small magnitudes. As a  typical C compiler cannot manipulate numbers of magnitudes smaller than (10^(-302)), underflow may occur. For example, if m=800, (Subscript[p, 0])=(e^(-800)) will be treated as zero,  if a function written in C is involved in its computation. This can invalidate all the subsequent computations. However, the Mathematica language can manipulate numbers much smaller (the user is referred to the Mathematica system variable $MinNumber). For this reason, SALVADOR offers several functions whose names end with the numeral 0. Such a "zero function"  uses the same algorithm as its counterpart does, but is implemented solely in the Mathematica language. For instance, newtonLD0 is a counterpart of newtonLD.

pmfLD0[5.7,15]

newtonLD0[demerec,InitialM->3,ShowIterations->True]

Methods of Historical Interest

The maximum likelihood method is widely regarded as the method of choice for estimating mutation rates; see, for example, the simulation studies of  Rosche and Foster (2000).  However, there are several methods of great historical interest; these relatively simple methods are still popular and some are actually adopted by SALVADOR to yield initial guesses for the maximum likelihood estimators.

Note that the (Subscript[P, 0]) method is not always applicable.

?MethodOfP0

MethodOfP0[demerec]

Note that Luria and Delbruck's method of the mean has several variants, one of which was suggested by Armitage.

?MethodOfMeans

Options[MethodOfMeans]

MethodOfMeans[demerec,ShowIterations->True]

To adopt Armitage's (1952) correction, chose β=3.46.

MethodOfMeans[demerec,ShowIterations->True,Correction->3.46]

Replace the sample mean by the sample median in Luria and Delbruck's method of the mean gives Drake's (1991) estimator.

?drakeMedian

Options[drakeMedian]

drakeMedian[demerec,ShowIterations->True]

This is different from the method of the median proposed by Lea and Coulson (1949).

?MethodOfMedians

MethodOfMedians[demerec, ShowIterations->True]

med=Median[demerec]

MethodOfMedians[med, ShowIterations->True]

Lea and Coulson (1949) also proposed a method that bears some resemblance to the maximum likelihood method.

?semiMLE

semiMLE[demerec]

It is noteworthy that Jones at al. (1994)   proposed a method based on the median that does not require an initial guess.

?jonesMedian

jonesMedian[demerec]

Haldane's Formulation

Under Haldane's formulation cell growth is synchronous. If a cell population is initiated by a single nonmutant cell, then the total number of cells (the sum of all mutants and all nonmutants) in the g-th generation is (2^(g )). Each division of a nonmutant cell can give rise to a mutant daughter cell with probability μ. The total number of mutants in the g-th generation is said to obey an Haldane distribution denoted by H(μ,g). Relevant algorithms for the Haldane formulation are discussed in Zheng (2007). 

To help visualize the Haldane formulation, the user can invoke simuHaldane to simulate the mutational process under Haldane's formulation.The function simuHaldane assumes that there is just one wild-type  cell at the beginning, that is, (Subscript[N, 0])=1.

?simuHaldane

simuHaldane[8,0.02,ShowGrowth->True]

Now we simulate a fluctuation experiment by assuming μ=0.002 and g=8.

datHal0=Table[simuHaldane[8,0.005],{75}]

Note that simuHaldane assumes  (Subscript[N, 0])=1. However, this is not an obstacle. For example, if (Subscript[N, 0]) =3 is the desired initial number of nonmutant cells, we can process the above data  as follows.

datHal=(Plus@@#)&/@ Partition[datHal0,3]

A maximum likelihood estimate of the mutation rate μ can be obtained by invoking newtonHaldane, which applies the Newton-Raphson method to calculate a estimate.

?newtonHaldane

newtonHaldane[datHal,8,ShowIterations->True, InitialCells->3]

newtonHaldane0[datHal,8,ShowIterations->True, InitialCells->3]

An asymptotic 95% confidence interval for the mutation rate can be obtained by invoking CIHaldane.

?CIHaldane

CIHaldane[datHal,8,ShowIterations->True, InitialCells->3]

Historically, the mutation rate μ was  estimated by the method of the mean. The method implemented in SALVADOR was that first appeared in Rossman et al. (1995).

?rossman

rossman[datHal,8,3]

Probability mass functions of the Haldane distribution can be computed as follows.

?pmfHaldane

ppH=pmfHaldane[10,0.01,200,1];

g1=ListPlot[ppH,Joined->True,PlotRange->All]

Now we find the corresponding Luria-Delbruck distribution and make a comparison.

phi=1-1/1024.

?haldaneMutations

mm=haldaneMutations[0.01,10]

Note in passing that the relevancy of Kendall's Assumptions in the present context is explained in Zheng (1999). The Galton-Watson formulation is explained in Kimmel and Axelrod (1994). Under either formulation the expected number of mutations would be computed as follows:

haldaneMutations[0.02,10]

Now we can compute the corresponding Luria-Delbruck distribution.

ppL=pmfLD0[10.23,phi,200];

g2=ListPlot[ppL,PlotRange->All,Joined->True,PlotStyle->Dashing[{0.01}]]

Now we can compare the two distributions.

Show[g1,g2,Frame->True, FrameLabel->{"Number of mutants","Probability"}]

As an illustration, we apply the Haldane model to analyze the data given by Demerec (1945).  To begin with, we estimate g as follows.

N0=90;Nt=1.9*10^8;
Log[Nt/N0]/Log[2]

The above result suggests that a reasonable estimate of g is 21. Under this assumption we can now proceed to obtain both point and interval estimates of the mutation rate.

newtonHaldane[demerec,21,InitialCells->N0,ShowIterations->True]

CIHaldane[demerec,21,InitialCells->N0,ShowIterations->True]

The Bartlett Formulation

This  formulation of the Luria-Delbruck model was proposed by Bartlett in 1952.  Under this formulation, nonmutant cells also grow stochastically, governed by a Yule birth process. For this reason, Bartlett called this formulation  a fully stochastic model, in contrast to the Lea-Coulson formulation under which nonmutant cells grow deterministically. The mutant distribution under this formulation is  called the Bartlett distribution, which can be computed as follows.

?pmfBartlett

pmfBartlett[0.00001,0.999,12]

The following function is for simulation under the Bartlett formulation.

?simuBartlett

simuBartlett[0.1,0.02,1,14.5,ShowGrowth->True]

Now we simulate a fluctuation experiment consisting of n=20 tubes.

datBart=Table[simuBartlett[1.5,0.002,5,4.5],{20}]

The mutation rate is known in advance.

truMu=0.002/1.5

We  use this data set as an illustrative example. First, let us extract the numbers of mutants from the simulated data.

simulated=#[[2]]&/@ datBart

Clearly, the parameter ϕ is

phi=1-Exp[-1.5*4.5]

Now we apply the maximum likelihood procedure (Zheng 2008) to estimate the mutation  rate.

?newtonBartlett

newtonBartlett[simulated,5,phi]

We can also construct an asymptotic confidence interval for the mutation rate.

?CIBartlett

CIBartlett[simulated,5,phi]

Now it is straightforward to estimate the mutation rate for the Demerec experiment.

demerec={33,18, 839,47,13,126,48,80,9,71,196,66,28,17,27,37,126,33,12,44,28,67,730,168,44,50,583,23,17,24};

phi=1-90/(1.9 10^8)

newtonBartlett[demerec,90,phi,ShowIterations->True]

CIBartlett[demerec,90,phi,ShowIterations->True]

Plating efficiency

Plating efficiency, denoted by ϵ,  refers to the probability that a mutant cell living in the tube actually forms a visible colony on a solid culture in a dish after being plated (transferred). Imperfect plating efficiency (ϵ<1) can be caused by partial plating or by the cell's imperfect ability to survive the plating process. An example of partial plating is Luria and Delbruck's experiment #16 (Luria and Delbruck 1943) . In this experiment a 0.08-ml portion of each of the 12 0.2-ml cultures was plated; thus the plating efficiency is ϵ=0.4.

e16=0.08/0.2

The original data of the experiment are as follows.

expt16={1,0,3,0,0,5,0,5,0,6,107,0,0,0,1,0,0,64,0,35}

The cell population size prior to plating  is (Subscript[N, T])=5.6×108.

n16=5.6  10^8

The only existing method to account for plating efficiency is to assume that the number of mutants in a tube prior to plating obeys the LD(m,1) distribution and that each mutant cell can survive the plating process with probability ϵ. The resulting mutant distribution can be calculated by invoking the function pmfLDPlating.

?pmfLDPlating

Assume that m = 37.2 and ϵ = 0.1.

dist=pmfLDPlating[37.2,0.1,150];dist[[Range[6]]]

ListPlot[Transpose[{Range[0,150],dist}],Joined->True,Frame->True,PlotRange->All]

The first method to account for plating efficiency is to adjust the sample mean and then use Luria and Delbruck's method of means. This is the way Luria and Delbruck actually analyzed their data in 1943. Because adjusting the sample mean is the same as adjusting the raw data, one can proceed as follows.

?MethodOfMeans

Options[MethodOfMeans]

m=MethodOfMeans[expt16/e16]

So the mutation rate is

m/n16

Note that the P0 method is different from the well-known one designed for ϵ=1.

?P0Plating

m=P0Plating[expt16,0.4]

Hence, the mutation rate is

m/n16

The maximum likelihood method (Zheng 2008a) can be applied  as follows.

?newtonLDPlating

m=newtonLDPlating[expt16,0.4,ShowIterations->True]

m/n16

An asymptotic 95% confidence intervals for m can be constructed as follows.

?CILDPlating

Options[CILDPlating]

ci=CILDPlating[expt16,0.4,ShowIterations->True]

Therefore, an CI for the mutation rate is

ci/n16

Note that SALVADOR also offers  Jones median method and  Lea and Coulson's median method as modified by Angerer.

?JonesMedianPlating

JonesMedianPlating[expt16,0.4]

LCAMedian[expt16,0.4]

A limiting distribution: a special topic

The Bartlett mutant distribution gives us an interesting two-parameter distribution whose probability generating function is of the form
                      G (z; A, k) = [ (((1 - A((z^(-1)) - 1) log(1 - z)])^(-k)),
 where A and k are positive real numbers.   This is an infinitely divisible distribution with divergent moments (Zheng 2008b).  The distribution is denoted by (B^0)(A,k).
 The probability mass function  can be calculated as follows.

?pmfBartzero

prob=pmfBartzero[1.3,5.8,100];prob[[Range[5]]]

ListPlot[prob,Joined->True]

The following sample was generated from the (B^0)(1.3,5.8) distribution.

simuSample={11,8,14,4,34,37,3,30,50,17,69,4,180,968,28,28,7,25,17,15,41,29,279,0,21,19,38,159,66,50};

Assuming that k=5.8 is known, we can find the maximum likelihood estimate of A as follows.

newtonBartzero[simuSample,5.8,ShowIterations->True,MaxIterations->33]

Similarly, an asymptotic 95% for the parameter A can be constructed.

?CIBartzero

CIBartzero[simuSample,5.8,ShowIterations->True]

There is also a P0 method for the (B^0)(A,k) distribution.

?p0Bartzero

p0Bartzero[simuSample,5.8]

A Concluding Example

Boe et al. (1994) meticulously conducted a large scale fluctuation experiment, producing the following 1104 observations.

boe={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,3,4,4,4,4,6,8,14,15,17,27,29,39,39,66,74,265,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,4,4,4,5,5,5,7,8,12,16,18,19,30,42,132,152,513,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,2,2,2,3,4,5,5,6,7,8,12,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,2,2,2,3,4,4,5,6,6,13,18,19,26,36,68,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,4,4,5,6,6,7,12,16,24,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,2,2,2,3,4,4,5,5,6,10,11,12,14,24,27,29,30,73,482,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,3,3,5,11,16,30,31,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,2,2,2,2,2,2,2,2,4,5,5,6,6,6,7,12,12,12,16,513,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,5,7,8,11,11,18,52,73,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,4,4,4,5,7,7,7,151,513,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,3,3,3,4,4,5,6,11,13,16,40,258,320,513,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,4,4,4,11,11,17,107,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,4,4,6,6,16,18,20,27,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,2,2,3,3,5,8,8,9,21,21,26,37,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,3,3,3,3,4,5,5,6,6,7,7,8,10,12,14,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,2,3,3,4,5,5,6,13,16,19,26,41,105,116,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,3,6,6,6,11,13,19,21,34,140,146,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,4,9,10,10,11,13,19,20,23,59,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,2,2,3,3,3,4,4,4,4,4,5,5,8,8,9,10,21,26,192,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2,3,3,4,4,6,9,9,10,14,25,26,29,57,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,5,6,7,8,10,14,16,28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,3,4,4,6,6,7,7,18,32,35,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,4,4,5,5,8,10,11,12,19,21};

Assuming equal growth rate we obtain the maximum likelihood estimate of m  as follows.

newtonLD[boe]

The fitting is less than perfect, but still an LD distribution catches the most prominent features of the data.

plotFit[boe,logPlot->True]

If we use the differential growth model, we can improve the fitting significantly.

newtonM[boe]

plotFit[boe,DifferentialGrowth->True,logPlot->True]

By setting options appropriately, one can use plotFit to produce informative graphics.

Options[plotFit]

Note that when you set the option logPlot to True, the horizontal axis is log(x+1) in stead of log(x),  to avoid log(0). The meaning of the two options InitialM and Initial R is the same as in newtonM and newtonLD. As a illustrative example, we now use graphics to compare the two models fitted to the Boe el al. (1994) data.

graphic1=plotFit[boe,fitStyle->{RGBColor[0,0,0],Dashing[{0.015}]},dataStyle->{RGBColor[0,0,0],PointSize[0.012]},logPlot->True]

graphic2=plotFit[boe, DifferentialGrowth->True,fitStyle->RGBColor[0,0,0],dataStyle->{RGBColor[1,1,1]},logPlot->True]

Show[graphic2,graphic1,FrameLabel->{"log of k+1", "prob. of no more than k mutants"}]

References

W. P. Angerer, An explicit representation of the Luria-Delbruck distribution.  J. Math. Biol. 42 (2001) 145-174.

 P.  Armitage, The statistical theory of bacterial populations subject to mutation, J. R. Statist. Soc. B 14 (1952) 1-33.

M.S. Bartlett, An Introduction to Stochastic Processes With Special Reference to Methods and Applications, third edition, Cambridge University Press, 1978.

L. Boe, T. Tolker-Nielsen, K.-M. Eegholm, H. Spliid, A. Vrang, Fluctuation analysis of mutations to nalidixic acid resistance in Escherichia coli. J. Bacteriol. 176 (1994) 2781-2787.

J. Cairns,  J. Overbaugh, S. Miller, The origin of mutants. Nature 335 (1988) 142-145.

K.S. Crump, D.G. Hoel, Mathematical models for estimating mutation rates in cell populations. Biometrika  61 (1974) 237-252.

M. Demerec, Production of staphylococcus strains resistant to various concentrations of penicillin,  Proc.  Nat.  Acad.   Sci. USA 31 (1945) 16-24.

J.W. Drake, A constant rate of spontaneous mutation in DNA-based microbes, Proc. Nat. Acad. Sci. USA 88 (1991) 7160-7164.

 M.E. Jones, S. M. Thomas and A. Rogers, Luria-Delbruck fluctuation experiments: design and analysis. Genetics  136 (1994) 1209-1216.

M. Kimmel, D.E. Axelrod, Fluctuation test for two-stage mutations: application to gene amplification, Mutation Res. 306 (1994) 45-60.

A.L. Koch, Mutation and growth rates from Luria-Delbruck fluctuation tests. Mutation Res. 95 (1982) 129-143.

E.A. Lea, C.A. Coulson, The distribution of the numbers of mutants in bacterial populations.  J. Genetics 49 (1949) 264-285.

 S.E. Luria, M. Delbruck, Mutations of bacteria from virus sensitivity to virus resistance. Genetics 28 (1943) 491-511.

W.T. Ma, G. vH. Sandri, S. Sarkar, Analysis of the Luria and Delbruck distribution using discrete convolution powers.  J. Appl. Prob. 29 (1992) 255-267.

B. Mandelbrot, A population birth-and-mutation  process, I: Explicit distributions for the number of mutants in an old culture of bacteria. J. Appl. Prob. 11 (1974) 437-444.

W.A. Rosche, P. L. Foster, Determining mutation rates in bacterial populations. Methods,  20 (2000) 4-17.

 T.G. Rossman, E.I. Goncharova, A. Nadas, Modeling and measurement of the spontaneous mutation rate in mammalian cells, Mutation Res. 328 (1995) 21-30.

S. Sarkar, Haldane's solution of the Luria-Delbruck distribution, Genetics 127 (1991) 257-261.

F.M. Stewart, D.M. Gordon, B.R. Levin, Fluctuation analysis: the probability distribution of the number of mutants under different conditions, Genetics 124 (1990) 175-185.

 Q. Zheng, Progress of a half century in the study of the Luria-Delbruck  distribution, Math. Biosci. 162 (1999) 1-32.

Q. Zheng, Statistical and algorithmic methods for fluctuation analysis with SALVADOR as an implementation, Math. Biosci. 176 (2002) 237-252. 

Q. Zheng, New algorithms for Luria-Delbruck fluctuation analysis,  Math. Biosci. 196 (2005) 198-214.

Q. Zheng,  On Haldane's formulation of Luria and Delbruck's mutation model, Math.  Biosci.  209 (2007) 500-513.

Q. Zheng, On Bartlett's formulation of the Luria-Delbruck mutation model, Math. Biosci. 215 (2008) 48-54.

Q. Zheng,  A note on plating efficiency in fluctuation experiments, to appear in Math. Biosci. (2008a).

Q. Zheng, On a limiting distribution derived from the Bartlett formulation of the Luria-Delbruck mutation model, preprint, (2008b).