A computational model of the epidermis with the deformable dermis and its application to skin diseases

The skin barrier is provided by the organized multi-layer structure of epidermal cells, which is dynamically maintained by a continuous supply of cells from the basal layer. The epidermal homeostasis can be disrupted by various skin diseases, which often cause morphological changes not only in the epidermis but in the dermis. We present a three-dimensional agent-based computational model of the epidermis that takes into account the deformability of the dermis. Our model can produce a stable epidermal structure with well-organized layers. We show that its stability depends on the cell supply rate from the basal layer. Modeling the morphological change of the dermis also enables us to investigate how the stiffness of the dermis affects the structure and barrier functions of the epidermis. Besides, we show that our model can simulate the formation of a corn (clavus) by assuming hyperproliferation and rapid differentiation. We also provide experimental data for human corn, which supports the model assumptions and the simulation result.

r r r i = −µ ∂ ∂ r r r i (U der +U md +U wall ) , i ∈ Ω der .
(1) d dt r r r i = −µ ∂ ∂ r r r i (U md +U mm +U mc +U wall ) , i ∈ Ω memb , Individual interaction functions are explained below.

The dermis
U der represents interactions between dermal particles, which has been modified from the previous work 29 : where ψ(x) = 1 12 x −12 − 1 6 x −6 .
The function ψ describes weak adhesion when two particles do not overlap (x > 1) and the excluded volume effect, i.e., the hard-core repulsion when two cells overlap (x < 1). The function u d (x) then describes short-range repulsive and long-range adhesive interactions between two dermal particles; when overlapping (x < 1), they first feel an elastic repulsion (1 − δ d /2R der ≤ x < 1) and then a hard-core repulsion (x < 1 − δ d /2R der ). For x ≥ 1, the two particles feel adhesive interaction with a cutoff range l * . The parameter δ d determines the compressibility of the dermis. The stiffness of the dermis can be controlled by varying the strength of the elastic interaction K sp . The functional form of u d (x) has been changed from the previous work 29 to incorporate this stiffness control.

The basement membrane
U mm represents the elastic interactions among the membrane particles, which are the same as the previous work 29 : where E m represents the set of links of the basement membrane network; u str and u bend describe the stretching and the bending energies assigned to a given link, respectively; The third term describes excluded volume interactions between two membrane particles not connected by a link, which becomes relevant to avoid overlapping in the case of large deformation. The functions u str and u bend are given by Here u str (x) is determined in such a way that the stretching elasticity is linear when the compression of the edge is small (δ m ≤ x ≤ 1 and nonlinear elasticity outside of this range, which prevents indefinite deformation by constant external forces. N N N (1) i j = (r r r i −r r r j )×(r r r (1) i j −r r r j ) (r r r i −r r r j )×(r r r i j are the two membrane particles linked to both r r r i and r r r j ; both of the normals are defined in such a way that they point outward (i.e., from dermis to epidermis).
U md represents interactions between a membrane particle and a dermal particle, which has been modified from the previous work 29 : The first term represents adhesion of dermal particle i to the nearest membrane particle nmp(i); δ is the Kronecker delta. The second term represents excluded volume interactions. The functions are given by where the adhesive interaction u md has a cutoff at r r r i − r r r j = l * memb (R i + R j ) and vanishes when the distance exceeds this range. The function u ex expresses only the repulsive interaction.

2/11
U wall represents the interaction between a dermal or membrane particle and a boundary wall placed beneath the dermis at z = 0 through the excluded volume interaction, as in the previous work 29 : where the position r r r i is a mirror image of r r r i regarding the boundary wall at z = 0.

Cell adhesion to the basement membrane
U mc represents interactions between basal cells and the basement membrane 29 : The first term describes excluded volume interactions between basal cells and membrane particles. The second and the third terms describe adhesion to the basement membrane for transit amplifying (TA) cells and stem cells, respectively, where χ stem i = 1 when cell i is a stem cell, and χ stem i = 0 for a transit amplifying (TA) cell; and The two functions u ad (x) and u bind (x) are defined in such a way that only TA cells can leave the basement membrane, whereas stem cells are bound to the basement membrane. For TA cell i, the adhesive force vanishes at r r r i − r r r nmp(i) = l * (R i + R memb ). At this point the TA cell is regarded as differentiated, and S i starts to increase (explained below).

Cell-cell interaction and cell division
U cc represents the interaction energy between two cells, which consists of two parts: where Ω c = Ω basal ∪Ω diff is the set of cells, and Ω dp is the set of a pair of cells in a cell division process. The first two terms in the bracket describe interactions between cells that are not involved with cell division, which are composed of excluded volume interactions and adhesive interactions. Reflecting its spheroidal shape, a differentiated cell i interacts with another cell j through the effective radii R eff i j and R eff ji , as defined below. The coefficient K ad is modified from the previous work 29 , and is given by the mean adhesion strength of two cells i and j: The differentiation-dependent adhesion strength κ i (t) is described below. Supplementary Figure S1 shows the functional form of the interaction potential between two cells i and j, K ex u ex (x) + κu ad (x) when two cells are of the same kind, κ i = κ j = κ. Note that, for x < 1, u ad (x) vanishes and the two curves reduce to u ex (x), the excluded volume interaction. The third term represents cell division. The cell division process is modeled in the same way as in the previous work 29 : Two dividing cells (i, j) are completely overlapped at the onset of division: they start to be separated by a spring force whose natural length l(t) increases in time from 0. Hence, u div and l(t) The division is complete when ∥r r r i (t) − r r r j (t)∥ > 1.8R cell .

Cell division cycle
We assume that stem cells divide infinite number of times and TA cells divide N div times. For both types of cells, a cell division cycle has two stages: in the deterministic stage, the cell division phase ϕ i (t) is governed by 29 where R(x) is the ramp function defined by R(x) = xΘ(x) and Θ is the Heaviside function. Starting from ϕ = 0, cell division occurs when ϕ (t) reaches 2π from 0; Without calcium ions, the period is given by T div . After that, the cell enters a stochastic process, where division occurs following the Poisson process with the division rate γ div . The average waiting time for division in the stochastic process is γ −1 div . We chose γ −1 div = T 0 /3, so that the overall cell division cycle on average is (4/3)T div . After the cell division, the phase is reset to ϕ = 0.

Cell differentiation
The degree of differentiation for cell i is expressed by the state variable S i (t): S i (t) = 0 for a basal cell, 0 < S i ≤ S SG for a spinous cell, S SG < S i ≤ S SC for a granular cell, and S SC < S i for a cornified cell.
Cell i is regarded as a basal cell when a nonzero adhesive force to the basement membrane exists, as explained above. When starting differentiation, S i (t) obeys where the differentiation speed is assumed to be accelerated by calcium ions c i (t) and a stimulantB(t). Dynamics of c i (t) are described below. Dynamics ofB(t) are governed by d dtB WithoutB(t), the model of S i is the same as the previous work 26 . Release ofB(t) is assumed to occur (χ(i) = 1) when cell i is a cornified cell that has at least one granular cell adhesively connected to it; otherwise, χ(i) is set to zero. Λ i is a set of all cells adhesive to cell i, through whichB(t) diffuses.

Cell deformation
When differentiated, cell i starts to flatten into an oblate spheroid with the flattening rate α i (t), with which the radius of the short and the long axes given by R which increases drastically when cell i becomes a granular cell (S i ≥ S SG ). For differentiated cells, the unit vector of the short axis n n n i (t) obeys 28 : n n n i (t + ∆t) = n n n i ∥ n n n i ∥ , where v v v i (t) = dr r r i (t)/dt. n n n i is normalized for each time step ∆t. The effective radius of cell i, when interacting with cell j, depends on n n n i and α i as where cos θ i j = ((r r r i − r r r j )/||r r r i − r r r j ||) · n n n i . Note thet R eff

Desquamation
In the basal, spinous, and granular layers, the adhesion strength κ i depends on differentiation state S i : In the cornified layer, we adopt the following desquamation model 49 . Cornified cells adhere to each other by corneodesmosomes. Corneodesmosomes are decomposed by degradation enzyme kallikreins (KLKs), which are expressed and released when cells differentiate from a granular cell to a cornified cell. The activity of KLKs is inhibited by binding with Kazal-type-related inhibitor (LEKTI), depending on the pH gradient in the cornified layer: binding is strong when pH is high as in the deep cornified layer and weak when pH is low as near the surface of the cornified layer. As increasing KLKs, corneodesmosomes decompose and cornified cells peel off from the epidermis.
Hence, in the cornified layer, we consider the following model for the adhesion strength κ i depending on the released KLKs, denoted by K Desquamation occurs when κ i ≤ K off and fewer than N conn cells are adhesive to cell i. The dynamics of K (r) i are described as follows. In the granular layer, we consider the produced and released KLKs, K where the constant K (p) max and L (p) max mean the maximum concentration of intercellular KLKs and LEKTI, respectively, and f prod and f rel describe the production and release of KLKs and LEKTI, respectively: Here the production and the release start to increase when a cell becomes a granular cell (S > S SG ). Normally, in the granular layer, calcium levels are small (c < c * ), so that the production is dominant; when c i exceeds c * due to calcium excitation occurs during cornification, the production stops and the release becomes dominant.
In the cornified layer, K , the concentration of the KLKs bound to LEKTI, are considered. The process of binding and decomposition of KLKs and LEKTI are described by is the set of cornified cells adhesive to cell i, and We have assumed that the diffusion of the KLKs-LEKTI complex is negligible and that, when a cell moves outward (as S i becomes larger) K i increase, due to the decrease of pH, which are described by f pH (x) and f pH (x). Therefore, in the deep cornified layer, the LEKTI-bound KLKs are dominant, whereas near the surface of the cornified layer, isolated KLKs become dominant, which decreases the adhesion strength κ i and induced desquamation.

Calcium dynamics
The model for calcium dynamics is the same as the previous work 26 . The dynamics of c i (t) (intercellular calcium), A(t, x x x) (the extracellular ATP concentration), P i (t) (the concentration of IP 3 ), h i (t) (the inactivation factor), w i j (t) (the activity of the gap-junctions between cell i and cell j), and B i (t) (the stimulant substance) is governed by

7/11
where G r r r, r r r i , Here ATP is released by calcium excitation described by the function G, and diffuses in the extra-cellular regions. IP 3 increase is caused by ATP as described in F P and diffuses through gap junctions, which appear when cells differentiate into granular cells. Calcium excitation is caused by the presence of IP 3 and a stimulant B and suppressed by the inactivation factor h i as in F c (P, c, h, B). The conductivity of the gap junctions w i j is determined by calcium concentrations. The stimulant B is released upon cornification (χ(i) = 1) and diffuses into adjacent cells. When calculating A(t, x), we impose periodic boundary conditions in the x and y directions and the Neumann boundary condition in the z direction.  26 . We compute the main dynamics with a constant calcium level. When 20 cells undergo cornification, we switch to the computation of calcium dynamics until reaching a steady state. Then we return to the main dynamics with a renewed calcium level. Parameter values used in the model equations are listed in Supplementary Table.1. Parameters that appear in the equations of motion for particles and cells are taken from Reference 29, and those in cell division cycle, cell differentiation, and calcium dynamics are taken from Reference 26. In these references, parameter values have been chosen so that numerical simulations could produce dermal protrusions with reasonable sizes and epidermal layers with well-defined boundaries. In this work, we have chosen smaller sizes of dermal particles and membrane particles, which resulted in the change of interaction strengths regarding the dermis and the basement membrane. To obtain numerical results consistent with the previous work, we modified parameter values in Reference 29. For the parameters that do not appear in these references, their values have been determined so that we could obtain reasonable numerical results. We should note that the parameter values in this work have not been directly verified by experiments.
We have chosen the non-dimensional unit of force by setting the mobility µ = 1, whereby the unit of other parameters is expressed by length [µm] and time [h].