Frequency of mutants in a very large number of observations

Question

During cellular division (generation of two cells from one) bacteria experience mutations, where a newly produced cell harbors altered genetic code. To measure the rate of this happening Max Delbrück and Salvador Luria have designed an experiment from which it is possible to measure such mutation rate.

Now, since mutation is a random event and it can arise any time during the growth of bacterial population the frequency of mutants can vary greatly between different populations when the population reaches its peak, because for example it exhausts the resources. Estimating the mutation rate solely on the frequency is thus not a sufficiently accurate measure.

Now my question is the following: let's suppose we start with a single cell that has zero mutations, and at each time step each cell divides in two. After each division, one of the daughter cells can gain a mutation with probability $\mu$, after gaining it, it doesn't lose it and as well each offspring has it. There is no death in such a population and each cell divides at the same rate (up until the exhaustion of resources, when it drops to zero).

If we had the possibility to look at the frequency of mutants in a very large number of different bacterial populations (approaching infinity), would the average frequency be dependent on $\mu$ and what would it be? If not, what would it be in that case?

A bit of background for better clarity: https://en.wikipedia.org/wiki/Luria%E2%80%93Delbr%C3%BCck_experiment

Answer 1

Since I was unable to analytically answer the question, I performed several simulations.

As it seems, by keeping the mutation rate constant, the frequency of mutants increases with increasing population limit (size). Conversely, by keeping the population limit (size) constant, the frequency of mutants decreases with decreasing mutation rate.

It thus looks like that the frequency is dependent on both the population limit (size) and the mutation rate. The frequency probably approaches 1 as the population size (limit) approaches infinity and limits towards 0 with the decreasing mutation rate.

The code below is a Python implementation of the above model. For each pair of the parameters 10000 simulations were performed.

import random


class Cell:
    """"Definition a cell object."""
    def __init__(self, mutation):
        self.mutation = mutation


def draw_mut():
    """Draw a mutation or not."""
    draw_m = random.choices(population=[True, False],
                            weights=[mut_rate, norm_rate], k=1)
    if draw_m[0]:
        mutation = True
    else:
        mutation = False
    return mutation


def cell_initialize():
    """Initialize the first cell."""
    mutation = draw_mut()
    cells.append(Cell(mutation))


def cell_division():
    """Divide cells."""
    for i in range(len(cells)):
        if cells[i].mutation:
            mutation = True
            cells.append(Cell(mutation))
        elif not cells[i].mutation:
            mutation = draw_mut()
            cells.append(Cell(mutation))


def average():
    """"Average frequency."""
    avg_freq = 0
    for frequency in frequencies:
        avg_freq += frequency
    return avg_freq / simulations


def run():
    """Run the simulations."""
    for n in range(simulations):
        mutants = 0
        global cells
        cells = []
        cell_initialize()

        while len(cells) < population_limit:
            cell_division()

        for cell in cells:
            if cell.mutation:
                mutants += 1
        frequencies.append(mutants / population_limit)
        if n + 1 == 2500:
            print(f"{round((n + 1) / simulations * 100)} % done.")
        elif n + 1 == 5000:
            print(f"{round((n + 1) / simulations * 100)} % done.")
        elif n + 1 == 7500:
            print(f"{round((n + 1) / simulations * 100)} % done.")
        elif n + 1 == 10000:
            print(f"{round((n + 1) / simulations * 100)} % done.")


populations = [2**4, 2**6, 2**8, 2**10, 2**12, 2**14, 2**16, 2**18]
scenarios = [0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001]

for pop in range(len(populations)):
    for sc in range(len(scenarios)):
        frequencies = []
        population_limit = populations[pop]
        mut_rate = scenarios[sc]
        norm_rate = 1 - mut_rate
        simulations = 10000
        run()

        print(f"Frequency: {average()}, mutation rate: {mut_rate}, "
              f"population: {population_limit}\n")

        with open("results.txt", "a") as f:
            f.write(f"Frequency: {average()}, mutation rate: {mut_rate}, "
                    f"population: {population_limit}\n")