Stochastic model of contact inhibition and the proliferation of melanoma in situ

Contact inhibition is a central feature orchestrating cell proliferation in culture experiments; its loss is associated with malignant transformation and tumorigenesis. We performed a co-culture experiment with human metastatic melanoma cell line (SKMEL- 147) and immortalized keratinocyte cells (HaCaT). After 8 days a spatial pattern was detected, characterized by the formation of clusters of melanoma cells surrounded by keratinocytes constraining their proliferation. In addition, we observed that the proportion of melanoma cells within the total population has increased. To explain our results we propose a spatial stochastic model (following a philosophy of the Widom-Rowlinson model from Statistical Physics and Molecular Chemistry) which considers cell proliferation, death, migration, and cell-to-cell interaction through contact inhibition. Our numerical simulations demonstrate that loss of contact inhibition is a sufficient mechanism, appropriate for an explanation of the increase in the proportion of tumor cells and generation of spatial patterns established in the conducted experiments.

In this paper we present theoretical and experimental frameworks designed to investigate the modulated contact inhibition and its role in the formation of a carcinoma in situ and to demonstrate that allelophilic properties of cancer cells is a key feature for their uncontrolled proliferation.

Keratinocytes and melanoma cells co-culture proliferation.
To evaluate the cell proliferation, the human metastatic melanoma (SK-MEL-147) and human immortalized keratinocytes (HaCaT) cell lines were selected for co-culture experiments. The choice of these cells allows us to mimic the interaction between the skin basal layer cells and the melanoma. Another reason for selecting these cell lines was to compare the co-culture development with patterns produced through a stochastic model dynamics. The latter involves a cell line that shows a distinctive degree of contact inhibition (a property of HaCaT) and another cell line that is highly tolerant, i.e., displays a loss of contact inhibition (which is a characteristic of SK-MEL-147). In the supplementary material we collect results of our experiments. At the post confluence stage the carrying capacity of HaCaT is at 1779.56 ± 130.47 cells/mm 2 while for SK-MEL-147 it equals 5043.51 ± 316.47 cells/mm 2 (see Data Fitting section and Fig. S1 at the Supplementary Information file). This demonstrates higher density levels achieved by melanoma cells (Fig. 1) confirming their distinctively lower degree of contact inhibition in comparison with keratinocytes. A similar phenomenon was observed in a different situation in ref. 4.
At the initial stage of the co-culture experiments, cells were seeded at 250 cells/mm 2 , at a proportion of keratinocytes to melanoma of 10:1, in a monolayer on a 24-well plate dish with coverslips. The co-culture was allowed to proliferate for eight days. The monolayer structure enabled us to investigate the role of contact inhibition in the cell proliferation at a quantitative level. After four days in the co-culture, cells reached confluence, and it was possible to observe the formation of growing melanoma clusters. These clusters are constrained by layers of keratinocytes cells, of density somewhat higher than normal (Fig. 1A). To evaluate the cell population growth, we counted the number of cells in images from 30 locations on the plate for each day of experiment. The obtained data were fitted by using the logistic growth model (Fig. 1B). The parameter ρ indicates the cell population growth rate, the maximum population density is denoted by K, and the time when the cell density equals K/2 is τ. The values for the fitting for each cell line are shown in Fig. 1D. Note that ρ and τ can be made (approximately) the same for both cell lines, while the ratio between the maximum densities is ~4. The change in time of the ratio between the two cell population densities is shown in Fig. 1C. One may also note that the proportion of HaCaT cells density decays from ~10:1 in the first-day measurement to ~4:1 at the eighth day. Although the logistic model gives a proper fit for the experimental data, it yields no insight on the mechanisms underpinning the increase on the proportion of melanoma cells. To this end we propose a more detailed model aiming at providing a qualitative description of the observed data.
A stochastic model of cell proliferation for a cell-to-cell interaction through the contact inhibition. A mathematical model explaining the above experimental results is related to the Widom-Rowlinson model [24][25][26] . Since the data refer to a monolayered co-culture arrangement, we can focus on a two-dimensional system composed of two cell types interacting as colored billiard balls of different sizes, where each cell type has two exclusion areas surrounding it (one for balls of the same color, the other for balls of different colors). Here type 1 are tumor/cancer cells (SK-MEL-147) and type 2 are healthy/normal cells (HaCaT). The cells are placed at the vertices of a two-dimensional grid Λ formed by unit squares. The distance between two vertices is represented by the minimal number of unit edges connecting them (a graph metric). An admissible configuration of this system satisfies the following rules: (i) each vertex is occupied by at most one cell; (ii) the distance between two cells i and j, with i, j = 1, 2, is a positive integer not less than a given value D(i, j), referred to as exclusion distance/ diameter between cell types i and j. In fact, D(i, j) represents the minimal allowed distance between sites occupied by cells of types i and j.
A formal mathematical description of the model and related results are contained in the supplementary information (see Theorem 1).
The experiments show a close relationship between the degree of contact inhibition among cells and the cell density. It indicates a possible biological interpretation of self-exclusion diameters D(i, j): they measure the degree of contact inhibition of melanoma and keratinocytes, respectively. A greater value of the exclusion diameter D(i, i) indicates a higher degree of contact inhibition. Intuitively, it can be understood in terms of dense-packing the grid Λ by a single cell type. As the self-exclusion diameter of the cell increases, the maximal proportion of occupied vertices reduces. In other words, the cell density within the grid is inversely proportional to D(i, i).
For the simulations, the minimal allowed distances were chosen to be D(1, 1) = 1, D(1, 2) = D(2, 1) = 3, and D(2, 2) = 2; this indicates that the smallest exclusion occurs between two tumor cells. Pictorially, a pair of healthy cells are twice less allelophilic and a healthy-tumor pair is three times less allelophilic than the pair of tumor cells. A diagram presenting features of an admissible configuration is shown on Fig. 2A: here normal (tumor) cells are represented as red (blue) colored circles. Note that our choice of D(1, 1) allows the tumor cells to occupy nearest-neighbor vertices. The Fig. 2B shows a magnification of an admissible configuration on a 25 × 35 grid as obtained in the course of a simulation. The spatial pattern is characterized by the formation of clusters of tumor cells surrounded by a sea of normal cells separated by a layer of empty vertices. Such a pattern resembles the keratinocytes and melanoma spatial arrangement observed experimentally, as shown in Fig. 1A. As was stated, the exclusion diameter D(i, j) is interpreted as a degree of contact inhibition between cell types i and j: the smaller D(i, j) is the lower the degree of contact inhibition.
Stochastic dynamics of cell proliferation: a Markov chain on admissible configurations. We consider three possible transitions for the i-th cell type: division, migration, and death, occurring at a rate α i , δ i , and ρ i , respectively, with i = 1, 2. For simplicity, quiescence of cell type i is considered to happen at rate α i but a general version of the model might take different values for that rate. Cell transformations from normal to tumor is neglected due to the time scale of the experiments, and because cell lines used in experiments are from different embryonic origins. We treat the randomness of cell proliferation in terms of a discrete-time Markov chain on admissible configurations on a square grid of size L × L where L = 200. A given vertex of the grid can assume a state 0, 1, or 2 indicating that it is empty (0), occupied by a tumor cell (1), or by a normal cell (2).
The transition from a current state towards the next one begins with selecting a vertex x ∈ Λ with probability L −2 and verifying its state as occupied or empty. (i) If occupied by a cell type i, vertex x may remain occupied (quiescence) by it with probability α i /Q or becomes empty with probability is the number of empty vertices at distance D(i, i) from x, which can receive the migrated cell in accordance with the admissibility rule. In case of migration (occurring with The red/blue circles indicate the normal/tumor cells placed at sites/vertices of a square grid. The four blue circles positioned at neighboring vertices represent a part of a densely-packed configuration for tumor cells. The four red circles surrounded by diagonal red squares represent a part of a checker-board densely-packed configuration for normal cells: their inhibition level prevents any further concentration. Each diagonal red square around a red circle marks an exclusion area for other red (healthy) cells. The larger purple region around the red circles covers the vertices forbidden for blue (tumor) cells in the vicinity of healthy ones. The gray circles mark two randomly selected empty vertices, located in a shadowed and non-shadowed region, respectively. The arrows from circles 1 and 2 point to possible states at the next transition: (i) vertex 1 may remain empty or become occupied by a red cell (but not blue), since it is located in a purple shadowed area; (ii) vertex 2 is in a non-shadowed area and may remain empty or become occupied by either a red or blue cell. (B) A simulated configuration: blue (red) dots indicate tumor (normal) cells. After a sufficiently long time, a spatial pattern is formed: normal cells surround tumor cell clusters, with empty separating layers. Note similarities with experimental images in , and it becomes occupied by the migrating cell. (ii) If vertex x is empty, it remains empty with probability (ρ 1 + ρ 2 )/R, and notice that this probability might be defined in terms of different rates and this choice was done aiming at avoiding introduction of more parameters. Next, a cell type i may occupy x with probability , and n i (x) stands for the number of vertices currently occupied by cell type i among the nearest allowed neighbors of x. The migrating cell is chosen with probability 1/n i (x) and its vertex becomes empty after migration. Further,  = i ( ) 1 or 0 when vertex x can or cannot be occupied by a cell of type i without violating the admissibility condition. If R = ρ 1 + ρ 2 , i.e. vertex x cannot be occupied by any type of cell, it remains empty.
The proposed structure of Markov chain transitions gives a rather coarse-grained (and simplified) approximation to possible dynamics of the co-culture cell population. In particular, the outcomes of cell division take into account geometric characteristics of the current cell configuration only partially. (Through numbers n x ( ) i e ( ) and n i (x)). It is not clear to what extent these simplifications change the overall picture of the co-culture evolution; we consider them as an initial step in the study of mathematical models of such an evolution. Nevertheless, even such a simplified model shows a similarity between curves in Figs 1 and 2.
In the stochastic dynamics, tumor cells increase their proportion due a low degree of contact inhibition. The results of simulations based on the stochastic dynamics are analyzed in terms of the i-th cell type density, that is, the proportion of vertices of the grid which are occupied by the i-th cell type. Figure 2C shows the dynamics of the normal and the tumor cell populations. The cell densities are shown on the vertical axis and the time steps of the Markov chain are marked in a logarithmic scale on the horizontal axis. Here one should not confuse the Markov time steps with experimental time scales. In the co-culture experiment there are multiple cell state transitions going on in parallel while in the simulation we take at most one transition per a single time step. An initial configuration is chosen with a the type 1 (tumor) cell density 0.01, type 2 (normal) cell density 0.09 and the total cell population density 0.1. The transition rates are made the same for both cell types, and only the self-exclusion diameters D(1, 1) = 1 and D(2, 2) = 2 are different, with normal cells having the self-exclusion diameter twice as long as tumor cells. (The choice of the cross-exclusion diameter D(1, 2) = D(2, 1) = 3 seems to play a lesser role here.) We avoid border effects by maintaining them empty during an entire simulation. The red (or blue) lines dots indicate the normal (or tumor) cell population density. During the time interval of the simulations there is a variation in the ratio of normal to tumor cell-densities, as presented on Fig. 2D. Since all proliferation rates of both cell types are the same, such a change should be attributable to the difference on the values of their exclusion diameters. In biological terms, it means that a low degree of contact inhibition favors the tumor cells population to increase its proportion within the total cell population.  The pictures on frames E and F, do not conform with the co-culture experiments presented here. Their presence has two purposes, the first is to show an asymptotic behavior of our model and how one might expect a co-culture experiment to behave when the only difference between the tumor and normal cells was their degree of allelophilia. A disagreement between experimental results and simulations can be explained by the fact that the experiments have been conducted for relatively short time-periods. Also, the proposed stochastic model might miss some important aspects of the co-culture development (e.g., because the two cell types are different on more aspects than their degree of contact inhibition). See the Discussion section of the paper.

Mutual contact inhibition promotes the formation of cell clusters.
A robustness analysis of prevalence of the most tolerant cell type. Figure 4 shows results of an analysis of prevalence of the most tolerant species after 3 billion iterations. Frame A indicates the average density at the final iteration after 100 repetitions on the y axis and multiple values for the ratio α 2 /α 1 , with error bars included for each value of the ratio. The average densities are approximately the same for the two cell types when α α ∼ / 6 2 1 and cell type two prevails for α 2 /α 1 = 10 at the final iteration. Frame B comprises three samples of the Markov chain simulation displaying the cell densities when α 2 /α 1 = 6. The three samples have a qualitatively similar behavior but each of them shows a difference in the degree of prevalence of the cell type 1 due to the randomness of simulations. The final spatial configurations for the three samples are shown in frames C, D and E. The different patterns can be explained by large error bars indicated on frame A. Frame C shows a coexistence-like pattern where type two cells surround cells of type one. Frame D displays a condition where type two cells still prevail while prevalence of type one cells is observed in frame E. Results in frames B and E indicate that for a larger amount of iterations one may expect the most tolerant cell type to prevail in a spatial domain; increasing the ratio α 2 /α 1 causes this phenomenon to take longer to happen. Frame F presents an analysis of melanoma cell clusters in the experiments and simulations. The top graphs show that the cluster areas in both simulations and experiments have similar box-plots. Here 1 pixel 2 of area of the heatmap obtained from the simulation is equated with ~93 μm 2 of the area in the experiment. The square root is used in the identification of the perimeter correspondence. Notice that there is a good agreement between the median of both aspect ratios and convex hull obtained from both (F) The analysis of cluster properties as obtained from the experimental data and simulations. We have obtained a good agreement between the box-plots for the cluster area and cluster perimeter when we take 1 pixel on a simulation heatmap to be ~10 μm. experiments and simulations. Furthermore, frame F shows that clusters have a diversity of shapes and sizes, randomly distributed over the plate and the simulation grid. This emphasizes stochasticity of the observed phenomenon.
The distance distribution quantifies the degree of allelophilia of the cells. To evaluate the degree of tolerance between cells of the same type, we performed an image-based analysis of the distribution of the distances between keratinocytes or melanoma cell lines in the co-culture. From the distribution of distances between cells we defined a minimal distance allowed for cells of the same type as the most frequent distance observed (a mode of the histogram). The distance distribution of keratinocytes assumes only distances appearing within the same range as those found in melanoma cell clusters. Three different distance distribution signatures have been observed (see Fig. 5). (a) A pre-confluence state where no significant difference has been detected between cells due to cell spread and low density (Fig. 5A). (b) At confluence, melanoma clusters are formed, and it is possible to detect a difference between histograms of distance distributions (Fig. 5B). (c) In a post-confluence state, a peak in the histogram for melanoma cells has been observed (Fig. 5C), showing a shorter distance more frequent than for keratinocytes. Tumors are denser than the surrounding normal tissue. This indicates a high degree of proximity and therefore a greater tolerance between cancer cells.
Keratinocyte density depends on the distance from melanoma. To evaluate melanoma proliferation restriction we quantified the density of keratinocytes in three distinct regions surrounding the tumor cells at 50, 100, and up to 150 μm and at the high density cancer cells domain (Fig. 5G). As one can see, the keratinocyte density is affected by the distance from the melanoma cell domain (Fig. 5H). This morphology of epithelial cells are observed when cells are at high contact of inhibition. The keratinocytes are affected by proximity to melanoma cells, and more distant cells as a transitional morphology towards more dispersed state. Such a pattern is repeated and is induced by the communication between the SK-MEL-147 and HaCaT cells. A change on the proliferation rates in co-culture experiments indicates that the interaction between the tumor and normal cells is affecting their growth capabilities. As our data show, cells rearrange their spatial patterns to ensure that an initially well-mixed population of two cell types turns into clusters of a single cell type, either normal or tumor.

Discussion
The novelty of this work is that quantitative co-culture experiments are combined with a specific mathematical model. The co-culture experiments demonstrate that the low sensitivity to contact inhibition is responsible for a decrease in proportion rates of keratinocytes to melanoma cells and formation of tumor-like melanoma clusters. The explanation of these features has been provided through a specific mathematical model.
Our aim was to understand the role of a low degree of contact inhibition on the development of cancer in situ. The melanocyte and keratinocyte interaction through a direct contact is an important mechanism controlling their proliferation 27 . We observed a decrease on the proportion of keratinocytes to melanoma in the course of time: see Fig. 1C. Furthermore, a development of spatial pattern was observed, progressing from a well-mixed configuration (Fig. 5A) towards dense melanoma clusters surrounded by spatially spread domains with a relatively high density of keratinocytes (Fig. 5C).
To explain the observed phenomena, we use a mathematical model with a particular form of interaction between cells of the same and different types through different degrees of contact inhibition. The melanoma cells are supposed to have a lower level of contact inhibition than the keratinocyte cells. A mathematical mechanism of a stochastic nature has been considered, based on a Widom-Rowlinson-type model and its modifications [24][25][26] . The dynamics is presented as a Markov chain on admissible cell configurations with specified exclusion diameters for the keratinocyte (healthy/normal) and melanoma (tumor) cells; in the simulations we use the same rates for supposed transitions for cells of both types: division, death, and migration. Numerical simulations have been performed, confirming a number of features observed in the co-culture experiments. Tumor cells have smaller exclusion diameters; according to our interpretation, they are less sensitive to contact inhibition than healthy ones.
A major question is about including mechanical aspects of the cell motions into the model. This approach have been employed in investigations of homeostasis and morphogenesis 28,29 , and has shaded light on how mechanic stress affects contact inhibition molecular mechanisms 30 . Such an implementation in our stochastic approach would need a serious re-vision of the model which could be a topic for future research. In the mechanical approach cell-to-cell interaction happens through'springs' which compression (dilation) controls cells growth rates 31,32 . Our approach corresponds to the limit for which the'spring' is fully compressed and it generates an infinite potential barrier, that is, the exclusion diameters are establishing a hard core potential for cell-to-cell interaction.
The proposed theoretical approach enabled us to reproduce a temporal increase on the melanoma proportion in a series of numerical simulations. A formation of spatial patterns on the two-dimensional grid has been also confirmed, albeit not at a desired level of accuracy. An existing agreement between the co-culture experiments and numerical simulations suggests that the development of a co-culture of melanoma and keratinocyte cells has a distinctive stochastic aspect which can be quantified. This may lead to a new methodology of prognosis for the development of carcinoma in situ, particularly formation of tumor-type clusters of melanoma cells and a collective reaction of keratinocyte cells to such a formation. Specific biochemical mechanisms of cellular interaction standing behind these phenomena need further studies. This would involve a variety of mathematical models modifying the present model. In particular, exclusion areas for healthy cells may vary and depend upon their distances from developing clusters of cancer cells.
Our characterization of the degree of contact inhibition suggests the use of the cells-to-cells distances to quantify sensitivity to contact inhibition. Tissues with tumors are denser than the normal ones [33][34][35] , and, thus, tumor cells are expected to be closer to each other forming geometrically packed arrangements. This is confirmed in Fig. 5F where we evaluate the distance distributions for keratynocytes and melanoma cells after confluence. The distribution of melanoma cell-to-cell distances is displaced to the left in comparison with the keratinocytes; this reinforces our interpretation of the exclusion diameters as a measure of the degree of contact inhibition of a cell.
The present experiments have been conducted with co-culture cell mono-layer. In the mathematical model it was interpreted as a stochastic dynamics of a two-dimensional admissible configuration on a grid. The passage to dimension three is mathematically possible (at a cost of an increase in the simulation time). We expect the qualitative features of the model presented here to be preserved. Such an extension would enable us to evaluate experiments on the formation of in vivo carcinoma in situ and to compare experimental results with simulations, aiming to establish a quantitative model to investigate carcinogenesis at its early stage. However, it would require new methods of marking and counting cells in the co-culture.
Contact inhibition and the necessity for its full characterization. The theoretical curve for the ratio between the tumor and normal cell density has a more pronounced shape in comparison with the experimental data (Figs 3B and 1C). This may be attributed to other effects controlling cell proliferation in the co-culture. A more precise correspondence between experimental and theoretical results might be established by engineering two cell lines with very similar proliferation, death, and migration rates but having different degrees of modulation of contact inhibition. This could have practical applications on the design of cancer treatment aiming to increase the degree of contact inhibition of tumor cells 4,36,37 .
The simulation based on a mathematical model of a realistic cell interaction is quite promising: the higher level of allelophilia of the tumor cells may be a key feature for their prevalence within tissues. This can be detected through exclusion diameters as measure of the degree of contact inhibition of cell proliferation. Other types of mutual influence could also be taken into account. The co-culture experiments show that the mechanism of a higher allelophilic degree might be due to molecular complexes surrounding the cells (see Fig. 5G and H). This shows the necessity of establishing a biochemical characterization of the exclusion diameters (or more factors) involved in the cell interaction.
A view in the general context of cancer research. Let us address some aspects of possible criticism of our approach. We recognise that an excessive use of analogies between biological phenomena taking place at different scales (e.g., tissue and ecological levels) is risky. On the other hand, the evolutionary paradigm is an established tool for understanding carcinogenesis, provided that a proper caution is taken 18,21,38 . In particular, while comparing ecological phenomena with those at a tissue level one must consider the specificity of cell-to-cell interactions and propose proper forms for quantitative modeling. An important mechanism here is the presence of membrane molecules causing the less immunogenic cancer cells to be selected for survival 39,40 . Another example is the development of cancer cells capable of performing glycolysis in hypoxic environment, when tumors reach a critical size interpreted as an adaptation 41,42 . Such a new feature of cancer cells causes local acidosis toxic to normal cells; this is considered as an analog of allelopathic interactions between plants or bacteriocins in microtube 43 . In contrast, the present study focuses on an impact of an allelophilic behavior upon the formation of a carcinoma in situ which, within the evolutionary paradigm, indicates that the selective advantage of cancer cells may happen because of their greater degree of allelophilia.

Methods
Cell culture. The Supplementary information file contains details on cell culture proliferation assay. Human immortalized keratynocites (HaCaT) and human metastatic melanoma (SK-MEL-147) cell lines were used in this study. Cells were grown in Dulbecco's Modified Eagle's Media (DMEM) supplemented with 10% fetal bovine serum (FBS) in tissue culture incubator (5% CO 2 at 37 °C).

CDH1 immunofluorescence.
HacaT and SK-MEL-147 cells were seeded over 24-well plate with coverslips and grown for eight days. Cells were fixed with ice-cold methanol: acetone (1:1) at 4 °C for 15 min. Cells were blocked with 5% bovine serum albumin (BSA) in phosphate-buffered saline for 1 h at room temperature. After blocking, cells were incubated with antibody against CDH1 (1:200 in 1% BSA, Transduction Laboratories) for 45 min at room temperature. After washing, the coverslips were then incubated with Alexa Fluor 546 secondary antibody (1:400, Invitrogen) and 50 μg/mL Hoechst 33258 (Invitrogen) in 1% BSA for 2 h at room temperature. Immunofluorescence images for each day of experiment (n = 30) were acquired with an Nikon Eclipse E600 microscope, under 20x objective. The cells were counted based on the number of nuclei present in an image and the distance between cells were measured based on the centroid point of the image of the cell nucleus. The nuclei count were made manually using Cell Counter plug-in in ImageJ software (See the Image Analysis section in supplementary information file for more details).

Cluster features.
To evaluate melanoma cell clusters features, we performed an analysis of 60 clusters at the 8th day of the experiment. Shape descriptions performed on images captured under 10x objective consist of area, perimeter, convex hull ((Area)/(ConvexArea)) and aspect ratio (Major axis/Minor axis). Features of the clusters of cells of type 1 presented by one hundred repetitions of the simulations were also measured in the same way in heatmaps of 10 5 iterations.
Data availability statement. The data generated during and/or analysed during the current study are available from the corresponding author on reasonable request.