{'fraction_non_alphanumeric': 0.0282763167533514, 'fraction_numerical': 0.010680536224452421, 'mean_word_length': 4.987893674883762, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Senescent fibroblasts can drive melanoma initiation and progression

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 enhance 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 co-culture experiments show that the senescent fibroblasts promote the growth and invasion of non-tumorigenic melanoma cells. We also observed increased proteolytic activity in stromal fields adjacent to melanoma lesions in human histology. Collectively, senescent fibroblasts create a pro-oncogenic environment that synergizes with mutations to drive melanoma initiation and progression and should therefore be considered as a potential future therapeutic target.


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 keratin-filled corneocytes that form a 10-30 cell thick water resistant protective layer that gives the skin its critical barrier function. Protection of the keratinocytes from the DNA-damaging 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 UV-induced DNA damage [4]. Keratinocytes in turn regulate the growth and behavior of melanocytes through a complex network of cell-cell 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 β (TGF β ), 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 mitogen-activated 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 oncogene-induced 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 full-blown malignant melanoma. There is some suggestion that activation of phosphoinositide 3-kinase (PI3K)/protein kinase B (PKB or AKT) due to loss of a tumor-suppressor 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 old-age, 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 metalloproteinase-3 (MMP-3) and Interleukin-6 (IL-6), 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 p16 IN K4A (a cell cycle inhibitor) expression within the dermis [24]. In other words, concentrations of p16 IN K4A increase dramatically as tissue ages. Photo-damage 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 early-stage melanoma cells by secreting the growth factor IGF-1 and that stromal remodeling through TGF β 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 reaction-diffusion equations, has been applied to model tumor-immune 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 tumor-stroma 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 super-physiological 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.

Mathematical model
We propose a systemic model of skin capable of growth control as well as self-repair. In this section, we describe the development of the vSkin model, which represents a cross-section 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 cross-section of human skin (6 mm × 1.24 mm), as shown in Figure 1B. Each grid point may also be occupied by up to five microenvironmental variables, EGF, bFGF, TGF β , 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 produce or promote inhibit Melanocyte TGFβ EGF Fibroblast Keratinocyte bFGF Figure 2. Cell interaction network for normal skin. Green lines represent promoting events, while red lines imply inhibiting effects. Note, cell coloration is used in all simulations to represent the equivalent cell. In the model, a fibroblast is assumed to produce three growth factors, Epidermal growth factor (EGF), Basic fibroblast growth factor (bFGF) and Transforming growth factor (TGF β ). EGF promotes growth of keratinocytes, whilst TGF β inhibits the growth. The growth of melanocyte is controlled by both bFGF and physical interaction (contact inhibition) with neighboring keratinocytes.
maintaining normal skin structure and function. The network was derived through lengthy discussion with our experimental co-authors 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 TGF β , 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 (G E ) at a node η ≡ (η x , η y ) ∈ Ω ⊂ Z 2 is production by fibroblasts and keratinocytes binding to EGF receptors of keratinocytes or melanocytes where X η = 1 if a node η is occupied by a cell type X at time step t, 0 otherwise, where a cell type X can be keratinocyte (K), melanocyte (M ) or fibroblast (F ). In Eq (1), [·] represents the concentration of a chemical at a node η at time step t (e.g., represents the concentration at the next time step, and δt represents a time step. The operator ∆ 2 defines a two-dimensional discrete Laplacian, where δh is the grid size. Equation (1) defines EGF as a diffusing chemical which is increased by a rate of α E or β E when a grid node η is occupied by a fibroblast (F ) or a keratinocyte (K), respectively. The concentration is decreased if a node is occupied by a keratinocyte or melanocyte (M ) at a rate of γ E and ζ E , respectively, due to uptake of EGF. Lastly, the concentration has a natural decay rate of λ per time step. Similarly, the dynamics of bFGF (G b ) is modeled as follows, production by fibroblasts, keratinocytes, abnormal fibroblasts and transformed melanocytes binding to bFGF receptors of keratinocytes, melanocytes, fibroblasts, transformed melanocytes, and abnormal fibroblasts where fibroblasts and keratinocytes produce bFGF at the rate of α b and β b , respectively, and transformed fibroblasts and transformed melanocytes produce bFGF at the rate of α s b and β s b , respectively. All cell types bind to bFGF at a rate of γ b , ζ b , ξ b , ζ s b , and ξ s b . Finally, TGF β (T β ) is produced by fibroblasts and keratinocytes at a rate of α T β and β T β , respectively, diffuses at a rate D T β , and binds to keratinocytes, melanocytes (or transformed melanocytes), fibroblasts (or abnormal fibroblasts) at rates of γ T β , ζ T β , ξ T β , respectively.
production by fibroblasts and keratinocytes binding to TGF β receptors of keratinocytes, melanocytes, fibroblasts, transformed melanocytes, and abnormal fibroblasts 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 κ E and ρ E only if there is any loss compared to the initial density of epidermis (E 0 ) and dermis (E 1 ), respectively. However, transformed fibroblasts always produce ECM with a rate of ρ s E . 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 non-dimensional units). ECM is degraded by MMP at a rate of µ. The governing equation for ECM is where E 0 represents initial density of epidermis and E 1 is that of dermis. Note that when a fibroblast is in contact with either a melanocyte or a keratinocyte, we set E 1 to be its maximal value 1.0 to model fibroblast's basement membrane generation. The function H is the Heaviside step function defined as MMPs (P) are produced by abnormal fibroblasts, diffuse at a rate D p , 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 β s p . The governing equation is

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 N η . 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 (t c ). 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.
Cell-cell adhesion: In order to incorporate one of the important mechanical aspects of skin into the model, we define an adhesion level (0 ≤ A η ≤ 1) at each node η ∈ Ω explicitly by using cell-cell adhesive coefficients C η,ξ , where ξ is one of its orthogonal neighbors (i.e. ξ ∈ N η ). 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 [3,55]. To model E-cadherin 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 C η,ξ for all cases is summarized in Figure 4 A-B. By averaging the four adhesion coefficients C η,ξ at each grid point, we define an adhesion level A η for a cell at that grid point η.  Migration: In vSkin, every cell has capacity to move to one of its orthogonal neighbors (N η ) regulated by its cell-cell adhesion restriction (i.e., Pr(migration) ∼ (1 − A η )), where A η is a cell-cell 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 N η . We first determine the cell-cell adhesion value (A η ). If there is only one empty neighbor in the neighborhood N η , the cell moves to this point with probability (1 − A η ). 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 4C-D.
In addition to the cell-cell 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 )the melanocyte is not allowed to move north and therefore can only move south. ECM gradient will always dictate the movement direction, provided one exists and their is sufficient space to move.
one of its orthogonal neighbors with regards to its cell-cell 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 cell-matrix interaction mechanism [56,57] and migrate via haptotaxis. Based on rule (iii), a melanocyte is allowed to be moved to another grid point ξ in N η 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 4E-F.
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 G E has to be sufficient (G E > B k ) for its division. Then, a keratinocyte proliferates with a probability of Pr(keratinocyte division) ∼ (1 − T β ). A melanocyte proliferates if (i) the melanin unit is disrupted and (ii) the concentration of bFGF is sufficient for cell division (G b > B m ). A fibroblast is allowed to proliferate if the growth factor bFGF (G b ) is greater than the threshold (B f ). 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 (α E , α b , α T β ) 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 (N η ). 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, S k for a keratinocyte and S m 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 α b , 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 (Pr(fibroblast death) ∼ [G b ]/age). 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.
1. Define a rectangular grid (Ω) that determines both microenvironmental concentrations and cell positions.
3. Place fibroblasts, melanocytes and keratinocytes on the grid.
5. Cell-cell adhesion level A η is determined for each η ∈ Ω.
6. The concentration of microenvironmental factors and adhesion level are coupled with each cell to determine its phenotype.
7. Action associated with the determined phenotype is realized and the grid is updated accordingly.

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 non-dimensionalized [45].
The diffusion rate of EGF (D E ) is estimated using the relation given in [58] with the molecular weight 133.07 kDa [59]. The estimated value is D E = 2.5 × 10 −8 cm 2 /s. The production rate of EGF by a fibroblast has not been measured. Thus, we estimate it to be a 0.5 (dimensional value = 1.2 fg day −1 cell −1 ) from multiple simulations varying this parameter within range of [0,1] to achieve normal growth (net growth ∼ 0) 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; non-dimensional). 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 [0, 1] to achieve skin homeostasis. We estimate the diffusion rate of TGF β using the relation in [58] and the molecular weight 44.3112 kDa [59]. The estimated value is D T β = 5.8 × 10 −8 cm 2 /s. The production rate of TGF β by a fibroblast is chosen to be 0.21333 fg day −1 cell −1 from in vitro study [60]. The rate for a keratinocyte is assumed to be significantly smaller than that of a fibroblast [61]. Since TGF β is known to mainly regulate the growth of keratinocytes [49][50][51], we assumed that keratinocytes bind to TGF β the most. The rates for a melanocyte and a fibroblast were estimated to achieve normal skin homeostasis. The consumption rates of TGF β are estimated to be 0.2, 0.6, 0.1 for a melanocyte, a keratinocyte, and a fibroblast, respectively. Using the same method [58] with molecular weight 17.1531 kDa, the diffusion rate D b is estimated to be 1.2×10 −7 cm 2 /s. The production rate of bFGF by a fibroblast is also estimated from multiple simulations varying these parameters within a range of [0,1] (non-dimensional). The estimated parameter is 0.05 (0.12 fg day −1 cell −1 ). 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 D p = 5.01 × 10 −8 cm 2 /s in [58,59]. The production rate of MMP by fibroblasts is varied within a range 0.024 -0.12 fg day −1 cell −1 . Estimates for the kinetic parameters λ and µ are not available since they are very difficult to obtain experimentally. Thus, we set non-dimensional 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.
The cell cycle time is set it to be a typical value t c = 24 h. The size of the grid is set to be 300 × 62. The grid size δh is set to be δh = 20 µm (a typical cell diameter), and thus the dimensions of the skin slice we simulate are 6 mm × 1.24 mm. The time step for the microenvironmental variables is set to be δt = 4.8 h. We assume that all cell types have the same volume 3.1416 × 10 −9 cm 3 with a radius 10 µm, and therefore the base number of cells is set to be 3.1831 × 10 8 . The background density of ECM is set to be 5.0 × 10 −4 g cm 3 [63]. The background concentration of microenvironmental variables is Dermis ECM density -0.7 assumed to be G 0 = 10 ng cm −3 [64]. We assume that the concentration of growth factors above which a cell proliferates in normal skin is set to the non-dimensional threshold value 0.02 for B m , 0.02 for B f and 0.01 for B k . Since the concentrations below which cells die from starvation are also variables, we set a starvation threshold to be 1% (S k ), 2% (S m ), and 0.1% of the background microenvironmental concentration G 0 for keratinocyte, melanocyte, and abnormal fibroblast respectively.

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).

Co-culture 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 pre-senescent fibroblasts were plated upon 10 cm tissue culture plates overnight before the addition of Adeno-GFP 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 Nikon-TS100 inverted fluorescence microscope.

Western blotting
Cells were plated and grown overnight. To detect ADAM-9 levels, ADAM-9 antibody (ab) was used for immunoblotting (IB). For all IB experiments, proteins were denatured prior to separation on 6-18% Tris-Glycine gels. Proteins were transferred to polyvinylidene fluoride (PVDF) and blocked with 5% non-fat milk in Tris-Buffered Saline Tween-20 (TBST). Blots were incubated overnight in primary ab diluted according to the manufacturers 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 ADAM-9 activity, proteins were not denatured prior to separation upon a casein gel. Cleavage of casein was imaged using an Epson V300 photo scanner.
Three-dimensional spheroid assays WM793 melanoma spheroids were prepared using the liquid overlay method [65]. Spheroids were treated with serum free RPMI, pre-senescent fibroblast conditioned RPMI or senescent fibroblast conditioned RPMI. Pictures of the invading spheroids were taken using a Nikon-TS100 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 Alexa-Fluor 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 A-B, 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 6-8 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 (f (t))." We quantify ratios of each cell compartment (ratio of keratinocyte to melanocyte and that of fibroblasts to melanocytes) at each time step (t) and compare with the ratios of the normal state (t 0 ). We also took changes of epidermal thickness into account to monitor skin fitness. Changes from time step t 0 to time step t (t 0 → t) 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, where r represents the sum of two ratios, the ratio of keratinocyte to melanocyte and that of fibroblast to melanocyte, ∆ represents the change from time step t 0 to t, h stands for a epidermal thickness, and w 1 and w 2 are the weights. In this study we choose the same weight (i.e.,w 1 = w 2 = 0.5).

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 (≥ 50) 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.  Figure 7. Wound healing. A: Positive correlation between healing time and wound size (depth) derived from a suite of wound healing simulations with a range to different wound depths (average of 100 runs). The healing time was defined as the time it takes for the basement membrane to be continuously connected. The mean healing time increases as wound size increases. B: Three snapshots from a wound healing simulation, with a wound that was generated by removing a triangular shape slice with height of 0.96 mm and width of 1.86 mm.

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, TGF β , 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 wound-healing 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. After the perturbed skin reaches its new equilibrium, the skin fitness was measured using (8). The EGF and TGF β impulses have a negligible effect on skin fitness. The bFGF and MMP have a greater impact on skin fitness, reducing it to 60% and 0%, respectively. However, whilst the skin structure is quite close to normal (see upper inset vSkin histology) after bFGF perturbation, after MMP over-expression, the skin has quite a different structure (see lower inset vSkin histology)characterized by a lower basement membrane, and the epidermal cells (keratinocytes and melanocytes) have populated the dermal layer.

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 super-physiological microenvironmental changes. Specifically, each growth factor was maximized in the domain during 3, 9, 15 and 21 days. The skin fitness function f (t) in equation (8) was used to quantify abnormality of vSkin after perturbation. Note that we measure skin fitness f (t) 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 TGF β impulses on vSkin were also negligible even though the growth of keratinocytes was inhibited by high TGF β 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 TGF β 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 TGF β 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, TGF β , 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 long-term 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, TGF β , 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 senescent-like 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  Figure 10. Cell interaction network for abnormal skin. Key interactions among normal cells (keratinocyte, melanocyte, fibroblasts), abnormal cells (senescent fibroblast and transformed melanocyte) and microenvironmental variables represented as network. Note, cell coloration is used in all simulations to represent the equivalent cell. Green lines represent promoting events, while red lines imply inhibiting effects. A transformed melanocyte can produce its own growth factor (bFGF) and MMP, while senescent fibroblasts can produce bFGF, MMP and extracellular matrix proteins. senescent fibroblasts compared to non-senescent 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 (f (t)) using the average number of each cell compartment and the average epidermal thickness from  Figure 11. Degree of senescence and its effects on melanocyte growth. Senescent fibroblasts were placed in the center of the dermal layer, among normal fibroblasts, only 10% were chosen to be senescent. We considered a range of senescence in terms of the production rate of bFGF and that of MMP. Specifically, the nondimensional parameter space for senescent fibroblasts was (α s b , α s p ) ∈ {(x, y)|0.05 ≤ x ≤ 0.8, 0.0 ≤ y ≤ 0.04}, where α s b represents bFGF production rate by senescent fibroblasts and α s p indicates MMP production rate by senescent fibroblasts. Multiple realizations (50) with each parameter set (α s b , α s p ) were performed, and the mean skin fitness was quantified at each time step using (8). The first row shows the skin fitness landscapes after two and six years. A2 and A6 represent skin fitness obtained from simulations with mildly senescent fibroblasts (bFGF production rate is 0.05 and MMP production rate is 0.01). B2 and B6 show the fitness with highly senescent fibroblasts (bFGF production rate is 0.8 and MMP production rate is 0.040). Lower panels show representative vSkin histology for each of these points and highlights how different they are from normal skin (see figure 10 for a key of color coded cell types). 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, A 2 ). When senescent fibroblasts have stronger secretory phenotypes (B 2 and B 6 ) the melanocyte population expands rapidly and skin fitness deteriorates to -0.4 (B 2 ). Interestingly, the fitness subsequently increases to 0.0 as time progresses (see the difference between B 2 and B 6 ). This is due to the rapid growth of melanocytes initially (B 2 ), 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.  Figure 12. Degree of melanocyte malignancy and its effects on melanoma progression. The degree of melanocyte transformation was varied by increasing production rate of bFGF and MMP. The parameter range was given as (β s b , β s P ) ∈ {(x, y)|0.01 ≤ x ≤ 0.05, 0.0 ≤ y ≤ 0.04}, where β s p represents bFGF production rate and β s p describes MMP production rate of a transformed melanocyte. The first row shows the mean skin fitness landscape after 0.25 and 0.5 years (averaged over 50 runs). A 0.25,0.5 represents the fitness from simulations with minimally transformed melanocytes (bFGF and MMP production rates of 0.01), while B 0.25,0.5 indicates the fitness achieved with more aggressive melanoma cells (bFGF production was 0.05 and MMP was 0.04. Lower panels show representative vSkin histology for each of these points (see figure 10 for a key of color coded cell types).

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 (A 0.25,0.5 → B 0.25,0.5 ). The mutant melanocytes disintegrate skin structure and completely take over as shown in a sample snapshot at B 0.5 . This vSkin transformation recapitulates a type of melanoma which progresses very rapidly.  Figure 13. The synergistic effects of senescent fibroblasts and minimally transformed melanocytes in driving melanoma development. The degree of senescence was varied as in the Figure 11. The transformed melanocytes (blue) are assumed to be only minimal transformed, i.e., loss of growth inhibition. Mildly senescent fibroblasts A 2,6 still manage to produce a small melanoma nodule (lower panel) but only impact skin fitness a little (upper panel). However, increasing the fibroblast senescence significantly impacts melanoma progression, destroying the skin (virtual histology, B 2,6 ) and reducing the fitness all the way to zero (upper panel). See figure 10 for a key of color coded cell types.

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 B 2,6 ) can transform essentially normal skin (with the exception of the initial mutant melanocyte) to a pathologic state (see density graph B 2,6 ). Therefore, senescent fibroblasts have the potential to propagate and maintain melanoma development, provided a minimally transformed melanocyte population already exists.  Figure 14. Senescent fibroblasts enhance melanoma cell growth and invasion. A: Senescent fibroblasts aid the growth of early stage (non-tumorigenic) melanoma cells. Equal numbers of GFP-tagged (green) WM35 melanoma cells were seeded on top of either primary skin human fibroblasts (pre-senescent) or fibroblasts treated with radiation and allowed to undergo senescence (senescent) and allowed to grow for 72 hrs. Graph shows the mean (of three independent experiments) relative growth rates (relative to WM5 cells grown alone) where melanoma cells were plated on top of pre-senescent or senescent fibroblasts. B: Poorly invasive melanoma cells (WM793) were grown as spheroids alone, in co-culture with pre-senescent fibroblasts or in co-culture with senescent fibroblasts. The spheroids were then implanted in collagen and allowed to invade. Mean increases in invasion were calculated using ImageJ.

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 co-culture assay.
We explored the prediction that senescent fibroblasts enhanced the growth and invasion of minimally transformed melanocytes in a series of in vitro fibroblast co-culture models. It was noted that plating early stage WM35 melanoma cells (which are non-tumorigenic in mice) on top of senescent fibroblasts enhanced their growth compared to non-senescent fibroblasts (Figure 14 A,B). To study invasion in a tissue-like context, we next generated spheroids in which poorly invasive melanoma cells (WM793) were co-cultured with either non-senescent 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 co-cultured with the senescent fibroblasts. Although not direct comparison with the vSkin results, our co-culture 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 (TGF β , 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 short-term 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 over-production 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 cell-cell 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 hyperpla-   Anti-MMP therapy effects the growth of all populations but is more effective for melanoma progression driven by senescent fibroblasts ((h)-(j)) than melanocyte transformation driven melanomas ((f) and (g)).
sia ( 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 co-operation between mutant melanocytes and senescent fibroblasts ( Figure 13). Given that most known driver mutations in melanoma are associated with entry into oncogene-induced 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 V600E-induced 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 host-tumor 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 MMP-1 and MMP-2 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, anti-MMP therapy might be a good treatment option. Indeed, there is already evidence that mutations in MMP-8, 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 anti-MMP 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 anti-MMP 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 anti-oncogenic 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.