Abstraction of the biological system

The microphysiological environment of the dermal-epidermal junction is the natural habitat of acquired nevi and the origin of most melanomas. The microanatomy of this area is governed by the characteristic shape of the basal membrane, a collagen type IV matrix interfacing both compartments [44, 45]. Normal melanocytes express CCN3, a matricellular protein that inhibits melanocyte proliferation and stimulates adhesion to collagen type IV. As a consequence, melanocytes are situated in the basal layer of the epidermis under physiologic conditions. In nevi, attachment of melanocytes to collagen IV is mediated through DDR1, which is under the control of CCN3 [46, 47]. Hence, we assume that the migration of nevus forming melanocytes in the basal layer can be described by mechanisms that let cell-agents glide along the membrane surface. Modeling the basement membrane in a differential geometric approach allows us to simplify the movement of nevus forming melanocytes as abstracted velocity vectors, while, at the same time, enables smooth and continuous migration along the complex geometric shape of the membrane.

Under physiologic conditions there exists a strong symbiosis between keratinocytes and melanocytes [16, 48, 49], which plays a crucial role in the proliferation and migration of melanocytes in the basal layer. In view of the feasibility and tractability of the individual-cell simulation approach, we integrate relevant interaction with keratinocytes in the abstracted behavior of melanocyte-agents. Thus, we simulate only nevus forming melanocytes but not keratinocytes or the background of normal melanocytes. Similarly, in our approach melanocyte-melanocyte interaction is not modeled on a biomolecular level [16] but abstracted in the movement dynamics of cell-agents (attraction, repulsion, collision). Hence, forces and collisions among cells in our model must not be understood as physical quantities (conservation of momentum, etc.), but as abstracted mechanisms that reproduce a larger set of complex processes [14, 15]. Detailed discussion of different configurations of cell motility are found in the following sections on the formation of the reticular and globular pattern.

Because our simulation approach is characterized by heavy abstraction of microscopic and cellular processes, also the parameters must be simplified and aggregate representations of biological quantities. For instance, we model the division of cells as discrete stochastic events that are parameterized with probabilistic likelihoods [30, 36]. In reality, cell division is a complex biological process that is not instantaneous but controlled by a large number of molecular reactions. To map the most crucial (extracellular) influences on the proliferation rate, we use heuristic damping factors that reproduce the biological effects we know or expect to exist in the biological system. This is, for instance, a reduction of the proliferation rate in areas with high cell density via contact inhibition [50] or induction of senescence of cells depending on the generation number (number of past divisions) due to telomere shorting [51].

Geometric model of the epidermal microanatomy

The spatial domain of simulated cell-agents is a cuboid section of virtual tissue that is separated into a dermal and epidermal compartment by a vaulted surface representing the basement membrane. We developed a detailed three-dimensional geometric model of the basement membrane to reproduce the typical microanatomy. For simplicity, skin appendages (hair follicles and eccrine glands) are not included. Dermal papillae are modeled as circular evaginations and their lateral S-shape is obtained from the revolution of parameterized shape curves. We formalize the basement membrane as a two-dimensional composite differentiable manifold.

In order to reproduce realistic tissue anatomy, we inferred statistical distributions of shape parameters and papilla density from measurements of histopathologic sections, and from in vivo images obtained by dermatoscopy and laser scanning confocal microscopy. Due to the absence of skin appendages, the number of dermal papillae per area is slightly greater in our model than in vivo. Histopathologic sections of excised tissue are distorted by shrinking [52], which has been taken into account in the parameterization of the geometric model (Fig 2).

Fig 2

Presentation of the three-dimensional geometric model of dermal papillae and the dermal-epidermal junction.

(A) Rendering of a computer generated basement membrane section of 1,000 μm × 1,000 μm. (B) Schematic lateral section of the geometric papilla model (location x, height H, radius R, scalar shape parameter p and derived shape curve S). (C, D) Routine hematoxylin and eosin stained histopathologic section of epidermis with dermo-epidermal junction measuring roughly 600 μm in width (C) and corresponding visual representation of computer generated tissue (D). (E, F, G) Horizontal views of the epidermal microenvironment as seen with dermatoscopy (E), in vivo confocal microscopy (F), and visualization of the geometric model (G), each measuring 500 μm × 500 μm. The lighter circular areas in (E) and (G) correspond to dermal papillae, which appear as hyporeflective areas surrounded by a rim of hyperreflective cells in confocal microscopy images (F). All three horizontal views show a square section with 500 μm side length. The horizontal line in (G) indicates the sectional plane in (D).

Agent-based model for migration and proliferation

We conceptualize melanocytes as three-dimensional spherical objects with internal states such as size and position that are subject to abstracted attractive, repulsive, and random movement forces. In concordance with the natural behavior and localization of melanocytes, the geometric model of the basement membrane is used to constrain the movement of simulated cells to the basal layer of the epidermis, and to prevent emigration from the epidermal compartment. Cellular processes such as division and differentiation are modeled as stochastic events. The temporal evolution of the simulated melanocyte population is obtained from an iterative algorithmic scheme.

To prevent two cells from occupying the same space, we use a technical mechanism that relaxes the force vectors of colliding spherical agents. Interaction with the manifold (basement membrane) is simulated by geometric transformations of velocity vectors and geodesic displacements. Fig 3A–3C presents a schematic of the movement model as a combination of movement forces, collision handling and geodesic translations. We introduce relatively unconstrained movement at a later point to simulate the scenario of nest formation. To modulate the movement forces and the likelihoods for cell division and differentiation on an individual basis, we use a measure of the local concentration of melanocytes and the cell generation number. Hence, the parameterization of the dynamic agent model consists of various stochastic quantities and control functions of the form

Fig 3

Schematic outline of the agent-based model.

(A) Abstracted intercellular forces (attraction and repulsion), random forces and external forces impose a velocity vector on each individual cell. (B) To minimize the overlap of spherical cells, a simplistic corrector algorithm is applied on the resulting velocity vector (for clarity shown in two dimensions). In all cases where the anticipated new position of a cell results in overlap with a neighboring cell, a certain (projected) component of the velocity, which we call the correction, is subtracted from the vector. Due to iterative application, this procedure ensures that the total overlap of simulated cells is minimized. (C) The corrected three dimensional velocity vector is projected into the coordinate space of the manifold model. By solving the geodesic differential equation, motion along the membrane surface can be simulated. (D) Outline of dynamic effects in the agent-based model. Large arrows (green and red) indicate effects that emerge naturally from the structural concept and the resulting internal logic of the model. Increase of the generation number of cells and of the local density result from cell division. Whereas attractive forces among melanocytes lead to higher density values, repulsion and random motion (diffusion) decrease the local concentration. Effects that do not emerge automatically, are controlled by functional influences (small arrows, blue).

Emergence of the reticular pattern

We show in a first experiment that our model can reproduce the connection of meso- and macroscopic traits in nevi with a reticular pattern. It is known that in vivo the reticular pattern reflects the microanatomy of the dermal-epidermal junction. A relative higher density of melanocytes and pigment on the lateral surfaces of dermal papillae produce ring shapes in horizontal dermatoscopic and microscopic images [6]. In our model we were able to reproduce this configuration with a small set of rules, which state that migration of cells is driven by random forces only. A selection of model parameterizations that produce reticular nevi with different features such as central hyperpigmentation, hypo-hyperpigmentation and homogeneous pigmentation as described by Hofmann-Wellenhof [53] are presented in Table B in S2 Text. Plots of the according functional damping factors as well as visual and quantitative simulation results are presented in Fig 4. A video showing the temporal evolution of a simulated reticular nevus is provided in S2 Video.

Fig 4

Parameter variation in simulated reticular nevi.

The results of four simulations runs with different parameter configurations are compared. Each scenario corresponds to a number (1-4) and a specific color (blue, orange, green, red). (A) Configuration of generation and density dependent damping factors. For each of the four scenarios, four parameter curves are displayed in the respective colors. The parameter curves describe the damping factors for proliferation and diffusivity depending on generation number and local density. (B, C, D) Visualization of the spatial cell generation profiles (B)—the color map refers to generation numbers between 0 and 60—, simulated dermatoscopic views (C) and histologic views (D). All images depict the state of the simulated lesions after 2,000 or 2,500 days, horizontal images display a region of 5,000 μm × 5,000 μm and white markers indicate the location of the histologic section. (E) Temporal evolution of the population size. (F) Temporal evolution of the average generation number in simulated melanocytes. Scenarios 2 and 4 use a tissue model with larger dermal papillae. An additional downwards force was applied in scenarios 2, 3 and 4. Different instances of growth arrest and corresponding radial generation profiles can be compared in scenario 1 and 3. The parameterizations of all four scenarios are listed in Table B in S2 Text.

We conclude that the reticular pattern results from a complex interplay of simple cellular dynamics (proliferation and diffusive migration) and the microanatomy of the dermo-epidermal junction (cell-matrix adhesion). A dynamic balance between proliferation and migration preserves a density equilibrium during growth. In low density areas, undamped proliferation increases the number of melanocytes, whereas in medium density areas, diffusive migration leads to the emigration of cell-agents. In high density areas, the packing of cell-agents impedes their displacement. Technically, the collision mechanism cancels out a large amount of the movement forces in individual cell-agents. The resulting freezing effect corresponds to increased melanocyte-melanocyte bonding [16]. However, to prevent unrealistic jittering of simulated cell-agents and to account for the presence of other cells that are not displayed by our model (e.g. keratinocytes), we explicitly dampen the initial movement forces in high density areas. In high density areas we also reduce the proliferation rate of melanocytes to mimic the biologic effect of contact inhibition [50, 54]. We further configure our model in such a way that proliferation and migration degrade with higher generation numbers, which corresponds to oncogene-induced senescence in vivo [55]. As a consequence, only a limited proliferative potential [29, 32] is available and the simulated nevus reaches a steady state and a final size as it is observed in real nevi. A certain proliferative potential in the confined center region cannot be exhausted because higher density prevents the proliferation of cells (compare the concentric density and generation profiles in Fig 4). Accordingly, the generation number is highest in the periphery of the nevus. In S3 Text we investigate mathematical models for the growth of simulated populations in abstracted non-spatial scenarios.

Different shapes and borders of simulated nevi, which mirror the variations seen in real nevi, are a consequence of stochastic events during simulation. However, the morphological appearance and quantitative characterization of simulated nevi is stable with respect to parameterization (S5 Text). That is, repeated simulation runs with the same parameterization produce adequately similar results.

In histologic sections of reticular nevi, melanocytes are preferentially situated at the base of the rete ridges and avoid the tips of dermal papillae. The reasons for this phenomenon are not completely resolved [56]. In our computer simulations, this pattern is automatically reproduced to a certain extent with above model configuration. This may indicate that this phenomenon emerges as a consequence of the microanatomy of the epidermis. A pronounced expression of this behavior in simulated melanocytes can be obtained by imposing a slight downward force on the cells. In this case, a higher density of adhesion molecules may be responsible for the preferential attachment of nevus forming melanocytes in this part of the epidermis. In S1 Text and S1 Video we investigate how the movement behavior of melanocytes and the particular geometric form of the basement membrane jointly produce this vertical inhomogeneity.

Nest formation and emergence of the globular pattern

The biological and molecular mechanisms that lead to the formation of nests are not fully resolved. However, it is known that the activation of certain adhesion molecules plays a pivotal role [15, 30, 31, 36]. Here, we interpret the emergence of globular nevi as the result of the occurrence of a nest-forming phenotype. In silico, commencing from the basic reticular configuration we introduce a certain cellular differentiation to simulate melanocytic nests as local strains of highly proliferative and slightly adhesive cells. This approach aligns with results from computational simulations that attribute the formation of melanocytic nests to proliferative activity rather than migration [30, 36].

Technically, differentiation into a globular strain occurs randomly during the division of melanocytes. The membership to a globular strain is inherited during subsequent cell divisions such that each nest stems from a single cell [29]. All members of a nest are allowed to detach from the basement membrane and the geodesic movement model (restriction to the basal layer) is replaced by unconstrained movement. As a consequence, adhesion among melanocytes in the same strain, which is simulated by local attractive forces, leads to spherical nests. Slight repulsive forces among cells, which are not in the same strain, prevent the fusion of different globules. External forces are used to prevent globules from entering the dermal compartment and also simulate cell-matrix adhesion in such a way that nests are predominantly located in the lower regions of the rete ridges [57]. Hence, in silico, the formation of nests is associated with a distinct differentiation and additional migratory dynamics whereas the behavior of single cells is analogous to the reticular scenario (compare the parameter configurations in Table B and C in S2 Text).

The size and temporal evolution of melanocytic nests was investigated in experiments in vitro and in silico [30, 36]. In our model the size and temporal evolution of nests is stochastically controlled through proliferation and emigration (cells leaving the nest). The likelihood of these events depends on the within-strain generation number of individual cells (number of divisions since the initial differentiation). Emigrated cells return to the behavior of the default population. If all cells of a globular strain emigrate, the nest is dissolved. In S4 Text and S3 Video according mathematical models for the size of simulated melanocytic nests are formalized and reproduced in simulations.

In Table C in S2 Text we present four slightly different parameterizations of this extended model. The parameterizations correspond to the four simulation scenarios in Fig 5 and demonstrate that a broad range of synthetic globular nevi can be generated by slight variation of a limited set of parameters. Furthermore, the formation, size and distribution of globules is similar to real nevi and repeated simulations with the same parameterization lead to quantitatively and visually similar results (S5 Text).

Fig 5

Parameter variation in simulated nevi with the globular pattern.

Results of four simulation runs with different parameter configurations are compared. Numbering and color coding is analogous to Fig 4. (A) Configuration of the emigration and proliferation factors for each parameter scenario. (B, C, D) Dermatoscopic visualization of the cell population (B) with enlarged region (C) and histologic views (D). All images depict the state of the simulated lesions after 1,000 days, horizontal images display a region of 5,000 μm × 5,000 μm and white markers indicate the position of the enlarged region and of the histologic section. (E) Temporal evolution of the total and nested population size. (F) Temporal evolution of the number of globules. Simulated dermal tissue with larger papillae was used in scenarios 1 and 3. A detailed overview on the parameterization is available in Table C in S2 Text.

If proliferation within nests is high (i.e. large within-strain generation numbers allowed) and emigration is low, globules are larger and the density of melanocytes in the nests is higher. An increased emigration factor allows melanocytes to steadily detach from nests and results in smaller globules and a larger population of single cells. The dermatoscopic and the histopathologic representations of scenarios 1 and 2 in Fig 5 reflect the expected behavior. If emigration does not decline with increasing generation number or exceeds a critical threshold with respect to proliferation, globules dissolve over time. In scenarios 3 and 4 in Fig 5 this behavior was reproduced. In all four scenarios, the nests reach their maximal size after approximately 100 days. Average nest sizes range between 200 and 500 cells.

In our simulations not only the size of nests but also the spatial distribution results from a balance between proliferation and emigration. If nests dissolve over time, the densely populated central part of the lesion, which has a low proliferative activity, is covered by a reticular pattern. The periphery, on the other hand, which has a relatively high proliferative activity, is typified by a globular pattern. This can be recognized in scenarios 3 and 4 in Fig 5 and corresponds to the active border of real growing nevi with the typical rim of globules in the periphery [5, 58] (see also Figs 1 and 6). As in real nevi, the peripheral globules indicate that the lesion is still growing.

Fig 6

Comparison of the growth dynamics of nevi as documented in clinical imagery and generated by our model.

(A, B, C, D) Comparison of the spatio-temporal evolution of the areal size of four simulated reticular (A) and globular (C) nevi as well as 13 reticular (B) nevi and 25 globular (D) nevi that were retrieved retrospectively from collected anonymized sets of digital dermatoscopic images. We used linear interpolation to quantify the growth dynamics in silico and in vivo. Black lines show the averaged linear growth dynamics. (E, F, G, H) Temporal evolution of the dermatoscopic exposition of real (E, G) and simulated (F, H) nevi with a peripheral rim of globules. Dissolving melanocytic nests in the center region yield a reticular pattern as seen in all four image series. In real nevi, additional to biological variations in the mesoscopic behavior of cells, tension forces such as relaxed-skin-tension, Blaschko and Langer lines can be the cause for morphological eccentricity and border morphology.

Temporal evolution and growth patterns

To further assess the plausibility of our simulation results, we compare the temporal evolution of simulated nevi with longitudinal clinical data. In particular, as a quantitative indicator, we correlate the evolution of the area occupied by simulated melanocyte populations with sequential area measurements of real nevi in retrospectively collected digital images of 25 globular nevi and 13 reticular nevi with more than 10 observations over time and an initial size of 30 mm² or less [5, 58–60]. The area of lesions in dermatoscopic images was obtained by image segmentation with a pretrained fully convolutional network where ResNet34 layers are reused as encoding layers of a U-Net style architecture [61] and subsequent conversion of pixels into square millimeters according to the resolution of the used videodermatoscope. The area of artificial lesions was measured from the positioning of all simulated melanocytes by calculating the area of the smallest enclosing disk and of the 0.9 quantile error disk. Due to the limited visible pigmentation effect of sparse outlying melanocytes, the 0.9 quantile error disk is a better approximation of the areal size as detected in dermatoscopic images. In Fig 6A–6D we demonstrate that the simulation of the growth dynamics of reticular and globular nevi is in accordance with biology. We further quantified the slopes of the growth curves by linear least-square interpolation. The obtained slopes for simulated and real reticular nevi were 2.271 × 10⁻³ mm² d⁻¹ (standard deviation 0.364 × 10⁻³) and 1.928 × 10⁻³ mm² d⁻¹ (standard deviation 1.889 × 10⁻³) and for simulated and real globular nevi 7.442 × 10⁻³ mm² d⁻¹ (standard deviation 0.876 × 10⁻³) and 6.005 × 10⁻³ mm² d⁻¹ (standard deviation 6.099 × 10⁻³). A unpaired t-test shows that globular nevi grow faster than reticular nevi in simulation (P = 0.398 × 10⁻³) and in vivo (P = 4.365 × 10⁻³). Futhermore, we see that for both patterns there is no significant difference between the linear growth approximations of real and simulated nevi (P = 0.546 in the reticular scenario and P = 0.277 in the globular scenario). The growth data is available in S2 Table.

In Fig 6E–6H we further qualitatively compare the mesoscopic and macroscopic dynamics of real nevi and simulated nevi composed of melanocytes that are able to form nests. When parameters are set in such a way that the probability of emigration (cells leaving the nest) is high, a rim of globules forms in the periphery and the reticular pattern is maintained in the center. This pattern is typical for growing nevi and can be reproduced in simulated lesions. If the probability for emigration is low, nests are maintained in the center. Videos showing the temporal evolution of a simulated globular nevus and a simulated nevus with peripheral globules are available in S4 and S5 Videos.

Differential geometric model of the basal membrane

A two-dimensional manifold is a topological space, locally homeomorphic to the Euclidean space $$ {\mathbb{R}}^2$$. That is, there exists a set of continuous charts with continuous inverse, that map surroundings (e.g. disks) in the manifold M onto surroundings in $$ {\mathbb{R}}^2$$. More specifically, there exists a covering collection of open sets {U_i}_i such that ⋃_iU_i = M and corresponding local charts φ_i: U_i → V_i where V_i are open sets in $$ {\mathbb{R}}^2$$. The collection {(U_i, φ_i)}_i is called an atlas. A manifold is k-times differentiable, if all chart-compositions $$ {\phi }_j^{-1}\circ {\phi }_i{|}_{U_i\cap U_j}$$ are k-times continuously differentiable. This compatibility property of the atlas is required for defining tangent spaces and vector fields on M. For a complete introduction to manifold theory, we refer to [65–67].

We start with a subset M₀ of the hyperplane $$ \left\{(x,y,z)\right|z=0\}\subset {\mathbb{R}}^3$$ representing a planar square section of the basal membrane without protruding dermal papillae. In combination with the projection mapping $$ \pi :{\mathbb{R}}^3\supset M_0\to {\mathbb{R}}^2:(x,y,0)\mapsto (x,y)$$, the flattened membrane is a smooth manifold with the atlas {(M₀, π)}.

In a next step, we regard dermal papillae as local evaginations from M₀ in $$ {\mathbb{R}}^3$$, formalized as charts (U_i, φ_i) where the disjoint open sets $$ U_i\subset {\mathbb{R}}^3$$ describe the surfaces of nonoverlapping papillae and the images φ_i(U_i) = V_i are disjoint open disks in $$ {\mathbb{R}}^2$$. The planar part of the basal membrane between dermal papillae (compare rete ridges) can be formalized by $$ N_0:=M_0\backslash {\bigcup }_i\overline{{\pi }^{-1}\left(V_i\right)}$$. The compatibility condition yields a differentiability requirement for the mappings φ_i and π at the interfaces between the papillae surfaces U_i and the planar membrane part N₀. With {(N₀, π)} ∪ {(U_i, φ_i)}_i we obtain an atlas for the (composite) manifold M₁: = N₀∪⋃_i U_i representing the basal membrane with protruding dermal papillae. Because π is the projection mapping, the representation of the basal membrane M₁ in $$ {\mathbb{R}}^2$$ appears distorted only within the base areas of papillae V_i.

In order to model coarser elevations and depressions, we additionally deform the membrane vertically by a mapping $$ {\chi }^{-1}:{\mathbb{R}}^3\to {\mathbb{R}}^3$$. The final manifold surface M₂: = χ⁻¹(M₁) is equipped with the atlas {(χ⁻¹(N₀), π ∘ χ)} ∪ {(χ⁻¹(U_i), φ_i ∘ χ)}_i and for convenience we introduce the notation $$ \psi :M_2\to {\mathbb{R}}^2$$ to refer to the collection of charts that identify locations on the basal membrane M₂ with locations in $$ {\mathbb{R}}^2$$.

A key feature of the formalization by differentiable manifolds is the identification of tangential vectors on the vaulted basal membrane with tangential vectors in $$ {\mathbb{R}}^2$$ by the linear mapping

Shape model for dermal papillae

We model each dermal papilla as a surface of revolution which is obtained by rotating a one-dimensional curve around a vertical axis in $$ {\mathbb{R}}^3$$. The parameterization of a dermal papilla is composed of height $$ H\in {\mathbb{R}}_{+}$$, base radius $$ R\in {\mathbb{R}}_{+}$$ and the shape curve S. According to the differentiability requirement, the surface of each papilla U_i must transition smoothly into the planar part N₀ of the basement membrane. To meet this requirement and to be able to reproduce the characteristic lateral S-shape of dermal papillae, we model shape curves as three-times continuously differentiable paths,

With the parameter tuple (x_i, R_i, H_i, r_i, h_i), where R_i is the radius of a disk V_i with center $$ x_i\in {\mathbb{R}}^2$$, each dermal papilla can be represented as a chart $$ {\phi }_i^{-1}:V_i\to U_i$$, with local polar representation

Fig 7

B-spline model for the shape of dermal papillae.

(A) Lateral section of an ensemble of randomly generated papillae (green lines). (B) Positioning of control vertices (red markers connected by lines) for a particular shape curve (black). The positions of the fourth and fifth vertex along the r-axis can be altered with the shape parameter p. (C, D) Radial and vertical component of the shape curve (blue) with corresponding derivatives (red) and auxiliary lines. Compare the differentiability requirements in Eq (2).

The transformation χ⁻¹ is modeled as a superposition of larger circular elevations and depressions. For the vertical components the function model for h can be reused.

Parameterization of the geometric model

We measured the dimensions of 100 dermal papillae from routine histopathologic sections of a retrospectively obtained convenience sample of 13 biopsies. Measured dimensions include the height $$ \overline{H}$$, radius (half of the diameter) at the base $$ {\overline{r}}_{h=0}$$, at half of the height $$ {\overline{r}}_{h=0.5H}$$ and at 10 μm below the tip $$ {\overline{r}}_{h=H-10}$$. The average measured height was 59.7 μm with a standard deviation of 18.2 and the mean base radius was 28.0 μm with a standard deviation of 10.1. The raw measurement data is provided in S1 Table. We accounted for dehydration by increasing all lengths by 20% [52] but still assumed a certain underestimation and distortion in the measurements due to inhomogeneous shrinking and deformation. Additionally, the histological section technique implies that papillae are usually not cut at their thickest position and that the sectional plane does not perfectly fit the vertical alignment of papillae.

We used optimization techniques to fit the mathematical model for shape curves (R, H, p) to each recorded papilla $$ (\overline{H},{\overline{r}}_{h=0},{\overline{r}}_{h=0.5H},{\overline{r}}_{h=H-10})$$ individually. Neglecting any statistical correlations, the solution ensemble (parameters R and p) as well as additional measurements were fitted with beta and log-normal distributions. Based on the resulting statistical model we were able to randomly generate dermal papillae and basal membranes.

From visual inspection and comparison with horizontal reflectance confocal microscopic images taken in vivo (compare Fig 2F and 2G), we came to the conclusion that simulation of realistic tissue anatomy requires approximately twice the radius (and hight) as we extracted from histopathologic measurements. Hence, artificial papillae are sampled with a base radius between 56 and 75 μm (standard deviation between 8.5 and 11) and a height between 70 and 135 μm (standard deviation between 13 and 21) depending on the simulated characteristics of skin tissue. We also recognized that the shape of measured papillae is equally distributed between the broad-shouldered and cylindrical types. In the parameterization of the geometric model we intentionally introduce an over-representation of the shouldered type (S-shape) in order to give enough space for the formation of melanocytic nests between (non-deformable) dermal papillae.

To validate the number of dermal papillae per skin area, we counted dermal papillae on 40 1 mm x 1 mm sectors of the inner lower arm of a 35-year male volunteer via horizontal in vivo reflectance confocal microscopy (VivaScope 3000; MAVIG GmbH, Munich, Germany; raw data in S1 Table). We found a mean of 48.8 papillae per mm² with standard deviation of 5.29 corresponding to values reported in literature [68–71]. Due to the absence of additional skin features like eccrine gland ducts and hair follicles in silico, the density of papillae in the computer model must be higher in order to obtain realistic spacing (corresponding to the dimensions of rete ridges). Also the perfect circular contour of artificial papillae introduces artifacts in the perception of the density. We came to the conclusion that the optimal result is obtained with maximum packing density, corresponding to a value between 50 and 160 papillae per mm² depending on the typical size of papillae and the simulated tissue characteristics. The model parameters obtained with this approach are presented in Table A in S2 Text.

From the papilla shape model (Fig 7) it follows that the largest diameter of dermal papillae is always attained at the base. As a consequence, the problem of generating dense spatial layouts of non-intersecting papilla bodies can be reduced to placing non-overlapping disks on a plane. To generate spatial layouts with maximum packing density, we align papillae on a hexagonal grid and use iterative local optimization of their center positions to account for the differences in disk radii.

Melanocyte population dynamics

We developed an iterative scheme that advances the melanocyte-agents in constant time steps Δt. A typical value for Δt is 1 d, but all algorithms and calculations were designed to scale correctly with different time increments (see S6 Video and also S3 and S4 Text). To reduce the effect of the initial conditions and the warm-up phase (high stochastic volatility due to small population) on the overall long-term result, simulations are usually initialized with a group of about 10 melanocytes located at the center of the basal membrane segment.

According to our conceptual model, we assume that certain intra- and extracellular conditions can impact the rate of proliferation [16, 54]. Hence, individual cell behavior is regulated by two agent-centric state values of which the first is the generation number and the second is the local density. The generation number g is an integer value that counts the number of cell divisions a particular cell has undergone since the start of a simulation run. We define the local density as the distance-weighted neighborhood measure

Overall, we expect exponential growth dynamics in the number of melanocytes $$ M^{\prime }\left(t\right)=\tilde{p}M\left(t\right)$$ with the proliferation rate $$ \tilde{p}$$ degrading with generation and density such that $$ \tilde{p}=\tilde{p}(g,\rho )=A\left(g\right)B\left(\rho \right)p_0$$ where A and B are polynomial functions with values in the interval [0, 1] and p₀ is the base proliferation rate. For instance, to limit the maximum allowed generation to 60, we let A(g)→0 with g → 60. The density-dependent damping factor B(ρ) is used to reduce population growth in already congested (high-density) areas. In the quantification of growth rates of melanocytic lesions in vivo different measures are possible [5, 72]. The deduction of cell-level proliferation rates from macroscopic data is difficult and we assume that division rates can vary greatly depending on environment, cell-type, and differentiation [33]. Experimental studies in vitro indicate values between 0.02 h⁻¹ and 0.04 h⁻¹ for melanoma cells [30, 36]. We heuristically set the base proliferation rate p₀ to 0.05 d⁻¹. To simulate elevated proliferation within melanocytic nests we use a proliferation rate p₀ of 0.1 d⁻¹. From the macroscopic exponential growth model we calculate the individual per-step cell-division probability

Our model neglects programmed cell death (apoptosis), only cells leaving the simulation domain are discarded. Cell division is assumed to be symmetric; if a division occurs in a cell, the generation number of both filial cells is increased by one. Various individual attributes of the parent such as size, location but also the base proliferation rate p₀ are inherited. This enables the simulation of ancestral strains of melanocytes with specific genetic expressions [13, 17, 29] (e.g. for nested melanocytes p₀ = 0.1). We distort the inherited values with additive Gaussian noise (e.g. σ = 0.01 for the base proliferation rate). During division of normal cells, differentiation into the nesting type occurs with a (conditional) probability q₀, which is set to 0 in scenarios with only the reticular pattern and 0.1 in globular scenarios. Nested cells are emitted with a base rate s₀ of 0.1 d⁻¹.

An overview on the parameterization of cell proliferation is presented in Table A in S2 Text. Visual outlines of the algorithms and routines as well as further implementation details are presented in S6 Text. Possible configurations of the generation dependent damping A(g) and an effects analysis in the aggregated exponential growth model can be found in S3 Text. A comparison of different within-nest proliferation dynamics can be found in S4 Text and S3 Video.

Cell migration model

The following iterative algorithm is used to simulate the migration of cells. For each cell-agent a new position at time t + Δt is calculated from the current position taking into account the cells generation, the local density, interactions and collisions with other cells and the physiological constraint. A visual schematic of the movement model was presented in Fig 3. Here, we present in detail the mathematical formulation of this model. An outline of the corresponding algorithm and implemented code is provided in S6 Text.

At the beginning of each movement phase, vectors pointing to neighboring cells n and the respective distances are calculated by

An initial velocity vector v⁽ⁱ⁾ is composed of external forces v_ext, such as a downwards moving trend, and diffusive noise. The noise term is a trivariate normal random variable scaled by the base diffusivity D and damping factors Q and R, such that

For each neighboring cell n, attractive and repulsive forces can modulate the initial velocity. Depending on the distance to the neighbor, on local cell density and on the generation number, the functions F_a and F_r in combination with the functional damping factors G_a,r and H_a,r determine the orientation and intensity of pairwise intercellular force vectors (compare Fig 3A).

The following corrector mechanism implements collisions between cell-agents. If the predicted velocity vector either leads to a collision with a neighbor n or if prevailing overlap will not be resolved automatically ($$ d_n<\frac{v_n}{\Vert v_n\Vert }·v^{\left(ii\right)}\Delta t$$), a fraction K of the respective projected velocity component is subtracted from the predicted vector (compare Fig 3B).

In order to confine cells to the basal layer (we assume the cell is currently located on the membrane surface, compare Fig 3C), the velocity vector v⁽ⁱⁱⁱ⁾ is projected (or rotated) into the tangential plane of the manifold M₂ at location x^(t). The resulting tangential vector v^(iv) is then transformed into a tangential vector $$ \overline{v}:=d{\psi }_{x^{\left(t\right)}}·v^{\left(iv\right)}$$ of the coordinate space $$ {\mathbb{R}}^2$$ at location $$ {\overline{x}}^{\left(t\right)}:=\psi \left(x^{\left(t\right)}\right)$$ using the linear mapping in Eq (1). By solving the geodesic differential equation with initial conditions $$ ({\overline{x}}^{\left(t\right)},\overline{v})$$ on the interval [t, t + Δt], a new parameter location $$ {\overline{x}}^{(t+\Delta t)}\in {\mathbb{R}}^2$$ is obtained. Finally, the new position on the basal membrane M₂ is the inverse transform $$ x^{(t+\Delta t)}={\psi }^{-1}\left({\overline{x}}^{(t+\Delta t)}\right)$$. Further details on the application of the geodesic differential equation can be found in S1 Text.

For cells that are part of a nest we skip the last step and calculate the new position as x^(t+Δt) = x^(t) + Δt(v⁽ⁱⁱⁱ⁾ + v_derm) where v_derm is non-zero only if the cell is located in the dermis (below the membrane). For simplicity v_derm is either a vertical vector if the cell is below the base level of the membrane or a horizontal vector if the cell is located within a dermal papilla.

A parameter overview is presented in Table A in S2 Text. For simulations we always assume the default case with absent intercellular forces F_a,r ≡ 0 and constant damping factors (≡ 1).

Visualization techniques

Dermatoscopic and histopathologic diagnosis methods are the current state-of-the-art for evaluation of melanocytic lesions. We developed rendering techniques that visualize the physiology and global status of simulated cell populations in ways mimicking the image styles of both microscopic techniques. Dermatoscopy, as a non-invasive screening technique [41, 42], enables physicians through cross-polarization or immersion fluid to have an en-face view of pigmented lesions on the dermo-epidermal junction. However, melanocytes are not directly visible, but morphological features of melanocytic lesions are indicated by the local pigmentation of epidermal tissue. To produce comparable visual representations, we calculate two-dimensional horizontally resolved histograms from the positioning of simulated melanocytes and apply Gaussian filters on the normalized data in order to approximate tissue pigmentation. As a consequence, all melanocytes are assumed to eradiate pigmentation in a radially symmetric fashion irrespective of microanatomic boundary conditions and depth. To obtain more realistic visual representations, we plan to use three-dimensional rendering techniques, where effects like the depth-dependent translucency of epidermal tissue, normalization of color-space and the basal membrane are taken into account. Histopathological analysis, on the other hand, is the current gold-standard for a final and confirmative diagnosis of skin tumors, albeit evaluation is not entirely independent from clinical information [74]. In typical hematoxylin and eosin stained biopsy samples, pigmentation can be recognized to occur mainly in keratinocytes; melanocytes are usually typified by a white halo around the nucleus. In the computer generated histologic image style we represent sliced melanocytes as brown disks.