Estimating the Rate of Mutation to a Mutator Phenotype

1. Introduction

Replication fidelity is paramount to the survival of species and to the origin of complex life [1,2]. The fidelity required for survival increases with the complexity of an organism [1], the length of its genome [1,3,4], and the genetic linkage within that genome [5]. That such high fidelity of replication is vital to survival is evidenced by the observation that, across the tree of life, large fractions of organisms’ genomes consist of genes dedicated to aspects of replication, proofreading and repair that insure high-fidelity replication [6,7,8]. Any mutation appearing within any of these many genes has the potential to affect replication fidelity. In this light, it is not surprising that variation in replication fidelity has been observed in many species [9,10,11]. The amount of variation in mutation rate, however, is not necessarily reflective of the rate of mutation to mutation-rate variants, as indirect selection on the mutation rate in nature is a major confounding factor that is variable. Additionally, the number of mutations occurring in replication and repair genes that actually affect replication fidelity is largely unknown. And of those mutations that do affect fidelity, the magnitude of that effect is also largely unknown. This deficits in our knowledge surrounding the mutation rate to mutation-rate variants is primarily due to methodological difficulties and certainly not to a lack of interest or importance. Filling these gaps in our knowledge would impact both basic and applied sciences, shedding new light on basic evolutionary processes and providing key parameters that have been missing from cancer biology, bacterial evolution, and viral evolution. The present manuscript develops a novel method for estimating these missing parameters. It takes advantage of an under-appreciated way to increase statistical power in such assays, dubbed the “Jones protocol” [12,13,14]. We test our methods against computer simulations of feasible laboratory experiments using biologically reasonable parameters, and we find that it may indeed be possible to empirically estimate rates of mutation to mutation-rate variants. We briefly discuss application of these methods to the study of organismal and somatic evolution.

The observation that all eukaryotes have costly mechanisms for proofreading and repair suggests that low mutation rates must be advantageous and perhaps even essential for the survival of species. Even the large RNA virus SARS-CoV-2 has primitive proofreading/repair machinery, providing a large mutational target for mutation-rate variation [15]. A fairly simple information-theoretic argument [1,16] predicts an inverse relationship between the length of an organism’s genome and the maximum tolerable mutation rate allowed (or minimum replication fidelity required) for the organism’s long-term survival [4,16,17,18,19,20,21].

One way organisms can reduce the effective length of their genome and thereby tolerate higher mutation rates is by introducing recombination [22]. The observation that asexual species on the tree of life are rare and relatively short-lived may be interpreted as evidence for the evolutionary advantage of sex and recombination [23,24]; but it can also be interpreted as evidence for the evolutionary advantage of low genomic, or more precisely, low effective-genomic, mutation rates. Under free recombination, the effective-genomic mutation rate is the genic mutation rate, i.e., the per-gene mutation rate, which can be much lower than the per-genome mutation rate.

A previous study [25] estimated the rate of mutation to mutator phenotype made possible by an observation that mutators would form colonies on plates that looked qualitatively different from non-mutators (or wildtype). This study depended, however, on the judgment of a human being peering into a microscope. Furthermore, the authors of this study employed one type of cell and made no claims that such distinguishable colony morphology might be observed with other cell types. The methodology we present here does not depend on cell type, colony morphology or on the accuracy of a possibly overworked lab technician.

In addition to basic evolutionary questions our experimental protocol can address, such as those outlined above, there are also significant biomedical applications that can be addressed, not least of which is cancer research [26,27,28,29,30,31,32,33,34]. The rate at which mutation rate variants are produced in populations of somatic cells is a key parameter for understanding and predicting somatic evolution leading to oncogenic transitions; despite its critical role, however, this parameter remains unknown. The work we present here may help to remedy this deficiency.

2. Materials and Methods

2.1. Definitions

$N_0,N$	initial and final population size
r	growth rate
$\mu$	marker mutation rate of wildtype
${\mu }^{\prime }$	marker mutation rate of mutator
${\mu }_m$	rate of mutation from wildtype to mutator
m	factor by which wildtype and mutator mutation rates differ: ${\mu }_m=m\mu$
$\varphi \left(t\right)$	recruitment rate of mutations; e.g. $\varphi \left(t\right)=\mu N_0{e}^{rt}$ for marker recruitment rate on the wildtype background
$p(k;t)$	probability that a mutation appearing at time t leaves k mutants in the final population
${\lambda }_k$	expected number of mutations of “type k”, i.e., that leave k mutants in the final population
$t_f$	time at which population attains its maximum size: $N=N\left(t_f\right)$
$\Lambda$	total number of marker mutations that occur on the wildtype background: $\Lambda =\mu (N-N_0)\approx \mu N$
${\Lambda }_m$	total number of mutator mutations that occur on the wildtype background: ${\Lambda }_m={\mu }_m(N-N_0)\approx {\mu }_{m}N$
$\phi (x;\Lambda )$	pgf for the total number of mutants observed in the final population, given that a total of $\Lambda$ mutations occurred during growth.
$\check{f}$	fast Fourier transform of function f.
${\mathcal{L}}_i$	likelihood function for the i-th growth cycle.

A “mutation” is an error that occurs during replication. A “mutant” is an individual that carries a “mutation”. Mutations of interest here are specific mutations (replication errors) that can be observed experimentally. For example, two ways in which specific mutations can be observed are: 1) by eliminating everyone in the population that does not carry the specific mutation (for example, applying an antibiotic to screen for antibiotic-resistant mutations), and 2) by causing a fluorescent marker to fluoresce. Such observable mutations are sometimes called “selectable markers”, and we will refer to such mutations as “marker mutations”, and to the rate at which such mutations occur as “marker mutation rate”. When a mutation occurs, it becomes the founding member of a new “mutant lineage” in which all descendants will carry that “mutation”. In what follows, we use the terms “mutation” and “mutant lineage” interchangeably. For example, we will introduce the concept of a “type-k mutation”, referring to a marker mutation that forms a mutant lineage that ultimately gives rise to k “marker mutants” in the final culture. (We could equivalently say “type-k mutant lineages”.)

By “mutator”, we refer to a genotype with some mutation in its replication/proofreading/repair mechanisms that causes it to replicate with reduced accuracy. And here we digress briefly to explain the title of our paper: “Estimating the rate of mutation to a mutator phenotype”. A more accurate title might be “Estimating the rate of mutation to a mutation-rate-variant phenotype”, as our methods could in theory detect a decrease in mutation rate as well as an increase. We have kept the former title, however, because mutator mutations (mutations that increase mutation rate) are typically loss-of-function mutations and are therefore overwhelmingly more common than mutations that decrease mutation rate (which are gain-of-function mutations).

2.2. Experimental Protocol

• Grow a bacterial population from initial size $N_0$ to a large (known) number, N, at time T.
• At time T:
  • Protocol A: Take two random samples of size $N_S$ and $N_B$, where $N_S\ll N_B\ll N$ (S for sample, and B for bottleneck). Use the sample of size $N_S$ to inoculate a number c of flasks with fresh media; these independent cultures grow and each is then screened for marker mutants. Record the number of marker mutants observed in each of the c cultures. Use the sample of size $N_B$ to inoculate fresh media in a single flask to start the next growth cycle.
  • Protocol B: Take one random sample of size $N_B$ to inoculate fresh medium to start the next growth cycle.
• Repeat Step 2 $n-1$ times (for n growth cycles)
• After the final $n^{th}$ growth cycle, take one random sample of size $N_S$. Use this sample to inoculate a number c of flasks with fresh media; these independent cultures grow and each is then screened for marker mutants. Record the number of marker mutants observed in each of the c cultures.

Figure 1. Schematic of the two experimental protocols. To maximize statistical power, the serially-passaged populations (i.e., the populations at each growth cycle) should be grown to the largest N possible, then diluted to obtain countable numbers of mutants on the selective plates [12]. This takes advantage of the increased statistical power conferred by the “Jones protocol” [12,13,14], an improvement on the standard fluctuation test [35,36,37,38,39].

Figure 1. Schematic of the two experimental protocols. To maximize statistical power, the serially-passaged populations (i.e., the populations at each growth cycle) should be grown to the largest N possible, then diluted to obtain countable numbers of mutants on the selective plates [12]. This takes advantage of the increased statistical power conferred by the “Jones protocol” [12,13,14], an improvement on the standard fluctuation test [35,36,37,38,39].

2.3. Analysis

For simplicity of presentation, we will assume that the population grows exponentially during the growth phase, such that the initial population size is $N_0$, and the final population size at time $t_f$ is given by $N=N_0{e}^{r{t}_f}$. We note that our results do not depend on this assumption: we only require that the population grows from size $N_0$ to N by means of some pure-birth process (meaning that the total number of replications is exactly $N-N_0$). We let $\varphi \left(t\right)$ denote "recruitment rate" of mutations at time t, i.e., $\varphi \left(t\right)=\mu N_0{e}^{rt}$, where $\mu$ is mutation rate per individual per unit time. We let $p(k;t)dt$ denote the probability that a mutation appearing in the small time interval $(t,t+dt)$ leaves k mutants in the final population.

Table 1. Simulated data from protocol B with three growth cycles. Parameters were $N_0=1$, $N=2^{33}$, $N_B=2^{20}$, $N_S=4295$

Numbers of marker mutants
	Growth cycle
Culture	1	2	3
1	1	6	9
2	1	9	13
3	4	8	19
4	23	8	15
5	5	8	13
6	3	5	14
7	12	3	9
8	1	7	7
9	8	8	13
10	4	5	7
11	1	7	8
12	3	10	12

Table 1. Simulated data from protocol B with three growth cycles. Parameters were $N_0=1$, $N=2^{33}$, $N_B=2^{20}$, $N_S=4295$

Numbers of marker mutants
	Growth cycle
Culture	1	2	3
1	1	6	9
2	1	9	13
3	4	8	19
4	23	8	15
5	5	8	13
6	3	5	14
7	12	3	9
8	1	7	7
9	8	8	13
10	4	5	7
11	1	7	8
12	3	10	12

We will refer to a mutation (or mutant lineage) as being of “type k” if it leaves k mutants in the final population. We define random variable $M_k$ to denote the number of mutations (or mutant lineages) of type k. We make the very weak assumption that the $M_k$ are independent, in which case, we have:

$$\mathbb{P}\{M_k=j\}=\frac{{\lambda }_k^j}{j!}{e}^{-{\lambda }_k}$$

(1)

where

$${\lambda }_k={\int }_0^{t_f}p(k;t)\varphi \left(t\right)dt.$$

The total number of mutants in the final population is:

$$M=\sum _{k=1}^{\infty }k{M}_k.$$

For a "pure birth" process of constant birth rate r, and assuming that the mutations have no effect on fitness, we know that:

$$p(k;t)=e^{-r(t_f-t)}{\left(1-e^{-r(t_f-t)}\right)}^{k-1}$$

([40], p 450). While we assume r to be a constant here, we note that the results we obtain are robust to this assumption; i.e., any pure-birth regime ($r=r\left(t\right)>0$) should in principle give the same results. This fact may be intuited by recalling that we implicitly make a very weak assumption (i.e., a very reasonable assumption from a biological perspective) that mutation occurs primarily during replication. Under this assumption, what counts is simply the number of replication events and not the dynamics of replication rates. This assumption can be violated for example, when ultra-violet radiation is the primary source of mutation, but we suppose that such scenarios are exceptional. From here we have:

$$\begin{array}{cc}\hfill {\lambda }_k& ={\int }_0^{t_f}\mu N_0{e}^{rt}p(k;t)dt\hfill \\ & =\frac{\mu N}{r}{\int }_{e^{-r{t}_f}}^{1}u{(1-u)}^{k-1}du\hfill \\ & \approx \frac{\mu N}{k(k+1)},\hfill \end{array}$$

given that $r\approx 1$, $e^{-r{t}_f}=\frac{N_0}{N}\to 0$, and

$${\int }_0^{1}u{(1-u)}^{k-1}du=\frac{1}{k(k+1)}.$$

Now that we have derived the overall expected number of mutant lineages of type k, ${\lambda }_k$, we return to equation 1 and write the pgf for the number of type k mutant lineages as:

$$f_k\left(x\right)=e^{{\lambda }_k(x-1)}.$$

Since each of these lineages contributes exactly k mutants to the final population, the pgf for the total number of mutants contributed by type k lineages to the final population is simply:

$$f_k\left(x^k\right)=e^{{\lambda }_k(x^k-1)}.$$

The pgf for the total number of mutants observed in the final population, from all types, is then given by

$$\phi \left(x\right)=\prod _{k=1}^{\infty }{e}^{{\lambda }_k(x^k-1)}.$$

Substituting the expression for ${\lambda }_k$ we find:

$$\phi (x;\Lambda )=e^{-\Lambda }exp\left(\sum _{k=1}^{\infty }\frac{\Lambda x^k}{k(k+1)}\right),$$

(2)

where $\Lambda =\mu (N-N_0)\approx \mu N$. The notation $\phi (x;\Lambda )$ emphasizes that the pgf for the total number of observed mutants at the final time depends on a single parameter, the product of the final population size and the mutation rate.

We can simplify this expression by considering the following lemma and theorem.

Lemma 1.Given $z\in [0,1]$ then

$$\begin{array}{cc}\hfill \sum _{j=1}^{\infty }& \frac{z^j}{j(j+1)}=\left\{\begin{array}{cc}1+\left(\frac{1-z}{z}\right)log(1-z),\hfill & \mathrm{if}z\in [0,1).\hfill \\ 1,\hfill & \mathrm{if}z=1.\hfill \end{array}\right.\hfill \end{array}$$

Proof. 

If $z=1$ we have that

$$\begin{array}{cc}\hfill \sum _{j=1}^{\infty }\frac{z^j}{j(j+1)}& =\underset{k\to \infty }{lim}\sum _{j=1}^k\left(\frac{1}{j}-\frac{1}{j+1}\right)\hfill \\ \hfill & =\underset{k\to \infty }{lim}\left(\frac{1}{1}-\frac{1}{2}\right)+\left(\frac{1}{2}-\frac{1}{3}\right)+\cdots +\left(\frac{1}{k}-\frac{1}{k+1}\right)\hfill \\ \hfill & =\underset{k\to \infty }{lim}1-\frac{1}{k+1}\hfill \\ \hfill & =1.\hfill \end{array}$$

For case $z\in [0,1)$, let’s recall that the Taylor series of $log(1+\zeta )$ is

$$\begin{array}{c}\hfill T\left(\zeta \right)=\sum _{j=1}^{\infty }\frac{{(-1)}^{j+1}}{j}{\zeta }^j,\zeta \in (-1,1).\end{array}$$

Which implies that

$$\begin{array}{cc}\hfill log(1-\zeta )& =T(-\zeta )\hfill \\ \hfill & =\sum _{j=1}^{\infty }\frac{{(-1)}^{2j+1}}{j}{\zeta }^j\hfill \\ \hfill & =-\sum _{j=1}^{\infty }\frac{{\zeta }^j}{j},\hfill \end{array}$$

(3)

for each $\zeta \in (-1,1)$. Therefore

$$\begin{array}{cc}\hfill \sum _{j=1}^{\infty }\frac{z^j}{j+1}& =\frac{1}{z}\sum _{j=1}^{\infty }\frac{z^{j+1}}{j+1}\hfill \\ \hfill & =\frac{1}{z}\left(z-z+\sum _{j=2}^{\infty }\frac{z^j}{j}\right)\hfill \\ \hfill & =\frac{1}{z}\left(-z+\sum _{j=1}^{\infty }\frac{z^j}{j}\right)\hfill \\ \hfill & =-1-\frac{1}{z}log(1-z).\hfill \end{array}$$

(4)

From Eq. (4) and Eq. (3) it follows that

$$\begin{array}{cc}\hfill \sum _{j=1}^{\infty }\frac{z^j}{j(j+1)}& =\sum _{j=1}^{\infty }\left(\frac{1}{j}-\frac{1}{j+1}\right)z^j\hfill \\ \hfill & =\sum _{j=1}^{\infty }\frac{z^j}{j}-\sum _{j=1}^{\infty }\frac{z^j}{j+1}\hfill \\ \hfill & =-log(1-z)+\frac{1}{z}log(1-z)+1\hfill \\ \hfill & =1+\left(\frac{1-z}{z}\right)log(1-z),\hfill \end{array}$$

as stated.    □

Theorem 1.Given $z\in [0,1]$ and under the convention that ${\displaystyle \underset{(x,y)\to (0,0)}{lim}{x}^y=1}$ then

$$e^{-\Lambda }exp\left(\sum _{j=1}^{\infty }\frac{\Lambda z^j}{j(j+1)}\right)={(1-z)}^{\frac{\Lambda (1-z)}{z}}.$$

(5)

Proof. 

In the limit as $z\to 1$, Eq. (5) holds by the convention.

Let $z\in [0,1)$ be a fixed number and let ${\displaystyle \psi =\sum _{j=1}^{\infty }\frac{z^j}{j(j+1)}}$. By Lemma 1, we have that

$$\begin{array}{cc}\hfill e^{\Lambda \psi }& =exp\left(\Lambda +\Lambda \left(\frac{1-z}{z}\right)log(1-z)\right)\hfill \\ \hfill & =e^{\Lambda }exp\left(\Lambda \left(\frac{1-z}{z}\right)log(1-z)\right)\hfill \\ \hfill & =e^{\Lambda }exp\left(log\left({(1-z)}^{\Lambda \left(\frac{1-z}{z}\right)}\right)\right)\hfill \\ \hfill & =e^{\Lambda }{(1-z)}^{\Lambda \left(\frac{1-z}{z}\right)}.\hfill \end{array}$$

Thus

$$e^{-\Lambda }{e}^{\Lambda \psi }={(1-z)}^{\Lambda \left(\frac{1-z}{z}\right)}.$$

   □

Applying Theorem 1 to Equation 2, we see that if a population grows from a single individual to final size N, creating de novo (neutral) mutant lineages at rate $\mu$ per individual per unit time, then the total number of mutant individuals in the population at the final time will be described by the pgf:

$$\phi (z;\Lambda )={(1-z)}^{\frac{\Lambda (1-z)}{z}},$$

(6)

which we will also write as:

$$\phi (z;\Lambda )=h{\left(z\right)}^{\Lambda },$$

(7)

where

$$h\left(z\right)={(1-z)}^{\frac{1-z}{z}}.$$

(8)

Finally, we note that Equation 7 gives the pgf for the number of mutant individuals in a final population of N individuals. Trivially, we re-write as $\phi (z;\Lambda )={\left(h{\left(z\right)}^{\mu }\right)}^N$; we find that $h{\left(z\right)}^{\mu }$ gives the pgf for numbers of mutants that a single replication event will leave in the final population. We will use this result in interpreting the expression for $F\left(z\right)$ at the end of the following section.

2.3.1. Incorporating mutators

The foregoing theory does not account for the fact that, as a population grows, variants carrying an elevated mutation rate (mutators) may appear as a consequence of mutations occurring in genes encoding any aspect of replication, proofreading or repair. We now modify the foregoing theory to account for mutators.

Let ${\mu }^{\prime }>\mu$ denote the mutation rate of the mutator to observable mutants (i.e., to the selectable phenotype). Consider a mutator lineage of type k (i.e., a clone that leaves k mutators in the final culture). Since this lineage grows to final size k and produces observable mutants at rate ${\mu }^{\prime }$, we can use Equation 6 directly. The number of observable mutants in the final population that were produced within a mutator clone of type k therefore has pgf

$$\phi (z;k{\mu }^{\prime })={(1-z)}^{k{\mu }^{\prime }(1-z)/z},$$

where this distribution depends again on the product of the final population size – in this case k – and the mutator mutation rate, ${\mu }^{\prime }$. If there are j mutator clones of type k, then the number of observable mutants produced by type-k clones has pgf:

$$\phi {(z;k{\mu }^{\prime })}^j={(1-z)}^{jk{\mu }^{\prime }(1-z)/z}$$

And the total number of observable mutants produced by type-k mutators has pgf:

$$\begin{array}{cc}\hfill \mathbb{E}\left[\phi (z;k{\mu }^{\prime })\right]& =\sum _{j=0}^{\infty }\phi {(z;k{\mu }^{\prime })}^j\mathbb{P}\{M_k=j\}\hfill \\ \hfill & =exp\left\{{\lambda }_k(\phi (z;k{\mu }^{\prime })-1)\right\}\hfill \end{array}$$

The total number of observable mutants produced by mutators of all types is therefore described by the pgf:

$$f\left(z\right)=\prod _{k=1}^{\infty }\mathbb{E}\left[\phi (z;k{\mu }^{\prime })\right]$$

which may be rewritten as:

$$logf\left(z\right)=\sum _{k=1}^{\infty }{\lambda }_k(\phi (z;k{\mu }^{\prime })-1)=\sum _{k=1}^{\infty }{\lambda }_k(h{\left(z\right)}^{k{\mu }^{\prime }}-1)$$

where

$$h\left(z\right)={(1-z)}^{(1-z)/z}$$

We define ${\Lambda }_m={\mu }_{m}N$ to be the recruitment rate of mutators in the population, where ${\mu }_m$ is the rate of mutation to mutators. Substituting the expression for ${\lambda }_k$, using Lemma 1 and letting $g\left(z\right)=h{\left(z\right)}^{{\mu }^{\prime }}$, we find:

$$logf\left(z\right)={\Lambda }_m\sum _{k=1}^{\infty }\frac{g{\left(z\right)}^k-1}{k(k+1)}=\frac{{\Lambda }_m(1-g\left(z\right))log(1-g\left(z\right))}{g\left(z\right)}$$

The pgf for the total number of observable mutants produced by mutators may thus be written as:

$$f\left(z\right)={(1-g\left(z\right))}^{{\Lambda }_m(1-g\left(z\right))/g\left(z\right)}=h{\left(h{\left(z\right)}^{{\mu }^{\prime }}\right)}^{{\Lambda }_m}.$$

Since the pgf describing the total number of observable mutants produced by the wildtype is given by $h{\left(z\right)}^{\Lambda }$ (Eq 6), the total number observable mutants produced in the culture (i.e., by either wildtype or mutator) therefore has pgf:

$$F\left(z\right)=h{\left(z\right)}^{\Lambda }h{\left(h{\left(z\right)}^{{\mu }^{\prime }}\right)}^{{\Lambda }_m}$$

This expression has a clear intuitive explanation. We know that $h{(·)}^{{\Lambda }_m}$ describes the total number of mutators produced by the wildtype in the final population. Since the pgf for the number of observable mutants produced by a single replication event within a mutator lineage is given by $h{\left(z\right)}^{{\mu }^{\prime }}$, the composition $h{\left(h{\left(z\right)}^{{\mu }^{\prime }}\right)}^{{\Lambda }_m}$ gives the $pgf$ associated with the total number of observable mutants in the final population produced by mutators. The product of this term with the pgf for the observable mutants produced by the wildtype, $h{\left(z\right)}^{\Lambda }$, then gives the pgf for the total number of observable mutants in the final population.

2.4. Model fitting and validation

Our overall goal is to use the analytical results described above to estimate the underlying mutation rates in populations that have been propagated as described in Section 2.2. To test this idea, we generate simulated data using three independent approaches, which are described in detail, with pseudocode, in Appendix A. In Section 3.1 we test the validity of these simulated data sets, using the simple case of a single bottleneck and no mutators as a test case. We then use all three methods to produce datasets including mutators. Finally, in Section 3.2, we describe and test the accuracy of three methods to estimate the mutation-to-mutator rate from the simulated data.

3. Results

3.1. Model reliability

Since our theoretical approaches can be regarded as numerical methods, we validated reliability by performing simulations. It is widely known that the Moran model can be used to simulate cell growth where mutations are allowed [41, Figure 5B]. Here, we employ the Moran model to model a growing population; this differs from common usage of the Moran model but is equally valid.

A cell population where only observable (marker) mutations are allowed, such that no mutators or mutant cells with mutator phenotypes appear, can also be modeled with the graphical model (described in Appendix A.1) and the quantile function model (described in Appendix A.2), stopping them at Eq. (A5) and Eq. (A18), respectively. Additional growth cycles can be simulated by implementing a random sampling of the desired size to found the next growth cycle.

For each of the final population sizes $N\in \{20000,40000,60000,80000,100000,200000,300000,400000,500000,600000,700000,800000,900000,1000000\}$, 100 simulations were generated by the Moran, graphical, and the quantile function models, respectively. Goodness-of-fit tests were carried out to determine whether the null hypothesis, $H_0$: the data sets are drawn from the same distribution, can be rejected. Pairwise two-sample Kolmogorov-Smirnov tests for discrete distributions [42] were performed, as used in Ycart (2013) [43].

As shown in Figure 2, the distributions associated with the simulations generated by Moran, graphical, and quantile function models were not significantly different $(p>0.05)$ for all final population sizes tested.

Furthermore, the underlying distribution of the data sets is known since the number of mutants produced by wildtype cells has the Luria-Delbrück distribution [44, Eq. 15]. For each data set with a final population size $N=1000000$, used for Figure 2, analysis on the web tool bz-rates [45] was performed.

The bz-rates platform provides visual verification of the fitting of data sets to the Luria-Delbrück distribution, an estimation of the mutation rate, the 95% confidence interval for the mutation rate, etc. In addition, a ${\chi }^2$ goodness-of-fit test is performed to determine whether the null hypothesis, $H_0$: the data set was drawn from a Luria-Delbrück distribution, can be rejected.

As shown in Figure 3, it was not possible to reject the null hypothesis that the data sets were drawn from a Luria-Delbrück distribution, with the Luria-Delbrück and data set distributions not being statistically different ($p>0.05$), see Table 2.

Figure 2. Statistical comparison of data set distributions. Given a final population size N, $N\in \{20000,40000,60000,80000,100000,200000,300000,400000,500000,600000,700000,800000,900000,1000000\}$, 100 simulations were generated with each of the models. For the three models, the parameters were $\mu ={10}^{-7}$, ${\mu }^{\prime }=0$, and ${\mu }_m=0$; hence only mutant cells produced by wildtype cells were allowed. The horizontal red dashed line shows the upper bound of the rejection region ($p=0.05$) for the null hypothesis $H_0$ which states that both data sets come from the same distribution. (a) p-values of the Kolmogorov-Smirnov test, per final population size, between the data sets of the Moran and graphical models. (b) p-values of the Kolmogorov-Smirnov test, per final population size, between the data sets of the Moran and quantile function models. (c) p-values of the Kolmogorov-Smirnov test, per final population size, between the data sets of the graphical and quantile function models.

Figure 2. Statistical comparison of data set distributions. Given a final population size N, $N\in \{20000,40000,60000,80000,100000,200000,300000,400000,500000,600000,700000,800000,900000,1000000\}$, 100 simulations were generated with each of the models. For the three models, the parameters were $\mu ={10}^{-7}$, ${\mu }^{\prime }=0$, and ${\mu }_m=0$; hence only mutant cells produced by wildtype cells were allowed. The horizontal red dashed line shows the upper bound of the rejection region ($p=0.05$) for the null hypothesis $H_0$ which states that both data sets come from the same distribution. (a) p-values of the Kolmogorov-Smirnov test, per final population size, between the data sets of the Moran and graphical models. (b) p-values of the Kolmogorov-Smirnov test, per final population size, between the data sets of the Moran and quantile function models. (c) p-values of the Kolmogorov-Smirnov test, per final population size, between the data sets of the graphical and quantile function models.

3.2. Estimations reliability

By inspection, it was noted that considering more than one growth cycle provides no estimation improvement. Thus, we present the estimates and confidence intervals for the first growth cycle only, obtained by three methods, described in detail in Section 3.2: 1) bootstrap confidence interval and point estimation [46, Ch. 13] applied to the direct estimation of ${\mu }_m$ (direct estimations bootstrap), 2) the likelihood-based confidence region ([47], Eq. 14; [48]) applied to the maximum likelihood estimation of ${\mu }_m$ (MLE), and 3) the likelihood-based confidence region where the mutator strength was fixed applied to the maximum likelihood estimation of ${\mu }_m$ (Two-variable MLE).

To standardize the estimation methods, given a mutation rate $\eta \in \{\mu ,{\mu }^{\prime },{\mu }_m\}$ we defined the interval $J\in \{I,I^{\prime },I_m\}$, as

$$\begin{array}{c}\hfill J=[{log}_{10}\left(\eta \right)-2,{log}_{10}\left(\eta \right)+2].\end{array}$$

(9)

Moreover, parameter $p^{*}$ was set as 20 and since only the first growth cycle will be considered $c^{*}=1$.

Three data sets were created to test the estimation methods, with $r^{*}=125$ replicates for the data set built with the Moran model and $r^{*}=150$ replicates for the data sets built with the graphical and quantile function models. The first data set was built with the Moran model (parameters: $N=2^{33}$, $t_f=33$, $\mu ={10}^{-7}$, ${\mu }^{\prime }=5\times {10}^{-5}$, ${\mu }_m={10}^{-4}$), the second one was built with the quantile function model (parameters: $N=2^{33}$, $t_f=33$, $\mu ={10}^{-7}$, ${\mu }^{\prime }=5\times {10}^{-4}$, ${\mu }_m={10}^{-5}$), and the last one was built with the graphical model (parameters: $N=2^{33}$, $t_f=33$, $\mu ={10}^{-7}$, ${\mu }^{\prime }=5\times {10}^{-5}$, ${\mu }_m={10}^{-6}$).

For the direct estimations bootstrap method, a total of 500 estimations were performed maintaining the dilution parameters as ${\delta }_B=2^{20}/2^{33}=0.0001220703$ for dilution due to bottlenecks between growth cycles and ${\delta }_P=0.00001$ for dilution used for plating and counting mutants.

For the MLE method, a single estimation was performed setting ${\delta }_B=2^{20}/2^{33}$.

Lastly, for the two-variable MLE, a single estimation was performed setting ${\delta }_B=2^{20}/2^{33}$ and ${\mu }^{\prime }$ was taken as $k\mu$, where $k>1$ denoting the mutator strength was chosen to be $k=500$ when data sets built with the Moran model and graphical method were used to perform the estimations, and $k=500$ when the data set built with the quantile method were used to perform the estimations.

Given a data set, a performance comparison between the three methods was conducted. For comparisons related to the data set built with the Moran model, the quantile function method, and the graphical method see Table 3, Table 4, and Table 5, respectively.

Figure 4. Contour graph of the log-likelihood function and likelihood function, respectively. Given the data set built with the Moran model used for Table 3, Figure 4a shows the log-likelihood function $log\left({\mathcal{L}}_1\right)$, where the mutator strength was fixed (two-variable MLE), was constructed; while Figure 4b shows the likelihood function multiplied by the scaling factor ${\left({\displaystyle \underset{\omega \in \Omega }{max}\left\{{\mathcal{L}}_1\left(\omega \right)\right\}}\right)}^{-1}$. All parameters remain unchanged except for $p^{*}=150$, instead of $p^{*}=20$, for better image resolution. The red dot shows where the maximum of the function was attained, the blue point displays the actual mutation rates, and the region between the black curves represents the 95% confidence region.

Figure 4. Contour graph of the log-likelihood function and likelihood function, respectively. Given the data set built with the Moran model used for Table 3, Figure 4a shows the log-likelihood function $log\left({\mathcal{L}}_1\right)$, where the mutator strength was fixed (two-variable MLE), was constructed; while Figure 4b shows the likelihood function multiplied by the scaling factor ${\left({\displaystyle \underset{\omega \in \Omega }{max}\left\{{\mathcal{L}}_1\left(\omega \right)\right\}}\right)}^{-1}$. All parameters remain unchanged except for $p^{*}=150$, instead of $p^{*}=20$, for better image resolution. The red dot shows where the maximum of the function was attained, the blue point displays the actual mutation rates, and the region between the black curves represents the 95% confidence region.

Figure 5. Contour graph of the log-likelihood function and likelihood function, respectively. Given the data set built with the quantile method used for Table 4, Figure 5a shows the log-likelihood function $log\left({\mathcal{L}}_1\right)$, where the mutator strength was fixed (two-variable MLE), was constructed; while Figure 5b shows the likelihood function multiplied by the scaling factor ${\left({\displaystyle \underset{\omega \in \Omega }{max}\left\{{\mathcal{L}}_1\left(\omega \right)\right\}}\right)}^{-1}$. All parameters remain unchanged except for $p^{*}=150$, instead of $p^{*}=20$, for better image resolution. The red dot shows where the maximum of the function was attained, the blue dot displays the actual mutation rates, and the region between the black curves represents the 95% confidence region.

Figure 5. Contour graph of the log-likelihood function and likelihood function, respectively. Given the data set built with the quantile method used for Table 4, Figure 5a shows the log-likelihood function $log\left({\mathcal{L}}_1\right)$, where the mutator strength was fixed (two-variable MLE), was constructed; while Figure 5b shows the likelihood function multiplied by the scaling factor ${\left({\displaystyle \underset{\omega \in \Omega }{max}\left\{{\mathcal{L}}_1\left(\omega \right)\right\}}\right)}^{-1}$. All parameters remain unchanged except for $p^{*}=150$, instead of $p^{*}=20$, for better image resolution. The red dot shows where the maximum of the function was attained, the blue dot displays the actual mutation rates, and the region between the black curves represents the 95% confidence region.

Figure 6. Contour graph of the log-likelihood function and likelihood function, respectively. Given the data set built with the graphical method used for Table 5, Figure 6a shows the log-likelihood function $log\left({\mathcal{L}}_1\right)$, where the mutator strength was fixed (two-variable MLE), was constructed; while Figure 6b shows the likelihood function multiplied by the scaling factor ${\left({\displaystyle \underset{\omega \in \Omega }{max}\left\{{\mathcal{L}}_1\left(\omega \right)\right\}}\right)}^{-1}$. All parameters remain unchanged except for $p^{*}=150$, instead of $p^{*}=20$, for better image resolution. The red dot shows where the maximum of the function was attained, the blue dot displays the actual mutation rates, and the region between the black curves represents the 95% confidence region.

Figure 6. Contour graph of the log-likelihood function and likelihood function, respectively. Given the data set built with the graphical method used for Table 5, Figure 6a shows the log-likelihood function $log\left({\mathcal{L}}_1\right)$, where the mutator strength was fixed (two-variable MLE), was constructed; while Figure 6b shows the likelihood function multiplied by the scaling factor ${\left({\displaystyle \underset{\omega \in \Omega }{max}\left\{{\mathcal{L}}_1\left(\omega \right)\right\}}\right)}^{-1}$. All parameters remain unchanged except for $p^{*}=150$, instead of $p^{*}=20$, for better image resolution. The red dot shows where the maximum of the function was attained, the blue dot displays the actual mutation rates, and the region between the black curves represents the 95% confidence region.

4. Discussion

Our results suggest that a feasible experimental protocol – a variant of the famous Luria-Delbruck experiment – for estimating the rate of mutation to mutator phenotypes may not be out of reach. To assess this possibility, we have explored a parameter space within which we know the relevant parameters lie. There is, however, one parameter that presents us with a caveat; somewhat vexingly, this parameter is precisely the one that we are most interested in, namely, the rate of mutation to mutator. The reason for this caveat lies in a fact that we point out at the very beginning of this manuscript: aside from hand-wavy gene-counting arguments, we have no a priori knowledge of what a range of plausible values might be for the rate of mutation to mutator. This is of course not a problem for assessing our methods by comparison with simulation. It is a problem, however, when assessing the practical feasibility of our proposed protocol in the lab: if this parameter is much lower in reality than what we are assuming, the population sizes required for its estimation may become unfeasible.

One parameter whose effect we have not explored is the bottleneck size: we have only employed a bottleneck size of $2^{20}$. With this bottleneck size, we found that applying more than one growth cycle did not yield much improvement in the estimates. It may be that for a different bottleneck size, estimates can be improved by performing more than one growth cycle.

Appendix A.1. Graphical model

Suppose we grow a bacterial population inside a flask, where diet and space requirements are fulfilled for up to N cells and no cell dies. Taking N as the carrying capacity and considering a constant growth rate r associated with each individual in the bacterial population, we have a pure birth process where the population size reaches the carrying capacity in $t_f$ units of time. Here, as in many biological contexts, time is measured in generations; i.e., $t_f\in \mathbb{R}$ is such that $2^{t_f}=N$. Additionally, let’s suppose ${\mu }_m$ is the mutation rate to mutator phenotypes and ${\mu }^{\prime }$ is the mutation rate of mutator cells to the marker phenotype.

With the above assumptions, we have a pure birth process called Yule process [49], and the expected number of mutations to the mutator phenotype occurring from the beginning of the experiment until time $t_f$ is

$$\begin{array}{c}\hfill N{t}_f{\mu }_m.\end{array}$$

(A1)

Moreover, if M is a random variable modeling the number of mutations occurring in the time interval $[0,t_f]$, M follows a Poisson distribution [50], so

$$\begin{array}{c}\hfill M\sim \mathrm{Poisson}\left(N{t}_f{\mu }_m\right).\end{array}$$

(A2)

Our framework leads to the following simplification, the bacterial population grows exponentially; this is, given $t\in [0,t_f]$ then $e^{rt}$, with $r=\frac{log\left(N\right)}{t_f}$, denotes the population size at time t.

Let $f:[0,t_f]\to [0,1]$ be a function such that

$$\begin{array}{c}\hfill f\left(t\right)=\frac{e^{rt}}{N}.\end{array}$$

(A3)

Hence, $f\left(t\right)$ denotes the population size at time t compared to the final population, as a fraction.

Let $C=(T,F)$ be a random vector, where $T\sim \mathrm{Uniform}(0,t_f)$ and $F\sim \mathrm{Uniform}(0,1)$ then C represents a mutant cell that appeared at time T and belongs to the fraction F of the final population. By considering a random sample of size M from C, $C_1=(T_1,F_1),C_2=(T_2,F_2),\dots ,C_M=(T_M,F_M)$, and scattering them within $[0,t_f]\times [0,1]$, we obtain feasible and not feasible mutants since, regarding f as a probability density function, all feasible mutations lie under the graph of f and disregard them otherwise.

Figure A1 shows an example of the method above described.

Figure A1. Graphical model for selecting feasible mutations. $M\sim \mathrm{Poisson}\left(N{t}_f{\mu }_m\right)$ mutations were generated as random points $C_1,C_2,\dots ,C_M$ on $[0,t_f]\times [0,1]$. The parameter for final population size was $N=2^{33}$, for wildtype mutation rate was ${\mu }_m={10}^{-7}$, and for the final time was $t_f={log}_2\left(N\right)\approx 29.89735$. The graph of $f\left(t\right)=\frac{1}{N}exp\left(\frac{log\left(N\right)t}{t_f}\right)$ (black curve) played a discriminant role since all points below the graph of f (red points) represent actual mutations while the rest of the points (blue points) were disregarded.

Figure A1. Graphical model for selecting feasible mutations. $M\sim \mathrm{Poisson}\left(N{t}_f{\mu }_m\right)$ mutations were generated as random points $C_1,C_2,\dots ,C_M$ on $[0,t_f]\times [0,1]$. The parameter for final population size was $N=2^{33}$, for wildtype mutation rate was ${\mu }_m={10}^{-7}$, and for the final time was $t_f={log}_2\left(N\right)\approx 29.89735$. The graph of $f\left(t\right)=\frac{1}{N}exp\left(\frac{log\left(N\right)t}{t_f}\right)$ (black curve) played a discriminant role since all points below the graph of f (red points) represent actual mutations while the rest of the points (blue points) were disregarded.

Let $M_m$ be the number of actual mutations from wildtype to mutator cells. Without loss of generality, let’s suppose $C_i=(T_i,F_i)$ is the mutant cell associated with the i-th mutation, $i=1,\dots ,M_m$. Let’s note that each mutant cell $C_i$, appearing at time $T_i$, has $t_f-T_i$ units of time to give rise to its clonal linage; i.e. a subpopulation of size $X_i\ge 1$. Now, in a Yule process the clone population size that starts with $j>0$ individuals, each one with growth rate $\lambda$, has a negative binomial ($\mathrm{NB}$) distribution with parameters j and $e^{-\lambda t}$ ([49],[pp. 377–378; [51], [Eq. 3.15; [52]).

Therefore, $X_i$ is a random variable such that

$$\begin{array}{c}\hfill X_i-1\sim \mathrm{NB}\left(1,e^{-r(t_f-T_i)}\right),\forall i\in \{1,\dots ,M_m\},\end{array}$$

(A4)

since $X_i-1$ denotes the number of cells in the linage of cell $C_i$, excluding the initial cell $C_i$. This is due to the support of the negative binomial distribution we are computationally working with [53, func. rnbinom], $\{0,1,2,\dots \}$. When the support of this distribution is $\{1,2,\dots \}$, as in ([49], pp. 377–378; [51], Eq. 3.15; [52]), Equation (A4) should be rewritten as

$$X_i\sim \mathrm{NB}\left(1,e^{-r(t_f-T_i)}\right),\forall i\in \{1,\dots ,M_m\}.$$

Furthermore, in both cases, the total number of mutators in the population at time $t_f$ is

$$\begin{array}{c}\hfill X=\sum _{i=1}^{M_m}{X}_i.\end{array}$$

(A5)

To incorporate mutant cells with mutator phenotypes, let’s suppose they appear as the result of mutations within the mutator cells subpopulation. Then, mutant cells with mutator phenotypes can appear only along the growth cycle of each mutant cell $C_i$, $i=1,\dots ,M_m$.

Let $i\in \{1,\dots ,M_m\}$ be a fixed index. We know there are exactly $X_i$ mutator cells that have been raised by mutator cell $C_i=(T_i,F_i)$, and if any mutant cell with a mutator phenotype is appearing will be in the time interval $[T_i,t_f]$. Setting $X_i$ as the carrying capacity, mutant cells with mutator phenotypes have $t_f-T_i$ units of time to appear, which happens with frequency ${\mu }^{\prime }$.

Consequently, under a similar reasoning for Equation (A2), we have that the number of mutations leading to mutant cells with mutator phenotypes is a random variable $M_i$ such that

$$\begin{array}{c}\hfill M_i\sim \mathrm{Poisson}\left(X_i(t_f-T_i){\mu }^{\prime }\right).\end{array}$$

(A6)

Let $C_{i,1}=(T_{i,1},F_{i,1}),\dots ,C_{i,M_i}=(T_{i,M_i},F_{i,M_i})$, with $T_{i,j}\sim \mathrm{Uniform}(0,t_f-T_i)$ and $F_{i,j}\sim \mathrm{Uniform}(0,1)$, a random sample of size $M_i$. Hence, each $C_{i,j}$, $j\in \{1,\dots ,M_i\}$ is a random point on $[0,T_i]\times [0,1]$, representing a potential mutant cell with a mutator phenotype.

Analogous to Equation (A3), function $f_i:[0,t_f-T_i]\to [0,1]$ defined by

$$\begin{array}{c}\hfill f_i\left(t\right)=\frac{e^{r_{i}t}}{X_i},\end{array}$$

(A7)

with $r_i=\frac{log\left(X_i\right)}{t_f-T_i}$, plays the discriminant role for the $M_i$ mutants cells with mutator phenotypes. Where, as before, all points except those below the graph of $f_i$ were disregarded since they represent actual mutant cells with mutator phenotypes.

Let $M_{m,i}$ be the number of actual mutant cells with mutator phenotypes that appeared along the growth cycle of mutator cell $C_i$. Without loss of generality, let’s say $C_{i,1},\dots ,C_{i,M_{m,i}}$ represent mutant cells with mutator phenotypes then, for each $j\in \{1,\dots ,M_{m,i}\}$, the mutant cell with a mutator phenotypes $C_{i,j}$ has $t_f-(T_i+T_{i,j})$ units of time to give rise to its offspring of size $X_{i,j}$, which, according to Equation (A4), $X_{i,j}$ is a random variable such that

$$\begin{array}{c}\hfill X_{i,j}-1\sim \mathrm{NB}\left(1,e^{-r_i(t_f-T_i-T_{i,j})}\right),\end{array}$$

(A8)

for every $j\in \{1,\dots ,M_{m,i}\}$.

By repeating this process for each $i\in \{1,\dots ,M_m\}$ we have that the total number of mutators in the population at time $t_f$ is

$$\begin{array}{c}\hfill Y=\sum _{i=1}^{M_m}\sum _{j=1}^{M_{m,i}}{X}_{i,j}.\end{array}$$

(A9)

Furthermore, we can study the distribution of the numbers of mutator cells and mutant cells with mutator phenotypes, as shown in Figure A2.

For a pseudocode description of the graphical model, see Algorithm 4.

Figure A2. One thousand simulations of graphical method. (a) Histogram of the number of mutants on log scale and (b) histogram of the number of mutators on log scale. Given the parameters in Figure A1 and a mutation rate to mutator phenotype ${\mu }^{\prime }={10}^{-5}$ we simulated 1000 cultures in 354 seconds (printing plots) and in 5 seconds (printing no plots) described method. The distributions of the numbers of mutants and mutators were heavy-tailed, as expected.

Figure A2. One thousand simulations of graphical method. (a) Histogram of the number of mutants on log scale and (b) histogram of the number of mutators on log scale. Given the parameters in Figure A1 and a mutation rate to mutator phenotype ${\mu }^{\prime }={10}^{-5}$ we simulated 1000 cultures in 354 seconds (printing plots) and in 5 seconds (printing no plots) described method. The distributions of the numbers of mutants and mutators were heavy-tailed, as expected.

Appendix A.2. Quantile function model

Let’s recall that given a function $g:A\to B$ we say it is invertible if there exists a function $g^{-1}:B\to A$ (called the inverse function of g) such that

$$\begin{array}{cc}\hfill & g^{-1}\left(g\left(a\right)\right)=a,\forall a\in A,\mathrm{and}\hfill \\ \hfill & g\left(g^{-1}\left(b\right)\right)=b,\forall b\in B.\hfill \end{array}$$

Definition A1

(Quantile function). Let $F:\mathbb{R}\to [0,1]$ be a cumulative distribution function. The quantile function $F^{-}:[0,1]\to [-\infty ,\infty ]$ of F is defined by

$$\begin{array}{c}\hfill F^{-}\left(x\right)=inf\{y\in \mathbb{R}:F\left(y\right)\ge x\},\end{array}$$

(A10)

$0\le x\le 1$.

Quantile functions generalize the concept of the inverse of a function g when g is a cumulative distribution function. Now, quantile functions have a remarkable application, often called the inversion method, which is used to simulate random variables [Proposition 2 [54], a tool used in the following method.

Under assumptions and using the notation of the graphical model, given $t\in [0,t_f]$, $e^{rt}$ denotes the population size at time t. Then ${\mu }_m{e}^{rt}$ denotes the expected number of mutant cells at time t, which, as a function, can be regarded as a probability density function. Furthermore, let $k=r/(N-1){\mu }_m$ then $h:[0,t_f]\to [k{\mu }_m,k{\mu }_{m}N]$ given by

$$\begin{array}{c}\hfill h\left(t\right)=k{\mu }_m{e}^{rt}\end{array}$$

(A11)

is a probability density function, where its cumulative distribution function is $H:[0,t_f]\to [0,1]$ defined by

$$\begin{array}{c}\hfill H\left(t\right)=k{\mu }_m\frac{e^{rt}-1}{r}.\end{array}$$

(A12)

Hence its quantile function is $H^{-}:[0,1]\to [0,t_f]$ given by

$$\begin{array}{c}\hfill H^{-}\left(u\right)=\frac{1}{r}log\left(\frac{ur}{k{\mu }_m}+1\right).\end{array}$$

(A13)

So if $U\sim \mathrm{Uniform}(0,1)$, then

$$\begin{array}{cc}\hfill V& =\frac{1}{r}log\left(\frac{Ur}{k{\mu }_m}+1\right)\hfill \\ \hfill & =\frac{1}{r}log\left(U(N-1)+1\right)\hfill \end{array}$$

(A14)

is such that $V\sim H$, by the inversion method.

Applying the inversion method as above, only actual mutations leading to mutator phenotypes and their respective appearing time are provided (see Figure A3). Thus, instead of asking for the expected number of mutations, as in Eq. (A1), it is sufficient to ask for the expected number of mutators in the population, which is

$$\begin{array}{c}\hfill N{\mu }_m.\end{array}$$

(A15)

Therefore, Equation (A2) is rewritten as

$$\begin{array}{c}\hfill M\sim \mathrm{Poisson}\left(N{\mu }_m\right)\end{array}$$

(A16)

Then we can get the appearing time of M mutator cells by generating a random sample from V of size M, $V_1,\dots ,V_M$. Note that with this framework there is no need to disregard any point. Therefore, the number of actual mutations $M_m$ coincides with M, which implies Equation (A4) is rewritten as

$$\begin{array}{c}\hfill X_i-1\sim \mathrm{NB}\left(1,e^{-r(t_f-V_i)}\right),\forall i\in \{1,\dots ,M\},\end{array}$$

(A17)

and then the total number of mutators in the population at time $t_f$ is

$$\begin{array}{c}\hfill X=\sum _{i=1}^M{X}_i.\end{array}$$

(A18)

To incorporate mutant cells with mutator phenotypes we proceed as in Eq. (A15). For each $i\in \{1,\dots ,M\}$, the expected number of mutants within the i-th mutator’s offspring of size $X_i$ is $X_i{\mu }^{\prime }.$ In consequence, the random variable, $M_i$, counting the number of mutations leading to mutant cells with the i-th mutator phenotype is such that

$$\begin{array}{c}\hfill M_i\sim \mathrm{Poisson}\left(X_i{\mu }^{\prime }\right).\end{array}$$

(A19)

Now, given $i\in \{1,\dots ,M\}$ a fixed index, let $k_i=r_i/(X_i-1){\mu }^{\prime }$ then $h_i:[0,t_f-V_i]\to [k_i{\mu }^{\prime },k_i{\mu }^{\prime }{X}_i]$ defined by

$$\begin{array}{c}\hfill h_i\left(t\right)=k_i{\mu }^{\prime }{e}^{r_{i}t}\end{array}$$

(A20)

is the probability density function whose quantile function $H_i^{-}:[0,1]\to [0,t_f-V_i]$ given by

$$\begin{array}{c}\hfill H_i^{-}\left(u\right)=\frac{1}{r_i}log\left(\frac{u{r}_i}{k_i{\mu }^{\prime }}+1\right).\end{array}$$

(A21)

Thus if $U_i\sim \mathrm{Uniform}(0,1)$, then

$$\begin{array}{cc}\hfill W_i& =\frac{1}{r_i}log\left(\frac{U_i{r}_i}{k_i{\mu }^{\prime }}+1\right)\hfill \\ \hfill & =\frac{1}{r_i}log\left(U_i(X_i-1)+1\right)\hfill \end{array}$$

(A22)

is such that $W_i\sim H_i$, where ${\displaystyle H_i\left(t\right)={\int }_0^t{h}_i\left(\tau \right)d\tau }$.

Let $W_{i,1},\dots ,W_{i,M_i}$ a random sample from $W_i$, then $V_i+W_{i,j}$ denotes the appearing time of the j-th mutant cell with the i-th mutator phenotype. Therefore, analogous to Eq. (A17), the size of their offspring is a random variable $X_{i,j}$ such that

$$\begin{array}{c}\hfill X_{i,j}-1\sim \mathrm{NB}\left(1,e^{-r_i(t_f-V_i-W_{i,j})}\right),\forall j\in \{1,\dots ,M_i\}.\end{array}$$

(A23)

Implying that the number of mutant cells with mutator phenotypes is

$$\begin{array}{c}\hfill Y=\sum _{i=1}^M\sum _{j=1}^{M_i}{X}_{i,j}.\end{array}$$

(A24)

For a pseudocode description of the quantile function method, refer to Algorithm 5.

Figure A3. Comparison between simulated data and actual probability density function. Parameters were the same as in Figure A1. Data for the histogram of the appearing time of mutant cells were simulated by the inversion method and Equation (A14). The graph of $h\left(t\right)=k{\mu }_m{e}^{rt}$ (black curve) was included for reference since simulated data has distribution h.

Figure A3. Comparison between simulated data and actual probability density function. Parameters were the same as in Figure A1. Data for the histogram of the appearing time of mutant cells were simulated by the inversion method and Equation (A14). The graph of $h\left(t\right)=k{\mu }_m{e}^{rt}$ (black curve) was included for reference since simulated data has distribution h.

Figure A4. (a) Histogram of the number of mutants on log scale and (b) histogram of the number of mutators on log scale. The same parameters of Figure A2 were used to simulate 1000 cultures in two seconds using the quantile function method.

Figure A4. (a) Histogram of the number of mutants on log scale and (b) histogram of the number of mutators on log scale. The same parameters of Figure A2 were used to simulate 1000 cultures in two seconds using the quantile function method.

Appendix A.3. Description of the Moran model

We started the simulation with a single individual referred to as a normal cell (neither mutator nor antibiotic resistance). In order to simulate a pure birth process, a single individual is chosen at random from the population to reproduce. This individual divides into two equal daughter cells, and the mother cell is removed from the population.

At each birth event, each normal daughter cell may mutate to the antibiotic resistant phenotype at mutation rate $\mu$, or to the mutator phenotype at rate ${\mu }_m$. Individuals with the mutator phenotype have a mutation rate that is k times higher than the normal mutation rate. When either mutation happens, all subsequent offspring inherit the new phenotype.

We note that the mutator phenotype does not have antibiotic resistance directly, but its offspring can develop antibiotic resistance in the future. In particular, daughters of a mutator cell develop antibiotic resistance with mutation rate ${\mu }^{\prime }=k\mu$, while daughters of an antibiotic resistant cell become mutators at rate ${\mu }_m$. Consequently, the population consists of four types of individuals: normal cells, antibiotic resistant cells, mutator cells, and individuals exhibiting both antibiotic resistance and the mutator phenotype.

In summary, the simulation proceeds by choosing an individual at random to reproduce, removing that individual, and determining which types of daughters are added in its place, as described above. This process is repeated until the population size reaches the desired final size (e.g. $2^{33}$). We can also simulate repeated growth phases and population bottlenecks by randomly sampling the final population to achieve a desired initial population size (e.g. $2^{20})$, and repeating the procedure above. See Algorithm 1, 2, and 3 for a detailed pseudocode description.

Algorithm 1: Moran model definitions

Input: $\mu ={10}^{-7}$: Mutation rate from wildtype to antibiotic resistance. ${\mu }_m={10}^{-6}$: Mutation rate from wildtype to mutator. ${\mu }^{\prime }={10}^{-5}$: Mutation rate from mutator to antibiotic resistance. $N=2^{33}$: Final population size. Ntotal: Current population size, total of all types. NAMB[0], NAMB[1], NAMB[2], NAMB[3] : Numbers of each of four types of individuals including normal cells, antibiotic-resistant cells, mutator cells, and individuals including both antibiotic-resistant and mutator phenotype respectively. nreps=1000: Number of replicate lines. cumsum: Cumulative frequency of each type in the population, set cumsum[3]=1. pinit = [0,20,20,20]: (see below) Ninit: initial population size for each growth phase, $=2^{pinit}$. bottle: Bottleneck sampling fraction. r, $r_2$, $r_3$: Random numbers. ngrowths = 4: number of growth phases.

Algorithm 2: Moran model pseudocode

Algorithm 3: Moran model pseudocode, continued.

Algorithm 4: Graphical model pseudocode.

Algorithm 5: Quantile function model pseudocode.

Appendix B.1. Direct estimations

Let X be a discrete random variable with support S and probability mass function $p\left(x\right)=\mathbb{P}(X=x)$, for all $x\in S$. Let’s denote by R the convergence range of the empirical probability-generating function of X

$$\begin{array}{c}\hfill G_X\left(\theta \right)=\mathbb{E}\left({\theta }^X\right)=\sum _{x\in S}p\left(x\right){\theta }^x.\end{array}$$

It is known that $R_0:=\{\theta \in \mathbb{C}:|\theta |\le 1\}\subseteq R$ because for any arbitrary and fixed complex number $\theta$ such that $\left|\theta \right|\le 1$, it holds that

$$\begin{array}{c}\hfill \sum _{x\in S}\left|p\left(x\right){\theta }^x\right|\le \sum _{x\in S}p\left(x\right)=1.\end{array}$$

Implying $G_X$ converges absolutely on $R_0$, which proves the statement.

However, when working with the Luria-Delbrück probability generating function, it is customary to consider the set $\{\theta \in \mathbb{R}:0\le \theta \le 1\}$ exclusively; i.e., the interval $[0,1]$. This is primarily due to biological considerations. Values close to zero, approaching from the right side, reflect the effect of observing small quantities of mutant cells, while values close to one, approaching from the left side, reflect the influence of large quantities of mutant cells in cell cultures [64].

In works such as [45] and [64], the study of the probability-generating function has been restricted to the points $\theta \in \Theta$, where $\Theta =\{0.1,0.8,0.9\}$, to generate estimates of $\mu$. However, to estimate ${\mu }_m$, we consider the points $\theta \in {\Theta }^{*}$, where ${\Theta }^{*}=\{0.4,0.5,0.6,0.7,0.8,0.9\}$.

This choice is made because we will be using the probability generating function for the sum of the random variable counting the number of mutant cells produced by wildtype cells and the one counting the number of mutant cells produced by mutator cells. Since we expect to observe small quantities of mutant cells sporadically. Hence, considering the set ${\Theta }^{*}$ allows us to reflect this framework.

Let’s assume that $r^{*}$ replicates of the experiment have been considered, each with $c^{*}$ growth cycles. Let $m_i\in {\mathbb{R}}^{r^{*}}$, where $i=1,\dots ,c^{*}$, be a vector whose j-th entry, $j=1,\dots ,r^{*}$, represents the sum of the number of mutant cells produced by wildtype cells and the number of mutant cells produced by mutator cells.

To incorporate this information into the dilution process of the experimental protocol, let’s define the vector $n_i\in {\mathbb{R}}^{r^{*}}$, where $i=1,\dots ,c^{*}$, such that its j-th entry, $j=1,\dots ,r^{*}$, corresponds to a realization of the random variable $S_j$

$$\begin{array}{c}\hfill S_j\sim \mathrm{Binomial}(m_{i,j},{\delta }_i),\end{array}$$

where $m_{i,j}$ and ${\delta }_i\in (0,1)$ represent the j-th entry of the vector $m_i$ and the dilution factor considered for the i-th growth cycle, respectively.

Before introducing the estimation method, we propose intervals $I=[l,u]$, $I^{\prime }=[l^{\prime },u^{\prime }]$ and $I_m=[l_m,u_m]$ such that

$$\begin{array}{cc}\hfill {log}_{10}\left(\mu \right)& \in I,\hfill \\ \hfill {log}_{10}\left({\mu }^{\prime }\right)& \in I^{\prime },\mathrm{and}\hfill \\ \hfill {log}_{10}\left({\mu }_m\right)& \in I_m.\hfill \end{array}$$

Given that, for computational purposes, the approach above provides greater accuracy in the estimation process.

To work with a finite number of possibilities, the intervals will be partitioned to consider only ${(p^{*}+1)}^3$ points within the set $\Omega :=I\times I^{\prime }\times I_m$. Thus, instead of working with the interval $[a,b]$, where $(a,b)\in \{(l,u),(l^{\prime },u^{\prime }),(l_m,u_m)\}$, we will work with the set of points

$$\begin{array}{c}\hfill \left\{a{\left(\frac{b}{a}\right)}^{q/p^{*}}:q=0,1,\dots ,p^{*}\right\}\subset [a,b].\end{array}$$

(A25)

Let $F_1:{\Theta }^{*}\times \Omega \to \mathbb{R}$ be defined as

$$\begin{array}{c}\hfill F_1(\theta ,{\omega }_1,{\omega }_2,{\omega }_3)=h{\left(\theta \right)}^{{\omega }_1}h{\left(h{\left(\theta \right)}^{{\omega }_2}\right)}^{{\omega }_3},\end{array}$$

(A26)

and let $F_i:{\Theta }^{*}\times \Omega \to \mathbb{R}$ be given by

$$\begin{array}{c}\hfill F_i(\theta ,{\omega }_1,{\omega }_2,{\omega }_3)=h{\left(F_{i-1}\left(\theta \right)\right)}^{{\omega }_1}h{\left(h{\left(F_{i-1}\left(\theta \right)\right)}^{{\omega }_2}\right)}^{{\omega }_3},\end{array}$$

(A27)

for $i=2,\dots ,c^{*}$.

In this way, the estimators ${\widehat{\Delta }}_i$, ${\widehat{{\mu }^{\prime }}}_i$, and ${\widehat{{\Delta }_m}}_i$ for $\Delta$, ${\mu }^{\prime }$, and ${\Delta }_m$ in the i-th growth cycle, respectively, are those that satisfy

$$\begin{array}{c}\hfill \sum _{\theta \in {\Theta }^{*}}{\left(F_i(\theta ,{\widehat{\Delta }}_i,{\widehat{{\mu }^{\prime }}}_i,{\widehat{{\Delta }_m}}_i)-G_i\left(\theta \right)\right)}^2=\underset{({\omega }_1,{\omega }_2,{\omega }_3)\in \Omega }{min}\sum _{\theta \in {\Theta }^{*}}{\left(F_i(\theta ,N{10}^{{\omega }_1},{10}^{{\omega }_2},N{10}^{{\omega }_3})-G_i\left(\theta \right)\right)}^2,\end{array}$$

for $i=1,\dots ,c^{*}$.

Appendix B.2. Maximum likelihood estimations

Applying maximum likelihood methods requires working directly with probability mass functions. In this case, it would involve dealing with the explicit expression of the Luria-Delbrück distribution, which admits no explicit analytic expression [47].

However, by considering the direct approach, outlined up to Equation (A27), sufficient information is provided to recover the probability mass states of the Luria-Delbrück distribution and implement maximum likelihood methods. This can be done by employing the fast Fourier transform as in [65] and [66, Eq. 5].

The Luria-Delbruck distribution, in its simplest form (i.e., the Lea-Coulson formulation that we employ here), has no finite moments, hence its moment-generating function is not analytic. Its characteristic function, on the other hand is analytic, and it is $F\left(e^{\mathtt{i}\theta }\right)$. We leverage this fact to obtain the probability mass function defined by the inverse Fourier transform of $F\left(e^{\mathtt{i}\theta }\right)$.

Let ${\displaystyle d_i=\underset{j\in \{1,\dots ,r^{*}\}}{max}{n}_{i,j}}$, where $n_{i,j}$ is the j-th entry of the vector $n_i$, for $i=1,\dots ,r^{*}$. In other words, $d_i$ is the maximum number of mutant cells present in the i-th growth cycle after dilution.

In this way, given $\omega =({\omega }_1,{\omega }_2,{\omega }_3)\in \Omega$, and defining ${\mathcal{F}}_i(\theta ;\omega ):=F_i(\theta ,N{10}^{{\omega }_1},{10}^{{\omega }_2},N{10}^{{\omega }_3})$, the real part of the fast Fourier transform of ${\mathcal{F}}_i$

$$\begin{array}{c}\hfill {\check{\mathcal{F}}}_i\left(k\right)=\frac{1}{d_i}\sum _{l=0}^{d_i-1}{\mathcal{F}}_i\left(\frac{2\pi \mathtt{i}l}{d_i};\omega \right)exp\left(-\frac{2\pi \mathtt{i}lk}{d_i}\right),\end{array}$$

(A28)

where ${\mathtt{i}}^2=-1$, approximates the value of

$$\begin{array}{c}\hfill p_{X_i}(k;\omega )=P(X_i=k;\omega ),\end{array}$$

(A29)

for $k=0,1,\dots ,d_i$. Where $X_i$ is the random variable whose probability-generating function is ${\mathcal{F}}_i$. In other words, Eq. (A28) provides an approximation to the probability of observing k mutants on the i-th growth cycle after the dilution.

Therefore,

$$\begin{array}{c}\hfill {\mathcal{L}}_i\left(\omega \right)=\prod _{j=1}^{r^{*}}{p}_{X_i}(n_{i,j};\omega )\end{array}$$

denotes the likelihood function associated with the i-th growth cycle.

This implies that the maximum likelihood estimators ${\widehat{\Delta }}_i^{*}$, ${\widehat{{\mu }^{\prime }}}_i^{*}$, and ${\widehat{{\Delta }_m}}_i^{*}$ for $\Delta$, ${\mu }^{\prime }$, and ${\Delta }_m$ in the i-th growth cycle, respectively, are those satisfying

$$\begin{array}{c}\hfill {\mathcal{L}}_i({\widehat{\Delta }}_i^{*},{\widehat{{\mu }^{\prime }}}_i^{*},{\widehat{{\Delta }_m}}_i^{*})=\underset{\omega \in \Omega }{max}{\mathcal{L}}_i\left(\omega \right).\end{array}$$

(A30)