SALVADOR 2.1 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
February 6, 2007

Introduction

SALVADOR  is a package intended as a tool for estimating mutation rates using data from fluctuation experiments that were pioneered by Luria and Delbruck (1941). The package is designed mainly for use in mutation research, but it can also be adopted for educational purposes.  For the most part, this manual adheres to the notation and terminology used in the review article by Zheng (1999), from which 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).  The current version (v 2.1) also contains estimators derived under the Haldane formulation. The Haldane formulation is explained in Sarkar (1991), and more thoroughly examined recently by Zheng (2007). Some of the methods contained in SALVADOR appear outmoded, but they are of historical interest and can play a useful role in education. 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 .

Install["salvadorV2"]

<< salvadorV2.1.m

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

Off[General::spell1]

Simulation

Five basic parameters must be specified before you simulate the number of mutants in a  culture. These five parameters are: 1) mutation rate (per cell per unit time) μ; 2) cellular  birth rate for nonmutants (Subscript[β, 1]); 3) cellular birth rate for mutants (Subscript[β, 2]); 4) initial population size (Subscript[N, 0]); and 5) 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.1, 12, ShowGrowth -> True]

Important Distributions

SALVADOR can compute the probability mass functions (p.m.f.) of four  important 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}], PlotJoined -> 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

SALVADOR concentrates on  estimation of the mean number of mutations denoted by m, because an estimate of the mutation rate (Subscript[μ, β])  can be easily obtained from that of m. The recommended estimation method is the maximum likelihood method. We use the experimental data of Demerec (1945) to illustrate the capabilities  of SALVADOR 2.1.  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 different starting points and see how the two converge to the same m.l.e. Note that if InitialM is set to a negative number, then SALVADOR  uses its own rules to determine an initial guess, which was documented in Zheng (2004).

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

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

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

?newtonM

newtonM[demerec]

Options of newtonM are similar to those of newtonLD.

Options[newtonM]

You can also use the slower mleM to accomplish the same task.

Interval Estimation of Mutation Rates

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

?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 algorithms described in Zheng (2004) to determine reasonable 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 unusually large observations. This type of data may lead to very large estimates of m, say m=800. A typical C compiler cannot process numbers whose magnitude is smaller than (10^(-302)).  If a function invoked to process such data requires a sub-function written in  C , then an underflow error will occur  (e.g., consider the magnitude of the probability of having no mutants, (Subscript[p, 0]) = (e^(-800))).  SALVADOR 2.1 provides three functions, LDprob0,  newtonLD0 and LRIntervalLD0, to deal with extremely large observations. These three functions are equivalent to their counterparts, LDprob, newtonLD and LTIntervalLD, but they rely exclusively on the Mathematica language. Note that these three functions are considerably slower than  their respective counterparts.

pmfLD0[5.7, 15]

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

LRIntervalLD0[demerec,InitialLowerM->3,InitialUpperM->16,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).  But there are several methods of great historical interest; these methods are still popular.

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). To help visualize the Haldane formulation, the user can invoke simuHaldane to simulate the mutational process under Haldane's formulation.

?simuHaldane

simuHaldane[0.02, 10, ShowGrowth -> True]

Now we simulate a fluctuation experiment consisting of 22 test tubes.

data = Table[simuHaldane[0.03, 8], {22}]

A maximum likelihood estimate of the mutation rate μ can be obtained by invoking newtonHaldane, which applies the Newton-Raphson method to calculate the estimate. The algorithm is proposed by Zheng (2007).

?newtonHaldane

newtonHaldane[data, 8, ShowIterations -> True]

An asymptotic 95% confidence interval for the mutation rate can be obtained by invoking CIHaldane. The algorithm is given  in Zheng (2007).

?CIHaldane

CIHaldane[data, 8, ShowIterations -> True]

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

?rossmanMean

rossmanMean[data, 8]

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

?pmfHaldane

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

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

Now we take a look at 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, PlotJoined -> 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 by invoking the function CIHaldane.  Note that this computational process can take several minutes to complete. If only a point estimate is required, the user can invoke the function newtonHaldane.

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

Another 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 InitalM and Inital 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.

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.  to appear (2007).