Three-dimensional experiments and individual based simulations show that cell proliferation drives melanoma nest formation in human skin tissue

Background Melanoma can be diagnosed by identifying nests of cells on the skin surface. Understanding the processes that drive nest formation is important as these processes could be potential targets for new cancer drugs. Cell proliferation and cell migration are two potential mechanisms that could conceivably drive melanoma nest formation. However, it is unclear which one of these two putative mechanisms plays a dominant role in driving nest formation. Results We use a suite of three-dimensional (3D) experiments in human skin tissue and a parallel series of 3D individual-based simulations to explore whether cell migration or cell proliferation plays a dominant role in nest formation. In the experiments we measure nest formation in populations of irradiated (non-proliferative) and non-irradiated (proliferative) melanoma cells, cultured together with primary keratinocyte and fibroblast cells on a 3D experimental human skin model. Results show that nest size depends on initial cell number and is driven primarily by cell proliferation rather than cell migration. Conclusions Nest size depends on cell number, and is driven primarily by cell proliferation rather than cell migration. All experimental results are consistent with simulation data from a 3D individual based model (IBM) of cell migration and cell proliferation. Electronic supplementary material The online version of this article (10.1186/s12918-018-0559-9) contains supplementary material, which is available to authorized users.


Background
Clusters of melanoma cells, called nests, are a prominent feature of melanoma progression [1,2]. Identifying the presence and characteristics of melanoma nests in human skin is an important diagnostic tool [3,4]. Furthermore, nest size is an important characteristic because larger melanoma nests are associated with more aggressive melanoma [3]. Recent 3D experimental work by Wessels et al. [5] suggests that melanoma nest formation in Matrigel is driven by cell migration. However, nest formation might be different in human skin, where melanoma cells are in contact with other cell types [1,6]. We hypothesise that two different mechanisms could lead to nest formation: (i) cell proliferation, where clusters of melanoma cells are formed primarily through mitosis (Fig. 1a); and (ii) cell migration, where clusters of adhesive melanoma cells form primarily through melanoma cell migration (Fig. 1b). Cell migration occurs over a short time scale of hours, whereas cell proliferation takes place over a much longer time scale of days. Since our work is focused on the role of proliferation, we perform experiments over a period of four days so that we are able to observe and quantify the role of cell proliferation. This choice of experimental time scale means that our experimental observations do not resolve the details of cell migration, which would require a much finer time resolution in the experiments.
We use a 3D human skin experimental model [7,8] to discriminate between these two conceptual models by performing a suite of experiments in which we systematically vary the initial density of proliferative melanoma cells placed on 3D human skin. This initial series of experiments allow us to examine the role of initial cell number in driving nest formation. All experiments are repeated using non-proliferative, gamma-irradiated melanoma cells. We find that higher initial numbers of melanoma cells lead to larger nests, and that cell proliferation leads to dramatically-larger nests. All experimental outcomes are consistent with a series of 3D simulations from an IBM [9]. These results provide insight into the mechanisms driving nest formation, showing that the mechanisms in 3D human skin are different to monoculture experiments performed in Matrigel.

Results and discussion
Confirmation that irradiated melanoma cells do not proliferate and are capable of migrating in a twodimensional barrier assay Experiments involving populations of proliferative melanoma cells are performed using non-irradiated SK-MEL-28 cells [10]. Experiments where melanoma cell proliferation is suppressed are performed using irradiated, but otherwise identical SK-MEL-28 cells [11,12]. The melanoma cells are gamma-irradiated to inhibit mitosis. It is possible that irradiation may have other impacts on cellular behaviour and could also influence DNA functioning [12,13]. We perform a series of live assays to show that irradiation does not affect the adherence or morphology of melanoma cells. These live cell assays involve placing populations of irradiated and nonirradiated melanoma cells on a two-dimensional tissue culture plate and making observations of cell numbers of a period of 24 h [14]. Therefore, these assays provide quantitative information about whether the populations of melanoma cells are capable of proliferating. Results confirm that irradiated melanoma cells do not proliferate. Furthermore, these assays show that irradiation does not cause the cells to die and does not affect cell morphology [see Additional file 1].
Two-dimensional (2D) barrier assays confirm that irradiated melanoma cells survive and migrate. Populations of irradiated melanoma cells are monitored over four days to confirm that irradiation does not impede the ability of cells to migrate. We use circular barrier assays to compare the spatial expansion of irradiated and nonirradiated melanoma cell populations. The leading edge of the spreading populations is detected using ImageJ [15], which also provides an estimate of the area The scale bar in each image is 3 mm. c Data shows the average diameter of the spreading populations as a function of time (n = 3). All data generated using non-irradiated melanoma cells is in blue, and data generated using irradiated melanoma cells is in red. Plots in (c) also show the variability. The error bars correspond to the sample standard deviation (n = 3) occupied by the spreading population of cells. Since the spreading populations of cells maintain an approximately circular shape, we convert the estimates of area into an equivalent diameter and we report data in terms of the diameter of the spreading population. Results are obtained in triplicate. Images in Fig. 2a-b show the increase in the diameter of the spreading cell populations for both irradiated and non-irradiated melanoma cells over four days. The upper row of images in Fig. 2ab, show increased spatial expansion of the population of non-irradiated cells compared to the population of irradiated melanoma cells in the lower row. Since irradiated melanoma cells do not proliferate, we expect that the size of the expanding population of irradiated cells will be smaller than the size of the expanding population of non-irradiated cells [16]. However, the area occupied by the population of irradiated melanoma cells increases over the four-day period, and this increase is due to cell migration alone. To confirm these visual observations, we also quantify the spatial spreading of irradiated and non-irradiated melanoma cell populations. Data in Fig. 2c shows the increase in diameter of both irradiated and non-irradiated melanoma cell populations over four days. At all times considered, the average diameter of the irradiated cell population is less than the average diameter of the non-irradiated cell population. This is expected because the irradiated melanoma cells do not proliferate, and it is known that proliferative populations of cells expand and invade the surrounding empty space faster than non-proliferative populations of cells [9,16]. Most importantly, the experiments initialised with irradiated melanoma cells show an increase in the diameter of the spreading population, confirming that irradiation does not prevent migration. All experiments are performed in triplicate and the averaged results are presented. We now use both, irradiated and non-irradiated melanoma cells in 3D experiments to identify the mechanism that drives melanoma nest formation.
Identifying the dominant mechanism driving melanoma nest formation Nests of melanoma cells are well-characterised histological features of melanoma progression. Early identification of these nests is critical for successful melanoma treatment. However, in addition to examining the presence of melanoma nests, it is important to identify the biological mechanisms that lead to nest formation as this information might be relevant to the development of new drugs. To examine these pathways we use a 3D experimental skin model. Irradiated and non-irradiated melanoma cells are cultured with primary keratinocytes and primary fibroblasts in the 3D experimental skin model for four days. From this point we refer to keratinocyte and fibroblast cells as skin cells. All cells are initially placed onto the 3D experimental skin model as a monolayer, as uniformly as possible. MTT (3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) assays highlight the metabolic activity of all cells, and show the spatial extent and spatial structure of cells on the top surface of the 3D experimental skin model. Images in Fig. 3a-b show prominent dark purple clusters on the surface of some 3D experimental skin models. Control studies, where 3D experiments are constructed without melanoma cells, show a complete absence of nests [see Additional file 1] suggesting that the dark purple clusters in Fig. 3a-b are melanoma nests. We make the natural assumption that higher densities of metabolically active cells are associated with darker purple colouration.
Images in Fig. 3a show that larger nests are associated with higher initial numbers of melanoma cells. To quantify this we measure the area of individual nests using ImageJ [15], and data in Fig. 3c confirms our visual observation. Interestingly, larger initial numbers of melanoma cells lead to a smaller number of larger nests [see Additional file 2]. This is consistent with smaller sized nests coalescing into a smaller number of larger nests over time. These results suggest smaller nests might coalescence into larger nests over time. To confirm this conjecture we would need to analyse our experiments using time-lapse imaging. Since our results show that cell number plays a critical role, we now examine the role of proliferation by suppressing mitosis.
We examine the role of cell proliferation by constructing 3D experimental skin models with irradiated melanoma cells. Images in Fig. 3b show that this leads to the formation of dramatically smaller nests. To quantify our results, the area of individual nests is measured using ImageJ [15] [see Additional file 2]. Data in Fig. 3d shows a similar trend to data in Fig. 3c as the nest area increases with initial cell number. However, comparing results in Fig. 3c-d shows that proliferation plays a dominant role in nest formation. For example, experiments initialised with 8500 proliferative melanoma cells leads to a median nest area of 0.15 mm 2 , whereas the median nest area is just 0.027 mm 2 when proliferation is suppressed. These measurements of nest area do not provide direct estimates of the number of cells present in each nest. However, it is reasonable to assume that larger nests contain more cells than smaller nests.
Our results are different to previous 3D studies that show melanoma nests are formed by cell migration [5]. We anticipate that the difference in our outcome could be due to: (i) differences between the melanoma cell lines used; (ii) the interaction of melanoma cells with the surrounding skin cells in the 3D experiments; or, (iii) differences in the material used to construct the 3D model described in [5] and the 3D model used in this study. Since our experiments are performed in 3D materials derived from human skin, and our experiments involve culturing melanoma cells together with primary human skin cells, we feel that our results are more realistic than examining nest formation in monoculture experiments in Matrigel. We now perform immunohistochemistry to confirm that irradiated melanoma cells survive in the 3D experimental human skin model over a period of four days.

Irradiated melanoma cells survive in a 3D experimental skin model
Here, we perform a series of experiments using a specific melanoma marker to provide additional evidence that nests observed on the 3D experimental human skin models are clusters of melanoma cells, and that irradiated melanoma cells survive in a 3D environment over four days. The 3D experimental skin models are constructed using both irradiated and non-irradiated melanoma cells. Vertical cross-sections through the 3D experimental skin models initialised with melanoma cells are stained using S100, which is a reliable melanoma cell marker [17]. Both irradiated and non-irradiated melanoma cells are found in the 3D experimental skin model after four days. Images in Fig. 4a-f show positive S100 staining of melanoma cells. In particular, Fig. 4b-f show positive S100 staining of irradiated melanoma cells after four days. This immunostaining confirms that irradiation does not alter the antigen properties of melanoma cells for this marker, and the irradiated melanoma cells survive in a 3D experimental skin model for four days. Our experimental results use skin cells and skin dermis from one donor. Additional results using cells and dermis from two other donors show little variability between them.

Variability between skin samples
We now examine whether there is any important variability in our results between skin samples from different donors. To examine this we perform additional experiments using dermis and primary skin cells from three different donors, which we denote as donor A, donor B and donor C. We show MTT assays on the 3D experimental skin models initialised with non-irradiated and irradiated melanoma cells in Fig. 5. The upper row of images in Fig. 5a-c show 3D experimental skin models initialised with 1250, 5000 and 8500 non-irradiated melanoma cells, respectively. In each case, we see that larger nests are associated with higher initial number of melanoma cells. A similar trend is observed for the images in  [19]. Simulations with suppressed melanoma proliferation use P p (m) = 0 the lower row of images in Fig. 5a-c where the experiments are initialised with an equivalent number of irradiated melanoma cells. However, regardless of whether we consider results from donor A, donor B or donor C, we always see that nest formation is dramatically reduced when we consider irradiated, non-proliferative melanoma cells.
Visual inspection of the images in Fig. 5 suggests that the size, shape and number of individual nests does vary slightly between the three donors. However, the influence of the initial cell number and the action of cell proliferation on nest formation remains consistent between the skin samples from the three different donors. That is, larger initial numbers of cells produces larger nests, and the action of cell proliferation leads to dramatically larger nests. To provide additional evidence we also measure the area of individual nests on skin samples from all donors using ImageJ [15]. Data provided [see Additional file 2] confirm that the relationship between initial cell number and the action of cell proliferation holds for all three donor samples.
The nests on the 3D experimental skin model initialised with 1250 irradiated melanoma cells are very small. Most experimental replicates of this particular experiment do not lead to any visually observable nests, as shown in the lower row of images in Fig. 5a. Therefore, data for nest area in these experiments is omitted [see Additional file 2]. We now use an IBM to verify our experimental outcomes.

Modelling melanoma nest formation using an individual based model
To corroborate our experimental findings, we use a random walk-based IBM to simulate the key features of the experiments. The IBM describes the spatial distribution of simulated cells on a 3D square lattice [18]. We use a 3D lattice of cross section 3 mm × 3 mm, and depth 2 mm, to represent the central region of each experimental 3D skin model (Fig. 6a). The lattice spacing is 20 μm. Simulated cells are called agents. We consider non-adhesive skin agents (green, Fig. 6b) and adhesive melanoma agents (blue, Fig. 6b). Note that the domain shown in Fig. 6b is a small subregion within the overall domain so that we visualise just the upper portion of the lattice where the majority of agents are located. It is well-known that it can be difficult to quantitatively calibrate stochastic IBMs to match complicated experimental data precisely [19,20]. Therefore we simply use parameters in the IBM that are adapted from previous work [9,21]. These previous studies report estimates of the proliferation rate of SK-MEL-28 melanoma cells, the proliferation rate of primary human fibroblast cells, the cell diffusivity of SK-MEL-28 melanoma cells and the cell diffusivity of primary human fibroblast cells [21]. We make a reasonable assumption that the proliferation rate and cell diffusivity of keratinocyte cells are the same as the proliferation rate and the cell diffusivity of the fibroblast cells, respectively. Our estimate of the strength of cell-to-cell adhesion is also adapted from a previous study where this parameter was determined using a series of two-dimensional barrier assays with a metastatic melanoma cell line [9]. This approach of using previously-reported parameter estimates allows us to focus on understanding the roles of the key underlying biological features, such as the role of cell migration and cell proliferation, without being distracted by the secondary task of obtaining precise parameter estimates. We achieve this by using previously determined parameter estimates and simply comparing simulation results where melanoma cell proliferation is present, with simulation results where melanoma cell proliferation is suppressed.
We initialise the IBM simulations to precisely mimic the way that cells are placed onto the upper surface of the 3D skin in the experiments. To initialise the simulations we randomly place a particular number of skin and melanoma agents onto the surface of the 3D lattice. The initial number of agents in each subpopulation is chosen to match to the initial cell density in the experiments. Figures 6b-d show smaller sub-regions of the 3D simulated skin to visualise the distribution of agents on the 3D lattice as clearly as possible. Results in Fig. 6b-c show that the IBM predicts the formation of clusters of adhesive melanoma agents on the surface of the 3D lattice. Results in Fig. 6d shows how the IBM predicts the downward movement of both skin and melanoma agents. Fig. 6d shows that skin agents move deeper into the 3D lattice than the melanoma agents, while nests of melanoma agents tend to remain on the surface. Overall, the spatial arrangement of skin and melanoma agents in the IBM (Fig. 6b-d) is similar to the spatial arrangement of cells in the 3D experiments (Figs. 3 and 4) [6].
To explore the role of initial melanoma cell number in nest formation, IBM results in Fig. 6e show that nests form on the surface of the 3D lattice, and that the trends in simulated nest area are qualitatively similar to those in the corresponding experiments. Therefore, the simulation outcomes in Fig. 6e confirm that initial melanoma cell number is an important factor in driving nest formation. We also explore the role of cell proliferation by repeating the simulations in Fig. 6e without any melanoma agent proliferation. Simulation results in Fig. 6f are comparable to the corresponding experimental results, as we observe similar trends in nest size and morphology. In conclusion, similar to the experiments, our 3D simulation results indicate that melanoma nest formation is driven by initial melanoma cell number, and that the presence of melanoma proliferation leads to dramatically-larger nests.
In addition to qualitatively visualising the trends in Fig. 6, we also use the IBM results to quantitatively examine trends in simulated nest size. Boxplots in Fig. 7 show data quantifying the size of nests predicted using the IBM under four different conditions. We measure the area of individual nests in the IBM using the Image Region Analyzer in MATLAB [22]. For model realisations where nests are not clearly defined we adjust the image manually by increasing the separation between neighbouring nests so that the Image Region Analyzer accurately measures nests separately. We exclude extremely small nests that are composed of less than four agents.
Results in Fig. 7a confirm that larger initial number of melanoma agents leads to larger simulated nests. Results in Fig. 7b show that suppressing proliferation in the IBM leads to dramatically smaller nests. These results in Fig. 7a-b correspond to the experimental results in Fig. 3c-d. Both the boxplots in Figs. 3c-d and 7a-b report nest area in the same units, therefore this is a direct comparison of the experimental observation and the prediction of the computational model.
In addition to using the IBM to replicate the experimental results, it is also straightforward to adjust the parameters in the IBM to make some simple predictions that have not been experimentally validated. Additional results in Fig. 7c show the distribution of simulated nest size when the proliferation rate of melanoma cells is reduced by half. Noting the difference in the vertical scale in Fig. 7a and b, we see that reducing the proliferation rate of melanoma cells by half leads to a reduction in simulated nest size by a factor of ten. Similarly, additional results in Fig. 7d show the distribution of simulated nest size when the strength of cell-to-cell adhesion for the melanoma cells is reduced by half. Again, noting the difference in vertical scale in Fig. 7a and d shows that reducing the strength of melanoma adhesion by half reduces the size of simulated nests be a factor of ten.

Conclusion
Our combined experimental and simulation findings demonstrate that cell proliferation plays the dominant role in melanoma nest formation. While it is wellaccepted that proliferation is important in the latter stages of tumour growth [23] and in the spatial spreading of cell populations [24], our work shows that proliferation is vitally important at the very earliest stages of melanoma progression. As far as we are aware, our work if the first to use a 3D experimental human skin model incorporating irradiated and non-irradiated melanoma cells and shows that cell proliferation is the dominant mechanism that drives melanoma nest formation.
Our results, pointing to the importance of cell proliferation, are interesting for a number of reasons: (i) previous monoculture experiments report that melanoma nests are formed by cell migration in Matrigel [5]. One potential explanation for this difference is that the Matrigel experiments are very different to our experiments since we study nest formation on 3D human tissues where melanoma cells are in contact with skin cells; (ii) some previous mathematical models of cluster/ nest formation focus on cell migration only, e.g. [25], whereas we find that cell proliferation plays the most important role; and (iii) our findings about the importance of cell proliferation in melanoma progression are consistent with the fact that many promising melanoma drugs aim to suppress proliferation [26][27][28].
Our suite of experimental results can be extended in many ways. For example, one limitation of our work is that we group the keratinocytes and fibroblast cells together, and refer to these cells as skin cells. It would be interesting to repeat our work and use specific markers to differentiate between these two populations of skin cells [29]. Another interesting extension of our experimental work could be to examine nest formation in 3D experiment using a mixture of irradiated and non-irradiated melanoma cells. This condition could mimic a partial reduction in proliferation, whereas our results correspond to a total inhibition of melanoma cell proliferation. Additionally in these experiments, cell proliferation can be blocked using drug treatments, such as mitomycin-C or some other commercially available proliferation inhibitor. Another relevant extension could be to perform a series of 3D skin experiments where the melanoma cells are treated so that they are non-migratory but maintain their ability to proliferate. Finally, it could also be interesting to repeat the 3D skin experiments as described here, and to image the formation of nests on a much shorter timescale that is comparable to the timescale of cell migration. If we had access to such time course data, it might then be possible to compare this kind of transient data from the experiment with transient information from a mathematical model [30].
Our suite of modelling results can also be extended in many ways. In this work we choose to work with a relatively simple mathematical model that represents just the key processes of interest, namely a population of motile, proliferative and adhesive melanoma cells, and a population of motile and proliferative skin cells. This model is parameterised using previouslydetermined parameter estimates [9,21]. While our model is useful in that it can both recapitulate our experimental results as well as generating new predictions that could be verified or challenged in future experimental studies, it would also possible to repeat all our simulation results using a more complicated mathematical model. For example, other modelling approaches such as continuous-space lattice-free models [31] or discrete models with force potentials between agents [32] could also be used in this context. While it is always tempting to use a more complicated mathematical model that incorporates additional biological detail, this approach is limited in that using more complex models requires additional parameters. We avoid this situation by always working with the simplest possible mathematical model that describes just the key features of interest.

Fibroblast isolation and culture
Human fibroblast cells are isolated following protocols in Haridas et al. [17]. Primary fibroblast cells are cultured at 37°C, in 5% CO 2 and 95% air.

Melanoma cell culture
The human melanoma cell line SK-MEL-28 is cultured as described in Haridas et al. [17]. SK-MEL-28 melanoma cells are kindly donated by Professor Brian Gabrielli (Mater Research Institute-University of Queensland). Cells are cultured at 37°C, in 5% CO 2 and 95% air.
A batch of SK-MEL-28 melanoma cells is irradiated to prevent cell proliferation. Approximately 1 × 10 7 melanoma cells are gamma-irradiated using a Gammacell 40 research irradiator (Australia) at approximately 0.8 Gy/ min for one hour resulting in a cumulative dose of 50 Gy. We refer to these non-proliferative cells as irradiated melanoma cells, and the proliferative cells as nonirradiated melanoma cells.
Identification of SK-MEL-28 cells is validated using short tandem repeat profiling (Cell Bank, Australia. January 2015).

Barrier assay
We perform circular barrier assays to observe and measure the spreading of populations of irradiated and non-irradiated melanoma cells. The protocol from Simpson et al. [16] is followed. Briefly, sterile stainless steel silicon barriers (Aix Scientific, Germany) are carefully placed in a 24-well tissue culture plate with 0.5 ml growth medium. The tissue culture plate containing cells is incubated for one hour at 37°C, in 5% CO 2 and 95% air. Viable cell suspensions of 20,000 cells/100 μl of irradiated and non-irradiated melanoma cells are carefully introduced into the barriers to ensure an even distribution of cells. The tissue culture plates containing cell suspensions are incubated for a further two hours to allow cells to attach to the plate. The barriers are removed and the cell layers are washed with serum-free medium (culture medium without foetal calf serum) and replaced with fresh growth medium. Plates are then incubated at 37°C, in 5% CO 2 and 95% air for zero, two and four days. We replace the growth medium after two days to replenish the nutrients. Each assay is performed in triplicate.

Crystal violet staining
We use the staining technique described by Simpson et al. [16] to analyse the barrier assays. In brief, cell monolayers are washed with phosphate buffered saline (PBS; Thermo Scientific, Australia) and fixed using 10% neutral buffered saline (United Biosciences, Australia) for 20 min at room temperature. The fixed cells are stained using 0.01% v/v crystal violet (Sigma Aldrich, Australia) in PBS for 20 min at room temperature. Excess crystal violet stain is removed using PBS, and the plates are air-dried. Images of irradiated and nonirradiated cell populations are acquired using a Nikon SMZ 800 stereo microscope fitted with a Nikon digital camera.

Establishing 3D experimental skin model with melanoma cells
We establish 3D experimental skin models using the skin collected from donors undergoing elective plastic surgery. The protocol for establishing the 3D skin equivalent model with melanoma cells is adapted from previous work [7]. In brief, sterile stainless steel rings (Aix Scientifics) with a radius of 3 mm are placed on the papillary side of the de-epidermised dermis in a 24-well tissue culture plate (Nunc®, Australia). We refer to the de-epidermised dermis as dermis. Single cell suspensions of primary keratinocyte cells (20000), primary fibroblast cells (10000) and non-irradiated melanoma cells (1250; 5000; 8500), are seeded onto the dermis in full Green's medium as uniformly as possible, and incubated at 37°C, in 5% CO 2 and 95% air for two days. We refer to the primary keratinocyte and fibroblast cells as skin cells. Subsequently, the stainless steel rings are removed and the dermis containing cells is submerged in full Green's medium for a further two days. After this four-day pre-culture period, the spatial distribution of cells in the 3D experimental skin model is analysed. We also perform a series of equivalent experiments using irradiated melanoma cells.
All experiments are performed in triplicate. Furthermore, all experiments are repeated using primary skin cells and dermis from three separate donors to account for variability between different donors.

MTT assay
An MTT (Thermo Scientific) assay is performed to check the metabolic activity of cells on the 3D experimental skin models. These assays are imaged with a stereo microscope (Nikon SMZ 800) fitted with a Nikon digital camera. We follow the protocol from Haridas et al. [7].

Immunohistochemistry on 3D experimental skin models with melanoma cells
We use immunohistochemistry to identify melanoma cells in the 3D experimental skin models. 10% neutral buffered formalin (United Biosciences, Australia) is used to fix the 3D experimental skin models. The tissue is divided through the centre of the MTT positive region using a sterile blade. The two smaller pieces of tissue are processed and embedded in paraffin. These samples are sectioned into 5 μm thick sections using a microtome. These sections are de-paraffinised, rehydrated and then subjected to heat-mediated antigen retrieval treatment using sodium citrate buffer (pH 6.0) in a decloaking chamber (Biocare Medical, USA) at 95°C for 5 min. Skin sections are washed in PBS followed by immunostaining using the MACH 4™ Universal HRP polymer kit (Biocare Medical). The primary antibody S100 (Dako, Australia) is diluted in DaVinci Green diluent (Biocare Medical) at 1:3000, and these sections are incubated with the primary antibody for one hour at room temperature. Positive immunoreactivity is visualized using 3,3-diaminobenzidine (DAB; Biocare Medical) and then counterstained with using Gill's haematoxylin (HD Scientific, Australia). The sections are dehydrated, and mounted on coverslips using Pertex® mounting medium (Medite, Germany). All stained sections are imaged using an Olympus BX41 microscope fitted with an Olympus digital camera (Micropublisher, 3.3RTV, QImaging; Olympus, Q-Imaging, Tokyo, Japan).

IBM simulation methods
We use a 3D lattice-based IBM, with adhesion between some agents, to describe the 3D experiments. In the IBM, cells are treated as equally sized spheres, and referred to as agents. These agents are restricted to reside on a 3D square lattice, with no more than one agent per site. The lattice spacing, Δ, represents the approximate size of each simulated agent, or the minimum spacing between agents. Here, we set Δ = 20 μm to match previous measurements [21]. We use a 3D lattice of dimension, 3 mm × 3 mm, and depth 2 mm, to represent the central region of each experimental skin model. This means that the number of lattice sites is 150 × 150 × 100. We choose the depth of the domain so that agents in the simulation never touch the bottom of the domain during the simulations. The parameters in the simulation model are adapted from previous studies [21]. Since we use the 3D lattice to represent the central region of the tissue, where cells are initialised uniformly across the surface, we apply periodic boundary conditions along all vertical boundaries. Since cells cannot leave the skin through the upper or lower surfaces, we apply no flux conditions on the upper and lower horizontal boundaries of the 3D lattice. We choose the depth of the 3D lattice to be large enough so that the agents never touch the bottom boundary of the lattice on the time scale of the simulations we consider.
To initialise simulations, we randomly place a particular number of simulated skin agents, N (s) (0), and a particular number of simulated melanoma agents, N (m) (0), onto the surface of the lattice. When the IBM is initialised we take care to ensure that no more than one agent occupies each lattice site. We always choose the initial number of agents in each subpopulation to match the equivalent initial density of cells in the experimental skin model. In the experiments, the initial populations of cells are uniformly placed inside a disc of radius 3 mm, whereas in the IBM the initial populations of agents are uniformly placed inside a square subregion of side length 3 mm. We set the initial number of skin agents to be N (s) (0) = 9549 to match the initial experimental population of 30,000 skin cells distributed in a disc of radius 3 mm. We vary the initial number of simulated melanoma agents to be N (m) (0) = 398, 1592 or 2706, to match the initial density of melanoma cells. This initial experimental density corresponds to 1250, 5000 and 8500 melanoma cells distributed in a disc of radius 3 mm. To match the experiments, the IBM simulations are run for four days.
At any time, t, there are N(t) = N (m) (t) + N (s) (t) agents on the lattice. In each discrete time step, of duration τ, we use a random sequential update method [34] to simulate motility and proliferation events. The algorithm involves executing two sequential steps: 1. N(t) agents are selected one at a time, with replacement and given the opportunity to move to a nearest neighbour lattice site with probability P m (s) ∈[0,1] and P m (m) ∈[0,1]. Here we can specify different motility probabilities for the skin cells and the melanoma cells, and this is important because previous work has shown that fibroblast cells are more motile than melanoma cells [21]. If the chosen agent is a melanoma agent, we incorporate adhesion into the model by examining the occupancy of the 26 nearest lattice sites in the 3D Moore neighbourhood. We count the number of those sites occupied by melanoma agents, a [18]. Potentially motile melanoma agents then attempt to move with a modified probability, P m * = (1 − q) a , which accounts for adhesion between neighbouring melanoma agents. The parameter q controls the strength of melanomamelanoma agent adhesion, with q = 0 corresponding to no adhesion, and increasing q leading to increased adhesion [18]. Setting q = 1 corresponds to maximal adhesion, and this would prevent any motility of melanoma agents that are in contact with other melanoma agents. We do not include any adhesion between skin agents as fibroblast cells are known to be mesenchymal and act as individuals rather than being strongly affected by adhesion [16]. If a movement event is successful, the agent attempts to move to a nearest neighbour lattice site from the six sites in the 3D von Neumann neighbourhood. To simulate crowding effects, potential motility events that would place an agent on an occupied site are aborted. 2. N(t) agents are selected one at a time, with replacement and given the opportunity to proliferate with probability P p (s) ∈[0,1] and P p (m) ∈[0,1]. Again, this framework allows us to specify different