Genetic composition of an exponentially growing cell population

TL;DR

This paper analyzes a DNA evolution model in growing cell populations, revealing that most sites mutate finitely many times and independence across sites generally holds, with some rare sites showing correlated mutations, challenging the infinite sites assumption.

Contribution

It introduces a model that relaxes the infinite sites assumption and provides insights into mutation independence and frequency in expanding cell populations.

Findings

Most sites mutate finitely many times in large populations.

Independence across sites generally holds except for rare highly mutated sites.

The model's predictions are applied to estimate mutation rates in lung cancer.

Abstract

We study a simple model of DNA evolution in a growing population of cells. Each cell contains a nucleotide sequence which randomly mutates at cell division. Cells divide according to a branching process. Following typical parameter values in bacteria and cancer cell populations, we take the mutation rate to zero and the final number of cells to infinity. We prove that almost every site (entry of the nucleotide sequence) is mutated in only a finite number of cells, and these numbers are independent across sites. However independence breaks down for the rare sites which are mutated in a positive fraction of the population. The model is free from the popular but disputed infinite sites assumption. Violations of the infinite sites assumption are widespread while their impact on mutation frequencies is negligible at the scale of population fractions. Some results are generalised to allow for cell death, selection, and site-specific mutation rates. For illustration we estimate mutation rates in a lung adenocarcinoma.

Equations

$$\mathbb{P}\left[V_i^{(r)}=\psi |(v_i)\right]=\{\begin{array}{ll}\mu /3,& \psi \ne v_i;\\ 1-\mu ,& \psi =v_i.\end{array}$$

$$\mathbb{P}[Y_1=j]=\frac{1}{j(j+1)}$$

$$B=\sum _{k=1}^K{Y}_k.$$

$$\mathbb{E}{z}^B=(1-z{)}^{c(z^{-1}-1)}.$$

$$\mathcal{T}={\cup }_{l=0}^{\infty }\{0,1{\}}^l.$$

$${\mathcal{T}}_t:=\left\{x\in \mathcal{T}:\sum _{y\prec x}{A}_y\le t<\sum _{y⪯x}{A}_y\right\}.$$

$$P_{x,t}:=\frac{|\{y\in {\mathcal{T}}_t:x⪯y\}|}{|{\mathcal{T}}_t|}.$$

$$\underset{t\to \infty }{\lim }{P}_{x,t}=P_x:=\prod _{\varnothing \prec y⪯x}{U}_y$$

$${\sigma }_n=\min \{t\ge 0:|{\mathcal{T}}_t|=n\},$$

$$B_i^{n,\mu }=\left|\{x\in {\mathcal{T}}_{{\sigma }_n}:V_i^{\mu }(x)\ne u_i\}\right|$$

$$n\approx {\mu }^{-1}\approx |\mathcal{S}|.$$

$$(B_i^{n,\mu }{)}_{i\in \mathcal{S}}\to (B_i{)}_{i\in \mathcal{S}}$$

$$\sum _{i\in \mathcal{S}}{\delta }_{B_i^{n,\mu }}.$$

$$\frac{1}{|\mathcal{S}|}\sum _{i\in \mathcal{S}}{\delta }_{B_i^{n,\mu }}\to \Lambda$$

$$\mathbb{P}[n^{-1}{B}_i^{n,\mu }>a]$$

$$\approx$$

$${\mu }^{-1}\mathbb{P}[n^{-1}{B}_i^{n,\mu }>a]\to 2(a^{-1}-1).$$

$$\mathbb{E}\sum _{i\in \mathcal{S}}{\delta }_{n^{-1}{B}_i^{n,\mu }}(a,1)\to 2\eta (a^{-1}-1).$$

$$\sum _{i\in \mathcal{S}}{\delta }_{n^{-1}{B}_i^{n,\mu }}\to \sum _{x\in \mathcal{T}\backslash \{\varnothing \}}{M}_x{\delta }_{P_x}$$

$$\mathbb{E}\left[\sum _{x\in \mathcal{T}\backslash \{\varnothing \}}{M}_x{\delta }_{P_x}(a,1)\right]=2\eta (a^{-1}-1),$$

$$\text{Var}\left[\sum _{x\in \mathcal{T}\backslash \{\varnothing \}}{M}_x{\delta }_{P_x}(a,1)\right]\ge 2\eta (a^{-1}-1).$$

$$\mathcal{S}={\mathcal{S}}_{\text{neut}}\cup {\mathcal{S}}_{\text{sel}},$$

$$\mathbb{P}\left[V_i^{(r)}=\psi |(v_i)\right]={\mu }_i^{v_i,\psi }.$$

$$\mu ={\left({\mu }_i^{\chi ,\psi }\right)}_{i\in \mathcal{S};\chi ,\psi \in \mathcal{N}}$$

$$B=\sum _{k=1}^K{Y}_k({\xi }_k).$$

$$\mathbb{E}{z}^B=\exp \left(c(b/a-1)F\left[1,\frac{\lambda }{a-b};1+\frac{\lambda }{a-b};\frac{b/a-z}{1-z}\right]\right)$$

$${\mathcal{G}}_i=\{v\in \mathcal{G}:\forall j\in \mathcal{S},(u_j\ne v_j⟺i=j)\},$$

$${\left(B_i^{n,\mu }\right)}_{i\in \mathcal{S}}\to {\left(\sum _{v\in {\mathcal{G}}_i}{X}_v\right)}_{i\in \mathcal{S}}$$

$$\left(\alpha (u^{\prime })-\beta (u^{\prime }),\alpha (v^{\prime }),\beta (v^{\prime }),\frac{2\alpha (u^{\prime }){\theta }_i^{u_i,v_i}}{\alpha (u^{\prime })-\beta (u^{\prime })}\right).$$

$${\mathcal{S}}_{\text{neut}}=\bigcup _{j\in J}\mathcal{S}(j),$$

$$\frac{1}{|\mathcal{S}|}\sum _{i\in \mathcal{S}}{\delta }_{B_i^{n,\mu }}\stackrel{p}{\to }\sum _{j\in J}q(j)\Lambda (j)$$

Full text

1. 1 Introduction

2. 2 Model

3. 3 Preliminaries

4. 3.1 Luria-Delbrück distribution

5. 3.2 Yule tree

6. 3.3 Mutation frequency notation

7. 3.4 Parameter regime

8. 4 Main results

9. 5 Generalisations

10. 5.1 Model and notation

11. 5.2 Generalised Luria-Delbrück distribution

12. 5.3 Results

13. 5.4 Open problems

14. 6 Mutations at finite numbers

15. 6.1 Counting genomes

16. 6.2 Construction

17. 6.3 Neglecting double mutations

18. 6.4 Convergence of genome counts

19. 6.5 Conditioning on reaching $n$ cells

20. 6.6 Law of large numbers

21. 7 Mutations at positive fractions

22. 7.1 Single site mutation frequencies and an upper bound

23. 7.2 Cell descendant fractions

24. 7.3 Mutation frequencies under the infinite sites assumption

25. 7.4 The infinite sites assumption approximation

26. 7.5 Mean and variance of the site frequency spectrum

27. 8 Infinite sites assumption violations

28. 9 Estimating mutation rates

29. 9.1 Diploid perspective

30. 9.2 A lung adenocarcinoma

1 Introduction

A population of dividing cells with a mutating DNA sequence is ubiquitous in biology. We study a simple model of this process. Starting with one cell, cells divide and die according to a supercritical branching process. As for DNA, we loosely follow classic models from phylogenetics [16, 25]. Each cell contains a sequence of the nucleotides A, C, G, and T, and each site (entry of the sequence) can mutate independently at cell division. We are interested in the sequence distribution when the population reaches many cells.

Let’s discuss a specific motivation. In recent years, cancer genetic data has been made available in great quantities. One especially common type of data consists of mutation frequencies in individual tumours. These data take the form of a vector $(x_{i})_{i\in\mathcal{S}}$, where $i\in\mathcal{S}$ denotes genetic sites and $x_{i}$ is the frequency of cells which are mutated at site $i$. To make sense of such data in terms of tumour evolution, simple mathematical models can be helpful.

Some important works on the topic are [27, 5, 26, 8]. They consider branching process and deterministic models of tumour evolution. They compare theory with data, estimating evolutionary parameters such as mutation rates. A central feature of their theory, and of countless other works, is the so-called infinite sites assumption (ISA). The ISA states that no genetic site can mutate more than once in a tumour’s lifetime. The assumption’s simplicity drives its popularity. However recent statistical analysis of single cell sequencing data [20] shows “widespread violations of the ISA in human cancers”.

For a ‘non-ISA’ model of a growing population of cells, there is in fact a famous example. Luria and Delbrück [22] modelled recurrent mutations in an exponentially growing bacterial population. Subsequent works [21, 18, 14, 17, 19] (and others) adapted Luria and Delbrück’s model to branching processes and calculated mutation frequencies. These works describe only two genetic states, mutated or not mutated, effectively restricting attention to a single genetic site. In [7] we offered an account of one such model, proving limit theorems for mutation times, clone sizes, and mutation frequencies. We then briefly studied an extension to a sequence of genetic sites. Now we offer a self-contained sequel to [7], slightly adapting the model, and aiming for a deeper understanding of the sequence distribution.

In [7] we studied several parameter regimes. In the present work by contrast, we study only one parameter regime which is the most biologically relevant. We take the final number of cells to infinity and the mutation rate to zero with their product finite. This limit is relevant because a detected tumour has around $10^{9}$ cells while the mutation rate per site per cell division is around $10^{-9}$ [15]. This limit is also standard in Luria-Delbrück-type models of bacteria.

Now we introduce our main results. The number of cells mutated at a given site (mutations are defined relative to the initial cell) converges to the Luria-Delbrück distribution. This recovers a well-known result of single site models [21, 14, 19, 17, 7]. So a site is mutated in only a finite number of cells, standing in contrast to the infinite total number of cells. Going beyond [21, 14, 19, 17, 7], we also study the rare event that a site is mutated in a positive fraction of cells. We show that, when appropriately scaled, this fraction of cells follows a power-law distribution.

Across sites, mutation frequencies are asymptotically independent. The independence leads to a many-sites law of large numbers. Specifically, the site frequency spectrum (empirical measure of mutation frequencies) converges to a deterministic measure concentrated at finite cell numbers. At positive fractions of cells, away from the mass concentration, independence breaks down and the site frequency spectrum converges to a Cox process. These results go beyond [10, 27, 5, 8, 7] who only give the expected site frequency spectrum, so our work contributes an appreciation of randomness.

Our results are not all at the same level of generality. For sites mutated in a positive fraction of cells, results are proven for a zero death rate and homogeneous division and mutation rates. For sites mutated in a finite number of cells, results are proven for sequence-dependent death, division, and mutation rates.

We also assess the infinite sites assumption’s validity. Our results say that for typical parameter values, the number of sites to violate the ISA is at least millions, or even billions, in a single tumour. Thus our work agrees with [20]’s statistical analysis of single cell sequencing data which says that ISA violations are widespread. It should be emphasised however that ISA violations do not neccessarily invalidate the ISA. One of our results says that ISA violations do not impact mutation frequencies viewed at the scale of population fractions. Bulk sequencing data, which is the majority of cancer genetic data [8], views mutation frequencies at the scale of population fractions. Therefore our work endorses analyses of bulk sequencing data which are reliant on the ISA, such as [27, 5, 26, 8].

Before commencing the paper, let’s note that there are a wealth of other works on mutations in branching processes. Especially common are infinite alleles models, for example [13, 6, 23, 9], where each individual in the population has an allele which can mutate to alleles never before seen in the population. In an infinite alleles model, a mutation always deletes an individual’s ancestral genetic information. In an infinite sites model on the other hand, a mutation never deletes ancestral genetic information; mutations simply accumulate. The DNA sequence model which we study sits between those extremes.

The paper is structured as follows. In Section 2, we introduce the model in its simplest form. In Section 3, we give notation and preliminary ideas. In Section 4, we present the paper’s main results. In Section 5, we give generalisations and open questions. In Section 6, we prove results on sites mutated in a finite number of cells. In Section 7, we prove results on sites mutated in a positive fraction of cells. In Section 8, we discuss the infinite sites assumption’s validity. In Section 9, we consider data from a lung adenocarcinoma and estimate mutation rates.

2 Model

Here the model is stated in its simplest form. It comprises two parts.

Population dynamics: Starting with one cell, cells divide according to the Yule process. That is, cells divide independently at constant rate. 2. 2.

Genetic information: The set of nucleotides is $\mathcal{N}=\{A,C,G,T\}$. The set of genetic sites is some finite set $\mathcal{S}$. The set of genomes (or DNA sequences) is $\mathcal{G}=\mathcal{N}^{\mathcal{S}}$. Each cell has a genome, i.e. is assigned an element of $\mathcal{G}$. Suppose that a cell with genome $(v_{i})_{i\in\mathcal{S}}\in\mathcal{G}$ divides to give daughter cells with genomes $(V_{i}^{(1)})_{i\in\mathcal{S}}$ and $(V_{i}^{(2)})_{i\in\mathcal{S}}$. Conditional on $(v_{i})$, the $V_{i}^{(r)}$ are independent over $i\in\mathcal{S}$ and $r\in\{1,2\}$, and

[TABLE]

It is also assumed that mutations occur independently for different cell divisions.

The model is generalised to cell death, selection, and nucleotide/site-specific mutation rates in Section 5.

3 Preliminaries

3.1 Luria-Delbrück distribution

Let $(Y_{k})_{k\in\mathbb{N}}$ be an i.i.d. sequence of random variables with

[TABLE]

for $j\in\mathbb{N}$. Let $K$ be an independent Poisson random variable with mean $c$. The Luria-Delbrück distribution with parameter $c$ is defined as the distribution of

[TABLE]

It is commonly seen in its generating function form (e.g. [21, 30])

[TABLE]

The connection between (1) and (2) is made explicit in [7] for example. Although the distribution is named after Luria and Delbrück (due to their groundbreaking work [22]), it was derived by Lea and Coulson [21]. See [30] for a historical review.

The Luria-Delbrück distribution’s power-law tail was derived in [24].

Lemma 3.1**.**

$\lim_{m\rightarrow\infty}m\mathbb{P}\left[B\geq m\right]=c.$**

3.2 Yule tree

The set of all cells to ever exist, following standard notation, is

[TABLE]

A partial ordering, $\prec$, is defined on $\mathcal{T}$. For $x,y\in\mathcal{T}$, $x\prec y$ means that cell $y$ is a descendant of cell $x$. That is, $x\prec y$ if

there are $l_{1},l_{2}\in\mathbb{N}_{0}$ with $l_{1}<l_{2}$ and $x\in\{0,1\}^{l_{1}},y\in\{0,1\}^{l_{2}}$; and 2. 2.

the first $l_{1}$ entries of $y$ agree with the entries of $x$.

Note that $\emptyset\in\mathcal{T}$ and that $\emptyset\prec x$ for any $x\in\mathcal{T}\backslash\{\emptyset\}$. So $\emptyset$ is the initial cell from which all other cells descend. For further notation, write $x\preceq y$ if $x\prec y$ or $x=y$. Also, write $x0$ and $x1$ for the daughters of $x\in\mathcal{T}$; precisely, if $x\in\mathcal{T}$ and $j\in\{0,1\}$, then $xj$ is the element of $\{0,1\}^{l+1}$ whose first $l$ entries are the entries of $x$ and whose last entry is $j$.

Let $(A_{x})_{x\in\mathcal{T}}$ be a family of i.i.d. exponentially distributed random variables with mean $1$. $A_{x}$ is the lifetime of cell $x$. The cells alive at time $t$ are

[TABLE]

The proportion of cells alive at time $t$ which are descendants of cell $x$ (including $x$) is

[TABLE]

Lemma 3.2**.**

For each $x\in\mathcal{T}$,

[TABLE]

almost surely, where

the $U_{y}$ are uniformly distributed on $[0,1]$; 2. 2.

for any $y\in\mathcal{T}$, $U_{y0}+U_{y1}=1$; 3. 3.

$(U_{y0})_{y\in\mathcal{T}}$* is an independent family.*

Lemma 3.2 will be proved in Section 7.

3.3 Mutation frequency notation

When DNA is taken from a tumour, the tumour’s age is unknown, but one may have a rough idea of its size. Therefore we are interested in the cells’ genetic state at the random time

[TABLE]

when the total number of cells reaches some given $n\in\mathbb{N}$.

Write $V^{\mu}(x)=(V_{i}^{\mu}(x))_{i\in\mathcal{S}}\in\mathcal{G}$ for the genome of cell $x\in\mathcal{T}$ (where $\mu$ is the mutation rate). So $(V^{\mu}(x))_{x\in\mathcal{T}}$ is a Markov-process indexed by $\mathcal{T}$ with transition rates given in Section 2. Write $V^{\mu}(\emptyset)=(u_{i})_{i\in\mathcal{S}}$ for the initial cell’s genome. A genetic site is said to be mutated if its nucleotide differs from that of the initial cell. Note that, according to this definition, a site which mutates and then sees a reverse mutation to its initial state is not considered to be mutated. Write

[TABLE]

for the number of cells which are mutated at site $i\in\mathcal{S}$ when there are $n$ cells in total. The quantity (3), and its joint distribution over $\mathcal{S}$, is the key object of our study.

3.4 Parameter regime

The number of cells in a detected tumour may be in the region of $n=10^{9}$, whereas the mutation rate is in the region of $\mu=10^{-9}$ [15]. The human genome’s length is around $|\mathcal{S}|=3\times 10^{9}$. Very roughly,

[TABLE]

Therefore we study the limits:

• •

$n\rightarrow\infty$, $\mu\rightarrow 0$, $n\mu\rightarrow\theta<\infty$;

• •

$n\rightarrow\infty$, $\mu\rightarrow 0$, $n\mu\rightarrow\theta<\infty$, $|\mathcal{S}|\rightarrow\infty$ (sometimes with $|\mathcal{S}|\mu\rightarrow\eta<\infty$).

Remark 3.3**.**

Taking the number of sites to infinity is not to be confused with the infinite sites assumption.

4 Main results

The first result shows that sites are typically mutated in only a finite number of cells, and that these numbers are independent across sites.

Theorem 4.1**.**

As $n\rightarrow\infty$ and $n\mu\rightarrow\theta\in[0,\infty)$,

[TABLE]

in distribution, where the $B_{i}$ are i.i.d. and have Luria-Delbrück distribution with parameter $2\theta$.

Remark 4.2**.**

Taking $|\mathcal{S}|=1$, Theorem 4.1 recovers results of single site models [21, 14, 17, 19, 7].

The site frequency spectrum is a popular summary statistic of genetic data. It is defined as the empirical measure of mutation frequencies:

[TABLE]

The site frequency spectrum sees a law of large numbers.

Theorem 4.3**.**

As $n\rightarrow\infty$, $n\mu\rightarrow\theta\in[0,\infty)$, and $|\mathcal{S}|\rightarrow\infty$,

[TABLE]

in probability, where $\Lambda$ is the Luria-Delbrück distribution with parameter $2\theta$. Convergence is on the space of probability measures on the non-negative integers equipped with the topology of weak convergence.

Theorems 4.1 and 4.3 teach us that almost every site is mutated in only a finite number of cells. What about the rare sites which are mutated in a positive fraction of cells? Heuristically, the Luria-Delbrück distribution’s tail gives the probability that site $i$ is mutated in at least fraction $a$ of cells:

[TABLE]

Approximation (4) is a hand-waving consequence of Theorem 4.1. Approximation (5) is due to Lemma 3.1. The next result offers rigour.

Theorem 4.4**.**

Let $i\in\mathcal{S}$ and $a\in(0,1)$. As $n\rightarrow\infty$ and $n\mu\rightarrow\theta\in[0,\infty)$,

[TABLE]

Theorem 4.4 and linearity of expectation yield the mean site frequency spectrum at positive fractions of the population.

Corollary 4.5**.**

Let $a\in(0,1)$. As $n\rightarrow\infty$, $n\mu\rightarrow\theta\in[0,\infty)$, and $|\mathcal{S}|\mu\rightarrow\eta\in[0,\infty)$,

[TABLE]

The next result gives the distribution of the site frequency spectrum at positive fractions of the population.

Theorem 4.6**.**

As $n\rightarrow\infty$, $n\mu\rightarrow\theta\in[0,\infty)$, and $|\mathcal{S}|\mu\rightarrow\eta\in[0,\infty)$,

[TABLE]

in distribution, with respect to the vague topology on the space of measures on $(0,1]$. That is, the measure applied to a finite collection of closed intervals in $(0,1]$ sees joint convergence. The random variables which appear in the limit are:

• •

$(M_{x})$* is a family of i.i.d. Poisson($\eta$) random variables;*

• •

$(P_{x})$* is from Lemma 3.2 and is independent of $(M_{x})$.*

Remark 4.7**.**

The mean site frequency spectrum, according to Theorem 4.6’s limit, is

[TABLE]

which recovers the limit of Corollary 4.5.

Remark 4.8**.**

The variance of the site frequency spectrum, according to Theorem 4.6’s limit, is bounded below by

[TABLE]

In particular, the coefficient of variation tends to infinity as $a\uparrow 1$.

The details of Remarks 4.7 and 4.8 are given in Section 7.5.

5 Generalisations

Motivated by biological reality, we introduce some generalisations: cell death, selection, and heterogeneous mutation rates.

5.1 Model and notation

Starting with one cell, the cell population grows according to a continuous-time multitype Markov branching process. The types are the genomes, elements of $\mathcal{G}=\mathcal{N}^{\mathcal{S}}$. It will be helpful to classify different types of genetic site. Partition the sites into neutral and selective sites:

[TABLE]

with $\mathcal{S}_{\text{sel}}\not=\emptyset$. For a genome $v=(v_{i})_{i\in\mathcal{S}}$, write $v^{\prime}=(v_{i})_{i\in\mathcal{S}_{\text{sel}}}$ for its restriction to the selective sites. Let $\alpha$ and $\beta$ be functions with domain $\mathcal{N}^{\mathcal{S}_{\text{sel}}}$ and range $[0,\infty)$. A cell with genome $v$ divides at rate $\alpha(v^{\prime})$ (to be replaced by two daughter cells) and dies at rate $\beta(v^{\prime})$.

The initial cell is said to have genome $u$, which is assumed to give a positive growth rate: $\alpha(u^{\prime})>\beta(u^{\prime})$.

Consider a cell with genome $(v_{i})_{i\in\mathcal{S}}$ dividing to give daughter cells with genomes $(V_{i}^{(1)})_{i\in\mathcal{S}}$ and $(V_{i}^{(2)})_{i\in\mathcal{S}}$. Conditional on $(v_{i})$, the $V_{i}^{(r)}$ are independent over $i\in\mathcal{S}$ and $r\in\{1,2\}$, and

[TABLE]

Slightly adapting previous notation, write

[TABLE]

for the collection of mutation rates. Now let’s state the notation for mutation frequencies (for brevity, unlike in Section 3.3, we shall do so in words). Write $B_{i}^{n,\mu}$ for the number of cells which are mutated at site $i$ when $n$ cells are first reached conditioned on the event that $n$ cells are reached.

5.2 Generalised Luria-Delbrück distribution

Let $(\xi_{k})_{k\in\mathbb{N}}$ be an i.i.d. sequence of exponentially distributed random variables with mean $\lambda^{-1}$. Let $(Y_{k}(\cdot))_{k\in\mathbb{N}}$ be an i.i.d. sequence, where $Y_{1}(\cdot)$ is a birth-death branching process with birth and death rates $a$ and $b$ respectively and initial condition $Y_{1}(0)=1$. Let $K$ be a Poisson random variable with mean $c$. The $(\xi_{k})$, $(Y_{k}(\cdot))$, and $K$ are independent. The generalised Luria-Delbrück distribution with parameters ($\lambda,a,b,c$) is defined as the distribution of

[TABLE]

Its generating function

[TABLE]

when $a>b$ is seen in [17, 19, 7]. Here $F$ is Gauss’s hypergeometric function.

Taking parameters $(\lambda,\lambda,0,c)$ recovers the Luria-Delbrück distribution with parameter $c$.

The generalised Luria-Delbrück distribution with parameters ($\lambda,\lambda a,\lambda b,c$), for $\lambda>0$ and $a,b,c\geq 0$, does not depend on $\lambda$. So one could define the distribution with $3$ rather than $4$ parameters. We choose $4$ for a cleaner interpretation of results.

5.3 Results

To begin, Theorem 4.1 is generalised. The genomes whose only difference from the initial cell’s genome is at site $i\in\mathcal{S}$,

[TABLE]

will play a crucial role.

Theorem 5.1**.**

Take $n\rightarrow\infty$ and $n\mu_{i}^{\chi,\psi}\rightarrow\theta_{i}^{\chi,\psi}\in[0,\infty)$ for all $i\in\mathcal{S}$ and $\chi,\psi\in\mathcal{N}$ with $\chi\not=\psi$. Then

[TABLE]

in distribution, where the $X_{v}$ are independent and have generalised Luria-Delbrück distributions with parameters

[TABLE]

In the next result, which generalises Theorem 4.3, we keep the number of selective sites finite while taking the number of neutral sites to infinity. For this limit, mutation rates require consideration. Partition the set of neutral sites:

[TABLE]

such that mutation rates and the initial genome’s nucleotides are homogeneous on $\mathcal{S}(j)$ ($J$ is just some indexing set). Write $\mu^{\chi,\psi}(j)=\mu_{i}^{\chi,\psi}$ for the mutation rates of the sites $i\in\mathcal{S}(j)$. Write $u(j)=u_{i}$ for the initial genome’s nucleotide at the sites $i\in\mathcal{S}(j)$.

Theorem 5.2**.**

Take $n\rightarrow\infty$, $n\mu^{\chi,\psi}(j)\rightarrow\theta^{\chi,\psi}(j)\in[0,\infty)$, $|\mathcal{S}_{\text{neut}}|\rightarrow\infty$, and $|\mathcal{S}(j)|/|\mathcal{S}_{\text{neut}}|\rightarrow q(j)$, for all $j\in J$ and $\chi,\psi\in\mathcal{N}$ with $\chi\not=\psi$. Then

[TABLE]

where the $\Lambda(j)$ are generalised Luria-Delbrück distributions with parameters

[TABLE]

Convergence is on the space of probability measures on the non-negative integers equipped with the topology of weak convergence.

5.4 Open problems

To generalise Theorem 4.6 to a non-zero death rate, selection, and heterogeneous mutation rates, we conjecture the following.

Conjecture 5.3**.**

Take $n\rightarrow\infty$, $n\mu^{\chi,\psi}(j)\rightarrow\theta^{\chi,\psi}(j)\in[0,\infty)$, and $\mu^{\chi,\psi}(j)|\mathcal{S}(j)|\rightarrow\eta^{\chi,\psi}(j)\in[0,\infty)$, for all $j\in J$ and $\chi,\psi\in\mathcal{N}$ with $\chi\not=\psi$. Then

[TABLE]

in distribution, where convergence is in the same sense as Theorem 4.6. The random variables which appear in the limit are:

• •

$(P_{x})_{x\in\mathcal{T}\backslash\{\emptyset\}}$* is distributed as in Lemma 3.2 (and Theorem 4.6), and $P_{\emptyset}=1$;*

• •

$(R_{x})_{x\in\mathcal{T}\backslash\{\emptyset\}}$* is an i.i.d. family of geometric random variables with parameter $(\alpha(u^{\prime})-\beta(u^{\prime}))/(\alpha(u^{\prime})+\beta(u^{\prime}))$, and $R_{\emptyset}$ is independent of $(R_{x})$ but with $R_{\emptyset}\overset{d}{=}R_{0}-1$;*

• •

$(P_{x})$* is independent of $(R_{x})$ if and only if $\beta(u^{\prime})=0$;*

• •

$(M_{x,r})_{x\in\mathcal{T},r\in\mathbb{N}}$* is an i.i.d. family of Poisson random variables with mean $\sum_{j}\sum_{\psi\not=u(j)}\eta^{u(j),\psi}(j)$, independent of $(P_{x},R_{x})$.*

See the Appendix for a heuristic derivation of Conjecture 5.3, which is based on a Yule spinal decomposition of the branching process.

Selection in cancer is a major research topic, and there have been attempts to infer selection from cancer genetic data [4, 26, 8]. Pertinently, Theorem 5.2 and Conjecture 5.3 suggest that selection may not be visible in mutation frequency data, which according to [27] is the case for around $1/3$ of tumours. However we have assumed that the number of selective sites is kept finite. According to [4], there are $3.4\times 10^{4}$ selective sites at which mutations can positively affect growth rate. Thus insight could be gleaned, for example, by taking $|\mathcal{S}_{\text{sel}}|\rightarrow\infty$ with $\mu^{\gamma}|\mathcal{S}_{\text{sel}}|\rightarrow\eta$ for $\gamma\in(0,1]$.

6 Mutations at finite numbers

In this section we prove results on mutations present in only a finite number of cells. In Subsections 6.1 to 6.5 we prove Theorems 4.1 and 5.1 (where $\mathcal{S}$ is finite). In Subsection 6.6 we prove Theorems 4.3 and 5.2 (where $|\mathcal{S}|$ tends to infinity).

6.1 Counting genomes

Assuming mutation rates $\mu=(\mu_{i}^{\chi,\psi})_{i\in\mathcal{S};\chi\psi\in\mathcal{N}}$, write

[TABLE]

for the number of cells with genome $v\in\mathcal{G}$ at time $t\geq 0$. (Recall that the initial condition is $X_{v}^{\mu}(0)=\delta_{u,v}$.) Write

[TABLE]

for the time at which $n\in\mathbb{N}$ cells are reached, and use the convention $\min\emptyset=\infty$.

Recall from (6) that $\mathcal{G}_{i}$ is the subset of genomes with exactly one mutation which is at site $i$. Write

[TABLE]

for the subset of genomes with at least two mutations.

Theorem 6.1**.**

Take $n\rightarrow\infty$ and $n\mu_{i}^{\chi,\psi}\rightarrow\theta_{i}^{\chi,\psi}\in[0,\infty)$ for $i\in\mathcal{S}$ and $\chi,\psi\in\mathcal{N}$ with $\chi\not=\psi$. Then

[TABLE]

in distribution, where the $X_{v}$ are independent and distributed according to:

• •

for $v\in\mathcal{G}_{i}$, $X_{v}$ has generalised Luria-Delbrück distribution with parameters

[TABLE]

• •

for $v\in\mathcal{G}_{\geq 2}$, $X_{v}=0$.

Theorem 6.1 says that cells with at least two mutated sites are non-existent. However simulations and biology tell the opposite story, that cells typically have many mutated sites. This apparent contradiction comes because, while the population size and mutation rate reciprocal converge to infinity, the number of sites is kept finite. So the result only makes sense if one is considering a small subset of the billions of sites.

The mutation frequencies are

[TABLE]

conditional on the event $\{\sigma_{n}^{\mu}<\infty\}$. Therefore Theorem 6.1, via the continuous mapping theorem, implies Theorems 4.1 and 5.1.

Theorem 6.1’s proof is rather lengthy. So, before jumping in with the technical details, let’s give an overview.

In Subsection 6.2 we present a construction of $(X_{v}^{\mu}(\sigma^{\mu}_{n}))_{v\in\mathcal{G}}$. The construction will ultimately illuminate the importance of various subpopulations and the mutations between them. Of particular importance is the primary subpopulation, which is defined as those unmutated cells which have an unbroken lineage of unmutated cells going back to the initial cell. The primary subpopulation, in the limit, grows deterministically and exponentially.

In Subsections 6.3 and 6.4 we show that several events are negligible: primary cells divide to give two mutant daughters; primary cells mutate at multiple sites at once; mutated cells receive further mutations, including backwards mutations. With these events neglected, the situation is pleasingly simplified. The primary subpopulation seeds, as a Poisson process with exponential intensity, single-site mutant subpopulations. The mutant subpopulations grow without further mutations, independently. This gives independent Luria-Delbrück distributions. Finally, in Subsection 6.5 we condition on the event that the population reaches $n$ cells.

Although the proof’s overview may sound simple, the details are less so. The reader will find that the random time $\sigma_{n}^{\mu}$ shoulders a large responsibility for complexity.

6.2 Construction

Additional notation to be used in the proof: for $v\in\mathcal{G}$,

[TABLE]

is the element of $(\mathbb{N}_{0})^{\mathcal{G}}$ denoting that there is one genome $v$ and zero other genomes.

Let

[TABLE]

be a sequence of mutation rates. Assume that

[TABLE]

for $\chi\not=\psi$.

Fix $n\in\mathbb{N}$. For $v,w\in\mathcal{G}$, write

[TABLE]

for the probability that a cell with genome $u$ which divides, gives daughters with genomes $v,w$ (which implies that we have assumed an ordering of the daughters - the first has genome $v$ and the second has genome $w$). Now the construction of $(X_{v}^{\mu_{n}}(\sigma^{\mu_{n}}_{n}))_{v\in\mathcal{G}}$ begins. For the foundational step, introduce the following random variables on a fresh probability space.

[TABLE]

is a birth-death branching process with birth and death rates

[TABLE]

and

[TABLE]

The initial condition $Z^{n}(0)=1$ is assumed. 2. 2.

For $j\in\mathbb{N}$,

[TABLE]

are $\{\emptyset\}\cup(\mathcal{G}\backslash\{u\})^{2}$-valued random variables, with

[TABLE]

and for $v,w\in\mathcal{G}\backslash\{u\}$

[TABLE] 3. 3.

For $v\in\mathcal{G}\backslash\{u\}$ and $j\in\mathbb{N}$,

[TABLE]

is a $(\mathbb{N}_{0})^{\mathcal{G}}$-valued Markov process, with the same transition rates as $(X^{\mu_{n}}_{x}(\cdot))_{x\in\mathcal{G}}$ (defined in (7)) and with the initial condition $\mathcal{Y}^{n}_{v,j}(0)=e_{v}$. 4. 4.

For $v,w\in\mathcal{G}\backslash\{u\}$ and $j\in\mathbb{N}$,

[TABLE]

is a $(\mathbb{N}_{0})^{\mathcal{G}}$-valued Markov process, with the same transition rates as $(X^{\mu_{n}}_{x}(\cdot))_{x\in\mathcal{G}}$ and with the initial condition $\mathcal{Y}^{n}_{v,w,j}(0)=e_{v}+e_{w}$. 5. 5.

For $v\in\mathcal{G}$,

[TABLE]

are Poisson counting processes with rate $1$.

The random variables

[TABLE]

are assumed to be independent ranging over $v,w,j$.

Let’s explain the meaning of the random variables introduced so far. $Z^{n}(\cdot)$ represents the ‘primary’ subpopulation - which we define as the type $u$ cells whose ancestors are all of type $u$. That is to say, there is an unbroken lineage of type $u$ cells between any primary cell and the initial cell. The rate, $\alpha_{n}$, that a primary cell gives birth to another primary cell, is simply the type $u$ division rate multiplied by the probability that no mutation occurs in either daughter cell. The rate, $\beta_{n}$, that a primary cell is removed, is the rate that a type $u$ cell dies plus the rate that a type $u$ cell divides to produce two mutant daughter cells.

The $E_{j}^{n}$ describe what happens at the $j$th downstep in the primary subpopulation trajectory. If $E_{j}^{n}=\emptyset$, then the downstep is a primary cell death. If $E_{j}^{n}=(v,w)$, then the downstep is a primary cell dividing to produce two mutant daughter cells of types $v$ and $w$.

Sometimes a primary cell divides to produce one primary cell and one mutant cell of type $v$. For the $j$th time that this occurs, $\mathcal{Y}^{n}_{v,j}(t)$ is the vector which counts the cells with each genome amongst the descendants of that type $v$ cell, time $t$ after its birth.

Sometimes a primary cell divides to produce two mutant cells of types $v$ and $w$. For the $j$th time that this occurs, $\mathcal{Y}^{n}_{v,w,j}(t)$ is the vector which counts the cells with each genome amongst the descendants of the two mutants time $t$ after their birth.

The $N_{v}(\cdot)$ will soon be rescaled in time to represent the times at which primary cells divide to produce one primary cell and one cell with genome $v$.

The random variables introduced so far, seen together in (9), provide all the necessary ingredients for the construction of $(X_{v}^{\mu_{n}}(\sigma_{n}^{\mu_{n}}))_{v\in\mathcal{G}}$. Now we build upon these founding objects, defining further random variables.

For $v\in\mathcal{G}\backslash\{u\}$ and $t\geq 0$,

[TABLE] 2. 7.

For $j\in\mathbb{N}$ and $v\in\mathcal{G}\backslash\{u\}$,

[TABLE] 3. 8.

[TABLE]

and then for $j>1$, recursively,

[TABLE]

(Here $Z^{n}(t^{-}):=lim_{s\uparrow t}Z^{n}(s)$.) 4. 9.

For $v,w\in\mathcal{G}\backslash\{u\}$,

[TABLE]

and then for $j>1$, recursively,

[TABLE] 5. 10.

For $v,w\in\mathcal{G}\backslash\{u\}$, and $t\geq 0$,

[TABLE]

Let’s explain the meaning of the new random variables. The $K^{n}_{v}(t)$ specify the number of times before time $t$ that primary cells have divided to produce one primary cell and one type $v$ cell. Let’s check that this interpretation makes sense. Conditioned on the trajectory of $Z^{n}(\cdot)$, $K_{v}^{n}(\cdot)$ is certainly a Markov process, and increases by $1$ at rate $2p_{n}(u,v)\alpha(u^{\prime})Z^{n}(t)$ - i.e. the rate at which primary cells divide multiplied by the probability that exactly one daughter cell is primary and one is type $v$.

$S^{n}_{j}$ is the time of the $j$th downstep of the primary subpopulation size. Then $T^{n}_{v,w,j}$ is the time of the $j$th primary cell division which produces two mutant cells of types $v$ and $w$. Note that a primary cell division which produces two mutant cells neccessarily coincides with a downstep in the primary subpopulation size. $K^{n}_{v,w}(t)$ is the number of primary cell divisions before time $t$ which produce cells of types $v$ and $w$.

The reader might question why we have decided to construct the ‘single mutation’ times and the ‘double mutation’ times so differently. The reason for the difference is that single and double mutations will play different roles in the limit, and require different techniques for the proof.

At last the construction reaches its dénouement.

For $t\geq 0$,

[TABLE] 2. 12.

[TABLE]

where $|\cdot|$ is the $l_{1}$-norm on $\mathbb{R}^{\mathcal{G}}$.

Note that $\mathcal{X}^{n}(\cdot)$ has the same distribution as $(X_{v}^{\mu_{n}}(\cdot))_{v\in\mathcal{G}}$; both objects are Markov processes on $(\mathbb{N}_{0})^{\mathcal{G}}$, whose initial conditions and transition rates coincide.

Next we will show that certain elements of the construction converge in distribution. Convergence will sometimes be in the Skorokhod sense. For notation, write $\mathbb{D}(I,R)$ for the space of càdlàg functions from an interval $I\subset[0,\infty)$ to a metric space $R$ (which will always be complete and separable). The space $\mathbb{D}(I,R)$ is equipped with the standard Skorokhod topology. Less standard, we will also consider the space $\mathbb{D}([0,\infty],R)$, which is defined by identification with $\mathbb{D}([0,1],R)$. Let’s be specific. Define $\omega:[0,1]\rightarrow[0,\infty]$ by

[TABLE]

For $f\in\mathbb{D}([0,1],R)$ define $\hat{w}(f)=f\circ\omega^{-1}$. Then the space $\mathbb{D}([0,\infty],R)$ is the image of $\hat{w}$ equipped with the induced topology. Note that according to this definition, for all $z\in\mathbb{D}([0,\infty],R)$, $\lim_{t\rightarrow\infty}z(t)$ exists. In fact we will only consider $z\in\mathbb{D}([0,\infty],R)$ with $\lim_{t\rightarrow\infty}z(t)=z(\infty)$.

Lemma 6.2**.**

As $n\rightarrow\infty$,

[TABLE]

in distribution, on the space $\mathbb{D}([0,\infty],\mathbb{R})$. Here

[TABLE]

is the growth rate of the primary cell population;

[TABLE]

is the large $n$ limit of $\lambda_{n}$; and $Z^{*}(\cdot)$ is a birth-death branching process with birth and death rates $\alpha(u^{\prime})$ and $\beta(u^{\prime})$.

Remark 6.3**.**

The processes of Lemma 6.2 are defined at $t=\infty$. For $n$ large enough that $\lambda_{n}>0$,

[TABLE]

and

[TABLE]

The limits $W^{n}$ and $W^{*}$ exist and are finite almost surely, which is a classic branching process result [3].

Proof of Lemma 6.2.

The transition probabilities of the $\left(e^{-\lambda_{n}t}Z^{n}(t)\right)_{t\in[0,\infty]}$ and $\left(e^{-\lambda t}Z^{*}(t)\right)_{t\in[0,\infty]}$ are well-known [3, 11, 7]. These transition probabilities depend continuously on the birth and death rates, so finite-dimensional convergence is given. To show tightness we shall use Aldous’s criterion [1]. Extend $\omega(s)$ to $s\in[0,2]$ by setting $\omega(s)=\omega(1)=\infty$ for $s\in[1,2]$. Write

[TABLE]

for $s\in[0,2]$. Let $(\rho_{n})$ be a sequence of $[0,1]$-valued stopping times with respect to $(M^{n}(\cdot))$. Let $(\delta_{n})$ be a positive deterministic sequence converging to zero. Then, writing $\mathcal{F}_{\rho_{n}}$ for the sigma-algebra generated by $M^{n}(\cdot)$ up to time $\rho_{n}$,

[TABLE]

where the last equality comes thanks to the fact that

[TABLE]

is a martingale. But

[TABLE]

Now,

[TABLE]

Take $n\rightarrow\infty$ to see that $M^{n}(\rho_{n}+\delta_{n})-M^{n}(\rho_{n})$ converges to zero in $L_{2}$ and hence in probability, thus satisfying Aldous’s criterion. ∎

Lemma 6.4**.**

As $n\rightarrow\infty$,

[TABLE]

in distribution, where $Y_{v,j}(\cdot)$ is a birth-death branching process with birth and death rates $\alpha(v^{\prime})$ and $\beta(v^{\prime})$ and initial condition $Y_{v,j}(0)=1$. Convergence is on the space $\mathbb{D}([0,\infty),\mathbb{R}^{\mathcal{G}})$.

Proof.

It is enough to note that the transition rates converge (see for example page 262 of [12]). ∎

Lemma 6.5**.**

As $n\rightarrow\infty$,

[TABLE]

in distribution, on the space $\mathbb{R}^{\mathbb{N}}$.

Proof.

Note that

[TABLE]

Then

[TABLE]

∎

Remark 6.6**.**

The number $3/2$ which appears in Lemma 6.5 is not special. It only matters that $3/2\in(1,2)$. The relevance of the result will be seen in Section 6.3.

Lemma 6.7**.**

As $n\rightarrow\infty$,

[TABLE]

converges in distribution to

[TABLE]

on

[TABLE]

where the $Z^{*}(\cdot)$ and $Y_{v,j}(\cdot)$ are independent.

Proof.

The convergence seen in Lemmas 6.2, 6.4, and 6.5 is joint convergence over the product space due to independence. ∎

We are yet to say how the random variables in (9) are jointly distributed over $n\in\mathbb{N}$. In fact, the choice of this joint distribution over $n\in\mathbb{N}$ has no relevance to the statement of Theorem 6.1. Hence the choice can be freely made, in a way that streamlines the proof. We assume that:

[TABLE]

almost surely, on the space $\mathbb{D}([0,\infty],\mathbb{R})$;

[TABLE]

almost surely, on the space $\mathbb{D}([0,\infty),\mathbb{R}^{\mathcal{G}})$, for $v\in\mathcal{G}\backslash\{u\}$ and $j\in\mathbb{N}$; and

[TABLE]

almost surely, on the space $\mathbb{R}^{\mathbb{N}}$.

To justify that it is possible to have constructed the random variables in such a way that (11), (12), and (13) hold, one can bring in Skorokhod’s Representation Theorem, to use with Lemma 6.7.

6.3 Neglecting double mutations

Call the event that a primary cell divides to produce two mutant cells a ‘double mutation’. Recall that double mutations are represented by the events $\{E_{j}^{n}=(v,w)\}$, which occur at the times $S_{j}^{n}$ when the primary cell population steps down in size. In order to comment on double mutations, we will first prove a rather crude upper bound for the number of downsteps in the primary cell population trajectory. Write

[TABLE]

for the time at which the primary cell population hits [math] or $n$. Write

[TABLE]

for the number of downsteps in the primary cell population before time $\tau_{n}$.

Lemma 6.8**.**

[TABLE]

almost surely.

Proof.

For each $n\in\mathbb{N}$, let $(R_{j}^{n})_{j\in\mathbb{N}}$ be a sequence of i.i.d. random variables with

[TABLE]

so

[TABLE]

is a random walk, whose distribution matches that of the discrete-time embedded chain of $Z^{n}(\cdot)$. Write

[TABLE]

for the number of steps until the walk hits $n$ or [math]. Then the number of downsteps before hitting $n$ or [math] is

[TABLE]

Therefore we can bound the tail of $D_{n}$’s distribution:

[TABLE]

But $\{\rho_{n}>n^{3/2}\}\subset\{1+\sum_{j=1}^{\lfloor n^{3/2}\rfloor}R_{j}^{n}<n\}$, so

[TABLE]

for some constant $c>0$. Inequality (6.8) holds for large enough $n$ and Inequality (16) is Chebyshev’s inequality. Finally, (17) gives that

[TABLE]

and the result is proven by Borel-Cantelli. ∎

Now it is to be seen that double mutations occurring before time $\tau_{n}$ can be neglected.

Lemma 6.9**.**

Let $v,w\in\mathcal{G}\backslash\{u\}$. As $n\rightarrow\infty$,

[TABLE]

almost surely.

Proof.

From Lemma 6.8, $C:=\sup_{n\in\mathbb{N}}n^{-3/2}D_{n}<\infty$. Then

[TABLE]

By (13) this converges to zero as $n\rightarrow\infty$. ∎

6.4 Convergence of genome counts

The purpose of this section is to show that $\mathcal{X}^{n}(\sigma_{n})$ converges when conditioned on the event $\{W^{*}>0\}$ ($W^{*}$ is defined in Remark 6.3). The times $\tau_{n}$ (defined in (14)) will play the role of a helpful stepping stone in the proof.

Lemma 6.10**.**

Condition on $\{W^{*}>0\}$. Then, almost surely,

there exists $n_{0}$ such that for all $n\geq n_{0}$, $Z^{n}(\tau_{n})=n$; and 2. 2.

$\lim_{n\rightarrow\infty}\tau_{n}=\infty$.

Proof.

To see the first statement, observe that there exists $n_{0}$ such that for all $n\geq n_{0}$, $W^{n}>W^{*}/2>0$. For such $n$, $\lim_{t\rightarrow\infty}Z^{n}(t)=\infty$, and hence $Z^{n}(\cdot)>0$. To see the second statement, suppose for a contradiction that there exists a bounded subsequence $(\tau_{n_{k}})\subset[0,C]$. Then, for large enough $k$,

[TABLE]

The left hand side of the inequality is unbounded over $k$. On the other hand, the right hand side, which does not depend on $k$, is finite thanks to (11). ∎

Lemma 6.11**.**

Condition on $\{W^{*}>0\}$. Suppose that $(a_{n})_{n\in\mathbb{N}}$ is a sequence of real-valued random variables on the same probabiity space as everything else, with

[TABLE]

[TABLE]

for each $n$, and

[TABLE]

almost surely. Then, almost surely,

[TABLE]

where

[TABLE]

Moreover, for $v\in\mathcal{G}_{i}$ and $j\in\mathbb{N}$,

[TABLE]

almost surely, where

[TABLE]

Proof.

Let $t\in\mathbb{R}$. Since $Z^{n}(\tau_{n})=n$,

[TABLE]

Thanks to (11):

for any sequence $(t_{n})$ which converges to infinity, $\lim_{n\rightarrow\infty}e^{-\lambda_{n}t_{n}}Z^{n}(t_{n})=W^{*}$ almost surely; and 2. 2.

$\sup_{n\in\mathbb{N}}\sup_{t\in[0,\infty]}e^{-\lambda_{n}t}Z^{n}(t)<\infty$ almost surely.

So, using dominated convergence,

[TABLE]

almost surely. Note also that

[TABLE]

Then

[TABLE]

Hence, recalling (10),

[TABLE]

converges almost surely to

[TABLE]

because $N_{v}(\cdot)$ is almost surely continuous at any fixed point.

Finally we check convergence of the $a_{n}-T^{n}_{v,j}$. Let $\epsilon>0$. For sufficiently large $n$,

[TABLE]

so

[TABLE]

or equivalently

[TABLE]

The argument is now repeated for an upper bound. For sufficiently large $n$,

[TABLE]

so

[TABLE]

or equivalently

[TABLE]

∎

Lemma 6.12**.**

Condition on $\{W^{*}>0\}$. Suppose that $(a_{n})$ satisfies the conditions of Lemma 6.11. Then, almost surely,

[TABLE]

where the $Y_{v,j}(\cdot)$ are from Lemma 6.4 and the $K_{v}^{*}(\cdot)$ and $T_{v,j}^{*}$ are from Lemma 6.11.

Remark 6.13**.**

By definition, $T^{*}_{v,j}\leq l$ for $j=1,..,K_{v}^{*}(l)$. So the limit in Lemma 6.12 is well defined.

Proof of Lemma 6.12.

Recall that

[TABLE]

The ‘double mutation’ term in (6.12) converges to zero, because

[TABLE]

which converges to zero by Lemma 6.9. As for the ‘single mutation’ term in (6.12), Lemma 6.11 says that the $K^{n}_{v}(a_{n})$ and $a_{n}-T^{n}_{v,j}$ converge to $K^{*}_{v}(l)$ and $l-T^{*}_{v,j}$, while (12) says that the $\mathcal{Y}_{v,j}^{n}(\cdot)$ converge to $e_{v}Y_{v,j}(\cdot)$. ∎

Lemma 6.14**.**

Condition on $\{W^{*}>0\}$.

[TABLE]

almost surely.

Proof.

By Lemma 6.10, for large enough $n$, $Z^{n}(\tau_{n})=n$. So $|\mathcal{X}^{n}(\tau_{n})|\geq n$, and hence $\sigma_{n}\leq\tau_{n}$. Therefore

[TABLE]

Suppose, looking for a contradiction, that

[TABLE]

Take a subsequence with

[TABLE]

Then, by Lemma 6.12,

[TABLE]

converges, and so must be a bounded sequence. However it is also true that, taking $k\rightarrow\infty$,

[TABLE]

which is unbounded. ∎

Lemma 6.15**.**

Condition on $\{W^{*}>0\}$.

[TABLE]

almost surely.

Proof.

By Lemma 6.14, $(\sigma_{n})=(a_{n})$ satisfies the conditions of Lemma 6.11 with $l=0$. Then Lemma 6.12 gives the result. ∎

Let’s look at the limit in Lemma 6.15. For $v\in\mathcal{G}_{i}$, $K^{*}_{v}(0)$ is Poisson distributed with mean $2\lambda^{-1}\alpha(u^{\prime})\theta_{i}^{u_{i},v_{i}}$. Conditional on $K^{*}_{v}(0)$, the times $(-T^{*}_{v,j})_{j=1}^{K^{*}_{v}(0)}$, unordered, are i.i.d. exponentially distributed random variables with mean $\lambda^{-1}$. So

[TABLE]

has generalised Luria-Delbrück distribution with parameters

[TABLE]

Therefore the limit of Lemma 6.15 is a vector of independent generalised Luria-Delbrück distributions:

[TABLE]

where the $X_{v}$ are as stated in Theorem 6.1. To complete the proof of Theorem 6.1 we need to show that conditioning on $\{W^{*}>0\}$ can be translated to conditioning on $\{\sigma_{n}<\infty\}$, which is the subject of the next subsection.

6.5 Conditioning on reaching $n$ cells

In order to connect $\{W^{*}>0\}$ and $\{\sigma_{n}<\infty\}$, the next result is the key. It states that these events are approximately the same for large $n$.

Proposition 6.16**.**

$\lim_{n\rightarrow\infty}\mathbb{P}[W^{*}>0,\sigma_{n}=\infty]=0$, and 2. 2.

$\lim_{n\rightarrow\infty}\mathbb{P}[W^{*}=0,\sigma_{n}<\infty]=0$.

Let’s break the proof of Proposition 6.16 into several lemmas; the idea is that the random variable $W^{n}$ be used as an intermediary.

Lemma 6.17**.**

[TABLE]

Proof.

If $W^{*}>0$, then there exists $n_{0}$, such that for all $n\geq n_{0}$

[TABLE]

So

[TABLE]

Therefore, by dominated convergence,

[TABLE]

∎

Lemma 6.18**.**

[TABLE]

Proof.

If $W^{n}>0$, then $\lim_{t\rightarrow\infty}X^{n}(t)=\infty$, and so $\sigma_{n}<\infty$. ∎

Proof of Part 1 of Proposition 6.16.

[TABLE]

as $n\rightarrow\infty$, by Lemmas 6.17 and 6.18. ∎

The structure for the proof of Part 2 of Proposition 6.16 is much the same as that of Part 1. However the details will require a little extra work.

Lemma 6.19**.**

[TABLE]

Proof.

Let $\epsilon>0$. If $W^{*}=0$, then there exists $n_{0}$ such that for all $n\geq n_{0}$

[TABLE]

So

[TABLE]

almost surely. Then by dominated convergence,

[TABLE]

Meanwhile for each $n$,

[TABLE]

(the distribution of $W^{n}$ is seen in [3, 7]). Therefore

[TABLE]

But $\epsilon>0$ was arbitrary, giving the result. ∎

Lemma 6.20**.**

[TABLE]

Proof.

If the primary population size never reaches $n$ and there are never any mutations, then the total population size never reaches $n$. That is, if $Z^{n}(\tau_{n})=0$, $K_{v}^{n}(\cdot)=0$ and $K^{n}_{v,w}(\cdot)=0$ for all $v,w\in\mathcal{G}\backslash\{u\}$, then

[TABLE]

which means that $\sigma_{n}=\infty$. Equivalently,

[TABLE]

where the equality relies on the fact that $\{Z^{n}(\tau_{n})=0\}\cup\{Z^{n}(\tau_{n})=n\}$ covers the whole probability space. It follows that

[TABLE]

We will show that each term of the right hand side of Inequality (19) converges to zero. Firstly,

[TABLE]

which is the probability that $Z^{n}(\cdot)$, if starting at size $n$, eventually goes extinct; this clearly converges to zero.

Secondly,

[TABLE]

because $p_{n}(u,v)\rightarrow 0$. Hence

[TABLE]

Lastly,

[TABLE]

But $\mathbb{P}[K^{n}_{v,w}(\tau_{n})\not=0]$ converges to zero by Lemma 6.9, while $\mathbb{P}[Z^{n}(\tau_{n})=0]$ converges to $\mathbb{P}[W^{*}=0]>0$ by (11). ∎

Proof of Part 2 of Proposition 6.16.

Just as for Part 1,

[TABLE]

as $n\rightarrow\infty$, by Lemmas 6.19 and 6.20. ∎

Corollary 6.21** (to Proposition 6.16).**

For any sequence of events $(H_{n})_{n\in\mathbb{N}}$,

[TABLE]

if the limit exists.

Proof.

Partition the event $\{H_{n}\cap(W^{*}>0\cup\sigma_{n}<\infty)\}$ in two ways to obtain

[TABLE]

and take $n\rightarrow\infty$. ∎

Finally we are in a position to prove Theorem 6.1.

Proof of Theorem 6.1.

For any $R\subset(\mathbb{N}_{0})^{\mathcal{G}\backslash\{u\}}$,

[TABLE]

where (20) is due to Corollary 6.21 and (21) is due to Lemma 6.15.

∎

6.6 Law of large numbers

Proof of Theorem 5.2 (and Theorem 4.3).

The expected site frequency spectrum is given by

[TABLE]

for $k\in\mathbb{N}_{0}$. The penultimate term of (22) vanishes:

[TABLE]

because $|\mathcal{S}_{\text{sel}}|/|\mathcal{S}|\rightarrow 0$. The last term of (22) can be written as

[TABLE]

where $\mathbb{P}[B^{n,\mu}(j)=k]=\mathbb{P}[B^{n,\mu}_{i}=k]$ for $i\in\mathcal{S}(j)$. But

[TABLE]

while Theorem 5.1 implies that

[TABLE]

Therefore the expected site frequency spectrum converges:

[TABLE]

The variance is

[TABLE]

Because $\mathcal{S}_{\text{sel}}$ and $J$ are finite sets and the random variables are exchangable over $\mathcal{S}(j)$, the maximum is taken over a finite set. Theorem 5.1 says that the covariances converge to zero. ∎

7 Mutations at positive fractions

In this section we return to the basic Yule process setting, proving results on mutations present in a positive fraction of cells. In Subsection 7.1 we prove Theorem 4.4 and also prove an upper bound for mutation frequencies. In 7.2 we prove Lemma 3.2 and another result concerning cell descendant fractions. In 7.3 we determine mutation frequencies under the infinite sites assumption. In 7.4 we show that the infinite sites assumption can offer an approximation for mutation frequencies, concluding the proof of Theorem 4.6. In 7.5 we give details of Remarks 4.7 and 4.8.

7.1 Single site mutation frequencies and an upper bound

Write

[TABLE]

for the cells which are mutated at site $i\in\mathcal{S}$, and

[TABLE]

for their descendants. Recall that $\mathcal{T}_{\sigma_{n}}$ are the cells alive when the total number of cells reaches $n$. Note the inequality

[TABLE]

The goal of this subsection is to prove Theorem 4.4 and the following closely related result, which will later play a crucial role in the proof of Theorem 4.6.

Proposition 7.1**.**

Let $i\in\mathcal{S}$ and $a\in(0,1)$. As $n\rightarrow\infty$ and $n\mu\rightarrow\theta\in[0,\infty)$,

[TABLE]

For this subsection we are always talking about a single site $i\in\mathcal{S}$; for convenience, let’s drop the subscript $i$ from the notation. To begin the proofs of Theorem 4.4 and Proposition 7.1, fix $\mu$, and observe that $(B^{r,\mu})_{r\in\mathbb{N}}$ is a Markov process on the nonnegative integers with transition probabilities

[TABLE]

Here, $j/r$ is the probability that one of the $j$ mutant cells divides, while $\mu/3$ is the probability that a mutant’s daughter reverts to the unmutated state. The process $(\hat{B}^{r,\mu})_{r\in\mathbb{N}}$ has transition probabilities

[TABLE]

The key idea of the proof will be to condition on the number of cells when the first mutant (with respect to site $i$) arises. For this purpose, introduce

[TABLE]

for the total number of cells when the first mutant cell arises. Let

[TABLE]

be the event that the first cell division to see a mutation gives $j$ mutant cells, for $j\in\{1,2\}$.

Lemma 7.2**.**

For $r\in\mathbb{N}$,

[TABLE]

Proof.

The probability that the first $r-2$ cell divisions give no site $i$ mutations multiplied by the probability that the $(r-1)$th cell division gives exactly one mutant daughter is

[TABLE]

Similarly

[TABLE]

Divide by $\mu$ and take $\mu\rightarrow 0$. ∎

The next result gives conditional mutation frequencies.

Lemma 7.3**.**

Let $a>0$. As $n\rightarrow\infty$ and $n\mu\rightarrow\theta\in[0,\infty)$,

[TABLE]

and

[TABLE]

Proof.

Calculating from the transition probabilities (24),

[TABLE]

So

[TABLE]

and hence

[TABLE]

For the rest of the proof we will condition on the event $\{\xi^{\mu}=r,\Xi^{\mu}_{1}\}$. That is, we will consider the processes $(B^{s,\mu})_{s\geq r}$ and $(\hat{B}^{s,\mu})_{s\geq r}$ conditioned on $B^{r,\mu}=\hat{B}^{r,\mu}=1$. Write

[TABLE]

for the conditional expectation. From (26), for $s\geq r$,

[TABLE]

Combining (27) with

[TABLE]

we have that, for $n\geq r$,

[TABLE]

Therefore, as $\mu\rightarrow 0$ and $n\mu\rightarrow\theta$,

[TABLE]

In just the same manner,

[TABLE]

Consider the single mutant cell present when the total number of cells is $r$. Write $D^{n}$ for the number of cells which have descended from this mutant cell when the total number of cells is $n\geq r$. The process $(D^{n},n-D^{n})_{n\geq r}$ is just Polya’s urn. So

[TABLE]

Moreover a well-known result (e.g. [10]) says that, as $n\rightarrow\infty$, $n^{-1}D^{n}$ converges to a Beta random variable. That is, for $a\in(0,1)$,

[TABLE]

which is exactly the limit we wish to show for $n^{-1}B^{n,\mu}$ and $n^{-1}\hat{B}^{n,\mu}$. To show that $n^{-1}D^{n}$, $n^{-1}B^{n,\mu}$, and $n^{-1}\hat{B}^{n,\mu}$ share the same limiting distribution, we will show that their differences converge to zero. The inequality

[TABLE]

gives that

[TABLE]

The inequality

[TABLE]

gives that

[TABLE]

∎

Lemma 7.4**.**

Consider $(\mu_{n})_{n\in\mathbb{N}}$ with $n\mu_{n}\rightarrow\theta\in[0,\infty)$. Then

[TABLE]

where $c>0$ does not depend on $r,j$.

Proof.

The first inequality is immediate. We prove the second. From the transition probabiities (25),

[TABLE]

So, for $s\leq n$,

[TABLE]

where $c_{1}>0$ is a constant which does not depend on $k,n,s$. Now let’s condition on $\{\xi^{\mu_{n}}=r,\Xi^{\mu_{n}}_{j}\}$, again writing $\mathbb{E}_{r,j}$ for the conditional expectation. From (29), for $s\in\{r,..,n\}$,

[TABLE]

This leads to, for $s\in\{r,..,n\}$,

[TABLE]

where $c_{2}$ is a constant which does not depend on $j,n,r,s$. (In the following, $c_{3},..,c_{6}$ will also be constants.) Calculating third moments from the transition probabilities,

[TABLE]

Then

[TABLE]

Hence

[TABLE]

which combined with (30) gives that

[TABLE]

Apply Markov’s inequality to conclude. ∎

Proof of Theorem 4.4 and Proposition 7.1.

[TABLE]

Lemmas 7.2, 7.3, and 7.4, along with the Dominated Convergence Theorem, show that the limit of (31) is

[TABLE]

The same argument works for $\hat{B}^{n,\mu}$. ∎

7.2 Cell descendant fractions

Here we are concerned with the $P_{x,t}$ (the fraction of cells alive at time $t\geq 0$ which are descendants of cell $x\in\mathcal{T}$).

Aldous [2], in a different language to ours, gives a similar result to Lemma 3.2. Rather than adapting his result, we now give a distinct proof of Lemma 3.2.

Proof of Lemma 3.2.

Write

[TABLE]

for the descendants of cell $x\in\mathcal{T}$, and write

[TABLE]

for the descendants of cell $x\in\mathcal{T}$ which are alive at time $t\geq 0$. Observe that

[TABLE]

Hence

[TABLE]

is measurable with respect to the sigma-algebra generated by $(A_{y})_{y\in\mathcal{D}_{x}}$, and has the same distribution as

[TABLE]

It follows that

[TABLE]

almost surely, where $W_{x}\sim$Exp$(1)$; moreover if $x,y\in\mathcal{T}$ are such that $\mathcal{D}_{x}\cap\mathcal{D}_{y}=\emptyset$, then $W_{x}$ and $W_{y}$ are independent. In particular, $W_{x0}$ and $W_{x1}$ are independent. Now,

[TABLE]

almost surely, and $U_{x0}+U_{x1}=1$. A standard calculation shows that $U_{x0}$ is uniformly distributed on $(0,1)$: for $u\in(0,1)$,

[TABLE]

It remains to show independence of the $U_{x0}$. Another standard calculation shows that

[TABLE]

is independent of

[TABLE]

for $(u,v)\in(0,1)\times(0,\infty)$,

[TABLE]

Now fix $l\in\mathbb{N}$. Because $U_{x0}$ and $W_{x0}+W_{x1}$ are measurable with respect to the sigma-algebra generated by $(A_{y})_{y\in\mathcal{D}_{x}\backslash\{x\}}$, we have that

[TABLE]

forms an independent family of random variables.

Finally we complete the proof by induction. Suppose that $(U_{x0})_{x\in\mathcal{T}:|x|<l}$ is an independent family. Observing that for any $x\in\mathcal{T}$,

[TABLE]

we have that $(U_{x0})_{x\in\mathcal{T}:|x|<l}$ is measurable with respect to the sigma-algebra generated by

[TABLE]

Then, thanks to the independence of (32), $(U_{x0})_{x\in\mathcal{T}:|x|\leq l}$ forms an independent family.∎

Next comes a technical result whose value will become apparent in the next subsection.

Lemma 7.5**.**

Let $\epsilon\in(0,1)$. The set

[TABLE]

is almost surely finite.

Proof.

For $t\geq 0$, let $\mathcal{F}_{t}$ be the sigma-algebra generated by $\left(\mathcal{T}_{s}\right)_{s\in[0,t]}$. For $t\geq s\geq 0$, conditional on $\mathcal{F}_{s}$, $\left(P_{y,t}\right)_{y\in\mathcal{T}_{s}}$ is exchangable. So for $y\in\mathcal{T}$,

[TABLE]

Now let $x\in\mathcal{T}\backslash\{\emptyset\}$. We have

[TABLE]

and hence

[TABLE]

That is, $(P_{x,t})_{t\geq 0}$ is a submartingale with respect to $(\mathcal{F}_{t})_{t\geq 0}$. Then by Doob’s inequality,

[TABLE]

But $P_{x}$ is simply a product of $|x|$ independent Uniform$(0,1)$ random variables, where $|x|\in\mathbb{N}$ is the generation of $x$ (that is, $x\in\{0,1\}^{|x|}$). So $\mathbb{E}[(P_{x})^{2}]=3^{-|x|}$. Hence

[TABLE]

Now

[TABLE]

so the Borel-Cantelli lemma concludes the proof. ∎

7.3 Mutation frequencies under the infinite sites assumption

Enumerate the elements of $\mathcal{T}$,

[TABLE]

in such a way that

[TABLE]

Let’s give an example of such an enumeration: map $(x(r))_{r=1}^{l}\in\{0,1\}^{l}\subset\mathcal{T}$ to $2^{l}+\sum_{r=1}^{l}2^{l-r}x(r)$.

Assuming a mutation rate $\mu$, write

[TABLE]

for the first cell (with respect to the enumeration) which sees a mutation at site $i\in\mathcal{S}$.

Remark 7.6**.**

$\phi_{i}^{\mu}$* has geometric distribution:*

[TABLE]

In this subsection we are concerned with $P_{\phi_{i}^{\mu},\sigma_{n}}$, which is the fraction of cells alive at time $\sigma_{n}$ (when $n$ total cells are reached) which are descendants of cell $\phi_{i}^{\mu}$. Phrased another way, $P_{\phi_{i}^{\mu},\sigma_{n}}$ is the fraction of cells alive at time $\sigma_{n}$ which are mutated at site $i$ under the infinite sites assumption.

Next we give an infinite-sites analog of Theorem 4.6.

Proposition 7.7**.**

As $n\rightarrow\infty$, $\mu\rightarrow 0$, and $|\mathcal{S}|\mu\rightarrow\eta\in[0,\infty)$,

[TABLE]

in distribution, with respect to the vague topology on the space of measures on $(0,1]$.

The proof of Proposition 7.7 will require us to count the number of sites which see their first mutation at cell $x\in\mathcal{T}\backslash\{\emptyset\}$ (with respect to the enumeration); write

[TABLE]

Lemma 7.8**.**

As $\mu\rightarrow 0$ and $|\mathcal{S}|\mu\rightarrow\eta\in[0,\infty)$,

[TABLE]

in distribution, where the $M_{x}$ are i.i.d. Poisson($\eta$) random variables.

Proof.

The initial cell is $x_{1}=\emptyset$. The number of sites which mutate in cell $x_{2}$, $M_{x_{2}}^{\mu,\mathcal{S}}$, is binomially distributed with parameters $\mathcal{S}$ and $\mu$. This converges to a Poisson($\eta$) random variable. Now, for induction, suppose that

[TABLE]

in distribution, where the $M_{x}$ are i.i.d. Poisson($\eta$) random variables. Then

[TABLE]

Due to the property (33) of the enumeration, $M_{x_{k+1}}^{\mu,\mathcal{S}}$ conditioned on the event $\left(M_{x_{j}}^{\mu,\mathcal{S}}\right)_{j=2}^{k}=(m_{j})_{j=2}^{k}$ is just a binomial random variable with parameters $|\mathcal{S}|-\sum_{j=2}^{k}m_{j}$ and $\mu$. Therefore (34) converges as required. ∎

Proof of Proposition 7.7.

Fix a sequence of sets of sites $(\mathcal{S}_{n})_{n\in\mathbb{N}}$ and a sequence of mutation rates $(\mu_{n})_{n\in\mathbb{N}}$ with $\mu_{n}|\mathcal{S}_{n}|\rightarrow\eta$. Apply Skorokhod’s Representation Theorem to Lemma 7.8 to obtain random variables $(M_{x}^{n})_{x\in\mathcal{T}\backslash\{\emptyset\},n\in\mathbb{N}}$ and $(M_{x}^{\prime})_{x\in\mathcal{T}\backslash\{\emptyset\}}$ which satisfy

$(M_{x}^{n})_{x\in\mathcal{T}\backslash\{\emptyset\}}\overset{d}{=}(M_{x}^{\mu_{n},\mathcal{S}_{n}})_{x\in\mathcal{T}\backslash\{\emptyset\}}$ for each $n\in\mathbb{N}$; 2. 2.

$(M_{x}^{\prime})_{x\in\mathcal{T}\backslash\{\emptyset\}}\overset{d}{=}(M_{x})_{x\in\mathcal{T}\backslash\{\emptyset\}}$; and 3. 3.

$\lim_{n\rightarrow\infty}(M_{x}^{n})_{x\in\mathcal{T}\backslash\{\emptyset\}}=(M^{\prime}_{x})_{x\in\mathcal{T}\backslash\{\emptyset\}}$ almost surely.

Put the $M^{n}_{x},M^{\prime}_{x}$ on the same probability space as the $P_{x,t},P_{x}$ so that the $M^{n}_{x},M^{\prime}_{x}$ are independent of the $P_{x,t},P_{x}$. Then

[TABLE]

Let $I_{1},..,I_{k}\subset(0,1]$ be closed intervals. Then

[TABLE]

Lemma 3.2 says that $P_{x,\sigma_{n}}$ converges to $P_{x}$; and $P_{x}$ does not lie on the boundaries of the $I_{j}$ with probability one, so the summands of (35) converge pointwise. Meanwhile Lemma 7.5 says that the sum is over a finite subset of $\mathcal{T}\backslash\{\emptyset\}$. ∎

Lemma 7.9**.**

Let $a\in(0,1)$. As $n\rightarrow\infty$, $\mu\rightarrow 0$, and $|\mathcal{S}|\mu\rightarrow\eta\in[0,\infty)$,

[TABLE]

Proof.

First,

[TABLE]

by Fatou’s lemma, Proposition 7.7, and Remark 4.7. Second,

[TABLE]

by the inequality $P_{\phi_{i}^{n},\sigma_{n}}\leq n^{-1}\hat{B}_{i}^{n,\mu}$ and then Proposition 7.1. ∎

7.4 The infinite sites assumption approximation

Lemma 7.10**.**

Let $a\in(0,1)$. As $n\rightarrow\infty$, $n\mu\rightarrow\theta\in[0,\infty)$, and $|\mathcal{S}|\mu\rightarrow\eta\in[0,\infty)$,

[TABLE]

Proof.

The inequalities

[TABLE]

imply that

[TABLE]

Hence

[TABLE]

which converges to zero thanks to Corollary 4.5, Proposition 7.1, and Lemma 7.9. ∎

Proof of Theorem 4.6.

Let $I_{1},..,I_{k}\subset(0,1]$ be closed intervals. Writing $||\cdot||$ for the $l_{1}$-norm on $\mathbb{R}^{k}$,

[TABLE]

due to Lemma 7.10. Therefore

[TABLE]

and

[TABLE]

share the same limiting distribution, if it exists. This limiting distribution, by Proposition 7.7, is

[TABLE]

as required. ∎

7.5 Mean and variance of the site frequency spectrum

Finally, let’s check Remarks 4.7 and 4.8. Note that

[TABLE]

is a sum of i.i.d. mean-$1$ exponentially distributed random variables, which is just a gamma random variable with parameters $|x|$ and $1$. Then

[TABLE]

So

[TABLE]

As for the variance of the site frequency spectrum,

[TABLE]

8 Infinite sites assumption violations

The infinite sites assumption (ISA) is a popular modelling assumption, stating that each genetic site can mutate at most once during the population’s evolution. There are influential and insightful analyses of tumour evolution which rely on the ISA, for example [27, 5, 26, 8]. However, recent statistical analysis of single cell sequencing data shows “widespread violations of the ISA in human cancers” [20]. Thus it is unclear to what extent [27, 5, 26, 8]’s analyses can be trusted. Our studied model of DNA sequence evolution does not use the ISA and invites a theoretical assessment of the ISA’s validity.

Let’s check the prevalence of ISA violations. For simplicity, consider the most basic version of the model, which was introduced in Section 2. Building upon notation of Section 3.2, write

[TABLE]

for the set of ancestors of those cells alive at time $\sigma_{n}$ (when the total number of cells reaches $n$). Write

[TABLE]

for the number of times that site $i$ mutates up to time $\sigma_{n}$. Observe that $X_{i}^{n,\mu}$ is binomially distributed with parameters $2n-2$ and $\mu$. Site $i$ is said to violate the ISA if $X_{i}^{n,\mu}\geq 2$, which occurs with probability

[TABLE]

Then the number of sites to violate the ISA,

[TABLE]

is binomially distributed with parameters $|\mathcal{S}|$ and $p(n,\mu)$. If parameter values are indeed in the region of $n=10^{9}$ and $\mu=10^{-9}$, then the expected proportion of sites to violate the ISA is in the region of $0.5$. This means that the expected number of sites to violate the ISA may be in the billions. Even if very conservative parameter estimates were plugged in, the number of violations is still massive. In fact violations are even more common if one considers cell death. Suppose that cells divide at rate $\alpha$ and die at rate $\beta$. Then to go from a population of $1$ cell to $n$ cells requires approximately $n\alpha/(\alpha-\beta)$ cell divisions, where the factor $\alpha/(\alpha-\beta)$ may be as large as $100$ [5].

Depite the apparent prevalence of ISA violations, our results suggest that their impact on mutation frequencies is negligible at the scale of population fractions. Importantly, bulk sequencing data is only sensitive on the scale of population fractions. Our theoretical work stands in support of the data-driven works [27, 5, 26, 8, 20].

Note however that our model only considers point mutations; it does not, for example, consider deletions of genomic regions, which are thought to be a significant cause of ISA violations [20].

9 Estimating mutation rates

In this section we wish to give the reader a light flavour of mutation frequency data and its relationship to the model. We estimate mutation rates in a lung adenocarcinoma.

9.1 Diploid perspective

Before presenting data, an additional ingredient needs to be considered: ploidy. Normal human cells are diploid. That is, chromosomes come in pairs. Therefore a particular mutation may be present zero, one, or two times in a single cell. It should be said that the story is far more complex in tumours, with chromosomal instability and aneuploidy coming into play. Even so, many tumour samples display an average ploidy not so far from two (for example see Figure (1a) of [28]). We imagine an idealised diploid world.

To illustrate the diploid structure, label the genetic sites as

[TABLE]

for some $L\in\mathbb{N}$. The first coordinate of a site $(i,j)\in\mathcal{S}$ states on which chromosome of a pair the site lies, and the second coordinate refers to the site’s position on the chromosome. Mutations at sites $(1,j)$ and $(2,j)$ are typically not distinguished in data. In the original model set up, mutations were defined as differences to the initial cell’s genome. Let’s slightly improve that definition. Now a genome $v\in\mathcal{G}$ is said to be mutated at site $(i,j)\in\mathcal{S}$ if $v_{i,j}\not=r_{j}$, where $(r_{j})_{j=1}^{L}$ is some reference. Then data is simplistically stated in the model’s language as

[TABLE]

for $j=1,..,L$. That is, the total number of mutations at position $j$ divided by the total number of chromosomes which contain position $j$.

9.2 A lung adenocarcinoma

The mutation frequency data of a lung adenocarcinoma was made available in [29] (499017, Table S2). The data is plotted in Figures 4 and 4. This is just one tumor to illustrate our results. A broader picture of data is seen in [27, 5]. They analysed hundreds of tumors. Around $1/3$ of the tumors were said to have a power-law distribution for mutation frequencies, resembling the one we consider.

Our method to estimate mutation rates is, to a large extent, inspired by [27, 5]. Their attention is restricted to a subset of mutations. They ignore mutations at frequency less than $0.1$, saying that their detection is too unreliable. They ignore mutations above frequency 0.25, in order to neglect mutations present in the initial cell (which are relatively few). We do the same.

Write

[TABLE]

for the number of mutations with frequency in $(a,b)\subset(0.1,0.25)$. Then, adapting Corollary 4.5 to (36), the expected number of mutations with frequency in $(a,b)$ is

[TABLE]

Under different models, [27, 5] derive the same approximation (37). They estimate the mutation rate $\mu$ by applying a linear regression to (37). We simplify matters even further. Our estimator for $\mu$ is

[TABLE]

which (37) says is asymptotically unbiased. Now let’s calculate $\hat{\mu}$ for the data example. The data shows mutations on the exome, which has rough size $|\mathcal{S}|=3\times 10^{8}$ [5]. And the number of mutations in the specified frequency range is $\mathcal{M}(0.1,0.25)=112$. This gives

[TABLE]

Next let’s consider mutation rate heterogeneity. Write $\mu_{\chi}$ for the rate that nucleotide $\chi\in\mathcal{N}$ mutates. Partition the genetic sites:

[TABLE]

where

[TABLE]

is the set of sites which are represented by nucleotide $\chi$ in the initial cell. Just as before,

[TABLE]

is an unbiased estimator for $\mu_{\chi}$. The data gives

[TABLE]

This method could easily be extended to offer more detail, for example to estimate the rate at which nucleotide $A$ mutates to $C$ or to estimate mutation rates on different chromosomes.

The just presented statistical analysis is of course simple and brief. We recommend [8] for a far more comprehensive statistical analysis of mutation frequency data. However their infinite-sites framework does not consider mutation rate heterogeneity.

Appendix

A heuristic ‘proof’ of Conjecture 5.3 is given.

First we argue that, in the conjecture’s limit, selection is unimportant. Write

[TABLE]

for the set of genomes which are mutated at a selective site. Write

[TABLE]

for the proportion of cells at time $t\geq 0$ whose genomes are mutated at a selective site. Then, according to Theorem 6.1,

[TABLE]

in probability. Therefore we neglect selection.

Cells divide and die at rates $\alpha(u^{\prime})$ and $\beta(u^{\prime})$, which we now abbreviate to $\alpha$ and $\beta$. Some cells have an ultimately surviving lineage of descendants. Other cells eventually have no surviving descendants. Name these cells immortal and mortal respectively. In a supercritical birth-death branching process, it is well-known (eg. [10]) that the immortal cells grow as a Yule process and the mortal cells grow as a subcritical branching process. An immortal cell divides to produce two immortal cells at rate $\alpha-\beta$, or it divides to produce one immortal and one mortal cell at rate $2\beta$. A mortal cell divides at rate $\beta$ to produce two mortal cells, or it dies at rate $\alpha$. Because the process is conditioned to reach a large population size, let’s assume that the initial cell is immortal.

The notation of Section 3.2, $\mathcal{T}=\cup_{l=0}^{\infty}\{0,1\}^{l}$ and its partial ordering $\prec$, will be used to represent the immortal cells. Let $(A_{x})_{x\in\mathcal{T}}$ be i.i.d. Exp($\alpha-\beta$) random variables, which represent the times for immortal cells to divide to produce two immortal cells. The immortal cells at time $t\geq 0$ are

[TABLE]

The immortal descendants of $x\in\mathcal{T}$ are

[TABLE]

The number of immortal descendants of cell $x$ at time $t$ is

[TABLE]

Let $\left((R_{x}(t))_{t\geq 0}\right)_{x\in\mathcal{T}}$ be i.i.d. Poisson processes with rate $2\beta$. Write $R_{x,r}=\min\{t\geq 0:R_{x}(t)=r\}$ for $r=1,..,R_{x}(A_{x})$. Then the seeding times of mortal cells are

[TABLE]

Each seeding event initiates a subpopulation of mortal cells; let $(Y_{x,r}(t))_{t\geq 0}$ be i.i.d. birth-death branching processes with birth and death rates $\beta$ and $\alpha$. Then the number of mortal descendants of $x$ at time $t$ is

[TABLE]

The number of descendants of $x$ at time $t$ is

[TABLE]

The next result shows the long-term proportion of a cell’s descendants which are immortal. The result is a basic consequence of classic branching process theory [3], and was mentioned in its specific form by [10].

Lemma A.1**.**

There is $c\in(0,\infty)$ with

[TABLE]

almost surely.

We use Lemma A.1 to see the number of descendants of a cell as a proportion of the total population.

Lemma A.2**.**

For $x\in\mathcal{T}\backslash\{\emptyset\}$,

[TABLE]

almost surely, where the $P_{x}$ are as in Lemma 3.2.

Proof.

By Lemma 3.2 and Lemma A.1,

[TABLE]

converges to the required limit. ∎

Let’s look at mutations. In the proof of Theorem 4.6 it was shown that the number of new mutations to arise at a cell’s birth is approximately Poisson. Here, with heterogeneous mutation rates, the number of new mutations to arise at a cell’s birth is approximately Poisson with mean

[TABLE]

Each $x\in\mathcal{T}\backslash\{\emptyset\}$ witnesses $1+R_{x}(A_{x})$ cell divisions, while $\emptyset$ witnesses $R_{\emptyset}(A_{\emptyset})$ cell divisions (one less because there is not a cell division associated to the initiation of $\emptyset$). So the number of new mutations to arise at $x$ is

[TABLE]

where the $M_{x,r}$ are i.i.d. Poisson random variables with mean $\eta$. In the proof of Theorem 4.6 it was also shown that a mutation which arises in cell $x$ will have approximate frequency $P_{x}$. Here, thanks to Lemma A.2, the situation appears identical. It only remains to discuss mutations arising in mortal cells. Any subpopulation of cells which descended from a mortal cell must eventually die out. Hence mutations arising in mortal cells are negligible when compared to the infinite total population size.

Declarations of interest

None.

Acknowledgements

We thank Michael Nicholson, Trevor Graham, and Marc Williams for inspiring discussions. We thank two anonymous referees for numerous helpful corrections and suggestions. David Cheek was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016508/01), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.