Mean number of mutants in a culture of size N

For calculating the mean number of mutants in the culture, we used the following approach. We determined the mean number of division events that a cell in the final culture has experienced on the path from the cell that seeded the culture to a cell in the culture of size N. We then calculated the probability that a mutation occurred along this path.

To determine the mean number of divisions, we first determined the growth rate in the culture as a function of the cell-cycle time distribution, and then we used it to determine the mean generation number of a cell in the final culture. I will first outline the derivation of the growth rate in the culture. The age at which any individual cell divides is a random variable, A. For a cell randomly chosen at birth, this age of division is described by the cumulative distribution function $\Phi$ defined by Prob(A < a) = $\Phi(a)$, and, equivalently, by its density function $\phi(a) = d\Phi(a)/da$.

At division, the parent is lost and two cells of age zero are created. Consider a population of such cells. If the density of cells of age a in the population is denoted y(a,t) where t is the absolute time, then the equation for loss of the parent cells by division is

$$(\partial_t + \partial_a) y(a,t) = -h(a) y(a,t),$$ (6.4)

where h(a) is given by

$$h(a) = \frac{\phi(a)}{1 - \Phi(a)}.$$ (6.5)

Note that we can write

$$h(a) = -\frac{d}{da} \log(1- \Phi(a)).$$ (6.6)

The equation for gain of new cells is derived by integrating Eq.() over a and demanding that the total rate of production be given by

$$\frac{d}{dt} \int_0 ^\infty y(a,t) \:da = \int_0 ^\infty h(a) y(a,t) \: da.$$ (6.7)

This results in the production equation

$$y(0,t) = 2 \int_0 ^\infty y(a,t)\: h(a)\: da$$ (6.8)

We seek solutions to Eq.() of the form $y(a,t) = \eta(a) \rho(t)$, the so-called separable solutions. Substituting this form into Eq. (), we obtain

$$\frac{d\rho(t)}{dt}\frac{1}{\rho(t)} = -h(a) - \frac{d\eta(a)}{da}\frac{1}{\eta(a)}.$$ (6.9)

Note that the right-hand side depends only on a, while the left-hand side depends only on t. Therefore, if this equation is to hold for all values of both t and a, then either side must be constant. If we call this constant $\alpha$, we have

$$\frac{d \rho}{dt} = \alpha \rho$$ (6.10)

and

$$\frac{d \eta}{da} = \left[- h(a) - \alpha \right] \eta.$$ (6.11)

The solution for $\rho$ is simply

$$\rho(t) = \rho(0) e^{\alpha t}$$ (6.12)

while for $\eta$ we have

$$\eta(a) = \eta(0) (1-\Phi(a)) e^{-\alpha a}$$ (6.13)

The condition on $\alpha$ is obtained by substituting Eq. () into Eq. () to give

$$\eta(0) = -2 \eta(0) \int (1 - \Phi(a))\: e^{-\alpha a} \left[\frac{d}{da} \log( 1 -\Phi(a))\right] \:da.$$ (6.14)

Now an integration by parts and use of Eq. () yields

$$\frac{1}{2} = \int_0 ^\infty da\: \phi(a)\: e^{-\alpha a}$$ (6.15)

This is the eigenvalue equation for $\alpha$ that we seek. Since $\phi(a)$ is the density function for cell-cycle time, we can write this last result as

$$\mbox{E}[e^{-\alpha a}] = \frac{1}{2}$$ (6.16)

where E denotes expectation with respect to $\phi$. If $\phi$ is a gamma distribution of shape parameter q and mean $\tau$ (Eq. ), then the growth rate is given by

$$\alpha = \frac{q}{\tau} \left( 2^{\frac{1}{q}} - 1 \right).$$ (6.17)

Note that in the limits that we understand a priori, we have agreement with our expectations. For q=1, corresponding to an exponential density function, we should recover a simple Markov model. In this case, we get $\alpha = 1/\tau$. For the limit as $q \rightarrow \infty$, we get a process describing a polymerase chain reaction, with all "cells" replicating at exactly equal times. For this we have $\alpha = \log(2)/\tau$. Both of these results conform to prior knowledge.

The calculation of the mean generation number requires us calculating the mean age of a cell at division. For this, we assume that the culture was in stationary growth from the beginning. That is, we assume that the age distribution in the culture is constant as a function of time. We then calculate the average age at division using the density function for age, $\phi(a)$, but weighted by the proportion of cells of age a in the culture. To determine this, observe that cells that divide, by chance, earlier than usual, will leave, on average, more offspring than those that divide later. If the growth rate is g, i.e., the number of cells grows like $N(t)=N(0)\exp(\alpha t)$, then two cells that divide $\Delta t$ time units apart will leave different numbers of offspring and the ratio in that number is $\exp(\alpha \Delta t)$. For the simplicity of notation, I will denote $N(0) \equiv N_0$. Following this argument, first given by Fisher (1930), we obtain the average age at division

$$\mbox{E}[a\vert N]=\frac{\mbox{E}[a e^{-\alpha a}]}{\mbox{E}[e^{-\alpha a}]}$$ (6.18)

but in light of the definition of $\alpha$, we get

$$\mbox{E}[a\vert N]=2\mbox{E}[a e^{-\alpha a}].$$ (6.19)

For the specific case of the gamma distribution, parameterized as above

$$\mbox{E}[a\vert N]=\tau 2^{-1/q}.$$ (6.20)

The mean number of divisions is given by

$$\mbox{E}[g\vert N] = \frac{t}{\mbox{E}[a\vert N]} = \frac{\log(N/N_0)}{(\alpha \mbox{E}[a\vert N])}$$ (6.21)

which again in the case of the gamma distribution gives

$$\mbox{E}[g\vert N] = \frac{\log(N/N_0)}{q (1 - 2^{-1/q})}.$$ (6.22)

For the exponential distribution, $E[g\vert N] = 2 \log(N/N_0)$ and for the PCR process, $E[g\vert N] = \log(N/N_0)/\log(2) = \log_2(N/N_0)$.

We may now calculate the mean number of mutants by assuming that at each cell division, each of the daughters has a probability $\mu$ of becoming mutant (and therefore all of her daughters as well). The probability that no mutation occurs in g divisions is $(1 - \mu)^g$, so the probability that at least one mutation occurs is $1 - (1 - \mu)^g$. Then the probability of a cell being a mutant is $\sum_g (1 - (1 - \mu)^g) p(g)$, where p(g) is the probability that the cell underwent g divisions. For small mutation rates, such that $\mu g \ll 1$, this expression can be approximate by $\mu \left<g\right>$, that is, the probability that an individual cell is mutant is $\mu \mbox{E}[g\vert N]$. Each of the N cells in the final culture has this chance of being a mutant, thus the mean number of mutants in the culture is $\left<M\right> = \mu N \mbox{E}[g\vert N]$. For a gamma distribution of cell-cycle times, this becomes

$$\left<M\right> = c \mu N \log(N/N_0)$$ (6.23)

where the correction factor c is given by

$$c = \frac{1}{q(1- 2^{-1/q})}.$$ (6.24)

We now have an analytical form for the mean number of mutants in a culture of size N, in which the cells have a gamma-distributed cell cycle time.

A similar derivation gives us the correction factor for the mean in the case of a 2-phase model of cell cycle time. Assume that the cell cycle time has a constant component, T_B, and an exponentially distributed part, of parameter $\lambda$. Thus, the cell-cycle time distribution is

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

The mean cell cycle time is $T_B + 1/\lambda$, and the variance in cell cycle time is $1/\lambda^2$. The coefficient of variation of this cell-cycle time distribution is

$$r = \frac{1/\lambda}{T_B + 1/\lambda} = \frac{1}{\lambda T_B + 1}$$

Let us derive the mean proportion of mutants in the culture with this new cell-cycle time distribution. Assuming that the culture is in stationary growth, with growth rate $\alpha$, the number of cells in the culture as a function of time is given by

$$N = N_0 e^{\alpha t}.$$

Conform Eq. , the eigenvalue equation for the growth rate $\alpha$ is

$$\frac{1}{2} = \int_0^\infty \phi(a) e^{-\alpha a} da,$$

where $\phi(a)$ is the distribution of the age of cells at division. For the shifted exponential cell-cycle time distribution,

$$E[e^{-\alpha a}] = \int_{T_B}^\infty p(t) e^{-\alpha t} dt \\... ...t} dt \\ = \frac{\lambda e^{-\alpha T_B}}{\lambda + \alpha}.$$

Thus $\alpha$ must satisfy

$$\frac{\lambda e^{-\alpha T_B}}{\lambda + \alpha} = \frac{1}{2}.$$ (6.26)

The solution of this equation can be given in terms of the Lambert's W function:

$$\alpha = -\lambda + \frac{W[2 \lambda T_B e^{\lambda T_B}]}{T_B}.$$ (6.27)

The average age at division of a cell on the lineage from the root to a leaf in the genealogical tree of the culture is given by Eq. :

$$E[a\vert N] = \frac{E[a e^{-\alpha a}]}{E[e^{-\alpha a}]}.$$

From the eigenvalue equation for $\alpha$, we find (Eq. )

$$E[a\vert N] = 2 E[a e^{-\alpha a}].$$

The average number of generations from the founder cell to a cell in the current culture is then given by

$$E[g\vert N] = \frac{t}{E[a\vert N]} = \frac{\log(N/N_0)}{\alpha E[a\vert N]} = \frac{\log(N/N_0)}{2 \alpha E[a e^{\alpha a}]}.$$

Now

$$E[a e^{-\alpha a}] = \ \int_{T_B}^\infty a \lambda e^{-\lambd... ...T_B}}{\lambda + \alpha} \ (T_B + \frac{1}{\lambda + \alpha})].$$

This can be expressed in a simpler form, given that the growth rate, $\alpha$ satisfies Eq. . Namely,

$$E[a e^{-\alpha a}] = \frac{\lambda e^{- \alpha T_B}}{(\lambda... ...} = \frac{1}{2} \left[T_B + \frac{1}{\lambda + \alpha}\right].$$

As before, if we write the mean number of mutants in the culture as

$$\left<M\right> = c \mu M \log(N/N_0),$$

where c is a function of the cell-cycle time distribution, for the shifted exponential we obtain:

$$c = \frac{\lambda + \alpha}{\alpha\left[T_B\left(\lambda + \alpha \right) + 1 \right]}.$$ (6.28)

As expected, if we set T_B = 0, we obtain c = 2, which is the correction factor for exponentially-distributed cell cycle time.

 

Table 6.1: Mean proportion of mutants in cultures in which cells have gamma-distributed cell cycle time.

 

N0	N	$\mu$	q	Observed mean (S.E.)	Predicted mean
1	104	10-4	1	0.001742 (0.000154)	0.001842
1	104	10-4	3	0.001417 (0.000134)	0.001488
1	104	10-4	10	0.001330 (0.00011)	0.001375
1	104	$3 \times 10^{-4}$	1	0.005128 (0.000207)	0.005526
1	104	$3 \times 10^{-4}$	3	0.004045 (0.000151)	0.004464
1	104	$3 \times 10^{-4}$	10	0.004355 (0.000204)	0.004126
1	104	10-3	1	0.017548 (0.000463)	0.018421
1	104	10-3	3	0.014822 (0.000366)	0.014882
1	104	10-3	10	0.013536 (0.000327)	0.013754
1	104	$3 \times 10^{-3}$	1	0.050005 (0.00068)	0.055262
1	104	$3 \times 10^{-3}$	3	0.043824 (0.000625)	0.044645
1	104	$3 \times 10^{-3}$	10	0.041336 (0.000586)	0.041261
1	105	10-4	1	0.002223 (0.000168)	0.002303
1	105	10-4	3	0.001702 (0.000084)	0.00186
1	105	10-4	10	0.001354 (0.000086)	0.001719
1	105	10-3	1	0.023307 (0.000524)	0.023026
1	105	10-3	3	0.017987 (0.000349)	0.018602
1	105	10-3	10	0.016823 (0.000315)	0.017191

 

Table 6.2: Mean proportion of mutants in cultures in which the cell cycle time is distributed as a shifted exponential.

 

N0	N	$\mu$	r	Observed mean (S.E.)	Predicted mean
1	104	10-4	1	0.001784 (0.000176)	0.001425
1	104	10-4	3	0.001338 (0.000107)	0.001353
1	104	10-4	9	0.00132 (0.00011)	0.001333
1	104	$3 \times 10^{-4}$	1	0.004549 (0.000221)	0.004268
1	104	$3 \times 10^{-4}$	3	0.00397 (0.000176)	0.00402
1	104	$3 \times 10^{-4}$	9	0.003866 (0.000146)	0.004
1	104	10-3	1	0.014188 (0.000339)	0.014225
1	104	10-3	3	0.013547 (0.000317)	0.013533
1	104	10-3	9	0.013415 (0.00033)	0.01333
1	104	$3 \times 10^{-3}$	1	0.041419 (0.000556)	0.042676
1	104	$3 \times 10^{-3}$	3	0.039641 (0.000548)	0.040599
1	104	$3 \times 10^{-3}$	9	0.038854 (0.00052)	0.039991
1	105	$3 \times 10^{-4}$	1	0.005331 (0.000178)	0.005335
1	105	$3 \times 10^{-4}$	3	0.005005 (0.000177)	0.005075
1	105	$3 \times 10^{-4}$	9	0.005172 (0.000168)	0.004999
1	105	10-3	1	0.017139 (0.000323)	0.017782
1	105	10-3	3	0.017509 (0.00036)	0.016916
1	105	10-3	9	0.016393 (0.00029)	0.016663

The mean proportion of mutants as obtained from simulations, together with the theoretical prediction is presented in Tables (for the gamma-distributed cell cycle time) and (for the shifted exponential). 10000 independent runs were performed for each of the parameter sets.

As can be seen from the tables, there is a good agreement between the means that I obtained from simulations, and the ones that I calculated. It is also apparent that if the cell cycle is exponentially-distributed, the mean proportion of mutants is higher, for the same mutation rate per division, than if the cell-cycle time distribution had a higher order. Turning the argument around, if we assume that the cell cycle time is exponentially-distributed, when in reality it is not, leads to underestimation of the mutation rate. Although this comes out of the expression for the mean, I would like to give an intuitive argument for how the cell-cycle time distribution enters into the mutant distribution.

Assuming that the cells have an exponentially-distributed cell cycle time is equivalent to assuming that they all have a constant probability of dividing per unit time. This seems reasonable at first, but of course, it is false. Any cell that has just divided will have very small probability of dividing again too soon. That this actually makes a difference in the distribution of mutants in a population of a given size is seen clearly when considering the case of a population of four cells arising from a single ancestor. There are only two topologically distinct genealogical trees (Fig. ) for the cells in this culture.

 

In the balanced tree (Fig. A), each of the four individuals has two divisions in their history and so has the same probability of being mutated: $2 \mu$ for $\mu$ small. In the second tree, which is skewed (Fig. B), the four individuals went through 1,2,3 and 3 division events, with probability $\mu$ of mutating at each division, for a mean mutant frequency of $9/4 \mu$. Also, there is clearly a larger variance compared to the balanced tree. To make the point even clearer, consider the polymerase-chain reaction. Here one starts with a small number, for simplicity say one, molecules of a nucleic acid, called template. By adding a polymerization enzyme, and energy-rich nucleotide monomers, the complementary strand will be synthesized for the initial template molecule. Then the complementary strands are dissociated from one another and the reaction is repeated for n cycles, to yield 2ⁿ molecules of nucleic acid. At each cycle, all the molecules in the vat act as templates, and the complementary strand is synthesized for each of them. It is clear that only the completely balanced tree will be realized in this type of reaction. Theoretically, the probability of obtaining a given genealogical tree can be computed given the distribution of cell cycle times. The proportion of mutants in the final culture will depend only on the relative probabilities of realizing different types of genealogical trees with the same number of leaves, and on the probability of mutation at cell division. The problem is that this computation becomes intractable for even very small trees.

The approach that we eventually designed for accounting for the cell cycle time in the growth of the culture was suggested by our findings that:

• the mean proportion of mutants depends on the growth rate of the culture, and
• that the growth rate is only a function of the cell cycle parameters.

Analyzing the empirical distributions of mutants that we obtained for various cell cycle parameters, we found that they closely resemble in shape the Luria-Delbrück distribution. This prompted us to attempt to generalize a variant of the Luria-Delbrück distribution for cell-cycle time distributions other than exponential. This variant is a continuum approximation of L-D, due to my collaborator, T. Kepler (Kepler & Oprea, in preparation). In the next section, I present a brief outline of the derivation of this distribution.