Multiscale modeling of layer formation in epidermis

The mammalian skin epidermis is a stratified epithelium composed of multiple layers of epithelial cells that exist in appropriate sizes and proportions, and with distinct boundaries separating each other. How the epidermis develops from a single layer of committed precursor cells to form a complex multilayered structure of multiple cell types remains elusive. Here, we construct stochastic, three-dimensional, and multiscale models consisting of a lineage of multiple cell types to study the control of epidermal development. Symmetric and asymmetric cell divisions, stochastic cell fate transitions within the lineage, extracellular morphogens, cell-to-cell adhesion forces, and cell signaling are included in model. A GPU algorithm was developed and implemented to accelerate the simulations. These simulations show that a balance between cell proliferation and differentiation during lineage progression is crucial for the development and maintenance of the epidermal tissue. We also find that selective intercellular adhesion is critical to sharpening the boundary between layers and to the formation of a highly ordered structure. The long-range action of a morphogen provides additional feedback regulations, enhancing the robustness of overall layer formation. Our model is built upon previous experimental findings revealing the role of Ovol transcription factors in regulating epidermal development. Direct comparisons of experimental and simulation perturbations show remarkable consistency. Taken together, our results highlight the major determinants of a well-stratified epidermis: balanced proliferation and differentiation, and a combination of both short- (symmetric/asymmetric division and selective cell adhesion) and long-range (morphogen) regulations. These underlying principles have broad implications for other developmental or regenerative processes leading to the formation of multilayered tissue structures, as well as for pathological processes such as epidermal wound healing.

Initially we begin with 50 basal cells randomly placed on the basement membrane, and each cell is set with a random initial time. The domain is set with periodic boundary condition in the horizontal direction, and no flux boundary condition on the basement membrane. For simplicity, we assume that the signal cannot diffuse beyond the top layer of the epidermal tissue. Rather than defining no flux boundary conditions on complex surfaces, we instead extend the computational domain beyond the domain containing the cells and assign the chemical diffusion coefficient to be D c =0 on the extended domain. In addition to simplifying boundary conditions, this also allows the use of a rectangle shaped domain, which allow efficient computation.
Inter-cellular potential =0.05, =4.5, cut-off distance is 10.0. ! = ! = 2 for Base and Asymmetric Division Model; ! = 4, ! = 1 for Signal Model. External force !"#!$%&' =0.001, cut-off distance is 5.0. Probability Use values from Model 1 and 2 of The 3D model contains stochastic components from the individual cell activity, and we computed the mean and standard deviation of observed quantities. We determined the number of runs of simulations when the values of mean and standard deviation remain relatively stable as the number of runs of simulations increases. We also considered balancing modeling performance (more runs of simulation provide a more accurate prediction) and computational costs (more runs of simulation is more expensive). Fig A below shows Sharpness Index of each model with 25 simulation runs, and the difference of SI between 20 simulation runs (Figs 5,6,9) and 25 simulation runs (Fig A) are small. We found the summary statistics based on 20 simulations provides a "convergent" estimate of the interested quantities, along with affordable computational cost.

D. Stratification measured using Ripley's K function
Ripley's K function is a good measure to analyze point patterns at multiple distances [6]. Here we use it as an additional measure to investigate the spatial pattern of the same type cells at multiple distances in each slice along the z direction to compare with the investigations using Sharpness Index and Isolation Ratio. We have plotted the Ripley's K function for a typical simulation of each model (Fig C-G). The results of Ripley's K function show consistency with the analysis using Sharpness Index and Isolation Ratio.
In the Base Model (Fig C), Ripley's K function calculations show that every layer along the z-axis is a mixture of three cell types. Basal stem cells form clusters at short cell distance, which is a natural result of self-proliferation. Spinous and granular cells distribute evenly across the tissue. This analysis shows consistency with the results of Sharpness Index and Isolation Ratio (Fig 3, 5AB).
In the Asymmetric Division Model (Fig D), Ripley's K function calculations show that basal stem cells distribute regularly in the first layer with attachment to the basement membrane, while spinous and granular cells distribute regularly in the tissue except in the first layer, which is similar to the observation based on Sharpness Index and Isolation Ratio (Fig 4B, 5CD).
In the Selective Adhesion Model, for the good scenario (Fig E) Ripley's K function calculations show that the mechanism yields the stratified pattern. However, in the bad scenario (Fig F), there is a cluster of spinous cells above the granular layer. Both cases are consistent with the analysis using Sharpness Index and Isolation Ratio for the Selective Adhesion Model ( Fig  6CD, B).
In the Signal Model (Fig G), Ripley's K function calculations show that the mechanism yields the stratified pattern, which is consistent with the analysis using Sharpness Index and Isolation Ratio for the Signal Model (Fig 9BCD).
We also checked the pair correlation function for each model, and they provide the similar results. For example, Fig H presents the pair correlation function for a typical simulation of the Signal Model, showing the tissue is stratified and the same type cells are distributed more regularly within their layer. The small peaks when cell pair distances equal 1.5 cell diameter indicates the same type cells tend to cluster at short cell distance, which is likely due to proliferation and selective cell adhesion.   Table A in S2 Text.  Table A in S2 Text.   Table A in S2 Text.  Table A in S2 Text.   Table A