SALVADOR 1.0 Manual

Qi Zheng
National Center for Toxicological Research
3900 NCTR Road, Jefferson, Arkansas 72079
November, 2001

Introduction

SALVADOR  is a package intended as a tool for exploring on the computer various aspects of fluctuation analysis. 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["salvador"]

<< salvador.m

For notation, nomenclature and other background information, the user can consult the 2 references listed at the end of this mannual and numerous additional references contained therein. The sole purpose of this mannual is to demonstrate the various functionalities of SALVADOR.

Simulation

Five basic parameters must be specified before you simulate the number of mutants in a sister 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.999, 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 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 are as follows.

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 using the maximum likelihood method.

?mleLD

mleLD[demerec]

It is also possible to fit an M(m,ρ)  distribution to this data set. Because two parameters are involve in the model, it takes considerably more time to compute the maximum likelihood estimates for both m and ρ.

?mleM

mleM[demerec]

Fisher's Information

Fisher's information for m is needed in constructing confidence intervals for mutation rates. Two methods of computing  Fisher's information are included in SALVADOR. For mathematical details about these two methods, the user is referred to Zheng (2001).

est = mleLD[demerec]

?FisherInfoLD

We can choose two different values of k to help determine how many terms are needed to get a sufficiently accurate approximation of I(m).

{FisherInfoLD[est, 1, 1000], FisherInfoLD[est, 1, 5000]}

If we assume ρ≠1, then we proceed differently.

{paraM, paraRho} = mleM[demerec]

?FisherInfoM

FisherInfoM[paraM, paraRho, 2000]

?FisherInfoRho

FisherInfoRho[paraM, paraRho, 2000]

The method based on the original idea of Lea and Coulson takes more time when a large number of terms are retained in the series. You should not see much difference if only a small number of terms are kept. Of course the two methods give identical results.

?FisherInfoMLea

FisherInfoMLea[paraM, paraRho, 2000]

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]

How to Get Results Quickly

No computer package can replace critical thinking. However, SALVADOR offers 2 functions that enable the user to get a feeling of his experimental data without thinking about mathematical details. Critical thinking is strongly recommended when you look at your preliminary results obtained in such a completely automated manner.  The first function allows the user to get a confidence interval for the mutation rate (Subscript[μ, β]) automatically.

?CIMutationRate

We now find a 95% confidence interval for (Subscript[μ, β]) using Demerec's experimental data.

CIMutationRate[demerec, 0.95, Nt, 2000]

\(We\ can\ also\ get\ a\ 90  %\ confidence\ interval\ for\ \(\(μ\_β\)\(.\)\)\)

CIMutationRate[demerec, 0.90, Nt, 2000]

If you choose too small a k, SALVADOR will complain.

CIMutationRate[demerec, 0.95, Nt, 50]

The second function enables you to visualize your data along with a fitted distribution in an automated manner.

?plotFit

Options[plotFit]

Note that the logPlot option should be set to false if some sister cultures contain no mutants, because log(0) is undefined.

g1 = 
      plotFit[demerec, 
        fitStyle -> {RGBColor[0, 0, 1], Thickness[0.006], 
            Dashing[{0.01, 0.02}]}, 
        dataStyle -> {RGBColor[1, 0, 0], PointSize[0.02]}, 
        logPlot -> True]

Here is another way to view the same data set.

g2 = 
      plotFit[demerec, logPlot -> False, 
        DifferentialGrowth -> True, 
        dataStyle -> {RGBColor[0, 0, 1], PointSize[0.01]}, 
        fitStyle -> {RGBColor[1, 0, 0]}]

Performance Testing

The following parameters and random seed are deliberately chosen to generate 60 sister cultures, among which one sister culture contains an unusually large number of mutants (3544). This simulated data will allow you to test the speed of SALVADOR. On the author's PC with 256 MB RAM and a CPU of 550 MHz, it took about 21 seconds to obtain a maximum likelihood estimate of m,  when ρ=1 is assumed. It takes about 39 seconds to find the maximum likelihood estimates of both m and ρ, when an M(m,ρ,1) model is assumed.

mu = 10^-8; Nt = 10^8; N0 = 10; b1 =  .08; b2 =  .072;  
    SeedRandom[5878]

Here are the true parameters.

{mu/b1 (Nt - N0), b2/b1}

simudata = Table[simuMutants[mu, b1, b2, N0, Nt], {60}]

t0 = AbsoluteTime[]; 
    est = mleM[simudata]; {AbsoluteTime[] - t0, est}

References

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

Q. Zheng (2001), Statistical and algorithmic methods for fluctuation analysis with SALVADOR as an implementation, submitted manuscript.