Determination of Mutation Rates with Two Symmetric and Asymmetric Mutation Types

Abstract

We revisit our earlier paper, with two of the coauthors, in which we proposed an unbiased and consistent estimator ${\widehat{\mu }}_n$ for an unknown mutation rate $\mu$ of microorganisms. Previously, we proved that the associated sequence of estimators $\left\{{\widehat{\mu }}_n\right\}$ converges to $\mu$ almost surely pointwise on a nonextinct set ${\Omega }_0$. Here, we show that this sequence converges also in the mean square with respect to conditional probability measure ${\mathbb{P}}_0\left(·\right)=\mathbb{P}\left(·\cap {\Omega }_0\right)/\mathbb{P}\left({\Omega }_0\right)$ and that, with respect to ${\mathbb{P}}_0,$ the estimator is asymptotically unbiased. We further assume that a microorganism can mutate or turn to a different variant of one of the two types. In particular, it can mean that bacteria under attack by a virus or chemical agent are either perishing or surviving, turning them to stronger variant. We propose estimators for their respective types and show that they are a.s. pointwise and $L^2$-consistent and asymptotically unbiased with respect to measure ${\mathbb{P}}_0$.

1. Introduction

In numerous biological applications related to the proliferation of microorganisms, such as bacteria, viruses, and even cancer cells, at some point of their division, they start mutating. New mutants are typically stronger and more virulent. In most common situations, bacteria divide in two forming a deterministic branching process. If a colony starts from a single organism, in the nth generation, its population size is $2^n$ (unless some progenies perish). Suppose at some point that one of its offspring mutates, say, with probability $\mu$ (which is typically small). Then, a new mutant divides like a nonmutant and it engenders a deterministic branching process. As the rate of reverse mutation on the same phenotype, called revertant, is extremely small [1,2], we do not consider this option in our paper. The mutation probability $\mu$ is called the mutation rate.

1.1. Our Approach Compared to Other Methods

The mutation rate is a common target in biomedical research and its applications. Ever since the seminal paper by Luria and Delbrück [3] in 1943, numerous works on determination of mutation rate have emerged (cf., [4,5,6,7,8,9,10,11,12,13,14]). They all have dealt with various estimates for $\mu$.

The deterministic branching process (a rooted tree) contained a branching subprocess (random subtree) of non mutants that turned out to be a Galton–Watson process, although the target process was that of new mutants. The new mutants and pre-existing mutants form non-Markov processes.

Why do we deem the branching processes approach more rigorous? In Luria and Delbrück [3], the authors assumed that the number of mutations $M\left(\tau \right)$ that occur in interval $\left(t,t+\tau \right)$ equals $\mu \tau N_t$, where $N_t$ is the total number of bacteria at time t. This assumption is inaccurate, because to determine the mutation rate $\mu$, they had to exclude from $N_t$ the number of pre-existing mutants, even if this number is small. We actually prove that even if the mutation rate $\mu$ is small, the number of pre-existing mutants can converge to infinity and thus cannot be neglected. Moreover, $M\left(\tau \right)$ would be better off to equal some $\mu \tau 2X_t$, where $X_t$ is the number of nonmutants at time t, which doubles in interval $\left(t,t+\tau \right)$ unless $\tau$ is too small. Stewart [15] also states that the Luria and Delbrück method is inaccurate.

Furthermore, if that would not be the case and $N_t$ would be approximately equal to the number of nonmutants at time t, then the Luria and Delbrück method is still problematic. Here is why: suppose $X_n$ is the number of nonmutants in generation n and suppose that each nonmutant can turn to a mutant with equal probability $\mu$. Thus, each nonmutant is a Bernoulli r.v. with parameter $\mu$. The number of newly born mutants $N{M}_n$ in generation n is thus binomial with parameters $\left(X_n,\mu \right)$, that is, $N{M}_n\in \left[\mathbb{B}\left(X_n,\mu \right)\right]$. Because $X_n$ is large and $\mu$ is small, the number $N{M}_n$ is approximately Poisson, with parameter ${\lambda }_n=N{M}_n\mu$, which carries unknown $N{M}_n$ and $\mu$. The probability that there is no mutant at generation n would have been $p_0=e^{-{\lambda }_n}$, which the authors claimed could be determined experimentally and is fairly simple. However, there is still a problem of finding $N{M}_n$, which is not equal to the total number cells $N_n$, nor does it equal the number of all mutants $M_n$ at generation n. This is the reason why we think this method is unpredictably inaccurate (although Zheng in [16,17] claims that the Luria/Delbrück method is the best). Precisely, $\mu \approx N{M}_n/2X_{n-1}$ (≠$M_n/N_n=M_n/2^n$ as the authors allude), provided that the latter is consistent in some sense and that n is sufficiently large. None of these requirements has ever been met in their work, nor in the work of their followers. This problem was also pointed out by Stewart [15].

Another problem in the statistical evaluation of the popular $p_0$-method is that in order to maintain a “nonmutant” plate (which is required to reasonably estimate $p_0$), it applies to counts rendered at early times or for cultures with extremely low mutation rates to warrant at least one plate with no mutants. However, this carries an adverse effect on statistical significance.

In various other alternatives, many formulas contain the ratio $M_n/N_n$ as a major component in their estimators (rather $M_{n+k}/N_{n+k}-M_n/N_n$) which has a similar shortcoming and it cannot be salvaged by logarithmic or other multipliers often occurring in the literature. We noticed that no work offers a general formula for the determination of $\mu$ (whether practical or not). All heuristic methods may work in some special cases, but none are claimed to be general. In the literature on mutation rates, we see that some methods claimed to be efficient when following certain heuristic instructions imposed on $\mu$ (which was to be determined) and the mean number of mutants. It is not clear how the results obtained for $\mu$ were validated. We compared several tables of the results produced from different methods and they all vary.

In Niccum et al. [18], we offered the first rigorous method of determining $\mu$ by a pointwise rapidly convergent sequence $\left\{{\widehat{\mu }}_n\right\}$ of estimators. At first, the original formula for ${\widehat{\mu }}_n$ we operated on looked impractical because it contained newly formed mutants that are difficult to tell from pre-existing mutants, but it was modified by replacing it with an identical formula requiring just the knowledge of the mutants in any two consecutive generations that turned out to be practical in real labs. That “it was still difficult to observe” as stated in Zheng [17] was not the situation with a series of our lab experiments [18] performed by us and our students (who coauthored the paper) and beyond. It also needs to be noted that the proposed estimator had a rapid speed of convergence that has been validated in numerous subsequent simulations rendered by our students in senior undergraduate and graduate projects with no single exception.

1.2. Synchronized Cultures

In our original mathematical model [18], we assumed that bacteria replicate synchronously. We retain this mathematical assumption in the present work as well. As explained in [18], the real count in the lab is a little different. This method has been rendered for decades and the count of mutants (that is compatible with our theoretical assumption) agrees with quite a few experiments made by different researchers, which obviously renders coincidences (and thus inaccuracies) probabilistically impossible.

Yet for the skeptics of our mathematical methodology, we would like to mention the common practice of synchronization of bacterial cultures known in the biological literature. First, a synchronized culture is a culture that contains cells in the same growth stage when they replicate, as opposed to asynchronous growth when the replication cycles are random. Most known cell divisions in a population are asynchronous. It is an observed, although unproven, assumption that random deviations of the cycle lengths are symmetric around the mean time of division. However, that needs to be rigorously validated. If this holds true, it is easy to deduce that a mere observations of a bacterial culture would give a pretty accurate identification of an underlying generation.

Alternatively, there are widely used synchronization approaches from 1957 (cf. Campbell [19]) and earlier and these continue to be described to the present time (cf. Chang et al. [20] of 2019). A popular technique, known as Helmstetter–Cummings [21,22] of 1968 and 2015 deals with an unsynchronized culture filtered through cellulose nitrate membrane filter. There are numerous other synchronization techniques by Anderson and Pettijohn [23], Kepes and Kepes [24,25], Kubitcheck [26], Ling and Chang [27], Noack, Klöden, and Bley [28], Shehata and Marr [29], and Wallden et al. [30], to name a few. A rigorous description of all those methods lies outside the scope of this paper.

Note that experimentally cells may not remain synchronized during growth complicating analysis. As mentioned above, it would be interesting to conduct lab experiments and learn about the statistical impact of unsynchronized replication. It seems that even if we assume the replication to be synchronous while in facts it is not, the deviations will be symmetric and not large, and as the result, the incorrect assumption will not corrupt the counts, because so far, the observed numbers were not conflicting with our experiments, and they agreed with experiments conducted in other labs. This conjecture is by no means intended to replace a rigorous investigation, which we plan to render.

1.3. Mutation or No Mutation?

There are many anecdotal stories about antibiotics causing antibiotic resistance and breeding new forms of virulent bacteria, in quite a few cases offering little to no remedy. One such notorious form is called MRSA, which stands for methicillin-resistant Staphylococcus aureus (S. aureus), a type of bacteria that has become resistant to potent antibiotics used to treat ordinary staph infections. It thus means that it is tougher to treat compared to other strains of Staphylococcus.

MRSA most commonly causes skin infections and can sometimes trigger pneumonia and other conditions. The symptoms largely depend on which part of the body is affected. Untreated, MRSA can become severe, leading to sepsis. MRSA infections are usually seen in individuals who were previously or currently in hospitals or other health care facilities such as dialysis centers or nursing homes. Many professionals are alarmed by the spread of MRSA as due to its resistance to many types of drugs; it is becoming known as a “super bug”.

On the other hand, there is a general consensus that while bacteria mutate, this process is neither impacted nor exacerbated by the use of antibiotics or other chemical agents. They claim that mutations occur periodically without any external force. However, that natural selection allows these to thrive over other bugs that do not have resistance. (It must be noted that antibiotics are natural products of other organisms that also try to keep bacteria away.)

Before we elaborate it, we wonder what a mutation is referred to in cellular/molecular biology. A mutation is formally defined as an alteration to the gene sequence of an organism. A “natural mutation” occurs spontaneously and can happen during bacterial DNA replication where enzymes choose a wrong base to pair with the original DNA, or due to the environmental UV light, modifying the structure of the original DNA strand that allows it to pair with bases that are not normally compatible.

While the entire bacterial genome is susceptible to mutation, certain mutations are more easily detected than others because of selection. An example is the trpE gene in the Escherichia coli (E. coli) WP2 strain [31]. This high mutation region allows researchers to detect mutation on the strain using the so-called reverse mutation assay, in which the bacteria are grown in an amino-acid-deficient culture to test if the mutation occurred ($tr{p}^{+}$) or not ($tr{p}^{-}$) [32].

The natural mutation rate for E. coli in a rich medium mutation rate is ${10}^{-3}$ per genome per generation [33]. The mutation rate of methicillin-resistant Staphylococcus aureus (S. aureus)—MRSA strain is similar to the reported mutation rate of wild-type E. coli [34]. The elevated mutation rate is suspected to cause the increased rate of bacteria vancomycin-resistant genes, especially for the vancomycin-intermediate-resistant S. aureus strains—VRSA [35].

1.4. Other Means of Genetic Change

Other means of genetic change exist, in addition to classical spontaneous or genotoxicant-induced mutation, and, in some cases, these resulting genetic changes may also be covered by our model in as much as they may represent stochastic events. While the mechanisms responsible for these changes are well understood and these processes are not considered as classical mutations, they do result in genetic change in the affected organism. In some cases, the probabilities of these events are rare and, as such, the resulting genetic change may be described by probabilistic models. For the purposes of a more comprehensive discussion of other modes genetic changes, these alternative modes of change will be briefly described.

Transformation. Transformation is a process in which a bacteria takes up exogenous (free) DNA from its surrounding environment. This DNA can come from dead (lysed) bacteria. Transformation can happen artificially (such as a heat shock/electroporation), or naturally. A natural transformation occurs when there is a cooperative expression of multiple genes under specific conditions called competence [36]. This state is induced by a limited nutrition (usually amino-acid-deficient environment), associated with the stationary phase. For some other bacteria, transformation occurs most efficiently at the end of exponential growth approaching the stationary phase.

Transformation is one process that can lead to Horizontal Gene Transfer (HGT). There is a hypothesis about the HGT facilitated by the transformation of antibiotic-induced cell-wall-deficient bacteria [37]. This hypothesis states that before the bacteria die, they generate a cell-wall-deficient form and increase the uptake rate of exogenous DNA, because the DNA from dying bacteria release DNA. They are not doing the uptake as they are dead and could not then restore life to propagate this. Combining with the fact that patients who were treated with diverse antibiotics are more prone to get infected with MRSA makes this hypothesis more credible.

In this process, two bacterial cells attach, and plasmids and free intracellular DNA fragments travel from the donor cell to the target cell. In the target cell, plasmids can begin to replicate, whereas DNA fragments may become incorporated in the target cell’s genome. A study on MRSA conjugation stated that this mechanism could be one of the pathways by which methicillin resistance is transmitted among S. aureus strains [38].

Transduction. Another mechanism for bacteria genome variation is called transduction. When bacteria are infected with a virus (usually a bacteriophage), the virus uses the cell’s materials to produce new viruses and kill the bacteria when it reaches some threshold. During the production of new viruses, some fragments of the original bacteria’s DNA are taken up into the virus’ DNA. When this virus infects other cells, it also introduces the DNA to the new cell. If the infection goes to the latent phase and the fragments contain antibiotic genes, the bacteria become antibiotic-resistant [39]. Virus latency (or viral latency) is the ability of a pathogenic virus to lie dormant (latent) within a cell, denoted as the lysogenic part of the viral life cycle. If virus becomes latent, it neither reproduces nor kills the host cell.

So, as we see it, there is a variety of genetic alterations to the bacteria that are technically not mutations, but they often produce virulent species. While it is a common conjecture that bacteria do not mutate under chemical agents such as antibiotics, the latter does cause formidable DNA changes turning some strains to dangerous forms that can be even life-threatening. Most of these transformations are detectable and they can be observed and analyzed mathematically to determine the odds of their changes whether or not we call them mutations. Therefore, it stands to reason to consider those phenomena and place them in the same category as mutations and apply the same mathematical tools to establish rigorous estimators of their genetic transformations. It is also possible that an epigenetic change can occur that may be detectable by selection.

1.5. The Results

In the present paper, we show that the estimator proposed in [18] is not only almost surely pointwise consistent but also consistent in the $L^2$-norm giving us yet another mathematical confirmation of the goodness of the proposed estimator. We also prove the convergence with respect to other probability measures. We further assume that a microorganism can mutate in one of the two types, which in particular, can model the spread of virulent antibiotics-resistant bacteria (referred to in the media as “superbugs”) that evolve under attacks by chemical agents so that they either die or mutate to survive the attacks. Consequently, mutated or altered bacteria’s variants that survive are referred to as type 1 mutants and other mutants are called type 2 mutants. We also propose estimators for either type of mutants. These estimators turn out to be unbiased and consistent in every common sense, that is, almost surely pointwise and in the $L^2$-norm. (Some other forms of consistency are discussed throughout the forthcoming sections.)

We note that while death of bacteria is present in most bacterial cultures, it is rarely used in mathematical modeling. We address this matter in the present article and even generalize our model allowing replicated offspring to turn to one of the two mutant types with probabilities p and $1-p$, respectively. It makes sense to call such division symmetric when $p=\frac{1}{2}$. In the event that this represents bacterial death, that, as mentioned, takes place under the use of a chemical agent (such as antibiotic) or bacteriophage attack.

2. Types of Stochastic Convergence

To make our paper self-contained, we present a short background on types of convergence of random variables (r.v.’s) related to the consistency of our estimators. (Cf. Dshalalow [40].)

Definition 1

(Types of convergence for a sequence of r.v.’s). Let $Z,Z_1,Z_2,\dots$ be a sequence of r.v. on a probability space $\left(\Omega ,\mathcal{F},\mathbb{P}\right)$. We say that the sequence $\left\{Z_n\right\}$ converges to r.v. Z

(i)

in probability$\left(\mathit{in}\mathit{notation}{Z}_n\stackrel{\mathbb{P}}{\to }Z\right)$if

$$\underset{n\to \infty }{lim}\mathbb{P}\left\{\left|Z_n-Z\right|>\epsilon \right\}=0,\mathit{for}\mathit{each}\epsilon >0,$$

(ii)

in the p-norm (or in the $L^p$-norm or in the p-th mean)$\left(\mathit{in}\mathit{notation}{Z}_n\stackrel{L^p}{\to }Z\right)$if

$$\underset{n\to \infty }{lim}\mathbb{E}\left[{\left(Z_n-Z\right)}^p\right]=0,p\ge 1,$$

(iii)

in the mean (or in the$L^1$-norm)$\left(\mathit{in}\mathit{notation}{Z}_n\stackrel{L^1}{\to }Z\right)$if

$$\underset{n\to \infty }{lim}\mathbb{E}[Z_n-Z]=0,$$

(iv)

almost surely (a.s.) (or$\mathbb{P}$-almost surely)$\left(\mathit{in}\mathit{notation}{Z}_n\stackrel{a.s.}{\to }Z\right)$if

$$\mathbb{P}\left\{\omega \in \Omega :\underset{n\to \infty }{lim}\left|Z_n\left(\omega \right)-Z\left(\omega \right)\right|=0\right\}=1.$$

The latter means that$Z_n$converges to Z pointwise for almost all$\omega \in \Omega$.

Theorem 1.

The following relations hold true:

$$Z_n\stackrel{a.s.}{\to }Z\Rightarrow Z_n\stackrel{\mathbb{P}}{\to }Z\mathrm{and}{Z}_n\stackrel{L^2}{\to }Z\Rightarrow Z_n\stackrel{L^1}{\to }Z\Rightarrow Z_n\stackrel{\mathbb{P}}{\to }Z$$

Definition 2.

If${\widehat{\theta }}_n$is an estimator of some parameter θ, then${\widehat{\theta }}_n$is a statistic (Borel measurable function) of a random sample$X_1,\dots ,X_n\in \left[X\right]$. In general, suppose$X_1,X_2,\dots$is a sequence of r.v.’s on a filtered probability space$\left(\Omega ,\mathcal{F},\left({\mathcal{F}}_n\right),\mathbb{P}\right)$adapted to$\left({\mathcal{F}}_n\right)$. Notice that the r.v.’s$X_1,X_2,\dots$need not be independent nor identically distributed. Let

$${\widehat{\theta }}_n:=\delta \left(X_1,\dots ,X_n\right)$$

be a statistic of the sample of the first n r.v.’s. Then,$\left\{{\widehat{\theta }}_n\right\}$is a sequence of estimators associated with$\left\{X_n\right\}$and induced by a Borel measurable function δ. Let θ be a real number referred to as a parameter of$\left\{X_n\right\}$.

(i)

The sequence$\left\{{\widehat{\theta }}_n\right\}$is unbiased if$\mathbb{E}{\widehat{\theta }}_n=\theta$holds for all n.

(ii)

The sequence $\left\{{\widehat{\theta }}_n\right\}$ is asymptotically unbiased if ${lim}_{n\to \infty }\mathbb{E}{\widehat{\theta }}_n=\theta$

(iii)

The sequence$\left\{{\widehat{\theta }}_n\right\}$is consistent

(a)

in probability if${\widehat{\theta }}_n\stackrel{\mathbb{P}}{\to }\theta$

(b)

in the pth mean (or in the$L^p$-norm) if${\widehat{\theta }}_n\stackrel{L^p}{\to }\theta$for some$p\ge 1$

(c)

almost surely if${\widehat{\theta }}_n\stackrel{a.s.}{\to }\theta$

We say that ${\widehat{\theta }}_n$ is a consistent estimator in probability, pth mean, or almost surely, respectively, if so is the sequence $\left\{{\widehat{\theta }}_n\right\}$.

3. A Background on Stochastic Estimators of Mutation Rate

Consider a population of microorganisms, such as bacteria, stemming from a single parent. Assume that the bacteria replicate by the division in exactly two progeny. At generation n, there are exactly $2^n$ bacteria in the population forming a deterministic branching process, provided that the bacteria do not die. Suppose that a generation began with a single parent bacterium (at generation zero) and beginning from generation 1, each bacteria, independently of the others, mutates with probability $\mu \in \left(0,\frac{1}{2}\right)$. The parent bacterium of the generation is thus assumed to be a nonmutant organism. Let $X_n$ be the total number of nonmutants at generation n, $N{M}_n$ be the total number of new mutants at generation n, and $M_n$ be the total number of all (preexisting and new) mutants at generation n. We assume that $\mu$ (referred to as the mutation rate) is constant at every generation. To determine $\mu$ experimentally, we proposed in [18] the estimator

$${\widehat{\mu }}_{n+1}=N{M}_{n+1}/2X_n=\left(2X_n-X_{n+1}\right)/2X_n=1-\frac{X_{n+1}}{2X_n}$$

(1)

which turned out to be unbiased and consistent a.s. on the non extinction set ${\Omega }_0$ (defined below) such that $\mathbb{P}\left({\Omega }_0\right)=1-{\mu }^2/{\nu }^2$, where $\nu =1-\mu$. (For various bacteria, such as E. coli, this probability is extremely close to 1.)

Denote

$$C_n=\left\{\omega :X_n\left(\omega \right)>0\right\}\mathrm{and}{\Omega }_0=\bigcap _{n=1}^{\infty }{C}_n.$$

(2)

Since the sequence $\left\{C_n\right\}$ is monotone decreasing, the notation ${\Omega }_0$ in (2) makes sense and ${\Omega }_0$ is referred to as the nonextinction set. As mentioned above, $\mathbb{P}\left({\Omega }_0\right)=1-\frac{{\mu }^2}{{\nu }^2}=1-\varphi$.

In light of our next applications, it is more convenient to turn to the confined space $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right),$ where

$${\mathbb{P}}_0\left(·\right)=\mathbb{P}\left(·\cap {\Omega }_0\right)/\mathbb{P}\left({\Omega }_0\right)=\mathbb{P}\left(·|{\Omega }_0\right).$$

(3)

If $X:{\Omega }_0\to \overline{\mathbb{R}}$ is a r.v. on the confined space $\left({\Omega }_0,\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$, then (cf., Dshalalow [41]), the expected value ${\mathbb{E}}_0$ (with respect to measure ${\mathbb{P}}_0$) is

$${\mathbb{E}}_{0}X=\int Xd{\mathbb{P}}_0=\frac{1}{\mathbb{P}\left({\Omega }_0\right)}\underset{{\Omega }_0}{\int }Xd\mathbb{P}=\frac{1}{\mathbb{P}\left({\Omega }_0\right)}\mathbb{E}\left[{\mathbf{1}}_{{\Omega }_0}X\right].$$

(4)

4. Preliminaries

Let $W_n=\frac{X_n}{{\left(2\nu \right)}^n},n=0,1,\dots$, be the normalized sequence of the respective nonmutants populations. This sequence is a martingale that almost surely on ${\Omega }_0$ converges to a r.v. $W_{\infty }$ (cf. [42]). Furthermore, Theorem 2 [Section 6, Chapter I] of [42] reveals the nature of this r.v., also stating that the convergence of $W_n$ to $W_{\infty }$ is also in the mean square.

Note that Theorem 2 below is borrowed from [42] originally formulated and proved as Theorem 2. It could be more rigorously crafted under the specification of the confined probability space $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$ introduced in (3) and (4). We proceed with those key assumptions throughout the rest of the paper.

Theorem 2

([42]). If $\nu >\frac{1}{2},$ then

(i)

${∥W_n-W_{\infty }∥}_{L^2}\to 0$.

(ii)

$\mathbb{E}{W}_{\infty }=1,\mathrm{Var}{W}_{\infty }=\frac{1}{2\nu }·\frac{{\sigma }^2}{2\nu -1}=\frac{\mu }{\nu -\mu }$.

(iii)

$\mathbb{P}\left\{W_{\infty }=0\right\}=\varphi =1-\mathbb{P}\left({\Omega }_0\right)$.

From Theorem 2 ($iii$), we have that ${\mathbb{P}}_0\left\{W_{\infty }>0\right\}=1$.

Assume that $\nu >\frac{1}{2}$. Since $W_n\to W_{\infty }:\Omega \to \left\{0\right\}\cup (0,\infty ]$ with the distribution $\varphi$ and $1-\varphi$, respectively, the sequence $\left\{X_n\right\}$ converges to

$$X_{\infty }=\underset{n\to \infty }{lim}{W}_n{\left(2\nu \right)}^n$$

where $X_{\infty }$ is defined on $\Omega$ and valued in $\left\{0\right\}\cup \left\{\infty \right\}$, with the distribution $\varphi$ and $1-\varphi =\mathbb{P}\left({\Omega }_0\right)$, respectively. Therefore, on the confined probability space, $\left({\Omega }_0,\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$, $X_{\infty }$ equals ∞${\mathbb{P}}_0$-a.s. Hence, the r.v. $\frac{1}{X_{\infty }}$ on $\left({\Omega }_0,\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$ equals $0,{\mathbb{P}}_0$-a.s., which implies that

$${\mathbb{E}}_0\left[\frac{1}{X_{\infty }}\right]=0$$

(5)

Now, because $X_n\ge 1$ on $C_n=\left\{X_n>0\right\}$, $\frac{1}{X_n}\le 1$ on $C_n$. Therefore,

$$\frac{1}{X_n}{\mathbf{1}}_{C_n}\le {\mathbf{1}}_{C_n}.$$

Because ${\mathbf{1}}_{C_n}$ is $L^1\left(\mathbb{P}\right)$-integrable, using the Lebesgue dominated convergence theorem (LDCT),

$$\begin{array}{ccc}\hfill \underset{n\to \infty }{lim}\mathbb{E}\left[\frac{1}{X_n}{\mathbf{1}}_{C_n}\right]& =\mathbb{E}\left[{\lim }_{n\to \infty }\frac{1}{X_n}{\mathbf{1}}_{C_n}\right]\hfill & \hfill =\mathbb{E}\left[{\lim }_{n\to \infty }\frac{1}{X_n}·{\lim }_{n\to \infty }{\mathbf{1}}_{C_n}\right]\\ \hfill & =\mathbb{E}\left[\frac{1}{X_{\infty }}{\mathbf{1}}_{{\Omega }_0}\right]=0\hfill \end{array}$$

(6)

Thus, we proved the following lemma,

Lemma 1.

${\mathbb{E}}_0\left[\frac{1}{X_{\infty }}\right]=\mathbb{E}\left[\frac{1}{X_{\infty }}{\mathbf{1}}_{{\Omega }_0}\right]=0$, where ${\mathbb{E}}_0$ and $\mathbb{E}$ are the integrals on respective spaces $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$ and $\left(\Omega ,\mathcal{F},\mathbb{P}\right)$.

Next, recall [18] that $N{M}_{n+1}$ is the number of new mutants at generation $n+1$ equal to

$$N{M}_{n+1}=2X_n-X_{n+1}.$$

Thus,

$$\frac{N{M}_{n+1}}{{\left(2\nu \right)}^{n+1}}=\frac{1}{\nu }{W}_n-W_{n+1}\to \frac{\mu }{\nu }{W}_{\infty }\mathrm{a}.\mathrm{s}.\mathrm{pointwise}.$$

Similarly, because $N{M}_{n+1}\ge 1$ a.s. on $A_{n+1}=\left\{N{M}_{n+1}>0\right\},$ it follows that $\frac{1}{N{M}_{n+1}}\le 1$ on $A_{n+1}$ that implies

$$\frac{1}{N{M}_{n+1}}{\mathbf{1}}_{A_{n+1}}\le {\mathbf{1}}_{A_{n+1}}.$$

By the LDCT and because $A_{n+1}\subseteq C_n,$

$$\begin{array}{cc}\hfill \underset{n\to \infty }{lim}\mathbb{E}\left[\frac{1}{N{M}_{n+1}}{\mathbf{1}}_{A_{n+1}}\right]=\mathbb{E}\left[{\lim }_{n\to \infty }\frac{1}{N{M}_{n+1}}{\mathbf{1}}_{A_{n+1}}\right]& \le \mathbb{E}\left[{\lim }_{n\to \infty }\frac{1}{N{M}_{n+1}}{\mathbf{1}}_{C_n}\right]\hfill \\ \hfill & =\mathbb{E}\left[\frac{1}{N{M}_{\infty }}{\mathbf{1}}_{{\Omega }_0}\right]=0.\hfill \end{array}$$

(7)

Using the same arguments as in the proof of Lemma 1, we easily conclude that

$${\mathbb{E}}_0\left[\frac{1}{N{M}_{\infty }}\right]=0.$$

In a nutshell,

Lemma 2.

${\mathbb{E}}_0\left[\frac{1}{N{M}_{\infty }}\right]=\mathbb{E}\left[\frac{1}{N{M}_{\infty }}{\mathbf{1}}_{{\Omega }_0}\right]=0$, where ${\mathbb{E}}_0$ and $\mathbb{E}$ are the integrals on respective spaces $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$ and $\left(\Omega ,\mathcal{F},\mathbb{P}\right)$.

5. The ${\mathbf{L}}^{\mathbf{2}}\left({\mathbb{P}}_{\mathbf{0}}\right)$-Convergence of ${\widehat{\mathbf{\mu }}}_{\mathbf{n}+\mathbf{1}}$

Theorem 3.

Let $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$ be the confined space introduced in Section 3 and $\left(\Omega ,\mathcal{F},\mathbb{P}\right)$ be the original space. The estimator ${\widehat{\mu }}_{n+1}$ is ${\mathbb{P}}_0$-asymptotically unbiased and it converges to μ in the $L^2\left({\mathbb{P}}_0\right)$-norm.

Proof. 

We use the formula

$${\mathrm{Var}}_0{\widehat{\mu }}_{n+1}={\mathbb{E}}_0{\widehat{\mu }}_{n+1}^2-{\left[{\mathbb{E}}_0{\widehat{\mu }}_{n+1}\right]}^2=\frac{1}{\mathbb{P}\left({\Omega }_0\right)}\mathbb{E}\left[{\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{{\Omega }_0}\right]-{\left\{\frac{1}{\mathbb{P}\left({\Omega }_0\right)}\mathbb{E}\left[{\widehat{\mu }}_{n+1}{\mathbf{1}}_{{\Omega }_0}\right]\right\}}^2.$$

(8)

(i)

${\widehat{\mu }}_{n+1}{\mathbf{1}}_{{\Omega }_0}\le {\widehat{\mu }}_{n+1}$ (dominating sequence that is $\mathbb{P}$-integrable with $\int {\widehat{\mu }}_{n+1}d\mathbb{P}=\mathbb{E}{\widehat{\mu }}_{n+1}=\mu$). Thus, by the LDCT,

$$\underset{n\to \infty }{lim}\mathbb{E}\left[{\widehat{\mu }}_{n+1}{\mathbf{1}}_{{\Omega }_0}\right]=\mathbb{E}\left[\underset{n\to \infty }{lim}{\widehat{\mu }}_{n+1}{\mathbf{1}}_{{\Omega }_0}\right]=\mu \mathbb{P}\left({\Omega }_0\right).$$

implying that

$${\mathbb{E}}_0{\widehat{\mu }}_{n+1}\to \mu$$

(9)

which means that the estimator ${\widehat{\mu }}_{n+1}$ is ${\mathbb{P}}_0$-asymptotically unbiased.

(ii)

Next, we show that ${\widehat{\mu }}_{n+1}$ is $L^2\left({\mathbb{P}}_0\right)$-integrable. First, from (1),

$$\mathbb{E}\left[\left.{\widehat{\mu }}_{n+1}^2\right|X_n\right]=\mathbb{E}\left[\left.{\left(\frac{N{M}_{n+1}}{2X_n}\right)}^2\right|X_n\right]=\frac{1}{4X_n^2}\mathbb{E}\left[\left.N{M}_{n+1}^2\right|X_n\right].$$

Given the $\sigma$-algebra $\sigma \left(X_n\right)$, the r.v. $N{M}_{n+1}\in \left[\mathbb{B}\left(2X_n,\mu \right)\right]$, i.e., $N{M}_{n+1}$ is conditionally binomial with parameters $2X_n$ and $\mu$. Thus,

$$\mathbb{E}\left[\left.{\widehat{\mu }}_{n+1}^2\right|X_n\right]=\frac{1}{4X_n^2}\left[2X_n\mu \nu +4X_n^2{\mu }^2\right]=\frac{1}{X_n}\frac{1}{2}\mu \nu +{\mu }^2$$

implying that

$$\mathbb{E}\left[{\widehat{\mu }}_{n+1}^2\left.{\mathbf{1}}_{C_n}\right|X_n\right]=\frac{1}{X_n}\frac{1}{2}\mu \nu {\mathbf{1}}_{C_n}+{\mu }^2{\mathbf{1}}_{C_n}.$$

(10)

and that

$$\mathbb{E}\left[{\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{C_n}\right]=\frac{1}{2}\mu \nu \mathbb{E}\left[\frac{1}{X_n}{\mathbf{1}}_{C_n}\right]+{\mu }^2\mathbb{P}\left(C_n\right).$$

From (6) and (10), we have

$$\mathbb{E}\left[{\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{C_n}\right]\to {\mu }^2\mathbb{P}\left({\Omega }_0\right).$$

Since ${\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{{\Omega }_0}\le {\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{C_n}$ and since ${\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{C_n}$ is $\mathbb{P}$-integrable, so is ${\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{{\Omega }_0}$. Furthermore, because ${\widehat{\mu }}_{n+1}^2\le 1$ and because ${\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{{\Omega }_0}\to {\mu }^2{\mathbf{1}}_{{\Omega }_0}$ a.s., by the LDCT,

$$\mathbb{E}\left[{\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{{\Omega }_0}\right]\to {\mu }^2\mathbb{P}\left({\Omega }_0\right)$$

implying that

$${\mathbb{E}}_0{\widehat{\mu }}_{n+1}^2=\frac{1}{\mathbb{P}\left({\Omega }_0\right)}\mathbb{E}\left[{\widehat{\mu }}_{n+1}^2{\mathbf{1}}_{{\Omega }_0}\right]\to {\mu }^2.$$

(11)

(iii)

Using (9) and (11),

$$\begin{array}{cc}\hfill & \underset{n\to \infty }{lim}{∥{\widehat{\mu }}_{n+1}-{\mathbb{E}}_0{\widehat{\mu }}_{n+1}∥}_{L^2\left({\mathbb{P}}_0\right)}^2=\underset{n\to \infty }{lim}{\mathrm{Var}}_0{\widehat{\mu }}_{n+1}\hfill \\ \hfill & =\underset{n\to \infty }{lim}\left\{{\mathbb{E}}_0\left[{\widehat{\mu }}_{n+1}^2\right]-{\left({\mathbb{E}}_0{\widehat{\mu }}_{n+1}\right)}^2\right\}={\mu }^2-{\mu }^2=0\hfill \end{array}$$

(12)

or

$$=\frac{1}{\mathbb{P}\left({\Omega }_0\right)}{\mu }^2\mathbb{P}\left({\Omega }_0\right)-\frac{1}{{\mathbb{P}}^2\left({\Omega }_0\right)}{\mu }^2{\mathbb{P}}^2\left({\Omega }_0\right)=0.$$

(iv)

Notice that while ${\widehat{\mu }}_{n+1}$ is an unbiased estimator of $\mu$ with respect to measure $\mathbb{P},$ it need not be unbiased with respect to measure ${\mathbb{P}}_0$. As mentioned, however, ${\widehat{\mu }}_{n+1}$ is asymptotically unbiased with respect measure ${\mathbb{P}}_0$. To establish the $L^2\left({\mathbb{P}}_0\right)$-convergence of ${\widehat{\mu }}_{n+1}$ to $\mu$ we use (9) and (12):

$$\begin{array}{cc}\hfill & {∥{\widehat{\mu }}_{n+1}-\mu ∥}_{L^2\left({\mathbb{P}}_0\right)}\le {∥{\widehat{\mu }}_{n+1}-{\mathbb{E}}_0{\widehat{\mu }}_{n+1}∥}_{L^2\left({\mathbb{P}}_0\right)}+{∥{\mathbb{E}}_0{\widehat{\mu }}_{n+1}-\mu ∥}_{L^2\left({\mathbb{P}}_0\right)}\hfill \\ \hfill & =\sqrt{{\mathrm{Var}}_0{\widehat{\mu }}_{n+1}}+\left|{\mathbb{E}}_0{\widehat{\mu }}_{n+1}-\mu \right|\to 0\mathrm{as}n\to \infty .\hfill \end{array}$$

Therefore, we proved that ${\widehat{\mu }}_{n+1}\stackrel{L_2\left({\mathbb{P}}_0\right)}{\to }\mu$. □

6. The Speed of Convergence

In [18], we validated the pointwise convergence of estimator ${\widehat{\mu }}_{n+1}$ to $\mu$ by simulation and most importantly noticed that the speed of convergence was fast. Considering the replication time of various bacteria (20–40 min), it generally takes not more than 12 to 15 h of lab experiments (or 25 to 30 generations) to attain to a high accuracy approximation of $\mu$. We wondered if similar qualities of ${\widehat{\mu }}_{n+1}$ apply to the $L^2$-convergence (mean square) in some form. The considerations below are carried out in the confined space $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right)$.

Let

$${\xi }_i^{\left[n\right]}=\left\{\begin{array}{cc}1,& i\mathrm{th}\mathrm{non}\mathrm{mutant}\mathrm{in}\mathrm{generation}n\mathrm{does}\mathrm{not}\mathrm{mutate}\\ 0,& \mathrm{otherwise}.\end{array}\right.$$

(13)

Then, ${\eta }_i^{\left[n\right]}=1-{\xi }_i^{\left[n\right]}$ is 1 if the ith nonmutant mutates and 0 otherwise.

We drop the superscript n in ${\xi }^{\prime }s$, and ${\eta }^{\prime }s$ assuming that they are iid Bernoulli r.v.’s with parameters $\nu$ and $\mu$, respectively. Then,

$$X_{n+1}=\sum _{i=1}^{2X_n}{\xi }_i,n=0,1,\dots ,X_0=1$$

and

$$N{M}_{n+1}=2X_n-X_{n+1}=2X_n-\sum _{i=1}^{2X_n}{\xi }_i=\sum _{i=1}^{2X_n}{\eta }_i$$

implying that

$${\widehat{\mu }}_{n+1}=\frac{{\displaystyle \sum _{i=1}^{2X_n}}{\eta }_{i·}}{2X_n},\mathrm{where}{\eta }_i\in \mathbb{B}\left[1,\mu \right].$$

Setting $X_n=\mathfrak{X}$ (constant), ${\widehat{\mu }}_{n+1}$ is the sample mean of $2\mathfrak{X}$ Bernoullir.v.’s and thus with $\mathfrak{X}$ large, by the CLT (central limit theorem), ${\widehat{\mu }}_{n+1}$ has a Gaussian distribution with parameters $\mu$ (mean) and $\mu \nu /2\mathfrak{X}$ (variance). Now, with n being just 10 or larger and $\mu$ small, $X_n=\mathfrak{X}$ must be a large number. Thus, with simulation, because $X_n$ is obtained empirically it is a constant and as noticed ${\widehat{\mu }}_{n+1}\in N\left[\mu ,\mu \nu /2\mathfrak{X}\right]$, good for all $n\ge 10$. The variance $\mu \nu /\mathfrak{X}$ of ${\widehat{\mu }}_{n+1}$ is small and must approach 0.

Suppose we “observe” bacterial replication on M plates where each colony starts with a single parent. Denote ${\widehat{\mu }}_{n+1}^{\left[k\right]}$ the estimator associated with kth plate, $k=1,\dots ,M$. For notational brevity, we will use the same symbol ${\widehat{\mu }}_{n+1}^{\left[k\right]}$ for its estimate (i.e., observed value). Denote by ${\overline{\mu }}_{n+1}$ the sample mean of M estimators (estimates) ${\widehat{\mu }}_{n+1}^{\left[k\right]},k=1,\dots ,M$. Thus,

$${\widehat{\sigma }}_{n+1}^2=\frac{1}{M-1}\sum _{k=1}^M{\left({\widehat{\mu }}_{n+1}^{\left[k\right]}-{\overline{\mu }}_{n+1}\right)}^2$$

is the unbiased sample variance that is known to be consistent and weakly convergent to $\mu \nu /2\mathfrak{X}$ as $M\to \infty$. However, we expectthat ${\widehat{\sigma }}_{n+1}$ rapidly converges to 0 with just reasonable n and M. This would be a heuristic feel of the true convergence speed of the variance of ${\widehat{\mu }}_{n+1}$ (cf. Equation (12)).

Example 1.

To validate the results above, we will simulate the process for several fixed values of μ and compute the sample estimates of the mean and variance of μ, i.e., ${\widehat{\mu }}_n$ and ${\widehat{\sigma }}_n^2$, for generations $n=1,2,\dots ,25$ for different numbers of plates ranging from $M=1$ to $M=20$. The goals are twofold: (1) to validate the accuracy of the estimators empirically and (2) to experimentally examine the convergence of the estimators for various choices of generations n and plates M. The estimators are more promising for practical use if convergence occurs for relatively small n and M.

To compute estimates, we must simulate M paths of the process for each experiment. For each sample path, we start with $X_0=1$ nonmutant at the 0th generation and, for each generation $n\ge 1$, we assume the number of nonmutants is a binomial random variable $X_n\in \mathbb{B}\left(2X_{n-1},1-\mu \right)$ for a preselected mutation rate $\mu$. We generated M paths from $n=1$ to $n=30$. The estimated mutation rate is computed at each generation $n=0$ to $n=29$ as ${\widehat{\mu }}_{n+1}=N{M}_{n+1}/2X_n$ (rewritten as $1-X_{n+1}/2X_n$ for convenience) and the mean ${\overline{\mu }}_n$ is taken for each experiment with M plates. The sample variance will also computed. Ideally, the estimate ${\widehat{\mu }}_n$ will be near $\mu$ and the estimate ${\widehat{\sigma }}_n^2$ will approach 0 quickly for relatively small M.

We follow this approach in several examples. In Figure 1, we let $\mu =0.01$ and simulated the experiment for 25 generations for each selection of plates $M\in \{2,5,10,20\}$ (We did not use $M=1$ as we evaluate sample variance with degree of freedom of 1). There are several intuitive occurrences here.

Since ${\widehat{\mu }}_{n+1}=N{M}_{n+1}/2X_n$, the estimate is 0 if there are no new mutants in generation $n+1$. As we see in Figure 1, the estimates are indeed 0 before the first mutant appears in each experiment. An interesting pattern emerges here: for $M=20$ and $M=10$ plates, the first mutant appears in generation 1 and 2, respectively; for $M=5$ plates, the first mutant appears at generation 3; and for $M=2$, the first mutant appears at generation 4. The time of the first mutation is inversely related to the number of plates, which makes sense because there are more opportunities for mutations with more plates and, therefore, the estimate is nonzero sooner with more plates.

Second, the estimates behave erratically between the initial escape from 0 until roughly generation 7–8. During the early generations, there are smaller nonmutant populations in the plates, meaning intuitively that the random mutations have larger impacts on the estimates.

As such, variance is higher when the prior generation’s nonmutant population is smaller, and this nonmutant population grows over time in expectation since the mutation rate $\mu \ll 0.05$, so we should expect to high variances in early generations. After 7 or 8 generations, the variance seems to dissipate enough to allow the estimates start to approach the true mutation rate $\mu =0.01$.

Continuing beyond generation 8, estimates in all experiments approach $\mu =0.01$, with the experiments with fewer plates, especially $M=2$, taking a longer time to converge—again, this is expected as the conditional variance of ${\overline{\mu }}_{n+1}$ is inversely proportional to $M^2$, and a shrinking variance should result in convergence of the unbiased estimator.

Clearly, and unsurprisingly, the estimates have higher quality for larger generation n and larger number of plates M, but more interestingly, what does this mean for experimental use of the estimator? First, we define a modified relative error for the estimate

$$e_{0.01}=max\left(\left|\frac{{\overline{\mu }}_n-\mu }{\mu }\right|,0.01\right),$$

which is the relative error or 0.01, whichever is larger. Figure 2 is a heatmap of this modified relative error for each pair n and M from our simulations.

Thresholding the relative error as $e_{0.01}$ allows us to see clearly in Figure 2 that relative errors of less than 0.01 occurred for all generations $n\ge 20$ regardless of the number of plates, but for more plates, near $M=20$, these small errors occur reliably earlier, at generations $n<15$. Thus, even low numbers of plates can be used if the experiment can continue for more generations. This is a promising finding for biology lab protocol, as growing more than 20–25 dishes per experiment is highly infeasible. After 15–20 generations (starting with 1 starting cells), we have the mean converged.

Beyond the behavior of the mean mutation rate themselves, we can consider the estimated variance of the mutation rate, ${\widehat{\sigma }}_{n+1}^2$. Figure 3 shows the estimated sample variance for several experiments with $M=2,5,10,$ and 20 plates.

As we see in Figure 3, the variance converges to 0 within 10–12 generations regardless of the number of plates, as prior analysis suggested should happen.

7. Estimator ${\widehat{\mathit{p}}}_{\mathit{n}}$

We now assume that once a cell or microbe mutates, it turns into one of the two mutant types, with probabilities p and $q=1-p,$ respectively. The type i mutant then divides in accordance with a deterministic branching process throughout the generations and it does not mutate any further, nor does it alter its type. In particular, we can also interpret mutation types as a continual existence upon mutation (type 1) with probability p or the death (type 2) with probability $q=1-p$. Note that our analysis in this section is interchangeably carried out in the original and confined probability spaces $\left(\Omega ,\mathcal{F},\mathbb{P}\right)$ and $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0\right),$ respectively. In all cases, the underlying functionals can be distinguished by the absence or presence of subscript 0.

Denote by

$${\widehat{p}}_n=\frac{N{M}_n^{\left[1\right]}}{N{M}_n}\mathrm{and}{\widehat{q}}_n=\frac{N{M}_n^{\left[2\right]}}{N{M}_n}$$

(14)

the proportions of type 1 and 2 new mutants at the nth generation from among the $N{M}_n$ new mutants, respectively. We will show that the estimator ${\widehat{p}}_n$ is unbiased and consistent (and so is ${\widehat{q}}_n$).

Recall that $A_n=\left\{N{M}_n>0\right\}$. Denote $g_n\left(z\right)=\mathbb{E}{z}^{X_n}$. Below, we observe an interesting and unexpected result concerning the sequence $\left\{A_n\right\}$. Whereas it obviously is not monotone, as the following two assertions show, the sequence $\left\{\mathbb{P}\left(A_n\right)\right\}$ of their measures under a weak and natural assumption, is strictly monotone increasing (we recollect that the sequence $\left\{\mathbb{P}\left(C_n\right)\right\}$ with respect to nonmutants is monotone nonincreasing) and furthermore, like $\left\{\mathbb{P}\left(C_n\right)\right\}$, $\left\{\mathbb{P}\left(A_n\right)\right\}$ converges to $\mathbb{P}\left({\Omega }_0\right)=1-\varphi$.

Theorem 4.

The following two assertions are valid.

(i)

$\mathbb{P}\left(A_n^c\right)=g_{n-1}\left({\nu }^2\right).$

(ii)

If $\mu <{\nu }^2<1,$ the sequence $\left\{\mathbb{P}\left(A_n\right)\right\}$ is strictly monotone increasing.

Proof. 

(i)

Let $h_n\left(z\right)=\mathbb{E}{z}^{N{M}_n}$. Then, because $N{M}_{n+1}\in \left[\mathbb{B}\left(2X_n,\mu \right)\right]$ given $\sigma \left(X_n\right)$,

$$\begin{array}{cc}\hfill & h_{n+1}\left(z\right)=\mathbb{E}\left[\mathbb{E}\left[z^{N{M}_{n+1}}{X}_n\right]\right]=\mathbb{E}\left[{\left(\mu z+\nu \right)}^{2X_n}\right]=g_n\left[{\left(\mu z+\nu \right)}^2\right].\hfill \\ \hfill & \Rightarrow \mathbb{P}\left(A_{n+1}^c\right)=h_{n+1}\left(0\right)=g_n\left({\nu }^2\right).\hfill \end{array}$$

(ii)

We show that $\left\{g_n\left({\nu }^2\right)\right\}$ is strictly monotone decreasing.

First off, $p\left(z\right)={\left(\nu z+\mu \right)}^2$ is strictly monotone increasing for $z\in \left[0,1\right]$. Now since ${\nu }^2<1$ and $p\left(1\right)=max\left\{p\left(z\right):z\in \left[0,1\right]\right\}=1$,

$$p\left({\nu }^2\right)=g_1\left({\nu }^2\right)={\left({\nu }^3+\mu \right)}^2<1.$$

Now, we show that

$$g_{n+1}\left({\nu }^2\right)<g_n\left({\nu }^2\right)$$

or that

$$g_{n+1}\left({\nu }^2\right)=p\left[g_n\left({\nu }^2\right)\right]<g_n\left({\nu }^2\right).$$

Let $g_n\left({\nu }^2\right)=z$. Consider $p\left(z\right)\diamond z$, where ⋄ stands for one of the three relations, $<,>,\mathrm{or}=$. We have

$${\left(\nu z+\mu \right)}^2-z\diamond 0$$

or

$${\nu }^2{z}^2-\left(1-2\nu \mu \right)z+{\mu }^2\diamond 0$$

The discriminant of the equation ${\nu }^2{z}^2-\left(1-2\nu \mu \right)z+{\mu }^2=0$ is

$$\begin{array}{cc}\hfill D& ={\left(1-2\nu \mu \right)}^2-4{\nu }^2{\mu }^2=1-4\mu \nu =1+{\nu }^2-4\nu ={\left(2\nu -1\right)}^2\hfill \\ \hfill & =\left\{\begin{array}{cc}=0,& \mu =\nu =\frac{1}{2}\\ >0,& \nu \ne \frac{1}{2}\end{array}\right..\hfill \end{array}$$

Furthermore, after some algebra,

$${\nu }^2{z}^2-\left(1-2\nu \mu \right)z+{\mu }^2={\nu }^2\left(z-z_1\right)\left(z-z_2\right),$$

where $z_1=\frac{{\mu }^2}{{\nu }^2}$ (is the extinction probability $\varphi$) and $z_2=1$. Therefore, ${\nu }^2{z}^2-\left(1-2\nu \mu \right)z+{\mu }^2<0$ (and ⋄ is <) for all $z\in \left(\frac{{\mu }^2}{{\nu }^2},1\right)$. In conclusion,

$$\mathrm{if}z={\nu }^2,\mathrm{then}{\nu }^2>\frac{{\mu }^2}{{\nu }^2}\mathrm{whenever}{\nu }^2>\mu .$$

Therefore, if ${\nu }^2>\mu$, $p\left({\nu }^2\right)=g_1\left({\nu }^2\right)={\left({\nu }^3+\mu \right)}^2<g_0\left({\nu }^2\right)={\nu }^2$. Next, since $p\left(z\right)$ is strictly monotone increasing in $z\in \left[0,1\right]$ and because $p\left({\nu }^2\right)<{\nu }^2$, it follows that

$$g_2\left({\nu }^2\right)=p\left[p\left({\nu }^2\right)\right]<g_1\left({\nu }^2\right)=p\left({\nu }^2\right).$$

Finally, since $p\left(z\right)$ is strictly monotone increasing, by the principle of mathematical induction,

$$g_{n+1}\left({\nu }^2\right)=\underset{n+1p’\mathrm{s}}{\underbrace{p\left(p\left(p\dots \left({\nu }^2\right)\right)\right)}}<g_n\left({\nu }^2\right)=\underset{np’\mathrm{s}}{\underbrace{p\left(p\dots \left({\nu }^2\right)\right)}}$$

that proves the assertion that $\left\{g_n\left({\nu }^2\right)\right\}$ is a strictly monotone decreasing sequence and so is the sequence $\left\{\mathbb{P}\{N{M}_n=0\}=\mathbb{P}\left(A_n^c\right)\right\}$. □

Corollary 1.

${lim}_{n\to \infty }\mathbb{P}\left(A_n\right)={lim}_{n\to \infty }\mathbb{P}\left\{N{M}_n>0\right\}=1-\varphi =\mathbb{P}\left({\Omega }_0\right)=1-\frac{{\mu }^2}{{\nu }^2}$.

Proof. 

Let ${\alpha }_n=g_n\left({\nu }^2\right)$. Then, from the equation

$$g_{n+1}\left({\nu }^2\right)=p\left[g_n\left({\nu }^2\right)\right]\mathrm{holding}\mathrm{for}n=0,1,\dots ,$$

we have

$${\alpha }_{n+1}=p\left({\alpha }_n\right),n=0,1,\dots$$

(15)

Furthermore, since $\left\{{\alpha }_n\right\}$ is monotone decreasing, $\alpha :={lim}_{n\to \infty }{\alpha }_n$ exists and it is unique. Because $p\left(z\right)$ is continuous, this limit from (15) can be found from the equation

$$\alpha =p\left(\alpha \right).$$

From the theory of branching processes, we know (cf. Kimmel and Axelrod [43]) that the extinction probability $\varphi$ is the smallest positive root of the equation $z=p\left(z\right)$. In our case, $p\left(z\right)={\left(\nu z+\mu \right)}^2$ and the extinction probability $\varphi <1$ if the mean $2\nu >1$ which is routinely met when $\mu <\frac{1}{2}$. The equation $z={\left(\nu z+\mu \right)}^2$ has two roots: 1 and the other one $\varphi =\frac{{\mu }^2}{{\nu }^2}$ which is the smallest of two. Consequently,

$$\alpha =\underset{n\to \infty }{lim}\mathbb{P}\left\{N{M}_n=0\right\}=\varphi =\frac{{\mu }^2}{{\nu }^2}$$

ends up being equal to the extinction probability. Therefore,

$${\displaystyle \underset{n\to \infty }{lim}\mathbb{P}\left\{N{M}_n>0\right\}=\mathbb{P}\left({\Omega }_0\right)=1-\frac{{\mu }^2}{{\nu }^2}.}$$

□

Since ${\widehat{p}}_n=\frac{N{M}_n^{\left[1\right]}}{N{M}_n}$, given $N{M}_n$, $N{M}_n^{\left[1\right]}\in \left[\mathbb{B}\left(N{M}_n,p\right)\right]$ and therefore,

$$\mathbb{E}\left[\left.{\widehat{p}}_n\right|N{M}_n\right]=\frac{1}{N{M}_n}\mathbb{E}\left[\left.N{M}_n^{\left[1\right]}\right|N{M}_n\right]=\frac{1}{N{M}_n}N{M}_{n}p=p\Rightarrow \mathbb{E}{\widehat{p}}_n=p.$$

Furthermore,

$$\begin{array}{cc}\hfill & \mathbb{E}\left[{\widehat{p}}_n{\mathbf{1}}_{A_n}N{M}_n\right]=p{\mathbf{1}}_{A_n}\Rightarrow \mathbb{E}{\widehat{p}}_n{\mathbf{1}}_{A_n}=p\mathbb{P}\left(A_n\right)\hfill \\ \hfill \Rightarrow & {\mathbb{E}}_n{\widehat{p}}_n=p.\hfill \end{array}$$

Consequently, ${\widehat{p}}_n$ is $\mathbb{P}$- and ${\mathbb{P}}_n$-unbiased, where ${\mathbb{P}}_n\left(·\right)=\mathbb{P}\left(·\cap A_n\right)/\mathbb{P}\left(A_n\right)$.

Theorem 5.

If $\nu >\frac{1}{2}$, ${\widehat{p}}_n$ converges to p in the $L^2\left({\mathbb{P}}_n\right)$-norm, $L^2\left({\mathbb{P}}_0\right)$-norm, and a.s. pointwise. Furthermore, ${\widehat{p}}_n$ is asymptotically ${\mathbb{P}}_0$-unbiased.

Proof. 

(i)

We start off with

$$\mathbb{E}\left[\left.{\widehat{p}}_n^2{\mathbf{1}}_{A_n}\right|N{M}_n\right]=\frac{1}{N{M}_n^2}{\mathbf{1}}_{A_n}\left[N{M}_{n}pq+N{M}_n^2{p}^2\right]=\left[p^2+\frac{1}{N{M}_n}pq\right]{\mathbf{1}}_{A_n}$$

(16)

implying

$$\mathbb{E}\left[{\widehat{p}}_n^2{\mathbf{1}}_{A_n}\right]=pq\mathbb{E}\left[\frac{1}{N{M}_n}{\mathbf{1}}_{A_n}\right]+p^2\mathbb{P}\left(A_n\right).$$

(17)

Because $\frac{1}{N{M}_n}{\mathbf{1}}_{A_n}\le 1$, ${\widehat{p}}_n^2{\mathbf{1}}_{A_n}$ is $\mathbb{P}$-integrable. By Lemma 2 (cf. Equation (7)),

$${\mathrm{Var}}_n{\widehat{p}}_n=pq\frac{1}{\mathbb{P}\left(A_n\right)}\mathbb{E}\left[\frac{1}{N{M}_n}{\mathbf{1}}_{A_n}\right]+p^2-p^2=pq\frac{1}{\mathbb{P}\left(A_n\right)}\mathbb{E}\left[\frac{1}{N{M}_n}{\mathbf{1}}_{A_n}\right]\to 0$$

Thus, because ${\widehat{p}}_n$ is ${\mathbb{P}}_n$-unbiased,

$${∥{\widehat{p}}_n-p∥}_{L^2\left({\mathbb{P}}_n\right)}={∥{\widehat{p}}_n-{\mathbb{E}}_n{\widehat{p}}_n∥}_{L^2\left({\mathbb{P}}_n\right)}=\sqrt{{\mathrm{Var}}_n{\widehat{p}}_n}\to 0.$$

(ii)

Now we show that

$${\widehat{\mu }}_{n+1}^{\left[1\right]}=N{M}_{n+1}^{\left[1\right]}/2X_n\to \mu p\mathrm{a}.\mathrm{s}.$$

The latter is true due to the following: As per our assumption, every new mutation takes place with a constant probability $\mu$. The new mutants form process $\left\{N{M}_n\right\}$. The estimator ${\widehat{\mu }}_n$ of $\mu$ is consistent in the sense that ${\widehat{\mu }}_n\to \mu$ a.s. pointwise. If a new mutant is formed from a nonmutant, in our present model, it becomes a type 1 new mutant with a constant probability p, i.e., a type 1 new mutant emerges from a nonmutant with probability $p\mu$. We then focus on $X_n$ and $N{M}_{n+1}^{\left[1\right]}$ of nonmutants in generation n and type 1 new mutants in generation $n+1$, respectively. We can then reduce this model to the original one with $p\mu$ being a new mutation rate. Thus, we define the estimator ${\widehat{\mu }}_{n+1}^{\left[1\right]}$ of $p\mu$ as $N{M}_{n+1}^{\left[1\right]}/2X_n$, and the latter must converge to $p\mu$ a.s. pointwise [18].

Consequently, from Theorem 3, ${\widehat{\mu }}_{n+1}^{\left[1\right]}$ converges to $p\mu$ in the $L^2\left({\mathbb{P}}_0\right)$-norm.

(iii)

Now, we turn to ${\widehat{p}}_{n+1}$ defined as $N{M}_{n+1}^{\left[1\right]}/N{M}_{n+1}$. Dividing the numerator and denominator by $2X_n$, we get

$${\widehat{p}}_{n+1}=\frac{N{M}_{n+1}^{\left[1\right]}/2X_n}{N{M}_{n+1}/2X_n}=\frac{{\widehat{\mu }}_{n+1}^{\left[1\right]}}{{\widehat{\mu }}_{n+1}}\to \frac{p\mu }{\mu }=p,\mathrm{a}.\mathrm{s}.\mathrm{on}{\Omega }_0.$$

More specifically, if the numerator converges to $p\mu$ pointwise on ${\Omega }_1$ and the denominator converges pointwise to $\mu$ on ${\Omega }_2$ such that $\mathbb{P}\left({\Omega }_1\right)=\mathbb{P}\left({\Omega }_2\right)=\mathbb{P}\left({\Omega }_0\right)$,then ${\widehat{p}}_{n+1}\to p$ pointwise on set ${\Omega }_1\cap {\Omega }_2$ and clearly $\mathbb{P}\left({\Omega }_1\cap {\Omega }_2\right)=\mathbb{P}\left({\Omega }_0\right)$. Thus, the convergence of ${\widehat{p}}_n$ to p is a.s. pointwise on ${\Omega }_0$.

Consequently,

$${\widehat{p}}_n^2\to p^2\mathrm{a}.\mathrm{s}.\mathrm{pointwise}\mathrm{on}{\Omega }_0$$

(18)

and

$${\widehat{p}}_n^2{\mathbf{1}}_{{\Omega }_0}\to p^2{\mathbf{1}}_{{\Omega }_0}\mathrm{a}.\mathrm{s}.$$

Because ${\widehat{p}}_n^2\le 1$ for all $\omega$’s, by the LDCT and Equation (18),

$$\underset{n\to \infty }{lim}\mathbb{E}{\widehat{p}}_n^2{\mathbf{1}}_{{\Omega }_0}=\mathbb{E}\left[{\lim }_{n\to \infty }{\widehat{p}}_n^2{\mathbf{1}}_{{\Omega }_0}\right]=p^2\mathbb{P}\left({\Omega }_0\right).$$

(iv)

${\mathrm{Var}}_0{\widehat{p}}_n={\mathbb{E}}_0{\widehat{p}}_n^2-{\left[{\mathbb{E}}_0{\widehat{p}}_n\right]}^2=\frac{1}{\mathbb{P}\left({\Omega }_0\right)}\mathbb{E}\left[{\widehat{p}}_n^2{\mathbf{1}}_{{\Omega }_0}\right]-\frac{1}{{\mathbb{P}}^2\left({\Omega }_0\right)}{\mathbb{E}}^2\left[{\widehat{p}}_n{\mathbf{1}}_{{\Omega }_0}\right]$

${\widehat{p}}_n{\mathbf{1}}_{{\Omega }_0}\le {\widehat{p}}_n\mathrm{and}{\widehat{p}}_n\mathrm{is}\mathbb{P}\text{-}\mathrm{integrable}\left(\int {\widehat{p}}_{n}d\mathbb{P}=p\right)$

Thus by the LDCT,

$$\underset{n\to \infty }{lim}\mathbb{E}{\widehat{p}}_n{\mathbf{1}}_{{\Omega }_0}=p\mathbb{P}\left({\Omega }_0\right)$$

and

$${\mathbb{E}}_0{\widehat{p}}_n\to p,$$

so that ${\widehat{p}}_n$ is asymptotically ${\mathbb{P}}_0$-unbiased. Consequently,

$${\mathrm{Var}}_0{\widehat{p}}_n\to p^2-p^2=0.$$

(v)

${∥{\widehat{p}}_n-p∥}_{L^2\left({\mathbb{P}}_0\right)}\le {∥{\widehat{p}}_n-{\mathbb{E}}_0{\widehat{p}}_n∥}_{L^2\left({\mathbb{P}}_0\right)}+{∥{\mathbb{E}}_0{\widehat{p}}_n-p∥}_{L^2\left({\mathbb{P}}_0\right)}$

$=\sqrt{\mathrm{Var}{\widehat{p}}_n}+\left|{\mathbb{E}}_0{\widehat{p}}_n-p\right|\to 0.$

□

Remark 1.

(i)

If the second mutation type of bacteria are dead, the formula ${\widehat{p}}_{n+1}=\frac{N{M}_{n+1}^{\left[1\right]}}{N{M}_{n+1}}$ is impractical for a lab utility, since the dead cells are not counted. Here is a way out.

First, notice that in this case, $N{M}_{n+1}^{\left[1\right]}=M_{n+1}-2M_n$ (i.e., the same formula as for $N{M}_{n+1}$ with only one mutation type). This is exactly because the dead bacteria cannot be counted. Next, we have

$$2X_n=X_{n+1}+N{M}_{n+1}^{\left[1\right]}+N{M}_{n+1}^{\left[2\right]},$$

(19)

so that the “deficit” of the count is

$$N{M}_{n+1}^{\left[2\right]}=2X_n-X_{n+1}-N{M}_{n+1}^{\left[1\right]}=2X_n-X_{n+1}-M_{n+1}+2M_n.$$

(20)

Therefore,

$${\widehat{p}}_{n+1}=\frac{N{M}_{n+1}^{\left[1\right]}}{N{M}_{n+1}}=\frac{N{M}_{n+1}^{\left[1\right]}}{N{M}_{n+1}^{\left[1\right]}+N{M}_{n+1}^{\left[2\right]}}=\frac{M_{n+1}-2M_n}{2X_n-X_{n+1}}.$$

(21)

(ii)

If the second mutation type are not dead bacteria, then $N{M}_{n+1}^{\left[1\right]}\ne M_{n+1}-2M_n$, but

$$N{M}_{n+1}^{\left[1\right]}=M_{n+1}^{\left[1\right]}-2M_n^{\left[1\right]}.$$

(22)

However, Equation (19) still applies, whereas Equation (20) can be modified as

$$N{M}_{n+1}^{\left[2\right]}=2X_n-X_{n+1}-N{M}_{n+1}^{\left[1\right]}=2X_n-X_{n+1}-M_{n+1}^{\left[1\right]}+2M_n^{\left[1\right]}.$$

(23)

Obviously, Equation (21) can be modified as follows

$${\widehat{p}}_{n+1}=\frac{M_{n+1}^{\left[1\right]}-2M_n^{\left[1\right]}}{2X_n-X_{n+1}},$$

(24)

where superscripts $\left[1\right]$ and $\left[2\right]$ stand for the first and second type mutants, respectively.

In a nutshell, if type one bacteria die upon mutation, the formula for the estimator ${\widehat{p}}_n$ of p is simple and it does not require one to distinguish the bacteria types. This is not the case with two surviving types, as per Equation (24).

8. Simulation of Estimators for Types 1 and 2

To validate the convergence and its speed, we render simulation analogous to that in Section 6. After generation of nonmutant paths similar to Section 6, we evaluate new mutant (all types) at each generation using $N{M}_{n+1}=2X_n-X_{n+1}$. For each new mutant generation, we run a binomial random number generator with distribution $\mathbb{B}\left(N{M}_{n+1},p\right)$ for mutation type 1, assuming type 2 is the rest of the mutation population with distribution $\mathbb{B}\left(N{M}_{n+1},q=1-p\right)$, because $N{M}_{n+1}=N{M}_{n+1}^{\left[1\right]}+N{M}_{n+1}^{\left[2\right]}$. This new mutant type i generation is added to the current population of type i mutant (previous generation mutant doubled): $M_{n+1}^{\left[i\right]}=N{M}_{n+1}^{\left[i\right]}+2·M_n^{\left[i\right]}$, with the assumption that $M_0^{\left[i\right]}=0$ (initially, we have 0 mutants). The algorithm is looped for n generations and we store the result in a new $(n+1)\times M$ matrix with first row being zero row. Type 1 mutation rate is estimated by Equation (21): ${\widehat{p}}_{n+1}=\frac{M_{n+1}^{\left[1\right]}-2M_n^{\left[1\right]}}{2X_n-X_{n+1}}$ and the output was an $n\times M$ matrix. Mean mutation and variance of the mutation rate were evaluated using the sample variance equation in Section 6, knowing the convergence is established with Theorem 5.

For this simulation, we let the probability that a newly mutated bacterium is of type 1 be $p=0.25$. Figure 4 is the simulation results—specifically, the mean of the type 1 mutation rate estimate.

Similar to simulation in Section 6, we see a pattern here: the estimate is 0 if there are no new mutants. In this simulation, first mutation does not occur until generation 3 for $M=20$, generation 4 for both $M=10$ and $M=5$, and generation 7 for $M=2$. The time of the first type 1 mutant is also inversely related to the number of plates, due to the appearance of type 1 mutation depends on the time of first mutation, which is because of, as discussed previously, more opportunities for mutations to occur with more plates.

The estimates stabilize after 9–10 generations, which is a couple of generations lag after first mutation occurs, as established above. This can be explained due to the probability of getting type 2 mutation is larger than type 1 and, due to this type 1 mutation, take more generation to appear. The following modified relative error heatmap (Figure 5) also describes this phenomenon.

Figure 5 shows the similar pattern as Figure 2: for all generations $n\ge 20$, regardless of the number of plates, the relative error is less than 0.01; the more plates there are, the less occurrence of large relative errors in earlier generations. For all number of plates, the threshold for small relative error is $n\ge 16$; which is 1 to 2 generation behind the observed threshold from Figure 2 ($n\ge 14$). Similar to Section 6, this result is a promising finding for biology researchers as the number of dishes is feasibly small to be made in laboratory settings.

As described in Figure 6, similar to simulation in Section 6, the estimated variance of ${\widehat{p}}_n$ quickly converges to 0 after 13–17 generations (which is also a couple of generations behind the convergence of variance discussed in Section 6), regardless of the number of plates.

9. Summary

In 2012, Niccum et al. [18] offered a simple stochastic estimator ${\widehat{\mu }}_n$ of $\mu$ in a most rigorous fashion proving that it was unbiased and consistent. The latter meant that the sequence $\left\{{\widehat{\mu }}_n\right\}$ was almost surely pointwise convergent to $\mu$ on a non extinct set. Best of all, the pointwise convergence proved to be fast, as it was validated by countless simulation runs for various values of $\mu$.

In this paper, we studied the stochastic process describing bacterial mutation that begins at some point of their division. In our earlier paper [18] (by two of the present coauthors), we proposed a stochastic estimator ${\widehat{\mu }}_n$$\left(n=1,2,\dots \right)$ of the mutation rate $\mu$ (or more precisely, mutation probability) and showed that the sequence $\left\{{\widehat{\mu }}_n\right\}$ converges to $\mu$ a.s. pointwise on the non extinct set ${\Omega }_0\subseteq \Omega$, where $\left(\Omega ,\mathcal{F},\mathbb{P}\right)$ is the probability space on which all processes were considered. $\mathbb{P}\left({\Omega }_0\right)=1-\frac{{\mu }^2}{{\nu }^2}$, where $\nu =1-\mu$, and thus, $\mathbb{P}\left({\Omega }_0\right)\approx 1$ if $\mu$ is small which is the case in many practical situations. The stochastic estimator ${\widehat{\mu }}_n$ obeys Equation (1),

$${\widehat{\mu }}_{n+1}=\frac{N{M}_{n+1}}{2X_n},$$

where $X_n$ is the number of nonmutants in generation n and $N{M}_{n+1}$ is the number of newly formed mutants in generation $n+1$.

One of the subjects of interest was to establish a different type of convergence, namely $L^2$. We formed the confined space $\left({\Omega }_0,{\mathcal{F}}_0=\mathcal{F}\cap {\Omega }_0,{\mathbb{P}}_0=\mathbb{P}/\mathbb{P}\left({\Omega }_0\right)\right)$ and established the following results.

1.

The sequence of estimators $\left\{{\widehat{\mu }}_n\right\}$ is ${\mathbb{P}}_0$-asymptotically unbiased (that is, ${\mathbb{E}}_0{\widehat{\mu }}_{n+1}\to \mu$) and it converges to $\mu$ in the $L^2\left({\mathbb{P}}_0\right)$-norm.

2.

In our earlier paper [18], even though we did not rigorously established the speed of pointwise convergence of $\left\{{\widehat{\mu }}_n\right\}$ to $\mu$, we conducted simulation as well as numerous times thereafter, showing that the convergence was rapid. With the assumed $\mu =7.6\times {10}^{-6}$ (the number produced in our lab results), a good accuracy was observed in the 30th generation. We also wanted to test the $L^2$-speed of convergence in the context of our paper. We used ${\widehat{\mu }}_{n+1}$ in the form

$${\widehat{\mu }}_{n+1}=\frac{{\displaystyle \sum _{i=1}^{2X_n}}{\eta }_{i·}}{2X_n},\mathrm{where}{\eta }_i\in \mathbb{B}\left[1,\mu \right]$$

setting $X_n$ as a constant, due to the step-by-step simulation procedures. With generation n relatively small, like 10 or greater, because $X_n$ was large, ${\widehat{\mu }}_{n+1}\in N\left[\mu ,\mu \left(1-\mu \right)/2X_n\right]$ (Gaussian) using the CLT argument. The variance of ${\widehat{\mu }}_{n+1}$ is thus $\mu \left(1-\mu \right)/2X_n$ being virtually zero. Yet, we formed the unbiased sample variance

$${\widehat{\sigma }}_{n+1}^2=\frac{1}{M-1}\sum _{k=1}^M{\left({\widehat{\mu }}_{n+1}^{\left[k\right]}-{\overline{\mu }}_{n+1}\right)}^2$$

where M is the number of plates, with bacterial colonies on each and with ${\overline{\mu }}_{n+1}$ being the sample mean of M estimators (estimates) ${\widehat{\mu }}_{n+1}^{\left[k\right]},k=1,\dots ,M$. The conclusion thus was that ${\widehat{\sigma }}_{n+1}^2$ will converge to zero with M not very large. Indeed, in Example 1, we fully validated it with $\mu =0.01$, $M=2,5,10,20$ and with $n\ge 10$. Granted, the convergence would be slower for a smaller $\mu$, but the pointwise convergence of ${\widehat{\mu }}_n$ to $\mu \approx 7.6\times {10}^{-6}$ was around $n=30$ and with $M=10$ or 20.

3.

We further generalized our model by allowing bacteria mutate into one of the two mutant types, with constant probabilities p and q. One of the interpretation of this model is when the second type could be a dead bacterium (that could be attacked by viruses or chemical agent). As mentioned in the introduction, it can depose its DNA to be picked up by other bacteria. We proposed two respective estimators

$${\widehat{p}}_n=\frac{N{M}_n^{\left[1\right]}}{N{M}_n}\mathrm{and}{\widehat{q}}_n=\frac{N{M}_n^{\left[2\right]}}{N{M}_n}$$

and proved that they were unbiased and consistent. Here NM${}^{\left[i\right]}$ is the number of newly formed type 1 mutants. Namely, we proved in Theorem 5 that ${\widehat{p}}_n$ converged to p in the $L^2\left({\mathbb{P}}_n\right)$-norm$,L^2\left({\mathbb{P}}_0\right)$-norm, and a.s. pointwise, and ${\widehat{p}}_n$ was asymptotically ${\mathbb{P}}_0$-unbiased. To prove this result, we had to show that the sequence $\left\{\mathbb{P}\left(A_n\right)\right\}$ was strictly monotone increasing and converged to $\mathbb{P}\left({\Omega }_0\right)=1-\frac{{\mu }^2}{{\left(1-\mu \right)}^2}$. As in the above case, to test the speed of the $L^2$ convergence, we ran simulation with $\mu =0.01,p=0.25,$ and the number of plates $M=2,5,10,20$ showing a good speed for the mean and variance for $n\approx 20$ and $15,$ respectively.