SALVADOR 2.0 Manual

Qi Zheng
Department of Epidemiology and Biostatistics
School of Rural Public Health,  Bryan, Texas 77802
qzheng@srph.tamhsc.edu
February, 2005

Introduction

SALVADOR  is a package intended as a tool for exploring  various quantitative aspects of fluctuation experiments. For reasons of computational efficiency, the package was written in a mixture of the Mathematica and the C language. Therefore, before using SALVADOR, two programs must be install from within Mathematica.

Install["salvadorV2"]

<< salvadorV2.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]

For notation, nomenclature and other background information, the user can consult Zheng(1999,2002,2005) and numerous additional references contained therein. The sole purpose of this manual is to demonstrate the various functionalities of SALVADOR 2.0.

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}]

Important Distributions

SALVADOR can compute the probability mass function (p.m.f.) of four most important distributions in fluctuation analysis: the Poisson distribution, the LD or Luria-Delbruck  distribution, the M or Mandelbrot distribution, and the convolution of a Poisson and an LD distribution. The last distribution was  popularized in fluctuation analysis by Cairns et al. in 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 let's look at both 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  (based on an algorithm of Lea and Coulson) 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]

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.0.  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 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 critiria or to specify the maximum number of iterations. Thus it is helpful to familiarize yourselve 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 2.0 uses its own rules to deterimine 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.0 is the reciporacol 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 melM to accomplish the same task.

Interval Estimation of Mutation Rates

The most important new feature of SALVADOR 2.0 is its capability of computing (profile) likelihood ratio-based confidence intervals.  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 Alpah should be set to α, not α/2, because critical values are determined by  2 chi-squared distribution.

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 intialUpperM, are statrting 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 reasosnalbe 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 boudary 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.

?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 simular 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 informtion 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 unit is not contained in the above CI, the equal growth model may be sufficient.

Sometimes you may want to use some of the options when constructing a likelihood based CI.

Options[profileIntervalM]

For example, occasionally you may have to try a few starting points in constructing a confidence interval.

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

Methods of Historical Interest

The maximum likelihood method is obviously the method of choice for estimating mutation rates. But there are several methods of great historical interest; these methods are still popular.

?MethodOfP0

MethodOfP0[demerec]

?MethodOfMeans

MethodOfMeans[demerec]

?MethodOfMedians

MethodOfMedians[demerec]

?semiMLE

semiMLE[demerec]

It is noteworthy that Jones has proposed a method based on the median that does not require an initial guess.

JonesMedian[demerec]

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 have 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 carefully setting options, one can use plotFit to draw publication quality graphics.

Options[plotFit]

Note that when you set 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

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

M. Demerec (1945), production of staphylococcus strains resistant to various concentrations of penicillin, proceedings of the National Academic of Science, 31: 16-24.

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

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

Q. Zheng (2005), New algorithms for Luria-Delbruck fluctuation analysis, Mathematical Biosciences , to appear.