Stochastic modeling indicates that aging and somatic evolution in the hematopoietic system are driven by non-cell-autonomous processes

Conceptual basis of the model

We assume that all possible random somatically heritable genetic and epigenetic changes, collectively referred to hereafter as mutations, have a theoretical distribution of fitness effects (DFE) per cell division. Neither the shape nor the variance of this DFE is known for HSCs. At the organismal population level, most mutations have either no or negligible phenotypic effects or are lethal (reviewed in [50]). Mutations that reduce fitness are much more frequent than those that confer a selective advantage. However, as somatic cells are subject to selection at both the animal population and tissue levels, the net DFE of somatic mutations that affects inter-cellular competition and somatic evolution at the tissue level could differ from DFEs typically observed at the animal population level.

To infer the nature of DFE in HSC pools, we applied the logic that mutation fitness effects will affect competition and selection processes by creating a fitness differential in a population of cells. A narrow DFE with mostly neutral mutations will result in a largely drift-driven population, where the chance of survival is random (Fig. 1; upper panels). A wide DFE, where a larger proportion of mutations have fitness effects, will engender a broader fitness differential among cells within a population. This fitness differential lends itself to natural selection, wherein cells with disadvantageous phenotypes are less likely to survive into subsequent generations. As most functional mutations are likely to decrease fitness (residing in the negative tail of the DFE), we reasoned that variance of DFE should have an effect on the slope of mutation accumulation in HSCs, whereby wide variance (more functional mutations) should act to purge mutated cells with decreased fitness from the pool and thus act to lower the speed of mutation accumulation in the population (Fig. 1; lower panels).

Mutation accumulation in the Tier 3 genome can be used as a marker for overall mutation accumulation rates. Tier 3 represents (to best approximations) the non-conserved, non-coding, and non-repetitive sequences of the genome, spanning ~1 billion base pairs [51]. The vast majority of Tier 3 DNA should therefore not be under selection (particularly in somatic cells), and mutations in Tier 3 should (with rare exceptions) not alter cellular phenotype. But the accumulation of these passenger mutations will depend on the intensity of purifying selection acting on functional mutations (including epigenetic alterations) as shown in Fig. 1, and can be used to infer the intensity of selection and thus the parameters of mutation DFE acting on HSC. Basically, increased accumulation of mutations with detrimental phenotypic effects will confer a fitness cost to cell division. We will use our model to determine the parameters of mutation DFE and rate increase over a human lifetime that can replicate experimentally determined mutation accumulation in the Tier 3 genome of HSC.

Since the most common leukemias appear to initiate in HSC [52-56], the initial clonal expansions of an oncogenically-initiated HSC should be rate-limiting for the occurrence of successive oncogenic events in the clonal context and eventually for the development of hematopoietic malignancies. We will thus also ask if the estimated realistic parameters for mutation DFE are concomitantly permissive for age-dependent exponential increases in the extent of clonal expansions in the HSC pool. The rates of somatic evolution (extent of mutation-driven clonal expansions) is a major factor defining the ultimate odds and frequency of leukemia over time/age, as can be inferred from the following equation integrating the probability that a set of n successive oncogenic drivers {d1, …, dn} will happen in any single cell in a stem cell pool over a lifetime (Eq. 1): $P_{d1\dots dn}\left(t\right)=\text{D}\left(t\right)\times {\displaystyle \underset{0}{\overset{t}{\int }}\left({\displaystyle \prod _{i=1}^n{p}_i}\right)}\left(t\right)dt$(1) where P_d1…_dn(t) is the probability of acquiring n drivers in one cell by time t, D(t) is the total number of cell divisions that have happened in the pool by time t, and p_i is the probability of acquiring a driver i ∈ {d₁, …, d_n} per cell per division as a linear function of the total number of mutations per cell per division.

The chances of acquiring a set of n drivers in any given cell in the pool over time depend on two factors: mutation rate and the total number of cell divisions by the time in question. From Eq. 1, it is evident that cell proliferation driven by acquisition of a fitness advantage will have a much more dramatic effect on the final probability of the whole set of drivers occurring within one genetic context compared to mutation rate. While a change in mutation rate can lead to a linear increase in this probability, the expansion of a selectively advantageous clone will elevate the probability of occurrence of subsequent driver mutations within a cell of this clone/genetic context exponentially. As Eq. 1 defines the probability density function of n drivers happening in any cell in the pool over time, clones making up a greater share of the pool will harbor proportionally more dividing cells and have proportionally higher chances of this set of drivers happening in a cell within the clone. Based on this logic, we assume that the shape of the age-dependent incidence of leukemia will mostly be determined by the age-dependent magnitude of selection-driven clonal expansions possible under given parameters for mutation DFE and rate. Therefore, we asked what mutation parameters are compatible with both the reported slope of mutation accumulation in the Tier 3 genome and with exponential increases in the magnitude of clonal expansions (increased positive selection) that replicate the shape of the age-dependent leukemia incidence curve.

Architecture of the model

To fully investigate the many parameters governing somatic evolution in HSC pools, we designed a stochastic model to replicate HSC population dynamics, to simulate the impact of mutations in HSCs over human lifetimes, and to model the effects of tissue microenvironment on selection and clonal expansion within the HSC pool. This model is a stochastic, discrete time continuous parameter space model realized in a Monte Carlo experiment run in the Matlab programming environment (The MathWorks, Inc., Natick, Massachusetts). A chart of cell fate decisions in the simulated HSC pool during a model's run is shown in Fig. 2. The model starts with a matrix of the initial number of HSC, and each cell's state is updated on a “weekly” basis throughout the simulated lifespan of 85 years (4420 weeks). The weekly update included stochastic cell fate decisions to divide or stay dormant based on estimated HSCdivision rates throughout life (modeled based on published data; Fig. 3), and to stay in the pool or leave for whatever reason (such as death or differentiation) based on niche space availability, current number of competing cells at different ages (modeled based on published data for HSC pool size; Fig 3), and each cell's relative fitness. Cell fitness changed after each division stochastically, initially based only on mutation DFE. Cells that diverged in fitness more than a certain threshold from their parental cells upon division were designated as new clones, thus replicating functional (clonal) divergence of HSC in the pool with age (Fig. S1A-B). The code and detailed parameter description and justification are presented in Supplemental Methods.

Fixed fitness effects of mutations cannot explain age-related functional decline and somatic evolution in HSC pools

We found mutation DFE variance (in standard deviations, denoted as σ) to influence the slope of mutation accumulation within the range σ= 5×10⁻¹–5×10⁻⁶ and mutation rate increase from stable to up to 8-fold over lifetime (Fig. 4). As described above and outlined in Fig. 1, higher σwill suppress mutation accumulation, including in the Tier 3 genome, by essentially imposing a penalty on cell division. An increasing mutation rate effectively widens DFE variance per cell division following the rule of distribution sums [57] (Fig. S1C-D). We used these ranges for mutation rate increase and DFE variance over lifetime as the X and Y axes of a plane, with the Z axis represented by the model output of interest (more details on generating the surface plots can be found in Supplemental Methods). In this manner, we built surfaces of model outputs of interest (such as slope of mutation accumulation or magnitude of clonal expansions over lifetime) within the effective range of mutation DFE variance and mutation rate increase. Mutation DFE was always centered on zero, with the majority of mutations being neutral. Within this principal design of measuring the simulated model's output, we also tested a range of mutation DFE shapes (DFE tail ratios), with 0, 1, 10, 33, and 50% of mutations in the positive tail. A larger positive tail means a greater frequency of beneficial mutations. Note that even when 0% of mutation DFE is in the positive tail, fitness differential will still build up within the population (with some cells “less negatively” affected). We assessed three basic model outputs: 1) the Tier 3 mutation accumulation slope; 2) the maximum extent of somatic evolution possible at any given age measured as the percentage of the pool occupied by the most successful clone; and 3) the dynamics of average cell fitness in the pool, to see if the model replicates the general stem cell functional decline over lifetime. The reference Tier 3 mutation accumulation slope with age for AML (0.09162; 95% confidence interval, CI: 0.03759-0.1457; Fig. 4A; The Cancer Genome Atlas Research Network, 2013) was used to delineate the range of mutation DFE and rate parameters that reproduced the slope of DNA mutation accumulation in the Tier 3 genome of HSC. As argued by others [58], mutation accumulation in AML is thought to reflect mutation accumulation in individual HSC (see Supplemental Methods). We will hereafter call this range the plausible range, depicted as shaded areas in Fig. 4F and 5A,B, to signify the ranges of mutation DFE variance and rate parameters that are most likely to enclose real parameters. Notably, the slope of DNA methylation accumulation in a set of CpG sites in blood cells ([19]; 95% CI: 0.061-0.062; Fig. 4B) is well within the CI for Tier 3 mutations (Fig. S2), indicating that mutation accumulation rates for genetic and epigenetic changes may be similar (both being consistent with cell division kinetics; Fig. 3A). We registered maximum rates of somatic evolution throughout the lifespan to see what parameters of the mutation process are permissive for age-dependent exponential increases in the magnitude of mutation-driven clonal expansions that would resemble the age-dependent increase in leukemia incidence. As evident from Eq. 1, the frequency of leukemia should closely follow the extent of somatic evolution, and somatic evolution rates, in turn, are dependent on the frequency of functional phenotype-altering mutations, i.e. on the parameters of the mutation rate and DFE. Typical age-dependent curves of clonal expansions generated by the model under different mutation DFE are shown in the shape-match plots with leukemia incidence in Fig. 5C. Lastly, the average fitness decline was monitored to test if the experimentally observed ~2-3-fold decline in HSC fitness over lifetime [5, 59, 60] was replicated in the model output.

Results shown in Fig. 5B (see also Table S1A, mutation-only conditions) reveal a substantial discrepancy between conditions required to replicate re-alistic mutation accumulation (shaded areas in Fig. 5A,B) and those that replicate clonal expansions (rates of somatic evolution) resembling the age-dependent shape of leukemia incidence. We used 0-1 normalization of leukemia incidence and the simulated clonal expansion curves to compare their shapes by the mean root square error (MRSE) method, with shape similarity between 0 (no similarity) and 1 (perfect match) (see color code in Fig. 5C).

A maximum of ~5% overlap was observed between the plausible range of mutation DFE and rate parameters replicating the empirically known slope of mutation accumulation (shaded areas in Fig. 5A,B) and the DFE range permissive for leukemia incidence-like shape of clonal expansions (red areas in Fig. 5B) under a liberal cutoff of ≥0.7 of shape similarities by MRSE. See Fig. 6 and Table S1A, “Mutations alone” conditions, for quantitative data. This overlap (the common set of mutation parameters needed to replicate both reference phenomena) only occurs in the middle range of mutation DFE variance (Y axis in plots in Fig. 5B) and requires a ≥4-fold mutation rate increase over lifetime (X axis in plots in Fig. 5B). As shown in Fig. 4D, under such combinations of mutation parameters, mutation accumulation is impeded or even reversed at later portions of life. The increasing mutation rate effectively increases DFE variance per cell division (see explanation Fig. S1A-B) and thus intensifies selection against mutant cells. Such a late-life slowdown in mutation accumulation or mutation purging has not been observed experimentally, and therefore we consider these combinations of mutation rate increase and DFE variance to be unrealistic.

Surprisingly, both the above-mentioned area of overlap and the magnitude of clonal expansions are drastically suppressed by an increased positive tail (more beneficial mutations) of the mutation DFE (Fig. 5B, Fig. 6, and Table S1A, “Mutations-alone” conditions). This suppression likely arises from the phenomenon known from bacterial populations as clonal interference, whereby an increasing number of high fitness clones mutually suppress the expansion of each other and that of other phenotypes [61, 62]. In our simulations, the strongest fitness gain (clonal expansions) over time based on sole action of mutations was possible under a zero to minute (1%) positive tail of the mutation DFE, consistent with the idea that stem cells reside at a local fitness peak (low probability of advantageous mutations), at least in young healthy individuals [34]. Therefore, the plausible range of mutation parameters is not permissive for significant exponential age-dependent increases in the rates of somatic evolution (clonal expansions).

The average cell fitness drop in the pool throughout lifespan also did not replicate the experimentally reported 2-3-fold decline in fitness for HSC in old age [5, 59, 60] and was restricted to a maximum of a few percent drop over lifetime throughout the whole range of conditions (Fig. 5D). Despite the mostly negative fitness effects of mutations, purifying selection appears to buffer the general cell fitness decline by purging mutation-affected cells from the pool. This effect of selection is particularly evident under wider mutation DFEs (Fig. 5D). Thus, accumulation of cell-autonomous damage is insufficient to account for HSC fitness decline in old age.

Mutation fitness effects modified by microenvironment explain higher late-life rates of somatic evolution and fitness decline in HSC pools

Since we were unable to define a common set of mutation parameters that would recapitulate the experimentally observed mutation accumulation rates, somatic evolution/leukemia incidence, and HSC fitness decline with age, we tested an alternative, evolutionary model. For this purpose, we developed a model of bicomponent fitness effects exerted by the tissue microenvironment, containing a uniform and a randomly distributed part. The uniform component derives from the loss of tissue integrity with age, which should impact all cells in a tissue, as tissue microenvironment moves from its optimum towards a more degraded state. Roughly, this component dictates that cells generally have a lower fitness in a degraded (aged) microenvironment compared to the optimal (young) one. However, the degree of this influence should vary among cells based on the mutation-generated phenotypic diversity of cells. This variation is the distributed part of the microenvironmental effect that adds to selection processes by contributing to the fitness differential buildup in the cell population. We propose that this part of the environmental effect is distributed independently from the mutation DFE, so that any given phenotype that has a fitness advantage in an optimal microenvironment will not necessarily have this advantage in an altered environment (and vice versa), thus making the composite fitness effect of any given mutation context-dependent, just as it is in natural populations of organisms. Such a bicomponent effect of the tissue microenvironment is based on inference from the Sprengel-Liebig's system of limiting factors, initially known as Liebig's Law of the Minimum (summarized in [15]).

The principle for how the Sprengel-Liebig's system defines the fitness of a given phenotype is shown in Fig. 7. In aggregate, any given phenotype has a certain degree of adaptation to each environmental factor so that the factor has its optimal intensity (e.g. optimal concentration) and extreme intensities, also called pessima (Fig. 7A). Fitness decreases as the intensity of a factor in the environment changes from the optimum to either of the pessima for a given phenotype, a phenomenon known as Shelford's Law of Tolerance [63]. Selection leads to improved adaptation so that populations consist of phenotypes for which a given environment is optimal. We argue that selection at the germline level should lead to co-evolution of stem cells and tissue microenvironment to optimize performance during pre- and reproductive periods (Fig. 7B). In a complex, multi-factorial environment a phenotype will have different degrees of adaptation to each particular factor, but the phenotype's net fitness will be limited by its adaptation to the factor it is least adapted to, called the limiting factor, following Sprengel-Liebig's Law of the Minimum (Fig. 7C) ([15, 64]). General adaptation to a complex environment leads to the evolution of phenotypes that have optimal net fitness (“evolved phenotype” in Fig. 7C), and this process reduces the likelihood that new phenotypes that arise from random mutations will improve fitness relative to the population as it becomes better adapted.

Environmental alterations lead to changes in factor intensities, which results in perturbation of the relative fitness of different phenotypes within a population and leads to selection for minority phenotypes that happen to have a better net fitness in an altered environment (Fig. 7D). Alteration of an environment reduces the average fitness of the population (Fig. 7E,F), revealing a uniform fitness effect on the population. However, in a phenotypically diverse population the exact degree of fitness change for any particular phenotype will depend on the phenotype's set of adaptations to particular environmental factors as shown in Fig. 7C,D, revealing a stochastic component of environmental effects which redistributes the relative fitness of particular phenotypes in a population independently from their initial distribution in an optimal environment. This general principle is the mechanism for how phenotypic effects of mutations are translated by environment into the phenotype's fitness.

To model this principle, we first generated a proposed curve of fitness decline, which describes a general non-linear ~3-fold fitness drop as a function of age using the following equation: $F\left(A\right)=F_{\mathit{max}}-\left(\frac{F_{\mathit{max}}-F_{\mathit{min}}}{1+5200A{e}^{-0.0031}}\right)$(2) Where F is average HSC pool fitness, F_max is maximum initial fitness equal to 1, F_min is minimum end-life fitness equal to 0.3, and A is age in weeks. This curve was not intended to exactly replicate the natural fitness decline, as its shape is unknown, but reflects the general principle shown in Fig. 7B [65]. The overall 3-fold fitness reduction was chosen based on the 2-3 fold reductions in hematopoietic output per HSC observed for older humans and mice [5, 59, 60]. We used the first order derivative of this function to calculate the average fitness drop (ΔF) for any given discrete period. This amount was subtracted from each simulated cell in the pool at each weekly model update and represented the uniform part of the microenvironmental effects of cell fitness. To reproduce the distributed part of environmental effects, we then corrected the resulting fitness of each cell by an amount drawn from a normal distribution centered on zero and having variance σ (t) = ΔF(t)/2. As the chemical composition and physical conditions of tissue microenvironment are complex, their change with age will impact the existing mutation-generated functional diversity through an independent stochastic component, contributing to diversification of the relative fitness of different cells in a stochastic manner (Fig. 7F). In this way, mutations accumulated in a simulated cell lineage at any age are re-evaluated at each ”weekly” update of the model run; for example, mutations occurring during ontogeny could have an impact on HSC fitness later in life that differs greatly from their impact upon first occurrence (as the result of microenvironmental perturbations).

We applied this additional bicomponent modification of mutation fitness effects by microenvironment and tested the same series of outputs in a composite model, varying the proportion of the microenvironment-imposed DFE (ENV DFE) that resides in the positive tail. The composite model results in an extensive overlap between the plausible range (shaded areas in Fig. 8A,B) and conditions allowing for exponential age-dependent increase in clonal expansions (red areas in Fig. 8B; see color bar in Fig. 8C). Moreover, the magnitudes of these clonal expansions are substantially higher in the composite versus the mutations-only model (Fig. 6; Table S1). Both the extent of overlap and the magnitude of clonal expansions are maximal for a small (0-1%) positive tail for microenvironment-imposed DFE. As shown in Fig. 8D, the overwhelming cumulative effect of the microenvironmental uniform component results in the expected substantial average fitness drop, but without promoting excessive purifying selection to buffer fitness decline. In this way, the bi-component model of tissue microenvironment effects, by interacting with the mutation pre-generated cellular diversity, creates a range of mutation DFE variance and rate increase that are concurrently permissive for realistic mutation accumulation, exponential age-dependent increase in rates of somatic evolution (clonal expansions), and realistic cell fitness decline.

Given different estimates of the size of adult human HSC pool, we also tested if the results shown in Fig. 5 and 8 are sensitive to pool size. Both with a larger HSC pool kept at a stable size through adult ages and a pool size increasing over the whole life, the composite mutation/microenvironment model demonstrates a much greater overlap between DFE parameters that replicate leukemia incidence-like clonal expansions in the simulated pool and the slope of mutation accumulation in the Tier 3 genome (Fig. 6B and S3; Table S1B), suggesting that the principal discrepancy between the mutations-only (Fig. 5B) and composite (Fig. 8B) models is independent of pool size. We also tested the same model output under a lower adult cell cycling speed, as the lowest published estimates based on telomere attrition dynamics suggest ~0.6 divisions per year for human HSC [21]. The results in Fig. 6B and Table S1 demonstrate that the discrepancy in the mutations-only model still persists and is corrected by the composite model.