# Single-cell mutational burden distributions in birth-death processes

Christo Morison1 1School of Mathematical Sciences, Queen Mary University of
London, London, E1 4NS United Kingdom. 2Group of Theoretical Biology, The
State Key Laboratory of Biocontrol, School of Life Sciences, Sun Yat-sen
University, Guangzhou, 510060 China.
∗Corresponding author<EMAIL_ADDRESS>Dudley Stark1 & Weini Huang1,2,∗
1School of Mathematical Sciences, Queen Mary University of London, London, E1
4NS United Kingdom. 2Group of Theoretical Biology, The State Key Laboratory of
Biocontrol, School of Life Sciences, Sun Yat-sen University, Guangzhou, 510060
China.
∗Corresponding author<EMAIL_ADDRESS>
###### Abstract

Genetic mutations are footprints of cancer evolution and reveal critical
dynamic parameters of tumour growth, which otherwise are hard to measure in
vivo. The mutation accumulation in tumour cell populations has been described
by various statistics, such as site frequency spectra (SFS) from bulk or
single-cell data, as well as single-cell division distributions (DD) and
mutational burden distributions (MBD). An integrated understanding of these
distributions obtained from different sequencing information is important to
illuminate the ecological and evolutionary dynamics of tumours, and requires
novel mathematical and computational tools. We introduce dynamical matrices to
analyse and unite the SFS, DD and MBD based on a birth-death process. Using
the Markov nature of the model, we derive recurrence relations for the
expectations of these three distributions. While recovering classic exact
results in pure-birth cases for the SFS and the DD through our new framework,
we also derive a new expression for the MBD as well as approximations for all
three distributions when death is introduced, confirming our results with
stochastic simulations. Moreover, we demonstrate a natural link between the
SFS and the single-cell MBD, and show that the MBD can be regenerated through
the DD. Surprisingly, the single-cell MBD is mainly driven by the
stochasticity arising in the DD, rather than the extra stochasticity in the
number of mutations at each cell division.

Keywords: Birth-death process, single-cell mutational burden distribution,
stochastic dynamics of mutation accumulation, site frequency spectrum,
infinite sites approximation, cancer evolution.

## Introduction

Somatic mutations are important for the evolution of biological systems with
clonal reproduction, including the development from healthy tissues to cancer
[1, 2]. While less is known about the somatic mutation rates in clonal species
such as plants and corals, they have been studied extensively in human
tissues. Healthy cells may accumulate in the order of $1$ to $2$ mutations per
cell per division, which is directly observable in early development [3, 4,
5]. The mutational rate of tumour cells is often thought to be higher, which
can be caused for example by genomic instability [6, 7, 8]. This large number
of mutations accumulated in tumours serves as a genetic footprint to reveal
their evolutionary history. Since the majority of these mutations are neutral
[9], not impacting the fitness of a cell compared to its parental cell,
neutral theory has been used to explain mutational patterns in many patient
samples across different tumour types [10, 11, 12]. These measurements often
demonstrate an early expansion of tumour populations, wherein driver mutations
are clonal and the intratumour heterogeneity arises from neutral passenger
mutations accumulated after cancer initiation. Although clonal interference,
where cells carrying different sets of driver mutations intercompete, is a
likely alternative scenario especially in large populations [13, 14], here we
focus on a further understanding of mutation accumulation under neutral
selection.

Distributions of genetic heterogeneity under neutral selection have been
studied in population genetics for over half a century [15, 16]. One such
measure is the site frequency spectrum (SFS), which describes the frequencies
of mutations in a population [17]. Because the SFS deals with population-level
information, it can be compared to bulk genomic data [12, 18]. For an
exponentially-growing population, a rescaling of the SFS, the cumulative
variant allele frequency spectrum, has been shown to follow a $1/f^{2}$ power
law, for $f$ the frequency of a mutation in the population [19, 20, 18, 12].
More recently, exact expressions for the SFS were found under the assumption
of neutral evolution [21].

The advent of single-cell data [22, 23] opens the door for combining bulk and
single-cell data to understand mutation accumulation [24, 25]. There is great
need for new mathematical and computational machinery to cope with single-cell
data, which provides different mutational distributions beyond the SFS. The
number of unique mutations in the population, also known as the overall tumour
mutational burden (TMB) [26], has been studied both in a single tumour [21]
and (its distribution) between tumours [27, 28]. However, the distribution of
mutational burdens between cells, the so-called single-cell mutational burden
distribution (MBD), has only recently been experimentally observable through
single-cell sequencing. Understanding the MBD may further help in inferring
important evolutionary parameters, determining the growth history of the
tumour and the level of selection at play, with neutral selection as a
baseline with which to compare. Using data from healthy haematopoietic stem
cells and oesophageal epithelial cells, Moeller et al. showed that analysis of
single-cell and bulk data complement each other and narrowed down the
parameter inference of the mutation rate and stem cell population size [25].
More specifically, the mean and variance of the MBD for a growing population
were derived and used to estimate the underlying mutation rate [25]. However,
the exact analytical shape of the MBD has not yet been explicitly found.

Another distribution that can be constructed through single-cell data is the
division distribution (DD) of individual cells, describing the frequencies of
cells having a given number of divisions in their history. It is perhaps best
seen through a network theory lens: birth-death processes can be viewed as
growing trees, where cell division is represented by a leaf bifurcating into
two leaves, and cell death is the removal, or pruning, of a leaf. A pure-
birth, or Yule, process can then generate phylogenetic trees [29], where
results from graph theory are applied to understand speciation models. The DD
has been solved for the Yule process [30]. There is a natural connection
between the DD and MBD; indeed, using a birth-death model, we recover this DD
result and use it to better understand the somatic evolution of single cells
within a population via the MBD.

We introduce a new framework via dynamical matrices to investigate mutation
accumulation in a birth-death process and explain how key quantities such as
the SFS, DD and MBD are obtained from these mutational matrices. This
framework allows us to derive exact solutions of these distributions by
recurrence relations in the pure-birth case, as well as first-order
approximations when death is introduced, which hold in the low-death and
large-population limits. We compare our solutions for the SFS and DD to known
results in population genetics and all three distributions to our stochastic
simulations. We find that the MBD can be generated via the DD and the mean
mutation rate per cell division, independently of the stochasticity in the
number of mutations per cell division.

## Model

### Mutation accumulation in a birth-death process

In a birth-death process where either a cell divides with probability $\beta$
or dies with probability $\delta=1-\beta$, the population size $N_{i}$ at time
step $i$ can be described by a discrete-time Markov chain (Figure 1a). The
state space in this Markov chain is the finite integer set $\\{0,\dots,N\\}$ ,
where $N$ is the largest possible population size. In some cases, we are
interested in the limit $N\to\infty$.

Often, the stochastic birth-death process explicitly involves a continuous
time parameter $t$ instead of a discrete step count $i$. This allows for rates
to be considered instead of probabilities; however, as long as events are
assumed not to be simultaneous, these two schemes can be mapped to one another
by choosing a distribution of times between events. Most often, events are
assumed to be exponentially distributed, and thus their frequency grows with
the population size.

Figure 1: (a) Discrete-time Markov chain description of the population size.
(b) Growing binary tree representation of an example realisation of the birth-
death process with mutations described in the main text, with birth
probability $\beta=2/3$, death probability $\delta=1/3$, mutational mean
$\mu=2$, initial population $N_{0}=1$ and final population $N_{12}=5$. Cells
that are crossed out have died. Edges are labelled by the number of new
mutations occurring during that division. Leaves (living cells) are labelled
by their mutational burden, which is equal to the sum of the edges that
connect them to the root, or the mutation-free progenitor cell. (c) The
population size $N_{i}$ (solid line) and number of unique mutations $M_{i}$
(dashed line) versus the step count $i$ for the example realisation in (b).

Here, we focus on the growing-population case $\beta>\delta$ and assume
$N_{0}=1$ unless otherwise mentioned. In this case, the birth-death process is
a growing rooted binary tree, where the root is the lone progenitor cell
(assumed to be mutation-free, as any of its mutations will be clonal in the
population), leaves are living cells, and pruned leaves are dead cells. See
Figure 1b for an example realisation of such a process. Novel mutations are
accumulated during cell divisions and old mutations may be lost in cell death.
When a cell divides, its daughter cells inherit all mutations carried by the
mother cell and acquire a random number of new mutations drawn from a Poisson
distribution with mean $\mu$: $U_{1},U_{2}\sim\text{Pois}(\mu)$, where the
indices refer to the two daughter cells. These mutations are unique under the
infinite sites approximation, where the probability of two point mutations
occurring at the same location along the large genome is vanishingly small
[31]. The point mutations occurring during each duplication of the genome are
thought of independent from each other and thus their number follows a Poisson
distribution [32, 33].

We are most interested in three key quantities of this birth-death process
with mutations at each step $i$: (i) the site frequency spectrum (SFS)
$\\{S_{j,i}\\}_{j}$, whose elements $S_{j,i}$ denote the number of mutations
which occur $j$ times in the population [17]; (ii) the single-cell mutational
burden distribution (MBD) $\\{B_{k,i}\\}_{k}$, whose elements $B_{k,i}$ are
the number of cells having a mutational burden of $k$ [25]; and (iii) the
division distribution (DD) $\\{D_{\ell,i}\\}_{\ell}$, whose elements
$D_{\ell,i}$ indicate the number of cells having undergone $\ell$ divisions
during the process. Note that $D_{\ell,i}$ is also the number of leaves lying
at a distance of $\ell$ edges from the root in the growing tree framework.

While the importance of the SFS and the DD have been investigated in growing
populations [29, 19, 20, 12, 18, 21], we are interested in the relationship
between them and how it can help us understand the MBD. Here, we introduce a
novel discrete-time framework to demonstrate the symmetry between these
distributions. The number of unique mutations at step $i$ is
$M_{i}=\sum_{j}S_{j,i}$, and the population size is
$N_{i}=\sum_{k}B_{k,i}=\sum_{\ell}D_{\ell,i}$, both of which are plotted in
Figure 1c for the example found in Figure 1b. Thus, the MBD and the DD form
partitions of the number of cells in a way similar to the SFS partitioning the
number of unique mutations. Next, we introduce dynamical matrices to connect
those quantities arising from the same population of individual cells.

### Dynamical matrices to unite the SFS, DD and MBD

We consider a collection of matrices $Y_{i}$, where the rows refer to cells
and the columns refer to mutations, known as genotype matrices or SNP (single
nucleotide polymorphism) matrices in bioinformatics [34]. Our matrices are
dynamical in that their entries are updated at each step by a binary filling
in the following manner: the $(n,m)$th entry of the matrix is equal to $1$ if
the $n$th cell possesses the $m$th mutation at step $i$ and equal to $0$
otherwise. When a cell dies, its row is removed from the matrix. Figure 2a
shows an example of the matrix $Y_{i}$ associated to the tree example in
Figure 1b. We extend the concept of genotype matrices by marking mutations
arising during a single (past) division by grey shaded areas. Note that if no
mutations arise during a division, the corresponding $0$ entry would still be
shaded in grey, as this shading tracks the division burden. These matrices are
also how mutational data can be stored in stochastic simulations.

Figure 2: The matrix framework described in the main text, where $i$ refers to
the Markov step count (in this example, $i=12$). (a) The matrix $Y_{i}$
corresponding to the example realisation of the birth-death process depicted
in Figure 1b, where entry $(n,m)$ is $1$ if cell $n$ possesses mutation $m$
and $0$ otherwise. Grey shaded mutational entries arose during the same
division and thus always occur together in the descendants of their
progenitor. For example, mutations $m-1$, $m$ and $m+1$ in cell $n$ (shown
outlined in black) were generated in the same past division. (In Figure 1b, we
can determine by inspection this ancestor cell to be the one that later
divided into two cells with mutational burdens of $6$.) The row sum of the
entries of $Y_{i}$ is shown in orange, the column sum of the number of grey
areas is in pink and the column sum of the entries of $Y_{i}$ is in blue. (b)
The site frequency spectrum (SFS) $\\{S_{j,i}\\}_{j}$: a histogram of the
bottommost (orange) vector of (a). (c) The division distribution (DD)
$\\{D_{\ell,i}\\}_{\ell}$: a histogram of the middle (pink) vector of (a). (d)
The single-cell mutational burden distribution (MBD) $\\{B_{k,i}\\}_{k}$: a
histogram of the rightmost (blue) vector of (a), or equivalently of the sums
of the weights of the edges in Figure 1b from the root to the leaves.

We can obtain the distributions of our key quantities, the SFS, DD and MBD
(Figures 2b–d), from our dynamical mutation matrices (Figure 2a). For each
mutation (column), the number of cells carrying this mutation is the row sum
of the entries of $Y_{i}$ (orange vector in Figure 2a). Thus, the histogram of
this vector is the SFS at step $i$. For each cell (row), the number of
divisions that the cell has undergone is the column sum of the number of grey
areas (pink vector in Figure 2a), and the number of mutations in this cell is
the column sum of the entries of $Y_{i}$ (blue vector in Figure 2a).
Correspondingly, the histograms of these two vectors lead to the other
distributions obtainable from single-cell information: the DD and the MBD.

The symmetry provided by this mutation matrix $Y_{i}$ gives rise to the
following relationship between the site frequency spectrum and the single-cell
mutational burden distribution:

$\sum_{j=1}^{N_{i}}jS_{j,i}=\sum_{k=0}^{M_{i}}kB_{k,i}.$ (1)

We call this quantity the number of mutational occurrences: that is, the sum
of the entries of $Y_{i}$.

Symbol |  Description
---|---
$\beta$ |  Birth probability (Markov transition rate)
$\delta$ |  Death probability (Markov transition rate)
$N_{i}$ |  Population size at step $i$
$N$ |  Maximal population size
$\mu$ |  Mutational mean: mean number of mutations acquired per division per daughter cell
$M_{i}$ |  Number of unique mutations (tumour mutational burden, TMB) at step $i$
$\\{S_{j,i}\\}_{j}$ |  Site frequency spectrum (SFS), where $S_{j,i}$ is the number of mutations occurring $j$ times in the population at step $i$
$\\{B_{k,i}\\}_{k}$ |  Single-cell mutational burden distribution (MBD), where $B_{k,i}$ is the number of cells with $k$ mutations at step $i$
$\\{D_{\ell,i}\\}_{\ell}$ |  Division distribution (DD), where $D_{\ell,i}$ is the number of cells having undergone $\ell$ divisions at step $i$
$Y_{i}$ |  Mutational matrix at step $i$: entry $(n,m)$ is $1$ if cell $n$ possesses mutation $m$, $0$ otherwise
$\genfrac{[}{]}{0.0pt}{}{i}{\ell}$ |  Unsigned Stirling number of the first kind with indices $1\leq\ell\leq i$
Table 1: Notation used in this manuscript.

## Results

Our primary approach for deriving the distributions of our key quantities from
the discrete-time model is as follows. We use the law of total expectation
($\mathbb{E}\left[X\right]=\mathbb{E}\left[\mathbb{E}\left[X\,\middle|\,Y\right]\right]$,
for any random variables $X$ and $Y$) to equate an expected quantity at step
$i+1$ to a conditional expectation. This usually is a function of the quantity
at step $i$ conditional on knowledge at step $i$, as earlier knowledge is
never needed due to the Markov nature of the model. From this, we derive a
recurrence relation for the expected values of our desired quantity, which can
be solved.

We first note that conditioning on the survival of the entire population plays
a role in all of our subsequent expected values. In the pure-birth case, the
population is deterministic and equal to $N_{i}=i+1$. Once death is included,
however, the population size becomes a random variable. For the birth-death
chain of Figure 1a, its expected value at step $i$, both conditioned on
survival and not, can be exactly calculated, which is done in Proposition S2
of the Supplementary Information. Figure S1 shows that the expected population
size conditioned on survival can be linearly approximated by
$\mathbb{E}\left[N_{i}\right]\simeq(\beta-\delta)i+1$, valid for low death,
since in this limit the expected gain in population in one step is
$\beta-\delta$. All of our ensuing expectations are conditioned on survival
and the initial condition $N_{0}=1$, which we will omit from our notation for
brevity.

### Site frequency spectrum

With the recurrence relation method outlined above, we can formally derive the
pure-birth ($\beta=1$) site frequency spectrum. When there is no death,
$N_{i}=i+1$, and thus the expected number of $j$-abundant mutations in the
dividing cell at step $i$ is $jS_{j,i}/(i+1)$. After the division in step $i$,
any $j$-abundant mutations in the dividing cell will become $(j+1)$-abundant
and thus no longer contribute to the $j$-site. Similarly, $(j-1)$-abundant
mutations in the dividing cell now contribute to the $j$-site. We therefore
have

$\mathbb{E}\left[S_{j,i+1}\right]=\mathbb{E}\left[\mathbb{E}\left[S_{j,i+1}\,\middle|\,\\{S_{j^{\prime},i}\\}_{j^{\prime}}\right]\right]=\mathbb{E}\left[S_{j,i}-\frac{jS_{j,i}}{i+1}+\frac{(j-1)S_{j-1,i}}{i+1}+\left(U_{1}+U_{2}\right)\delta_{1,j}\right],$
(2)

for $\delta_{\cdot,\cdot}$ the Kronecker delta symbol, whose source term
arises from the new $1$-abundant mutations occurring during division, which
are drawn from a Poisson distribution: $U_{1},U_{2}\sim\text{Pois}(\mu)$. We
make the change of variables $Q_{j}=\mathbb{E}\left[S_{j,i}\right]/(i+1)$,
whose lack of dependence on $i$ can be computed by brute force from equation
(2), absorbing the source term into the boundary condition $Q_{1}=\mu$. Using
the linearity of expectation, equation (2) becomes simply
$(j+1)Q_{j}=(j-1)Q_{j-1}$, which telescopes to obtain the known result

$\mathbb{E}\left[S_{j,i}\right]=\frac{2\mu(i+1)}{j(j+1)}.$ (3)

An identical procedure can be applied in the birth-death case to recover the
large-population expected SFS, though a difficulty here is that now $N_{i}$
becomes a random variable itself. Since to first order
$\mathbb{E}\left[A/B\right]\simeq\mathbb{E}\left[A\right]/\mathbb{E}\left[B\right]$
(when $A$ and nonzero $B$ are close to their expected values), our first-order
approximation (whose derivation is found in Proposition S4 of the
Supplementary Information), which matches the $N\to\infty$ result from [21],
is

$\mathbb{E}\left[S_{j,i}\right]\simeq\sum_{j^{\prime}=0}^{\infty}\frac{2\mu(\delta/\beta)^{j^{\prime}}}{(j+j^{\prime})(j+j^{\prime}+1)}\mathbb{E}\left[N_{i}\right],$
(4)

where all expectations are conditioned on non-extinction of the whole
population [21]. In the limit of low death ($\delta\ll\beta$), this
approximation is sound, as then the variance in population size is small (see
Proposition S6 for details).

### Division distribution

The expected division distribution in the pure-birth case can be obtained in a
similar manner as the site frequency spectrum. The probability of selecting a
cell with $\ell$ divisions in its history is $D_{\ell,i}/(i+1)$, which will
then no longer contribute to $D_{\ell,i+1}$, dividing into two cells with one
more division in their history than before. The law of total expectation then
becomes

$\mathbb{E}\left[D_{\ell,i+1}\right]=\mathbb{E}\left[\mathbb{E}\left[D_{\ell,i+1}\,\middle|\,\\{D_{\ell^{\prime},i}\\}_{\ell^{\prime}\leq
i}\right]\right]=\mathbb{E}\left[D_{\ell,i}-\frac{D_{\ell,i}}{i+1}+\frac{2D_{\ell-1,i}}{i+1}\right],$
(5)

which can be solved to recover the result from [30]:

$\mathbb{E}\left[D_{\ell,i}\right]=\genfrac{[}{]}{0.0pt}{}{i}{\ell}\frac{2^{\ell}}{i!},$
(6)

where the unsigned Stirling numbers of the first kind
$\genfrac{[}{]}{0.0pt}{}{i}{\ell}$ are defined by the relation

$\genfrac{[}{]}{0.0pt}{}{i+1}{\ell}=i\genfrac{[}{]}{0.0pt}{}{i}{\ell}+\genfrac{[}{]}{0.0pt}{}{i}{\ell-1}\quad\text{
for }1\leq\ell\leq i,$ (7)

with boundary conditions $\genfrac{[}{]}{0.0pt}{}{1}{1}=1$ and
$\genfrac{[}{]}{0.0pt}{}{i}{\ell}=0$ if $\ell>i$ or $\ell=0$. Results from
stochastic simulations using a Gillespie algorithm agree well with this
prediction (see Figure 3).

Figure 3: Average (solid dark pink line) of $200$ simulation realisations
(representatives in solid pale pink lines) of the division distribution (DD)
for a pure-birth process up to final population size $N=10^{3}$, along with
the predicted expected DD obtained from equation (6) (dashed black line).

Similarly to with the SFS, in the $\delta>0$ case we can obtain a first-order
approximation (see Proposition S5) for the DD:

$\mathbb{E}\left[D_{\ell,i}\right]\simeq\frac{\genfrac{[}{]}{0.0pt}{}{i}{\ell}2^{\ell}\left(1-\delta/\beta\right)^{-\ell}}{\sum_{\ell^{\prime}=1}^{i}\genfrac{[}{]}{0.0pt}{}{i}{\ell^{\prime}}2^{\ell^{\prime}}\left(1-\delta/\beta\right)^{-\ell^{\prime}}}\mathbb{E}\left[N_{i}\right],$
(8)

wherein it is evident that the division distribution partitions the population
size (since summing the fraction on the right-hand side over $\ell$ gives
unity) and that this partitioning is orchestrated by the functions
$f(i,\ell)=\genfrac{[}{]}{0.0pt}{}{i}{\ell}2^{\ell}\left(1-\delta/\beta\right)^{-\ell}$.

### Mutational burden distribution

The single-cell mutational burden distribution differs from the division
distribution because there is an additional stochasticity at each cell
division due to mutational (Poisson) distributions. To obtain an expected MBD
from a DD, we can employ a procedure to introduce this stochasticity as
follows: each cell contributing to a division burden $D_{\ell,i}$ will have
undergone $\ell$ divisions, so will have acquired $\sum_{p=1}^{\ell}U_{p}$
mutations, where $U_{p}\sim\text{Pois}(\mu)$ represents the number of
mutations acquired during the cell’s $p$th division. Since the Poisson
distribution is additive, this sum is in turn a Poisson-distributed random
variable with mean $\ell\mu$. The left-hand side of Figure 4a qualitatively
depicts the elements of the DD being converted into Poisson probability mass
functions associated to these sums of Poisson variables. These are then summed
to obtain the MBD, shown on the right-hand side of Figure 4a, in the following
manner.

Writing $U_{p,q,\ell}\sim\text{Pois}(\mu)$ for the number of mutations
acquired during the $p$th division of the $q$th cell (for some labelling of
cells $1\leq q\leq D_{\ell,i}$) having undergone $\ell$ divisions, we can sum
over the elements of the DD labelled by $\ell$ to obtain

$\mathbb{E}\left[B_{k,i}\right]=\mathbb{E}\left[\mathbb{E}\left[B_{k,i}\,\middle|\,\\{D_{\ell^{\prime},i}\\}_{\ell^{\prime}\leq
i}\right]\right]=\mathbb{E}\left[\mathbb{E}\left[\sum_{\ell=0}^{i}\sum_{q=1}^{D_{\ell,i}}\mathbbm{1}_{\left\\{\sum_{p=1}^{\ell}U_{p,q,\ell}=k\right\\}}\,\middle|\,\\{D_{\ell^{\prime},i}\\}_{\ell^{\prime}\leq
i}\right]\right],$ (9)

where we have used the indicator function $\mathbbm{1}_{A}$ to be $1$ on the
set $A$ and $0$ elsewhere.

Figure 4: An illustration of the conversion of a division distribution (DD)
into a single-cell mutational burden distribution (MBD). (a) The elements of a
DD are translated into Poisson distributions with means $\ell\mu$, weighted by
$D_{\ell,i}$ (such that the sums over the teal and cyan distributions are
$D_{\ell^{\prime},i}$ and $D_{\ell^{\prime\prime},i}$, respectively, for
example), and then summed to obtain the corresponding MBD. Note that if the
mean of the DD is $\langle\ell\rangle$, then the mean of the resulting MBD is
$\langle\ell\rangle\mu$. (b) Average (solid dark pink line) of $200$
simulation realisations (representatives in solid pale pink lines) of the DD
for a pure-birth process up to final population size $N=10^{4}$ with
mutational mean $\mu=10$. (c) Average (solid dark blue line) of the MBD for
the same simulation realisations as (b) (representatives in solid pale blue
lines), along with the MBD obtained from converting the average DD as
explained in (a) and the main text (dashed red line).

Now, using the linearity of expectation, that $i$ is fixed and $\ell$ just an
index, and the independence of the random variables $D_{\ell,i}$ and
$U_{p,q,\ell}$, the right-hand side of equation (9) becomes

$\sum_{\ell=0}^{i}\mathbb{E}\left[D_{\ell,i}\mathbb{E}\left[\mathbbm{1}_{\left\\{\sum_{p=1}^{\ell}U_{p,q,\ell}=k\right\\}}\right]\right]=\sum_{\ell=0}^{i}\mathbb{E}\left[D_{\ell,i}\right]\mathbb{E}\left[\mathbbm{1}_{\left\\{\sum_{p=1}^{\ell}U_{p,q,\ell}=k\right\\}}\right].$

Finally, substituting the expression in equation (6) for the expected DD and
the probability mass function for the Poisson distribution with mean
$\ell\mu$, we find the pure-birth expected MBD:

$\mathbb{E}\left[B_{k,i}\right]=\sum_{\ell=0}^{i}\genfrac{[}{]}{0.0pt}{}{i}{\ell}\frac{2^{\ell}}{i!}\frac{\mathrm{e}^{-\ell\mu}(\ell\mu)^{k}}{k!}.$
(10)

Figures 4b–c verify the conversion from a DD to a MBD described in the
previous discussion with simulations.

The same conversion procedure can be implemented in the birth-death case.
Again, working with first-order approximations, the expression in equation (8)
for the expected birth-death DD can be used instead of the pure-birth
expression in equation (6) during the final step to obtain a first-order
approximation of the expected birth-death MBD.

Finally, consider the number of mutational occurrences: that is, the sum of
the entries of the mutational matrix $Y_{i}$ or, equivalently, either side of
equation (1). If this quantity is divided by the number of mutations $M_{i}$,
we obtain the mean of the SFS; and if it is divided by the population $N_{i}$,
we obtain the mean of the MBD. We can derive the expected number of mutational
occurrences using our recurrence relation approach, from which we deduce that
this mean, representing the expected mutational burden of a cell, grows
logarithmically with the step count $i$ (see Propositions S7 and S8). In the
pure-birth case, it is simply a rescaling of the harmonic numbers.

## Discussion

The distribution of genetic mutations in cell populations has been studied
both in the cases of constant [17, 35, 34, 25] and growing populations [20,
12, 36, 37, 38, 39, 40, 41]. With the advent of single-cell sequencing
technologies, exploration of more precise information in single cells is sure
to follow in the footsteps of population-level research [22, 23, 42]. At the
population level, both site frequency spectra (SFS) and overall tumour
mutational burden (TMB) have been investigated analytically [19, 20, 18, 12,
21, 25]. Here we focus on the single-cell distribution of the latter (the
single-cell mutational burden distribution, or MBD), and use the foundation of
the SFS to better understand the MBD analytically.

A new framework uniting the SFS and the MBD is presented, relying on a simple
procedure: dynamical matrices store the mutational information of a population
of cells, whose size is dictated by a birth-death process. Our approach of
encoding the data in binary matrices, where the entry $(n,m)$ is $1$ when cell
$n$ has mutation $m$ and $0$ otherwise, naturally emerges from the (neutral)
evolution-motivated idea wherein a cell is identified by its mutation load
[34]. Two different ways of partitioning the entries of this mutational matrix
provide definitions of both the SFS and the MBD as histograms of the row- and
column-sums, respectively, as shown in Figure 2. With this symmetry in mind,
which gives rise to equation (1), an identical analytical approach depending
on the discrete-time Markov nature of the model can be applied to both cases,
along with an intermediary case of the division distribution (DD), to obtain
recurrence relations for the distributions of interest: we employ the law of
total expectation to write the expected value of a quantity of interest in
terms of expected values at the previous time step. These recurrences are
solved exactly in the pure-birth case and approximately in the birth-death
case, giving rise to analytical predictions for the SFS, DD and MBD, which are
compared to stochastic simulations as well as previous work on the SFS and the
DD.

Indeed, our model recovers the expected values of the SFS and TMB (via the sum
of the SFS) derived by Gunnarsson et al. [21], as found in Propositions S7 and
S8 of the Supplementary Information. Our stochastic time first-order
approximation in equation (4) matches theirs from the stochastic population
scenario with a fixed elapsed time in the large-population limit, where the
regimes coincide according to their convergence analysis [21]. Our derivation
for the pure-birth DD in equation (6) recovers a result from previous work on
phylogenetic trees produced by Yule processes: combinatorics results relating
to binary search trees [30] were then applied to the phylogenetic context
[29].

Intuitively, we would think that the explicit single-cell MBD results from
both the DD and the extra stochasticity arising from the mutational
distribution at each past cell division (the internal nodes of the single-cell
phylogenetic tree). Surprisingly, we found that the the latter nodal
stochasticity does not play a large role in the MBD. While there is certainly
higher variance in the MBD than in the DD, as evidenced by Figures 4b–c, the
shapes of the two distributions remain similar and we can construct the MBD
based on the DD and $\mu$, the mean value of number of mutations acquired per
cell in each past cell division. The derivation from equation (9) to (10)
demonstrates that only the mean of the mutational distribution matters when
obtaining the MBD, rather than its higher moments. We further tested this
conclusion by applying other mutational distributions than the Poisson in
stochastic simulations, which lead to the same predicted MBD as shown in
Figure S2. Employing the binning procedure described in equation (S28) of the
Supplementary Information allows us to retrace our steps from the MBD to the
DD, which reinforces that it is only the mean of the mutational distribution
that is of critical importance to the shape of the MBD, not the exact form of
the distribution.

We showed that the expected mutational burden for an arbitrary cell in a
population (the mean of the MBD) increases logarithmically with the step count
$i$ in our model (see Propositions S7 and S8). In Moeller et al.’s continuous-
time framework, this mean is shown to be the product of the expected number of
divisions in the cell’s past and the mutational mean $\mu$ [25], much as we
have argued in Figure 4 for our conversion from the DD to the MBD. Under their
intuitive assumption of mutation burdens arising from a compound Poisson
distribution, the variance of the MBD is dependent on the means and the
variances of the DD and the mutational (Poisson) distribution [25], whereas
our derivations and simulations show that only the mean of the mutational
distribution plays a significant role, not its higher moments.

While single-cell sequencing is still in its adolescence, grappling with
hurdles such as trade-offs between sequencing noise, sample size and cost [43,
44], there is a growing need for mathematical and computational machinery to
handle the vast quantities of data being produced [22, 42]. Our model serves
as a new framework to integrate single-cell and bulk information, and shows
how various distributions of accumulated mutations are linked through the same
stochastic process.

### Acknowledgements

This research was supported by the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie EvoGamesPlus grant
number 955708. We thank Sabin Lessard, Nathaniel Mon Père and Alexander Stein
for fruitful discussions.

### Competing interests

The authors declare no competing interests.

## References

  * [1] Robert A Weinberg  “The Biology of Cancer”  Garland Science, 2013  DOI: 10.1201/9780203852569
  * [2] Thorsten B.H. Reusch, Iliana B. Baums and Benjamin Werner  “Evolution via somatic genetic variation in modular species” Publisher: Elsevier  In _Trends in Ecology & Evolution_ 36.12, 2021, pp. 1083–1092  DOI: 10.1016/j.tree.2021.08.011
  * [3] Taejeong Bae et al.  “Different mutational rates and mechanisms in human cells at pregastrulation and neurogenesis” Publisher: American Association for the Advancement of Science  In _Science_ 359.6375, 2018, pp. 550–555  DOI: 10.1126/science.aan8690
  * [4] Henry Lee-Six et al.  “Population dynamics of normal human blood inferred from somatic mutations”  In _Nature_ 561.7724, 2018, pp. 473–478  DOI: 10.1038/s41586-018-0497-0
  * [5] Benjamin Werner et al.  “Measuring single cell divisions in human tissues from multi-region sequencing data”  In _Nature Communications_ 11.1, 2020, pp. 1035  DOI: 10.1038/s41467-020-14844-6
  * [6] Steven A. Frank and Martin A. Nowak  “Problems of somatic mutation and cancer” Publisher: John Wiley & Sons, Ltd  In _BioEssays_ 26.3, 2004, pp. 291–299  DOI: 10.1002/bies.20000
  * [7] Natalia L. Komarova  “Cancer, aging and the optimal tissue design”  In _Somatic Evolution of Cancer Cells_ 15.6, 2005, pp. 494–505  DOI: 10.1016/j.semcancer.2005.07.003
  * [8] Rebecca A. Burrell, Nicholas McGranahan, Jiri Bartek and Charles Swanton  “The causes and consequences of genetic heterogeneity in cancer evolution”  In _Nature_ 501.7467, 2013, pp. 338–345  DOI: 10.1038/nature12625
  * [9] Ivana Bozic et al.  “Accumulation of driver and passenger mutations during tumor progression” Publisher: Proceedings of the National Academy of Sciences  In _Proceedings of the National Academy of Sciences_ 107.43, 2010, pp. 18545–18550  DOI: 10.1073/pnas.1010978107
  * [10] Shaoping Ling et al.  “Extremely high genetic diversity in a single tumor points to prevalence of non-Darwinian cell evolution” Publisher: Proceedings of the National Academy of Sciences  In _Proceedings of the National Academy of Sciences_ 112.47, 2015, pp. E6496–E6505  DOI: 10.1073/pnas.1519556112
  * [11] Andrea Sottoriva et al.  “A Big Bang model of human colorectal tumor growth”  In _Nature Genetics_ 47.3, 2015, pp. 209–216  DOI: 10.1038/ng.3214
  * [12] Marc J Williams et al.  “Identification of neutral tumor evolution across cancer types”  In _Nature Genetics_ 48.3, 2016, pp. 238–244  DOI: 10.1038/ng.3489
  * [13] Su-Chan Park and Joachim Krug  “Clonal interference in large populations” Publisher: Proceedings of the National Academy of Sciences  In _Proceedings of the National Academy of Sciences_ 104.46, 2007, pp. 18135–18140  DOI: 10.1073/pnas.0705778104
  * [14] Kasper Karlsson et al.  “Deterministic evolution and stringent selection during preneoplasia”  In _Nature_ 618.7964, 2023, pp. 383–393  DOI: 10.1038/s41586-023-06102-8
  * [15] W.. Ewens  “The pseudo-transient distribution and its uses in genetics” Edition: 2016/07/14 Publisher: Cambridge University Press  In _Journal of Applied Probability_ 1.1, 1964, pp. 141–156  DOI: 10.2307/3212065
  * [16] Motoo Kimura  “Genetic variability maintained in a finite population due to mutational production of neutral and nearly neutral isoalleles” Edition: 2009/04/14 Publisher: Cambridge University Press  In _Genetics Research_ 11.3, 1968, pp. 247–270  DOI: 10.1017/S0016672300011459
  * [17] Y.X. Fu  “Statistical Properties of Segregating Sites”  In _Theoretical Population Biology_ 48.2, 1995, pp. 172–197  DOI: 10.1006/tpbi.1995.1025
  * [18] Ivana Bozic, Jeffrey M. Gerold and Martin A. Nowak  “Quantifying Clonal and Subclonal Passenger Mutations in Cancer Evolution” Publisher: Public Library of Science  In _PLOS Computational Biology_ 12.2, 2016, pp. e1004731  DOI: 10.1371/journal.pcbi.1004731
  * [19] Richard Durrett  “Population Genetics of Neutral Mutations in Exponentially Growing Cancer Cell Populations” Publisher: Institute of Mathematical Statistics  In _The Annals of Applied Probability_ 23.1, 2013, pp. 230–250  DOI: 10.1214/11-aap824
  * [20] Richard Durrett  “Branching Process Models of Cancer”  In _Branching Process Models of Cancer_ Cham: Springer International Publishing, 2015, pp. 1–63  DOI: 10.1007/978-3-319-16065-8˙1
  * [21] Einar Bjarki Gunnarsson, Kevin Leder and Jasmine Foo  “Exact site frequency spectra of neutrally evolving tumors: A transition between power laws reveals a signature of cell viability”  In _Theoretical Population Biology_ 142, 2021, pp. 67–90  DOI: 10.1016/j.tpb.2021.09.004
  * [22] Ehud Shapiro, Tamir Biezuner and Sten Linnarsson  “Single-cell sequencing-based technologies will revolutionize whole-organism science”  In _Nature Reviews Genetics_ 14.9, 2013, pp. 618–630  DOI: 10.1038/nrg3542
  * [23] Yong Wang and Nicholas E. Navin  “Advances and Applications of Single-Cell Sequencing Technologies”  In _Molecular Cell_ 58.4, 2015, pp. 598–609  DOI: 10.1016/j.molcel.2015.05.005
  * [24] Federico Abascal et al.  “Somatic mutation landscapes at single-molecule resolution”  In _Nature_ 593.7859, 2021, pp. 405–410  DOI: 10.1038/s41586-021-03477-4
  * [25] Marius E. Moeller, Nathaniel V. Père, Benjamin Werner and Weini Huang  “Measures of genetic diversification in somatic tissues at bulk and single cell resolution suggest sources of unknown stochasticity”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2023  DOI: 10.1101/2022.11.07.515470
  * [26] Zachary R. Chalmers et al.  “Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden”  In _Genome Medicine_ 9.1, 2017, pp. 34  DOI: 10.1186/s13073-017-0424-2
  * [27] Evan M. Fernandez et al.  “Cancer-Specific Thresholds Adjust for Whole Exome Sequencing–Based Tumor Mutational Burden Distribution” Publisher: Wolters Kluwer  In _JCO Precision Oncology_ , 2019, pp. 1–12  DOI: 10.1200/PO.18.00400
  * [28] Elizabeth Martínez-Pérez, Miguel Angel Molina-Vila and Cristina Marino-Buslje  “Panels and models for accurate prediction of tumor mutation burden in tumor samples”  In _npj Precision Oncology_ 5.1, 2021, pp. 31  DOI: 10.1038/s41698-021-00169-0
  * [29] Mike Steel and Andy McKenzie  “Properties of phylogenetic trees generated by Yule-type speciation models”  In _Mathematical Biosciences_ 170.1, 2001, pp. 91–112  DOI: 10.1016/S0025-5564(00)00061-4
  * [30] William C. Lynch  “More Combinatorial Properties of Certain Trees”  In _The Computer Journal_ 7.4, 1965, pp. 299–302  DOI: 10.1093/comjnl/7.4.299
  * [31] Motoo Kimura  “The number of heterozygous nucleotide sites maintained in a finite population due to steady flux of mutations”  In _Genetics_ 61.4, 1969, pp. 893–903  DOI: 10.1093/genetics/61.4.893
  * [32] J Lederberg  “Replica plating and indirect selection of bacterial mutants: isolation of preadaptive mutants in bacteria by sib selection.”  In _Genetics_ 121.3, 1989, pp. 395–399  DOI: 10.1093/genetics/121.3.395
  * [33] Qi Zheng  “Mathematical Issues Arising From the Directed Mutation Controversy”  In _Genetics_ 164.1, 2003, pp. 373–379  DOI: 10.1093/genetics/164.1.373
  * [34] Roy Ronen, Nitin Udpa, Eran Halperin and Vineet Bafna  “Learning Natural Selection from the Site Frequency Spectrum”  In _Genetics_ 195.1, 2013, pp. 181–193  DOI: 10.1534/genetics.113.152587
  * [35] Richard Durrett  “Probability Models for DNA Sequence Evolution”, Probability and Its Applications  Springer New York, NY, 2008  DOI: 10.1007/978-0-387-78168-6
  * [36] Benjamin D. Simons  “Deep sequencing as a probe of normal stem cell fate and preneoplasia in human epidermis” Publisher: Proceedings of the National Academy of Sciences  In _Proceedings of the National Academy of Sciences_ 113.1, 2016, pp. 128–133  DOI: 10.1073/pnas.1516123113
  * [37] Lawrence A. Loeb et al.  “Extensive subclonal mutational diversity in human colorectal cancer and its significance” Publisher: Proceedings of the National Academy of Sciences  In _Proceedings of the National Academy of Sciences_ 116.52, 2019, pp. 26863–26872  DOI: 10.1073/pnas.1910301116
  * [38] Caroline J. Watson et al.  “The evolutionary dynamics and fitness landscape of clonal hematopoiesis” Publisher: American Association for the Advancement of Science  In _Science_ 367.6485, 2020, pp. 1449–1454  DOI: 10.1126/science.aay9333
  * [39] Gladys Y.. Poon, Caroline J. Watson, Daniel S. Fisher and Jamie R. Blundell  “Synonymous mutations reveal genome-wide levels of positive selection in healthy tissues”  In _Nature Genetics_ 53.11, 2021, pp. 1597–1605  DOI: 10.1038/s41588-021-00957-1
  * [40] Hwai-Ray Tung and Richard Durrett  “Signatures of neutral evolution in exponentially growing tumors: A theoretical perspective” Publisher: Public Library of Science  In _PLOS Computational Biology_ 17.2, 2021, pp. e1008701  DOI: 10.1371/journal.pcbi.1008701
  * [41] Monika K. Kurpas and Marek Kimmel  “Modes of Selection in Tumors as Reflected by Two Mathematical Models and Site Frequency Spectra”  In _Frontiers in Ecology and Evolution_ 10, 2022  DOI: 10.3389/fevo.2022.889438
  * [42] Heyrim Cho, Ya-Huei Kuo and Russell C. Rockne  “Comparison of cell state models derived from single-cell RNA sequencing data: graph versus multi-dimensional space”  In _Mathematical Biosciences and Engineering_ 19.8, 2022, pp. 8505–8536  DOI: 10.3934/mbe.2022395
  * [43] Samantha L. Goldman et al.  “The Impact of Heterogeneity on Single-Cell Sequencing”  In _Frontiers in Genetics_ 10, 2019  DOI: 10.3389/fgene.2019.00008
  * [44] Bora Lim, Yiyun Lin and Nicholas Navin  “Advancing Cancer Research and Medicine with Single-Cell Genomics”  In _Cancer Cell_ 37.4, 2020, pp. 456–470  DOI: 10.1016/j.ccell.2020.03.008
  * [45] William Feller  “An Introduction to Probability Theory and Its Applications”  John Wiley & Sons, 1968
  * [46] Maria Deijfen and Mathias Lindholm  “Growing networks with preferential deletion and addition of edges”  In _Physica A: Statistical Mechanics and its Applications_ 388.19, 2009, pp. 4297–4303  DOI: 10.1016/j.physa.2009.06.032
  * [47] Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik  “Table of Integrals, Series, and Products”  Academic Press, 2007

## Supplementary Information

We first derive results about the expected population size described by the
Markov chain depicted in Figure 1a and its probability of extinction, which
allow us to make rigorous our main text statement that all of our expectations
are conditioned on survival. We then prove several results from the main text.
Next, we make additional statements about the mutational occurrences $C_{i}$
defined by (1) and the number of unique mutations $M_{i}$ using our recurrence
relation approach. Finally, we comment on what happens if other (than Poisson)
mutational distributions are used, and we define a conversion from a single-
cell mutational burden distribution (MBD) to a division distribution (DD).
Table S1 summarises the new notation for Supplementary Information.

Symbol | Description
---|---
$\gamma_{z}$ | Extinction probability for a birth-death process starting at $N_{0}=z$
$u_{z,n}(p,q;a)$ | Random walk absorption probability; see Proposition S3 for full definition
$T$ | Markov chain transition matrix, with entries $T_{m,n}=\mathbb{P}\left[N_{i}=n\,\middle|\,N_{i-1}=m\right]$
$\zeta_{k,i}$ | coefficient equal to $1-(i-k)/(k+1)$ for $\lfloor(i+1)/2\rfloor\leq k\leq i$
$C_{i}$ | Number of mutational occurrences at step $i$: the sum of the entries of $Y_{i}$
Table S1: Notation used in Supplementary Information, in addition to those in
Table 1.

### Extinction probabilities and expected population sizes

Consider the Markov chain on a finite state space depicted by Figure 1a.

######  Proposition S1.

The extinction probability for the birth-death process starting at $N_{0}=z$
with $\beta\neq\frac{1}{2}$ is

$\mathbb{P}\left[N_{i}=0\text{ for some }i\geq
0\,\middle|\,N_{0}=z\right]=\gamma_{z}=\frac{\left(\delta/\beta\right)^{N}-\left(\delta/\beta\right)^{z}}{\left(\delta/\beta\right)^{N}-1}.$
(S1)

###### Proof.

Conditioning on the previous event (birth or death), we have for $1\leq z\leq
N-1$:

$\gamma_{z}=\beta\gamma_{z+1}+\delta\gamma_{z-1},$ (S2)

with boundary conditions $\gamma_{0}=1$ and $\gamma_{N}=0$. We note that (S2)
is a linear homogeneous second-order recurrence relation, so if we find two
particular solutions (in this case $1$ and $\left(\delta/\beta\right)^{z}$),
the general solution is a linear combination of them. The coefficients $A$ and
$B$ of this linear combination are found via the boundary conditions:

$A+B=1\quad\text{ and }\quad A+B\left(\delta/\beta\right)^{N}=0,$

which solve to give (S1). ∎

######  Remark.

By symmetry, we note that the probability to absorb into the $N$ state (having
started at $N_{0}=z$) is equal to $1-\gamma_{z}$, as can be seen by making the
change of variables

$\\{\beta,\delta,z\\}\to\\{\delta,\beta,N-z\\},$ (S3)

which allows us to view the Markov chain in reverse, in a sense. Therefore we
conclude that only on a set of zero probability does this Markov chain not
reach an absorbing state. Sometimes, we will consider the large-population
limit $N\to\infty$, where then the extinction probability is simply
$(\delta/\beta)^{z}$. For large enough $i$, we then obtain the approximation
for the survival probability
$\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]\simeq 1-\delta/\beta$.

######  Remark.

A second interpretation of this Markov chain is of a biased random walk on the
integers $\\{0,\dots,N\\}$. This specific set-up with absorbing boundaries is
known as the gambler’s ruin problem; the change of variables in (S3) is
analogous to putting ourselves in the other gambler’s shoes.

We can now describe the expected value of this Markov chain at step $i$. We
assume that $N$ is large enough; that is, that the maximum population size
$i+1$ (having started at $N_{0}=1$) is less than $N$.

######  Proposition S2.

The expected population at step $i$ of the birth-death process in Figure 1a is

$\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1\right]=\sum_{\begin{subarray}{c}n=1\\\
i-n\text{
odd}\end{subarray}}^{i+1}\frac{n^{2}}{i+1}\binom{i+1}{\frac{i-(n-1)}{2}}\beta^{\frac{i+(n-1)}{2}}\delta^{\frac{i-(n-1)}{2}},$
(S4)

and when conditioned on survival, we have
$\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]=\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1\right]/\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]$.

###### Proof.

Consider the transition matrix $T$ (independent of $i\geq 1$ by the
homogeneity of our Markov chain):

$T=\begin{pmatrix}1&0&0&\cdots&0&0\\\ \delta&0&\beta&\cdots&0&0\\\
0&\delta&0&\cdots&0&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\cdots&0&\beta\\\ 0&0&0&\cdots&0&1\end{pmatrix}\quad\text{ for entries
}\quad T_{m,n}=\mathbb{P}\left[N_{i}=n\,\middle|\,N_{i-1}=m\right].$

A property of stochastic Markov matrices is that the entries of the powers of
the matrix represent the multi-step probabilities:
$T^{i}_{m,n}=\mathbb{P}\left[N_{i}=n\,\middle|\,N_{0}=m\right]$ (this can be
simply proven, or concluded from the Chapman-Kolmogorov equations). We are
interested in the first row of these matrices $T^{i}$: that is, the case
$m=1$.

Note that the survival probability
$\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]$ is thus

$\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]=\sum_{n=1}^{i+1}T^{i}_{1,n}=1-T^{i}_{1,0}.$
(S5)

This sum as $i\to\infty$ is $1-\delta/\beta$, as noted in the remark following
Proposition S1.

For $i\geq 1$ the entries of the matrix powers satisfy the following
recurrence relation

$T^{i+1}_{1,n}=T_{n-1,n}T^{i}_{1,n-1}+T_{n+1,n}T^{i}_{1,n+1}=\beta
T^{i}_{1,n-1}+\delta T^{i}_{1,n+1}\quad\text{ for }\quad 1<n<N-1,$ (S6)

with boundary conditions $T^{i+1}_{1,1}=\delta T^{i}_{1,2}$,
$T^{i+1}_{1,N-1}=\beta T^{i}_{1,N-2}$, $T^{i+1}_{1,0}=T^{i}_{1,0}+\delta
T^{i}_{1,1}$ and $T^{i+1}_{1,N}=T^{i}_{1,N}$.

For $n>0$, consider the transition probability $T^{i}_{1,n}$: undergoing
births or deaths, $i$ steps are taken to progress from state $1$ to state $n$.
So, the sum of powers of $\beta$ and $\delta$ should be $i$ and their
difference should be $n-1$. The integer coefficients turn out to form the
Catalan triangle (see the proof of Corollary S2 for an interpretation of this
fact), so we obtain the ansatz

$T^{i}_{1,n}=\frac{n}{i+1}\binom{i+1}{\frac{i-(n-1)}{2}}\beta^{\frac{i+(n-1)}{2}}\delta^{\frac{i-(n-1)}{2}}\quad\text{
for }\quad 1\leq n\leq i+1,$ (S7)

where the binomial coefficients are defined to be zero if $i$ and $n$ have the
same parity, since it is impossible to reach an odd (respectively even)
destination $n$ in an even (respectively odd) number of steps having started
from $1$. It is straightforward to show that the expression (S7) satisfies the
recurrence (S6):

$\displaystyle\beta T^{i}_{1,n-1}+\delta T^{i}_{1,n+1}$
$\displaystyle=\beta\frac{n-1}{i+1}\binom{i+1}{\frac{i-(n-2)}{2}}\beta^{\frac{i+(n-2)}{2}}\delta^{\frac{i-(n-2)}{2}}+\delta\frac{n+1}{i+1}\binom{i+1}{\frac{i-n}{2}}\beta^{\frac{i+n}{2}}\delta^{\frac{i-n}{2}}$
$\displaystyle=\frac{i!\beta^{\frac{i+n}{2}}\delta^{\frac{i-n}{2}+1}}{\left(\frac{i-n}{2}\right)!\left(\frac{i+n}{2}\right)!}\left(\frac{n-1}{\frac{i-n}{2}+1}+\frac{n+1}{\frac{i+n}{2}+1}\right)$
$\displaystyle=\frac{n(i+1)!}{\left(\frac{i-n}{2}+1\right)!\left(\frac{i+n}{2}+1\right)!}\beta^{\frac{i+n}{2}}\delta^{\frac{i-n}{2}+1}$
$\displaystyle=T^{i+1}_{1,n}.$

Using (S7), the expected population size non-conditioned on survival is then

$\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1\right]=\sum_{n=0}^{i+1}n\mathbb{P}\left[N_{i}=n\,\middle|\,N_{0}=1\right]=\sum_{\begin{subarray}{c}n=1\\\
i-n\text{
odd}\end{subarray}}^{i+1}\frac{n^{2}}{i+1}\binom{i+1}{\frac{i-(n-1)}{2}}\beta^{\frac{i+(n-1)}{2}}\delta^{\frac{i-(n-1)}{2}},$

as desired. Conditioning on survival, we find

$\displaystyle\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]$
$\displaystyle=\sum_{n=0}^{i+1}n\mathbb{P}\left[N_{i}=n\,\middle|\,N_{0}=1,N_{i}>0\right]$
$\displaystyle=\sum_{n=0}^{i+1}n\frac{\mathbb{P}\left[N_{i}=n\cap
N_{i}>0\,\middle|\,N_{0}=1\right]}{\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]}$
$\displaystyle=\sum_{n=1}^{i+1}n\frac{\mathbb{P}\left[N_{i}=n\,\middle|\,N_{0}=1\right]}{\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]}$
$\displaystyle=\frac{\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1\right]}{\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]}.\qed$

######  Corollary S1.

For $k$ the number of births out of $i\leq 2k$ trials, we have

$T^{i}_{1,k-(i-k)+1}=\zeta_{k,i}\mathbb{P}\left[X=k\right],$ (S8)

where $X\sim\text{Binom}(i,\beta)$ is a binomially-distributed random variable
with probability $\beta$ of success and $\zeta_{k,i}$ is defined for
$\lfloor(i+1)/2\rfloor\leq k\leq i$ by

$\zeta_{k,i}=1-\frac{i-k}{k+1}\in(0,1].$ (S9)

###### Proof.

Write $k$ for the number of births out of the $i$ trials, so that
$n=k-(i-k)+1$. Then, (S7) becomes

$T^{i}_{1,n}=\frac{2n}{i+n+1}\mathbb{P}\left[X=\frac{i+(n-1)}{2}\right]=\left(1-\frac{i-k}{k+1}\right)\mathbb{P}\left[X=k\right].\qed$

######  Remark.

By itself, $\mathbb{P}\left[X=k\right]$ gives the probability of $k$ successes
out of $i$ trials, where the successes can occur anywhere. The inclusion of a
factor of $\zeta_{k,i}$ in (S8) thus describes the case where the number of
successes (births) must outweigh the number of failures (deaths) after any of
the first $1\leq i^{\prime}\leq i$ trials had the experiment ended there,
since we require the population to be non-extinct at each step. For example,
we want the ordered sequence $bbdd$ of births ($b$) and deaths ($d$) to
contribute to $T^{4}_{1,1}$, but we don’t want to count the sequence $bddb$,
as the final $b$ in the latter isn’t sensible, since the population is extinct
after the third step. To formalise this, we state the following definitions
and ensuing corollary.

A $(p,q)$-birth-death sequence is a word made from $p$ letters $b$ and $q$
letters $d$, and a non-empty word obtained by removing (possibly zero) letters
from the end of the original word is called a truncation. For example, $bbdbd$
and $bbdb$ are both truncations of the $(3,2)$-birth-death sequence $bbdbd$. A
$(p,q)$-birth-death sequence is called surviving if none of its $p+q$
truncations contain strictly more $d$s than $b$s. We then have the following
interpretation of the coefficients $\zeta_{k,i}$.

######  Corollary S2.

The proportion of surviving $(k,i-k)$-birth-death sequences is $\zeta_{k,i}$.

###### Proof.

By observation, the binomial coefficient $\binom{i}{k}$ is the total number of
$(k,i-k)$-birth-death sequences. On the other hand, another interpretation of
the Catalan triangle coefficients $\zeta_{k,i}\binom{i}{k}$ arising in (S7) is
the number of $i$-step walks from $(0,1)$ to $(i,k-(i-k)+1)$, where each step
goes from $(x,y)$ to $(x+1,y\pm 1)$ and the walk stays in the positive
quadrant. Therefore, they count the number of surviving $(k,i-k)$-birth-death
sequences. ∎

######  Remark.

Indeed, the definition (S9) of $\zeta_{k,i}$ resembles the large-$i$ survival
probability $1-\delta/\beta$, since $i-k$ is the number of deaths and $k+1$ is
one more than the number of births. Thus the expression (S8) can be
interpreted as this ‘survival probability’ $\zeta_{k,i}$ multiplied by a
binomial distribution, which doesn’t see the absorbing boundary at $0$.

Now, we can use our new expression (S8) for the transition probabilities to
compare the exact expected population with an intuitive linear approximation.

######  Remark.

The sums in the expected population (S4) and the survival probability (S5)
only involve terms $1\leq n\leq i+1$ where $i-n$ is odd. When $n=k-(i-k)+1$,
this translates to summing over integers $\lfloor(i+1)/2\rfloor\leq k\leq i$.
Using (S5) and (S8), the expected population size conditional on survival is
then

$\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]=\frac{\sum_{k=\lfloor\frac{i+1}{2}\rfloor}^{i}(k-(i-k)+1)\zeta_{k,i}\mathbb{P}\left[X=k\right]}{\sum_{k^{\prime}=\lfloor\frac{i+1}{2}\rfloor}^{i}\zeta_{k^{\prime},i}\mathbb{P}\left[X=k^{\prime}\right]}.$
(S10)

######  Remark.

Consider a naïve derivation of the expected population size, conditioning on
step $i-1$

$N_{i}=\beta(N_{i-1}+1)+\delta(N_{i-1}-1)=N_{i-1}+\beta-\delta=\cdots=(\beta-\delta)i+1,$
(S11)

and taking expectations of both sides. This is an appropriate linear
approximation in the limit of low death, but doesn’t correctly condition on
survival, as shown by the discrepancy for low $\beta$ (see Figure S1).

Figure S1: Plot of the expected population size conditioned on survival
$\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]$ versus the birth
probability $\beta$ for different numbers of steps $i=100$ (pale grey),
$i=200$ (dark grey) and $i=400$ (black). The dotted lines show the exact
expression from Proposition S2, and the solid lines are the linear
approximation $(\beta-\delta)i+1$, showing excellent agreement for low death
$\delta\ll\beta$ regardless of $i$.

######  Remark.

A heuristic derivation of the linear approximation (S11) from the exact
expression (S10) is the following. For large $i$, the probability mass
function of the binomial distribution tends to a function with mean
$i\beta=k$, the expected number of births in $i$ steps, and variance
$i\beta(1-\beta)$. When computing the sums in (S10) for large $i$, consider
the dominant $k=i\beta$ term (where we gloss over whether or not $i\beta$ is
an integer) as the primary contribution to the sum:

$\frac{\sum_{k}(k-(i-k)+1)\zeta_{k,i}\mathbb{P}\left[X=k\right]}{\sum_{k^{\prime}}\zeta_{k^{\prime},i}\mathbb{P}\left[X=k^{\prime}\right]}\simeq\frac{\left(i\beta-(i-i\beta)+1\right)\zeta_{i\beta,i}}{\zeta_{i\beta,i}}=(2\beta-1)i+1,$

recovering (S11). It is also worth noting that
$\zeta_{i\beta,i}=1-(i-i\beta)/(i\beta+1)$ coincides with the large-$i$
survival probability $1-\delta/\beta$, supporting our intuition that
$\zeta_{k,i}$ can be thought of as a survival probability of sorts, as
anticipated by Corollary S2 and the remark that follows it.

The exact survival probability at step $i$ is given explicitly by (S5) and
(S7). A related result from Feller [45] gives the distribution of times until
the Markov process goes extinct.

######  Proposition S3 (Feller).

Consider a random walk on the integers $\\{0,\dots,a\\}$ (with absorption on
both end points) beginning at $z$ with probability $p\in(0,1)$ of moving right
and $q=1-p$ of moving left. The probability of absorption on the $0$ boundary
at step $n$ is

$u_{z,n}(p,q;a)=2^{n}p^{(n-z)/2}q^{(n+z)/2}\frac{1}{a}\sum_{k=1}^{a-1}\cos^{n-1}\frac{\uppi
k}{a}\sin\frac{\uppi k}{a}\sin\frac{\uppi zk}{a}.$ (S12)

###### Proof.

See Chapter XIV,4 of [45]. ∎

######  Remark.

Sending $a\to\infty$ in (S12) provides a simplification [45]:

$\lim_{a\to\infty}u_{z,n}(p,q;a)=\frac{z}{n}\binom{n}{\frac{n-z}{2}}p^{\frac{n-z}{2}}q^{\frac{n+z}{2}},$
(S13)

where the binomial coefficient implicitly assumes that $n$ and $z$ have the
same parity. This is the case of only a single absorbing boundary at $0$, as
in our considered birth-death case. The similarity between (S7) and (S13)
allows to state that $u_{1,i+1}=\delta T^{i}_{1,1}$, which is reasonable as
the former is the probability of going extinct on step $i+1$ and the latter is
the probability of being at $1$ on step $i$ multiplied by the death
probability $\delta$. We can then use (S13), the previous argument and
Corollary S1 to find alternate expressions for the survival probability, where
$X\sim\text{Binom}(i,\beta)$ as usual:

$\mathbb{P}\left[N_{i}>0\,\middle|\,N_{0}=1\right]=1-\sum_{\begin{subarray}{c}n=1\\\
n\text{
odd}\end{subarray}}^{i}u_{1,n}(\beta,\delta;N)=1-\delta\sum_{\begin{subarray}{c}m=0\\\
m\text{
even}\end{subarray}}^{i-1}T^{m}_{1,1}=1-\delta\sum_{\ell=0}^{\lfloor(i-1)/2\rfloor}\frac{\mathbb{P}\left[X=\ell\right]}{\ell+1}.$

### Proofs of recurrence relations results

Our approximate expected distributions depend on Taylor expansions of expected
values of functions of random variables, along with the expected population
size low-death approximation (S11).

######  Lemma S1.

Consider a function $f\in C^{3}$ of two random variables $A$ and $B$ with
finite expectations; the second-order Taylor expansion of
$\mathbb{E}\left[f(A,B)\right]$ about
$(\mathbb{E}\left[A\right],\mathbb{E}\left[B\right])$ is

$\mathbb{E}\left[f(A,B)\right]=f(\mathbb{E}\left[A\right],\mathbb{E}\left[B\right])+\frac{1}{2}f_{AA}\operatorname{var}(A)+f_{AB}\operatorname{cov}(A,B)+\frac{1}{2}f_{BB}\operatorname{var}(B)+O^{3},$
(S14)

where the second derivatives $f_{ij}=\partial_{i}\partial_{j}f$ are evaluated
at $(\mathbb{E}\left[A\right],\mathbb{E}\left[B\right])$ and $O^{3}$ is big-O
notation for any cubic terms.

######  Remark.

The first-order expansion of $\mathbb{E}\left[f(A,B)\right]$ about
$(\mathbb{E}\left[A\right],\mathbb{E}\left[B\right])$ is
$f(\mathbb{E}\left[A\right],\mathbb{E}\left[B\right])$.

######  Corollary S3.

For $\mathbb{E}\left[B\right]\neq 0$ and $f(A,B)=A/B$, to second order we have

$\mathbb{E}\left[\frac{A}{B}\right]\simeq\frac{\mathbb{E}\left[A\right]}{\mathbb{E}\left[B\right]}-\frac{\operatorname{cov}(A,B)}{\mathbb{E}\left[B\right]^{2}}+\frac{\mathbb{E}\left[A\right]\operatorname{var}(B)}{\mathbb{E}\left[B\right]^{3}}=\frac{\mathbb{E}\left[A\right]}{\mathbb{E}\left[B\right]}-\frac{\mathbb{E}\left[AB\right]}{\mathbb{E}\left[B\right]^{2}}+\frac{\mathbb{E}\left[A\right]\mathbb{E}\left[B^{2}\right]}{\mathbb{E}\left[B\right]^{3}}.$

We can now derive the birth-death site frequency spectrum, conditioned on
survival.

######  Proposition S4.

Starting from one mutation-free progenitor cell, the birth-death process with
new mutations arising with expectation $\mu$ has approximate expected site
frequency spectrum

$\mathbb{E}\left[S_{j,i}\,\middle|\,N_{0}=1,N_{i}>0\right]\simeq\sum_{j^{\prime}=0}^{\infty}\frac{2\mu(\delta/\beta)^{j^{\prime}}}{(j+j^{\prime})(j+j^{\prime}+1)}\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right].$
(S15)

###### Proof.

As in the main text, we use the law of total expectation to write a recurrence
relation for the SFS at step $i+1$ conditional on the SFS at step $i$:

$\displaystyle\mathbb{E}\left[S_{j,i+1}\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[S_{j,i+1}\,\middle|\,\\{S_{j^{\prime},i}\\}_{j^{\prime}\leq
i}\right]\right]$
$\displaystyle=\mathbb{E}\left[S_{j,i}+\beta\left(-\frac{jS_{j,i}}{N_{i}}+\frac{(j-1)S_{j-1,i}}{N_{i}}+2\mu\delta_{1,j}\right)+\delta\left(-\frac{jS_{j,i}}{N_{i}}+\frac{(j+1)S_{j+1,i}}{N_{i}}\right)\right],$

for $\delta_{\cdot,\cdot}$ the Kronecker delta symbol, and where we have
omitted the conditioning on survival and the initial condition $N_{0}=1$ for
brevity. The first term in each set of parentheses in the second line arises
because any $j$-abundant mutations in the dividing (respectively dying) cell,
the number of which is $jS_{j,i}/N_{i}$, will become $(j\pm 1)$-abundant and
thus no longer contribute to the $j$-site. The second terms depict the
corresponding contribution from any $(j\mp 1)$-abundant mutations in the
dividing (respectively dying) cell, which will become $j$-abundant. The
Kronecker delta source term represents the $2\mu$ new (thus $1$-abundant)
mutations arising in any division.

######  Remark.

For $N_{0}=1$, by observation, the expected population size and SFS at step
$i=1$ are $\mathbb{E}\left[N_{1}\right]=2\beta$ and
$\mathbb{E}\left[S_{1,1}\right]=2\mu\beta$, with the $\beta$ factors dropping
if we condition on survival.

In the pure-birth case, the missing telescoping calculation (for
$Q_{j}=\mathbb{E}\left[S_{j,i}\right]/(i+1)$) in the main text is

$Q_{j}=\frac{j-1}{j+1}Q_{j-1}=\frac{j-1}{j+1}\frac{j-2}{j}Q_{j-2}=\cdots=\frac{2}{j(j+1)}Q_{1}=\frac{2\mu}{j(j+1)},$

which is indeed (S15) when $\delta=0$ (as only the $j=0$ term in the infinite
sum remains). A proof similar to this was presented in [25].

Now, for the birth-death case: if next we consider only the first-order
approximation from (S14)

$\mathbb{E}\left[\frac{S_{j,i}}{N_{i}}\,\middle|\,N_{0}=1,N_{i}>0\right]\simeq\frac{\mathbb{E}\left[S_{j,i}\,\middle|\,N_{0}=1,N_{i}>0\right]}{\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]}=X_{j},$

where the final equality is the same ansatz as before (reasonable given our
knowledge from [21] that
$\mathbb{E}\left[S_{j,i}\,\middle|\,N_{0}=1,N_{i}>0\right]\propto\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]$).
Again, we have a second-order non-homogeneous linear recurrence relation for
$X_{j}$:

$X_{j}\mathbb{E}\left[N_{i+1}\,\middle|\,N_{0}=1,N_{i}>0\right]-X_{j}\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]=-jX_{j}+\beta(j-1)X_{j-1}+\delta(j+1)X_{j+1}+2\mu\beta\delta_{1,j}.$

For low death, where our first order approximation holds, from Figure S1 and
(S11) we note that the expected gain in population size in one time step is
$\beta-\delta$. Thus the previous expression becomes

$\delta(j+2)X_{j+2}-(\beta-\delta+j+1)X_{j+1}+\beta jX_{j}=0\quad\text{ for
}\quad j\geq 1,$ (S16)

where we have absorbed the source term into the boundary conditions $X_{1}=C$
and $X_{2}=\beta(C-\mu)/\delta$.

######  Remark.

One way to proceed is to use a shift operator $S$ (such that $SX_{j}=X_{j+1}$)
and factorise (S16), solving first for one root of the ensuing quadratic
equation in $S$ and then the other.

Instead, we use Laplace’s method [46] for solving a second-order recurrence
relation of the form

$\left(\alpha_{2}(j+2)+\beta_{2}\right)X_{j+2}+\left(\alpha_{1}(j+1)+\beta_{1}\right)X_{j+1}+\left(\alpha_{0}j+\beta_{0}\right)X_{j}=0,$
(S17)

which, in our case (S16), we have
$\\{\alpha_{i}\\}_{i=0}^{2}=\\{\delta,-1,\beta\\}$ and
$\\{\beta_{i}\\}_{i=0}^{2}=\\{0,\delta-\beta,0\\}$. To solve (S17) with
appropriate boundary conditions we must find an interval $[a,b]$ and a
function $h(t)$ such that

$\left[t^{j}h(t)\left(\alpha_{2}t^{2}+\alpha_{1}t+\alpha_{0}\right)\right]_{a}^{b}=0\quad\text{
and
}\quad\frac{h^{\prime}(t)}{h(t)}=\frac{\beta_{2}t^{2}+\beta_{1}t+\beta_{0}}{t(\alpha_{2}t^{2}+\alpha_{1}t+\alpha_{0})},$

where $h(t)$ is defined up to a multiplicative constant determined by the
boundary conditions. The solution to the recurrence relation (S17) is then

$X_{j}\propto\int_{a}^{b}t^{j-1}h(t)\mathrm{d}t.$ (S18)

Following [46], we find $[a,b]=[0,1]$ and

$h(t)=\frac{1-t}{\beta/\delta-t},$

which is valid for $\delta>0$. From a table of integrals [47] we have

$\int_{0}^{1}x^{\lambda-1}(1-x)^{\mu-1}(1-\xi
x)^{-\nu}\mathrm{d}x=B(\lambda,\mu)\,{}_{2}F_{1}(\nu,\lambda;\lambda+\mu;\xi)$
(S19)

for $\real(\lambda)>0$, $\real(\mu)>0$ and $|\xi|<1$, where we have defined
the beta function $B(w,y)$ as

$B(w,y)=\int_{0}^{1}s^{w-1}(1-s)^{y-1}\mathrm{d}s,$

and the hypergeometric function ${}_{2}F_{1}(k,l;m;z)$ as

${}_{2}F_{1}(k,l;m;z)=\sum_{n=0}^{\infty}\frac{(k)_{n}(l)_{n}}{(m)_{n}}\frac{z^{n}}{n!},$

for the rising Pochhammer symbol $(q)_{n}$

$(q)_{n}=\begin{cases}1&\text{ if }n=0,\\\ q(q+1)\cdots(q+n-1)&\text{ if
}n>0.\end{cases}$

Substituting $\\{\lambda,\mu,\nu,\xi\\}\to\\{j,2,1,\delta/\beta\\}$ into (S19)
and momentarily ignoring a constant factor of $\delta/\beta$ since $h(t)$
defined up to a constant, (S18) becomes

$\displaystyle X_{j}$
$\displaystyle\propto\int_{0}^{1}t^{j-1}\frac{1-t}{1-(\delta/\beta)t}\mathrm{d}t$
$\displaystyle=B(j,2)\,{}_{2}F_{1}(1,j;j+2;\delta/\beta)$
$\displaystyle=\left(\int_{0}^{1}t^{j-1}(1-t)\mathrm{d}t\right)\sum_{n=0}^{\infty}\frac{(1)_{n}(j)_{n}}{(j+2)_{n}}\frac{(\delta/\beta)^{n}}{n!}$
$\displaystyle=\left(\frac{1}{j}-\frac{1}{j+1}\right)\sum_{n=0}^{\infty}\frac{j(j+1)}{(j+n)(j+n+1)}\left(\frac{\delta}{\beta}\right)^{n},$

which is (S15) after cancelling the factor $j(j+1)$ and multiplying by $2\mu$
(so that $h(t)$ matches the boundary condition
$X_{1}=\mathbb{E}\left[S_{1,1}\right]/\mathbb{E}\left[N_{1}\right]=\mu$). ∎

######  Remark.

The integral (S18) is found via different methods in [21] for the continuous
time (stochastic population) case, though with bounds $[0,1-1/N]$. Thus in the
large population limit these coincide, which aligns with our intuition that in
the large population limit the continuous and discrete time cases should
converge, as discussed in [21].

We now derive the division distribution in the pure-birth and birth-death
cases.

######  Remark.

For the expected pure-birth DD, the main text is simply missing an observation
that the solution (6) solves the recurrence relation (5):

$\genfrac{[}{]}{0.0pt}{}{i+1}{\ell}\frac{2^{\ell}}{(i+1)!}=\left(i\genfrac{[}{]}{0.0pt}{}{i}{\ell}+\genfrac{[}{]}{0.0pt}{}{i}{\ell-1}\right)\frac{2^{\ell}}{(i+1)!}=\frac{i}{i+1}\genfrac{[}{]}{0.0pt}{}{i}{\ell}\frac{2^{\ell}}{i!}+\frac{2}{i+1}\genfrac{[}{]}{0.0pt}{}{i}{\ell-1}\frac{2^{\ell-1}}{i!}.$

######  Proposition S5.

For the birth-death process described in Figure 1a, the expected division
distribution is approximated to first order by

$\mathbb{E}\left[D_{\ell,i}\,\middle|\,N_{0}=1,N_{i}>0\right]\simeq\frac{\genfrac{[}{]}{0.0pt}{}{i}{\ell}2^{\ell}\left(1-\delta/\beta\right)^{-\ell}}{\sum_{\ell^{\prime}=1}^{i}\genfrac{[}{]}{0.0pt}{}{i}{\ell^{\prime}}2^{\ell^{\prime}}\left(1-\delta/\beta\right)^{-\ell^{\prime}}}\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right].$
(S20)

###### Proof.

We employ our same approach as for the pure-birth case:

$\displaystyle\mathbb{E}\left[D_{\ell,i+1}\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[D_{\ell,i+1}\,\middle|\,\\{D_{\ell^{\prime},i}\\}_{\ell^{\prime}\leq
i},N_{i+1}>0\right]\right]$
$\displaystyle=\mathbb{E}\left[D_{\ell,i}+\beta\left(-\frac{1}{N_{i}}D_{\ell,i}+\frac{2}{N_{i}}D_{\ell-1,i}\right)+\delta\left(-\frac{1}{N_{i}}D_{\ell,i}\right)\right]$
$\displaystyle\simeq\left(1-\frac{1}{\mathbb{E}\left[N_{i}\right]}\right)\mathbb{E}\left[D_{\ell,i}\right]+\frac{2\beta}{\mathbb{E}\left[N_{i}\right]}\mathbb{E}\left[D_{\ell-1,i}\right],$

up to first order, where in the second line the impact of a cell dying is
simply to reduce its contribution to the division distribution by its expected
amount. As before, we have omitted the conditioning on survival and the
initial condition $N_{0}=1$ for brevity.

Next, consider the ansatz

$\mathbb{E}\left[D_{\ell,i}\right]_{\text{az}}=\frac{\genfrac{[}{]}{0.0pt}{}{i}{\ell}(2\beta)^{\ell}(\beta-\delta)^{i-\ell}}{\prod_{i^{\prime}=1}^{i-1}\mathbb{E}\left[N_{i^{\prime}}\right]},$
(S21)

which satisfies the previous recurrence relation by noting that, as argued in
the remark following Proposition S2 (as well as by (S11) and Figure S1),
$\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]-1\simeq(\beta-\delta)i$
in the low-death limit:

$\displaystyle\mathbb{E}\left[D_{\ell,i+1}\right]_{\text{az}}$
$\displaystyle=\frac{\genfrac{[}{]}{0.0pt}{}{i+1}{\ell}(2\beta)^{\ell}(\beta-\delta)^{i+1-\ell}}{\prod_{i^{\prime}=1}^{i}\mathbb{E}\left[N_{i^{\prime}}\right]}$
$\displaystyle=\left(i\genfrac{[}{]}{0.0pt}{}{i}{\ell}+\genfrac{[}{]}{0.0pt}{}{i}{\ell-1}\right)\frac{(2\beta)^{\ell}(\beta-\delta)^{i+1-\ell}}{\mathbb{E}\left[N_{i}\right]\prod_{i^{\prime}=1}^{i-1}\mathbb{E}\left[N_{i^{\prime}}\right]}$
$\displaystyle=\frac{\mathbb{E}\left[N_{i}\right]-1}{\mathbb{E}\left[N_{i}\right]}\genfrac{[}{]}{0.0pt}{}{i}{\ell}\frac{(2\beta)^{\ell}(\beta-\delta)^{i-\ell}}{\prod_{i^{\prime}=1}^{i-1}\mathbb{E}\left[N_{i^{\prime}}\right]}+\frac{2\beta}{\mathbb{E}\left[N_{i}\right]}\genfrac{[}{]}{0.0pt}{}{i}{\ell-1}\frac{(2\beta)^{\ell-1}(\beta-\delta)^{i-(\ell-1)}}{\prod_{i^{\prime}=1}^{i-1}\mathbb{E}\left[N_{i^{\prime}}\right]}.$

It remains to show that the ansatz (S21) can be rewritten as (S20). In this
pursuit, we employ the following lemma.

######  Lemma S2.

For $x\neq 0$, the unsigned Stirling numbers defined by (7) satisfy

$\prod_{k=1}^{i}\left(\frac{k}{x}+1\right)=x^{-i}\sum_{k=1}^{i}\genfrac{[}{]}{0.0pt}{}{i}{k}(1+x)^{k}.$
(S22)

###### Proof of Lemma S2.

Note that both expressions coincide to $1+1/x$ when $i=1$; suppose they do for
$1\leq i\leq I$. For $i=I+1$, the right-hand side of (S22) is

$\displaystyle
x^{-(I+1)}\sum_{k=1}^{I+1}\genfrac{[}{]}{0.0pt}{}{I+1}{k}(1+x)^{k}$
$\displaystyle=x^{-(I+1)}\sum_{k=1}^{I+1}\left(I\genfrac{[}{]}{0.0pt}{}{I}{k}+\genfrac{[}{]}{0.0pt}{}{I}{k-1}\right)(1+x)^{k}$
$\displaystyle=x^{-(I+1)}\left(\sum_{k=1}^{I}I\genfrac{[}{]}{0.0pt}{}{I}{k}(1+x)^{k}+\sum_{k=1}^{I}\genfrac{[}{]}{0.0pt}{}{I}{k}(1+x)^{k+1}\right)$
$\displaystyle=\frac{1}{x}\left(I+1+x\right)x^{-I}\sum_{k=1}^{I}\genfrac{[}{]}{0.0pt}{}{I}{k}(1+x)^{k}$
$\displaystyle=\left(\frac{I+1}{x}+1\right)\prod_{k=1}^{I}\left(\frac{k}{x}+1\right),$

where we have used the definition (7) and the boundary conditions of the
unsigned Stirling numbers of the first kind, resulting the left-hand
expression of (S22), concluding the proof by induction. ∎

Applying Lemma S2 with $x=1/(\beta-\delta)$, (S21) transforms into

$\frac{\genfrac{[}{]}{0.0pt}{}{i}{\ell}(2\beta)^{\ell}(\beta-\delta)^{i-\ell}}{\prod_{i^{\prime}=1}^{i-1}\mathbb{E}\left[N_{i^{\prime}}\right]}=\frac{\genfrac{[}{]}{0.0pt}{}{i}{\ell}(2\beta)^{\ell}(\beta-\delta)^{i-\ell}}{(\beta-\delta)^{i}\sum_{\ell^{\prime}=1}^{i}\genfrac{[}{]}{0.0pt}{}{i}{\ell^{\prime}}(1+1/(\beta-\delta))^{\ell^{\prime}}}\mathbb{E}\left[N_{i}\right],$

which is the desired expression (S20) after the identifications

$\frac{2\beta}{\beta-\delta}=1+\frac{1}{\beta-\delta}=2\left(1-\frac{\delta}{\beta}\right)^{-1}.\qed$

######  Remark.

When the distributions of $D_{\ell,i}$ and $N_{i}$ are not sufficiently
clustered around their means (that is, when $\delta$ isn’t small), however,
the first-order approximation for the division distribution is not sound. As
both diverge as $i\to\infty$, we require their variances to grow more slowly
than their means do.

######  Proposition S6.

The variance of the population size $N_{i}$ for the birth-death process with
$N_{0}=1$, using the low-death approximation
$\mathbb{E}\left[N_{i}\,\middle|\,N_{0}=1,N_{i}>0\right]\simeq(\beta-\delta)i+1$,
satisfies

$\mathrm{var}\left(N_{i}\right)\simeq\left(1-(\beta-\delta)^{2}\right)i,$

and thus tends to $0$ as $\delta\to 0$ (since then $\beta\to 1$), as expected.

###### Proof.

We omit the conditioning on survival and the initial condition $N_{0}=1$ for
brevity. Using our usual approach, we have

$\mathbb{E}\left[N_{i+1}^{2}\right]=\mathbb{E}\left[\mathbb{E}\left[N_{i+1}^{2}\,\middle|\,N_{i}\right]\right]=\mathbb{E}\left[N_{i}^{2}\right]+2(\beta-\delta)\mathbb{E}\left[N_{i}\right]+1\simeq\mathbb{E}\left[N_{i}^{2}\right]+2(\beta-\delta)\left((\beta-\delta)i+1\right)+1.$

Solving the ensuing recurrence, we find

$\mathbb{E}\left[N_{i}^{2}\right]\simeq(\beta-\delta)^{2}i^{2}+\left(-(\beta-\delta)^{2}+2(\beta-\delta)+1\right)i+1,$
(S23)

from which we subtract
$\mathbb{E}\left[N_{i}\right]^{2}\simeq\left((\beta-\delta)i+1\right)^{2}$:

$\mathrm{var}\left(N_{i}\right)\simeq(\beta-\delta)^{2}i^{2}+\left(-(\beta-\delta)^{2}+2(\beta-\delta)+1\right)i+1-\left((\beta-\delta)i+1\right)^{2}=\left(1-(\beta-\delta)^{2}\right)i.\qed$

######  Remark.

There is an equivalent recurrence relation for the division distribution when
expanded to second order using Corollary S3, though it cannot be solved by
elementary methods:

$\displaystyle\mathbb{E}\left[D_{\ell,i+1}\right]$
$\displaystyle=\mathbb{E}\left[D_{\ell,i}\right]-\mathbb{E}\left[\frac{D_{\ell,i}}{N_{i}}\right]+2\beta\mathbb{E}\left[\frac{D_{\ell-1,i}}{N_{i}}\right]$
$\displaystyle\simeq\mathbb{E}\left[D_{\ell,i}\right]-\left(\frac{\mathbb{E}\left[D_{\ell,i}\right]}{\mathbb{E}\left[N_{i}\right]}-\frac{\mathbb{E}\left[D_{\ell,i}N_{i}\right]}{\mathbb{E}\left[N_{i}\right]^{2}}+\frac{\mathbb{E}\left[D_{\ell,i}\right]\mathbb{E}\left[N_{i}^{2}\right]}{\mathbb{E}\left[N_{i}\right]^{3}}\right)$
$\displaystyle\quad+2\beta\left(\frac{\mathbb{E}\left[D_{\ell-1,i}\right]}{\mathbb{E}\left[N_{i}\right]}-\frac{\mathbb{E}\left[D_{\ell-1,i}N_{i}\right]}{\mathbb{E}\left[N_{i}\right]^{2}}+\frac{\mathbb{E}\left[D_{\ell-1,i}\right]\mathbb{E}\left[N_{i}^{2}\right]}{\mathbb{E}\left[N_{i}\right]^{3}}\right)$
$\displaystyle=\mathbb{E}\left[D_{\ell,i}\right]-\frac{\mathbb{E}\left[D_{\ell,i}\right]}{\mathbb{E}\left[N_{i}\right]}+2\beta\frac{\mathbb{E}\left[D_{\ell-1,i}\right]}{\mathbb{E}\left[N_{i}\right]}+\frac{\mathbb{E}\left[N_{i}^{2}\right]}{\mathbb{E}\left[N_{i}\right]^{3}}\left(2\beta\mathbb{E}\left[D_{\ell-1,i}\right]-\mathbb{E}\left[D_{\ell,i}\right]\right)$
$\displaystyle\quad+\frac{1}{\mathbb{E}\left[N_{i}\right]^{2}}\sum_{\ell^{\prime}=1}^{i}\left(\mathbb{E}\left[D_{\ell,i}D_{\ell^{\prime},i}\right]-2\beta\mathbb{E}\left[D_{\ell-1,i}D_{\ell^{\prime},i}\right]\right),$

where we used the linearity of expectation to simplify and have omitted all
conditioning on survival and the initial condition $N_{0}=1$ for brevity, as
usual. The first three terms in the final expression are those from the first-
order expansion of Proposition S5, and the next term can be simplified using
(S23). Despite this, the equation above for
$\mathbb{E}\left[D_{\ell,i+1}\right]$ not only remains complicated but
includes terms of the form
$\mathbb{E}\left[D_{\ell,i}D_{\ell^{\prime},i}\right]$, which prevent it from
being solved using our usual recursive methods.

### Mutational occurrences

Recall the definition of the number of mutational occurrences from (1); we
will write $C_{i}$ for this quantity. Note that this is also equal to (up to
renormalisation by the number of unique mutations $M_{i}$ or the population
$N_{i}$) the mean of the site frequency spectrum $\\{S_{j,i}\\}_{j}$ or the
single-cell mutational burden distribution $\\{B_{k,i}\\}_{k}$, respectively.
In the pure-birth case, we have the following result.

######  Proposition S7.

In a pure-birth process with mutation rate $\mu$ starting from a single
mutation-free progenitor cell, the expected number of unique mutations and
mutational occurrences at step $i$ are given by

$\mathbb{E}\left[M_{i}\right]=2\mu i\quad\text{ and
}\quad\mathbb{E}\left[C_{i}\right]=2\mu(i+1)\left(H_{i+1}-1\right),$

for $H_{n}=\sum_{k=1}^{n}k^{-1}$ the $n$th harmonic number.

###### Proof.

With each birth, an expected $2\mu$ new mutations are added to the population,
so the result for $\mathbb{E}\left[M_{i}\right]$ is clear by observation.
Taking the sum over $j$ of the expected pure-birth SFS (3) as in the
definition $M_{i}=\sum_{j}S_{j,i}$, and using the linearity of expectation and
$j^{-1}(j+1)^{-1}=j^{-1}-(j+1)^{-1}$ to telescope terms, we also obtain the
desired expression $\mathbb{E}\left[M_{i}\right]=2\mu i$.

Using the usual law of total expectation approach, the recurrence relation for
$C_{i}$ is

$\mathbb{E}\left[C_{i+1}\right]=\mathbb{E}\left[\left(1+\frac{1}{N_{i}}\right)C_{i}+2\mu\right],$
(S24)

since not only are $2\mu$ new mutations arising with each division event, but
all of the mutations in the dividing cell are duplicated. The proportion of
cells with $k$ mutations is $B_{k,i}/N_{i}$, so the expected number of
mutations in the dividing cell is found by multiplying this proportion by $k$
and summing over $1\leq k\leq M_{i}$. This sum is exactly $C_{i}/N_{i}$ by
(1), giving rise to the parenthetical term in (S24).

Recalling that in the pure-birth case we have $N_{i}=i+1$, (S24) can be solved
recursively in a similar manner to the pure-birth case of Proposition S4: we
divide both sides of (S24) by $i+2$, make a change of variables to
$G_{i}=C_{i}/(i+1)$ such that $G_{1}=C_{1}/2=\mu$, telescope the $G_{i}$ terms
to find

$G_{i+1}=G_{1}+2\mu\sum_{k=1}^{i}\frac{1}{k+2}=\mu+2\mu\left(H_{i+2}-1-\frac{1}{2}\right)=2\mu\left(H_{i+2}-1\right),$

which, after multiplying by $i+2$ and relabelling indices, is the desired
expression. ∎

######  Proposition S8.

For a birth-death process with mutation rate $\mu$, the expected number of
unique mutations $\mathbb{E}\left[M_{i}\right]$ to first order has the
following form:

$\mathbb{E}\left[M_{i}\right]\simeq
2\mu\left(\beta-\frac{\beta(\beta-\delta)}{\delta}\log\left(1-\frac{\delta}{\beta}\right)(i-1)\right)$
(S25)

and the expected total mutational occurrences $\mathbb{E}\left[C_{i}\right]$
to first order obey

$\mathbb{E}\left[C_{i}\right]\simeq
2\mu\beta\mathbb{E}\left[N_{i}\right]\sum_{i^{\prime}=1}^{i}\frac{1}{\mathbb{E}\left[N_{i^{\prime}}\right]},$
(S26)

where all expectations are conditioned on survival and the initial condition
$N_{0}=1$.

###### Proof.

First we will find the expression (S25) for $\mathbb{E}\left[M_{i}\right]$.
Note that the previous contributions during birth events of the pure-birth
case are simply weighted by $\beta$; the death events mean that $1$-abundant
mutations found in the dying cell (whose expected number is
$\mathbb{E}\left[S_{1,i}/N_{i}\right]$, which we will expand to first order as
before) vanish from the count at step $i+1$. Thus, the recurrence relation for
$\mathbb{E}\left[M_{i}\right]$ is

$\mathbb{E}\left[M_{i+1}\right]=\mathbb{E}\left[M_{i}+2\mu\beta-\delta\frac{S_{1,i}}{N_{i}}\right],$

which is solved by summing over $i$ and telescoping to find

$\mathbb{E}\left[M_{i}\right]\simeq 2\mu\beta
i-\delta\sum_{i^{\prime}=1}^{i-1}\frac{\mathbb{E}\left[S_{1,i^{\prime}}\right]}{\mathbb{E}\left[N_{i^{\prime}}\right]}=2\mu\beta
i-\delta(i-1)\sum_{n=0}^{\infty}\frac{2\mu(\delta/\beta)^{n}}{(n+1)(n+2)},$
(S27)

where we have taken advantage of the fact that (under our approximations)
$\mathbb{E}\left[S_{1,i}\right]/\mathbb{E}\left[N_{i}\right]$ is independent
of $i$ and substituted (S15). A crude approximation would be this expression
for $i^{\prime}=1$:
$\mathbb{E}\left[S_{1,1}\right]/\mathbb{E}\left[N_{1}\right]=\mu$. However,
our approximations hold better for large $i$, so instead we compute

$\sum_{n=0}^{\infty}\frac{w^{n}}{(n+1)(n+2)}=1-\left(1-w\right)\sum_{m=0}^{\infty}\frac{w^{m}}{m+2}=1+\frac{1-w}{w^{2}}-\frac{1-w}{w^{2}}\sum_{m=1}^{\infty}\frac{w^{m}}{m}=\frac{1}{w}+\frac{1-w}{w^{2}}\log(1-w),$

which we substitute into (S27) with $w=\delta/\beta$ to obtain

$\mathbb{E}\left[M_{i}\right]\simeq 2\mu\beta
i-2\mu\delta(i-1)\left(\frac{\beta}{\delta}+\left(1-\frac{\delta}{\beta}\right)\frac{\beta^{2}}{\delta^{2}}\log\left(1-\frac{\delta}{\beta}\right)\right),$

which can be rearranged to find (S25).

To find the expression (S26) for the total mutational occurrences, note that
as before the birth contributions are weighted by $\beta$; the death
contributions to $C_{i+1}$ include the previous quantity $C_{i}$, minus the
number of mutations found in the dying cell, whose expected number is
$\mathbb{E}\left[C_{i}/N_{i}\right]$, as argued in the proof of Proposition
S7. Thus, the recurrence relation satisfied by $\mathbb{E}\left[C_{i}\right]$
is

$\mathbb{E}\left[C_{i+1}\right]=\mathbb{E}\left[\left(1+\frac{\beta-\delta}{N_{i}}\right)C_{i}+2\mu\beta\right].$

We notice that by (S11) the parenthetical term expands to first order to
$\mathbb{E}\left[N_{i+1}\right]/\mathbb{E}\left[N_{i}\right]$, so by dividing
by $\mathbb{E}\left[N_{i+1}\right]$ and solving for
$\mathbb{E}\left[C_{i}\right]/\mathbb{E}\left[N_{i}\right]$ (by the same
telescoping trick as in the proof of Proposition S7), we obtain

$\frac{\mathbb{E}\left[C_{i+1}\right]}{\mathbb{E}\left[N_{i+1}\right]}\simeq\frac{\mathbb{E}\left[C_{1}\right]}{\mathbb{E}\left[N_{1}\right]}+2\mu\beta\sum_{i^{\prime}=1}^{i}\frac{1}{\mathbb{E}\left[N_{i^{\prime}+1}\right]}=\mu\left(1+2\beta\left(\sum_{i^{\prime}=1}^{i+1}\frac{1}{\mathbb{E}\left[N_{i^{\prime}}\right]}-\frac{1}{2\beta}\right)\right)=2\mu\beta\sum_{i^{\prime}=1}^{i+1}\frac{1}{\mathbb{E}\left[N_{i^{\prime}}\right]},$

which gives the desired expression after multiplying both sides by
$\mathbb{E}\left[N_{i+1}\right]$ and relabelling indices. ∎

### Other mutational distributions

As noted in the main text, when other mutational distributions from the
Poisson distribution are used in simulations, the conversion from the division
distribution (DD) to the single-cell mutational burden distribution (MBD)
described by Figure 4 still holds under some conditions. In Figure S2, the
mutational distributions used are $\text{Unif}\\{1,\dots,2\mu-1\\}$ and
$\text{Geom}(1/\mu)$, both of which have mean $\mu$ and are recaptured by the
predicted distribution.

Figure S2: Conversions from the division distribution to the single-cell
mutational burden distribution as in Figure 4 of the main text, using other
(than Poisson) mutational distributions, with means $\mu$: (a) uniform
distribution over the set $\\{1,\dots,2\mu-1\\}$ and (b) geometric
distribution with mean $\mu$. Average (solid dark grey line) of $200$
simulation realisations (representatives in solid pale grey lines) and the
predicted MBD distribution converted from the DD (dashed red line).

If, however, the number of new mutations $U_{1}$ and $U_{2}$ (where the
indices refer to daughter cells) are drawn from a distribution with small
support, such as the uniform distribution $\text{Unif}\\{\mu-1,\mu+1\\}$, then
it is unsurprising that the smoothing procedure from a DD to a MBD of Figure 4
doesn’t hold, since fewer possible mutational burdens are obtainable for
cells. In the extreme case of the delta distribution $\text{Delta}(\mu)$, we
simply obtain a scaling of the DD, where all cells have mutational burden
equal to the product of the mutational mean and their division burden.

On the other hand, given an expected MBD, it is possible to recover the
expected DD via a binning procedure. By histogramming the MBD with bins of
width $\mu$ centred at integer multiples of $\mu$ (which results in the
distribution supposing that exactly $U_{1}=U_{2}=\mu$ mutations were acquired
with each division), we can then rescale by $\mu$ to obtain the DD.
Symbolically,

$\mathbb{E}\left[D_{\ell,i}\right]=\sum_{k\in\left[\left(\ell-\frac{1}{2}\right)\mu,\left(\ell+\frac{1}{2}\right)\mu\right)}\mathbb{E}\left[B_{k,i}\right].$
(S28)

In Figure 4a, this involves summing $\mu$ adjacent bins of the MBD, and then
rescaling the $x$-axis by $\mu$ to obtain the DD. This is exactly the
aforementioned distribution, obtained by using $\text{Delta}(\mu)$ as the
mutational distribution.