Senescent fibroblasts can drive melanoma initiation and progression.
Eunjung Kim, Vito Rebecca, Inna V. Fedorenko, Jane L. Messina, Rahel Mathew, Silvya S. MariaEngler, David Basanta, Keiran S.M. Smalley, Alexander R.A. Anderson
1 Integrated Mathematical Oncology Department, H. Lee Moffitt Cancer Center & Research Institute, Tampa, FL, USA.
2 Molecular Oncology, H. Lee Moffitt Cancer Center & Research Institute, Tampa, FL, USA.
3 College of Medicine Pathology and Cell Biology, University of South Florida, Tampa, FL, USA.
4 Department of Clinical Chemistry & Toxicology, School of Pharmaceutical Sciences, University of So Paulo, So Paulo, Brazil.
Email: Eunjung.Kim@moffitt.org
Abstract
Skin is one of the largest human organ systems whose primary purpose is the protection of deeper tissues. As such, the skin must maintain a homeostatic balance in the face of many microenvironmental and genetic perturbations. At its simplest, skin homeostasis is maintained by the balance between skin cell growth and death such that skin architecture is preserved. This study presents a hybrid multiscale mathematical model of normal skin (vSkin). The model focuses on key cellular and microenvironmental variables that regulate homeostatic interactions among keratinocytes, melanocytes and fibroblasts, key components of the skin. The model recapitulates normal skin structure, and is robust enough to withstand physical as well as biochemical perturbations. Furthermore, the vSkin model revealed the important role of the skin microenvironment in melanoma initiation and progression. Our experiments showed that dermal fibroblasts, which are an important source of growth factors in the skin, adopt a phenotype that facilitates cancer cell growth and invasion when they become senescent. Based on these experimental results, we incorporated senescent fibroblasts into vSkin model and showed that senescent fibroblasts transform the skin microenvironment and subsequently change the skin architecture by enhancing the growth and invasion of normal melanocytes as well as early stage melanoma cells. These predictions are consistent with our experimental results as well as clinical observations. Our coculture experiments show that the senescent fibroblasts promote the growth and invasion of nontumorigenic melanoma cells. We also observed increased proteolytic activity in stromal fields adjacent to melanoma lesions in human histology. This leads us to the conclusion that, senescent fibroblasts create a prooncogenic environment that synergizes with mutations to drive melanoma initiation and progression and should therefore be considered as a potential future therapeutic target. This study suggests a potential link between aging in the skin microenvironment and the development of melanocytic neoplasms.
Author Summary
Skin homeostasis depends upon the complex interplay of skin cells as well as interactions between cells and the microenvironment. Here, we generated a virtual skin (vSkin) model in order to test our hypothesis that dysregulation of cellmicroenvironment interaction leads to aberrant skin structure and furthermore recapitulates pathologic conditions of the skin. To this end, we used the hybrid cellular automata method. The model couples a cellular automata that describes biological rules for cell types with partial differential equations that describe continuous microenvironmental factors. Our simulation recapitulated normal skin growth as well as selfrepair. It also showed that microenvironmental factors produced by stromal cells (fibroblasts) play an important role in maintaining normal skin homeostasis. We find that the transdifferentiation of fibroblasts due to senescence changes the skin microenvironment to drive melanoma initiation and progression. Our model gives new insights into the processes of melanoma initiation and progression, and suggests a novel strategy for treatment.
Introduction
Skin is the largest organ of the body. It is a physical barrier that limits the flow of water and electrolytes, while providing protection from ultraviolet radiation, microorganisms and toxic substances. Human skin consists of three layers, the epidermis, dermis, and subcuits. The epidermis consists of keratinocytes, melanocytes, Langerhans cells and Merkel cells. The dermis is composed of connective tissue, fibroblasts, blood vessels and lymphatics. The subcutis is made up of loose connective tissue and insulating fat cells. Keratinocytes are the main constituent of the epidermis, composing around 95% of the total cell number [1]. Within the epidermis the keratinocytes form a stratified squamous epithelium in which the keratinocytes are tightly connected to each other through desmosomes. Like other epithelial tissue, the epidermis continuously renews itself. In the epidermis, cell proliferation occurs at the basal layer followed by a stepwise upward migration and maturation of the cells. As the keratinocytes mature they lose their nuclei and turn into keratinfilled corneocytes that form a 1030 cell thick water resistant protective layer that gives the skin its critical barrier function. Protection of the keratinocytes from the DNAdamaging effects of ultraviolet radiation is mediated by the pigment producing cells, melanocytes, which locate to the basement membrane of the skin in a 1:5 ratio with the basal keratinocytes [2, 3]. Upon exposure to UV radiation, melanocytes produce melanin which is transferred to the surrounding keratinocytes via an active transport process. Once taken up into keratinocytes, melanosomes orientate themselves over the nucleus in a protective “hat” that shields the nuclei of the basal keratinocytes from UVinduced DNA damage [4]. Keratinocytes in turn regulate the growth and behavior of melanocytes through a complex network of cellcell adhesion proteins and secretory factors [3]. In addition to these interactions, the structure and strength of skin is also critically dependent upon the extracellular matrix (ECM) and the fibroblasts in the dermis. Dermal fibroblasts secrete ECM components and are responsible for mechanical strength of the dermis [5] as well as producing growth factors such as epidermal growth factors (EGF), transforming growth factor (), and basic fibroblast growth factor (bFGF) [6]. These growth factors promote growth and facilitate the constant renewal of the epidermis.
Melanoma is the most devastating form of skin cancer [7, 8], arising from the malignant transformation of melanocytes. The first phenotypic change in melanocytes is the disruption of growth controls and development of melanocytic nevi, a benign mole. At the molecular level, the growth is stimulated by the abnormal activation of the mitogenactivated protein kinase (MAPK) signaling pathway which can arise through activating mutations in the NRAS (15% of melanomas) or BRAF (50% of melanomas) [9, 10, 11, 12]. Despite BRAF mutations being important for melanoma development, the majority of benign nevi also harbor BRAF mutations [13] and rarely undergo full malignant transformation. Instead, benign nevi cease proliferation and remain quiescent for decades, entering a state of oncogeneinduced senescence (OIS) [14, 15]. This suggests that nevi must acquire additional intrinsic (molecular) or extrinsic (microenvironmental) damages to free themselves from the OIS state and develop into fullblown malignant melanoma. There is some suggestion that activation of phosphoinositide 3kinase (PI3K)/protein kinase B (PKB or AKT) due to loss of a tumorsuppressor gene, phosphatase and tensin homologue (PTEN) may be required for escape from senescence [16]. The suggested mechanism is that the inactivation of PTEN fails to attenuate the PI3K/AKT level, and consequently increases PI3K/AKT expression that facilitates cell proliferation and survival. However, the precise molecular events that underlie the transformation of nevus cells develop into melanoma are not well understood.
Cancer is typically a disease of oldage, and there is increasing evidence that senescence within the stroma, particularly in the fibroblast compartment can drive tumor development. A number of reports have shown that senescent stromal fibroblasts stimulate premalignant and malignant epithelial cells to grow in cell culture and to form tumors in mice [17, 18, 19, 20]. Mechanistically, this seems to involve the secretion of factors from senescent fibroblasts, such as Matrix metalloproteinase3 (MMP3) and Interleukin6 (IL6), that in turn remodel the microenvironment, alter epithelial differentiation, promote endothelial cell motility and stimulate cancer cell growth both in vitro and in vivo [21, 22]. It is well known that aged skin tends to harbor large numbers of senescent fibroblasts [23], and there is a direct correlation between age and the expression of (a cell cycle inhibitor) expression within the dermis [24]. In other words, concentrations of increase dramatically as tissue ages. Photodamage following exposure to UV radiation can also increase the level of senescence within the skin through the effects of oxidative damage upon the telomeres shortening [25, 26]. Although the exact role of the dermal microenvironment in melanoma development has been little studied, there is evidence that the mildly hypoxic environment of the skin contributes to melanocyte transformation [27]. It is also known that fibroblasts contribute towards the survival of earlystage melanoma cells by secreting the growth factor IGF1 and that stromal remodeling through overexpression can enhance melanoma survival in mice [28].
In order to gain a better insight into the mechanisms underlying normal skin homeostasis, and ultimately the role of the stroma in melanoma initiation and progression, we developed a virtual skin model. Several computational studies of skin have recently been developed. Deterministic ordinary differential equations have been used to model immune signaling in human skin, providing a quantitative strategy that distinguishes healthy from pathologic inflammatory responses [29]. A similar method has been employed to model mouse embryonic melanoblast proliferation dynamics [30]. Another continuous modeling framework, a system of deterministic reactiondiffusion equations, has been applied to model tumorimmune interactions [31]. A continuum fluid mechanics approach has been adopted to analyze the stability of the interface between epithelia and stroma [32]. A mixture theory has been used for structure stability analysis of melanoma growth[33] and microstructure patterning [34]. More discrete based modeling frameworks have also been applied to model skin. Agent based methods have been used to model epidermal keratinocyte behavior in normal skin [35, 36, 37], epidermal homeostasis control [38], keratinocyte originated skin disease [39], the interactions between keratinocytes and fibroblasts [40], and melanocyte distribution in epidermis [41]. However, none of these previous studies have considered normal skin homeostasis and its disruption as a route to melanoma development.
In this study, we have developed a hybrid multiscale mathematical model that simulates normal skin homeostasis. We named the model virtual skin (vSkin). We employed a hybrid cellular automata (HCA) approach [42, 43, 44, 45]. The general modeling approach is to derive a set of rules for discrete cells and couple this discrete cellular population with a suite of continuous microenvironmental factors. Previously, we have applied this method to investigate tumorstroma interactions in prostate cancer progression [46]. In this study, we developed a minimal set of cell rules for melanocytes, keratinocytes, and fibroblasts, which determine each cell’s life cycle. Then, we developed a cell interaction network that defines local interactions between cells and their microenvironment. These local interactions lead to the emergence of normal skin structure that recapitulates the in vivo structure of skin. Using our model, we observed that the vSkin model is robust enough to withstand perturbations. The vSkin model not only recovers from massive loss of its constituents but also finds an equilibrium state after superphysiological perturbations of microenvironmental factors. The proposed model also reveals the possible involvement of these microenvironmental factors in conjunction with senescent stroma in driving melanoma initiation and progression.
The paper is organized as follows. In the first section, we give an overview of the vSkin model development process. In the second section, we initially present results showing the stability of vSkin. Then, we discuss the robustness of the vSkin as well as application of the model to melanoma initiation via melanocyte mutation and fibroblast transdifferentiation. We conclude with a discussion focusing on the implications of the model predictions in terms of a novel strategy for treatment.
Material and methods
Mathematical model
We propose a systemic model of skin capable of growth control as well as selfrepair. In this section, we describe the development of the vSkin model, which represents a crosssection of skin as shown in Figure 1. We first present key microenvironmental variables involved in normal skin structure and function. Then we move on to discuss the implementation of key cell phenotypes.
Hybrid cellular automata model
In order to keep the complexity of the model to a minimum, we only consider three key cell types critical for skin function, melanocytes, keratinocytes, and fibroblasts. These cells are defined as points on a twodimensional grid that represents a crosssection of human skin (), as shown in Figure 1B. Each grid point may also be occupied by up to five microenvironmental variables, EGF, bFGF, , MMPs, and basement membrane/ECM. The discrete cells are coupled with these microenvironmental concentrations to define the hybrid cellular automata (HCA) model. The cells on the lattice influence these microenvironmental concentration fields through production and consumption, but are in turn affected by the concentrations, as they serve as regulators of cell behavior. To better understand how these cells interact with one another and their microenvironmental variables we developed a cell interaction network (Figure 2). This network shows the key interactions that are believed to be important in maintaining normal skin structure and function. The network was derived through lengthy discussion with our experimental coauthors and an extensive literature review, detailed throughout the following sections.
Microenvironmental factors
Keratinocytes, melanocytes and fibroblasts require a large number of different growth factors and adhesion signals to promote and maintain their normal growth and maintenance. However, here we restrict ourselves to either one or two chemical variables for each cell type in the model for the sake of simplicity. For a keratinocyte, we consider EGF and , because these two factors are known to be the most effective regulators for keratinocyte growth control [47, 48, 49, 50, 51]. We consider bFGF for the growth control of both melanocytes [52] and fibroblasts [53]. Since fibroblasts produce protease when they are activated [5, 6], we include MMP in the model as well. Finally, ECM concentration is incorporated, since ECM is required for structural support, and it separates the dermal and epidermal layers [1, 5] in the form of basement membrane. The dynamics of the five microenvironmental variables are defined by a system of partial differential equations that describe how each of them evolves in space and time. Since these continuous variables interact with discrete cells, for clarity we present only the discretized version of each equation.
The discretized governing equation for EGF () at a node is
(1) 
where
where a cell type can be keratinocyte (), melanocyte () or fibroblast (). In Eq (1), represents the concentration of a chemical at a node at time step (e.g.,, represents the concentration at the next time step, and represents a time step. The operator defines a twodimensional discrete Laplacian,
where is the grid size. Equation (1) defines EGF as a diffusing chemical which is increased by a rate of or when a grid node is occupied by a fibroblast () or a keratinocyte (), respectively. The concentration is decreased if a node is occupied by a keratinocyte or melanocyte () at a rate of and , respectively, due to uptake of EGF. Lastly, the concentration has a natural decay rate of per time step. Similarly, the dynamics of bFGF () is modeled as follows,
(2) 
where fibroblasts and keratinocytes produce bFGF at the rate of and , respectively, and transformed fibroblasts and transformed melanocytes produce bFGF at the rate of and , respectively. All cell types bind to bFGF at a rate of , , , , and . Finally, () is produced by fibroblasts and keratinocytes at a rate of and , respectively, diffuses at a rate , and binds to keratinocytes, melanocytes (or transformed melanocytes), fibroblasts (or abnormal fibroblasts) at rates of , , , respectively.
(3) 
ECM (E) is produced by both fibroblasts or keratinocytes depending on the local concentration of ECM. In other words, keratinocytes and fibroblasts produce ECM with a rate of and only if there is any loss compared to the initial density of epidermis () and dermis (), respectively. However, transformed fibroblasts always produce ECM with a rate of . Whenever fibroblasts are in contact with keratinocytes or melanocytes, fibroblasts can make much denser matrix which will create new basement membrane [54]. Basement membrane (BM) is defined as a continuously connected curve of grid points, each of whose ECM densities is maximal (taken to be 1 in nondimensional units). ECM is degraded by MMP at a rate of . The governing equation for ECM is
(4) 
where represents initial density of epidermis and is that of dermis. Note that when a fibroblast is in contact with either a melanocyte or a keratinocyte, we set to be its maximal value 1.0 to model fibroblast’s basement membrane generation. The function is the Heaviside step function defined as
(5) 
MMPs (P) are produced by abnormal fibroblasts, diffuse at a rate , and are depleted as they degrade the ECM (E). Note that when melanocytes are transformed, they can produce MMPs as well. In the model, we assume that transformed melanocytes produce MMPs with a rate of . The governing equation is
(6) 
All boundary conditions for all microenvironmental variables are noflux on the surface and the bottom of the domain. Periodic boundary conditions were imposed on the left and right sides of the domain.
Cell phenotypes
The behavior of each cell depends on its neighbors and the microenvironmental concentrations in which it resides. In the vSkin model, we consider the von Neumann neighborhood of range 1 (i.e., four orthogonal neighbors) for each grid point , denoted by . In vSkin, each cell has three possible phenotypes: proliferative, migratory and dead. In this section, we describe how these phenotypes are implemented. Biological rules for each cell were abstracted from a literature review and are summarized as flowcharts in Figure 3. For the sake of simplicity, each cell is allowed to execute only one phenotype, proliferative, migratory, or dead, at each cell time step (). In other words, if the cell chooses one phenotype at current time step, it exits current life cycle and waits until next time step. If a cell does not execute any of these three phenotypes, it remains as a quiescent cell.
Cellcell adhesion: In order to incorporate one of the important mechanical aspects of skin into the model, we define an adhesion level () at each node explicitly by using cellcell adhesive coefficients , where is one of its orthogonal neighbors (i.e. ). The adhesive coefficient between the same cell types is defined as 0.5. A value zero is assigned when there is no known adhesion mechanism between two cells, such as between a fibroblast and a keratinocyte. The coefficient for two melanocytes is based on current literature [55, 3]. To model Ecadherin mediated adhesion between a keratinocyte and a melanocyte, we impose a higher value (1) between these cell types. To recapitulate the anchorage of a melanocyte to the basement membrane, we also assign a higher value (1). The adhesive coefficient for all cases is summarized in Figure 4 AB. By averaging the four adhesion coefficients at each grid point, we define an adhesion level for a cell at that grid point .
(7) 
Migration: In vSkin, every cell has capacity to move to one of its orthogonal neighbors () regulated by its cellcell adhesion restriction (i.e., Pr(migration) ), where is a cellcell adhesion level obtained from equation (7). Since the migration rules are slightly different for each cell type we will give a more detailed description in the following paragraphs.
Both keratinocytes and fibroblasts move only if empty grid points exist in . We first determine the cellcell adhesion value (). If there is only one empty neighbor in the neighborhood , the cell moves to this point with probability (). If there are more empty grid points, the new position for the cell is chosen randomly from these grid points. An example of keratinocyte migration is shown in Figure 4CD.
In addition to the cellcell adhesion restriction, we use more complex rules for melanocyte migration based on the current literature [2, 3]. Melanocytes locate to the basement membrane in a 1:5 stable ratio with basal keratinocytes. Typically, one melanocyte is associated with approximately 36 keratinocyte neighbors. This symbolic structural relationship is known as the melanin unit. In the vSkin model, this melanin unit is modeled as a number of keratinocytes in the upper half of the Moore neighborhood of range for each melanocyte. The value of is set to 4 in vSkin. In vSkin, a melanocyte is moved to one of its orthogonal neighbors with regards to its cellcell adhesion level, provided that (i) the melanin unit of the melanocyte is disrupted, (ii) there is an ECM gradient towards one of its neighbors, and (iii) at least one of its neighbors is either empty or only occupied by a keratinocyte. In order to satisfy (i), a check for the number of keratinocyte neighbors is made every cell time step. Rule (ii) is implemented because melanocytes are known to have a delicate cellmatrix interaction mechanism [56, 57] and migrate via haptotaxis. Based on rule (iii), a melanocyte is allowed to be moved to another grid point in even if the point is occupied by a keratinocyte. This additional rule is implemented to recapitulate observed melanocyte high motile activity [56]. If there exist one neighbor () satisfying rule (ii) and (iii), the melanocyte is moved to the new position . If more than one such neighbors exist, the new position for the melanocyte is chosen randomly from the neighbors. An example of melanocyte migration is provided in Figure 4EF.
Proliferation: In the vSkin model, we assume that each cell has capacity for proliferation and will produce two daughter cells provided that (i) the cell specific growth factor is sufficient for cell growth and (ii) there is sufficient space surrounding the parent cell for the two daughter cells to occupy. For a keratinocyte to divide, the concentration of has to be sufficient () for its division. Then, a keratinocyte proliferates with a probability of . A melanocyte proliferates if (i) the melanin unit is disrupted and (ii) the concentration of bFGF is sufficient for cell division (). A fibroblast is allowed to proliferate if the growth factor bFGF () is greater than the threshold (). In addition to the bFGF level, we also consider a space constraint that limits the number of fibroblasts in vSkin domain. A fibroblast can divide only if at most one fibroblast neighbor exists in its neighborhood. This is not an unreasonable assumption in that the main role of fibroblasts in the model is to supply microenvironmental factors and the number of fibroblasts and the production rate of the microenvironmental factors are scalable. In other words, we can first choose a fixed number of normal fibroblasts in the model and then scale the production rates () accordingly. In order to satisfy (ii), we assumed that one daughter cell replaces the parent cell and the other daughter cell will move to any one of the parent cell’s four empty orthogonal neighbors (). If more than one of the neighboring grid points is empty, then the new cell position is chosen randomly from these points.
Death: For any cell to survive, it requires sufficient growth factors. We assume that the concentration has to drop to starvation levels, for a keratinocyte and for a melanocyte, before death can occur. For a normal fibroblast, we use a different constraint that incorporates both cell age and the concentration of bFGF. Since every fibroblast is continuously producing bFGF with the rate , it is never deprived of growth factor in normal skin conditions. We make another assumption that as a fibroblast ages, it will need more growth factors. In other words, the probability of fibroblast death is positively correlated with age but has a reciprocal relationship to the concentration of bFGF (). The space that dead cells occupy immediately becomes available to new cells.
HCA Algorithm
We now summarize the computational algorithm that integrates the discrete cells and the continuous microenvironmental variables.

Define a rectangular grid () that determines both microenvironmental concentrations and cell positions.

Set initial concentration of EGF, bFGF, and ECM at each node .

Place fibroblasts, melanocytes and keratinocytes on the grid.

Cellcell adhesion level is determined for each .

The concentration of microenvironmental factors and adhesion level are coupled with each cell to determine its phenotype.

Action associated with the determined phenotype is realized and the grid is updated accordingly.

Go to step 4.
Parameterization
This vSkin model inevitably contains a large number of parameters as it has three different cell types and five microenvironmental factors. Unfortunately, many of the model parameters were difficult, if not, impossible to obtain, and therefore need to be estimated. However, due to the interdependent nature of the variables, precise parameterization may not be as important, rather the interactions between the variables is the real driving force. Homeostasis is our primary goal and as such served as central driver for parameterization. We do not doubt that other parameter sets could achieve similar outcomes. Note that all parameters are nondimensionalized [45].
The diffusion rate of EGF () is estimated using the relation given in [58] with the molecular weight 133.07 kDa [59]. The estimated value is . The production rate of EGF by a fibroblast has not been measured. Thus, we estimate it to be a 0.5 (dimensional value = ) from multiple simulations varying this parameter within range of [0,1] to achieve normal growth (net growth ) of keratinocytes. Fibroblasts are the main growth factor supplier in skin [5, 6]. Thus, we make an assumption that a keratinocyte produces EGF with a significantly smaller rate (0.005; nondimensional). Growth factor consumption rates are difficult to measure, especially at the level of single cells. We consider a normal activity rate of growth factor receptors. For example, we use 0.01 for a normal activation rate of EGF receptors of keratinocytes, meaning that a keratinocyte is able to bind 1% of EGF at the grid point (per cell) per time step. The rate is set to be a smaller value (0.001) for melanocytes since EGF is known to stimulate keratinocyte growth mainly [47, 48]. Note that we varied these consumption rates within a range to achieve skin homeostasis. We estimate the diffusion rate of using the relation in [58] and the molecular weight 44.3112 kDa [59]. The estimated value is . The production rate of by a fibroblast is chosen to be 0.21333 from in vitro study [60]. The rate for a keratinocyte is assumed to be significantly smaller than that of a fibroblast [61]. Since is known to mainly regulate the growth of keratinocytes [49, 50, 51], we assumed that keratinocytes bind to the most. The rates for a melanocyte and a fibroblast were estimated to achieve normal skin homeostasis. The consumption rates of are estimated to be for a melanocyte, a keratinocyte, and a fibroblast, respectively. Using the same method [58] with molecular weight 17.1531 kDa, the diffusion rate is estimated to be . The production rate of bFGF by a fibroblast is also estimated from multiple simulations varying these parameters within a range of [0,1] (nondimensional). The estimated parameter is 0.05 (0.12 ). The rate for keratinocytes is assumed to be significantly smaller, since keratinocytes appear to produce significantly less bFGF than fibroblasts [62]. Since both melanocytes and fibroblasts use bFGF as a growth promoting signal and there is no evidence showing which cell binds to bFGF more, we use the same value (0.02) for the bFGF binding rate. The rate for keratinocytes is set to be significantly smaller (0.005). The diffusion rate of MMP is estimated to be in [58, 59]. The production rate of MMP by fibroblasts is varied within a range 0.024  0.12 . Estimates for the kinetic parameters and are not available since they are very difficult to obtain experimentally. Thus, we set nondimensional values that best produced homeostasis. We know that ECM density in the dermis is significantly higher than that in the epidermis (Figure 1), we therefore set the epidermal ECM density to 0.2 and the dermal ECM density is 0.7. All of the parameter values described here are summarized in Table 1.
Parameter  Description  Value  Normalized value 

Diffusion rate of  [58, 59]  0.0838  
Fibroblast production rate  0.21333 [60]  0.0873  
Keratinocyte production rate  0.004  
Keratinocyte binding rate  0.6  
Melanocyte binding rate  0.2  
Fibroblast binding rate  0.1  
Diffusion rate of EGF  [58, 59]  0.0367  
Fibroblast EGF production rate  0.5  
Keratinocyte EGF production rate  0.005  
Keratinocyte EGF binding rate  0.01  
Melanocyte EGF binding rate  0.001  
Diffusion rate of bFGF  [58, 59]  0.1708  
Fibroblast bFGF production rate  0.05  
Keratinocyte bFGF production rate  0.01  
Fibroblast bFGF binding rate  0.02  
Keratinocyte bFGF binding rate  0.005  
Melanocyte bFGF binding rate  0.02  
Decay rate of all growth factors    0.01  
Diffusion rate of MMP  [58, 59]  0.0724  
Keratinocyte ECM production rate  
Fibroblast ECM production rate  
ECM decay rate    
Epidermis ECM density    
Dermis ECM density   
The cell cycle time is set it to be a typical value . The size of the grid is set to be . The grid size is set to be (a typical cell diameter), and thus the dimensions of the skin slice we simulate are . The time step for the microenvironmental variables is set to be . We assume that all cell types have the same volume with a radius , and therefore the base number of cells is set to be . The background density of ECM is set to be [63]. The background concentration of microenvironmental variables is assumed to be [64]. We assume that the concentration of growth factors above which a cell proliferates in normal skin is set to the nondimensional threshold value 0.02 for , 0.02 for and 0.01 for . Since the concentrations below which cells die from starvation are also variables, we set a starvation threshold to be 1% (), 2% (), and 0.1% of the background microenvironmental concentration for keratinocyte, melanocyte, and abnormal fibroblast respectively.
Experimental methods
Cell culture
Melanoma cells were kindly provided by Dr Meenhard Herlyn (The Wistar Institute, Philadelphia, PA) and were cultured in RPMI medium supplemented with 5% fetal calf serum (FCS).
Coculture growth
Senescence was induced in primary fibroblasts following irradiation with 10 Gy followed by 7 days of recovery as described in Campisi and colleagues [21]. Senescent and presenescent fibroblasts were plated upon 10 cm tissue culture plates overnight before the addition of AdenoGFP tagged melanoma cells (a gift from Dr Amer Beg, Moffitt Cancer Center). Images of the growing cells were taken after 24, 48 and 72 hours with a NikonTS100 inverted fluorescence microscope.
Western blotting
Cells were plated and grown overnight. To detect ADAM9 levels, ADAM9 antibody (ab) was used for immunoblotting (IB). For all IB experiments, proteins were denatured prior to separation on 618% TrisGlycine gels. Proteins were transferred to polyvinylidene fluoride (PVDF) and blocked with 5% nonfat milk in TrisBuffered Saline Tween20 (TBST). Blots were incubated overnight in primary ab diluted according to the manufacturerÕs instructions. Secondary ab’s conjugated to the enzyme horseradish peroxidase (HRP) were used for detection by chemiluminescense.
Zymography
Cells were plated and allowed to grow overnight. To detect functional ADAM9 activity, proteins were not denatured prior to separation upon a casein gel. Cleavage of casein was imaged using an Epson V300 photo scanner.
Threedimensional spheroid assays
WM793 melanoma spheroids were prepared using the liquid overlay method [65]. Spheroids were treated with serum free RPMI, presenescent fibroblast conditioned RPMI or senescent fibroblast conditioned RPMI. Pictures of the invading spheroids were taken using a NikonTS100 inverted fluorescence microscope. Images were converted to gray scale and the spheroid frames were outlined using ImageJ (Image Processing and Analysis in Java) by modifying the threshold. The spheroid frames were than inverted and absolute intensity of the invading cells was quantified using Adobe Photoshop.
Human specimen procurement
A representative tumor sample was obtained from a patient who underwent surgical resection for metastatic melanoma and was prospectively consented and accrued to an existing melanoma tissue procurement protocol approved by the Moffitt Cancer Center Scientific Review Committee and The University of South Florida Institutional Review Board. After deparaffinization, the slides were stained for ADAM9 (Chemicon) before a thorough washing and incubation with a AlexaFluor 647 secondary antibody. The slide was also stained with DAPI to indicate the nuclei and analyzed using a Leica confocal microscope. Areas of stroma, tumor and leading edge were delineated through examination of corresponding H & E stained slide by a dermatopathologist.
Initialization
Using the HCA algorithm and parameterization described above, we first ran an initial simulation to obtain the starting configuration of the domain (Figure 1 B) for all subsequent simulations. We follow the experimental steps in the in vitro 3D organotypic skin reconstruct experiment Figure 5 AB, as if we develop a virtual organotypic skin reconstruct. This 3D organotypic culture system has served as a foundation for many basic science studies as well as a skin transplantation model, and it is known to be very stable and homeostatic (see review in [66]). The initial simulation starts with a dermal layer mixed with fibroblasts that we simulate for two weeks until dermal fibroblasts produce enough growth factors and ECM for keratinocyte and melanocyte growth. Then a mixed population of keratinocytes and melanocytes are placed on the top of this matured dermis. As soon as keratinocytes are placed, they rapidly grow until they fill almost the entire domain. Soon after the most of keratinocytes become quiescent and start to die. Finally, the population becomes stabilized. The melanocyte population also rapidly increases initially but soon finds its equilibrium. This initial simulation shows the key interactions in the model are well balanced and manifest a normal skin homeostasis. Figure 5C shows how the skin structure develops over the period of a year as well as the concentrations of the microenvironemental factors, the age of all cells and the adhesion field after one year. After a period of around 68 weeks, vSkin has already reached a homeostatic configuration. Reassuringly, the timescales of vSkin formation and the organotypic skin reconstruct are similar. We will utilize the resulting skin structure that emerges after a year as our initial condition in all subsequent simulations, since the skin structure that naturally emerged from this initial simulation is morphologically close to real skin.
Skin fitness
The degree of abnormality of vSkin was characterized by a defined metric, “skin fitness ().” We quantify ratios of each cell compartment (ratio of keratinocyte to melanocyte and that of fibroblasts to melanocytes) at each time step () and compare with the ratios of the normal state (). We also took changes of epidermal thickness into account to monitor skin fitness. Changes from time step to time step () in both compartmental ratios and epidermal thickness were weighted equally to determine the final skin fitness. The skin fitness is scored from 1 to +1, where +1 represents the maximal fitness and a normal skin state, and is defined as follows,
(8) 
where represents the sum of two ratios, the ratio of keratinocyte to melanocyte and that of fibroblast to melanocyte, represents the change from time step to , stands for a epidermal thickness, and and are the weights. In this study we choose the same weight (i.e.,).
Results
The vSkin model maintains a stable homeostasis
We already know that our vSkin model can recapitulate normal skin structure. As a typical CA method inevitably contains stochastic components, multiple realizations () were carried out for each simulation. Our simulations show that vSkin can reach a stable cell net growth rate and total cell number in the domain. Keratinocytes rapidly grow for the first two months and then reach a stable state (Figure 6 A and B), while both melanocytes and fibroblasts retain their initial population and the net growth remains stable (Figure 6 A,C, and D). These data indicate that our vSkin model rapidly achieves and maintains a stable skin homeostasis.
The vSkin model is robust to physical perturbation
Once the stability of the vSkin model has been established, we next examined its robustness. Two different types of perturbations were applied to determine the robustness of our model system. First, we tested if our vSkin model can withstand massive loss of its constituents. To this end, a triangular shaped injury, which mimics an in vivo puncture wound to the skin, was created in the center of domain. As a new space was introduced into vSkin, cells nearby the space have an increased tendency to move toward the new space due to gradient in the force distribution. In the dermal layer, fibroblasts migrate into the injury site and start to produce growth factors (EGF, , and bFGF) and extracellular matrix proteins. These newly produced growth factors promote growth of keratinocytes in epidermis, thereby starting to fill the void (healing procedure). A new basement membrane is generated by interactions between fibroblasts and either keratinocytes or melanocytes. The recovery time was defined as the time it took the wound site to form a complete basement membrane.
We also changed the size of the wound by increasing injury depth, and quantified the relationship between healing time and injury size. The recovery time from wounding tends to increase with wound size (see Figure 7 A). One hundred realizations were carried out and the mean and standard errors are reported. Interestingly, the variation among realizations increases with wound size. This is presumably due to the stochastic nature of the cell migration, thereby in basement membrane generation. In other words, the development of new basement membrane is dependent on chance interactions between fibroblasts and epidermal cells (keratinocytes or melanocytes). An example snapshot of the woundhealing process is depicted in Figure 7 B. These wound healing simulations indicate that our vSkin model is robust in the face of significant physical perturbation. Another implication is that the vSkin model can recapitulate some basic in vivo skin dynamics.
The vSkin model returns to equilibrium after biochemical perturbations
The second perturbation involves manipulation of microenvironmental factors to determine if the vSkin system is robust enough to withstand superphysiological microenvironmental changes. Specifically, each growth factor was maximized in the domain during 3, 9, 15 and 21 days. The skin fitness function in equation (8) was used to quantify abnormality of vSkin after perturbation. Note that we measure skin fitness when vSkin reached a new equilibrium.
When EGF impulses were given to vSkin, keratinocytes responded very quickly. Keratinocytes proliferated until their population saturated the domain and then became quiescent. After most of EGF has either decayed or been consumed by keratinocytes and melanocytes, the quiescent keratinocytes deprived of EGF start to die. The resulting skin histology appears to be normal as it preserves normal epidermal thickness and the normal ratio of melanocytes to keratinocytes. In Figure 8, the blue line shows what little impact the EGF impulse made on vSkin, even for the longest perturbation time.
The effects of impulses on vSkin were also negligible even though the growth of keratinocytes was inhibited by high concentration. This is presumably due to fast turn over rate of keratinocytes. The light blue line in Figure 8 depicts how the skin fitness changes over different periods of perturbation.
Upon applying bFGF impulses, melanocytes proliferated faster and started to occupy more space than keratinocytes. The overgrowth of melanocytes resulted in inhibition of growth of keratinocytes due to a lack of space. When overgrown melanocytes have consumed the excess bFGF, they started to commit apoptosis. Unlike EGF or impulses, vSkin obtained a new equilibrium state with a slight increase of melanocyte number, particularly near the basement membrane. The overall fitness score is 0.6 as depicted by the orange line in Figure 8.
Microenvironmental factors can transform normal skin to pathologic skin
Unlike EGF, , or bFGF, the overexpression of MMP had profound impact on vSkin fitness. Upon applying MMP, the extracellular matrix and the basement membrane begin to degrade. Destruction of the basement membrane increased the probability of melanocytes and keratinocytes to migrate into the dermal layer. In particular, melanocytes had a higher tendency to migrate into the dermal layer (a new microenvironment for them) where they can find more growth factors and loss of control from keratinocytes. This new microenvironment stimulates rapid growth of melanocytes. As a result of MMP overexpression, the vSkin became thinner with both keratinocytes and melanocytes moving into the dermal layer (downward). vSkin was able to recover from small changes induced by short term MMP overexpression (3 days), but fails to compensate for longterm MMP overexpression (9, 15, and 21 days), as shown in Figure 8.
Collectively, the results indicate our vSkin model is able to return to its original homeostatic state or to find a new equilibrium in the face of significant microenvironmental deviations; however, this is dependent upon the perturbations only being to those variables present in normal skin (e.g., EGF, , and bFGF). When vSkin was disturbed by a factor not typically present in normal skin, its robustness is dependent on the duration of exposure. With long term exposure of MMP, vSkin transforms to a pathologic state. This implies that the regulation of microenvironmental factors contributes to vSkin transformation from a normal to an abnormal, pathologic state.
Homeostatic disruption can lead to melanoma initiation
Our in silico perturbations have shown that vSkin is a robust system which recapitulates normal skin form and function. These experiments have also shown that bFGF and MMP are important microenvironmental regulators in skin homeostasis, particularly in melanocyte homeostasis. Melanocytes are the cell of origin for one of the most malignant forms of skin cancer, melanoma. A key to understanding melanoma initiation is determining how transformed melanocytes disrupt skin homeostasis. As we have readily shown in the previous section this homeostasis is crucially dependent on both keratinocyte regulation of melanocytes as well as microenvironmental regulation. Fibroblasts are the primary source for many of the microenvironmental factors that regulate skin homeostasis. Therefore, disruption or transformation of fibroblasts should also have a significant impact on skin homeostasis. Intriguingly, the major characteristic of transformed melanocytes that distinguishes them from normal melanocytes is their ability to produce bFGF and MMP. Similarly, we now know that senescent fibroblasts show altered patterns of ECM expression, growth factor release and protease activity [21].
To investigate this secretory phenotype further we undertook our own experimental investigation of the changes that occur in fibroblasts when they become senescent. Irradiation (10 Gy) of human primary skin fibroblasts led to their entry into a senescentlike state, as demonstrated by increased staining for galactosidase (Figure 9A,B). The onset of senescence in the fibroblasts was associated with a marked change in their mRNA profiles with increases noted in ECM constituents, growth factors and proteases (Figure 9C). Of these, ADAM9 is a protease thought to be critical for melanoma/fibroblast interactions [67]. Western blotting demonstrated that ADAM9 expression was highly upregulated in senescent fibroblasts compared to nonsenescent fibroblasts and that this was associated with increased matrix degrading capacity (Figure 9D,E). Therefore fibroblasts produce protease, increased levels of bFGF, and extracellular matrix proteins as they become senescent.
In order to explore the impact of modifying these cell types on homeostasis, we incorporated transformed melanocytes and senescent fibroblasts into the vSkin model. Figure 10 shows how these altered cell phenotypes modify the cell interaction network, now incorporating normal skin cells, transformed melanocytes and senescent fibroblasts.
Senescent fibroblasts aid melanocyte proliferation and invasion
To tease apart the relative contributions of the two abnormal cell types we initially look at each of them separately. First, we studied the role of senescent fibroblasts by examining how the degree of senescence in a subpopulation of the fibroblasts transforms normal skin structure (Figure 11). The degree of senescence specifically refers to how much bFGF and MMP the fibroblasts produce, so the most senescent phenotype produces maximal bFGF and MMP. For each parameter set (bFGF production rate, MMP production rate), we carried out multiple realizations and quantified skin fitness at each time step () using the average number of each cell compartment and the average epidermal thickness from these realizations. We observed that any degree of senescence for the fibroblasts significantly enhanced growth of melanocytes and promoted invasion into the dermis. From the skin fitness landscapes (upper panel, Figure 11) it is clear that the intensity of skin disruption is directly dependent on the degree of fibroblast senescence. When senescent fibroblasts have weak secretory phenotypes (i.e., close to normal), the fitness level of the skin they produce is almost one (i.e., normal). The skin structure, however, is not completely normal as there is some accumulation of matrix near the senescent fibroblasts (see virtual histology Figure 11, ). When senescent fibroblasts have stronger secretory phenotypes ( and ) the melanocyte population expands rapidly and skin fitness deteriorates to 0.4 (). Interestingly, the fitness subsequently increases to 0.0 as time progresses (see the difference between and ). This is due to the rapid growth of melanocytes initially , driven by the senescent fibroblasts, but this excessive growth leads to deprivation of the growth factors for both fibroblasts and keratinocyes (see Figure 10 for the interdependence of cell types and growth factors) resulting in cell death. Whilst not cancer, excessive melanocyte production is a key feature of melanoma growth. It is important to note that the melanocytes here are normal and are only responding to the abnormal microenvironmental cues being produced by the senescent fibroblasts. However, if these melanocytes were to become transformed then this would be true melanoma, therefore, senescent fibroblasts may play a critical role in creating an environment ripe for melanoma initiation.
Transformation of melanocytes drives melanoma initiation
So far, our vSkin simulations have shown that senescent fibroblasts are ideal for disrupting the normal skin environment to enhance melanocyte proliferation and invasion. However, we know that melanocytes are significantly transformed in melanoma and so we also need to consider how such mutations will affect normal skin structure. We therefore integrated phenotypes of mutant melanocytes into the vSkin model in order to study the role of melanocyte transformation as a driver of melanoma initiation. Note that we assume all fibroblasts and keratinocytes are normal for these simulations.
Just as we considered a spectrum of senescence in the fibroblast population, we examine a range of melanocyte transformation. Minimally transformed melanocytes have only lost growth inhibition and will proliferate even in the absence of bFGF. As malignancy increases, melanocytes gain the ability to produce bFGF and MMP at ever increasing levels. Therefore degree of malignancy specifically refers to how much bFGF and MMP the melanocytes produce, so the most malignant phenotype produces maximal bFGF and MMP. Interestingly, minimally mutated melanocytes (i.e. those that have only lost proliferative control) remained largely quiescent due to competition for space and growth factors as well as their inability to invade (results not shown). However, melanocytes with a greater mutational load proliferated rapidly and invaded taking over the dermal layer of the skin, as shown in Figure 12. and quantified by the change in fitness landscape at 0.25 years and 0.5 years. Since the more transformed melanocytes have autocrine signaling, they don’t suffer from growth factor deficiency even with loss of fibroblasts. With any degree of malignancy, mutant melanocytes significantly lower skin fitness and do so in a very short time frame, contrast the 0.5 year landscape of Figure 12 to the 6 year landscape of Figure 11. The effects become more dominant as the degree of transformation increases (). The mutant melanocytes disintegrate skin structure and completely take over as shown in a sample snapshot at . This vSkin transformation recapitulates a type of melanoma which progresses very rapidly.
Transformed melanocytes exploit senescent fibroblasts to drive melanoma progression
The next logical step was to integrate both transformed melanocytes and senescent fibroblasts together in the vSkin model, to investigate how interactions between these abnormal populations drive melanoma initiation and progression. In order to highlight the impact of these interactions we only considered minimally transformed melanocytes in the following simulation. It is worth noting that these minimally transformed melanocytes have only lost control over proliferation and are incapable of producing cancer in the vSkin model independently since subsequent transformation would be needed for that to occur (i.e. resulting in increased bFGF and MMP production).
Figure 13 highlights how senescent fibroblasts help transformed melanocytes proliferate more rapidly and invade into the dermis. Even when minimally transformed melanocytes are combined with the least senescent fibroblasts in vSkin we see a positive impact on the melanoma growth, with a nodule of melanocytes (Figure 13A ) forming after 6 years. Increasing the degree of senescence in the fibroblast population results in a more striking effect with the transformed melanocytes taking over the dermal layer, especially those nearby senescent fibroblasts (Figure 13B). From these results we can conclude that strong secretory phenotypes (e.g. the fibroblasts in Figure 13 ) can transform essentially normal skin (with the exception of the initial mutant melanocyte) to a pathologic state (see density graph ). Therefore, senescent fibroblasts have the potential to propagate and maintain melanoma development, provided a minimally transformed melanocyte population already exists.
In vitro experiments show senescent fibroblasts enhance melanocyte/melanoma growth and invasion
So far, we have presented vSkin model predictions that senescent fibroblasts change the skin microenvironment and produce an aberrant skin structure that results from the overgrowth and invasion of normal melanocytes and early stage melanoma cells. In this section, we compare the model predictions with an in vitro coculture assay.
We explored the prediction that senescent fibroblasts enhanced the growth and invasion of minimally transformed melanocytes in a series of in vitro fibroblast coculture models. It was noted that plating early stage WM35 melanoma cells (which are nontumorigenic in mice) on top of senescent fibroblasts enhanced their growth compared to nonsenescent fibroblasts (Figure 14 A,B). To study invasion in a tissuelike context, we next generated spheroids in which poorly invasive melanoma cells (WM793) were cocultured with either nonsenescent or senescent human skin fibroblasts in a 3D collagen gel (Figure 14 C,D). It was noted that the WM793 cells invaded the collagen minimally when grown alone and that their invasive behavior was markedly increased when cocultured with the senescent fibroblasts. Although not direct comparison with the vSkin results, our coculture cell growth assay highlights one of predictions that the vSkin model makes, which is that senescent fibroblasts promote the growth of early stage melanoma cells. The potential clinical relevance of our model predictions and experimental findings was suggested by the observation that human melanoma specimens stained positively for the protease ADAM9 at the leading edge where the melanoma cells and fibroblasts interact (Figure 15). In contrast, little ADAM9 staining was noted in the tumor core where stromal infiltration was lacking.
Discussion
The regulation of skin homeostasis involves a complex interplay between different cell types and their microenvironment. To better understand how these interactions produce skin homeostasis, we developed a mathematical model that incorporates the key cellular (keratinocytes, melanocytes, and fibroblasts) and microenvironmental (, EGF, bFGF and ECM) components of normal skin. Our virtual skin model (vSkin) successfully recapitulated normal skin structure, skin homeostasis and wound healing responses. Surprisingly the homeostasis that emerged from this model was remarkably robust, being able to recover from both physical and biochemical perturbations. However, the recovery was altered by both the perturbation type and period. With shortterm perturbations inducing increased population growth, followed by a return to stable homeostasis. More persistent microenvironmental perturbations (specifically with bFGF) led to an altered but stable homeostatic state characterized by an increased level of melanocytes. One key perturbation that the system wasn’t robust to was overproduction of MMP and even then it was only when MMP levels were elevated for significant periods of time. Although, this isn’t too surprising as MMPs are not normally produced in normal skin. These results indicate the importance of microenvironmental factors in maintaining skin homeostasis and highlight both bFGF and MMP as key factors in disrupting this homeostasis.
The dialogue between our three cell types, keratinocytes, melanocytes, and fibroblasts is primarily mediated by growth factors produced by fibroblasts and cellcell contact. One of the primary reasons for incorporating fibroblasts into vSkin was to understand their role in normal skin homeostasis, disruption and cancer initiation. An important experimental observation that motivated this line of inquiry was that when fibroblasts take on a secretory phenotype (such as that seen when undergoing senescence) they produce increased levels of bFGF and matrix degrading proteases (Figure 9), the very factors that normal skin seems to be most sensitive to. By introducing a population of senescent fibroblasts into vSkin we made the unexpected prediction that such secretory cells may play a central role in the initiation and progression of melanoma. The least senescent fibroblasts caused a temporary disruption of homeostasis but as we increased the level of senescence (i.e. increasing the levels of bFGF and MMP production) we were able to permanently disrupt the skin structure and produced melanocytic hyperplasia i.e. a mole like structure. However, if the senescent fibroblasts were to be removed, normal skin homeostasis would be reintroduced since the melanocytes are normal and only reacting to the altered microenvironment created by the secretory cells.
Our vSkin model is the first model showing all the steps in skin carcinogenesis, from normal skin followed by melanocytic hyperplasia, and ending with melanoma. The activation of stroma due to senescence induces excessive growth of normal melanocytes, resulting in the development of melanocytic hyperplasia (Figure 11). The vSkin model also recapitulated melanoma development driven only by melanocyte transformation (Figure 12). In addition, the vSkin model proposes alternative routes to melanoma initiation, by showing that melanoma can occur as a result of cooperation between mutant melanocytes and senescent fibroblasts (Figure 13). Given that most known driver mutations in melanoma are associated with entry into oncogeneinduced senescence (OIS), our model suggests the possibility that signals from the senescent stroma may aid melanocyte transformation by allowing them to exit OIS. There is already good evidence that increased PI3K/AKT signaling can overcome BRAF V600Einduced OIS in melanocytes [16] and significantly most of the growth factors and altered ECM interactions we describe here are known to exert their effects through the PI3K/AKT pathway [68].
Biochemical analysis shows that the alterations in senescent fibroblasts included increased expression of multiple growth factors and proteases. These alterations then stimulated the growth of normal melanocytes and early stage melanoma cells. The increased level of protease in senescent fibroblasts enhanced the invasiveness of late stage melanoma cells. Among many proteases, ADAM9 was found to be a putative mediator of matrix remodeling in senescent dermal fibroblasts. Importantly, the analysis of ADAM9 expression in melanoma samples from patients revealed that the expression was increased at the hosttumor interface, but not within the tumor core (Figure 15). Our findings mirror previous studies implicating a role for ADAM9 expression in melanoma cell/fibroblast interactions at the invasive front that is required for MMP1 and MMP2 mediated matrix degradation [67]. Whilst not conclusive, the ADAM9 staining implicates the stroma is producing matrix degrading enzymes. This clinical observation along with all of the experimental results we presented were consistent with our model predictions that senescent fibroblasts contribute to melanoma initiation and early melanoma progression.
Our vSkin model proposes multiple routes to melanoma initiation and progression, summarized in the upper panel of Figure 16. The obvious and rapid route shown by the blue lines is due to the accumulation of mutations in melanocytes to drive melanoma development. The alternative route involves cooperation of mutant melanocytes with senescent fibroblasts, Figure 16 (red lines). When senescent fibroblasts stimulate and maintain the excessive growth of minimally transformed melanocytes we observe the slowest melanoma growth. As senescent fibroblasts adopt stronger secretory phenotypes, melanoma progresses more rapidly.
If our predictions stand up to future validation, the fibroblasts could represent an important target for melanoma chemoprevention. The vSkin transformation to melanoma suggests a possible novel treatment, by negating the influence malignant phenotypes it may be possible to normalize the skin microenvironment and return (or maintain) skin homeostasis. Since proteases are one of key contributing factors driving melanoma initiation and progression, both from the senescent fibroblast and transformed melanocyte perspective, antiMMP therapy might be a good treatment option. Indeed, there is already evidence that mutations in MMP8, leading to an inactivation of protease activity that shifts the pattern of matrix degradation and growth factor availability, is an important event in the genesis of 23% of all melanomas [69]. As a purely hypothetical exercise we ran a suite of simulations to test this idea, almost as if we applied antiMMP cream to the vSkin model such that MMP expression was blocked completely. Specifically, the MMP level is completely knocked down to zero uniformly in the whole domain, when the size of initial melanocyte population had doubled. This knock down is maintained until the simulation is finished. Our in silico antiMMP treatment experiments show that the treatment was far more effective for the stroma dependent melanoma in that the population growth rates were significantly decreased (see Figure 16 (lower panel)).
This study emphasizes the importance of normal skin homeostasis and its disruption in melanoma initiation and progression. The possible role of MMP inhibition as a chemoprevention strategy for melanoma will require extensive future experimental validation. MMPs and other related proteases constitute a large family of enzymes with distinct pro and antioncogenic activities. The effective use of MMP blocking drugs is likely to require an exquisite sensitivity profile that current inhibitors do not yet offer. Our integrated approach highlighted two key unifying factors in disrupting skin homeostasis, MMP and bFGF production. Whether these factors were produced by transformed melanocytes or senescent fibroblasts led to similar skin degeneration. The most important implication is that we cannot only consider melanoma as a transformation of the melanocyte population but must also consider the fibroblast population. It is perhaps no accident then that melanoma really is a disease of the aged [70, 71, 72], where fibroblast senescence may naturally occur as we age and aid minor melanocyte mutations (due to a lifetime of skin sun exposure) to become full blown melanoma.
Acknowledgments
We thank Dr. David Basanta for his initial input on the cell interaction network and crucially for providing his code, used in [46], that was used as the foundation for the vSkin code.
Author Contributions
Conceived and designed the computational model: EK and AA. Contributed initial implementation of the computational model: DB. Performed computational simulations: EK. Conceived and designed the experiments: KS. Performed the experiments: VR and IF. Analyzed the data: EK, JM, KS and AA. Contributed reagents/materials: JM, RM, and SE. Wrote the paper: EK, KS and AA.
References
 1. Rook A, Burns T (2010) Rook’s textbook of dermatology. Chichester, West Sussex, UK: WileyBlackwell, 8th edition.
 2. Fitzpatrick TB, Breathnach AS (1963) The epidermal melanin unit system. Dermatol Wochenschr 147: 4819.
 3. Haass N, Herlyn M (2005) Normal human melanocyte homeostasis as a paradigm for understanding melanoma. Journal of Investigative Dermatology Symposium 10: 15363.
 4. Costin G, Hearing V (2007) Human skin pigmentation: melanocytes modulate skin color in response to stress. FASEB J 21: 97694.
 5. Sorrell JM, Caplan AI (2004) Fibroblast heterogeneity: more than skin deep. J Cell Sci 117: 66775.
 6. Kalluri R, Zeisberg M (2006) Fibroblasts in cancer. Nat Rev Cancer 6: 392401.
 7. Chin L (2003) The genetics of malignant melanoma: lessons from mouse and man. Nat Rev Cancer 3: 55970.
 8. Miller AJ, Mihm MC Jr (2006) Melanoma. N Engl J Med 355: 5165.
 9. Albino AP, Nanus DM, Mentle IR, CordonCardo C, McNutt NS, et al. (1989) Analysis of ras oncogenes in malignant melanoma and precursor lesions: correlation of point mutations with differentiation phenotype. Oncogene 4: 136374.
 10. Davies H, Bignell GR, Cox C, Stephens P, Edkins S, et al. (2002) Mutations of the BRAF gene in human cancer. Nature 417: 94954.
 11. Omholt K, Platz A, Kanter L, Ringborg U, Hansson J (2003) NRAS and BRAF mutations arise early during melanoma pathogenesis and are preserved throughout tumor progression. Clin Cancer Res 9: 64838.
 12. Houben R, VetterKauczok CS, Ortmann S, Rapp UR, Broecker EB, et al. (2008) PhosphoERK staining is a poor indicator of the mutational status of BRAF and NRAS in human melanoma. J Invest Dermatol 128: 200312.
 13. Pollock PM, Harper UL, Hansen KS, Yudt LM, Stark M, et al. (2003) High frequency of BRAF mutations in nevi. Nat Genet 33: 1920.
 14. Michaloglou C, Vredeveld LCW, Soengas MS, Denoyelle C, Kuilman T, et al. (2005) BRAF V600Eassociated senescencelike cell cycle arrest of human naevi. Nature 436: 720–4.
 15. Dhomen N, ReisFilho JS, da Rocha Dias S, Hayward R, Savage K, et al. (2009) Oncogenic BRAF induces melanocyte senescence and melanoma in mice. Cancer Cell 15: 294303.
 16. Vredeveld LCW, Possik PA, Smit MA, Meissl K, Michaloglou C, et al. (2012) Abrogation of BRAF V600Einduced senescence by PI3K pathway activation contributes to melanomagenesis. Genes Dev 26: 105569.
 17. Krtolica A, Parrinello S, Lockett S, Desprez PY, Campisi J (2001) Senescent fibroblasts promote epithelial cell growth and tumorigenesis: a link between cancer and aging. Proc Natl Acad Sci U S A 98: 120727.
 18. Krtolica A, Campisi J (2002) Cancer and aging: a model for the cancer promoting effects of the aging stroma. Int J Biochem Cell Biol 34: 140114.
 19. Bhowmick N, Neilson E, Moses H (2004) Stromal fibroblasts in cancer initiation and progression. Nature 18: 3327.
 20. Parrinello S, Coppe JP, Krtolica A, Campisi J (2005) Stromalepithelial interactions in aging and cancer: senescent fibroblasts alter epithelial cell differentiation. J Cell Sci 118: 48596.
 21. Coppé J, Patil C, Rodier F, Sun Y, Muñoz D (2008) Senescenceassociated secretory phenotypes reveal cellnonautonomous functions of oncogenic RAS and the p53 tumor suppressor. PLoS Biol 6: 285368.
 22. Laberge RM, Awad P, Campisi J, Desprez PY (2012) Epithelialmesenchymal transition induced by senescent fibroblasts. Cancer Microenviron 5: 3944.
 23. Campisi J (1998) The role of cellular senescence in skin aging. J Investig Dermatol Symp Proc 3: 15.
 24. Ressler S, Bartkova J, Niederegger H, Bartek J, ScharffetterKochanek K, et al. (2006) p16INK4A is a robust in vivo biomarker of cellular aging in human skin. Aging Cell 5: 37989.
 25. von Zglinicki T (2000) Role of oxidative stress in telomere length regulation and replicative senescence. Ann N Y Acad Sci 908: 99110.
 26. von Zglinicki T (2002) Oxidative stress shortens telomeres. Trends Biochem Sci 27: 33944.
 27. Bedogni B, Welford SM, Cassarino DS, Nickoloff BJ, Giaccia AJ, et al. (2005) The hypoxic microenvironment of the skin contributes to AKTmediated melanocyte transformation. Cancer Cell 8: 44354.
 28. Li G, Satyamoorthy K, Meier F, Berking C (2003) Function and regulation of melanomastromal fibroblast interactions: when seeds meet soil. Oncogene 22: 316271.
 29. Valeyev NV, Hundhausen C, Umezawa Y, Kotov NV, Williams G, et al. (2010) A systems model for immune cell interactions unravels the mechanism of inflammation in human skin. PLoS Comput Biol 6: e1001024.
 30. Aylaj B, Luciani F, Delmas V, Larue L, De Vuyst F (2011) Melanoblast proliferation dynamics during mouse embryonic development. Modeling and validation. J Theor Biol 276: 8698.
 31. Eikenberry S, Thalhauser C, Kuang Y (2009) Tumorimmune interaction, surgical treatment, and cancer recurrence in a mathematical model of melanoma. Plos Comput Biol 5: e1000362.
 32. Basan M, Joanny JF, Prost J, Risler T (2011) Undulation instability of epithelial tissues. Physical Review Letters 106: 158101.
 33. Ciarletta P, Foret L, Ben Amar M (2011) The radial growth phase of malignant melanoma: multiphase modelling, numerical simulations and linear stability analysis. Journal of the Royal Society Interface 8: 345–368.
 34. Chatelain C, Balois T, Ciarletta P, Ben Amar M (2011) Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture. New Journal of Physics 13: 115013.
 35. Grabe N, Neuber K (2005) A multicellular systems biology model predicts epidermal morphology, kinetics and ca2+ flow. Bioinformatics 21: 3541–7.
 36. Suetterlin T, Huber S, Dickhaus H, Grabe N (2009) Modeling multicellular behavior in epidermal tissue homeostasis via finite state machines in multiagent systems. Bioinformatics 25: 2057–2063.
 37. Adra S, Sun T, MacNeil S, Holcombe M, Smallwood R (2010) Development of a three dimensional multiscale computational model of the human epidermis. PLoS ONE 5: e8511.
 38. Schaller G, MeyerHermann M (2007) A modelling approach towards epidermal homoeostasis control. J Theor Biol 247: 55473.
 39. Grabe N, Neuber K (2007) Simulating psoriasis by altering transit amplifying cells. Bioinformatics 23: 1309–12.
 40. Sun T, McMinn P, Holcombe M, Smallwood R (2008) Agent based modelling helps in understanding the rules by which fibroblasts support keratinocyte colony formation. PLoS ONE 3: e2129.
 41. Thingnes J, Lavelle TJ, Hovig E, Omholt SW (2012) Understanding the melanocyte distribution in human epidermis: an agentbased computational model approach. PLoS One 7: e40377.
 42. Anderson ARA, Sleeman BD, Young IM, Griffiths BS (1997) Nematode movement along a chemical gradient in a structurally heterogeneous environment. 2. Theory. Fundam Appl Nematol 20: 165172.
 43. Anderson AR, Chaplain MA (1998) Continuous and discrete mathematical models of tumorinduced angiogenesis. Bull Math Biol 60: 85799.
 44. Anderson A, Pitairn A (2003) Application of the hybrid discretecontinuum technique, in polymer and cell dynamics, eds. W. Alt, M. Chaplain, M. Griebel, J. Lenz. Birkhauser : 261279.
 45. Anderson ARA (2005) A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. Math Med Biol 22: 163–86.
 46. Basanta D, Strand DW, Lukner RB, Franco OE, Cliffel DE, et al. (2009) The role of transforming growth factorbetamediated tumorstroma interactions in prostate cancer progression: an integrative approach. Cancer research 69: 7111–20.
 47. Green MR, Couchman JR (1985) Differences in human skin between the epidermal growth factor receptor distribution detected by EGF binding and monoclonal antibody recognition. J Invest Dermatol 85: 23945.
 48. Nanney LB, Stoscheck CM, Magid M, King LE Jr (1986) Altered [125I] epidermal growth factor binding and receptor distribution in psoriasis. J Invest Dermatol 86: 2605.
 49. Krane JF, Murphy DP, Carter DM, Krueger JG (1991) Synergistic effects of epidermal growth factor (EGF) and insulinlike growth factor I/somatomedin C (IGFI) on keratinocyte proliferation may be mediated by IGFI transmodulation of the EGF receptor. J Invest Dermatol 96: 41924.
 50. Hannon GJ, Beach D (1994) p15INK4B is a potential effector of TGFbetainduced cell cycle arrest. Nature 371: 25761.
 51. Ewen ME (1996) p53dependent repression of cdk4 synthesis in transforming growth factorbetainduced G1 cell cycle arrest. J Lab Clin Med 128: 35560.
 52. Shih IM, Herlyn M (1993) Role of growth factors and their receptors in the development and progression of melanoma. J Invest Dermatol 100: 196S203S.
 53. Makino T, Jinnin M, Muchemwa FC, Fukushima S, KogushiNishi H, et al. (2010) Basic fibroblast growth factor stimulates the proliferation of human dermal fibroblasts via the ERK1/2 and JNK pathways. Br J Dermatol 162: 71723.
 54. Andriani F, Margulis A, Lin N, Griffey S, Garlick JA (2003) Analysis of microenvironmental factors contributing to basement membrane assembly and normalized epidermal phenotype. J Invest Dermatol 120: 92331.
 55. Haass N, Smalley K, Li L (2005) Adhesion, migration and communication in melanocytes and melanoma. Pigment Cell Melanoma Res 18: 1509.
 56. Nikolas K Haass LL Kerian S M Smalley, Herlyn M (2005) Adhesion, migration and commnication in melanocyte and melanoma. Pigment Cell Res 18: 1509.
 57. Morelli JG, Yohn JJ, Zekman T, Norris DA (1993) Melanocyte movement in vitro: role of matrix proteins and integrin receptors. J Invest Dermatol 101: 6058.
 58. Swabb EA, Wei J, Gullino PM (1974) Diffusion and convection in normal and neoplastic tissues. Cancer research 34: 2814–22.
 59. The UCSD Signaling Gateway. URL http://www.signalinggateway.org/.
 60. Partridge M, Green MR, Langdon JD, Feldmann M (1989) Production of TGFalpha and TGFbeta by cultured keratinocytes, skin and oral squamous cell carcinomas–potential autocrine regulation of normal and malignant epithelial cell proliferation. Br J Cancer 60: 5428.
 61. Amjad SB, Carachi R, Edward M (2007) Keratinocyte regulation of TGFbeta and connective tissue growth factor expression: a role in suppression of scar tissue formation. Wound Repair Regen 15: 74855.
 62. Halaban R, Langdon R, Birchall N, Cuono C, Baird A, et al. (1988) Basic fibroblast growth factor from human keratinocytes is a natural mitogen for melanocytes. J Cell Biol 107: 16119.
 63. Kim Y, Friedman A (2010) Interaction of tumor with its microenvironment: A mathematical model. Bulletin of Mathematical Biology 72: 1029–1068.
 64. Tavakkol A, Varani J, Elder JT, Zouboulis CC (1999) Maintenance of human skin in organ culture: role for insulinlike growth factor1 receptor and epidermal growth factor receptor. Arch Dermatol Res 291: 64351.
 65. Smalley KSM, Haass NK, Brafford PA, Lioni M, Flaherty KT, et al. (2006) Multiple signaling pathways must be targeted to overcome drug resistance in cell lines derived from melanoma metastases. Mol Cancer Ther 5: 113644.
 66. Brohem CA, Cardeal LBdS, Tiago M, Soengas MS, Barros SBdM, et al. (2011) Artificial skin in perspective: concepts and applications. Pigment Cell Melanoma Res 24: 3550.
 67. Zigrino P, Nischt R, Mauch C (2011) The disintegrinlike and cysteinerich domains of ADAM9 mediate interactions between melanoma cells and fibroblasts. J Biol Chem 286: 68017.
 68. Dey JH, Bianchi F, Voshol J, Bonenfant D, Oakeley EJ, et al. (2010) Targeting fibroblast growth factor receptors blocks PI3K/AKT signaling, induces apoptosis, and impairs mammary tumor outgrowth and metastasis. Cancer Res 70: 415162.
 69. Palavalli LH, Prickett TD, Wunderlich JR, Wei X, Burrell AS, et al. (2009) Analysis of the matrix metalloproteinase family reveals that MMP8 is often mutated in melanoma. Nat Genet 41: 51820.
 70. Chao C, Martin RCG 2nd, Ross MI, Reintgen DS, Edwards MJ, et al. (2004) Correlation between prognostic factors and increasing age in melanoma. Ann Surg Oncol 11: 25964.
 71. Lachiewicz AM, Berwick M, Wiggins CL, Thomas NE (2008) Epidemiologic support for melanoma heterogeneity using the surveillance, epidemiology, and end results program. J Invest Dermatol 128: 13402.
 72. Howlader N, Noone A, Krapcho M, Neyman N, Aminou R, et al. (2011) SEER Cancer Statistics Review, 19752009 (Vintage 2009 Populations) based on November 2011 SEER data submission, posted to the SEER web site, April 2012. National Cancer Institute. Bethesda, MD.