Royal Society Publishing

The radial growth phase of malignant melanoma: multi-phase modelling, numerical simulations and linear stability analysis

P. Ciarletta, L. Foret, M. Ben Amar


Cutaneous melanoma is disproportionately lethal despite its relatively low incidence and its potential for cure in the early stages. The aim of this study is to foster understanding of the role of microstructure on the occurrence of morphological changes in diseased skin during melanoma evolution. The authors propose a biomechanical analysis of its radial growth phase, investigating the role of intercellular/stromal connections on the initial stages of epidermis invasion. The radial growth phase of a primary melanoma is modelled within the multi-phase theory of mixtures, reproducing the mechanical behaviour of the skin layers and of the epidermal–dermal junction. The theoretical analysis takes into account those cellular processes that have been experimentally observed to disrupt homeostasis in normal epidermis. Numerical simulations demonstrate that the loss of adhesiveness of the melanoma cells both to the basal laminae, caused by deregulation mechanisms of adherent junctions, and to adjacent keratynocytes, consequent to a downregulation of E-cadherin, are the fundamental biomechanical features for promoting tumour initiation. Finally, the authors provide the mathematical proof of a long wavelength instability of the tumour front during the early stages of melanoma invasion. These results open the perspective to correlate the early morphology of a growing melanoma with the biomechanical characteristics of its micro-environment.

1. Introduction

Skin cancer has possibly first been described by the Greek physician Hippocrates around 400 BC [1], who mentioned in his writings about black herpetic type lesions as μɛ̀λανoζ- (‘dark’) -ωμα (meaning literally ‘inhuman, hard’, and later used in Medicine as a suffix for ‘tumour’). In Western medical literature, the first skin tumour was reported by the physician John Hunter in 1787, who made the diagnosis on a skin sample which is preserved inside the Hunterian museum in London [2]. Skin cancers are classified into two general types, non-melanoma and melanoma, having an overall world incidence which has been increasing over the past decades. Currently, the World Health Organization estimates that between 2 and 3 million non-melanoma skin cancers and about 130 000 melanoma skin cancers occur globally each year.1 The American Cancer Society reports that one in every three cancers diagnosed is a skin cancer,2 meaning that one in five Americans will develop skin cancer in the course of a lifetime [3]. Although the mortality rate from non-melanoma skin cancer is relatively low, significant morbidity occurs from tissue destruction, functional impairment, and disfigurement [4]. Worldwide incidence rates for cutaneous melanoma have risen faster than those for any other malignancy in Caucasian populations over the last 30 years [5]. Furthermore, even though cutaneous melanoma accounts for about 3 per cent of skin cancer cases, it causes more than 75 per cent of skin cancer-related deaths. Melanoma has proved to be disproportionately lethal despite centuries of awareness of the disease, its relatively low incidence compared with non-melanoma skin cancers and, above all, its potential for cure in the early stages [6]. In fact, the depth of a diagnosed melanoma lesion (invasion thickness according to Breslow) is an indicator of its prognosis. Early diagnosis of melanoma is of utmost importance because if it is made at an early stage, when the tumour is still confined to the epidermis (or less than 1 mm thick) the prognosis is excellent and the 10-year survival rate is estimated to be between 90 and 100 per cent. If the diagnosis is made at a more advanced stage the 5-year survival rate drops to 10–15% [7]. The clinical diagnosis of cutaneous melanoma is based on the morphological features of the skin lesion, according the ‘ABCDE rule’ (asymmetry, border irregularity, colour variegation, diameter and evolution over time). This clinical algorithm is now used worldwide, having sensitivity of 65–80%, but with a tendency to fail in melanoma smaller than 5 mm in diameter [8]. Using dermoscopy (also called epiluminescence microscopy), clinical dermatologists have access to subsurface structures of the skin for in vivo examination, with an increase in the diagnostic accuracy from 65–80% to 70–95%, depending on clinician and/or dermoscopist training, respectively [9]. Nevertheless, some pigmented skin lesions, especially early melanomas, may lack specific dermoscopic features and are difficult to diagnose even with dermoscopy. In other words, as we go back in tumour development we are approaching the limitations of dermoscopy [10].

In order to develop more accurate diagnostic techniques, it is fundamental to identify the interacting mechanisms which control the early growth of skin cancer. Various classes of models have been proposed in order to develop mathematical frameworks dealing with different aspects of tumour evolution at scales ranging from molecular to tissue level [11,12]. While discrete models are particularly useful for studying carcinogenesis, genetic instability and mutual interactions of individual cells, continuum methods are employed at larger scales of investigation, where cancer-related variables can be described as continuous fields by means of partial differential and integro-differential equations [13,14]. A number of recent researches in cancer modelling has focused on a multi-phase continuum description of solid tumour evolution, accounting for structural heterogeneities, non-uniform micro-environment and spatial limitation of available nutrients [1517]. Although such theoretical frameworks have proved their consistency for describing some crucial aspects concerning the growth of avascular solid tumours, a conceptual framework for modelling the peculiar micro-environment of cutaneous melanoma is lacking. The scientific aim of this article is to foster understanding on the role of micro-structural interactions on the initial progression of melanoma, and to study the occurrence of morphological changes in diseased skin at a mesoscopic level. This article will focus in particular on the identification of the biomechanical parameters determining the dynamics of tumour evolution and the occurrence of shape instabilities.

In the following, §2 will be devoted to a comprehensive description of the biomechanical function and structure of human skin, both in healthy and pathological conditions. Section 3 will be concerned with the theoretical formalism and the biological hypothesis for defining a multi-phase model of the radial growth phase for a cutaneous melanoma. The numerical results of such a biomechanical model will be collected and discussed in §4. The issue of contour instability during the early growth of skin cancer will be treated in §5, performing a linear stability analysis of the perturbed tumour front in the multi-phase model. Finally, in §6 the authors analyse the main results of the proposed multi-phase theory, evidencing the novelties and the limitations of the modelling approach.

2. Biomechanical functions and structure of human skin

Skin is a very complex organ both from the structural and the biochemical viewpoints. Protecting the human body as its outer covering, skin represents about 20 per cent of total body weight and serves many other physiological functions, such as thermoregulation, insulation, pain sensation, touch and chemical conversion. As shown in figure 1a, three different compartments can be defined in the layered structure of human skin [19]: the epidermis, an outer epithelium of stratified cells; the dermis, an intermediate cushion of vascularized connective tissue; and the hypodermis, the lowermost layer made of loose tissue and adipose cells.

Figure 1.

Schematic of the layered structure of human skin (a), and of the cellular organization in the lower epidermis (b). The epidermal structure includes keratinocytes (Kc), melanocytes (Mc) and malignant melanoma cells (mm), the latter shown both in the proliferating (dark grey) and in the necrotic state (grey). The reticular dermis is attached to the lower epidermis through the basal laminae, being characterized by a dispersion of fibroblasts (Fb) inside the ECM network of collagen type I fibres (Coll). The epidermal structure on the right side has been partially drawn using the schematics employed by [18].

The epidermis is an avascular, innervated layer made of four different cellular types [20]: keratinocytes, melanocytes, Langerhans cells and Merkel cells. The keratinocytes are the most abundant cellular type in this epithelium, representing about 80 per cent of the epidermal cells. Migrating from the basal layer to the outer skin surface over a period of about 26 days, the keratinocytes undergo significant morphological changes while ensuring three fundamental functions: providing cohesion to the epithelium, forming a barrier from the exterior and protecting against excessive light exposure. The melanocytes are the second most abundant cellular species of the epidermis, representing about 5–10% of the cell population. They are localized at the basal level, having the function of synthesizing melanin, a pigment serving as a photoprotector. The Langerhans cells represent about 4–8% of the cellular population in the epidermis, being dendritic cells with the function of processing microbial bodies to gain full functionality of antigen-presentation to T-lymphocytes. The Merkel cells represent the minority of the cell types in the epidermis, with the primary function of serving as mechanoreceptors, being associated with the sense of touch discrimination of shapes and textures.

The dermis layer of human skin is composed of collagen, elastin fibres, an amorphous extracellular matrix, and fibroblasts. It is separated from the epidermis by the basement membrane, or basal laminae, which provides both the adhesion and the connection network between the two layers. The dermis contains both blood and lymphatic vessels, together with mechanoreceptor/nerve endings and a number of cutaneous adnexal structures (e.g. sebaceous and apocrine glands). Its architectural organization provides the fundamental mechanical functions to support and to provide strength to the whole skin [21]. Collagen accounts for more than 70 per cent of the dry weight of the skin, with a predominance of collagen type I triple helices, providing a significant increase in the mechanical stiffness, together with collagen types III and V [22]. Elastin represents about 4 per cent of the dry weight of skin and forms a network of small fibres between collagen, having the role of mechanical recoil after deformation. While these two elastic fibres are responsible for strength and energy recovery of the skin, the ground substance in the extra-cellular matrix (ECM) determines the macroscopic rheological properties [23]. In fact, such a component is formed by polysaccharides and proteins, which form macromolecules with a marked capacity of holding water, accounting for about 20 per cent of the dry weight of skin. From a structural viewpoint, the dermis can be divided into two regions, according to a different architectural organization. The reticular dermis, or stratum papillare, is the upper, superficial region made of loose connective tissue, which consists of numerous small, vascular eminences, and the papillae, which rise perpendicularly from its surface. The deeper reticular dermis, or, stratum reticulare, is made of dense connective tissue which may contain some adipocytes. Although thicker than the reticular dermis, the extent of its thickness varies significantly among different parts of the body [24]. The hypodermis, or subcutis, is found deep in the reticular dermis, without a clear separation structure, with the purpose to attach to the underlying skeletal muscle (excluding eyelids, ears). It is mainly constituted of globules of adipose structures, with the main function of providing the insulation for the human body.

2.1. Cell differentiation in the layered structure of epidermis

Epidermal thickness varies in the range from 0.1 to 1 mm (for palms/soles) depending on the anatomical location and on the presence of hair. The differentiation of keratinocytes from the inner to the outer layers is responsible for the stratified morphology of epidermis, as visible from optical microscopy [25]. Five different structural organizations of the keratinocytes can be distinguished, as indicated in figure 1a. The outer epidermis is composed of three layers: a stratum corneum, which is the outer portion made of multiple layers of flake-like cells filled with keratin and without nuclei; a stratum lucidum, which is a homogeneous zone with cells containing the remains of ruptured keratohyaline granules; and a stratum granulosum, with 1–3 layers of flattened cells which contain keratohyaline and lamellar granules. The lower epidermis, also known as the Malphigian layer, is made of a thick arrangement of spinous prickle-like cells, the stratum spinosum, which lie on a single row of mitotically active cells, forming the stratum germinativum, or basal layer. Cells in the basal layer, whose shape varies from columnar to cuboidal, are undifferentiated and can be considered as the epidermis stem cells. Once formed through mitosis, cells move up in the upper layers producing the above cited characteristic structures as they change shape during their differentiation process. Melanocytes, which are neural crest-derived cells, are only situated in the basal layer, where they interact with adjacent basal keratynocytes and the anchoring sites of the basement membrane, as shown in figure 1b. The characteristic cell appearance in the stratum spinosum is owing to the presence of desmosomal attachments between keratinocytes. Such desmosomes (maculae adhaerentes) are formed by two plaques, composed of desmoplakin, which are joined together through transmembrane glycoproteins (desmogleins and E-cadherins), while tonofilaments attach each plaque to the cytoplasm of the associated cell [26]. During keratinization, such intermediate filaments are crosslinked by filaggrin to form large aggregates concentrated at the cell periphery, with crucial functions of structural integrity and stability of the epithelium. Basophilic granules, called keratohyalines, are produced by the Golgi apparatus of keratinocytes in the stratum granulosum: they are rich in sulphated amino acids and contain lamellar bodies consisting of glycolipids. Such granules later migrate towards the cytoplasmic membrane with which they fuse at the interface between the stratum corneum. During this fusion process, polar lipids and proteins are massively discharged in the intercellular space, localizing at the lower stratum corneum, where they form a barrier to water loss [27]. At this stage, this process involves breakdown of the nucleus and organelles, as well as thickening of the plasma membrane. Finally, Merkel cells are located uniquely within the basal layer, while Langerhans cells in the epidermis are found in the stratum spinosum.

2.2. Basement membrane: ultrastructure and anchoring functions

Basement membranes are present in most tissues as molecular composites of collagens, proteoglycans and non-collagenous glycoproteins with the function of forming a supporting structure where epithelial or endothelial cells may grow. The basement membrane of skin, also known as basal laminae, provides the primary function of anchoring down the basal cells of the epidermis to the loose connective tissue of the papillary dermis. Such adherent connection is achieved by cell–matrix adhesion performed through cell adhesion molecules, having the schematic diagram shown in figure 2. The basal laminae appears in optical microscopy as a wavy band-like structure about 0.5–1 µm thick, whose very complex ultrastructure has been identified and extensively investigated in the last decades [28]. Some of its components, in fact, are targeted by autoimmune reactions, facilitating the isolation of the autoantigens that allows study of the normal structure and the functions of this membrane. In electronic microscopy, the skin basement membrane can be divided into four areas, each having a different ultrastructural organization: the hemidesmosomes/upper lamina lucida, the lower lamina lucida, the lamina densa and the sublamina densa. The basal keratinocytes of the epidermis interact with the basal laminae through hemidesmosomes, which are exactly half of the desmosomal complex observed between two adjacent keratinocytes [29]. The intermediate filaments in basal keratinocytes, made from keratins 14 and 5, insert into the electron dense plaque of the hemidesmosomes which is situated at the cell undersurface. These keratin tonofilaments attach to the Bullous Phemphigoid Antigen 1 (BPAg1), a protein with sequence homology with desmoplakin and plectin, at the intracellular side of the hemidesmosome, which is positioned just beside the plasma membrane that overlies the lamina lucida. Two other constituents of the hemidesmosomes, the transmembranous protein BPAg2 and the integrin α6β4, have extracellular domains that project down into the lamina lucida [30]. This lamina appears lucent under electron microscopy, although it is traversed by a number of filamentous, anchoring structures extending from the hemidesmosomes of basal keratinocytes directly to the lamina densa. Proteins associated with such anchoring structures include a group of heterotrimers that are synthesized and secreted by keratynocytes. In particular, laminin-5, also called kalinin, is thought to play a key role in the attachment between basal keratinocytes and the basement membrane, binding to the integrins in the upper part and forming a complex with laminin-6 in the lower part [31]. At this level is also present a salt-split zone, which approximately represents the line of separation that is possible to achieve with chemical treatment (e.g. 1 M NaCl) of an intact skin sample. The lamina densa, owing its name to its electron-dense appearance under electron microscopy, consists of a lattice network of structural proteins, having an overall thickness of about 35–45 nm [32]. The major structural protein is considered type IV collagen, which provides a scaffold for the other macromolecules, such as laminin-1, nidogen. The presence of perlecan, i.e. heparin sulphate proteoglycan, provides clusters of negatively charged groups in this area. The lamina densa is connected to the dermis through anchoring fibrils that are attached to plaques in the sublamina dense region. Anchoring fibres consist mainly of laterally associated dimers of type VII collagen, secreted by both keratinocytes and fibrolasts; minor fibrous anchoring structures are oxytalan and elaunin [33]. The amino terminal portion of type VII collagen binds type IV collagen in lamina densa or anchoring plaques in the sublamina densa. In this way, type VII collagen forms a looping structural network through which collagen fibres of the papillary dermis may interweave [34], forming the lowermost section of the adhesive structure between basal keratinocytes and the dermis. The ultrastructural properties of the components forming the basal laminae in skin suggest many functional roles for this membrane. Its major role of providing an adherent connection between the epidermis and the papillary dermis is confirmed by several blistering diseases characterized by defective protein production. Moreover, as its adhesive molecules regulate the contacts between keratinocytes and the lamina densa, the cell signalling of these molecules is involved in the regulation of cell–matrix interaction. For skin cancers, in particular, it has been observed that the first stages of melanoma progression are characterized by a loss of adhesiveness to the basal laminae, which constitutes a selective barrier against tumour invasion of the reticular dermis. In the following, a biomechanical model based on mixture theory will be presented in order to describe the initial growth phase of a malignant melanoma.

Figure 2.

Schematic of the anchoring complex between a basal keratinocyte and the ECM network of the reticular dermis. The ultrastructure of the basal laminae of human skin is depicted, with indication of the three structural layers, as visible in electronic microscopy, and of the main molecular constituents.

3. Multi-phase modelling of the radial growth phase of malignant melanoma

The development of skin cancer possesses the typical characteristics of solid tumours, starting with the generation of uncontrolled cellular proliferation, which is followed by a change in the cellular differentiation pattern that may lead to tissue invasion and, finally, to the colonization of distant organs, i.e. metastasis. In physiological terms, skin cancer progression initiates with the deregulation of homeostasis between the cellular components of the tissue. Focusing our analysis on the development of cutaneous melanoma, we are interested in all the cellular processes that alter the balance existing between the cellular components of the epidermis, and which may eventually lead to an uncontrolled proliferation of melanocytes. In normal epidermis, there exists an epidermal melanin unit denoting the symbiotic relationship in which one melanocyte has the task to carry melanosomes to approximately 36 associated keratynocites [35]. The transportation of such specialized organelles, where melanin synthesis occurs, is allowed by the cellular structure of the melanocytes, which are able to extend their dendrites to the upper parts of the epidermis. The transfer of melanin protects keratinocytes from damage by UV radiation, which in turn control differentiation and proliferation of melanocytes in the basal layer of the epidermis: in health conditions, they keep a life-long stable ratio with basal keratinocytes between 1 : 8 and 1 : 5 [18]. In order to maintain this balance, melanocytes undergo a series of cellular movements that guarantee their homeostasis in normal skin. Before performing the cellular division, they need to detach both from the basal laminae and from keratinocytes, retracting their dendrites to proliferate. Mitogenesis, which begins with upregulation of growth factors by either keratinocytes and fibroblasts, is followed by migration along the basement membrane for repositioning. Once melanocytes are anchored again to the basal laminae and to adjacent keratinocytes, they extend again their dendrites to reform the characteristic epidermal melanin unit. This normal sequence made at different steps of melanocyte proliferation is kept balanced until such homeostasis is deregulated by the development of melanoma cells and their uncontrolled growth. According to Clark et al. [36], six steps of melanoma progression can be identified: (i) common acquired/congenital naevus without cytogenetic abnormalities; (ii) melanocytic naevus with aberrant differentiation; (iii) dysplastic naevus with architectural atypia; (iv) early radial growth phase of primary melanoma; (v) advanced vertical growth phase of primary melanoma, with dermal invasion; and (vi) metastatic melanoma. Although there exist rare cases of melanoma that arise directly in the dermis (e.g. some nodular melanomas), it is widely accepted the stepwise progression mechanism where primary melanoma arise in situ in the epidermis, before invading the dermis passing through the basal laminae. Recalling that the stage of melanoma has a dramatic prognostic significance, we focus in this paper on most common forms of malignant melanoma (e.g. superficial spreading, lentigo maligna, acral lentiginous types) whose early development consists of a radial growth phase. By histological criteria, in this stage primary melanoma cells are confined primarily to the epidermis but, although invasive and common, have little or no metastatic competence. The mechanisms through which the primary melanoma escapes the control of basal keratinocytes and passes through the basal lamina are still unknown, but they initiate a phenotypic transition to the more aggressive vertical growth phase, with high competence for metastasis [37]. In the following, we will discuss how a mathematical modelling, derived in the context of the multi-phase framework of mixture theory, might help in studying the growth patterns derived by the different cell-interaction scenarios proposed for homeostasis deregulation in radial growth phase of primary melanoma. Such a pathologic growth pattern of melanocytes results from an alteration of the homeostasis regulations in the normal epidermis, which are governed either by endocrine/paracrine factors or by intercellular communications. The theoretical modelling will focus on investigating the importance of such deregulations on the surface spreading of melanoma cells over an elastic substrate, which represents the passive response of the basal laminae. Two major switching mechanisms of cell interactions that are involved in early melanoma development will be taken into account [38]. First, the loss of adhesiveness and of anchorage between melanoma cells and the basal laminae is investigated within a mixture model with three components, where the altered expression of adhesion molecules is modelled as a viscoplastic interaction force. Second, the alterations in the keratinocytes–melanocytes interactions during the radial growth phase are modelled through a four-component mixture, where a change of cadherin class during melanoma development is taken into account.

3.1. Loss of adhesion to the basal laminae components: lower epidermis as a three-component mixture

3.1.1. Definition of the model and balance equations

Let us consider the early evolution of a skin cancer, its radial growth phase, which is characterized by an increase in the population of melanoma cells at the expense of the interstitial liquid, through which the chemical factors diffuse. During the growth process, the basement membrane acts as an elastic substrate, whose adhesive properties are determined by the mechanical behaviour of the integrin-mediated hemidesmosome complex that anchors the epidermal cells. Under these assumptions, the multicellular structure of the lower epidermis will be modelled as a mixture of an interstitial liquid (subscript l), a phase of melanoma cells (subscript M) and an elastic substrate representing the basal laminae (subscript m). The basic theoretical assumption of the mixture theory is that the infinitesimal volume of the mixture is filled co-jointly by its different constituents, so that each point in the space is assumed to consist of variables and properties associated with a particle of each constituent. In other terms, mixture theory in our problem considers the skin constituents dense enough that they can be homogenized in the same volume of the mixture, each one as a continuum having its own motion and its own physical variables. At this mesoscopic scale, the evolution of the diffusing nutrients can be described through the expression of their concentration n. We neglect the mechanical interaction with keratinocytes, which will be considered in the following paragraph. Defining the jth volumetric fraction ϕj of each phase as the relevant state variables of the growth problem, the saturation condition for the mixture imposes thatEmbedded Image 3.1

The system of mass balance equations for the mixture constituents can be written asEmbedded Image 3.2where vj indicates the velocities of the generic phase j, ϱj the relative mass densities and Γj the source/sink terms of each constituent, which may depend both on nutrient concentration and on ϕk. In the following, we make the assumption that the mass densities of the constituents are all equal to the mass density of water ϱ, so that ϱj = ϱ=constant.

Recalling that the basement membrane is considered as an elastic solid medium, it is useful to rewrite its mass balance equation as a function of its relative deformation gradient tensor Fm:Embedded Image 3.3where Jm = det Fm is the Jacobian of the relative volume transformation [39].

If the mixture is closed, i.e. each constituent may be transformed into another one but the total mass of the mixture is constant, a closure of the mass balances in equation (3.2) can be imposed:Embedded Image 3.4where Vc = ∑j=M,l,m (ϕj vj) is the average velocity of the mixture. Equation (3.4) represents the classical Eulerian relation of the incompressibility constraint for the whole mixture.

The mechanistic approach in mixture theory considers the relative motion of the constituents as driven by physical forces, so that the balance of linear momentum for each phase j can be written as follows [40]:Embedded Image 3.5where fij = −fji is the interaction force between the phases i and j, p is the Lagrange multiplier arising from the incompressibility constraint in equation (3.4) and having the physical dimensions of a hydrostatic pressure for the whole mixture, σj is the excess Cauchy stress tensor of the phase j and bj is the generic body force acting on it. Considering tumour growth as a slow process which is dominated by diffusion, both the inertial effects and the body forces (e.g. chemotactic and haptotactic) in equation (3.5) can be neglected compared with viscous effects. The general balance of linear momentum for the whole mixture can be obtained summing up the three equations expressed in equation (3.5), so thatEmbedded Image 3.6where the viscous (partial) stress of the interstitial liquid is missing because a perfect fluid is assumed. Finally, the evolution equation for the nutrients must be taken into account in order to close the system of equations governing the initial growth phase of the primary melanoma. We will make use of a simple advection–diffusion model [41], where nutrients are transported by the fluid phase and are uptaken by the melanoma cells:Embedded Image 3.7where δn and kn are the absorption and diffusion coefficients, respectively. For matters of completeness, we recall that the effects of taxis-inducing chemical and molecular species and the surface effects owing to gradients in concentrations could be taken into account by using the diffuse-interface theory of mixtures [42]. Although a promising theoretical framework for bridging the gap from molecular to tissue level, such an approach will not be considered in the following, as it would be necessary to impose a number of drastic simplifications to make it tractable for meaningful simulations.

3.1.2. Constitutive equations for the partial stresses

Let us focus on the definition of the constitutive behaviour of the melanoma cells and of the basal lamina, which is their elastic substrate. For what concerns a cellular phase in a mixture, the elastic fluid model provides a simple constitutive description of the cell-to-cell interactions during growth [43]. This model gives the excess Cauchy stress in the cellular phase as a diagonal tensor whose components only depend on its volume fraction within the mixture:Embedded Image 3.8where the function Σ(ϕM) is defined through heuristic arguments. Since cells must undergo compression when densely packed, Σ(ϕM) must be a monotonically increasing positive function for ϕM greater than a threshold volume fraction ϕM0, diverging to infinity when ϕM → 1. Few mathematical functions with such characteristics have been proposed for Σ(ϕj) in the study of tumour growth [44]. The choice of an elastic fluid model is demonstrated to be suitable for simulation purposes, although reductive of a more complex rheological behaviour of living cells. Anyway, a number of viscoelastic fluid models might be chosen, whose complexity would confer a gain of stability to the solution.

From a biomechanical viewpoint, the basal laminae acts as a viscoelastic substrate upon which the anchored melanoma cells may slide. Considering that this biological membrane is made mainly of a water ground substance, in which a loose network of elastic macromolecules is immersed, it may be modelled as a nearly incompressible medium undergoing finite deformation. The choice of a solid-like model is allowed by the fact that Γm = 0, and both a reference and an actual configuration can be defined. Assuming that the growth characteristic time of the tumour is much greater than the viscoelastic time of the substrate, an hyperelastic strain energy function can be defined. In order to account for the nearly incompressibility of the body in numerical applications, let us consider the following multiplicative decomposition of the relative deformation gradient Fm:Embedded Image 3.9where Fvol = J1/3I is the volume-changing (dilatational) part, and Fiso = J−1/3 F is a volume-preserving (distortional) part; so that det Fm = det Fvol and det Fiso = 1.

Bm being the left Cauchy tensor relative to the elastic substrate, let us define a strain energy function (per unit of volume) ψm in isothermal conditions asEmbedded Image 3.10where ψvol (Jm) and ψiso(Bisom) are the dilatational and distortional strain energy functions, respectively.

From equation (3.10), it is straightforward to derive the constitutive equations for the excess Cauchy stress tensor σm:Embedded Image 3.11whereEmbedded Image 3.12andEmbedded Image 3.13A useful expression for the volumetric strain energy function has been proposed by Ogden [45] for compressible rubber-like materials (Ogden), having the following form:Embedded Image 3.14where κ is the bulk modulus in the reference configuration and β an empirical coefficient. In the case of a nearly incompressible hyperelasticity, the bulk modulus κ plays the role of a penalty coefficient in numerical applications.

For what concerns the isochoric strain energy function for the elastic substrate, few elements of skin biomechanics must be considered. The mathematical expression of this strain energy must take into account the presence of both geometrical and mechanical non-linearities. The basal laminae possesses, in fact, a wavy undulation pattern which depends on the width characteristics of the dermis and of the epidermis (very ordered in the palms of hands and feet, mostly isotropic in their soles). Furthermore, the elastic behaviour of this substrate is strongly dependent on the mechanical properties of the papillary dermis. In fact, the Young modulus of the basal laminae has been experimentally evaluated at about 200 Pa [46], being negligible compared with the elastic stiffness of human dermis which, while strongly dependent on age, ranges from 0.5 to 50 MPa [47]. Moreover, the directional orientations of collagen type I and elastin fibres in this layer confer a significant level of anisotropy to the mechanical response. In general, the isochoric strain energy function can be decoupled in an isotropic term, function of invariants I1, I2 of Bisom, and an anisotropic term, depending on the pseudo invariants Iη, (η = 4, … , 8) [48], derived from Bisom and the two structural tensors defining the principal directions of fibre-reinforcement in the papillary dermis:Embedded Image 3.15

Using the isochoric strain energy function defined in equation (3.15) allows to account for material anisotropy, where the principal directions of stress are related to the Langer's cleavage lines of human skin [49]. Such cleavage lines are used to define the direction within the human skin along which the tissue has the least flexibility, having a great importance for surgical incisions. A correlation exists between Langer lines and the anisotropic mechanical behaviour of the elastic substrate [50] opening the possibility to investigate the variability of growth patterns for different anatomical origin of the primary melanoma.

3.1.3. Interaction forces

The simplest expression of the interaction force between the interstitial liquid and the cellular and elastic phase is given proportional to the relative velocity between the two phases:Embedded Image 3.16where Πki is a positive material coefficient, which is generally dependent on the volume fraction of the phases [51]. Equation (3.16) corresponds to a Darcy-like law expressing the porosity characteristics of the intercellular fluid. In fact, using equation (3.16) in combination with the momentum balance for the liquid in equation (3.6), one obtains the following relation:Embedded Image 3.17which allows to express the velocity of the fluid as a function of the other velocities of the mixture. From equation (3.31), the global incompressibility condition for the mixture can be expressed asEmbedded Image 3.18On the other hand, the interaction force between the elastic membrane and the cellular phase is determined by the adhesion properties of the anchoring hemidesmosomes. Experimental observations have demonstrated that the adhesion of melanocytes to the basal laminae is reduced during the uncontrolled proliferation in primary melanomas [18]. In particular, a group of integrins (e.g. α3β1, α5β1, ανβ3) and the immonoglobulin superfamily (e.g. ICAM-1) are implicated in the initial stages of melanocytic tumour progression. Furthermore, Rabinovitz et al. [52] have demonstrated that the integrin-binding complex α6β4 can transmit forces without the need to engage other adhesion molecules. The connection between melanocytes and basal laminae is more likely to have a viscoplastic behaviour, and the existence of a shear yield tension σyield should be considered. In modelling terms, if the interaction force between the two phases is below this limit, the cells remain anchored to the membrane, so that vM = vm. If the yield limit is reached, a simple visco-plastic constitutive behaviour can be postulated [53], so that:Embedded Image 3.19

Finally, the expression of fm→M may be recovered from equations (3.5) and (3.6), as follows:Embedded Image 3.20where it has been assumed for the basal laminae that the viscous interaction with the interstitial liquid may be neglected if compared with the adhesive force with the cells, i.e. (ΠmMΠml).

3.1.4. Reduced equations

The saturation condition of equation (3.1), together with the balance principles in equations (3.2) and (3.5), constitute a well-posed mathematical problem. However, it has been proven that this system of equations allows one to derive directly the expression of the velocity fields in function of stresses and volume fraction if few simplifications are considered [54]. In particular, an analytical solution has been reported under the assumption thatEmbedded Image 3.21in this case; in fact, equations (3.16) and (3.17) allow one to recover vl and ∇p independently from the system of balance equations for the cell and for the basal laminae, which gives vm and vM.

In our numerical simulations, we prefer to use a different argument for getting an analytical description of the velocity fields. If we consider that the whole mixture is not subjected to external forces and to surface constraints over the boundary, we can make the hypothesis that it has a null average velocity, Vc = 0. For an initially circular tumour, this condition is forced by the fulfilment of the incompressibility condition. Then we derive:Embedded Image 3.22

Summing up the linear momentum balance equation for the melanocytes and the basement membrane, we obtainEmbedded Image 3.23where it has been assumed that ΠmlΠMl. From equation (3.20), the relation between vM and vm is given byEmbedded Image 3.24

Putting together equations (3.22), (3.23) and (3.24), the velocities of the cellular phase and of the basement membrane can be expressed as follows:Embedded Image 3.25Embedded Image 3.26if |fm→M| ≥ σyield. Equations (3.22), (3.25) and (3.26) represent the closure of the equation system, giving the analytical expression of the velocity fields of the mixture components. Such velocity fields are a function both of the volumetric fractions and of the excess stresses of the phases, thus depending on the deformation of the elastic membrane. Let us define the displacement field u in the spatial description, i.e. Eulerian form, as the difference between the current position x of the mixture at time t and the referential position X(x, t), as follows:Embedded Image 3.27

Recalling that the spatial velocity for the elastic substrate is given by vm = du/dt, we may derive a very useful relationship for calculating the displacement components from the velocity field:Embedded Image 3.28where the connection between the deformation tensor Fm and the volumetric fraction of the phase m is given by equation (3.3), where Γm = 0.

Together with constitutive relations and the mass balance equations, equations (3.22), (3.25), (3.26) and (3.28) provide a suitable mathematical framework for numerical applications.

3.2. Cadherin switching in keratinocyte–melanocyte communications: lower epidermis as a four-component mixture

Under normal conditions, the epidermal melanin unit is kept in homeostasis through the control mechanisms exerted by the adhesion molecules between melanocytes and keratinocytes. In particular, a predominant role for E-cadherin in the regulation of the proliferation and the differentiation processes of melanocytes has been demonstrated. Melanoma cells show a downregulation of E-cadherin, which plays a critical role in tumour development, allowing them to escape from the control of adjacent keratinocytes [55]. A switching of cadherin subtype is observed during the initial stages of tumour development: melanoma cells lose E-cadherin expression and instead have upregulation of N-cadherin [56]. Such a switching mechanism determines a significant change in the extent of cell-to-cell communication within the epidermis during tumour progression. In fact, melamoma cells are allowed to interact with each other, conferring them new adhesive and communication properties with the basal laminae. Moreover, melanoma cells can interact with other N-cadherin expressing cells, such as vascular endothelial cells or fibroblasts, possibly influencing tumour invasion and stromal interactions in the vertical growth phase.

The aim of this section is to model the initial growth phase of melanoma cells, where a loss of adhesion with the basal laminae is accompanied by a reduced connection with surrounding keratinocytes, as observed during in vitro experiments [57]. In the following, the radial phase of growth of melanoma in the epidermis is modelled within the mixture theory formalism. The lower epidermal layer is considered as a mixture of four different phases: an interstitial liquid (subscript l), one cellular phase made of keratinocytes (subscript K), one cellular phase made of melanoma cells (subscript M), and an elastic phase representing the basal laminae (subscript m). Recalling the arguments of the previous section, the set of equations for the mixture is given by the integrity condition, the mass and linear momentum balances and the evolution equation for the nutrients:Embedded Image 3.29with i, j = (l,M,K,m). The general balance of linear momentum for the whole mixture readsEmbedded Image 3.30while the general incompressibility condition can be written asEmbedded Image 3.31where Vc = ∑j=M,l,m,K (ϕjvj) is the average velocity of the four-component mixture.

For what concerns the interaction forces between the phases, the following simplifications can be assumed:

  • the friction between the interstitial liquid and the basal laminae is much smaller than the interaction forces between the cellular phases and the interstitial liquid (fl→mfl→M, fl→K);

  • the friction between the basal laminae and the interstitial liquid is much smaller than the interaction forces between the basal laminae and the cellular phases (fm→lfm→M, fm→K);

  • the melanoma cell phase can be modelled having a liquid-like behaviour (fM→j = ΠMj · (vjvM)), as growing melanocytes have been shown to loose adhesion to the basal laminae and to have reduced anchoring activity with adjacent keratinocytes;

  • the complex keratinocytes-basal laminae determine a viscoplastic-like behaviour for fK→m, as such a hemidesmosome-mediated link has a binding functionality, which is not altered in this stage of melanoma development.

Under these assumptions, the equations for the balance of linear momentum can be rewritten as:Embedded Image 3.32Embedded Image 3.33Embedded Image 3.34andEmbedded Image 3.35

With the aim to get an analytical expression for the velocity fields, as done in the previous section, we make the assumption that the whole mixture is not subjected to external forces, so that a null average velocity can be imposed. This simplification allows recovery of the velocity field of the intercellular fluid, given byEmbedded Image 3.36which allows, together with equation (3.30), one to get a sub-set of reduced equations which do not depend explicitly on p and vl. At the beginning of the tumour development, the keratinocytes are anchored to the basal laminae until the intensity of their interaction force reaches a certain yield level σK. From momentum balance equations, if vK = vm, the force fm→K can be expressed asEmbedded Image 3.37Embedded Image 3.38Embedded Image 3.39andEmbedded Image 3.40

Under the assumption that |fm→K| ≤ σK, the reduced equations of the problem express the velocity fields of the mixture as follows:Embedded Image 3.41Embedded Image 3.42withEmbedded Image 3.43Embedded Image 3.44

If |fm→K| > σK, the keratinocytes detach from the basal laminae, and they can slide upon this elastic layer exchanging a visco-plastic interaction force with it. From a mathematical point of view, this problem is analogous to three components, but further considerations need to be done. In particular, the value of |fm→K| is dependent on fm→M, so that the expression of |fm→K| is dependent explicitly on the velocities of the cellular phases. Thus, the easiest modelling for the keratinocyte movement along the basal laminae is to introduce a threshold value over which cells detach from such an elastic substrate and start sliding with a constant value of ΠKm.

From momentum balance equations, if |fm→K| > σK the velocities of the constituents can be expressed as:Embedded Image 3.45with:Embedded Image 3.46Embedded Image 3.47Embedded Image 3.48Embedded Image 3.49andEmbedded Image 3.50Embedded Image 3.51withEmbedded Image 3.52Embedded Image 3.53Embedded Image 3.54andEmbedded Image 3.55

Equations (3.36), (3.41), (3.42), (3.45), (3.50) and (3.51) constitute the set of the reduced equations for the four-component mixture, where the melanocytes are allowed to slide through the cellular components of the basal epidermis, and a visco-plastic behaviour is considered for the adhesion between the keratinocytes and the basal laminae.

4. The radial growth phase of primary melanomas: numerical simulations

The radial growth of primary melanomas is characterized by a surface spreading of melanoma cells at the level of the basal epidermis. In the following, we will model this stage of tumour development through a series of two-dimensional numerical simulations of the mixture theories presented in the previous section. Such an important assumption is justified both from structural considerations on the biology of the lower epidermis and from histological findings on melanoma cell distribution. From a structural viewpoint, if on the upper side water loss is avoided by the lipid granules of the lower stratum corneum, the presence of proteoglycans in the lamina densa of the basal laminae forms a negatively charged barrier, which repels water. Moreover, the overall architecture of a malignant melanoma generally possesses a necrotic centre and a peripheral vital part [58]. The histological formation of a necrotic core is compatible with a diffusion of the nutrients which arises uniquely from the horizontal contour of the cancerous mass, and not from the vertical boundaries.

In the following, we will consider that the melanoma cells have a volume fraction which is confined within a material interface xM(t), which determines a free-boundary problem over the time t. Such a tumour interface moves over time according to the following boundary condition [41]:Embedded Image 4.1where n is the unit vector which is normal to the interface in xM. The other boundary conditions at the interface are given by classical arguments of continuum mechanics, which impose the continuity of the stress and of the nutrient flux across the boundary of the tumour.

The coupling between the mass increase in the melanoma cellular phase and the concentration n of the advected/diffused nutrients has been described by the following relation [43]:Embedded Image 4.2where next is the nutrient concentration outside the tumour, while γ and δ are the material parameters determining the production and the apoptosis of melanoma cells, respectively. In equation (4.2), the existence of a nutrient concentration nmin, which is the inferior limit for growth saturation to occur, i.e. it has been postulated that cellular production is inhibited for n < nmin < next. The material parameter μ determines a coupling between stress and growth, the latter decreasing with increasing compressive stress on melanoma cells. The mathematical expression for the compressive stress Σ(ϕj) in the elastic fluid model for the cellular component j is the following:Embedded Image 4.3which is a continuous and derivable function respecting the cited requirements of monotonicity and divergence for ϕj → 1. The last term on the right-hand side of equation (4.3) is a smoothe step function centred in ϕj0, with α > 1.

For what concerns the elastic behaviour of the basal laminae, we consider the Simo & Miehe [59] model for the volumetric strain energy function, which consists in putting β = −2 in equation (3.14), where κ plays the role of a penalty coefficient. For what concerns the isotropic term of the isochoric strain energy function, a neo-Hookean constitutive model is assumed:Embedded Image 4.4where c1 is a stress-like parameter. Defining the unit vector a along the principal direction of anisotropy, the last term in equation (3.15) is assumed to have an exponential stiffening [60], as follows:Embedded Image 4.5which is a typical stiffening behaviour used for fibre-reinforcement in a preferred direction.

4.1. Definition of dimensionless variables

In order to run numerical simulations for evaluating the radial growth patterns of the theoretical models over a plane (x,y), the reduced equations for the defined mixtures need to be put in function of dimensionless variables. The time must be scaled with the characteristic time of nutrient absorption, the space with the defined diffusion length-scale, and the nutrient concentration with its value outside the tumour border. Indicating with an upperscore the dimensionless variable, we can define the following relations:Embedded Image 4.6where the same kind of notation will be used for the operator ∇̄, which will be referred to the dimensionless variables and ȳ.

In order to solve numerically the theoretical model, we still need to compute the stress components of the nearly incompressible substrate, which are given by the combination of equations (3.11), (3.12), (3.13) with equations (3.14), (4.4) and (4.5), from the left Cauchy strain tensor Bm of the elastic phase.

The numerical solution of the proposed mixture models for an initially round tumour are calculated on a two-dimensional grid using Matlab software, with the aim to evaluate the influence of the material parameters on the pattern of surface spreading of the melanoma cells.

4.2. Numerical simulations of the mixture model with three components

Numerical simulations of the mixture model with three components have been performed considering that melanoma cells can grow only at the expense of interstitial liquid, so that ΓM = −Γl. The reduced equations for the velocity fields of the mixture phases have been computed in their dimensionless form, as follows:Embedded Image 4.7

In the case of |f̄mM| ≥ σ̄yield, the dimensionless velocity fields are given by:Embedded Image 4.8where σ̄m = σm/c1, σ̄M = σM/χM andEmbedded Image 4.9

The mass balance and the evolution equations for the nutrients may be rewritten asEmbedded Image 4.10andEmbedded Image 4.11withEmbedded Image 4.12where γ̄ = (γnext)/δn, μ̄ = μχM, Σ̄(ϕM) = Σ(ϕM)/χM and δ̄ = δ/δn.

The evolution equation for the interface of the melanoma phase is given in dimensionless form as:Embedded Image 4.13while at the borders of the numerical grid, which are considered sufficiently far away from the growing tumour, we impose that = 1 and ∇ = 0.

Numerical simulations have been performed over a square two-dimensional grid, using the constitutive parameters in table 1. The evolution over time of an initially round tumour mass, modelled as an axisymmetric distribution of the volume fraction of melanoma cells, is depicted in figure 3. The simulation evidences an initial growth phase with increase in volume fraction for melanoma cells, followed by a saturation of nutrients inside the mass, which causes the appearance of a necrotic core. At later stages, the tumour mass has a peripheral vital part that seems to move as a travelling wave, as visible from figure 4, consuming the nutrients which diffuse across the border. Furthermore, the results of this simulation report a sort of contour instability for the tumour mass for (t) = |M (t)|≫ (t = 0) = 0. The occurrence of a shape instability for the proposed multi-phase model of the radial growth phase of primary melanoma will be examined in the next section from a theoretical perspective. In the following, we aim at examining the effect of the biomechanical parameters of the model on the growth properties of the round tumour mass. Examining the evolution of the volume fraction of the basal laminae during melanoma progression in figure 5, we observe that the presence of a yield stress provokes a spatial concentration of the sliding movements between the phases inside a narrow annulus at the interface of the tumour mass. In figure 6, the time-evolution of the round tumour contour is depicted for different values both of the elastic modulus of the basal laminae and of the friction coefficient between melanoma cells and such elastic substrate. In both cases the range of variability of the growth velocity is very limited, and the diffusion of nutrients seems to dominate elasticity as a mechanism of growth control. On the other hand, the loss of adhesiveness of melanoma cells on the basal laminae seems to be a crucial mechanism for the initial stages of melanoma progression. As shown in figure 7, the change of adhesive properties of the anchoring complex between melanoma cells and basement membrane, modelled as a decreasing yield stress threshold in the visco-plastic behaviour, is able to determine a relevant variation in the growth process.

View this table:
Table 1.

Numerical values of the dimensionless parameters used in the simulation of the three-component multi-phase model of melanoma growth.

Figure 3.

Evolution of the volume fraction of melanoma cells ϕM over the dimensionless time in a two-dimensional grid (260 × 260). The material parameters used for this numerical simulation are Ā = = 1, Ā1 = 1 = 1.11, 0 = 6, κ/c1 = 20 and σ̄yield = 0.

Figure 4.

Spatial evolution over of the volume fraction ϕ̄M of the melanoma cells (a) and of the nutrient concentration (b), shown at different dimensionless times . The depicted values of the two-dimensional simulation refer to an axial symmetric initial distribution of melanoma cells, and are shown in a mid-section over the axis. The initial growth phases are depicted at times = 0, 1, 2, 3 (dashed lines), while later evolution is shown at time intervals = 5, 10, 15 (solid lines) and = 20 (bold line). The material parameters used for this numerical simulation are Ā = = 1, Ā1 = 1 = 1.11, κ/c1 = 20 and σ̄yield = 0.

Figure 5.

Evolution of the volume fraction ϕ̄m of the basal laminae for σ̄yield = 0 (a) and for σ̄yield = 0.06 (b), over the dimensionless axis. The volume fractions are depicted at times = 0,0.375, 3.75, 7.5, 11.25, 15. The material parameters used for this numerical simulation are Ā = = Ā1 = 1 = 1, and κ/c1 = 20.

Figure 6.

Evolution of the ratio /0 for the tumour border over the dimensionless time . The variability of such a ratio is shown in function of c1/χM ((a) ΠmM/ΠlM = 1, = 1 = 1), and of ΠmM/ΠlM ((b) c1/χM = 1, Ā = = 1). The material parameters used for this numerical simulation are σ̄yield = 0, and κ/c1 = 20 ((a) solid orange line, c1 = χM = 1; solid blue line, c1 = χM = 3; solid brown line, c1 = χM = 5; solid black line, c1 = χM = 7; solid pink line, c1 = χM = 9. (b) Solid orange line, ΠmM/ΠlM = 0.1; solid blue line, ΠmM/ΠlM = 0.3; solid brown line, ΠmM/ΠlM = 0.5; solid black line, ΠmM/ΠlM = 0.7, solid pink line, ΠmM/ΠlM = 0.9.

Figure 7.

Evolution of the ratio /0 for the tumour border over the dimensionless time . The variability of such a ratio is shown in function of σ̄yield ((a) 0 = 8: solid orange line, σ̄yield = 0; solid blue line, σ̄yield = 0.02; solid brown line, σ̄yield = 0.04; solid black line, σ̄yield = 0.06; solid pink line, σ̄yield = 0.08), and of 0 ((b) σ̄yield = 0: solid orange line, 0 = 8; solid blue line, σ̄yield = 10; solid brown line, σ̄yield = 12; solid black line, σ̄yield = 14). The material parameters used for this numerical simulation are Ā = Ā1 = = 1 = 1, and κ/c1 = 20.

4.3. Numerical simulations of the mixture model with four components

Following the variable scaling in equation (4.6), we need to write in dimensionless form the reduced model, expressed by equations (3.36), (3.41), (3.42), (3.45), (3.50) and (3.51), in order to solve the problem of the growth of the melanoma phase in a four-component mixture. In the case that |m → K| ≤ σ̄K, the reduced equations for the mixture may be rewritten asEmbedded Image 4.14andEmbedded Image 4.15with M* = [ΠlK (ΠKM + ΠMm + ΠlM) + ΠlM (ΠKM + ΠMm)], andEmbedded Image 4.16

In the case of |m → K| > σ̄K, the reduced equations for the mixture take the following form:Embedded Image 4.17Embedded Image 4.18andEmbedded Image 4.19whereEmbedded Image 4.20

Once defined the dimensionless velocity fields from the reduced equations, it is possible to derive a numerical solution of the theoretical model, whose mass balance equations are given byEmbedded Image 4.21andEmbedded Image 4.22and the evolution of the nutrient concentration has the same form as in equation (4.11). For what concerns the growth terms, we make the hypothesis that ΓM = −ΓK, in order to create a sharp front of melanoma cells which grow at the expense of the surrounding keratinocytes.

Numerical simulations for the four-component mixture model have been performed over a square two-dimensional grid, using the constitutive parameters in table 2. The growth pattern of an initially round tumour mass is depicted in figure 8, where the volume fractions of melanoma cells and of keratinocytes are shown at different stages of surface spreading. As for the numerical results of the three-component model, this simulation evidences the formation of a necrotic core after saturation of nutrients, having a peripheral vital part moving like a travelling wave, as shown in figure 9. While the elasticity of the basal substrate is found not to be a determinant for imposing the velocity of the tumour interface (figure 10a), the cadherin switch in melanoma cells, and the consequent reduced adhesion to adjacent keratinocytes, is found to be an effective mechanism for speeding up the radial growth of the primary melanoma. As shown in figure 10b, decreasing the friction coefficient between melanoma cells and keratinocytes causes a dramatic increase in the velocity of the melanoma front. Finally, the latter is found independent on the initial size of the tumour mass, as depicted in figure 11.

View this table:
Table 2.

Numerical values of the dimensionless parameters used in the simulation of the four-component multi-phase model of melanoma growth.

Figure 8.

Evolution of the volume fraction of melanoma cells ϕM over the dimensionless time in a two-dimensional grid (160 × 160). The material parameters used for this numerical simulation are c1/χM = c1/χK = 1, 0 = 6, κ/c1 = 20, and ΠKM = 100ΠMl = 100ΠKl = 50ΠMm = 25ΠKm.

Figure 9.

Spatial evolution over of the volume fraction ϕ̄M of the melanoma cells (a) and of the keratinocytes ϕ̄K (b) shown at different times . The depicted values of the two-dimensional simulation (grid 120 × 120) refer to an axial symmetric initial distribution of melanoma cells, and are shown in a section over the axis. The growth phases are depicted at times = 0 (bold, dark red), 4.6875 (dashed, red), 37.5 (bold, red), 56.25 (dashed, purple), 93.75 (bold, purple),131.25 (dashed, blue), 168.75 (dashed, bold, blue), 206.25 (bold, blue). The material parameters used for this numerical simulation are c1/χM = c1/χK = 1, κ/c1 = 20, and ΠKM = 20ΠMl = 20ΠKl = 10ΠMm = 5ΠKm.

Figure 10.

Evolution of the ratio /0 for the tumour border over the dimensionless time . The variability of such a ratio is shown in function of c1/χM ((a) ΠKM = 10ΠMm: solid orange line, C1/χM = 1; solid blue line, C1/χM = 10; solid brown line, C1/χM = 102; solid black line, C1/χM = 103; solid pink line, C1/χM = 5·103; solid yellow line, C1/χM = 104), and of ΠKM/ΠmM ((b) c1/χM = 1: solid orange line, ΠKM/ΠmM = 10; solid blue line, ΠKM/ΠmM = 15; solid brown line, ΠKM/ΠmM = 20; solid black line, ΠKM/ΠmM = 25). The material parameters used for this numerical simulation are 4ΠMl = 4ΠKl = 2ΠMm = ΠKm, σ̄yield = 1000, 0 = 9 and κ/c1 = 20.

Figure 11.

Evolution of the ratio /0 for the tumour border over the dimensionless time . The variability of such a ratio is shown in function of 0. The material parameters used for this numerical simulation are c1/χM = c1/χK = 1, κ/c1 = 20, and ΠKM = 20ΠMl = 20ΠKl = 10ΠMm = 5ΠKm (solid orange line, 0 = 6; solid blue line, 0 = 9; solid brown line, 0 = 13; solid black line, 0 = 16).

5. Linear stability analysis of the growing tumour front

In the previous sections we have proposed a theoretical model in the multi-phase framework, and we have presented some numerical results for the simulation of the radial growth phase of primary melanoma. In the following, we are interested in performing a linear stability analysis in order to determine the stability features of a small perturbation at the interface of the tumour mass.

As the elasticity of the basal membrane is found not to influence the growth characteristic of the melanoma ensemble, let us consider a simplified model of primary melanoma as a mixture of a melanoma cell phase M and an interstitial liquid l carrying the nutrients n. The main equations of such a model are the following:Embedded Image 5.1andEmbedded Image 5.2indicating the mass balance and the reaction–diffusion equations, respectively. We recall that ϕM + ϕl = 1, so that from equation (3.22) it can be stated thatEmbedded Image 5.3where vi indicate the velocity fields of the two phases.

The use of such constitutive assumptions for the mixture allow to derive the following relation typical of porous materials:Embedded Image 5.4where K(ϕM) is a porosity coefficient that allows to recover Darcy law from such modelling framework.

The goal of this section is to analyse the stability of a growing tumour front described by equations (5.1) and (5.2). Let us consider, for matter of simplicity, a one-dimensional tumour mass that is moving in the x-direction with a certain velocity S, having in mind to study the stability features after perturbations over the transversal direction y. Let us assume that the coupling tumour–nutrient involves the following expression of the mass transformation function Γ(ϕM, n):Embedded Image 5.5where γ and δ are the rates of tumour cell production and apoptosis, respectively, and next is the nutrient concentration outside the tumour mass. Also for simplicity, we consider travelling wave solutions arising from −∞ and shape invariant in the moving frame. We fix the tumour front at xint = L(t) = St. Such solutions mimic our numerical simulations when the tumour has reached its regime of constant growth rate. Note that a true travelling-wave solution exhibits a necrotic core ϕM = 0 at −∞ and that a necessary condition for the existence of such solution is then ∫−∞L(t) Γ(ϕM, n)dx = 0.

Let us rewrite the governing equations in terms of the dimensionless variables:Embedded Image 5.6

The aim of this section is to study the stability of the tumour front int by assuming a perturbation of its position given in the moving frame byEmbedded Image 5.7

We can express the variations on ϕM and in function of a moving frame Z = (), as follows:Embedded Image 5.8andEmbedded Image 5.9where ϕM0 (Z) and 0 (Z) are the solutions of the unperturbed travelling wave solution moving along x.

Defining F(Z) = K(∂Σ/∂ϕ)f(Z), with all the quantities dependent on ϕ being taken for ϕ = ϕM0, a potential function W for the perturbed cellular velocity can be expressed asEmbedded Image 5.10so thatEmbedded Image 5.11with 0 (Z) = −K(ϕM0)(∂ Σ(ϕM0))/(∂ϕM0)ϕM0′(Z), where ′ will indicate the derivative with respect to Z, the change of coordinates allowing to derive a system of ODEs in Z from a system of ODEs in (, ȳ, ).

Substituting such perturbations into equation (5.1), one obtains the following expressions at the first order:Embedded Image 5.12andEmbedded Image 5.13where we recall that = χM/(κnΠlM), γ̄ = γnext/δn and δ̄ = δ/δn.

If < 1, then diffusion is the dominant process involved in this kind of reaction–diffusion equations, so that we are allowed to neglect advection in equation (5.2). In this case, the following expressions hold:Embedded Image 5.14andEmbedded Image 5.15One boundary condition at the interface Zint imposes the continuity of the nutrient concentration, being (0) = 1, and a stress-free tumour contour, being ϕM (0) = ϕMext (0) = constant. Such boundary conditions impose for arbitrary ȳ the following relations:Embedded Image 5.16so thatEmbedded Image 5.17The other boundary condition at the interface of the tumour Zint, expressed by equation (4.1), can be written asEmbedded Image 5.18which isEmbedded Image 5.19

If λ = κ= 0, it can be easily shown that a solution of the problem is given by f0 (Z) = −ϕM0′ (Z) and g0 (Z) = −0′ (Z). It is important to remind that such a solution (a simple translation) respects all the boundary conditions at the interface Zint, being (ϕM0″ (0) + f0 (0)) = 0.

Let us consider the eigenvalue problem asymptotically for small values of κ, with λ = λ1·κ2 + O(κ4), assuming the following asymptotic expressions for f(Z) and g(z) [61]:Embedded Image 5.20

Using equations (5.12) and (5.14), the expressions in equations (5.13) and (5.15) can be rewritten asEmbedded Image 5.21Embedded Image 5.22The interface conditions at Zint impose that f1(0) = 0 and · F1′(0) = −λ1.

In order to make few assumptions, let us consider a few properties of the unperturbed solution. First, if the tumour is growing, than we have ϕM0′ (0) < 0. Integrating equation (5.12), the following relation holds:Embedded Image 5.23where ∫0x Γ (ϕM0, n0)dx is a positive function having a maximum at x̃ = x0. Considering that a necrotic core has appeared (satisfying the necessary condition that δγ), then ϕM0 (Z) has a global maximum close to Zint, and it is almost constant when Z → −∞ (being a very slow monotonically increasing function from equation (5.23)). According to this consideration, if we consider a sufficiently large growing mass we may consider that the perturbed solution decays at Z = −∞, so that it could be treated as another interface of the tumour mass. Under this assumption, g1 (−∞) = g1′(−∞) = 0 as well as f1(−∞) = f1′ (−∞) = 0.

Integrating equation (5.22) from 0 to ∞ and recalling that f(Z) and ϕM0 (Z) are defined for −∞ < Z ≤ 0, we obtain:Embedded Image 5.24Substituting equation (5.24) into equation (5.21), we have an analytical expression for the eigenvalue λ1:Embedded Image 5.25

The analysis of the asymptotic behaviour of equations (5.12,5.13) imposes that ϕM0 (−∞) = 0, indicating the formation of a necrotic core inside the tumour mass. Assuming that δ̄ is negligible, and recalling that instability occurs if λ1 ≥ 0, from equation (5.25) the following necessary condition can be derived:Embedded Image 5.26which is related to the constitutive equation for the cellular elastic stress as follows:Embedded Image 5.27Let us consider, as an example, the following constitutive equations:Embedded Image 5.28where K0 > 0 is a dimensionless permeability coefficient of a Darcy-like model, and ℵ is a positive definite constitutive variable (ℵ ≥ 1). The constitutive relation for the elastic stress in equation (5.28) describes an attractive force between cells for 0 < ϕ < ϕext, and a repulsive interaction for ϕ > ϕext, diverging at infinity for ϕ → 1, as shown in figure 12a. Using the constitutive assumptions in equations (5.28), the necessary condition for long wavelength instability in equation (5.26) can be simplified to the following expression:Embedded Image 5.29

Figure 12.

Different shapes of the elastic stress function Σ(Φ), shown for ϕext = 0.8 (a), and instability diagram in function of ϕext (b) for different values of the constitutive parameter ℵ (ℵ = 2, solid line: ℵ = 3, dashed line; ℵ = 5, dot-dashed line; ℵ = 10, x-marked line).

The instability diagram derived from equation (5.29) is depicted in figure 12b as a function of ϕMext for different values of the constitutive parameter ℵ. Note that the stress-free volumetric fraction ϕMext has been measured from experimental tests on transvascular and interstitial fluid movements in solid tumours [62].

In conclusion, if the condition in equations (5.27) is satisfied, a long wavelength perturbation of the tumour front appears unstable, causing the propagation of a front instability in the tumour mass, which would otherwise tend towards a stationary state. To the authors' knowledge, this is the first analytical demonstration of a contour instability during the planar growth of a tumour mass. Such an unstable front may occur only after the formation of a necrotic core, depending both on the elastic and on the diffusion properties of the tumour cells, as well as on the growth/absorption properties of the diffusing nutrients.

6. Discussion

The authors have proposed a biomechanical analysis of the radial growth phase of primary melanoma, in order to investigate the role of intercellular and stromal connections on the initial stages of tumour progression within epidermis. After a preliminary introduction to skin biology and of the structural properties of the epidermal–dermal junction, the radial growth phase of a primary melanoma has been modelled within the multi-phase theory of mixtures. The purpose of this approach aimed at taking into account the effects of those cellular processes that have been experimentally observed to disrupt homeostasis in normal epidermis, and to promote skin cancer invasion. A three-component mixture model has been proposed in order to represent the loss of adhesiveness of melanoma cells with respect to the basal laminae. The elastic substrate has been considered as a quasi-incompressible soft material, whose anisotropic constitutive energy represents the collagen fibre orientation in the reticular dermis. Taking into account the passive behaviour of a hyperelastic substrate is an advancement with respect to existing multi-phase models of tumour growth. A viscoplastic constitutive behaviour for the interacting forces between the melanoma cells and the basal laminae has been considered. This choice is motivated by the reported alteration of the binding functionalities of the adherent connections in the epidermal–dermal junction [63]. Although such a homeostasis deregulation can occur at different levels and through different mechanisms, the loss of adhesiveness to components of basal laminae is regarded as one of the main events determining the initial phase of melanoma progression [64]. Using an elastic fluid model for the melanoma cells and an advection–diffusion evolution equation for the nutrient concentration, a set of reduced equation for the velocity fields of the phases have been derived. After having rewritten the reduced equations in their dimensionless form, two-dimensional simulations have been performed to mimic the surface spreading of an initially round tumour mass. The numerical results are collected in figures 37. The initial growth of the simulated tumour mass is characterized by the formation of a necrotic central core, inside which nutrients become saturated. A peripheral vital annulus of melanoma cells expands with a front velocity that depends mostly on the viscoplastic characteristic of the adherent connections with the basal laminae. In particular, the presence of a yield value for this interaction force provokes a progressive sliding of melanoma cells which is concentrated in proximity of the tumour interface. This result is in accordance with the experimental observation that active cellular connections (e.g. nidogen, integrins, collagen type IV) are mostly concentrated at the periphery of the spreading tumour [65]. The loss of adhesiveness to the basal laminae, which is modelled as a decreased threshold value for this yield stress, causes in numerical simulation a relevant increase in the velocity of the tumour front. This quantitative finding provides an important validation of the biomechanical role played by the deregulation mechanisms of adherent junctions for promoting tumour initiation [66].

The role of keratinocytes in controlling the abnormal proliferation of melanoma cells during their radial growth phase has been studied through the definition a four-component mixture model. The theoretical framework has been presented as a extension of the three-component mixture, and the generalized reduced equations are derived using a viscoplastic behaviour for the interaction between keratinocytes and basal laminae. Proliferating melanoma cells are considered to have a fluid-like reduced adhesion with the basal laminae, whereas a selective adhesion to keratinocytes is maintained [67]. Once derived the reduced equations in their dimensionless form, two-dimensional simulations have been performed, whose numerical results are collected in figures 811. In particular, the essential role of keratinocytes in the regulation of melanocyte growth and differentiation has been demonstrated. The interaction force between keratinocytes and melanocytes is found to control the front velocity of the spreading melanoma. Downregulation of E-cadherin, a critical molecule for the intercellular communications [68], is modelled as a reduced friction coefficient between the two cellular phases, allowing melanoma cells to escape from adjacent keratinocytes.

In the last section, the authors performed a linear stability analysis of the tumour interface. As the elastic properties of the basement membrane in numerical simulations do not seem to influence the growth patterns of the tumour mass, a simplified two-phase model has been taken into account. Neglecting advection compared with diffusion as a transport mechanism, the authors provided important analytical results: a large tumour mass may become unstable to long wavelength perturbations of its expanding interface. In particular, the interface can become unstable only when a necrotic core has been formed, and the volume fraction inside the tumour is smaller than its interface value. When a melanoma lesion border is studied carefully, two types of irregularities may be noticed, referred to as texture and structure [69]. Texture irregularities are the fine variations along the lesion border: their measurement is more subject to noise from the hardware imaging devices, and thus they are not suitable for an evidence tool of the disease. In contrast, structure irregularities, which are general undulations of the perimeter, may infer the abnormal important histological signs, and a higher correlation with melanomas [70]. Therefore, the reported proof of a long wavelength instability of the tumour contour allows one to describe the formation of structural irregularities with the diffusion-driven growth properties of the expanding mass. This result opens the perspective to correlate the morphological feature of a growing melanoma with the biomechanical characteristics of its micro-environment. The reported linear stability analysis for a planar front interprets as instabilities the contour undulations for the radial tumour in figures 3 and 8. Using the finite difference method on a Cartesian grid, it should be recalled that numerical artefacts could arise in such simulations. Further work is needed to assess the stability characteristics of the primary melanoma in the radial geometry in function of the physical parameters.

7. Conclusion

The authors have defined a theoretical framework for investigating the radial growth phase of cutaneous melanoma. Compared with existing mixture models, a novelty of the proposed multi-phase theory is the possibility to study the biomechanical effects on melanoma evolution both of the elastic properties of the basal laminae and of its adhesive properties with the different cell populations in the epidermis. Numerical simulations have been performed to quantify the interrelationships between the epidermal architecture and the growth, invasion and cellular/microenvironmental characteristics of melanoma cells. The results of such simulations have demonstrated that the loss of adhesiveness of the melanoma cells both to the basal laminae, caused by a deregulation mechanism of adherent junctions, and to adjacent keratynocytes, consequent to a downregulation of E-cadherin, are fundamental biomechanical features for promoting tumour initiation. A methodological novelty of the article is represented by the development of a linear stability analysis on a mixture model of tumour progression. At the authors' knowledge, the article reports the first mathematical proof of a long wavelength instability of the tumour front, whose occurrence depends both on the elastic and on the mobility diffusion properties of the tumour cells, as well as on the growth/absorption properties of the diffusing nutrients.

In conclusion, the extent of such results suggests that short-term follow-up of pigmented skin lesions may be a promising approach to recognize malignant lesion as early as possible. Theoretical modelling and simulations may provide powerful tools in clinical dermatology, with the aim to improve the efficiency of early diagnosis. Taking into account the anisotropic behaviour of skin and the principal directions of its cleavage lines, the modelling approach has the potential to move towards a personalized screening method for cutaneous melanoma. Future research activities will be focused on merging together mathematical modelling, bioimaging and biomechanical methodologies for improving the ability to differentiate skin cancers from its simulators. Ultimately, these advances may have an impact on public health, preventing unnecessary biopsies (increased specificity) while increasing the sensitivity of the existing diagnostic methods.


The authors are thankful to Dr P. Guitera (Institut Curie, Paris) and Dr H. Kittler (University of Wien) for many helpful discussions about the structural biology and the clinical diagnosis of melanomas. The contributions of O. Antonie, for the early development of the numerical code, and of C. Chatelain, for pointing out some critical points in the linear stability analysis, are acknowledged. Dr Ciarletta was funded by Marie Curie EIF grant 042069 ‘CancerBioMechanics’ under the FP6 programme.


  • Received May 26, 2010.
  • Accepted June 30, 2010.


View Abstract