Roughness and dynamics of proliferating cell fronts as a probe of cell–cell interactions

Juxtacellular interactions play an essential but still not fully understood role in both normal tissue development and tumour invasion. Using proliferating cell fronts as a model system, we explore the effects of cell–cell interactions on the geometry and dynamics of these one-dimensional biological interfaces. We observe two distinct scaling regimes of the steady state roughness of in-vitro propagating Rat1 fibroblast cell fronts, suggesting different hierarchies of interactions at sub-cell lengthscales and at a lengthscale of 2–10 cells. Pharmacological modulation significantly affects the proliferation speed of the cell fronts, and those modulators that promote cell mobility or division also lead to the most rapid evolution of cell front roughness. By comparing our experimental observations to numerical simulations of elastic cell fronts with purely short-range interactions, we demonstrate that the interactions at few-cell lengthscales play a key role. Our methodology provides a simple framework to measure and characterise the biological effects of such interactions, and could be useful in tumour phenotyping.


Pharmacological modulator titration
To determine which concentration of the inhibitors would be most effective without significantly affecting cell mortality, separate titrations were conducted for each of the compounds on fully confluent healthy Rat1 cell cultures. Starting with values commonly used in our previous experiments, the applied concentration was increased in pseudo-logarithmic steps (1x, 2x, 5x, 10x, 20x, etc), up to a killing concentration determined by visual inspection of the cells after 24 hours of exposure, in comparison with cells prepared at the same time under identical conditions without inhibitor exposure.
Inhibitor effects on cell growth, morphology, and mortality were evaluated based on the presence and extension of lamelipodia, the distance between cell nuclei, loss of confluence, and of course the amount of dead cells. For colchicine, the initial concentration of 0.32 µM did not affect cell mortality, but morphology effects of a more rounded and globular cell shape were beginning to be apparent. At higher concentrations, these morphology changes were more significant, but we also observed significant mortality increase above 0.625 µM, leading us to conclude that the original concentration of 0.316 µM was already at the optimal value. For the cytochalasin B, starting from a concentration of 0.1 µM, cell mortality increased sharply from less than 5% to more than 30% for 0.5 to 1 µM, and the colony was seriously depopulated (50% mortality at 2 µM). We therefore chose an effective concentration of 0.5 µM. The toxicity of meclofenamic acid was likewise easily decided, the mortality again increasing rapidly from less than 5% to more than 80% between 50 and 100 µM, leading us to chose an effective concentration of 50 µM. Finally, with forskolin, we did not observe significant effects on cell mortality. At 5 µM we began to see apparent elongation of the cells. These effects appeared to saturate at 10 µM, which was therefore the concentration chosen for the experiments.We generally observed that at concentrations approaching the killing concentration, cells appear to shrink and we begin to observe loss of confluence with increasing amounts of empty space between the cells. The cells also detached readily from the substrate surface, dying and floating above the rest of the colony.

Image processing
Each of the 1000 x 1000 px images acquired using the Leica DM6000 were indexed by the x, y position of the motorised sample stage (the vertical height z was kept constant once the desired focal plane was established), allowing us to stitch the individual images in the frame matrix into a full panorama of the proliferating front, as shown in Fig. 1(a), using an in-house developed Python algorithm. To fully discriminate the cell colony from the background the original greyscale fluorescence images are binarised using the Otsu algorithm and the brightness histogram for the panorama, as seen in Fig. 1(b), with the colony in white. As a result of varying GFP expression between cells, the binarisation still leaves some dark areas behind the cell front, necessitating a "hole filling" step carried out with the Scipy function, "binary_fill_holes" , which leads to the image shown in Fig. 1(c). In this image, we identify the cell culture as the largest group of connected white pixels, obtained using the Scikit "label" function. The set of x, y coordinates of the front position is finally obtained using the Scikit function "find_contours", and can be seen in Fig. 1(d), overlaid on the original fluorescence microscopy image the cell colony.

Monovalued approximation
The extracted front position follows the contours of the individual cells, and shows small overhang features where a single cell or small clusters of a few cells at the edge of the colony orient somewhat diagonally. However, the framework of disordered elastic systems, with an analysis of the correlation function B(r) of relative displacements ∆u(r), is only applicable in the case of a univalued interface, meaning that elements such as overhangs, holes, and islands cannot be treated within this approach. We therefore approximate the extracted front position by a univalued function. This correction is carried out using 3 different approaches, as shown in Fig. 2: 'outer' contour (yellow), choosing the highest value of u(z) if multiple points of the front are present at a given lateral coordinate z; 'inner' contour (green), choosing the lowest value of u(z) if multiple points of the front are present at a given lateral coordinate z; and 'sum' contour (blue), adding all the pixels of the colony at a given lateral coordinate z to obtain an effective univalued u(r).
As can be seen in Fig. 2(a), the roughness B(r) calculated for each of these corrected fronts shows small differences at sub-cell lengthscales, but these become insignificant by the time the behaviour at few-cell lengthscales is examined. Since the sum approach effectively smooths the interface, to preserve at least partially the effect of the overhangs, which are clearly a feature of the biological system, we chose to use the outer contour for all the analysis presented in this study.
This results in abrupt but still univalued 'jumps' in the front position. In order to estimate the impact of these jumps on the front geometry, we subtracted the outer and inner contours of the front, and extracted the statistics of the resulting deviations from a straight line. As can be seen in Fig. 3(a,b), the jumps, reflecting the original overhangs, are generally limited to at most a width of 8 µm (i.e. less than the characteristic cell size) and an average height of 200 µm (i.e. approximately 10-12 cells).

Rounding of the initial configuration
During the lift-off patterning, the flexibility of the silicone insert used to demarcate the initial position of the cell front in a confluent cell culture can lead to a rounded arc geometry rather than a straight line configuration. This arc corresponds to a segment of a circle with a radius of approximately 1 m, which can be quite precisely overlaid over the initial front position, as shown in Fig. 4(a), and contributes an identifiable artefact to the calculation of the front roughness B(r). As can be seen in Fig. 4, B(r) for such a 1 m radius circle arc shows a power law growth, with a roughness exponent ζ = 1. However, the initial roughness levels are over 4 orders of magnitude lower than that of the cell front itself, starting at approximately 2 × 10 −4 µm 2 at r = 1 µm. The roughness increases to 3 × 10 4 at r = 30mm, the highest available lengthscale of the analysis. While this artefact, when present, becomes the dominant contribution to the experimental observations above r = 2 mm, and affects to a lesser degree the roughness of the cell front above 500 µm, it does not affect our analysis of the double power law growth behaviour with very different roughness exponents at subcell (4-18 µm) and few-cell lengthscales (80-180 µm).

Reproducibility of roughness analysis
Our previous studies of physical interfaces have shown that sufficiently long individual fronts, or sufficient different iterations of shorter fronts are necessary to accurately extract the value of the roughness exponent ζ 1, 2 . The relatively large cell front panoramas obtained in our experiments (length L ∼ 2 16 px or ∼ 5cm) provide ample statistics for the analysis of the roughness B(r), restricted to r ≤ L 2 , and thus a high level of confidence in its accuracy. As a further test, and to check the reproducibility of the experimental protocol, we repeated each experiment three times under control conditions to obtain the time evolution of three separate cell front panoramas over 40 hours. For each of these realizations i (i = 1, 2, 3), we obtain the functions u i,t (r) describing the position of the cell-front after an evolution time t = 0, 4, 8, · · · 40 hours. We compute the roughness B i,t (r) for each of these fronts, and the average B t (r) = 1 3 ∑ 3 i=1 B i,t (r). over the three realizations i at the same time t. As shown in Fig. 4, the roughness does not appear to change significantly from one realisation to the other (apart from the rounding artefact present in one of the realisations, as discussed in Sect. 4), and after a long evolution time, the roughness obtained for the three realisations appears to converge to the average value.

Evolution to steady state roughness under control conditions
During our observations of their proliferation the cell fronts are maintained in-situ under the microscope, without changing the medium or adding nutrients. We therefore limited our discussions of the time evolution of the fronts in the main paper to 40 hours, to ensure that chemical depletion did not play a role in the observed cell behaviour. However, to establish when significant cell death and deterioration of the colony does occur, we also followed a smaller set of fronts for longer periods. In pharmacologically modulated cell fronts, this occurs already at 60-70 hours, while under control conditions cell fronts appear healthy and proliferate beyond 100 hours.
The results of our longest experiments can be seen in Fig. 6(a,b), following the average position and roughness, respectively, of a single front for 112 hours under control conditions. The front proliferates with decreasing speed over the full period of observation, approximately following a logarithmic time dependence, although increased acceleration is observed at 12-18 hours, and very little proliferation appears to occur beyond 100 hours. The roughness B(r) calculated at subcell lengthscales (r = 10 µm, blue) remains essentially constant throughout the measurement, suggesting a steady state which is already reached by the time the initial configuration of the experiment is established, and which could be expected for a roughness determined largely by intracellular interactions. At few-cell lengthscales (r = 100 µm, red), on the other hand, the roughness increases over the first 24 hours from an initial lower value reflecting the artificially flat configuration induced by the removal of the insert. After approximately 24-30 hours, the value of B(r) appears to also reach a steady state, with no significant increase, despite the fact that the front continues to proliferate with no noticeable deterioration or cell death.

Extracting the roughness exponent
One of the most important features of the roughness function B(r) is that it provides information about the interactions present in a system. To capture this information, two correlated elements are necessary. First, the spatial regions where the function The roughness B(r) calculated at r = 10 µm (subcell) and r = 100 µm (few-cell) lengthscales, in blue and red, respectively. A steady state with little increase in the roughness is observed already from the initial configuration at subcell lengthscales, and beyond approximately 30 hours at few-cell lengthscales.
follows a power law need to be identified, and second, the fit of the function with a power law should be done over this region to extract the roughness exponent describing the interactions at this lengthscale.
In the theory of disordered elastic systems, where, for example, a system evolves as result of the competition between the elasticity of the interface c, the disorder or pinning, and thermal fluctuations proportional to the temperature T , it is well known that the roughness function at short lengthscales is dominated by thermal fluctuations, and follows a power law behavior T c r 2ζ th , where the roughness exponent is given by the thermal exponent ζ th = 0.5. However, at large lengthscales B(r) is dominated by disorder, if present, and depending on the disorder type follows a different power law proportional to ∼ r 2ζ dis . For example, for random bond disorder the scaling is governed by ζ dis = 2 3 . In the case of biological systems, the nature of the interactions between cells and how they affect the dynamics and statics of a cell front is not yet well established, and is what motivates the current study. In order to apply the theory of disordered elastic systems to effectively capture the information about the interactions in a growing colony of cells, we identify the regions where B(r) follows a power law behavior as follows: we perform a weighted fit of log B(r) with a linear function of r over the region [log r i , log(r i + w i )], where r i is the starting point and w i is the length of the fitting region chosen according to the initial value r i , in order to always obtain a window of equivalent size on a logarithmic scale for each starting point (see Fig. 7). The weights used for the fit are computed as the uncertainties ∆y associated to each point of y = log B(r). These uncertainties can be estimated as ∆y = 1 B(r) ∆B(r), where ∆B(r) = std(all observed values) √ number of observations-1 , and the number of observations decreases when r is increased (observations=total interface length -r). For each fit we compute the goodness of the fit through R 2 .
As can be seen in Fig. 8, which shows R 2 vales for fits of the roughness B(r) at t = 0 hours averaged for the three fronts under control conditions, two regions emerge where the power law fits with B(r) ∼ r ζ are highly pertinent. At short lengthscales (for small r i ) we obtain very high R 2 values above 0.9998, even for relatively large window sizes w i . R 2 then abruptly decreases, and a second local maximum with values around 0.9992 is reached at ∼ 80µm. This behaviour suggests two separate power law scaling regimes are present in B(r).
To choose the best window size for each region we complemented the mathematical consideration of R 2 with additional biological arguments. For r i = 4µm, even if the goodness of the fit is high up to 50µm, we restricted the fit to not exceed the average length of a cell ∼ 20 µm, in order to capture the behavior of B(r) at subcell lengthscales, thus defining a first region I = [4,16]µm. For the second region, we considered that important mechanisms of intercell communication, in a wide range of mammalian epithelial cells, have a range of at most a few cell lengths -for example the cell networks reported for gap junction communication in rat kidney cells 3 as well as mechanical interactions via the cytoskeleton and force transfer via the substrate, which were found to decrease linearly with distance 4 . We note here that the inhomogeneous distribution of cell division we observed, which seemed to be more prevalent in cells close, but not directly at the colony edge, can also be considered as an    interaction at this range. We therefore restricted the second region to II = [80, 180]µm. For all the experimental cell fronts, fits of B(r) were performed over these two regions I and II to obtain ζ I and ζ II . As a further check to ensure that the second region of power law scaling was not a spurious artefact of crossover between a single power law scaling regime and a flat function, we also compared our experimental results to simulations of specifically such a situation, an elastic line evolving in a thermal bath from an ideal flat configuration to an equilibrium rough configuration in the Edwards-Wilkinson model (E-W).
The E-W elastic line model describes the evolution of an interface u(y,t) with elasticity c subject to a thermal bath of temperature T with friction η: η∂ t u(y,t) = c∂ 2 y u(y,t) + ξ (y,t).
For an initially flat elastic line the correlations in its geometry will evolve in time as a result of the competition between the interface elasticity and the thermal fluctuations 5 . For finite times, a memory of the initial conditions remains in the roughness B(r,t). As t grows correlations spread along the interface 6 .

B(r,t) =
Tr The roughness of the elastic line is then characterized by two well-defined regimes: a power law T c r 2ζ th at short lengthscales, and a flat function at large scales. In order to check if the crossover between these two regimes could lead to a spurious second power law regime, we analyze the roughness of an elastic line described by the arbitrary chosen values c = 1, T = 50, η = 1. In Fig. Fig. 9 we show the time evolution of the roughness at time intervals of 1000 (in arbitrary units). For comparison we also show the same data as presented in Fig. 1(c) of the main paper for the roughness of a cell colony under control conditions.
We analyze the roughness of the elastic E-W line using the same approach as described above for the cell fronts. The goodness of fit R 2 for the elastic line at t = 1000 is shown in Fig.10, with a clear maximum observed at short lengthscales, corresponding to the expected region of power law scaling. The goodness of fit then drops when approaching the crossover between the power law scaling at short lengthscales and the flattening at large ones. Finally, when the fitting window covers only large lengthscales at which the configuration of the E-W line is completely flat, R 2 rises again. We note that since there are no fluctuations in our analytical data, at large lengthscales R 2 takes infinite values (there is no difference between the real data and the predicted data). We therefore forced the R 2 values to one instead of infinity to emphasise that in this interval there is a perfect fit by a power law with exponent 0.
Finally, an alternate procedure is to extract the roughness exponent ζ in reciprocal space, from the power law scaling of the structure factor S(q) S(q,t) =ũ(q,t)ũ(−q,t), is the Fourier transform of the displacement field u(z,t) defining the interface position. Formally, the structure factor S(q) and the roughness correlation B(r) contain the same geometrical information and are related through While this approach is more prone to noise, where sufficient statistics are available, it provides generally more reliable estimates of ζ than the real space autocorrelation functions 7 , and can fully address aspects of anomalous scaling not accessible via B(r) 2,[8][9][10][11][12] In the structure factor, as shown in Fig. 11, we again observe two regions of power law scaling. Extracting the roughness exponent values from intervals corresponding to the real space regions I (subcell lengthscale, dark green) and II (few-cell lengthscales, red) we obtain ζ I = 0.56 ± 0.01 and ζ II = 0.29 ± 0.08, consistent with the values obtained from B(r).

Probability distribution function of relative displacements
As a check of the nature of the interface and the symmetry of its relative displacements, we carried out a multiscaling analysis 13,14 to evaluate the probability distribution function (PDF) of relative displacements ∆u(r, z) and its characteristic scaling properties reflected in the behavior of its central moments 15 : σ n (r) = |∆u(r)| n ∼ r nζ n , where ζ n are the associated scaling exponents for the nth moment. For 1-dimensional monoaffine interfaces, the PDF is Gaussian [16][17][18] , and the roughness function B(r) ≡ σ 2 (r) ∼ r 2ζ is sufficient to fully characterize the scaling, with a single-valued exponent ζ n = ζ ∀n.
As shown in Fig. 12, the PDFs obtained for different r values and over the time of evolution of the front are highly symmetric for r above 100 µm, although the PDF at lowest r = 10µm values shows slight asymmetry. Consistent with the evolution of B(r) discussed in the main paper, we see that for small lengthscales r, the PDF remains essentially unchanged with time. For larger lengthscales r = 10000, the initial PDF falls off from a Gaussian behaviour for small ∆u(r) to a steady background level, related to the effective flattening of the initial front configuration at high lengthscales. As the front evolves and roughens also at the higher lengthscales, the PDF for r = 10000 becomes more Gaussian.
The third (skewness) and fourth (kurtosis) moments show two regimes of power law scaling at subcell and few-cell lengthscales in regions I and II as defined in the main paper, from which the characteristic exponents ζ 3 and ζ 4 can be extracted. As shown in Figs. 13, 14, the value of the scaling exponents in region II appear to be universal, but diminish for region I.

Effective flattening by filling of cavities
During the measurements, we occasionally found large cavities at the edge of the colony, as shown in Fig. 15(a). The presence of these features significantly affected the dynamics of the colony, with cells rapidly moving towards the cavity and progressively filling it, as can be seen in Fig. 15(b,c), with the inflowing cells aligned radially towards the cavity center. Only when the cavity was filled did this section of the front continue to move outward again, as previously observed by Kim et. al. 19 . While we excluded affected segments of the fronts containing such features from the statistical analysis of roughness and dynamics, we believe that they are potentially related to the steady-state saturation of B(r) at higher lengthscales.

Numerical simulations
The numerical simulations of proliferating cell fronts were performed by using the epicell package, developed at the University of Geneva 20 , adapted to the modeling of mechanical properties of two-dimensional cell sheets 21,22 .
The epicell model used is based on solving the following dimensionless Hamiltonian on a regular grid of two-dimensional hexagonal cells: A 0 , H the Hamiltonian, A 0 the preferred cell area, K the cell area elasticity, Γ the cell perimeter contractility, Λ the inter-cell adhesion represented by the cell-cell edge tension, A α the area of cell α, L i, j the length of the inter-cell edge e i, j . The re-normalization of the Hamiltonian with the constraint that all cells are of the same fixed type, elasticity and preferred size, yields a two-variable system, with reduced parameters Γ and Λ. Additionally, the model was modified to include a non-trivial cell division descriptor, favoring divisions to occur just behind the cell front, in line with out empirical observations. This was implemented through setting a weight prefactor to the division probability following W div = 0 at the cell front (d = 0) and W div = exp − d−1 3 for cells with distance d behind the front.

13/17
Parameters and execution Seventeen parameters are necessary to run the freeBoundary application included in the epicell package. They are listed below with their description and values taken for the current study: • modelId: The model dictates which cells are selected for division in priority. Here, set to 3, so that the largest and less compressed cells would divide first.
• nRep: Number of replicates of the cell front, here set to 1.
• myseed: Seed for the random number generator, here set to 1.
• K: The cell area elasticity parameter K. Set to 2.256 · 10 9 N/m 3 following previous studies on mechanical properties of cell sheets.
• A0: Preferred cell area A 0 . Set to 3.14 · 10 −10 m 2 , following the experimental observations on the average size of the Rat1 cell line.
• ncols: Number of columns of cells within the cell sheet. Set to 20 here to have a thin but long cell sheet.
• nrows: Number of rows of cells within the cell sheet. Set to 10 3 here to have a thin but long cell sheet.
• xmin: Left hard boundary position. Set to 0 here so that the system imposes boundaries automatically after the initial relaxation step.
• xmax: Right hard boundary position. Set to 0 here so that the system imposes boundaries automatically after the initial relaxation step.
• ymin: Bottom hard boundary position. Set to 0 here so that the system imposes boundaries automatically after the initial relaxation step.
• nbDivMax: The maximum number of cell divisions before the simulation terminates. Set here to 10 5 to allow a large number of divisions compared to the number of cells in the sheet.
• elasticCoef: Elastic coefficient of the hard boundaries, used in this work to mimic an infinite in-compressible cell sheet on three of the sides as to promote growth along one direction. Set to 10 here to provide a hard boundary.
• sdLambda: Standard deviation of the inter-cell adhesion Λ. Set to 0 in the current study.
• sdGamma: Standard deviation of the cell perimeter contractility Γ. Set to 0 in the current study.
• boundaryLambda: Normalized line tension of the cell sheet edges, used for simulating an interaction of the cells with the external media. Set to 0 in the current study.
The freeBoundary code was executed on the University of Geneva cluster, Baobab, with a time constraint before termination of 12 hours. The 4000 Λ, Γ parameter pairs were computed in parallel and 881 simulations were successful. Unsuccessful simulations were either due to nonphysical parameters (1411 simulations), impossible convergence at the initial relaxation stage (1050 simulations), nonphysical mechanical behaviour of vertices (651 simulations), and finally unexpected simulation terminations due to specific pathological parameter sets (7 simulations). The outcome of all parameter pairs is shown in Fig. 16, with successful simulations labeled as OK.  Figure 16. Complete simulation space with 4000 Λ, Γ parameter pairs. Out of this, 881 simulations were successful, labeled as OK. The unsuccessful simulations either occurred due to nonphysical parameters (1411 with infinite forces and 651 with nonphysical vertices) or due to non-convergence during the initial relaxation step (1050 simulations). A minority, 7 simulations, ended unexpectedly due to specific pathological parameter sets.

Phase diagram investigation
Our study focused on the successfully simulated region of the phase diagram, as shown in Fig. 17. The phase diagrams for the asymptotic roughness exponent ζ , quality of fit R 2 for the asymptotic exponent, total front displacement from its initial position, velocity of the front over the final 100000 iterations, and the number of iterations for each simulation are shown in Fig.  17(a,c,e,g,i) respectively. As discussed in the main text, the analysis of simulation parameters enables the extraction of two branches of physical observables, with the lower branch exhibiting pathological behaviour including quasi-nonexistent front displacement/velocity, periodic features in the roughness, and large error during extraction of the roughness exponent ζ . Through the combination of the phase diagrams, a region of stability was defined visually as a red bounding polygon, with edges located at (Λ, Γ) = [(0, 0.03), (0, 0.06), (−1.975, 0.49), (−1.975, 0.27)]. This filtering enables the removal of pathological scenarios, and allows for a better selection of the appropriate models for comparison with in-vitro experiments. For instance, the overall average roughness exponent ζ = 0.69 as shown in the distribution Fig. 17(b) is highly asymmetric, whereas sampling ζ values in the region of stability yields an average roughness exponent ζ = 0.74 with a smaller error and a more symmetric distribution.
It is also noteworthy that whilst ζ and the number of iterations following the initial relaxation step (Fig. 17(i,j)) seem to be concentrated around specific values, the total displacement and final velocity in Fig. 17(e,f,g,h) are more evenly distributed, and in fact tend to increase with decreasing inter-cell adhesion and cell contractility.