We present an integrated study to understand the key role of senescent fibroblasts in driving melanoma progression. Based on the hybrid cellular automata paradigm, we developed an in silico model of normal skin. The model focuses on key cellular and microenvironmental variables that regulate 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 model predicted the important role of the skin microenvironment in melanoma initiation and progression. Our in vitro experiments showed that dermal fibroblasts, which are an important source of growth factors in the skin, adopt a secretory phenotype that facilitates cancer cell growth and invasion when they become senescent. Our coculture experiments showed that the senescent fibroblasts promoted the growth of nontumorigenic melanoma cells and enhanced the invasion of advanced melanoma cells. Motivated by these experimental results, we incorporated senescent fibroblasts into our 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. The interaction between senescent fibroblasts and the early-stage melanoma cells leads to melanoma initiation and progression. Of microenvironmental factors that senescent fibroblasts produce, proteases are shown to be one of the key contributing factors that promoted melanoma development from our simulations. Although not a direct validation, 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 may create a prooncogenic skin microenvironment that cooperates with mutant melanocytes to drive melanoma initiation and progression and should therefore be considered as a potential future therapeutic target. Interestingly, our simulations to test the effects of a stroma-targeting therapy that negates the influence of proteolytic activity showed that the treatment could be effective in delaying melanoma initiation and progression. Cancer Res; 73(23); 6874–85. ©2013 AACR.

Major Findings

Our integrated approach merging mathematical modeling with experimental and clinical data suggests that the transdifferentiation of fibroblasts due to senescence creates a prooncogenic skin microenvironment that can aid melanoma initiation and progression. Furthermore, our model suggests that the normalization of microenvironment (in terms of proteolytic activity) may be effective in delaying melanoma initiation and progression.

Quick Guide to Model and Major Assumptions

Hybrid cellular automata model

We consider three key cell types: melanocytes, keratinocytes, and fibroblasts. These cells are defined as points on a two-dimensional grid that represents a cross section (6 mm ×1.24 mm, 300 × 62 cells; cell diameter: 20 μm) of human skin (Figs. 1 and 2A–B). Each grid point may be occupied by up to five microenvironmental variables: epidermal growth factor (EGF), basic fibroblasts growth factor (bFGF), TGF-β, matrix-degrading enzyme (MDE), and extracellular matrix (ECM)/basement membrane. The discrete cells are coupled with these microenvironmental concentrations to define the hybrid cellular automata (HCA) model. The cells on the lattice influence the microenvironmental variable concentrations through production and consumption, but are in turn themselves affected by these concentrations, as they serve as regulators of cell behavior.

Figure 1.

Cell–cell interaction network. The solid lines represent interaction among keratinocytes, melanocytes, and fibroblasts in a normal skin condition. The dashed lines represent additional interactions in an abnormal skin conditions.

Figure 1.

Cell–cell interaction network. The solid lines represent interaction among keratinocytes, melanocytes, and fibroblasts in a normal skin condition. The dashed lines represent additional interactions in an abnormal skin conditions.

Close modal
Figure 2.

A, hematoxylin and eosin (H&E) stained cross section of human normal skin. The epidermis contains a basal layer, melanocytes, and keratinocytes. The dermis contains fibroblasts and extracellular matrix. B, model domain with its key cell types. Continuously connected nodes (black) represent the basement membrane. C and E, flowcharts of the life cycle for fibroblasts, keratinocytes, and melanocytes. C, fibroblasts produce bFGF, EGF, and TGF-β, and proliferate if there is enough bFGF. Low concentrations of bFGF trigger fibroblast death. D, keratinocytes produce bFGF and proliferate when there is enough EGF and a low concentration of TGF-β. An insufficient level of EGF triggers keratinocyte death. E, melanocytes start their life cycle by checking the melanin unit and then decide to migrate or not. Sufficient bFGF stimulates melanocyte proliferation, whereas insufficient bFGF triggers cell death.

Figure 2.

A, hematoxylin and eosin (H&E) stained cross section of human normal skin. The epidermis contains a basal layer, melanocytes, and keratinocytes. The dermis contains fibroblasts and extracellular matrix. B, model domain with its key cell types. Continuously connected nodes (black) represent the basement membrane. C and E, flowcharts of the life cycle for fibroblasts, keratinocytes, and melanocytes. C, fibroblasts produce bFGF, EGF, and TGF-β, and proliferate if there is enough bFGF. Low concentrations of bFGF trigger fibroblast death. D, keratinocytes produce bFGF and proliferate when there is enough EGF and a low concentration of TGF-β. An insufficient level of EGF triggers keratinocyte death. E, melanocytes start their life cycle by checking the melanin unit and then decide to migrate or not. Sufficient bFGF stimulates melanocyte proliferation, whereas insufficient bFGF triggers cell death.

Close modal

The interaction network (Fig. 1) defines the key interacting elements of the model and is based on our current biologic understanding. Each fibroblast is assumed to produce three growth factors, EGF, bFGF, and TGF-β. EGF promotes growth of keratinocytes, whereas TGF-β inhibits the growth. The growth of a melanocyte is controlled by both bFGF and physical interaction (contact inhibition) with neighboring keratinocytes. In the abnormal skin model, we include two types of abnormal skin cells, transformed melanocytes and senescent fibroblasts. Each transformed melanocyte produces bFGF and MDE, whereas each senescent fibroblast produces bFGF, MDE, and extracellular matrix proteins. The properties of senescent fibroblasts were defined on the basis of our microarray analysis (Table 1 and Supplementary Fig. S1).

Table 1.

Microarray profile of senescent fibroblasts

 Extracellular matrix proteins Collagen IV, collagen XV, collagen XXI, fibronectin III, and laminin Paracrine growth factors Fibroblast growth factor, IGF binding protein V Angiogenic proteins Platelet-derived growth factor D, VEGF Proteolytic enzyme ADAM metallopeptidase V
 Extracellular matrix proteins Collagen IV, collagen XV, collagen XXI, fibronectin III, and laminin Paracrine growth factors Fibroblast growth factor, IGF binding protein V Angiogenic proteins Platelet-derived growth factor D, VEGF Proteolytic enzyme ADAM metallopeptidase V

Microenvironmental factors

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. Because these continuous variables interact with discrete cells, for clarity, we present only the discretized version of each equation.

The governing equations for the three microenvironmental variables (EGF, bFGF, and TGF-β) have the same form. The factors diffuse (diffusion coefficient D), are produced by fibroblasts (F) at a rate of α, consumed by keratinocytes (K) at a rate of γ and/or melanocytes (M) at a rate of ξ, and decayed naturally at a rate of λ. Thus, we define a general form for these three variables (see Supplementary Section 3.1 for a complete set of equations). The general form for the variables (EGF, bFGF, and TGF-β), denoted by G at a lattice point η ≡ (ηx, ηy), is

where δt is the time step. The operator |$\Delta _2$| defines a two-dimensional discrete Laplacian,

where δh is the grid size.

ECM (E) is produced by both fibroblasts and keratinocytes depending on the local concentration of ECM. In other words, keratinocytes and fibroblasts produce ECM with a rate of |$\kappa _E$| and |$\rho _E$| only if there is a loss compared with the initial epidermal (E0) or dermal (E1) densities, respectively. However, senescent fibroblasts |$F^{s}_{\eta}$| always produce ECM with a rate of |$\rho _E^s$|⁠. Whenever fibroblasts are in contact with keratinocytes or melanocytes, they can make a much denser matrix that creates and maintains the new basement membrane (1, 2). Basement membrane is defined as a continuously connected curve of grid points, each of whose ECM densities is maximal (1.0 in nondimensional units). ECM is degraded by MDE at a rate of μ. The governing equation for ECM is

where E0 represents initial density of epidermis and E1 is that of dermis. When a fibroblast is in contact with either a melanocyte or a keratinocyte, we set E1 to be its maximal value 1.0 to model basement membrane generation by the fibroblast. The function H is the Heaviside step function, used as a switch for ECM production.

MDE (P) is produced by senescent fibroblasts, diffuses at a rate Dp, and is depleted as it degrades the ECM (E). Note that when melanocytes are transformed |$(M^s_{\eta})$|⁠, they can also produce MDE. We assume that they do so with a rate of |$\beta _p^s$|⁠. The governing equation is

Boundary conditions for all microenvironmental variables are no-flux on the top 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 concentration of microenvironmental variables in which it resides. Each cell has three possible phenotypes, proliferative, migratory, and dead. Biologic rules for each cell are summarized as flowcharts in Fig. 2C–E. Each cell is allowed to execute only one phenotype at each cell time step. If a cell does not execute any of these three phenotypes, it remains as a quiescent cell (Supplementary Section 3.2). We assume the lower boundary of the domain allows cell flux for all cellular variables and keep track for the numbers of each cell type that leaves. A Dirichlet boundary condition was imposed on the top of the domain for keratinocytes. No-flux boundary conditions were imposed for both melanocytes and fibroblasts on the top of the domain. Periodic boundary conditions were imposed on the left and right sides of the domain.

Melanoma is the most devastating form of skin cancer (3, 4), arising from the malignant transformation of melanocytes, the pigment-producing cells of the skin. The acquisition of activating mutations in the kinase BRAF is an important initiating event in 50% of melanomas (5). Despite BRAF mutations being important for melanoma development, the majority of benign nevi that harbor BRAF mutations (6) rarely undergo full malignant transformation. Instead, benign nevi cease proliferation and remain quiescent for decades, entering a state of oncogene-induced senescence (OIS; refs. 7, 8). Although there is some suggestion that activation of phosphoinositide 3-kinase (PI3K)/protein kinase B (PKB or AKT) signaling, particularly in the context of loss of the tumor suppressor PTEN may be required for escape from senescence (9), the precise molecular events that underlie the transformation of nevus cells 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 excellent studies have shown that senescent stromal fibroblasts stimulate premalignant and malignant epithelial cells to grow in cell culture and to form tumors in mice (reviews in refs. 10, 11, and 12–15). Mechanistically, this seems to involve the secretion of factors from senescent fibroblasts, named as senescence-associated secretory phenotypes by Coppé and colleagues (12). The secretory factors include matrix metalloproteinase (MMP)-3 and interleukin (IL)-6, which in turn remodel the microenvironment, alter epithelial differentiation, promote endothelial cell motility, and stimulate cancer cell growth both in vitro and in vivo. Following ref. 12, we also undertook our experimental investigation of senescent fibroblasts and observed senescent fibroblasts show altered patterns of ECM expression, growth factor release, and protease activity. Using a suite of in vitro fibroblast coculture experiments, we also learned that these changes help the growth and invasion of melanoma cells.

To gain a better insight into the role of the stromal cells in skin cancer progression in human skin structure, we developed a hybrid multiscale mathematical model that simulates normal skin homeostasis (virtual skin model, vSkin). A number of studies presented excellent computational and mathematical models of skin (16–30). To the best of our knowledge, however, little attention has been paid to normal skin homeostasis and its disruption as a route for melanoma development. This study presents a minimal model of skin, including the dermis, and aims to capture homeostasis first before simulating melanoma formation. We employ a hybrid cellular automata (HCA) approach (31–35). We show that the vSkin model is robust enough to withstand both physical and biochemical perturbations, and still maintains a functional homeostasis. vSkin also suggests the involvement of microenvironmental factors in conjunction with senescent fibroblasts as an important driver of melanoma initiation and progression.

For clarity of composition, we have placed detailed explanations of the model development, parameterization, simulation process, and experimental methods in the Supplementary Material.

### Cell lines and human specimen procurement

Melanoma cells were kindly provided by Dr. Meenhard Herlyn (The Wistar Institute, Philadelphia, PA). GFP adenovirus was kindly provided by Dr. Amer Beg (Moffitt Cancer Center, Tampa, FL).

Representative tumor samples were 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.

### Initialization

We first ran an initial simulation to obtain the starting configuration of the domain (Fig. 3) for all subsequent simulations. We followed the experimental steps in the in vitro three-dimensional (3D) organotypic skin reconstruct experiment (Fig. 3A and 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, it is known to be very stable and homeostatic (36), and contains essentially the same elements as those of vSkin. The initial simulation starts with a dermal layer mixed with fibroblasts that simulate for 2 weeks until 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. Figure 3C and Supplementary Video S1 show how the skin structure develops over the period of a year.

Figure 3.

A, the generation of in vitro organotypic skin reconstructs. A stromal layer of collagen mixed with fibroblasts is placed on top of the acellular collagen, and the layer is incubated to allow the fibroblasts to constrict the collagen. The epithelial cells are layered on top of the fibroblast-containing stromal layer. Twenty-four hours after plating and after the contraction of collagen gel, the whole structure is passed to a steel grid to form air–liquid interface. B, photomicrograph of a mature organotypic culture of a normal skin. C, snapshots of vSkin initialization. D, mean populations dynamics of keratinocytes, melanocytes, and fibroblasts (out of 50 realizations). Sample net growth more than 1 year; E, keratinocytes; F, melanocytes; and G, fibroblasts.

Figure 3.

A, the generation of in vitro organotypic skin reconstructs. A stromal layer of collagen mixed with fibroblasts is placed on top of the acellular collagen, and the layer is incubated to allow the fibroblasts to constrict the collagen. The epithelial cells are layered on top of the fibroblast-containing stromal layer. Twenty-four hours after plating and after the contraction of collagen gel, the whole structure is passed to a steel grid to form air–liquid interface. B, photomicrograph of a mature organotypic culture of a normal skin. C, snapshots of vSkin initialization. D, mean populations dynamics of keratinocytes, melanocytes, and fibroblasts (out of 50 realizations). Sample net growth more than 1 year; E, keratinocytes; F, melanocytes; and G, fibroblasts.

Close modal

### Normal skin: dynamic homeostasis and robustness

As a typical HCA method inevitably contains stochastic components, any quantification of simulations outcomes is the average of multiple realizations (≥50). Our simulations show that vSkin can reach a stable cell net growth rate and total cell number in the domain. Initially, keratinocytes rapidly grow for the first 2 months then as growth factor (EGF) is consumed the keratinocytes become deprived of EGF (near the surface) and switch to a quiescent phenotype or start to die. As keratinocytes near the surface are continuously shed due to either growth factor deprivation or falling off the upper boundary, empty space is created allowing keratinocytes to migrate upward (closer to the surface), leading to free space for the basal keratinocytes (above basement membrane) to proliferate into, if the growth factor conditions are satisfied. The keratinocyte population maintains this dynamic but stable state (Fig. 3D and E). The melanocyte population also rapidly increases initially but soon finds its equilibrium (Fig. 3D and F). Fibroblasts retain their initial population levels (Fig. 3D and G).

To determine the robustness of the vSkin system, two different types of perturbations were applied. First, we tested whether 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 (Fig. 4A snapshot and Supplementary Video S2). As a new space was introduced into vSkin, cells nearby the space have an increased tendency to move toward the new space. In the dermal layer, fibroblasts migrate into the injury site and start to produce growth factors and extracellular matrix proteins. These newly produced growth factors promote growth of keratinocytes in epidermis, thereby starting to fill the void (healing procedure). We changed the size of the wound by increasing injury depth, and quantified the relationship between healing time and injury size. The healing time from wounding tends to increase with wound size (Fig. 4A).

Figure 4.

A, positive correlation between healing time, defined as the time it took the wound site to form a complete basement membrane, and wound size (depth) derived from a suite of wound healing simulations with a range to different wound depths (average of 100 runs). B, skin fitness from maximizing (1.0) the concentration of each growth factor over the entire domain during 3, 9, 15, and 21 days (averaged over 50 runs).

Figure 4.

A, positive correlation between healing time, defined as the time it took the wound site to form a complete basement membrane, and wound size (depth) derived from a suite of wound healing simulations with a range to different wound depths (average of 100 runs). B, skin fitness from maximizing (1.0) the concentration of each growth factor over the entire domain during 3, 9, 15, and 21 days (averaged over 50 runs).

Close modal

The second perturbation involves manipulation of microenvironmental factors. Specifically, each growth factor was maximized in the domain for different periods of time (i.e., 3, 9, 15, and 21 days). The degree of abnormality of vSkin after perturbation was characterized by a defined metric, “skin fitness” [see equation (5s) in the Supplementary Material], which measured how close each of cellular compartments was to normal skin after perturbation. The skin fitness is scored from −1 to +1, where +1 represents the maximal fitness and a normal skin state.

Both EGF and TGF-β impulses did not change skin fitness (Fig. 4B) or structure (Supplementary Video S3). Basic fibroblast growth factor (bFGF) impulses decreased skin fitness by 40% due to increased melanocyte population (Fig. 4B). vSkin was able to recover from small changes induced by short-term matrix-degrading enzyme (MDE) overexpression (3 days), but failed to compensate for long-term (≥9 days) MDE overexpression (Fig. 4B). Taken together, the results show that vSkin is remarkably robust to microenvironmental perturbations that are present in normal skin. However, when vSkin was insulted by a factor not typically present (e.g., MDE), its robustness was dependent on the duration of exposure. With long-term exposure of MDE, vSkin transformed to a pathologic state, implying that regulation of key microenvironmental factors can contribute to vSkin transformation to an abnormal state.

### Senescence in fibroblasts

Fibroblasts are the primary source for many of the microenvironmental factors that regulate skin homeostasis. Therefore, disruption of fibroblasts should also have a significant impact on skin homeostasis. To this end, we generated senescent fibroblasts by following ref. 12. The onset of senescence in the fibroblasts was associated with a marked change in ECM constituents, growth factors, and proteases (Table 1 and Supplementary Fig. S1).

Next, we examined the effects of senescent fibroblasts in the growth and invasion of 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 with presenescent fibroblasts (Fig. 5A). To study invasion in a tissue-like context, we next generated spheroids in which poorly invasive melanoma cells (WM793) were cocultured with either presenescent or senescent human skin fibroblasts in a 3D collagen gel (Fig. 5B). 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. Taken together, these in vitro experiments suggest that senescent fibroblasts promote the growth of early-stage melanoma cells (nontumorigenic) and the invasion of advanced melanoma cells.

Figure 5.

A, equal numbers of GFP-tagged (white) WM35 melanoma cells were seeded on top of either primary skin human fibroblasts (presenescent) or fibroblasts treated with radiation and allowed to undergo senescence (senescent) and allowed to grow for 72 hours. Graph shows the mean (of three independent experiments) relative growth rates (relative to WM5 cells grown alone) in which melanoma cells were plated on top of presenescent or senescent fibroblasts. B, poorly invasive melanoma cells (WM793) were grown as spheroids alone, in coculture with presenescent fibroblasts, or in coculture with senescent fibroblasts. The spheroids were then implanted in collagen and allowed to invade. Mean increases in invasion were calculated using ImageJ.

Figure 5.

A, equal numbers of GFP-tagged (white) WM35 melanoma cells were seeded on top of either primary skin human fibroblasts (presenescent) or fibroblasts treated with radiation and allowed to undergo senescence (senescent) and allowed to grow for 72 hours. Graph shows the mean (of three independent experiments) relative growth rates (relative to WM5 cells grown alone) in which melanoma cells were plated on top of presenescent or senescent fibroblasts. B, poorly invasive melanoma cells (WM793) were grown as spheroids alone, in coculture with presenescent fibroblasts, or in coculture with senescent fibroblasts. The spheroids were then implanted in collagen and allowed to invade. Mean increases in invasion were calculated using ImageJ.

Close modal

### Benign nevi

Motivated by our experiments, we incorporated senescent fibroblasts into the vSkin model to explore the impact of the senescent fibroblasts on normal melanocyte growth and skin structure. The degree of senescence was varied to test whether the degree differently transforms normal skin structure. The degree of senescence specifically refers to how much bFGF and MDE the fibroblasts produce, so the most senescent phenotype produces maximal bFGF and MDE. Specifically, the nondimensional parameter space for senescent fibroblasts was |$(\alpha _b^s,\alpha _p^s) \in \{ (x,y\left| {0.05 \le x \le 0.8,0.0 \le y \le 0.04} \right.)\}$|⁠, where |$\alpha _b^s$| represents bFGF production rate by senescent fibroblasts and |$\alpha _p^s$| indicates MDE production rate by senescent fibroblasts. From the skin fitness (average over 50 realizations for each degree) analysis, we learned that any degree of senescence for the fibroblasts significantly enhanced the growth of melanocytes and promoted invasion (left; Fig. 6A). Furthermore, it is clear that the intensity of skin disruption is dependent on the degree of fibroblast senescence. When senescent fibroblasts have weak secretory phenotypes, 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 (Fig. 6A, snapshot a2; Supplementary Video S4). When senescent fibroblasts have stronger secretory phenotypes (Fig. 6A, b2 and b6; Supplementary Video S5) the melanocyte population expands rapidly and skin fitness deteriorates to −0.4 (Fig. 6A, b2). Interestingly, the fitness subsequently increases to 0.0 as time progresses (see the difference between Fig. 6A, b2 and b6). This is due to the rapid growth of melanocytes initially, driven by the senescent fibroblasts, but this excessive growth leads to over consumption of the growth factors for both fibroblasts and melanocytes, resulting in cell death, followed by a return to a more normal structure and improved skin fitness. These simulations suggest that senescent fibroblasts may play a critical role in creating an environment ripe for growth and invasion of normal melanocytes and cause the creation nevi-like structures.

Figure 6.

Key of color-coded cell types: yellow (keratinocytes), brown (melanocytes), pink (fibroblasts), red (senescent fibroblasts), and blue (transformed melanocytes). A, mildly senescent fibroblasts (bFGF rate: 0.05 and MDE rate: 0.01) did not change skin fitness and skin structure after 2 years (a2), but produced a very small mole after 6 years (a6). Highly senescent fibroblasts (bFGF rate: 0.8 and MDE rate: 0.04) significantly reduced skin fitness (b2 = −0.3) and promoted normal melanocyte invasion (b2) after 2 years. After 6 years, the skin fitness is increased to 0.0 and skin structure becomes closer to normal skin (b6). B, mildly senescent fibroblasts reduced skin fitness 50% (c2,6) and produced a small melanoma nodule after 6 years (c6). Increasing the senescence significantly reduces the fitness all the way to its minimum (−1.0) and destroys the skin (d6). C, the total population of transformed melanocytes was quantified as a measure of tumor burden. Representative vSkin snapshot is also shown (left). The population of minimally mutated melanocytes (barely visible green line) remains constant. Senescent fibroblasts with weaker secretory phenotypes lead to a slower melanocyte proliferation and invasion [(i) and (ii)]. The most senescent fibroblasts [(iii)] help minimally mutated melanocytes to proliferate rapidly and progress to full melanoma. The mutated melanocytes that produce bFGF alone form melanoma in situ in which the melanoma cells are constrained to the epidermis [C, (iv)]. The most mutated melanocytes produce both bFGF and MMP and progress the most rapidly [C, (v)].

Figure 6.

Key of color-coded cell types: yellow (keratinocytes), brown (melanocytes), pink (fibroblasts), red (senescent fibroblasts), and blue (transformed melanocytes). A, mildly senescent fibroblasts (bFGF rate: 0.05 and MDE rate: 0.01) did not change skin fitness and skin structure after 2 years (a2), but produced a very small mole after 6 years (a6). Highly senescent fibroblasts (bFGF rate: 0.8 and MDE rate: 0.04) significantly reduced skin fitness (b2 = −0.3) and promoted normal melanocyte invasion (b2) after 2 years. After 6 years, the skin fitness is increased to 0.0 and skin structure becomes closer to normal skin (b6). B, mildly senescent fibroblasts reduced skin fitness 50% (c2,6) and produced a small melanoma nodule after 6 years (c6). Increasing the senescence significantly reduces the fitness all the way to its minimum (−1.0) and destroys the skin (d6). C, the total population of transformed melanocytes was quantified as a measure of tumor burden. Representative vSkin snapshot is also shown (left). The population of minimally mutated melanocytes (barely visible green line) remains constant. Senescent fibroblasts with weaker secretory phenotypes lead to a slower melanocyte proliferation and invasion [(i) and (ii)]. The most senescent fibroblasts [(iii)] help minimally mutated melanocytes to proliferate rapidly and progress to full melanoma. The mutated melanocytes that produce bFGF alone form melanoma in situ in which the melanoma cells are constrained to the epidermis [C, (iv)]. The most mutated melanocytes produce both bFGF and MMP and progress the most rapidly [C, (v)].

Close modal

### Melanoma

We integrated both transformed melanocytes and senescent fibroblasts together in the vSkin model, to investigate how interactions between these abnormal populations contribute to the melanoma development. To highlight the impact of these interactions, we only considered minimally transformed melanocytes in the following simulations. Note that these minimally transformed melanocytes have only lost control over proliferation and are incapable of producing cancer in the vSkin model independently. Figure 6B highlights how senescent fibroblasts help minimally transformed melanocytes proliferate more rapidly and invade into the dermis. Weakly senescent fibroblasts produced a nodule of minimally transformed melanocytes (Fig. 6B, c2,6; Supplementary Video S6) 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 and decreasing skin fitness to below zero (Fig. 6B, d2,6; Supplementary Video S7). From these results, we conclude that senescent fibroblasts can transform essentially normal skin (with the exception of the initial mutant melanocyte) to a pathologic state, a melanoma-like state. Furthermore, the melanoma progression rate is directly dependent on the degree of senescence. Weakly senescent fibroblasts facilitated the slowest melanoma growth [Fig. 6C, snapshot (i) and orange line]. As senescent fibroblasts adopt stronger phenotypes, they produce faster progressing melanomas [Fig. 6C, (ii) and (iii)].

In addition, the vSkin model recapitulates melanoma development driven only by melanocyte transformation. We considered three types of melanoma cells (i.e., most aggressive, mildly aggressive, and minimally aggressive melanoma cells). Minimally aggressive types are the minimally transformed cells we used in combination with the senescent fibroblasts. Mildly aggressive melanoma cells are assumed to produce bFGF, whereas the most aggressive ones are assumed to produce both bFGF and MDE. The minimally aggressive cells produce no cancerous growth (Fig. 6C, green line), and the mildly aggressive ones produce melanoma in situ [Fig. 6C, snapshot (iv) and light blue line], whereas the most aggressive ones grow more rapidly and produce invasive melanoma [Fig. 6C, snapshot (v) and dark blue line].

### Clinical relevance

From vSkin simulations, we observed that MDE production by senescent fibroblasts promoted melanoma invasion (independent of significant melanocyte transformation). Figure 7A shows a representative vSkin snapshot as well as the corresponding MDE expression levels. Note that we present MDE expression levels as if we stained a vSkin snapshot for MDE and melanoma cells. We observe that the MDE level is high (yellow) at the interface between tumor and senescent fibroblasts, and low at the tumor core (MDE expression: black).

Figure 7.

A, increased protease (MDE) expression in the interface (green dashed line) between tumor and senescent fibroblasts (left; MDE expression level, yellow) compared with tumor core (right; MDE expression level, black). Note that the dark blue dots in the vSkin MDE expression panel represent nucleus of all cell types (transformed melanocytes, fibroblasts, and senescent fibroblasts) B, representative hematoxylin and eosin-stained section of a nodule of metastatic melanoma, stained with anti-ADAM9. The yellow line demarcates the tumor–stroma interface. There is high ADAM9 expression at the tumor/host interface in the melanoma specimen but weaker, more diffuse expression in the tumor core. C, (i)–(iii), anti-MDE treatment of melanomas driven by senescent fibroblasts [(i) weakly, (ii) mildly, and (iii) strongly senescent fibroblasts]. Both the sizes of the nodule in the dermis and population growth rate are reduced (compare with Fig. 6C). C, (iv) and (v), anti-MDE treatments of melanoma driven by mildly [(iv) and more (v)] aggressive melanoma cells. A small reduction in both invasion and growth is observed.

Figure 7.

A, increased protease (MDE) expression in the interface (green dashed line) between tumor and senescent fibroblasts (left; MDE expression level, yellow) compared with tumor core (right; MDE expression level, black). Note that the dark blue dots in the vSkin MDE expression panel represent nucleus of all cell types (transformed melanocytes, fibroblasts, and senescent fibroblasts) B, representative hematoxylin and eosin-stained section of a nodule of metastatic melanoma, stained with anti-ADAM9. The yellow line demarcates the tumor–stroma interface. There is high ADAM9 expression at the tumor/host interface in the melanoma specimen but weaker, more diffuse expression in the tumor core. C, (i)–(iii), anti-MDE treatment of melanomas driven by senescent fibroblasts [(i) weakly, (ii) mildly, and (iii) strongly senescent fibroblasts]. Both the sizes of the nodule in the dermis and population growth rate are reduced (compare with Fig. 6C). C, (iv) and (v), anti-MDE treatments of melanoma driven by mildly [(iv) and more (v)] aggressive melanoma cells. A small reduction in both invasion and growth is observed.

Close modal

To find the potential clinical relevance of our model predictions, we explored protease expression in human melanoma specimens in similar regions. Because ADAM family proteases, cell surface proteins implicated in both proteolysis and cell adhesion, were found to be upregulated in senescent fibroblasts in both our microarray and in vitro experiments (Supplementary Fig. S1C–S1E), we stained a small cohort of human melanoma specimens for ADAM9. The specimens stained positively for the protease ADAM9 at the leading edge in which the melanoma cells and fibroblasts interact (Fig. 7B). Strong ADAM9 staining was also noted in the outer layers of the tumor in which the influence of the stromal environment was the strongest (Supplementary Fig. S2). In contrast, little ADAM9 staining was noted in the tumor core in which the stromal influence was lacking (Fig. 7B and 2S). It is worth noting that ADAM 9 is a protease thought to be critical for melanoma–fibroblast interactions (37).

### Targeting senescent fibroblasts

So far, we have shown that proteases are a key contributing factor in driving melanoma progression. Thus, anti-MDE therapy might be a good treatment option. Indeed, there is already evidence that mutations in MMP-8, lead 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 approximately 23% of all melanomas (38). To this end, we ran a suite of simulations with anti-MDE treatments, almost as if we applied anti-MDE cream to the vSkin model such that MDE expression was blocked completely. Specifically, the MDE concentration was completely knocked down to zero uniformly in the whole domain, when the size of initial melanocyte population had doubled, until the simulation was finished. The simulations show that the treatment was effective in decreasing both growth and invasion of melanoma cells (Fig. 6C vs. Fig. 7C). Interestingly, the treatment was far more effective in melanomas driven by senescent fibroblasts and minimally transformed melanocytes [(i)–(iii) vs. (iv)–(v)].

The regulation of skin homeostasis involves a complex dialogue between keratinocytes, melanocytes, and fibroblasts, which 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, the very factors that normal skin seems to be most sensitive to.

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 melanocytes, resulting in the development of melanocytic hyperplasia (Fig. 6A). When the melanocytes were transformed even minimally, their cooperation with senescent fibroblasts could lead to melanoma initiation and progression [Fig. 6B and Fig. 6C (i)–(iii)]. Given that most known driver mutations in melanoma are associated with entry into 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 (9) 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 (39). The most important implication is that we can not only consider melanoma as a transformation of the melanocyte population but must also consider the fibroblast population. It is tempting to speculate that melanoma may be a disease of the aged (40–42), when progressive age-related fibroblast senescence could cooperate with accumulating UV-induced mutations in melanocytes to drive melanoma development.

J.L. Messina is a consultant/advisory board member of GlaxoSmithKline. No potential conflicts of interest were disclosed by the other authors.

Conception and design: E. Kim, K.S.M. Smalley, A.R.A. Anderson

Development of methodology: E. Kim, S.S. Maria-Engler, D. Basanta, A.R.A. Anderson

Acquisition of data (provided animals acquired, and managed patients, provided facilities, etc.): E. Kim, V. Rebecca, I.V. Fedorenko, J.L. Messina, R Mathew, K.S.M. Smalley

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): E. Kim, R Mathew

Writing, review, and/or revision of the manuscript: E. Kim, K.S.M. Smalley, A.R.A. Anderson

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): E. Kim

Study supervision: A.R.A. Anderson

The authors thank Dr. David Basanta for his initial input on the cell interaction network and crucially for providing his code that was used as the foundation for the vSkin code.

This work was funded by Moffitt Cancer Center PS-OC NIH/NCI 1U54CA143970-01 and R01 CA 161107-01.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Andriani
F
,
Margulis
A
,
Lin
N
,
Griffey
S
,
Garlick
JA
.
Analysis of microenvironmental factors contributing to basement membrane assembly and normalized epidermal phenotype
.
J Invest Dermatol
2003
;
120
:
923
31
.
2.
Lee
DY
,
Cho
KH
.
The effects of epidermal keratinocytes and dermal fibroblasts on the formation of cutaneous basement membrane in three-dimensional culture systems
.
Arch Dermatol Res
2005
;
296
:
296
302
.
3.
Chin
L
.
The genetics of malignant melanoma: lessons from mouse and man
.
Nat Rev Cancer
2003
;
3
:
559
70
.
4.
Miller
AJ
,
Mihm
JMC
.
Melanoma
.
N Engl J Med
2006
;
355
:
51
65
.
5.
Davies
H
,
Bignell
GR
,
Cox
C
,
Stephens
P
,
Edkins
S
,
Clegg
S
, et al
Mutations of the BRAF gene in human cancer
.
Nature
2002
;
417
:
949
54
.
6.
Pollock
PM
,
Harper
UL
,
Hansen
KS
,
Yudt
LM
,
Stark
M
,
Robbins
CM
, et al
High frequency of BRAF mutations in nevi
.
Nat Genet
2003
;
33
:
19
20
.
7.
Dhomen
N
,
Reis-Filho
JS
,
da Rocha Dias
S
,
Hayward
R
,
Savage
K
,
Delmas
V
, et al
Oncogenic BRAF induces melanocyte senescence and melanoma in mice
.
Cancer Cell
2009
;
15
:
294
303
.
8.
Michaloglou
C
,
Vredeveld
LCW
,
Soengas
MS
,
Denoyelle
C
,
Kuilman
T
,
van der Horst
CMAM
, et al
BRAF V600E-associated senescence-like cell–cycle arrest of human naevi
.
Nature
2005
;
436
:
720
4
.
9.
Vredeveld
LCW
,
Possik
PA
,
Smit
MA
,
Meissl
K
,
Michaloglou
C
,
Horlings
HM
, et al
Abrogation of BRAF V600E-induced senescence by PI3K pathway activation contributes to melanomagenesis
.
Genes Dev
2012
;
26
:
1055
69
.
10.
Bhowmick
N
,
Neilson
E
,
Moses
HL
.
Stromal fibroblasts in cancer initiation and progression
.
Nature
2004
;
18
:
332
7
.
11.
Krtolica
A
,
Campisi
J
.
Cancer and aging: a model for the cancer promoting effects of the aging stroma
.
Int J Biochem Cell Biol
2002
;
34
:
1401
14
.
12.
Coppé
J-P
,
Patil
CK
,
Rodier
F
,
Sun
Y
,
Muñoz
DP
,
Goldstein
J
, et al
Senescence-associated secretory phenotypes reveal cell-nonautonomous functions of oncogenic RAS and the p53 tumor suppressor
.
PLoS Biol
2008
;
6
:
2853
68
.
13.
Krtolica
A
,
Parrinello
S
,
Lockett
S
,
Desprez
PY
,
Campisi
J
.
Senescent fibroblasts promote epithelial cell growth and tumorigenesis: a link between cancer and aging
.
Proc Natl Acad Sci U S A
2001
;
98
:
12072
7
.
14.
Parrinello
S
,
Coppe
J-P
,
Krtolica
A
,
Campisi
J
.
Stromal-epithelial interactions in aging and cancer: senescent fibroblasts alter epithelial cell differentiation
.
J Cell Sci
2005
;
118
:
485
96
.
15.
Laberge
R-M
,
P
,
Campisi
J
,
Desprez
P-Y
.
Epithelial-mesenchymal transition induced by senescent fibroblasts
.
Cancer Microenviron
2012
;
5
:
39
44
.
16.
S
,
Sun
T
,
MacNeil
S
,
Holcombe
M
,
Smallwood
R
.
Development of a three dimensional multiscale computational model of the human epidermis
.
PLoS ONE
2010
;
5
:
e8511
e
.
17.
Aylaj
B
,
Luciani
F
,
Delmas
V
,
Larue
L
,
De Vuyst
F
.
Melanoblast proliferation dynamics during mouse embryonic development. Modeling and validation
.
J Theor Biol
2011
;
276
:
86
98
.
18.
Basan
M
,
Joanny
J-F
,
Prost
J
,
Risler
T
.
Undulation instability of epithelial tissues
.
Phys Rev Lett
2011
;
106
:
158101
.
19.
Chatelain
C
,
Balois
T
,
Ciarletta
P
,
Ben Amar
M
.
Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture
.
New J Phys
2011
;
13
:
115013
.
20.
Ciarletta
P
,
Foret
L
,
Ben Amar
M
.
The radial growth phase of malignant melanoma: multi-phase modelling, numerical simulations and linear stability analysis
.
J R Soc Interface
2011
;
8
:
345
68
.
21.
Eikenberry
S
,
Thalhauser
C
,
Kuang
Y
.
Tumor-immune interaction, surgical treatment, and cancer recurrence in a mathematical model of melanoma
.
Plos Comput Biol
2009
;
5
:
e1000362
e
.
22.
Grabe
N
,
Neuber
K
.
A multicellular systems biology model predicts epidermal morphology, kinetics and Ca2+ flow
.
Bioinformatics
2005
;
21
:
3541
7
.
23.
Grabe
N
,
Neuber
K
.
Simulating psoriasis by altering transit amplifying cells
.
Bioinformatics
2007
;
23
:
1309
12
.
24.
Murphy
KE
,
Hall
CL
,
Maini
PK
,
McCue
SW
,
McElwain
DLS
.
A fibrocontractive mechanochemical model of dermal wound closure incorporating realistic growth factor kinetics
.
Bull Math Biol
2012
;
74
:
1143
70
.
25.
Murphy
KE
,
Hall
CL
,
McCue
SW
,
Sean McElwain
DL
.
A two-compartment mechanochemical model of the roles of transforming growth factor β and tissue tension in dermal wound healing
.
J Theor Biol
2011
;
272
:
145
59
.
26.
Schaller
G
,
Meyer-Hermann
M
.
A modelling approach towards epidermal homoeostasis control
.
J Theor Biol
2007
;
247
:
554
73
.
27.
Suetterlin
T
,
Huber
S
,
Dickhaus
H
,
Grabe
N
.
Modeling multi-cellular behavior in epidermal tissue homeostasis via finite state machines in multi-agent systems
.
Bioinformatics
2009
;
25
:
2057
63
.
28.
Sun
T
,
McMinn
P
,
Holcombe
M
,
Smallwood
R
.
Agent based modelling helps in understanding the rules by which fibroblasts support keratinocyte colony formation
.
PLoS ONE
2008
;
3
:
e2129
e
.
29.
Thingnes
J
,
Lavelle
TJ
,
Hovig
E
,
Omholt
SW
.
Understanding the melanocyte distribution in human epidermis: an agent-based computational model approach
.
PLoS ONE
2012
;
7
:
e40377
e
.
30.
Valeyev
NV
,
Hundhausen
C
,
Umezawa
Y
,
Kotov
NV
,
Williams
G
,
Clop
A
, et al
A systems model for immune cell interactions unravels the mechanism of inflammation in human skin
.
PLoS Comput Biol
2010
;
6
:
e1001024
e
.
31.
Anderson
AR
,
Chaplain
MA
.
Continuous and discrete mathematical models of tumor-induced angiogenesis
.
Bull Math Biol
1998
;
60
:
857
99
.
32.
Anderson
AR
,
Pitairn
A
.
Application of the Hybrid Discrete-Continuum Technique, In Polymer and cell dynamics - multiscale modeling and numerical simulations
.
Alt
W
,
Chaplain
M
,
Griebel
M
,
Lenz
J
, editors.
Birkhauser
;
2003
. p.
261
79
.
33.
Anderson
ARA
.
A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion
.
Math Med Biol
2005
;
22
:
163
86
.
34.
Anderson
ARA
,
Sleeman
BD
,
Young
IM
,
Griffiths
BS
.
Nematode movement along a chemical gradient in a structurally heterogeneous environment. 2. Theory
.
Fundam Appl Nematol
1997
;
20
:
165
72
.
35.
Basanta
D
,
Strand
DW
,
Lukner
RB
,
Franco
OE
,
Cliffel
DE
,
Ayala
GE
, et al
The role of transforming growth factor-beta-mediated tumor-stroma interactions in prostate cancer progression: an integrative approach
.
Cancer Res
2009
;
69
:
7111
20
.
36.
Brohem
CA
,
Cardeal
LB
,
Tiago
M
,
Soengas
MS
,
Barros
SB
,
Maria-Engler
SS
.
Artificial skin in perspective: concepts and applications
.
Pigment Cell Melanoma Res
2011
;
24
:
35
50
.
37.
Zigrino
P
,
Nischt
R
,
Mauch
C
.
The disintegrin-like and cysteine-rich domains of ADAM-9 mediate interactions between melanoma cells and fibroblasts
.
J Biol Chem
2011
;
286
:
6801
7
.
38.
Palavalli
LH
,
Prickett
TD
,
Wunderlich
JR
,
Wei
X
,
Burrell
AS
,
Porter-Gill
P
, et al
Analysis of the matrix metalloproteinase family reveals that MMP8 is often mutated in melanoma
.
Nat Genet
2009
;
41
:
518
20
.
39.
Dey
JH
,
Bianchi
F
,
Voshol
J
,
Bonenfant
D
,
Oakeley
EJ
,
Hynes
NE
.
Targeting fibroblast growth factor receptors blocks PI3K/AKT signaling, induces apoptosis, and impairs mammary tumor outgrowth and metastasis
.
Cancer Res
2010
;
70
:
4151
62
.
40.
Chao
C
,
Martin
RC
II
,
Ross
MI
,
Reintgen
DS
,
Edwards
MJ
,
Noyes
RD
, et al
Correlation between prognostic factors and increasing age in melanoma
.
Ann Surg Oncol
2004
;
11
:
259
64
.
41.
N
,
Noone
AM
,
Krapcho
M
,
Neyman
N
,
Aminou
R
,
Altekruse
SF
, et al
SEER Cancer Statistics Review, 1975–2009 (Vintage 2009 Populations) based on November 2011 SEER data submission, posted to the SEER web site, April 2012
.
Bethesda, MD
:
National Cancer Institute
;
2011
.
42.
Lachiewicz
AM
,
Berwick
M
,
Wiggins
CL
,
Thomas
NE
.
Epidemiologic support for melanoma heterogeneity using the surveillance, epidemiology, and end results program
.
J Invest Dermatol
2008
;
128
:
1340
2
.