Patterns in melanocytic lesions: impact of the geometry on growth and transport inside the epidermis

In glabrous skin, nevi and melanomas exhibit pigmented stripes during clinical dermoscopic examination. They find their origin in the basal layer geometry which periodically exhibits ridges, alternatively large (limiting ridges) and thin (intermediate ridges). However, nevus and melanoma lesions differ by the localization of the pigmented stripes along furrows or ridges of the epidermis surface. Here, we propose a biomechanical model of avascular tumour growth which takes into account this specific geometry in the epidermis where both kinds of lesions first appear. Simulations show a periodic distribution of tumour cells inside the lesion, with a global contour stretched out along the ridges. In order to be as close as possible to clinical observations, we also consider the melanin transport by the keratinocytes. Our simulations show that reasonable assumptions on melanocytic cell repartition in the ridges favour the limiting ridges of the basal compared with the intermediate ones in agreement with nevus observations but not really with melanomas. It raises the question of cell aggregation and repartition of melanocytic cells in acral melanomas and requires further biological studies of these cells in situ.


Introduction
Melanomas are the most deadly skin cancer, being responsible for 75% of the mortality due to skin cancer according to the American Cancer Society. 1 Unlike cancers affecting other organs, these tumours are directly observable since the primary tumour appears as a pigmented lesion at the surface of the skin. Early detection is therefore made possible by simple skin examination, eventually performed by the patient himself. When a melanoma is detected at an early stage, it can be treated by simple excision and the 10 year survival rate is higher than 99%. However, the survival rate drops to less than 50% when it penetrates deeply into the dermis. In the last decades, many efforts have been made to improve the methods of differential diagnosis in order to classify malignant and benign melanocytic lesions based on morphological criteria. Empirical studies from collections of clinical cases have led to the identification of shapes and microstructures, especially thanks to more and more precise tools like dermoscopy [1,2]. But the mechanisms generating these structures as the morphological differences between malignant and benign tumours remain largely unknown. We recently proposed physical mechanisms controlling the contour regularity of melanocytic tumours [3][4][5] and explaining the apparition of microstructures such as pigmented dots and globules [6]. However, the models that have been developed until now suppose that the various layers of the skin are horizontal with the avascular growth occurring in a thin epidermis with a simple flat geometry. In certain regions of the body, the skin has a more pronounced geometry, often associated with specific microstructures, such as dermal papillae. We discuss here how the skin geometry can influence the pattern formation in melanocytic tumours. The effect of the geometry on growth has already been demonstrated for morphogenesis experimentally [7] and theoretically [8].
Human skin can be divided into three layers, epidermis, dermis and hypodermis. The epidermis is the superficial layer of the skin and is mainly composed of keratinocytes and to a less extent of melanocytes (figure 1). Keratinocytes proliferate in a basal monolayer attached to the basement membrane forming the dermal-epidermal junction, and migrate towards the skin surface during their differentiation. The stratum corneum is the outermost part of the epidermis and is made of fully differentiated keratinocytes and of non-living corneocytes, in a lipid-rich matrix regulating skin permeability. In healthy tissue, each melanocyte remains connected to neighbour keratinocytes and to the basement membrane. Its main role consists of producing the pigment of the skin called melanin. Melanin is enclosed in vesicles and is then transported by neighbour keratinocytes via endocytosis and exocytosis. Melanocytic lesions such as nevi and melanoma originate from a dysregulation of melanocytes leading to the invasion of the surrounding tissue. The structure of human skin varies significantly depending on the location on the body. We can define two main types of skin: non-glabrous and glabrous skin. Non-glabrous skin is characterized by the presence of hair follicles and a thin epidermis (100 mm thick typically) almost planar except near the hair follicles. Hairless skin, which is found on the palms and soles, is instead characterized by a thick epidermis (1 mm thick typically), including a dense and thick stratum corneum. At the skin surface, dermatoglyphs are formed by a periodic alternation of ridges and furrows [9]. At the dermal-epidermal junction, an epidermal ridge lies under each surface furrow and each surface ridge: the limiting ridges (crista profunda limitans) and intermediate ridges (crista profunda intermedia), respectively. Intermediate ridges are narrower than the limiting ridges [10] and host the eccrine sweat glands (acrosyringium) that appear by dermoscopy as white dots at the centre of each surface ridge [11]. This complex structure plays an important role in the perception of touch and induces a strong mechanical coupling between the dermis and epidermis in palms and soles allowing the skin to support high mechanical stresses [12]. Among non-white populations, glabrous skin is the most common location for melanoma. For instance, in Japanese population 50% of detected melanomas are found in these locations [13]. Clinical research has shown that nevi and melanomas developing on glabrous skin have specific shapes, certainly influenced by the geometry of the skin in these locations, and requiring appropriate diagnostic criteria.
On glabrous skin, nevi and melanomas in situ are generally associated with parallel pigmented stripes. For nevi, these stripes are usually located along the furrows of the skin surface, parallel furrow patterns (40-45% of acral nevi [14]), unlike melanomas which are associated with similar pigmented parallel patterns but located along the ridges of the skin surface, parallel ridge patterns (83% of melanoma in situ) [14][15][16][17]. Saida et al. [14] reported a sensitivity of 86% and a specificity of 99% for the early diagnosis of acral melanomas using an algorithm based on these parallel patterns. These shapes suggest a strong influence of the geometry of the epidermis on the distribution of melanocytes and melanin in agreement with clinical observations [10,18]. For instance, using electron microscopy, Nagashima & Tsuchida [10] have observed the geometry of the dermal-epidermal junction at different sites of the feet and have shown that on each site, the pigmentation pattern of nevi follows the structure of epidermal ridges.
A tissue section of acral nevus performed perpendicularly to skin fingerprints frequently shows the presence of melanin  columns extending from the limiting ridge of the dermalepidermal junction towards the skin surface [19]. A preferential proliferation of nevus cells in the limiting ridges has been proposed to explain the presence of these columns. However, recent investigations by Palleschi et al.
[20] and by Saida et al. [17] show that this explanation is not sufficient. In some acral nevi, these authors show that aggregates of nevus cells are present in both limiting and intermediate ridges, but columns of melanin are found only above limiting ridges. Inhibition of melanin synthesis in intermediate ridges could explain this phenomenon.
In this paper, we aim to model such clinical observations. As dermascopes detect the melanin pigment and not the melanocytes, we consider the proliferation of tumour cells on a distorted basement membrane and also the transport of the originated melanin in the same geometry. For this purpose, we develop two mathematical models. Both models take into account the existence of epidermal ridges. These models provide a geometrical explanation for the apparition of parallel patterns in melanocytic lesions of glabrous skin and for the localization of melanin columns above limiting ridges. In §2, we revisit the avascular growth model for melanomas [3][4][5][6] taking into account the geometry of the epidermis. In §3, we propose a new model for melanin transport and melanin column distribution inside the epidermis. In §4, we discuss the physical results which correctly describe the nevi stripes, but not completely the pigmented stripes of acral melanoma. In appendix A, we give a more detailed description of the methods used in this study, a list of symbols is given in table 1 and a list of parameter values is given in table 2.

Model of melanoma cell proliferation on curved surfaces
In order to understand the clinical observations of melanocytic lesions on glabrous skin, we first focus on cell proliferation and tumour growth. Several approaches have been developed to model tumour growth [21][22][23]. We develop a continuous model based on the theory of mixtures [24,25]. Even though there have been studies of diffusion on a curved surface [26] and proliferation of curved epithelium [27], no cancer models take into account the curvature of the epithelium. For simplicity sake, this mixture model takes into account only two phases [4,6]: a cancerous proliferative phase with a concentration f c and a second healthy phase containing the interstitial fluid, the dead cells and the keratinocytes with a concentration f h ¼ 1 2 f c . The concentration represents the percentage of each kind of cell at a given point inside the lesion, but an average is achieved on a scale larger than the cell scale. Both phases need nutrients to ensure cell functions. In this section, we present briefly the chosen geometry of glabrous skin and we detail the governing equations induced by the curvature effects of the limiting basement membrane. These equations can be treated only numerically and the simulations will be compared to clinical observations.

Geometry and governing equations
We consider the growth of a thin layer with thickness e on a curved surface S B . S B represents the basement membrane, and the layer is constituted of the proliferative zone containing melanocytes, so basically the spinous and basal layer of the epidermis. Its order of magnitude is of a few cells size, so few times 6 mm. Therefore, we can consider the layer as two dimensional and the fraction of cancerous phase f c (x, y) at the point of coordinates (x, y), once averaged does not depend on z. In the following, we use the Monge representation of the surface S B , which means that the vertical coordinate of one point of the surface is z B ¼ h(x, y) (figure 2). The classical equations of tumour growth [3][4][5] are modified by this geometry. The induced modifications are rather technical requiring the tools of differential geometry. It turns out that it has never been treated before in the literature of tumour modelling perhaps because of its specificity to skin cancers. A complete description of curvature effects on differential mathematical operators is given in appendix A. Inside the thin layer, the small nutrient molecules of concentration n diffuse and are consumed by the cells. The timescale for diffusion being much shorter than the uptake time, the diffusion is treated at equilibrium, and we derive the following three-dimensional equation: where D is the diffusion matrix, with coefficients D ¼ D jj when the diffusion occurs along the surface S B , and D ¼ D ? when the diffusion occurs perpendicularly. d n (resp. k) is the nutrient consumption rate by the cancerous phase (resp. by the healthy phase). The thickness e being small, an averaged description of all quantities is enough for our purpose restoring a twodimensional description. It is why we integrate this equation along the z-axis giving where N(x, y) is the surface concentration at (x, y); G, g and H are metric tensors given in appendix A. D 0 is the Laplace-Beltrami operator (see equation (A 15)); n 1 and n 2 are typical values of the nutrient concentration outside the thin layer (table 2). The first term of equation (2.3) represents the diffusion along the curved surface, whereas the second and third terms represent the orthogonal diffusion. Note that the third term due to the curvature will change sign along the surface. Although small, it is positive for ridges and increases the orthogonal flux Figure 2. Monge representation of the surface S B and representation of the thin layer of thickness e. A point on the surface S B is represented by , and the normal vector to the surface S B at x(x, y) is written as n(x, y). The thin layer is represented by the set of rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20140339 and the local nutrient concentration as well. The two last terms represent nutrient consumption. Compared to previous works, the novelty of such equations concerns the real treatment of the basal geometry. Note that the averaging process allows recovery of a two-dimensional tractable equation, at least for simulations, without loss of pertinent information.
As we have assumed the same mass density for the two phases f c and f h , the mass exchanges between them are the same, and we can focus on the mass balance equation for f c where v c is the velocity of the cancerous phase, and G(f c , N) is the mass exchange rate. Using a variationnal principle, v c is related to the derivative of the free-energy between cancer cells Df c , C(f c ) being the free-energy between melanocytic cells. As expected for a cell-cell interaction, the volumic contribution f(f c ) in P is weak at low concentration, attractive at intermediate concentration and becomes repulsive. Thus, it has a simple representation reminiscent of the Lennard-Jones potential (figure 3). A surface term e 2 c Df c is added to P (f c ) penalizing large concentration gradients.
For the mass exchange rate G(f c , N), a simple linear law is assumed: G(f c , N) ¼ g c f c (N/n s 2 d) with death and proliferation rate, d and g c , being constant. Therefore, the final twodimensional equation for the cell concentration is given by with the Laplace -Beltrami operator (appendix A) noted with a superscript 0 given by þ @ @y f c K(1 À f c ) 2 g 1=2 g yy @S @y þ g yx @S @x : (2:7) We rewrite equations (2.3) and (2.6) with dimensionless quantities: x ¼ x/l n and l n ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi D k /d n p being the nutrient penetration length,t ¼ tg c n s andN ¼ N/n s . Finally, dropping hats and adapting all coefficients accordingly, we get the two coupled master equations to solve The first mathematical symbol of equation (2.8) represents the lateral diffusion on a curved surface (equation (A 15)). The second (resp. the third) term accounts for the nutrient consumption of the melanocytic cells (resp. the healthy cells, mostly keratinocytes). The fourth term is the main term of the transversal diffusion, and the fifth term is the curvature contribution to the transversal diffusion. In equation (2.9), we recover the terms of the mass balance equation (equation (2.4)), the velocity in the convection term is calculated on a curved surface (equation (2.7)) and the third term expresses the cell proliferation and death. Only numerical simulations will allow recovery of some features of clinical observations. Note that table 2, which gives biomechanical quantities measured in the literature, makes the model quantitative as well.

Numerical results
Numerical simulations of equations (2.8) and (2.9) were performed in a two-dimensional N Â N lattice, with f c ; 0 and n ; n 0 at the lattice border. The initial conditions are chosen to follow parabolic distribution for the concentration with addition of an initial white noise. For the simulations, we choose f in P as where p 1 1, in order to represent the form of figure 3. Here, we fix p ¼ 2 and f e ¼ 0.6. The simulations are executed in the following steps: -calculate the value of P on the grid and -implement the cell density (resp. the nutrient) value using a forward time centred space scheme to represent equation (2.9) (resp. equation (2.8)). A time dependence with a different timescale has been introduced in equation (2.8) to ensure that growth is slower than diffusion.
The program was implemented on a graphic processor unit, a GTX-580 NVIDIA graphic card, using the CUDA parallel programming language, in order to reduce computational time (appendix A.3.1). The value of parameters is chosen inside the range of biological values reported in table 2. The form of the basement membrane S B is taken regular and sinu- An example of the simulations is given in figure 4. The simulated patterns follow the ridges of the basal pattern and are elongated along its direction although initial conditions are circular. The pattern of figure 4h resembles the pattern of the acral nevus shown by dermoscopy in figure 4a. As expected, the nutrient distribution plays an important role in this matter, as growth is driven by nutrient diffusion in the avascular phase and only N is modified by the concavity of the surface S B . We focused on a small range of parameters in order to avoid the effect of a phase separation process [6], corresponding to a small value of D and an important value for e c . We chose a value of the height of the basal perturbation h 0 ¼ 1.0 l n (roughly 100 mm which is around the in vivo value). The parameters d, b 1 and k regulate the concentration of the melanocytic cells inside the lesion [28].
In this section, we have proposed a model and simulations exhibiting the importance of the geometry of the wavy basement membrane. Considering the small thickness of the epidermis, a two-dimensional description remains valid when the geometry is accurately taken into account. The fact that the equations are written in two-dimensional space dimensions is of utmost importance for simulations. Indeed, the basal geometry modifies mostly the nutrient repartition and by consequence the cell proliferation, with a higher cell density in the epidermal ridges. The cancerous cell repartition follows the basement geometry as shown in our numerical study and also clinically. However, one needs to remember that the signal of the dermascope used in clinical patient examination is not sensitive to the tumour cell concentration but to the melanin repartition in the whole depth of the epidermis. This melanin repartition also results from the basement geometry, where it is initially produced. It is why we also consider the transport of melanin from the basement membrane up to the upper layers of the epidermis.

Melanin transport in the epidermis
Several models both experimental and mathematical have been developed to describe the transport of molecules through the skin. To our knowledge, nothing has been done for the transport of melanin. Our model is based on the transport of the melanin vesicles convected by the displacement of keratinocytes. It means that we need to consider the dynamics of keratinocytes first.  with p the hydrostatic pressure, v the local rate of migration of keratinocytes proportional to the pressure gradient r p, K ¼ k/m, k the permeability m the viscosity of the tissue.

Model of keratinocyte migration and melanin transport in the epidermis
To demonstrate the generality of the observed mechanisms, we have considered in a second step a migration process described by a Stokes flow, which is more adapted to highly viscous and incompressible flow, and is given by the equations and From the basal layer S B , keratinocytes migrate along the normal to this surface (apical migration) imposing the boundary conditions for the two different flow descriptions v ¼ Vn on S B , ( 3 :5) where n is the normal to the surface of S B , and V the migration rate from the basal layer. The concentration of melanin in the basal layer c ¼ c 0 depends on the rate of synthesis of the melanin by melanocytes in this layer. It diffuses actively between the cells [29] and is convected with the migration of keratinocytes described by the velocity field v. Experimental studies suggest that melanin degradation occurs at a constant rate d m and is almost complete when reaching the skin surface SC [30,31]. The melanin concentration can therefore be described by a convection-diffusion equation This type of flow is well suited to the description of a viscous flow in a porous medium and has the advantage that it can be treated analytically [34]. The velocity field v can be represented by a velocity potential f satisfying the following equations: and with C (0) the average pigmentation and C (1) . 0 the amplitude of the corrections due to geometrical constraints on melanin transport. Equation (3.17) shows that C(x) is in phase opposition with z B (x), indicating that even in the case where the production of melanin is uniform in the basal layer a parallel pigmented pattern appears with a stronger pigmentation above the epidermal ridges. This result shows that the patterns observed in melanocytic lesions of glabrous skin can be explained by geometric constraints on melanin transport. The keratinocyte migration is faster above the epidermal ridges due to the local pressure (figure 6), as a consequence the vertical advection and dispersion of melanin is therefore increased.
To demonstrate the generality of the mechanism for appearance of the pigmented parallel pattern illustrated above, we consider now a keratinocyte migration described by a Stokes flow.

Stokes flow
We assume a similar form for S B and the same boundary as previously v ¼ Vn on S B , i.e. the same apical migration con- such that the incompressibility constraint can be satisfied. A numerical simulation of discretized equations (3.3) and (3.4) has been performed in Python using a relaxation method (appendix A.3.2). The results are shown in figure 7. As in the case of Darcy flow, the keratinocyte velocity is more important above the epidermal ridges than above the furrows. This difference of velocity causes a faster advection of melanin above these regions and the emergence of parallel patterns in skin pigmentation, even for a uniform production of melanin in the basal layer. In the case of the thick skin of hairless areas, H ) V/d m 400 mm, the concentration of intact melanin in the upper layers of the epidermis is negligible, the result presented here is independent of the form chosen for the interface SC and the results are qualitatively the same for other values of the parameters. The numerical studies confirm that the emergence of parallel grooves can be explained only by the geometric constraints imposed on the transport of melanin as shown in figure 7.

Localization of melanin column and stress inhibition of keratinocyte proliferation
The dermal-epidermal junction of the hairless skin consists of a periodic alternance of narrow ridges (intermediate ridges) and wide ridges (limiting ridges). Even if nevus cells tend to form aggregates in these two types of ridges, only the melanin produced in the limiting ridges is transported to the upper layers of the epidermis eventually making melanin columns (figure 8). We adapt here the geometry and the boundary conditions of our model to provide an explanation for this phenomenon.
To take into account the accumulation of melanocytes in epidermal ridges as suggested by the study of §2, the following boundary conditions on S B is chosen for the concentration of melanin: with Z nest a constant describing the vertical extension of the melanocyte nest. Many experiments have highlighted the important regulatory effect of mechanical stress on cell proliferation (see [35] for instance). We take now into account the inhibition of keratinocyte proliferation by the mechanical pressure assuming a homeostatic pressure p 0 in the basal layer. A local pressure p . p 0 leads to a decrease in proliferation and migration velocity, and a pressure p , p 0 stimulates proliferation and increases the rate of migration [36]. Boundary conditions on S B for the Stokes flow are therefore modified as ( 3 :20) and the boundary conditions on SC are unchanged. The numerical solution for melanin concentration obtained in a realistic epidermis geometry is shown in figure 8. In the narrow ridges, fast keratinocyte migration induces a high pressure on the basal layer, due to the non-slip assumption on this layer, leading to an inhibition of the proliferation and a decrease of the migration rate. In a steady state, migration rates are therefore much lower in the narrow intermediate ridges, than in the wide limiting ridges. As discussed previously, the vertical dispersion of melanin increases with the migration rate of keratinocytes, which explains why our model predicts the apparition of melanin columns only above limiting ridges.

Discussion
In most cases, melanocytic lesions in hairless skin have a notable structure and exhibit pigmented parallel stripes, usually attributed to the specific geometry of the skin in these areas. However, there are differences in the position of the stripe network between nevi and melanomas. In nevi, it is believed that the pigmented stripes coincide with the surface furrows while it would be the contrary for melanomas, giving perhaps a possible diagnosis. The mechanisms explaining the appearance of these structures, however, remain largely mysterious [20]. Here, we present two biomechanical models which explain physically such structures: one is based on tumour growth modelling, the other on melanin repartition inside the epidermis, both of them being strongly influenced by the basal geometry.
The first model considers the avascular growth of the lesion due to cancer cell proliferation in layers strongly distorted by the basal geometry. Focusing on a periodic undulated basal at the origin of fingerprints on soles and palms, we establish first a model able to accurately take into account such structure, then we perform simulations. It turns out that the nutrient concentration which decreases in the lesion due to consumption by the proliferative cells is more important in the valleys compared with the crests, increasing the proliferation in the ridges but also the activity of the tumour cells like perhaps melanin production. This effect is locally dependent on the curvature, so increases in narrow ridges compared with wide ridges.
Our second model focuses more on the transport of melanin in such glabrous areas of the skin. For that we take into account the specific transport of this pigment synthesized in the basal layer, advected to the upper layers of the epidermis by the apical migration of keratinocytes and degraded at a constant rate. The model has been solved analytically and numerically in the case of uniform melanin production on an undulated basal layer. It shows that the geometry of the basal surface influences the migration rate of keratinocytes which is predicted to be faster above epidermal ridges and leads to the apparition of a parallel pattern in the apparent skin pigmentation, with stronger pigmentation above the rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20140339 epidermal ridges. This corresponds to the clinical observation for nevi and shows that this pigmentation pattern can be explained only by the constraints imposed by the geometry of the epidermis on the keratinocyte migration. However, we note that histopathological studies show that melanocytes tend to aggregate and the assumption of uniform melanin synthesis is therefore probably not true as suggested by our first model. Parallel patterns along ridges are due to the combination of aggregation of melanocytic cells and melanin transport as suggested by our two models.
In nevi, nevus cells can be found in both intermediate and limiting ridges [17,20], but it was noted that only the melanin synthesized in limiting ridges is transported to the upper layers of the epidermis to form eventually a melanin column forming parallel furrow pattern. Using our model of melanin transport, we have proposed an explanation of this phenomenon based on the inhibition of keratinocyte proliferation by mechanical stress. The epidermis being treated as a viscous medium for migration, keratinocyte migration in narrow ridges lead therefore to significant constraints in the basal layer, attached to the dermal-epidermal junction, which inhibits cell proliferation and reduces the cellular migration rate as well as the advection of melanin. This explanation has been validated numerically and the results correspond to the melanin distribution observed in histopathology. As for melanomas, melanin transport should obey the same rules as in our model and as for nevi.
Therefore, the position of parallel ridges for melanoma cannot be explained by physical reasons alone. Only biological facts can justify them and perhaps be included in our study in the future. However, they are still missing to our knowledge. For instance, one questionable hypothesis would be the boundary condition on melanin distribution in equations (3.18) and (3.19). We can argue that the increase of nutrients in intermediate ridges would result in a higher cell activity and melanin production, but it is not sure that this physical argument is enough to justify the dichotomy between melanomas and nevi. However, more research on the biology of melanocytes in situ is necessary. Indeed, the melanocytic precursor cells originate from the dermis [37,38]. A possibility would be a preference of melanoma stem cells to migrate in intermediate ridges. These models enlighten the dominating mechanisms behind the parallel patterns, but they need to be improved by more biological information on the location and activities of the different melanocytic cells in order to be able to accurately predict the behaviour of a lesion and fully explain the significance of the current diagnostic tools. Using the same notations as in §2.1, we describe the metric allowing us to differentiate on a curved surface. Following Ogawa [26], we write (G ij ) the three-dimensional metric tensor at first order in e G xy ¼ @X @x Á @X @y ¼ @x @x Á @x @y þ z @x @x Á @n @y þ @x @y Á @n @x , Therefore, we can write (G ij ) as and with g ij is the metric tensor of the surface S B as and k ij is the second fundamental tensor, At the first order in e, we find with H the mean radius of curvature of S B , In our case, we only consider a dependence on x. So the expression of (g ij ), g and H is as follows:

A.2. Nutrient diffusion
Here, we detail the integration along the z-axis that allows a two-dimensional description of nutrients diffusion. The amount of nutrient contained in a section of surface S 0 of the thin layer is written as (A 10) Therefore, we use the average nutrient concentration in two-dimensional N(x, y) at (x, y) We choose a nutrient concentration depending weakly on the z-axis but with sharp variations n(x, y, z) ¼ n 0 (x, y) þ e 2 n 1 x, y, z In order to find a two-dimensional description of equation (2.1), we multiply it by (G/g) 1/2 and integrate it over z between 0 and e. The first term gives us G 1=2 G xx @n @x þ G xy @n @y þ @ @y G 1=2 G yy @n @y þ G yx @n @x dz |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} horizontal diffusion , (A14) D z perpendicular diffusion constant 10 4 d 21 [5] n s , N 1 , N 2 oxygen partial pressure in the skin 3320-10 400 Pa [46] l n nutrient penetration length 0.04-0.18 mm where (G ij ) ¼ (G ij ) 21 . The horizontal diffusion is equal to D 0 n, where D 0 is the standard Laplace-Beltrami operator.

(A 18)
A.3. Numerical scheme A.3.1. Scheme of melanocytic cells growth The numerical simulation are done by discretizing equations (2.8) and (2.9) on a two-dimensional lattice. In the case of the simple form of the surface chosen (h(x, y) ¼ h(x)), we consider the nutrient field n k i,j , the cells field f k i,j , the interaction field S k i,j and the derivation fields ffiffi ffi g p k i,j ( ffiffi ffi g p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi det(g ij ) p ) and H k i,j (the mean radius of S B ), where k is the time index, and i and j are the grid position. The first step of the simulation is to calculate the value of S k i,j , we use the following scheme for the Laplace -Beltrami operator: In order to express equation (2.8), we use the same scheme to express D 0 and we introduce a time dependence with a different timescale to ensure that growth is slower than diffusion n kþ1 i,j À n k i,j Dt ¼ D 0 n k i,j À k(1 À f k i,j )n k i,j À n k i,j f k i,j þ b 1 (n 1 À n k i,j ) þ b 2 H k i,j (n 2 À n k i,j ), (A 20) we choseDt ¼ 5Dt. We use the same approach to implement f kþ1 i,j using Dt instead ofDt. These schemes were written in CUDA and the results were drawn using Python.

A.3.2. Scheme of the Stokes flow model
We use a relaxation method to solve the stationary Stokes equation on the domain V with the following algorithm written in Python: -initialization of p i,j , v x,i,j and v y,i,j respecting the boundary conditions; p i,j is calculated to satisfy the Laplace equation Dp ¼ 0 using the relaxation method; v x,i,j and v y,i,j are determined by solving Dv x ¼ @ x p/m and Dv y ¼ @ y p/m using the relaxation method and -calculation of the incompressibility error dp ¼ r . v by finite differences. If the error does not satisfy the convergence criterium, i.e. S V dp 2 i,j e incomp , we implement p, pÃ i,j ¼ p i,j þ dp i,j and return to the third step. Otherwise the resolution is over.