Introduction

Intra-tissue genetic heterogeneity emerges from a multitude of dynamical processes [1, 2]. Cells randomly accumulate mutations, they self-renew, differentiate or die and clones of varying fitness may compete and expand [3–6]. These core evolutionary principals form the basis for describing the aging of somatic tissues, tumour initiation and tumour progression [7–12]. Naturally, there is great interest in understanding precisely how these dynamical processes act. This typically involves hypothesizing models of cell behavior, and deriving quantitative estimates of their underlying physical parameters [3, 5, 13–17]. These could by quantities such as mutation rates per cell division, active stem cell numbers, symmetric and asymmetric division rates, etc.. Some such parameters are already well characterised in certain tissues. We know for example that during early development healthy somatic cells accumulate on average between 1.2 and 1.3 mutations per genome per division [5, 18]. Other parameters, such as the number of stem cells in a tissue type, are well known in some tissues but harder to quantify in others.

In general, a lack of time resolved information complicates evolutionary inferences, and we must often rely on indirect information, for example the patterns of tissue-specific genetic heterogeneity within and across individuals [19]. These observed patterns of genetic heterogeneity emerge from the underlying stochastic processes and can strongly depend on the measurement technique [20]. Both enforce specific limitations on the inferability of the underlying dynamics. Different types of genomic data – such as bulk and single cell whole genome sequencing – may contain different information on the system, and each comes with their own limitations. When bulk sequencing data is used for inference, some evolutionary parameters such as the population size and proliferation rates of stem cells are entangled and cannot be obtained in isolation [3, 16]. If on the other hand single cell data is employed, some additional information becomes available, for example the individual cell mutational burdens and the co-occurrence of mutations. However, often the sampling size is sparse compared to the presumed number of long-term proliferating cells in many tissues [18, 21, 22], which in turn introduces other constraints. Here we show that combining such data obtained from different resolutions (bulk or single cell) can help to overcome limitations and further narrow possible ranges of inferred parameters. Here, we develop concrete mathematical and computational models for extracting evolutionary parameters from bulk and single cell information and apply these methods to whole genome bulk sequencing data in healthy oesophagus [10] and whole genome single cell sequencing data in healthy haematopoiesis [18].

Results

A stochastic model of mutation accumulation in healthy somatic tissue

We model the stem cell dynamics of healthy tissue as a collection of individual cells that divide, differentiate, and die stochastically at predefined rates (see Fig.1). Novel mutations can occur with every cell division, each daughter acquiring a random number drawn from a Poisson distribution with rate μ. We explicitly include symmetric and asymmetric stem cell divisions at potentially different rates, the former resulting in two stem or two differentiated cells, and the latter in one stem and one differentiated cell. Since differentiated cells are lost from the stem cell pool they are functionally equivalent to cell death. While only symmetric divisions allow variants to change in frequency, both division types introduce novel variants into the population and thus both contribute to the evolving heterogeneity.

Figure 1.

Mutation accumulation occurs through all developmental stages from conception to death [23]. We therefore consider three phases of demographic changes of the stem cell population. In the first early developmental phase the stem cell population rapidly expands from a single cell exclusively by symmetric divisions at rate γ. This is followed by a growth and maintenance phase wherein the population continues to expand but also undergoes turnover through asymmetric divisions at rate ϕ, as well as cell removal (due to differentiation or death) and replacement at rate ρ. In the final mature phase, cell turnover continues but the size of the total stem population N remains constant.

We employ two independent implementations of our theoretical model. One is a direct stochastic simulation utilizing a Gillespie algorithm. This approach creates a simulated dataset containing all single cell or bulk sequencing properties, however it can be computationally expensive. In addition, we compute the time-dynamical expected value of the distribution of variant allele frequencies (the VAF spectrum) directly, which we find obeys the partial differential equation

where κ = f N(t) denotes the number of cells sharing a variant (the variant frequency f times the total population size N), δ(x) is the Dirac impulse function, ∂_t and ∂_κ are the partial derivatives with respect to time and variant size, and

This set of equations allows for computationally efficient numerical solutions (SI 1.3).

Transition of developmental growth and constant population size signatures in bulk whole genome sequencing data of healthy adult oesophagus

We first discuss established properties of the VAF spectrum relating to somatic tissues. In certain model systems it is given by a power law ∝ f ^−α with a critical exponent α. The value of the exponent depends on the demographic dynamics of the population. For a well-mixed exponentially growing population without cell death the VAF spectrum v(f) is given by 2μ/(f + f ²) (a f ⁻² power law) and is independent of time [6]. In contrast, for a population of constant size – i.e. where birth and death rates are equal – the spectrum obeys v(f) ∝ 2μ/ f [24] (a f ⁻¹ power law; see also SI 1.5.1), though this solution is only valid at sufficiently long times. In the following, we focus on the dynamics of genetic diversification in healthy tissues and show that the expected VAF spectrum contains both contributions of growing and constant population dynamics.

Healthy adult somatic tissues are thought to initially expand rapidly during ontogenic growth, continue to expand during infancy and childhood, and reach homeostasis by adulthood. It is thus natural to ask what the expected VAF spectrum would look like in tissues experiencing homeostasis after a period of growth. We first investigated the underlying VAF dynamics in a minimal theoretical model where the population switches from pure exponential growth to constant size once a maximal cell number is reached. Numerical solutions of Eq. (1) show that the expected VAF distribution exhibits a gradual transition from the f ⁻² (growing population) to the f ⁻¹ (constant population) power law (Fig. 2a). These transitional states themselves do not adhere to some intermediate power-law (e.g. f ^−α for 1 < α < 2), but instead present a sigmoidal shape, with the low frequency portion following f ⁻¹ and the high frequencies f ⁻². Over time the shape changes as a wavelike front traveling from low to high frequency, with the constant-size equilibrium establishing earliest at the lowest frequencies and moving to higher frequency over time. Interestingly, the convergence towards equilibrium slows down over time – for evenly-spaced observation times the solutions lie increasingly closer together – further decreasing the speed at which the high frequency portion of the spectrum approaches equilibrium. More complex forms of such demographic changes, for example the inclusion of a mixed growth-and- maintenance phase or a logistically growing population, result in qualitatively similar behavior (see Supplementary Fig. 5).

Figure 2.

To test whether adult human tissues show such transitional signatures, we grouped the VAFs of healthy adult oesophagus samples from Martincorena and colleagues [10] into three age groups (young, middle and old, Fig. 2b). In accordance with our theoretical expectations, the averaged VAF spectrum of the young donor group is closest to the expected f ⁻² distribution of ontogenic growth, with low frequency variants only starting to approach the f ⁻¹ homeostatic scaling. The averaged spectra of the middle and old age groups on the other hand notably shift towards the expected f ⁻¹ homeostatic line. Interestingly, while there is a clear separation between the spectra of the youngest and middle age groups, those of the middle and older groups are much closer. This agrees with our prediction that the speed of convergence of the spectra towards homeostatic scaling slows down with age.

VAF based mutation rate inferences in healthy tissues

The theoretical prediction for the VAF spectrum of an exponentially growing population, which is recovered in many cancers given sufficient sequencing depth [20], contains information on the effective mutation rate μ [25]. Similarly, the VAF spectrum of a constant population includes the same mutation rate μ (SI). Thus we can in principle infer μ from the previously described oesophagus data, for example by applying a regression approach. However, since the scaling of the VAF spectrum of healthy tissues can be shown to be age dependent (Fig. 2a-b), we use two reference shapes for a growing and a constant population, respectively. The fitting is performed by a linear regression least-squares approach. The general solution is computed from the VAF spectrum of the data, V_d(μ), and the reference shape V_r, which assumes μ = 1, as follows:

If V_d(μ) has the shape of V_r, then E(V_d(μ)) = μ E(V_r), which can be used to prove equation 2:

The VAF shapes for the growing and constant population give close estimates, as can be seen in fig. 2, which compares the age of all individuals with corresponding inferred effective mutation rates. The estimates are in the range of 1 to 2 mutations per genome per cell division and agree with previous observations based on early developmental stem cell divisions [18]. Estimates based on a growing population are slightly higher compared to estimates based on a constant population, but differences are small.

Surprisingly, we observe a clear trend of a linearly increasing effective mutation rate with age. This trend cannot be explained by the transition of variants from a growing into a constant population alone. An actual increase of the true mutation rate per cell division with age would be a direct explanation. However, this seems unlikely: It is well established that the mutational burden across individuals and tissues increases linearly with age [9, 11, 26], and an increasing mutation rate would in contrast result in an accelerated increase of mutational burden with age, which is not observed experimentally.

An alternative explanation is a continued slow linear increase of the stem cell population with age (Fig. 2d). In such a scenario, the effective mutation rate would appear to increase linearly with age as well, although the true mutation rate per cell division remains unchanged. At first glance, a linearly increasing stem cell population appears unnatural. However, such a linear increase would emerge from a small bias of stem cells towards self-renewal (growth) and concomitantly a slowdown of the stem cell proliferation rate proportional to the change of population size [27]. The latter is a natural feedback mechanism ensuring a constant output of differentiated somatic cells. Direct in vivo imaging in multiple human tissues including Oesaphagus seem to support an increased density and slower proliferation of stem cells with age [28]. This explanation is attractive for another reason: A feedback of decreased cell proliferation with increasing population size would also function as a tumour suppressor mechanism, and may be another reason cancer driver mutations in healthy tissues are abundant but rarely progress to cancer [10, 29]. Progression would require additional genomic or environmental changes to overcome this regulatory feedback.

Single-cell mutational burden allows sample size independent inferences of stem cell proliferation and mutation rates

With the increasing availability of single-cell data, it becomes possible to investigate distributions of certain quantities for which bulk sequencing only provides average measurements. One such is the number of mutations present in individual cells. Recent work by Lee-Six and colleagues [18] showed that in a 59 year old donor haematopoietic stem cells had accumulated on average around 1000 mutations per cell; however there was significant variation between individual genomes, some carrying 900 mutations, while others up to 1200. This variation is the result of the stochastic nature of mutation accumulation and cell division, and can be exploited to estimate per division mutation and stem cell proliferation rates. To this end we look at the distribution of the number of mutations per cell, which we will from here on refer to as the mutational burden distribution. Formally, we describe the mutational burden m_j in a cell j by a stochastic sum over its past divisions , in which both the total number of divisions y_j in the cell’s past and the number of mutations per division u_ij are random variables. Even without knowing the actual distributions of y and u, general expressions for the expectation and variance of m_j that only depend on the expectation and variance of y and u exist and are given by E(m) = E(y) E(u) and Var(m) = E(y) Var(u) + E(u)² Var(y) (SI 1.7). These expressions allow direct estimates of the effective division and mutation rate. Assuming the mutation rate per daughter cell is Poisson distributed with expectation μ, the average number of divisions per lineage — expressed in terms of an effective division rate λ — becomes

In addition, if the mutation rate μ is unknown, it can similarly be obtained from

We first show that these estimators recover per division mutation and proliferation rates in stochastic simulations (Fig. 3a and b). More importantly, they only require information on comparably few single cells. Approximately, 100 cells are sufficient to reliably reconstruct the mutational burden distribution, which for example would constitute a small sample of the haematopoietic stem cell pool. In fact, the inference does not significantly improve if we increase data resolution beyond a few hundred single cells (Fig. 3a and b). This is in stark contrast to inferences based on single cell phylogeny: Although they contain more information in principle, important aspects of a single cell phylogeny, e.g. the number of observable branchings, are strongly affected by sampling and stem cell population size.

Figure 3.

We apply our estimators to single-cell mutational data in blood obtained from the study by Lee-Six et al. [18], in which 89 haematopoietic stem cells (HSCs) were extracted from a bone marrow aspirate of a 59-year-old individual (Fig. 4a). Using Eq. (4) to estimate the expected number of divisions per lineage from the sample burden (SI 1.8.1), we find a total division rate (including both symmetric and asymmetric divisions) in the constant population phase of λ = 10.6 (7.6-15.3) divisions per HSC per year, which is within ranges previously suggested [3, 30]. Notably, the single-cell mutational burden alone can not disentangle symmetric and asymmetric division rates. Applying Eq. (5) to the sample burden distribution and taking Var(y)/ E(y) = 1 (akin to assuming exponentially distributed division times, a common but simplified null model) we obtain an estimated mutation rate of per cell per division (Fig. 4a). This is significantly higher compared to the 1.2 mutations per division suggested previously [5, 18]. We posit different possible explanations for this discrepancy. Stem cell divisions and mutation accumulation are stochastic processes, therefore variation between individuals is expected. However, simulations of our model for a wide range of parameter values suggest the distribution of these inferred mutation rates is unlikely to present a wide enough variance to explain this difference (Fig. 7). Recently, Mitchell and colleagues [12] suggested that many mutations in somatic tissue are acquired independently of cell divisions. In fact, the previous estimates of 1.2 to 1.3 mutations per genome per division were derived from early developmental cell divisions where effects of aging would be negligible. The increased mutation rate inferred here could in part be explained by proliferation independent effects. However, it is important to note that the true distribution of divisions per lineage may also be over-dispersed compared to the Poisson model (i.e. Var(y) > E(y)), which would lead to an overestimation of μ in Eq. (5). Furthermore, the observed burden distribution is incompatible with a simple division-independent mutation model, which would lead to a Poisson distribution of the burdens that has a variance much smaller than what we observed (Fig. 4a). A combination of division independent mutations together with non-Poissonian cell divisions could possibly reproduce the result.

Figure 4.

Figure 5.

Figure 6.

Figure 7.

Figure 8.

Sparse sampling, single cell derived VAF spectra and evolutionary inferences

Current methods limit us to whole genome information of at most a few thousand single cells per tissue or tumour [22]. In many situations this will only constitute a small fraction of the underlying long term proliferating cell population. For example, although the exact number of HSCs remains unknown, most recent estimates suggest a population size of 10⁵ cells [3, 18]. In tumours, this number can be as large as 10¹⁰ cells [31]. As we have shown in the previous paragraph, the single-cell mutational burden distribution allows certain inferences even if sampling is sparse. If we however construct the VAF spectrum from such a small sample, the distribution of variants will be significantly transformed with respect to that of the total population (Fig. 3c). Since Eq. (1) gives the expected curve measured from the total population, one must include a correction, which is given by the expectation of hypergeometric sampling (i.e. without replacement). This can be obtained through the transformation (see Supplement 1.6)

with the VAFs observed in the sample and S the sample size. Furthermore, from stochastic simulations we note that the variance on the distribution increases with variant frequency f, making the lowest frequency state (i.e. 1/S) the best candidate for comparative fitting and evolutionary inferences (see Supplement Fig. 9).

Figure 9.

Taking the mutation rate from Werner et al. [5] and the total division rate obtained from the burden distribution as constant, we explored solutions of Eq. (1) (sampled by Eq. (6)) for a wide range of realistic values in the remaining parameter space {N, p, t_M, N_H }. Comparing the lowest sampled frequency state (f = 1/89) with the data (Fig. 4b) we identify a curve in the space of the numbers of HSCs N and p (fraction of divisions that are asymmetric) where theoretically predicted VAF distributions are identical (Fig. 4c). Interestingly, this curve in N-p space naturally identifies a maximal stem cell population size capable of producing the data, which occurs when stem cells divide exclusively symmetrically (p = 0). Variations of t_M (the time until maturity of the HSC population) and N_H (the HSC pool size at which maintenance begins to occur) had only a small effect on the this inference (Fig. 4d). This is partly because mutations arising during these phases are more likely to be found at higher frequencies upon measurement [13, 32]. We find the data to be congruent with an adult HSC pool of at most 200’000-300’000 cells, depending on the exact timings of the maturation phase (Fig. 4d). This estimate agrees with the original study and other independent inferences based on population data [3, 18]. However, this upper bound corresponds to the case of exclusively symmetric stem cell divisions. Accounting for the possibility of mutating asymmetric divisions reduces the estimated stem cell pool. For example, the data would be consistent with 50’000 stem cells if 90% (p = 0.9) of all stem cell divisions are asymmetric. In principle, even smaller stem cell pool sizes are possible, though decreasing an order of magnitude would imply an extremely large portion (> 99%) of asymmetric HSC divisions. We note that in the data non-developmental branchings are observed, which immediately implies that not all HSC divisions can be asymmetric (p < 1) and thus a scenario of an extremely small population of only a handful of HSCs maintaining heamatopoisis can be rejected. With an orthogonal population based method, Watson and colleagues estimated 25’000 as a reasonable lower bound for the number of HSCs. Based on our inference, this would imply 0.95 as an upper bound for the proportion of all stem cell divisions that are asymmetric.

Discussion

Here, we have shown that single cell and bulk whole genome sequencing can inform on different aspects of somatic evolution. Single cell information does support previous bulk based estimates of possible ranges of HSC numbers. Surprisingly, the same single cell data leads to a mutation rate inference during homeostasis approximately 4 times higher compared to previous developmental estimates.

Another open question is the role of selection and how it shapes intra-tissue genetic heterogeneity. Evidence is emerging that positively selected variants in blood are almost universally present in individuals above 60, while the effective observable dynamics in younger individuals is well described by neutral dynamics. How results presented here generalize or modify will critically depend on the model of selection realized in human hematopoiesis, e.g. a models of rare or frequent driver events. Details of the underlying biology are currently unknown.

An important recent study has suggested an amended model of mutation accumulation in somatic tissues, wherein the majority of mutations are acquired continuously over time throughout a cell’s lifespan rather than during DNA replication [12]. We note that for a time dependent rate of mutation accumulation that is constant across all cells at a given point in time, this model predicts Poisson distributed single-cell mutational burdens. In contrast, we find the single-cell mutational burden distribution in human HSCs to be highly over-dispersed compared to a Poisson model. A time-associated model of mutation accumulation with a constant rate alone cannot generate such a pattern. At least one additional source of stochasticity would be needed. This could be fluctuations in the time dependent mutation rates, non-constant cell proliferation rates, e.g. bursts of stem cell divisions potentially in reaction to injury or disease, or a combination of both effects. Either way, our observations suggest that our current theoretical models of somatic evolution in healthy tissues are incomplete. There is evidence for unknown processes contributing to intra-tissue genetic heterogeneity. The precise nature of these processes remains an open question.

Acknowledgements

The authors are grateful to Marc Williams and Christo Morison for discussions at different stages of the manuscript.

N.V.MP. acknowledges support by Télévie and the Barts Charity. B.W. is supported by a Barts Charity Lectureship (grant no. MGU045) and a UKRI Future Leaders Fellowship (grant no. MR/V02342X/1). This project has been facilitated by an Alan Turing Institute Data Science for Science Program grant TU/ASG/R-SPES-121.

1 Supplementary Information

1.1 Model system

We model a dynamical population of N(t) cells, each able to carry multiple mutations. It evolves by the occurrence of divisions — which replace cells in the pool with their (potentially mutated) daughter cells — and removals — which can correspond to differentiation or any form of cell inactivation. As outlined in Figure 1, we model different types of divisions affecting the system:

1. A symmetric (non-differentiating) division replaces a cell in the population with two daughter cells, each of which maintains the mutational profile of its parent and acquires a Poisson distributed number of new mutations.

2. A symmetric differentiation removes a cell from the population, as both of its daughter cells have differentiated.

3. an asymmetric division replaces a cell with a single daughter cell — as the other daughter has differentiated — and also acquires a Poisson distributed number of mutations.

The rates at which these divisions occur depends on what stage of development the population is in. A simplest approach is to consider two broad phases of evolution: an initial growth phase where N(t) increases according to some growth function, and a mature phase at a time t_M where N(t) is constant and stem cell divisions occur only to maintain homeostasis of the population. While the earliest stage of population expansion can reasonably be approximated as purely developmental — i.e. where little to no turnover occurs in the population — maintenance of the functional population (e.g. in the blood system) will at some point be required even as the stem cell pool continues grow. Thus we also include a mixed phase occurring before maturity but only after the stem cell pool has reached a size N_H, where stem cells can differentiate and homeostasis must be maintained. Note that this size is expected to be relatively small, as the earliest haematopoietic cells appear 16 days post conception [33], and peripheral blood is already driven mostly by functional HSC’s by 14 weeks post conception [23]. Constancy of the population size in the mature phase is ensured by applying a Moran model of divisions [34], where each non-differentiating symmetric division — which results in two daughter cells with the HSC phenotype — is simultaneously accompanied by a removal, either through a symmetric differentiation or some form of cell inactivation. Thus in our model we distinguish 3 types of events which occur independently with exponentially distributed waiting times:

1. Growth divisions, which are non-differentiating symmetric divisions (i) occurring at rate γ, and cause the population to increase in size by 1 cell;

2. Moran events, which consist of the simultaneous occurrence of a non-differentiating symmetric division (i) (which causes the population to increase in size by 1) and a cell removal (ii) (which causes the population to decrease in size by 1) either through symmetric differentiation or inactivation. These events occur at rate ρ;

3. Asymmetric divisions (iii), in which one of the daughter cells is differentiated and the other is not, occurring at rate ϕ.

In the early developmental phase there are only growth divisions; during the mixed phase, growth divisions occur alongside Moran events and asymmetric divisions; and in the mature phase, only the Moran and asymmetric division events occur. Note that under these dynamics the rates of the different types of divisions changes across phases: the rate of non-differentiating symmetric divisions (sometimes referred to as self-renewals) is given by γ during early development, γ + ρ during mixed growth, and finally ρ at maturity; while the rate of asymmetric divisions respectively jumps from 0 to ϕ in the mixed phase, and remains at ϕ in the mature phase. This can be summarized as:

We assume exponential growth (see section 1.4 for discussion of other possible growth curves) from the stem cell population’s most recent common ancestor (the zygote, or any later occurring stem cell), and make the simplifying assumption that all division rates remain constant in their respective range of operation. The population size then depends on the division rates as N(t) = exp (γ(t)t) with Θ(x) the Heaviside step function.

We introduce μ as the rate parameter of the Poisson distribution from which the number of mutations arising in a daughter cell is drawn, i.e. the mutation rate per daughter per division. Each new mutation is presumed to be distinct from all previously existing mutations in the system, an assumption known as the infinite sites model.

A summary of all parameters governing this model is given in table 1.

Table 1.

1.2 Stochastic simulation of the dynamical model

In direct simulations of the system we model individually distinct cells, each carrying a particular set of mutations. These are evolved according to the model described in Section 1.1 by means of a Gillespie algorithm [35]: Events are performed successively with exponentially distributed waiting times, and the acting cell in each event is selected uniform randomly from the population. Daughter cells arising from a division inherit all mutations from their parent cell, and acquire new mutations which are distinct from any existing mutations.

Simulations are initiated with a single cell carrying no mutations. In the early and mixed growth phases, taking the rate of growth divisions γ as fixed ensures exponential growth of the complete population. To investigate the effect of differing growth curves (see section 1.4), we also implemented a logistic growth phase, where cell division gradually decreases and cell death gradually increases until the population reaches its mature size. The average population growth then follows

The code for the simulations was written in the Julia computing language, and can be found at https://github.com/natevmp/hsc-vaf-dynamics.

1.3 The expected VAF distribution of the complete population

To obtain the dynamics of the VAF distribution, we first look at the probabilistic time evolution of a single variant’s size (i.e. the number of cells sharing that mutation) P_k(t), which is given by a Fokker-Plank equation. We first sketch the procedure for obtaining the result for a Moran model, where the population size is constant. Next, we show how to obtain the expression for a deterministically growing population. Finally we use these results to obtain the distribution of variants for a population that changes in size.

1.3.1 Probability density distribution of a single variant in a Moran process

The discrete probability distribution of variant’s size is governed by the transition probabilities ℙ{l| k}, which denotes the likelihood of the number of cells carrying this mutation changing from k to l. In a Moran model, a constant population size N is maintained by the simultaneous occurrence of one cell division and one cell death/removal. Thus, the possible transitions following the occurrence of such an event are given by the probabilities

Assuming events occur with exponentially distributed waiting times at rate ρ, P_k can be described through the master equation

When the population size N is large, this leads to a system of N differential equations, which is computational taxing to solve. Instead, a diffusion approximation can be made (on the assumption that N is large, i.e. that there are many states k = 0, 1, …, N), to describe p(κ), the probability density of variant’s size on the continuous support κ ∈ [0, N]. To this end a Fokker-Planck equation is constructed of the form

where the coefficients are given through the infinitesimal propagator t(κ + Δκ, t + Δt | κ, t):

The integrals are the first and second moments of a displacement Δκ, and can be obtained using the heuristic [36]:

with ⟨⟩_Δt the expectation taken over a time interval Δt. Typically a transformation to frequency space is first performed, f = ρ/N, where f is the frequency of cells carrying this mutation in the population. Thus, we obtain

and finally the well-known result

1.3.2 Probability density distribution of a single variant in a deterministic growing population

We now consider a growing population due to additional divisions occurring—at rate γ per cell— alongside the previously described Moran events. For simplicity we assume the population size N(t) increases deterministically according to the growth rate. The transition probabilities for the variant size are then given by

Using the previously described heuristic (13) the coefficients of the Fokker-Planck equation become

And the full equation

1.3.3 The expected VAF distribution in a varying population

For independently evolving variants, we denote the number of variants at size k (i.e. shared by k cells) as V_k(t). We now note that in the limit of infinite variants, the fluxes between adjacent k-states are identical to the above-derived transition probabilities. Thus the above-derived equations can similarly describe the expected dynamics of a an ensemble of variant sizes V(κ, t), where we once again introduce the parameter κ ∈ [0, N(t)] to represent the continuous space analog of the variant size k. In our system of interest we however also wish to include an injected flux of newly arising mutations, which is limited to the state κ = 1 and has amplitude

Note that because asymmetric divisions result in just a single daughter remaining in the system, they contribute only half as many novel mutations to the system as symmetric divisions. In the continuous domain we introduce this flux as a Dirac delta function 𝒞δ(κ−1). Note also that, while these diffusion equations are more commonly applied in the frequency domain f ∈ [0, 1], for the purpose of numerical simulation it will be convenient to remain in the space of absolute sizes κ = f N(t). Then finally, the time evolution of the VAF distribution is given by

which is the result stated in the main text.

1.3.4 Numerical approximation of the VAF distribution dynamics

While a solution to the fixed population Moran system (15) is known [37], it is unwieldy, and does not readily generalize to the more complex expressions (18) and (20). We therefore resorted to numerical approximations of their solution. We applied a method of lines approach with a finite difference discretization in the size coordinate κ. The incoming flux of variants at κ = 1 introduces a discontinuity in the derivative, which leads to stiffness in the system. This is the reason we eschewed the transformation to f = κ/N(t), since when working in frequency space the point of incoming flux would change position as the population grows. To maximize performance near the discontinuity, we implemented a variable step size smallest near κ = 1, with the distances in space given by

where i ∈ 1, …, n − 1, with n the number of discretized points, and

In this formalism the delta function in (1) was approximated as a step function with height 2/(Δκ₁ + Δκ₂).

1.4 Changing the growth function during the maturation phase

In the main text we have shown that, for a constant population, the past occurrence of a growth phase influences the shape of the observed VAF spectrum. Here we investigate to what extent the precise form of this growth function alters the observation. To this end we compare the expected VAF spectrum, calculated from Eq. (1), for two different growth models: One where the population grows exponentially until the maximal capacity N_M is reached, at which point it immediately stops expanding; and one where the population grows according to the logistic function, which is a common model for population growth in ecological systems. Since in the latter case the function only asymptotically approaches capacity, we take maturity to be the point where the population reaches 0.99N_M. The resulting growth functions and expected VAF spectra are shown in Fig. 5 for two time points measured after maturity. We observe that there is indeed a quantitative difference between the spectra of the two models. However, we note the same qualitative characteristics as described in the Results section of the main text: a pronounced move towards the fixed size shape in the lower frequencies, and a delayed response at higher frequency. The difference between the two models furthermore diminishes as the measurement occurs farther away from the end of the growth phase.

1.5 Analytical solutions for the expected VAF distribution in select cases

While in section 1.3 it was shown how the expected VAF distribution of our stochastic model can be calculated numerically, under certain (stricter) conditions some solutions can also be obtained analytically. We present here derivations for analytical solutions of the distribution in a constant population adhering to a Moran process, and a strictly growing population subject to the pure-birth process. Although the results derived in this section are themselves not novel (see e.g. [24], [25], and [6]), we present them here for completeness, with new derivations that facilitate a unified analysis and do not require readings of the previous literature. For convenience, a list of the parameters used in this section is given in Table 2.

Table 2.

1.5.1 The equilibrium VAF distribution in a Moran process

In a birth-death Moran process, one cell dies and one is born every time step. Then, let every birth coincide with an average of μ new unique neutral mutations. Every cell can hold any number of mutations, and since these are neutral, the cells that are dying/born are chosen completely at random. If two cells are chosen for death and birth, any particular mutation has a chance P_k to be in one of the cells but not in the other:

This was previously also shown in Eq. (9). If we want to know how the number of mutations changes, we need to know how many mutations increase or decrease in abundance on average. This simply means that we multiply P_k by the number of mutations with the given abundance k:

Note that because of the symmetry of the process, C_k(t) is both the average increase and decrease. Lastly, if we look at a certain abundance k, there are two ways to lose mutations and two ways to gain mutations: M_k increases if mutations of abundance k + 1 lose one cell or if mutations of abundance k + 1 gain one cell. M_k decreases if mutations of abundance k lose or gain an cell. Note that all of these things can happen at the same time. If we want to know the average change , we therefore need to incorporate all of these four possibilities. After we add the rate division rate ρ, we finally get the following equation for the average change per unit of time:

This is a variation of Eq. (10) applied to all mutations, as opposed to just a single one. For k = 1 the first term C₀ is simply replaced by a fixed mutation rate μ, for k = N, the middle term has to be removed (and the last term C_N(t) is always 0). Once we apply this to an entire distribution of mutations, it gives us an estimate of how, on average, this distribution changes every time step. Notably, it should approach some equilibrium state if we exclude the last state, which is absorbing. Let us consider the distribution from k = 1 to k = N−1. The only way for this distribution to lose mutations is by either a mutation of abundance 1 to die out, or a mutation of abundance N−1 to fixate. This implies an equilibrium at:

We can now apply this iteratively to the other C_k(t) (for convenience, we ignore the ρ, since it has no bearing on the equilibrium, merely on how fast we approach it):

Which can then be applied to C_N−1(t) itself:

This can now be used to solve for C₁(t), C₁(t) … C_N−2(t) iteratively, resulting in an equilibrium solution :

Finally, this can be transformed back into the V_k(t) to get the equilibrium solution V^*(k):

1.5.3 The VAF distribution of a growing population under a pure-birth process

In a pure-birth process only births, and no deaths, occur in the population. Suppose then a strictly growing population starting from a single cell increasing by 1 every time step, with a mutation rate of μ. We denote the number of mutations with abundance k at time t as V_k(t). Then we can set up a system for the average change per time step:

For simplicity we view the process in event time, so that every time step is exactly one division event. However, the process can be projected so that the amount of time δt that passes between every event is drawn at random from a given distribution. Since this projection is always possible in both directions now matter the chosen distribution, this proof is valid for any process, as long as births are not allowed to be simultaneous. The V_k(t) themselves are all strictly increasing every time step and hence approach infinity, but we can make a statement about the average shape of the distribution and the value at each specific point in time:

To prove that 2/(k ·(k + 1)) is indeed the expected distribution for this system, we will first show that the first several time steps follow this shape. Then we prove that every time a new cell and thus a new possible abundance k = t is added to the system, the new distribution at the new abundance will follow the expected distribution as well. Finally, we will show that if an abundance and the neighbouring abundance follow the distribution already, it will continue to follow the distribution.

Let’s show now that the first few time steps follow the claimed distribution:

Now, since the first steps seem to follow the distribution as expected, we can move on to showing that generally at time step t, the newly generated abundance k = t will always fit the distribution:

The only part remaining is proving that if the preceding time steps fit the distribution, the following time step also follows the distribution:

Note that unlike for the constant population case, the expected distribution does not need time to asymptotically approach the steady-state solution, as it is already at equilibrium from t = 1.

1.6 Sampling affects the observed VAF distribution

Eq. (1) describes the evolution of mutation frequencies with respect to the whole population, however, these are typically not known from sequencing data. Instead, an observation provides the frequencies of mutations as measured in the sequenced sample of the population. This distinction should be recognized when interpreting experimental data, as the statistical distribution for the sample’s VAF distribution differs from that of the true population (Fig. 3c). To obtain the expected distribution from a sample of S cells, we must consider how each mutation’s frequency in the true population can be altered in the sample, so that the resulting curve is given by the sum over all original mutation frequencies j/N weighted by their probability of being sampled to a different frequency i/S. If sampling is performed randomly this probability is given by the hypergeometric distribution (sampling without replacement), leading to the expression

Besides the anticipated loss of resolution, the effect of this operation is to shift the expected mutation frequencies upward with respect to their true frequency in the population (see Fig. 3c in the main text), which can lead to improper inferences if not correctly accounted for.

1.7 The mutational burden distribution

1.7.1 The expected mutational burden distribution

The expected distribution of mutational burdens can be obtained from the argument that the burden m_j of an individual cell at the time of observation must equal the sum over all mutations that occurred during each of its past cell divisions:

Without knowledge of the distributions of y and u other than the requirement that the u_ji are identically distributed and independently drawn, it is straightforward to show that the expectation and variance of m_j are given by E(m) = E(y) E(u) and Var(m) = E(y) Var(u) + E(u)² Var(y) (Supplement). With the number of novel mutations per division as Poisson distributed, the average number of divisions per lineage — expressed in terms of an effective division rate λ — can thus be estimated from a sample as

If the mutation rate μ is not known, it can similarly be obtained from

Note that while a common modeling practice is to take the number of divisions as Poisson distributed — i.e. assuming exponentially distributed times between divisions — the estimates (4) and (5) do not require this. In fact, the factor E(y)/ Var(y) in (5) can be interpreted as a measure for the under-or over-dispersion of divisions per lineage with respect to the Poisson model where E(y) = Var(y). In other words, if the Poisson model is assumed (then (40) describes the Compound Poisson distribution (CPD)) the mutation rate can be estimated directly from the data m, whereas if the mutation rate is known the accuracy of the Poisson model can be assessed. Finally, note that in the true system there are nonzero correlations between lineages – the early u_ji are shared amongst many cells – meaning that for a single population the m_j are not truly drawn from (40). While this does not affect the expected division rate (4), an inferred mutation rate will be influenced by a shared history, as it effectively reduces Var(y) and Var(u) (Fig. 3b).

1.7.2 Fluctuations of the mutational burden distribution

To investigate expected within- and between-patient variation of the mutational burden distribution we implemented individual-based stochastic simulations of single cell mutation acquisition. We start simulations with N cells and 0 mutations. Cells are picked randomly for division and acquire a number of mutations drawn from a Poisson distribution with rate parameter μ. In addition, cells may undergo symmetric or asymmetric divisions with a given probability p. If a cell divides asymmetrically, the novel mutations are added to that cell in the simulation. If a cell divides symmetrically, one cell in the population replaces another and novel mutations are added to the replaced and the dividing cell.

These simulation allow us to track the mutational burden distribution and its variance for arbitrary p, μ and large N over time. Realisations of these dynamics for μ = 1.3 per cell division are shown in Fig. 7 and are compared to the experimental observation of the mutational burden distribution in [18]. We observe significant stochastic variations, especially pronounced for small stem cell population sizes N and low asymmetric division probabilities p. However, the true biological parameters are more likely to be a combination of large N and large p and therefore a mutation rate of μ = 1.3 per cell division is unlikely to explain the observations in [18].

In Fig. 7 and Fig. 8 we show the likelihood that the ratio of Variance/Mean is within the experimentally observed interval 3 < θ < 5. Again we can see that for μ = 1.3 the ability to explain the experimental observations is limited for any combination of parameters. Increasing the mutation rate to μ = 3 in contrast can explain the experimental observations for a wide range of parameter combinations.

1.8 Data fitting pipeline

To estimate values of the model parameters from the data, we combine both prediction for the single-cell mutational burden as well as the mutation allele frequency distribution. The unknown parameters are those listed in Table 1, whereby the growth rate is implicitly defined from the mature population size and the maturation time as γ = log(N_M)/t_M.

1.8.1 Estimating division rates from the single-cell mutational burden distribution

While the single-cell mutational burden distribution allows for direct estimation of the average number of divisions that occurred per lineage through Eq. (4), this does not directly provide the division rates, as these combine in the integrand on the lhs. To perform the integration we must first determine the expected contributions of each type of division (growth, symmetric, or asymmetric) to the mutational burden of a single lineage, as they — perhaps counter-intuitively — do not influence the expected burden equally.

Consider the system of N cells which obey the dynamics described in Section 1.1. At any point in time each cell describes a single lineage, and we are interested in the expected number of mutations added per lineage over time. Since divisions occur as independent Poisson events, we may treat the distinct types individually and simply add their contributions together to obtain the total burden increase over time. When an asymmetric division occurs the situation is simple: The single stem-type daughter cell acquires on average μ new mutations, which increase the average burden per cell by μ/N. Since these dynamics are Markovian, and events occur at rate Nϕ in the population, we obtain the average number of mutations acquired per lineage in a time Δt as

When a symmetric division occurs (as performed in the Moran process described previously) the situation becomes slightly more complicated. Following a single division we discern two disjoint outcomes: either two different cells are selected for replication and removal, or the same cell is selected for both. The former case occurs with probability 1−1/N and results in the average burden per lineage increasing by 2μ/N (two daughter cells obtained on average μ new mutations), while the latter occurs with probability 1/N and causes the average burden to increase by μ/N (one of the daughter cells was removed).With these events occurring at rate Nρ the expected increase per lineage after a time Δt is:

with the final approximation obtained in the limit of large N. For the pure growth divisions occurring at rate Nγ we take a similar approach, though we now consider an infinitesimal time increment dt to account for the varying population size. After a single division the average burden increase is 2μ/(N(t) + 1), so that

Thus, we find that the growth and Moran divisions on average add twice as many mutations to a lineage as the asymmetric divisions. Note that we obtain the average increase in lineage (the number of divisions in a cell’s past) by setting μ = 1 in the equations above. Plugging these results into the integral of (4), we find for the expected number of divisions per lineage at time t

Taking into account the different phases of evolution we can separate the integral into three parts

Denoting λ = ρ + ϕ the total rate of homeostasis divisions within the stem cell compartment, performing the integration we solve this for

with R = E(m)/μ the average number of divisions per lineage estimated from the mutational burden distribution. With μ estimated from [5], the remaining unknown parameters are t_M, N_H, N_M, and p. Finding the standard error on E(m) as (with S the sample size) and sweeping the unknown parameters within reasonable limits (t_M ∈ [0, 10], N_H ∈ [0, 1000], N_M ∈ [20^′000, 500^′000]) gives an error interval for the expected λ.

1.8.2 Fitting to the lowest frequency of the VAF distribution

We must be cautious about attempting to directly fit realizations of Eq. (1) to an observed distribution in its entirety. Without knowledge of the variance across individual realizations for each possible allele frequency, there is no rigorous way to weight the errors from different frequencies to obtain a single best fit. Indeed, simulations suggest that the variance on the number of mutations per VAF grows with the frequency (Fig. 9), implying that a non-weighted fit would overemphasize fluctuations at one end of the spectrum. We therefore selected only the lowest frequency available in the sample – which is anticipated to be the best described by the expected value (Fig. 9b) – as the reference for parameter fitting.

For select values of t_M and N_H, we calculated the VAF distribution — by numerically solving Eq. (1) and applying the sampling operator Eq. (6) — at lattice points in the parameter space of N_M ∈ [5 · 10³, 5 · 10⁵]× p ∈ [0, 1]. For each point (N_Mi, p_i) we compared the lowest frequency of the theoretical expected distribution v_1/S with the value measured in the data, resulting in a 2D lattice of errors. From this lattice a continuous error function was obtained through spline interpolation, wherein we identified a zero-error line, which is the line of optimal fit shown in Fig 4c.