Interplay between epidermal stem cell dynamics and dermal deformations

We introduce a particle-based model of self-replicating cells on a deformable substrate composed of the dermis and the basement membrane and investigate the relationship between dermal deformations and stem cell pattering on it. We show that our model reproduces the formation of dermal papillae, protuberances directing from the dermis to the epidermis, and the preferential stem cell distributions on the tips of the dermal papillae, which the basic buckling mechanism fails to explain. We argue that cell-type-dependent adhesion strength of the cells to the basement membrane is crucial factors of these patterns.

Skin provides another example of such systems, where a sheet of proliferating cells (the basal layer) is attached via the basement membrane to a soft elastic substrate (the dermis). Its morphology has two distinct features: First, the interface between the basal layer and the dermis has many protuberances directing outward (namely from dermis to epidermis), which are called dermal papillae. Second, several experiments suggest that stem cells tend to be found on the tips of dermal papillae, while the transit amplifying cells occupy the rest of the dermal surface [17,18], which we have experimentally confirmed, as shown in Fig. 1. Although buckling mechanism is considered to be responsible for the formation of such dermal undulations [19], it alone does not explain why dermal papillae direct outward rather than inward, nor does it tell us why stem cells prefer the tips of dermal papillae.
Skin is an important organ that provides us with barriers such as protecting from damages and preventing dehydration [20,21], and its internal structure affects barrier functions: for example, diseases such as psoriasis are accompanied by altered structure of the epidermaldermal interface, where the germinative cell population increases and their activity is enhanced [22,23]. Al- though mathematical models have been proposed to understand skin barrier functions [24][25][26][27][28][29][30][31][32], the morphology of the epidermal-dermal interface and its relationship to stem cell patterning has been largely ignored, except for a speculation on the interplay between stem cell activities and the dermal structure [33].
In this work we introduce a mathematical model of cell division dynamics in the basal layer on a substrate composed of the basement membrane and the dermis, taking into account substrate deformability and cell type-dependent adhesion to the substrate. We employ particle-based modeling that has been used to study growing tissues [6,34] and the whole epidermal simulation [26,29,30]. By numerical experiments we show that our model can reproduce dermal papillae formation and preferential stem cell patterning on the tips of the dermal papillae. We argue that the outward deformation of the dermis is caused by the interplay of cell-membrane adhesion and membrane elasticity, and that the stem cell distribution is determined by the difference in adhesion strength between different types of cells.

II. THE MODEL
The inner structure of the skin is schematically shown in Fig. 2(a). The dermis and the basement membrane form an elastic substrate, and the basal layer cells are attached to the basement membrane. These cells repeat division on this substrate and, when differentiated, they migrate and join the suprabasal layer.
To investigate these cell dynamics and their possible effects on the substrate shape, we employ the following particle-based model, a schematic of which is shown in Fig. 2(b). Both the dermis and the basement membrane are approximated by different kinds of spheres, which we call dermal particles and membrane particles, respectively. Dermal particles, with radius R d , are considered as independent particles, whereas membrane particles, with radius R m , are connected by links to form a triangular lattice network. On this substrate we consider two kinds of cells, stem cells and transit amplifying (TA) cells, both represented as spheres with radius R c . We assume the following difference between the two: A stem cell divides into a stem cell and a TA cell, while a TA cell divides into two TA cells; stem cells can divide an infinite number of times, while TA cells can divide only a finite number of times; Stem cells are strongly bound to the basement membrane, while TA cells are weakly bound so that they are able to detach from the basement membrane when they undergo differentiation (conditions for differentiation is described below). For simplicity, we disregard the suprabasal layer in this model: differentiated cells are removed from the system.
Equations of motion for dermal particles, membrane particles, and the two kinds of cells are obtained by specifying interactions between them, together with the rules of cell division, which are explained below. defined by four network-connected membrane particles i, j, k, and l. (b) Cell division process. Cell i, initially atri, is dividing into two cells i and j. Division is completed when the natural length of the spring between the cells, which grows in time, exceeds 2Rc.

A. The dermis
A dermal particle can interact with (i) another dermal particle, (ii) a membrane particle, and (iii) the boundary wall, which is the lower boundary of the system. Let d be the distance between the two interacting entities. In the case (i), the interaction is strongly repulsive when d < Dδ 1 , weakly repulsive when Dδ 1 ≤ d < D, weakly adhesive when D ≤ d < ΛD, and zero when ΛD ≤ d, where D = 2R d is the contact distance of the two particles and the parameters δ 1 and Λ control the compressibility of the dermis and the cutoff range of the long-range interaction, respectively. This interaction is determined in such a way that the dermal particles as a whole exhibit plasticity and elasticity. In the case (ii), the interaction is the same as in the case (i), except that D = R d +R m and δ 1 = 1: there is no intermediate weak-repulsion range for the dermis-membrane interaction. In the case (iii), the interaction is repulsive only when the particle overlaps with the boundary and otherwise zero.
The total energy of dermal particles associated with these interactions is represented by U der , which is defined in Appendix B 1.

B. Membrane deformation
Elasticity of the basement membrane is introduced by assigning the stretching energy and the bending energy to individual links of the triangular membrane network. For a given link, the stretching energy is considered between the two particles with distance d connected by this link, which is harmonic in the range 2R m δ 2 < d < 2R m and steeply increases outside of this range, so that too much stretching and compressing are prevented.
The bending energy for a given link is considered between the two triangular regions sharing this link [see Here we assume that, as a function of the angle θ between the unit normals of the two triangular regions, the bending energy is proportional to θ 4 for small θ, which is softer than θ 2 as is usually used for elastic plates: This allows the substrate to exhibit local deformations, while preventing folding of the basement membrane for large deformations.
In addition, short-range repulsive interactions are required between membrane particles that are not networkconnected so that their overlapping is avoided when the membrane is largely deformed. The total energy associated with these interactions is represented by U mm , which is defined in Appendix B 2.

C. Cell-basement membrane adhesion and cell differentiation
Both stem cells and TA cells repulsively interact with membrane particles when they overlap. In addition, each cell adhesively interacts only with its nearest membrane particle. Note that, for a given cell, its nearest partner can change by the movement of the surrounding cells and membrane particles, so that cells can move along the surface of the basement membrane while keeping adhesiveness by changing their adhesive interaction partners. We assign different adhesive forces to stem cells and TA cells so that only TA cells can be detached from the basement membrane, whereas stem cells remain in the basal layer.
We identify differentiation of TA cells with mechanical detachment from the basement membrane: When a TA cell and its nearest membrane particle with distance d satisfies d > λ(R c + R m ), this cell is regarded as differentiated and removed from the system. Here λ is a parameter determining the differentiation threshold. We assume that non-proliferative TA cells, namely those which have lost the ability of cell division, lose adhesiveness to the basement membrane so that they are easily pushed away from the basal layer. Such a relationship between proliferation and adhesiveness has been experimentally observed [35,36]. The total energy of the cellbasement membrane interactions is represented by U mc , which is defined in Appendix B 3.

D. Cell division and cell-cell interactions
Cell division is modeled as a process of two initially overlapping cells being gradually separated [see Fig. 3(b)], which is different from our previous work [26] and similar to Ref. [6]. In this description, when cell i starts division, a new cell j is immediately introduced into the system, where the two cells are almost completely overlapped, with their initial positions given by r j =r i +ǫR c t and r i =r i −ǫR c t, wherer i is the position of i just before the division, ǫ a small parameter, and t a unit vector that is randomly oriented within a plane per-pendicular to the axis betweenr i and its nearest membrane particle, thus approximating a tangent vector of the basement membrane atr i . Let us call the two cells committing a division process a dividing pair. A dividing pair is connected with a spring whose natural length grows in time from the initial value 2ǫR c . The division is completed when the natural length reaches 2ǫR c ; At this point the spring interaction of the dividing pair is dissolved.
Interaction of two cells (with distance d) that do not form a dividing pair is repulsive when d < 2R c and adhesive when 2R c ≤ d < 2ΛR c . Note that each cell in a dividing pair behaves as an independent cell when it interacts with other cells or membrane particles. The total energy associated with cell division and cell-cell interactions is given by U cc defined in Appendix B 4. Modeling of cell cycle is the same as in the previous work [26] and is explained in Appendix C.

E. The equations of motion
The equations of motion are written down by using the above-introduced interaction energies as where Ω d , Ω m , and Ω c are the sets of dermal particles, membrane particles, and cells, respectively, and r i = (x i , y i , z i ) is the position of i in each set. Note that the number of TA cells in the basal layer changes in time due to cell division and differentiation, whereas the number of stem cells remains unchanged.
The model is considered in a three-dimensional region defined as 0 ≤ x < L x , 0 ≤ y < L y , and 0 ≤ z < L z , with periodic boundaries in the x and y directions. The boundary wall supporting the whole system is placed at z = 0. The system size is chosen as L x = L y = 50 and L z = 100. For the basement membrane, 200 × 230 membrane particles are arranged to form a triangular network on the plane z = 16R c . They are initially partially overlapped, with distance approximately equal to 2R m δ 2 . Dermal particles are placed between the basement membrane and the boundary wall by random packing. The radii are chosen as R c = 1.4, R m = 1.0, and R d = 1.4, so that neither cells or dermal particles can pass through the membrane network. The other parameters are set to δ 1 = 0.8, δ 2 = 0.25, Λ = 1.4, and λ = 1.4.

A. Formation of dermal papillae
First we perform numerical simulations of our model, starting with an initially flat basement membrane on which 64 stem cells are placed in a certain arrangement. Figure 4 shows a simulation result in the case of a ring arrangement of the stem cells. As time evolves, TA cells are continually produced first by stem cells and then by themselves, spreading over the surface of the basement membrane and soon entirely cover it to form the basal layer, while the basement membrane starts to deform in the regions where it is covered with TA cells. Finally many outward (upward in this case) protuberances are produced, which can be regarded as dermal papillae. During this process, stem cells also spread over the basement membrane. Noticeably, stem cells tend to be located on the tips of the dermal papillae. These observations show that our model reproduces basic features of the membrane shape and stem cell distributions observed in real skin.
We introduce the following statistical quantities to characterize these spatial patterns. First, the mean and Gaussian curvatures of the basement membrane at the position of cell i, denoted by κ H (r i ) and κ G (r i ), respectively, are evaluated from the membrane particle nearest to i and its link-connected membrane particles, using standard formulas of discrete differential geometry [37]. Then, to quantify the difference of these curvatures between stem cells and TA cells, we define the normalized curvature differences as follows: where Ω stem is the set of stem cells, N stem is the number of stem cells, and N c is the number of all basal layer cells (stem cells and TA cells). These quantities are the average membrane curvatures at the position of stem cells relative to all basal layer cells, non-dimensionalized by the cell radius R c . In addition, relative height of the stem cells (nondimensionalized by R c ) is defined as and the spatial inhomogeneity of the stem cell distribution in the x and y directions is defined as from which we measure σ = σ 2 x + σ 2 y . Figure 5 shows time evolution of these quantities for three different initial stem cell arrangements: the ring arrangement as shown in Fig. 4; a clustered arrangement, where all stem cells are gathered at the center; and a random arrangement. We find that the normalized curvature differences χ H and χ G and the normalized height difference ζ increase from zero and reach stationary values, which are significantly larger than zero, which means that stem cells tend to be located on the tips of the dermal papillae. Furthermore, χ H and χ G are faster than ζ to converge to the stationary states, indicating that small papillae are created at the position of stem cells first and then their sizes increase, not that the full-grown papillae are created first and then stem cells move toward their tips. Increase of χ H , χ G , and ζ starts at around t ∼ 300, at which the entire surface of the basement membrane is covered with TA cells, as shown in Fig. 4(b). This implies that dermal papillae formation is caused by the pressure created by continuous division of cells, which is weak when there is enough room for spreading of TA cells.
The quantities χ H , χ G , and ζ show no noticeable dependence on initial stem cell arrangements. In contrast, the spatial inhomogeneity measure σ shows this dependence, as shown in Fig. 5(d). Here the clustered arrangement, which is initially highly inhomogeneous, results in large σ values, whereas the random arrangement, which is initially already homogeneous, keeps small σ values, although in all three cases σ decreases in time and reaches stationary values at around t ∼ 300. This decrease in the early stage is because spreading of TA cells causes stem cells to separate from each other. In the later state, stem cells settle in the tips of the dermal papillae and resist further separation, resulting in the residual values of σ.

B. Numerical experiments by introducing non-dividing cells
The observed preference of stem cells to occupy the tips of the dermal papillae must be attributed to the difference between stem cells and TA cells, which in this model is either the way of adhesion to the basement membrane (strongly bound or detachable) or the ability of cell division (infinite or finite number of times). To explore this we introduce non-dividing (ND) cells, which have the same type of adhesion to the membrane as stem cells (strongly bound to and never detach from the basement membrane) but do not replicate at all; their dynamics are described by (3).
We place the same number of stem cells and ND cells (64 for each) on the basement membrane in four different initial arrangements, as shown in the top panels of Fig. 6. Results of numerical simulations at t = 1000 are shown in the bottom panels of Fig. 6. We find that both stem cells and ND cells prefer the tips of dermal papillae. Compared to stem cells, ND cells tend to share the same  (7)]. Results are shown for three initial stem cell arrangements: ring-shaped (corresponding to Fig. 4), clustered, or random.
papilla and sometimes form one large cluster as shown in Fig. 6(b).
In Fig. 7 we plot χ H , χ G , ζ, and σ for stem cells and ND cells: For the latter, Ωstem and N stem in these quantities are replaced by the set of ND cells and the total number of them, respectively. We find that significantly large positive values of χ H and χ G are observed for both stem cells and non-diving cells, indicating that these cells are located on the tips of dermal papillae. There is no noticeable difference in χ H and χ G between stem cells and ND cells, nor is there any difference between different initial arrangements.
Difference between stem cells and ND cells are captured in ζ and σ. The height difference ζ in ND cells is slightly smaller than that in stem cells in the early stage for all four initial arrangements, presumably because there are fewer cell division events in the vicinity of ND cells so that they feel smaller forces that can cause the membrane to locally deform. In the later stage, ζ in ND cells becomes larger than that in stem cells, especially in the case of Figs. 7(b) and (d). This tendency correlates with spatial inhomogeneity: The inhomogeneity measure σ in ND cells is significantly higher than that in stem cells. In particular, Figs. 7(b) and 7(d) show that ND cells have a stronger tendency to preserve initial inhomogeneity than stem cells, which reflects the observed clustering of ND cells in Figs. 6(b) and 6(d) (note the periodic boundaries). The correlation between ζ and σ, together with the slower increase of the former, indicates that clustering of ND cells enhances the formation of large outward protuberances.  Fig. 6.

IV. DISCUSSIONS
The results of the previous section indicate that the strength of adhesion to the basement membrane is crucial for the observed patterns. Now the following argument explains why protuberances of the basement membrane direct outwards. Let us consider a small membrane segment on which cells are attached. When a cell division occurs and the surface of the segment becomes crowded, either a cell has to detach from the segment or the segment has to deform so that it can accommodate all attached cells, depending on the energy cost for detachment and deformation. If deformation is energetically more preferable than detachment, the resulting shape of the segment will be convex so as to minimize the stretching energy. Thus from a local point of view every place has a tendency to create outward protuberances, resulting in dermal papillae.
This energetic argument also explains why stem cells prefer the tips of the dermal papillae: Cells feel the stronger pressure to push them out of the basal layer in the concave regions such as the bottom of the basement membrane than in the convex regions such as the top, which means that the larger cost is required to remain in the basal layer in the former case. In such situations the total cost can be minimized when strongly bound cells such as stem cells occupy the top and TA cells fill the remaining space.
A coarse-grained continuum description helps us understand this situation. For simplicity, let us consider a one-dimensional basement membrane, whose deviation (assumed to be small) from the flat surface is denoted by h(x, t), on which cells with different adhesion strength repeat division at a uniform rate. For a membrane segment in the region [x, x+∆x] with the length ∆s and the curvature κ, the adhesion energy between the cells and the segment is proportional to the length of the cell layer, which is estimated as ∆s(κ −1 + R c )/κ −1 , which reflects the fact that a convex shape can accommodate more cells than a concave shape of the same segment. Substituting ∆s = [1 + 1 2 (∂ x h) 2 ]∆x and κ = −∂ 2 x h at the lowest order of h and taking the limit ∆x → 0, we obtain the energy density associated with the adhesion, denoted by ǫ a : where K a (x) represents the difference in adhesion strength depending on cell types. Then the free energy functional can be written as F = (ǫ a + ǫ r )dx, where ǫ r represents the remaining interactions including elasticity of the substrate. By taking variation with regard to h, we arrive at where τ is a constant and R comes from ǫ r . Thus the cellsubstrate adhesion gives rise to the first two terms, both contributing to destabilizing the uniform state h = 0, which is balanced by the substrate elasticity included in R.
The first term is regarded as the inverse diffusion. This term is symmetric with respect to upward or downward deviations of h and its effect corresponds to buckling-like mechanism proposed in the previous works. In contrast, the second term is not symmetric in this sense and is responsible for both the formation of outward protuberances and the tip preference of stem cells: For example, when a strongly bound cell such as a stem cell is surrounded by weakly bound cells, in the continuum description K a (x) has a unimodal shape with a positive peak at the position of the strongly bound cell and K ′′ a < 0 holds near this peak. According to (9), this implies upward deviation of the substrate near a strongly bound cell. The same mechanism also works among TA cells when a cell division locally increases the adhesion cost at the position of the dividing TA cell, thereby creating the same unimodal shape of K a (x) there. This explains outward protuberances in the places where stem cells are absent, as observed in Fig. 4(a). Note that this second term appears only when we explicitly take into account the size of the cells, which vanishes when R c = 0.
It has been postulated [33] that the pressure created by reproducing cells is responsible for spatial patterning of stem cells. Surely this effect accounts for the spatial distancing of stem cells, clustering of ND cells by failing to create proliferating TA cells around them, and the accompanying deformation of the basement membrane between stem cells. However, it explains neither the growing direction of the dermal papillae nor preferential stem cell locations on their tips, which can be explained only when we take into account the differential adhesion strength between cells.

V. CONCLUDING REMARKS
We have shown that our model can reproduce basic features of stem cell distribution and the dermal structure, especially the formation of dermal papillae and preferential stem cell locations on their tips. We conclude that both the outward growth of dermal papillae and the tippreference of stem cells are the result of differential cell adhesion to the basement membrane. It has been recently reported that this tip-preference of stem cells could be experimentally reproduced on an elastic sheets with artificial papillae [38], which showed that no external signal from dermis was required for stem cell patterning, thus corroborating our adhesion-based patterning mechanism. Such a generic mechanism is expected to work in other pattern-forming systems.
Our results may shed light on photoaging, a physiological senescence process induced by ultraviolet exposure: Ultraviolet destroys elastic fibers in the dermis, leading to dermal stiffness. Photoaging causes changes in the dermal structures, such as flattening of the basement membrane and thinning of the epidermis [39], accompanied by the decrease in the number and the activity of stem cells [40]. This phenomenon is in accordance with our model, which predicts that when elasticity of the basement membrane is lost or the dermis is hardened by photoaging, undulations are suppressed and as a result the thickness of the whole epidermis is reduced because of the diminishing surface area of the basement membrane accommodating fewer TA cells. Numerical experiments of the whole epidermis with varying membrane stiffness would be possible by incorporating the present model into our previous epidermal model [26,27].
sections were stained without fixation by the following primary antibodies overnight at 4℃: anti-α6 integrin (BD Biosciences, San Diego, CA, USA; GoH3), anti-β1 integrin (Chemicon International, Billerica, MA, USA) and anti-anti-human type XVII collagen (TS39-3) [? ]. After washing in phosphate-buffered saline, the sections were incubated with secondary antibodies conjugated with Alexa488 for 1 hr at room temperature. All the stained sections were observed using a confocal microscope (Olympus Fluoview FV1000; Olympus, Tokyo, Japan). The institutional review board of Hokkaido University Graduate School of Medicine approved the procedures. Participants or their legal guardians gave their written informed consent.
Appendix B: Interaction potentials between the particles

The dermis
The total interaction energy of dermal particles with other dermal particles, membrane particles, and the boundary wall, is the following: where Ω m and Ω d are the sets of membrane particles and dermal particles, respectively, andr i is a mirror image of r i regarding the boundary wall at z = 0. The interaction functions u dm , u dd , and u ex are defined as where φ(x) = 1 12 Note that u dm = u dd for x ≥ 1: See Fig. 8(a) and (b). We choose K dd = 0.04, K dm = 0.04, δ 1 = 0.8, and α 1 = 0.125.

Membrane deformation
The total energy of the membrane deformation is given by where E m represents the set of links of the basement membrane network. The stretching energy is given by u str (x): where the constants A 0 = − α2(1−δ2 2 ) 2δ2 2 and A 1 = α2(1−δ2) δ2 2 are chosen so that the function is C 1 -continuous: See Fig. 8(c).
The bending energy is also considered in each network link: Link (i, j) uniquely determines a pair of two adjacent triangular regions, whose unit normals N Here the sign of the normals is chosen so that they have positive zcomponents at the initial membrane configuration. The interaction function is defined as We choose K str = 0.04, K bend = 0.5, δ 2 = 0.25, and α 2 = 1.25.

Cell-membrane adhesion
The energy between the membrane and the cells are given as follows: where Ω TA is the set of TA cells, Ω stem is the set of stem cells (including ND cells in Sec. III B), and nhd(i) ∈ Ω m is the nearest membrane particle to i. The adhesion function depends on cell type, where

Cell division and cell-cell interactions
The total cell-cell interaction energy is given by where Ω c is the set of all cells, Ω dp ⊂ Ω c × Ω c is the set of dividing pairs (two cells committing the same division process). The same adhesion function (B10) is used for cell-cell adhesion. For the interaction of a dividing pair, we adopt where t ij is the onset time of cell division between i and j: A pair of cells are stored in Ω dp when they start division and are removed when they complete division. Note that the natural length of the spring between a dividing pair is given by 2R c β 0 (t − t ij ).