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.
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.
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).
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 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 and only if there is a loss compared with the initial epidermal (E0) or dermal (E1) densities, respectively. However, senescent fibroblasts always produce ECM with a rate of . 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 , they can also produce MDE. We assume that they do so with a rate of . 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.
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.
Materials and Methods
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.
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.
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).
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.
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 , where represents bFGF production rate by senescent fibroblasts and 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.
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].
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).
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.
Disclosure of Potential Conflicts of Interest
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
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.
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.
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).
- Received June 14, 2013.
- Revision received August 29, 2013.
- Accepted September 19, 2013.
- ©2013 American Association for Cancer Research.