Cell culture
For experiments, LN229 and U138 GB cells and primary GB lines (GBM4, GBM10) were used. LN229 cells were purchased from the American Type Culture Collection (ATCC, CRL-261, Manassas, VA, USA) and U138 cells were obtained from Cell Lines Service (Cell Lines Service, 300363, Eppelheim, Germany). LN229 cells were cultured using 89% (v/v) Roswell Park Memorial Institute medium (Lonza, Basel, Switzerland, BE12-115 F), supplemented with 10% (v/v) fetal bovine serum (FBS, Gibco, Carlsbad, CA, USA, 10500-064) and 1% (v/v) penicillin/streptomycin (P/S, Gibco, Carlsbad, CA, USA, 15140-122). U138, GBM4 and GBM10 were cultured in 89% (v/v) Dulbecco’s Modified Eagle Medium (Invitrogen, Waltham, MA, USA, 41965-062), and supplemented with 10% (v/v) FBS and 1% (v/v) P/S.
Primary GB were derived from human biopsies as reported earlier38. All patients provided written informed consent. The study was conducted in accordance with the Declaration of Helsinki and was approved by the local Ethics Committee of the Martin Luther University Halle-Wittenberg (project reference number: 2015 − 144).
For experiments involving mitomycin C treatment to inhibit proliferation, cells were subjected to a non-lethal dose of the substance, corresponding to 0.1 µg/ml and were continuously incubated and imaged with mitomycin C39.
Single cell migration
For measuring and analyzing single cell migration 1,000 cells were seeded per well in untreated 12-well plates (Greiner, Kremsmünster, Austria) 24 h prior to the start of experiments. Individual cells were imaged every 15 min with a microscope (Leica DMi8, Leica, Wetzlar, Germany) equipped with CO2 (5% (v/v)) and temperature (37 °C) regulation. The average cell speed and directionality were calculated based on cell tracking as described before using the Sobel operator for edge detection40. The directionality was defined as the ratio of the total distance travelled and the sum of incremental distances the cell moved between successive frames.
Collective cell migration
For collective cell migration experiments 450,000 cells were seeded per well in untreated 12-well plates. Twenty-four hours afterwards cells were imaged with an inverted microscope (DMi 8, Leica, Wetzlar, Germany), in a fully humidified atmosphere with 5% CO2 (v/v) and at 37 °C. Images were acquired every 3 min for 60 h. For analyzing cell migration particle image velocimetry (PIV), based on PIVlab, was used with a final cross-correlation window size of 16 × 16 pixels (pixel size: 0.48 μm), as described elsewhere, to obtain local velocity fields directly from phase contrast images12,31,41. Briefly, PIV is a pattern matching technique, which divides an image into patches and identifies the most similar image parts in the subsequent image using cross correlation. The displacement of each patch from the current to the following image is used to estimate the local velocity.
To assess the effective motion of cells, the mean squared displacement (MSD) \(\: < \left| {\Delta \:\overrightarrow {{x\left( {\Delta t} \right)}} } \right|^{2} >\) and its scaling coefficient α(Δt) were evaluated, using the following equation:
$$\: < \left| {\Delta \overrightarrow {{x\left( {\Delta t} \right)}} } \right|^{2} > = K\,\,{\text{* }}\Delta t^{{\alpha \:\left( {\Delta t} \right)}}$$
With the generalized diffusion coefficient K. The scaling coefficient retrieves information about the type of motion observed inside the monolayer, with sub-diffusive motion for α < 1, diffusive motion for α ≈ 1 and super-diffuse motion for α > 1.
Identifying proliferation events and cell density dependent migratory properties
To automatically identify individual proliferation events in the confluent layer, a machine learning model was used, as described previously31. Briefly, individual proliferation events in a dense monolayer were identified using a three-step process, consisting of an initial U-Net based segmentation, to reduce the region of interest and speed up evaluation, followed by candidate extraction using 2D cross-correlation. Finally, candidate divisions were classified as true or false positives using a GoogLeNet classifier. To match division events occurring over multiple frames the Munkres global nearest neighbor algorithm was used for detection-to-track-assignment. As cell divisions can be incomplete or cells can be excluded from the layer, having a similar appearance as mitotically rounded cells, but with a significantly longer visibility, detections that were present for more than 60 min were discarded.
To assess the influence of individual proliferation events on the layer velocity, for each identified proliferation event the speed around the division was calculated as a function of distance to the division and normalized to the speed of the rest of the layer. The same procedure was applied to up to 20 images, corresponding to 60 min, before and after the detection of the mitosis event. This period was chosen to assess the effect of pre-mitotic contraction and post-mitotic expansion. To further assess the cell density in each field of view, the initial cell number at t = 0 h was manually counted in a sub-field of 144 × 144 μm² and extrapolated to the whole field of view. To calculate the cell density ρ over the whole time, for each field of view the number of proliferations was added to the initial cell density. Combined with the layer migration speed v this approach allowed to calculate the scaling coefficient α between both quantities as follows:
$$\:v\:\sim \:\rho \:^{{\alpha \:}}$$
Obtaining the number of proliferation events and the initial cell density allows the calculation of doubling times for each cell line. Therefore, the current cell number N(t) at time t in a field of view and the cell number at Δt = 12 h later were used to estimate the instantaneous doubling time td as follows:
$$\:t_{d} = \frac{{\Delta \:t{\text{*}}\:N\left( t \right)}}{{N\left( {t + \Delta \:t} \right) – N\left( t \right)}}$$
Estimating the contribution of proliferation events to monolayer speed
The contribution of division events to the overall layer speed was estimated by the temporal and spatial contribution of pre-mitotic contraction and post-mitotic expansion obtained before. This allows to distinguish between all points in the image that are affected by proliferation events up to 1 h before and after the current time point and those that are currently not affected, classifying each data point as either “unaffected by proliferation events” or as “affected by proliferation events”. The “unaffected regions” were used to define the baseline speed, and the ratio of the speeds of both regions was considered as the relative contribution of proliferation events to the overall layer speed. Notably, the first and last hour of measurement, and the edges of images were discarded, as no data were available about proliferation events outside of the field of view or before and after the measurement.
Analyzing self-organization in monolayers
To determine local orientation of cells in a layer the largest eigenvector of the structure tensor J of the image I was used. The structure tensor J was defined as42:
$$J_{{pq}} \left( {\vec{x}} \right) = \:\int \: w(\vec{x} – \vec{x}\prime \:)\left( {\frac{{\partial \:I\left( {\vec{x}\prime \:} \right)}}{{\partial \:\vec{x}\prime \:_{p} }}\:\frac{{\partial \:I\left( {\vec{x}\prime \:} \right)}}{{\partial \:\vec{x}\prime \:_{q} }}} \right)\:d^{2} \vec{x}\prime \:\:\:\; with \; p,q\: \in \:\:[x,y]\:$$
and the Gaussian window function w. From the cellular orientation, the (local) nematic order parameter S was calculated as:
$$\:S = \: < 2*{\text{cos}}^{2} \theta \: – 1 >$$
With θ being the angle between the cellular orientation and the mean orientation of all cells, either in the whole field of view or a defined sub-window. To evaluate the spatial dependence of the order parameter, a power law was fitted as \(\:S\left(l\right)\sim\:{l}^{a}\), with the decay exponent a and the length scale l35. For long range order a approaches 0 and S remains constant. For 0 > a>-1 the order is considered to be of a quasi-long range type, with the extreme case a = -1 corresponding to a random order of orientations. For short-range order, an exponential decay is expected, but for the data given here, a power law describes the data best.
As nematics contain topological defects, ± 1/2 defects were identified using the Matlab toolbox “defector find”43 on the previously calculated cellular orientation fields. This toolbox calculates the winding number around potential topological defects in order to classify them as ± 1/2 defects. Notably, this allows assessing the evolution of the number of topological defects as a function of time and the movement of topological defects.
To evaluate if the ordering might be of a polar type, a polar order parameter V was introduced32:
$$\:\:V = \:\frac{1}{N}\:\sqrt {\left( {\sum {\:_{{i = 1}}^{N} } {\text{sin}}\vartheta \:_{i} } \right)^{2} + \:\left( {\sum {\:_{{i = 1}}^{N} } {\text{cos}}\vartheta \:_{i} } \right)^{2} }$$
Here, N corresponds to the number of angles ϑ in the whole field of view or a defined sub-window. ϑ is the angle between cellular orientation and velocity at the same spatial position. Similar to the nematic order parameter, the decay of the polar order parameter as a function of coarse graining was assessed. Here, an exponential decay function fitted the data best:
$$\:V\left( l \right) = m*e^{{ – \frac{l}{{l_{0} }}}} + V_{0}$$
With the characteristic decay length l0 and the constant offset V0.
Cell exclusion assay
For the cell exclusion assay 50,000 cells were seeded per well in a confined ring structure in 12-well plates. Twenty-four hours after cell seeding, the cells reached confluence. Thereafter, the ring structure was removed, and cells were transferred to the microscope (Leica DMi8, Leica, Wetzlar, Germany, 5% CO2 (v/v), 37 °C) for monitoring cell expansion. Five regions of interest were imaged per well. Cell expansion was monitored for 60 h with images taken every 5 min.
To analyze cell expansion a U-Net was trained to segment the cell layer from the background. Local velocities and other quantities were calculated as described under “Collective Cell Migration”, with the restriction of the analyzed area to the cell covered area.
For assessment of cell density dependent effects, another U-Net was trained to identify cell-cell boundaries in phase contrast images and thus assess cell sizes inside the layer. For training 52 manually segmented sample images, containing 15,015 cells, were used. For training and later analysis, the original images of size 1280 × 960 pixels were split into batches of 256 × 256 pixels. From each training image, 64 random patches of the batch size were chosen and randomly augmented as follows: random shearing in x- and y- direction from − 15° to 15°, random rescaling by factors of 0.7 up to 1.3, and random rotations from 0 to 360°. The U-Net was trained for 50 epochs, with a learning rate drop factor of 0.8 and an initial learning rate of 0.5 *10− 4.
The confusion matrix, accuracy and F1-score for the pixel wise segmentation of the test data set were calculated to test the trained U-Net. In addition, an estimate of the goodness of cellular identification was performed. Therefore, the raw segmentation was post-processed as follows: disconnected objects smaller than 50 pixels were removed, followed by morphological closing and skeletonizing of the binary image. For skeletonizing a minimal branch length of 150 pixels was chosen. For matching, the centers of the segmented cells were compared to the cell centers obtained from the ground-truth measurements. If the distance between the centers of an automatically and manually segmented cell was less than 4.8 μm (5 pixels), the segmented object was declared as true positive. True positives were considered as correctly identified cells, false positives as segmented cells with no sufficiently close (4.8 μm) cell in the ground truth image and false negatives as cells in the ground truth image with no sufficiently close (4.8 μm) cell in the segmented image.
To further improve the detection of cells, too small (less than 120 μm²) and too large (more than mean size + 2 * median absolute deviation of size) segmented objects were discarded. Furthermore, segmented objects have been assigned to tracks via the Munkres global nearest neighbor algorithm and segmented objects that have not been assigned in at least 4 of the 5 following images were discarded too, to reduce the number of misdetections. The following energy function was used as input for the Munkres algorithm for the assignment of detection d to track t:
$$\:Cost\:\left( {d,t} \right) = \left( {1 + n} \right)*\Delta \:x\left( {d,t} \right)*(1 + |\Delta \:A\left( {d,t} \right)|)$$
Here, n corresponds to the number of frames the detection d was not assigned to any track t. Δx(d, t) was the normalized distance of the detection d to the last detection of track t. ΔA(d, t) corresponds to the normalized difference in cell size between the currently detected cell d and the last detection of track t. If the assignment cost for a detection was too high, a new track was created.
To calculate the gradient of the cell density in a field of view, the cell density in dependence of the distance to the cell cluster boundary was calculated and linearly fitted. The slope of the fit was used as cell density gradient estimate.
To assess the relation between the expansion rate of the monolayer and the cell density gradients, piece-wise linear fitting was employed. To identify segments, the method introduced by R. S. Bogarts was used44.
Statistics
Statistics were performed using the two-tailed ANOVA with Tukey post-hoc test or the two-sided sign test. Significance was defined for p < 0.05. All error-bars and shaded areas depict the standard error of the mean. Experiments were repeated at least three independent times. For live-cell imaging experiments assessing collective motion, 4–8 fields of views were measured per experiment and condition. Confidence intervals for fit parameters were obtained using the MatLab function confint.
