Homeostasis limits keratinocyte evolution

Abstract

Recent studies have revealed that normal human tissues accumulate many somatic mutations. In particular, human skin is riddled with mutations, with multiple subclones of variable sizes. Driver mutations are frequent and tend to have larger subclone sizes, suggesting selection. To begin to understand the histories encoded by these complex somatic mutations, we incorporated genomes into a simple agent-based skin-cell model whose prime directive is homeostasis. In this model, stem-cell survival is random and dependent on proximity to the basement membrane. This simple homeostatic skin model recapitulates the observed log-linear distributions of somatic mutations, where most mutations are found in increasingly smaller subclones that are typically lost with time. Hence, neutral mutations are "passengers" whose fates depend on the random survival of their stem cells, where a rarer larger subclone reflects the survival and spread of mutations acquired earlier in life. The model can also maintain homeostasis and accumulate more frequent and larger driver subclones if these mutations (NOTCH1 and TP53) confer relatively higher persistence in normal skin or during tissue damage (sunlight). Therefore, a relatively simple model of epithelial turnover indicates how observed passenger and driver somatic mutations could accumulate without violating the prime directive of homeostasis in normal human tissues.

Keywords: carcinogenesis; keratinocyte biology; mathematical modeling; somatic evolution.

Fig. 1.

Homeostatic epidermis model with high-resolution genomes. The homeostatic epidermis is composed of individual keratinocyte cells that are governed by a set of decisions (flowchart). Spatial structure of the 3D model and a growth factor (GF) diffusible gradient governs the loss/replacement in the stem-cell niche at the basal layer. The model is entangled with DNA sequence data by the Gattaca step, where base-pair mutations can be introduced with every in silico cell division. These mutations act as cell-fate markers and different subclones are represented by different cell colors. Therefore, cells in the simulated skin produce a patchwork of different-sized subclones, where each colored population differs by at least one mutation.

Fig. 2.

Homeostasis imposes log-linear clonal distributions. Neutral model dynamics from 3D simulations of various sizes where A is the cumulative clonal area frequencies of all patient biopsies (blue) and randomly chosen simulations for each of the patient-specific biopsies such that the number of simulations chosen is equal across the three simulated sizes and the sum of the number of biopsies equals that of the total biopsies for the patient (red). (Inset) The log-10 transformed first incomplete moment for the same random sampling of patient comparable simulations for PD20399 (error bars denote SD for 100 repeated samplings). (B) The difference from the Komlogrov–Smirnov test statistic (D_m,n) and critical value (D_α) for all patient biopsies to patient-specific model simulation’s first incomplete moment distributions for all simulation sizes. The red arrow denotes comparisons where the null hypothesis can be rejected.

Fig. 3.

Homeostasis, a priori, constrains clones in a functionally heterogeneous tissue, dictating an age-dependent clonal expansion. Clonal dynamics in homeostasis are a function of complex interactions between the microenvironment and the external environment. During homeostasis, every cell is equal relative to its neighbors, i.e., neutral. NOTCH1 disrupts neutral dynamics by offering a slight advantage via a blocking probability ($f_0$) which prevents neighboring cells from dividing into their positions (A), whereas TP53 mutations are not subject to UV damage, when a proportion ($\theta s$) of non-TP53 mutant cells may be killed by UV damage (B). (A and B) Next to each schematic are shown clonal expansions for a single simulation for NOTCH1-advantaged clones, blue, and TP53-advantaged clones, red, in a 1-mm² simulation to 58 y, respectively. The age dependence is shown in the density plots broken by mutation age for 100 replicate simulations up to 58 y for neutral (C), NOTCH1, with varying strengths ($f_0$) (D), and TP53, with varying strengths ($\theta s$) (E).

Fig. 4.

Combined effects reveal homeostatic recovery and recapitulate patient clone distributions. The primary, central plot reveals clone sizes for 1-mm² simulations for four different ages with five replicates across a manually defined scale of parameter pairs for NOTCH1 ($f0$) and TP53 ($\theta s$) advantages (n = 2,880 simulations). Density subplots show an idealized example of how each mutation type effects homeostasis (red and blue for TP53 and NOTCH1, respectively). The dashed boxes and lines show the evolutionary dynamics for a single simulation for a 58-y simulation (population frequencies greater than 0.001) for six parameter pairs and the corresponding overall population sizes over time (1–6 yellow circles, blue line is a locally weighted smoothing line). For a single parameter pair ($\theta s$ = 0.03, $f0$ = 0.05) a 0.4-mm² tissue is shown where only NOTCH1 (blue) and TP53 (red) clones are colored.

Fig. 5.

Combined effects adhere to age-dependent expansions. All 58-y 1-mm² paired parameter simulations from Fig. 3 are displayed where TP53 persistence ($\theta s$) is greater than NOTCH1 persistence ($f0$) (A) and $f0$ > $\theta s$ (B) (for $f0$ = $\theta s$ see SI Appendix). Individual clones are broken into three groups (passenger, NOTCH1, and TP53) with corresponding sizes of each clone.