Computational model of a growing culture of cells

The problem that I am trying to solve is to estimate the mutation rate in a culture of cells which is undergoing exponential growth and phenotypic mutations. The basic setup assumed by fluctuation analysis is the following. A culture is seeded with N₀ (generally 1) cells. The culture then grows exponentially up to N cells. When a wildtype cell divides, each of its progeny has a probability $\mu$ of undergoing a mutation that changes its phenotype. In bacteriological experiments, the change in phenotype generally means that the cell will be capable of using a nutrient that a wildtype cell could not metabolize. However, the assay for detecting the mutants is only performed after the growth of the culture. That means that, while the culture is growing, mutant cells do not have a selective advantage over wildtype cells. They grow at the same rate as wildtype cells. It is also assumed that mutants do not revert to the wildtype phenotype, and thus all progeny of a mutant cell will be mutants. When the culture has reached size N, we count the number of mutants, M. Note that lethal mutations are neglected in this analysis. In fact, all the mathematical treatments of this process neglect lethal mutations, as they probably affect only the growth rate of the culture, not the relative proportion of mutant and wildtype cells. Also, neutral mutations fall under the wildtype phenotype, so they do not need a separate mathematical treatment. The question is now, if we have a number of data points (N_i, M_i), from different cultures that all have presumably the same mutation rate, how do we estimate the mutation rate from these data?

I set up an event driven simulation of the process that I just described. Each cell is represented by an object characterized by its phenotype, wildtype or mutant, and a division time. I seed the system with one wildtype cell of age 0. When a new cell object is created, I assign it a cell cycle time. I explored both types of cell-cycle time distributions that have been used in the literature. The first is the gamma distribution of order q, and scale parameter $\theta$:

$$\phi[t\vert q,\theta] = \frac{\theta^q t^{q-1}}{\Gamma(q)} e^{-\theta t}.$$ (6.1)

The mean cell cycle time is $\left<t\right> = \frac{q}{\theta}$ and its variance $var(t) = \frac{q}{\theta^2}$. For simplicity, I set the scale parameter to 1 in all cases. This will only be reflected in the absolute values of the cell cycle time, not in their relative ordering. The second type of cell-cycle time distribution that I used I call shifted exponential. That is, there is a constant probability per unit time that the cell starts to divide, $\lambda$, but division takes a constant amount of time for all cells, T_B. Then the distribution of cell cycle time is

$$p(t) = \left\{ \begin{array}{ll} 0 & \mbox{if $t < T_B$ } ... ...xp(-\lambda (t - T_B)) & \mbox{otherwise} \end{array} \right.$$ (6.2)

The mean cell cycle time is in this case $T_B + 1/\lambda$, and the variance in cell cycle time is $1/\lambda^2$. I scaled these parameters such that the mean cell cycle time is 1. That is, only the ratio between the division time and the mean waiting time between two divisions will affect the mutant distribution in the culture, not the absolute values of these times. I will denote this parameter by

$$r = \frac{T_B}{\frac{1}{\lambda}} = \lambda T_B.$$ (6.3)

I maintain a priority queue of cell objects, the value used for determining the order of objects in the queue being the absolute value of time at which the cell divides. This in turn is the sum of the cell cycle time of the cell and the time at which the cell was born. At each step of the simulation, the object with the lowest time of division is removed from the queue. Two new objects are created, each having a division time which is the sum of the current time and a cell cycle time drawn from the gamma distribution. With probability $\mu$ each of the two daughter cells mutates. These operations are performed for a fixed number of steps, i.e., for N-1 duplications, if we are to achieve a culture of size N. This algorithm implements the dynamics of a culture that grows exponentially from 1 to N cells. Random deviates from a gamma distribution of a given order were generated using the standard Numerical Recipe function (1988). The special case of the shape parameter q = 1 gives us deviates from the exponential distribution. At the end of the simulation, I count the number of mutants among the N cell objects. I generate a large number (10⁴-10⁶) of replicates of this experiment for constructing the distribution of the proportion of mutants in the culture.