Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture

Clinical diagnosis of skin cancers is based on several morphological criteria, among which is the presence of microstructures (e.g. dots and nests) sparsely distributed within the tumour lesion. In this study, we demonstrate that these patterns might originate from a phase separation process. In the absence of cellular proliferation, in fact, a binary mixture model, which is used to represent the mechanical behaviour of skin cancers, contains a cell–cell adhesion parameter that leads to a governing equation of the Cahn–Hilliard type. Taking into account a reaction–diffusion coupling between nutrient consumption and cellular proliferation, we show, with both analytical and numerical investigations, that two-phase models may undergo a spinodal decomposition even when considering mass exchanges between the phases. The cell–nutrient interaction defines a typical diffusive length in the problem, which is found to control the saturation of a growing separated domain, thus stabilizing the microstructural pattern. The distribution and evolution of such emerging cluster morphologies, as predicted by our model, are successfully compared to the clinical observation of microstructural patterns in tumour lesions.


Introduction
The evolution of a tumour mass is often characterized by the emergence of inhomogeneous patterns or clusters of cancerous cells, such as the papillary or cribriform structures in ductal carcinoma [1]. This morphological feature has diagnostic relevance for skin cancers, where, for example, large (small) blue-grey nests (globules) and leaf-like structures are typically found in basal cell carcinoma [2], as shown in figure 1(A). Similarly, melanoma skin cancers often present black globules composed of pigmented cells (see figure 1(B)), which are key factors in the early detection criteria of malignant lesions [2,3]. In recent decades, the continuous increase in the melanoma incidence rate [4], together with the decrease in survival rates after late detection [5], has led to extensive clinical studies for improving the diagnosis algorithms in clinical diagnosis [6,7]. However, most of these efforts focused on empirical statistical studies [2] and machine learning approaches [8], and only a few of them focused on understanding the physical origin of the morphological patterns of tumour lesions. Recent works have given a first physical insight into contour irregularities in these lesions, demonstrating that the competition between cellular adhesion and nutrient diffusion can drive a boundary instability during tumour growth [9,10]. The goal of this paper is to investigate the mechanisms responsible for the occurrence and evolution of cancer cell clusters during the early development of skin cancers.
The spontaneous sorting out and aggregation inside a mixture of different cell types has been observed since 1955, with the experience of Townes and Holtfreter [11] on mixtures of amphibian embryonic cells of different types. Starting from an initially random distribution, such cells were found to aggregate, forming an organized tissue with well-separated domains containing a specific cell type. The differential adhesion hypothesis is today a widely accepted explanation for this phenomenon, as discussed in the review by Steinberg [12]. Based on the analogy with the classical spontaneous phase separation in immiscible fluids [13], such a hypothesis assumes that the driving mechanism is the differential cell-cell adhesion between the various cell species, resulting in an effective surface tension between the cell type domains. Extensive experimental [14] and theoretical [15] investigations have highlighted the importance of cell sorting by cell-cell adhesion in embryogenesis [16,17], in tissue homeostasis [18] and, more recently, in cancer progression [19][20][21], especially during metastasis [22].
In order to explain the formation of aggregated cellular clusters, many theoretical approaches have been performed using discrete models. Extended Potts models, for example, consider cells arranged on a grid and local rules for cell migration and replication [23]. Such models were used by Turner and Sherratt [24] to show numerically the influence of cell-matrix adhesion on tumour invasiveness, as well as the possible split of a tumour mass into two distinct clusters. Multiscale models often employ a discrete approach in which internal cell processes and diffusion on a continuous space can be taken into account, simulating a wide range of processes, as done by the software CellSys recently developed by Hoehme and Drasdo [25]. Ramis-Conde et al [21] used such a multiscale approach to show the influence of specific cell adhesion molecules (CAMs) in regulating the dynamics of cancer invasion. Nevertheless, discrete models are only adapted to relatively small systems due to long computational times, and the extent of their result is limited to numerical simulations. Although being simple, some 4 continuous models have been developed [26] in order to overcome this limitation, providing an analytical insight into the aggregation process. Recently, Klein et al [27] proposed a model based on the Cahn-Hilliard equation to explain both analytically and numerically the existence of clusters for stem cell repartition in the epidermis.
On the other hand, multiphase models provide a realistic framework based on the mixture theory and have been widely adopted for tumour growth modelling, as recently reviewed by Lowengrub et al [28] and Preziosi et al [29,30]. These models describe biological tissues as a mixture of multiple phases (e.g. different cell types, extracellular matrix and interstitial fluid) and are potentially valuable in the medical context for simulating cancer growth inside various tissues on large scales. However, very few analytical studies have been performed to study the mathematical properties of these models, as most of the work has been focused on numerical simulations in specific geometries [28,31,32]. In particular, Lowengrub and coworkers reported that early multiphase models include a dependence of cell-cell adhesion on the local cell volume fraction, leading to backward parabolic equations [33] with features similar to a spinodal decomposition in classical phase separation theory [13,34,35]. Wellposed models with governing equations of the Cahn-Hilliard type have been developed in order to correct previous inconsistencies [36], but nothing is known about the dynamics of the spinodal decomposition and the phase ordering in such models. Recalling the importance of cell clustering in skin cancer diagnosis, we aim at analysing the occurrence of spinodal decomposition and cell aggregation in a multiphase model of tumour growth.
The paper is organized as follows. In section 2, we provide a brief review of skin cancer biology before introducing a two-phase model for a skin tumour with a general form of cell-cell interaction. In section 3, we study the stability of the mixture, both analytically and numerically, first neglecting mass exchanges between the phase and then accounting for cell proliferation in a reaction-diffusion process. The phase-ordering dynamics and the resulting patterns are compared to the results of earlier models of spinodal decomposition in other fields. Finally, in section 4, we discuss the extent to which our results have possible application in skin tumour diagnosis.

Skin cancer biology
The skin is made of three main layers: the hypodermis, the dermis and the epidermis. In this work, we focus on the early stages of skin cancer development, occurring in the two upper layers. The dermis contains blood and lymphatic vasculature, hair follicles, glands and nerves. It is a very elastic tissue due to an extended network of collagen and elastin fibres. The epidermis is the outermost layer of the skin that acts as a protective barrier from the external environment. It is mainly constituted by a type of cells called keratinocytes (see figure 2). The keratinocytes proliferate on a basal membrane separating the dermis from the epidermis, and migrate to the skin surface while differentiating, ultimately creating a layer of dead cells and lipids, the stratum corneum. Keratinocytes adhere strongly to each other via desmosomes to provide impermeability to the tissue. Skin is protected from sun damage by another cell population, the melanocytes, which are located at the basal membrane, provide melanin to the surrounding cells and are responsible for skin pigmentation. Sun exposure [38] or virus infection [39] can cause deregulation of the epidermal cells that start proliferating Schematic representation of the two first layers of the dermis. The epidermis is the outermost layer and is mainly composed of 90% keratinocytes and 5% melanocytes responsible for skin pigmentation. The epidermis is separated from the dermis by a basement membrane, a thin sheet of fibres. Nutrient supply to these two layers is provided by diffusion from the dermis vasculature and, in the case of oxygen, from the atmosphere [37]. Skin cancer originates, in particular, from keratinocytes (basal cell carcinoma) and melanocytes (cutaneous melanoma) that proliferate abnormally and invade the epidermis and the dermis after penetration through the basement membrane.
abnormally, leading to uncontrolled invasion of the surrounding tissue. Most of the resulting tumours are classified into two types: non-melanoma (NMSC) and melanoma skin cancer. Basal cell carcinoma (BCC) is the most common NMSC and originates from the stem cells of the epidermis, lying on the basal membrane [40]. It is often clinically characterized by a lesion with large blue-grey nests, small blue-grey globules and leaf-like structures [2]. Although it has a low probability of creating metastasis, its high incidence rate is a major problem for health care organizations [41]. Squamous cell carcinoma (SCC) is another NMSC that originates from the stratified epithelium over the basal cells. Although less common than BCC, it has a much larger metastasis rate. Finally, melanoma skin cancers originate from melanocytes, sometimes from pre-existing benign lesions called nevi (as depicted in figure 1(C, D)). Melanoma usually develop first in the epidermis, undergoing a radial growth phase before vertically invading the dermis. They represent less than 5% of skin cancers but are responsible for more than 75% of skin cancer-related deaths. Irregular patterns and a non-uniform distribution of globules or clusters of pigmented cells are key features for melanoma diagnosis [3].
We are interested in understanding the influence of specific cell-cell adhesion and the tumour microenvironment on the occurrence of such microstructural patterns resulting from cell aggregation. For this purpose, we provide in the following a theoretical investigation of phase ordering phenomena in a multiphase tumour model.

Binary mixture model of a skin tumour
For the sake of simplicity, we consider a simple multiphase model of the tumour lesion: a mixture of a cancerous cellular phase (arising from keratinocytes or melanocytes) with volume fraction φ c and averaged local velocity v c , and a liquid phase (interstitial fluid) with volume fraction φ l and velocity v l also containing dead cells and eventually other cell types. These two phases fill all the available space imposing the saturation constraint φ c + φ l = 1. Assuming that the density of both phases is equal to the water density ρ, the mass conservation is given by the following equations: with ρ c = −ρ l the mass transfer rate from the liquid phase to the cellular phase, resulting from cell proliferation and death (or shrinkage). The saturation constraint also imposes the following incompressibility condition: Tissue reorganization is here only a result of tumour cell proliferation and mechanical interaction, under the saturation and incompressibility constraints. Cells adhere to each other via cell adhesion molecules (CAMs) located at their plasma membrane. In particular, normal melanocytes express proteins for adhering and communicating with keratinocytes (E-cadherin, desmoglein and connexin). These proteins are down-regulated during melanomagenesis, while proteins for melanocyte-melanocyte communication and adhesion are up-regulated (N-cadherin, Mel-CAM, integrin, ALCAM and connexin) [42]. In contrast, healthy keratinocytes strongly adhere to each other via desmosomes, which become strongly down-regulated in BCC cells, resulting in a less dense tissue [40]. Following Wise et al [36], we take into account these interactions between cells into the free energy of the cellular phase, assuming a weakly nonlocal dependence on the cellular phase density φ c where ψ is the free energy per unit of volume for a homogeneous tissue and the term in 2 represents a surface energy penalizing large gradients of cellular concentration [36]. We assume that the viscous interactions between the cells and the liquid phase are the main source of energy dissipation in the system, which therefore reads where M is a constant friction parameter. Measurements of the Darcy law in rat skin allow us to estimate this parameter to be in the range of M ∼ 963-11571 mm −2 Pa day [43,44]. Following Doi and Onuki [45], we use Rayleigh's variational principle indicating that the system dynamics can be obtained by minimizing the Rayleighian, defined by Doi and Onuki as R = W + dF c /dt, with respect to the velocities v c and v l . In order to account for the incompressibility constraint given in equation (3), we introduce a Lagrange multiplier p before minimizing the following functional: When integrating by part the variation of free energy (the second term of equation (6)) border terms have to be considered at infinity. The governing equations are therefore written: which can be interpreted as the force equilibrium in each phase, assuming a system with negligible inertia. Eliminating the pressure p in equations (8) and (7) the relative movement of the two phases is given by In the absence of external forces and considering a high viscosity for the mixture, symmetry requirements impose that the centre of mass of the mixture is not moving, i.e. φ c v c + φ l v l = 0, which is a particular incompressible motion. In this case, the velocity of the cell phase is given by a Darcy-like law, as follows: with the excess pressure exerted by the cells defined by In the following, f (φ c ) = ∂ψ/∂φ c denotes the pressure in a homogeneous system. For physical consistency, cell-cell interactions should be attractive at a moderate cell volume fraction ( f (φ c ) < 0 for φ c < φ e ), because of cell adhesion, and repulsive at a high volume fraction ( f (φ c ) > 0 for φ c > φ e ), as depicted in figure 3. Jain [46] reviewed the cell volume fraction in human epidermis, giving φ e = 0.63-0.87 assuming that the tissue is close to mechanical equilibrium. Moreover, adding equations (7) and (8) gives ∇ p = −∇ , showing that the interstitial fluid pressure p and the excess pressure generated by the cells have the same typical value χ. The review by Jain reports the interstitial fluid pressure in various tissues and tumours, giving the estimate χ ∼ 133-1330 Pa. The term / √ χ gives a characteristic distance of interaction between cells. An estimate of this interaction is given by the cell size, typically 6-20 µm for melanoma cells [47], but can eventually be larger when taking into account paracrine interactions. The dynamics of the cellular phase is given by Finally, we investigate the effect of nutrient limitation on cell growth in the system. In a first approximation, we assume that the cell proliferation c depends linearly on the concentration n of a diffusing nutrient (e.g. oxygen or glucose): where n s is a typical nutrient concentration inside the tissue. In particular, Creasey et al [48] measured the doubling time of primary melanoma cells in culture giving γ c n = 0.2 day −1 , while Bedogni and Powell [49] have reported the behaviour of normal and tumour cells in various hypoxia conditions, leading to typical values of δ c = 0.1-0.33. Nutrients are delivered from the vasculature (and from the atmosphere in the case of oxygen), diffuse into the interstitial fluid with a diffusion constant D n and are consumed by the cells at a constant rate, δ n . In particular, Johnson et al [50] measured oxygen diffusion in the healthy epidermis with D n = 39.7 mm 2 day −1 and Stucker et al [37] have reported skin oxygen consumption in the range of δ n = 1190-3310 day −1 . Due to the fast diffusion as compared to the growth velocity, we assume the following expression for the nutrient concentration at equilibrium: where S n is the diffusion constant from the source of nutrients, which can either be the blood supply or the atmosphere for an avascular skin tumour in the epidermis (e.g. radial growth phase of melanoma, see [51]), and n s is the nutrient concentration at the source. Stucker et al [37] measured the variation of oxygen partial pressure in healthy skin, giving the estimate n s = 3320-10400 Pa.

Dimensionless equations
With the aim both to provide an analytical study of the governing equations and to perform numerical simulations, we define the dimensionless quantities: x = x/l n ,n = n/n s ,t = tγ , where l n = √ D n /δ n is the nutrient penetration length. Dropping the hats and the index for the cell phase, the dimensionless governing equations can be written as An estimation of such dimensionless parameters is given in table 1 from the experimental data in the biomedical literature.

Linear stability analysis of the binary mixture and phase ordering
In this section, we consider an initially homogeneous binary mixture, with cancer cell volume fraction φ 0 , and we investigate the occurrence of a spinodal decomposition (spontaneous phase separation). We perform a linear stability analysis on our multiphase tumour model for an infinite system. We then compare the results with simulations realized on a system with periodic boundaries.

Linear stability analysis in the absence of mass exchanges
In the absence of mass exchanges between the phases ( c = 0), equation (19) can be simplified by re-scaling space and time byx = x/ andt = t D/ , and dropping the hats, we obtain the following relation: We consider an initially spatially uniform cell density φ 0 and an infinitesimal perturbation given by Substituting this expression into equation (21), at the order we obtain the growth rate λ of the perturbation: If f (φ 0 ) < 0 (i.e. for a cell volume fraction φ 0 < φ * as shown in figure 3), the binary mixture is always unstable in the long-wavelength domain k < √ − f (φ 0 )/ , and the system undergoes a spinodal decomposition for all values of parameters and D. Note that here the strength of the cell-cell interaction coefficient D plays a similar role as the temperature in the common phase separation models. In fact, the governing equation is equivalent to the phenomenological equation for phase-separating systems with an order-parameter-dependent mobility M(φ), which corresponds here to φ K (φ), as described in [52,53]. After a transitory regime of phase separation, these systems are composed of single-phase domains organized into self-similar patterns, which are characterized by a typical length L(t) increasing with time. Given such a similarity with our multiphase tumour model, we expect to observe the same patterns and growth law in our system. However, we note that in our multiphase model the equivalent of the mobility, φ K (φ), vanishes at φ = 0 in an unstable domain ( f (φ) < 0 for φ 1), which is not the situation studied in [52,53], so that we might expect some differences in the dynamics of our system. We anticipate, as shown in figure 10, that the typical length distribution in our model is different from the one obtained from the Cahn-Hilliard equation.

Linear stability analysis including cell death and proliferation
When taking into account cell death and proliferation, the order parameter φ is not conserved anymore and therefore we expect a different phase ordering dynamics [35], arising from the presence of mass exchanges in equation (19). The homogeneous and stationary solution of equations (19) and (20) is given by and corresponds to a homeostatic state where cell division (resp. nutrient consumption) equilibrates cell death (resp. nutrient delivery). An infinitesimal perturbation of wavevector k around this solution can be written as follows: n = n 0 + δn exp(λt) cos(kx).
Substituting this expression into equations (19) and (20), the growth rate of the perturbation is determined by the following relation: The coupling with the diffusing nutrients is responsible for the last term in equation (27), which stabilizes the mixture against long-wavelength perturbations, as depicted in figure 4. This term introduces a second typical length in the system, which has been identified as a typical mechanism of pattern selection in phase-ordering systems, such as block copolymers and reaction-controlled separating mixtures [54]. As it has been observed in those systems, a saturation of the domain growth at a finite size is expected, and the order parameter should be conserved at large scales. The dispersion relation depicted in figure 4 is indeed similar to that derived by Klein et al [27] for their Cahn-Hilliard model of stem cell clustering in the epidermis. In their case a second typical length is given through contact inhibition and stem cell differentiation. We find here similar properties in a different modelling framework and expect comparable behaviours. Assuming 1 in equation (27), the saturation length in our model is given by From equation (27), we can derive the value of the wavevector k max , which gives the maximum growth rate during phase separation, and define a function g relating all the parameters, whose sign controls the stability of the binary mixture. The domain where g > 0 indicates the occurrence of a phase separation (see diagram phases in figure 5), as follows: The sign of g gives the phase diagram depicted in figure 5, where the unstable region (red domain) becomes less important here compared to the case without mass exchanges (outside of the grey domain). Interestingly, the stability of the mixture seems to depend mostly on the initial  Phase diagrams for the stability of a homogeneous mixture with an initial equilibrium cancer cell volume fraction φ 0 and nutrient concentration n 0 , given as a function of the dimensionless parameters δ (cell death rate) and β (nutrient supply rate) for = 0.2 (left) and = 0.084 (right) and for D = 2 (top) and D = 20 (bottom). In the grey domains φ 0 > φ * , and the mixture is stable as in the absence of cancer cell proliferation. The coupling with the diffusing nutrients has the tendency to stabilize the mixture and the unstable area is restricted to the red domain. Interestingly, the limit of this domain seems to correspond to surfaces of the constant φ 0 .
cell concentration φ 0 , as the limit of the unstable region corresponds numerically to a surface of constant φ 0 .

Numerical simulations neglecting mass exchanges between phases
Numerical simulations of equation (21) were performed in a two-dimensional N × N periodic lattice for different homogeneous distributions of φ 0 and an initially uncorrelated white noise with amplitude 0.01φ 0 . We discretized the equations using centred second-order accurate differencing schemes in space and explicit first-order accurate forward schemes in time. Moreover, we use a phenomenological form of f in that respects the biological and physical constraints we discussed in section 2 (see figure 3) where p must be strictly greater than one. We choose p = 2 in the simulations and φ e = 0.6. The order-parameter-dependent mobility M(φ) = φ K (φ) vanishes for φ → 0 such that the value of φ remains always positive, ensuring the physical consistency of our model [53]. However, numerical fluctuations introduced by the discretization of the equations in the simulations can lead to domains with negative φ. In order to solve this problem, we set φ to 0 whenever it becomes negative and we use a fine mesh ( x = 0.15 and t = 10 −4 ) to ensure that the average value of φ is conserved. In order to improve the computational time, we implemented our program on a graphic processors unit, GPU, using the CUDA parallel programming language. We used a GTX-580 Nvidia graphic card, and our computational time was consistently reduced. In figure 6, we depict the evolution of a system with an initial cell density φ 0 = 0.25 and a small initial uncorrelated white noise. We can see that the system exhibits two kinds of domains after a transitory regime: one is empty of cells (φ = 0) and the other rich in cells (φ = 0.54 ≈ φ e ), forming maze-like patterns. Figure 7 shows the patterns obtained for a lower initial cell density φ = 0.15 with the appearance of isolated clusters of cells. Considering that the total number of cells in the system is conserved, the dependance of the domain geometry on the initial concentration of cell φ 0 is predicted in the classical theory of spinodal decomposition in such systems [35,52]: isolated circular clusters of cells for small φ 0 and isolated circular domains empty of cells for φ 0 ≈ φ * . Moreover, a critical concentration φ t ≈ φ e /2 is defined as the concentration at which cell domains occupy half of the space. For an intermediate initial concentration φ 0 < φ t (resp. φ 0 > φ t ), maze-like patterns appear during a transitory regime, before reorganizing into isolated circular clusters of cells (resp. domains empty of cells) at long times. The duration of the maze-like pattern regime diverges as φ 0 gets closer to φ t , explaining the structure observed in figure 6 with φ 0 = 0.25 and φ t ≈ φ e /2 = 0.3. In order to make a comparison with existing results on phase ordering, we calculate the spherically averaged time-dependent structure factor S(k, t) = φ (k, t)φ(k, t) * , whereφ is the Fourier transform of φ and the average is taken on the angles of the wave vector, which gives access to the length distribution in the system. In particular, this distribution provides a definition of the typical length of the system at time t [52] (e.g. the width of the maze paths in figure 6 or the diameter of the clusters in figure 7) As depicted in figure 8, the structure factor exhibits here a self-similar scaling [55], as in classical phase separation models, which reads where G is a master function. The interpretation of the dynamical scalings is as follows: the patterns maintain their morphology while their characteristic length L(t) grows over time. Interestingly, our master curve seems to exhibit a much narrower peak compared to that described in [52] for phase separation with order-parameter-dependent mobility. This effect might be related either to a lack of data or to the mobility coefficient that vanishes in an unstable domain of φ for the binary mixture. The time-dependent lengthscale L(t) exhibits a power-law evolution, as shown in figure 8, with an exponent 0.37 close to the Lifshitz-Slvozov growth law L(t) ∼ t 1/3 , which is found in late stage phase ordering in conserved order parameter systems [35]. Finally, these simulations show that multiphase models may undergo spinodal decompositions and that they share the same properties of well-known models as phase ordering.

Numerical simulations considering cell death and proliferation
Taking into account a binary mixture with cell death and proliferation, the linear perturbation analysis has shown that a homogeneous distribution of cells might be unstable given certain conditions. Numerical simulations have been performed in a two-dimensional N × N periodic lattice solving equations (19) and (20) with an initial uniformly distributed concentration of cancer cells and nutrient. The initial values of φ and n are set by the homogeneous and stationary solutions φ 0 and n 0 , defined in equation (24), with an initial uncorrelated white noise for φ of amplitude 0.01φ 0 . In order to solve equation (20) numerically, we have rewritten it as follows: using a very fast relaxation timed 1 (in the simulationsd is set to 5). Figure 9 shows the influence of the initial cell density φ 0 on the geometry of the system, as for the case with c = 0. In the first row, a system with a low initial density φ 0 = 0.05 forms isolated circular clusters of cancer cells, while in the second row a high initial density φ 0 = 0.3 gives rise to isolated circular domains empty of cancer cells. Dots and globular patterns described in figure 1 for skin tumours can therefore be reproduced by a spinodal decomposition in a multiphase model. Furthermore, we note that the existence of an additional typical length l s stabilizes the phase-ordering pattern in a regular inhomogeneous distribution of cells (the left column of figure 9). The growth saturation of the domain can stabilize the pattern into isolated circular clusters, as shown in the first two rows of figure 9, or into a maze-like pattern if the saturation occurs in an earlier stage (l s close to the typical length of the phase separation) or if the system is close to criticality, as depicted in the third row of figure 9. In particular, the last row of figure 9 shows an extreme case where the pattern is stabilized just after the phase separation without a transitory growth regime, due to a saturation length l s = 9.1 (derived from equation (28)) close to the typical length of the phase separation l max = 2π/k max = 6.89 (k max defined as the maximum of λ in equation (27)).
After phase separation, we observe a transitory regime where the domain grows as t α with an exponent α that depends on the model parameters. We found that α = 0.38 for the evolution depicted in the second row of figure 9 and α = 0.35 in the third row for a different cell-cell  The re-scaled structure factors collapse on a similar master function for our model without cell proliferation and with cell proliferation before the domain size saturates. Note the difference with the re-scaled structure factors for the standard Cahn-Hilliard equation. The difference could be due to a lack of data for our model or to the effective mobility that vanishes in an unstable domain φ = 0 in our multiphase model. interaction ( p = 3 and p = 2, respectively). However, we found that during the transitory growth regime the structure factor exhibits a self-similar behaviour with size distributions collapsing on a master curve, which is similar to that found for c = 0 (see figure 10). As predicted by the theoretical analysis, the domain growth saturates when reaching a size of the order of the nutrient penetration length, which is the dimensionless unit in our simulation. From a biological perspective, when the diameter of a cellular cluster reaches twice the nutrient penetration length, cell proliferation stops in the cluster core due to a lack of available nutrients.

Discussion and comparison with biological data
Skin cancers, both melanoma and non-melanoma, often present microstructural patterns such as clusters and dots (figure 1). Here we investigate the emergence of such patterns using a binary mixture model for the tumour, also accounting for a realistic cell-cell adhesion mechanism, and we explain them as a phase separation between two cellular types. Our theoretical model predicts the occurrence of structures with a typical size of l s ≈ 2l n = 0.2 mm (see table 1), which are compatible with the size of tumour cell clusters observed clinically. Being much larger than a typical cell size of 6-20 µm for melanoma cells [47], this characteristic length is therefore compatible with a multiphase modelling of the system. As classically found for binary fluid mixtures under a critical temperature, our model predicts that cell segregation occurs when adhesion between cells of the same type becomes large enough. This result corresponds to the cadherin switch clinically observed during melanomagenesis, where melanocyte-melanocyte interactions are up-regulated, while melanocyte-keratinocyte interactions are down-regulated [42].
Neglecting cell proliferation and death, our multiphase tumour model is very similar to the Cahn-Hilliard equation but leads to a different expression for the mobility [52,53]. After the phase separation, due for example to a phenotypic change inducing stronger cell-cell interactions, tumour cells form cluster (maze-like) patterns for an initially low (high) density. We found that these microstructures enlarge with time according to the Lifshitz-Slvozov growth law L(t) ∼ t 1/3 . Cell proliferation, controlled here by the local nutrient concentration, introduces a second typical length in the system other than the nutrient penetration length, preventing long-lengthscale variation of the cell density. When reaching this typical size, the microstructure stops growing and we found numerically that the patterns stabilize into quite a symmetrical organization (circular clusters or holes). This is compatible with the observations of Xu et al [3], reporting that regularly distributed symmetric clusters are typical of benign nevi, which are non-evolving structures. In contrast, asymmetric clusters, such as those observed in our numerical simulations during the transitory growth regime, are typical of early melanoma development. Our theoretical model predicts that the presence of such microstructural patterns indicates that the lesion is evolving, probably due to a recent change in cell phenotype or stromal interactions. Furthermore, we found that the distribution of cluster sizes during the transitory growth regime (before saturation) is self-similar. We expect that the introduction of a cell proliferation regulation by contact inhibition might induce a similar effect as that found in models of reaction-controlled separating mixtures [54]. These systems have been shown to exhibit the same patterns and growth dynamics as those found in our binary mixture model. In a similar way the introduction of a paracrine regulation of cell proliferation might introduce a typical length, given by the diffusion of some paracrine factors, as investigated in [51]. Indeed Klein et al [27] used a model based on the Cahn-Hilliard equation to understand clustering of epidermal stem cells with contact inhibition. We found also cluster patterns (figure 9) for a different cell proliferation regulation mechanism and in a model taking explicitly into account mixture saturation. This confirms that phase separation models can be useful in understanding pattern formation in biological systems.
In conclusion, we analysed the properties of spinodal decomposition in a binary mixture model of skin tumours, showing the analogies existing between our results and the phase separation phenomena in block copolymers and in reaction-controlled separating mixtures. Our results relate changes in cell phenotype and stromal interaction to the spatial organization of tumour cells, the statistical distribution of shapes and sizes of cell clusters and their dynamics in a skin lesion. This study allows us to foster an understanding of the emergence of microstructural patterns in skin tumours and might ultimately help us to improve accuracy in early skin-cancer diagnosis.