function mutate1()

mutate1= function( n.gen=50000, n.variables=2, A.init=0.00, mu.init=0.0001, nu.init = 0.00001, title = "Mutation: equilibrating allele frequencies")
{
    # n.samp = 1 means there is no sample population, we're computing the mean.
 
    # mu = mutation rate per generation A to a.
    # nu = mutation rate per generation a to A.
 
    # p = allele frequency of A(for all individuals and for all generations)
 
    # Parameters for environment
    mu = mu.init     # mutation rate per generation A to a.
    nu = nu.init     # mutation rate per generation a to A.
 
    # Initialize matrix
    k = matrix(-999, nr=n.variables, nc=n.gen) #matrix to store values of p
 
    for (i in 1:n.variables)
    {
        k[1,1] = A.init    #  initial freq of A
        k[2,1] = 1- A.init    #  initial freq of a
     
        for (j in 2:(n.gen) )   # n.gen generations
        {
            k[1, j] = (1-mu) *  k[1, j-1] + nu * (1 - k[1, j-1]) # freq of A
            k[2, j] = (1-mu) *  k[2, j-1] + nu * (1 - k[2, j-1]) # freq of a
        }
    }
    #######
    #
    # Plot allele frequency vs. generation
    graphics.off()
    plot(1,1, xlim=c(1,n.gen), ylim=c(0,1), ty='n', xlab='Generation', ylab='Allele Frequency' )
    mtext(line=0, cex=0.7, paste('gen=', n.gen,  'A.init=', A.init,  'a.init=', 1-A.init, 'mu.init=', mu.init, 'nu.init=', nu.init ))
    mtext(line=1, cex=1.0, title)
    #
    lines(1:n.gen, k[1, ], lwd=2, col=1, lty=1)  # allele A, black
    lines(1:n.gen, k[2, ], lwd=2, col=4, lty=2)  # allele a, blue
    legend(45000, 0.9, c("A","a"), cex=0.8,
    col=c(1,4), lty=c(1,1), lwd=2)
 
} ###### END