Disentangling collective motion and local rearrangements in 2D and 3D cell assemblies

Roberto Cerbino *a, Stefano Villa a, Andrea Palamidessi b, Emanuela Frittoli b, Giorgio Scita bc and Fabio Giavazzi *a
aUniversità degli Studi di Milano, Dipartimento di Biotecnologie Mediche e Medicina Traslazionale, 20090 Segrate, Italy. E-mail: cerbino@gmail.com; fabio.giavazzi@unimi.it
bIFOM-FIRC Institute of Molecular Oncology, Milano, Italy
cUniversità degli Studi di Milano, Dipartimento di Oncologia e Emato-Oncologia, Milano, Italy

Received 14th October 2020 , Accepted 3rd December 2020

First published on 8th December 2020

The accurate quantification of cellular motility and of the structural changes occurring in multicellular aggregates is critical in describing and understanding key biological processes, such as wound repair, embryogenesis and cancer invasion. Current methods based on cell tracking or velocimetry either suffer from limited spatial resolution or are challenging and time-consuming, especially for three-dimensional (3D) cell assemblies. Here we propose a conceptually simple, robust and tracking-free approach for the quantification of the dynamical activity of cells via a two-step procedure. We first characterise the global features of the collective cell migration by registering the temporal stack of the acquired images. As a second step, a map of the local cell motility is obtained by performing a mean squared amplitude analysis of the intensity fluctuations occurring when two registered image frames acquired at different times are subtracted. We successfully apply our approach to cell monolayers undergoing a jamming transition, as well as to monolayers and 3D aggregates that exhibit a cooperative unjamming-via-flocking transition. Our approach is capable of disentangling very efficiently and of assessing accurately the global and local contributions to cell motility.

1 Introduction

In the last decade, the characterization of multicellular motility, often in combination with the simultaneous, space-resolved measurement of other biophysical and biochemical parameters, such as traction forces,1 internal mechanical stresses,2 local rheological properties,3 or levels of biochemical signalling4 has played a pivotal role in improving the biophysical understanding of fundamental biological processes. Notable examples of these processes are different types of -taxes,2,5,6 embryonic morphogenesis,7,8 jamming/unjamming transitions,9,10 and cancer invasion.11,12 One of the most popular ways to extract the instantaneous displacement/velocity field from time-lapse microscopy experiments in multicellular assemblies is Particle Image Velocimetry (PIV). Originally introduced to measure flow fields in fluids,13 PIV requires seeding the fluid with small tracer particles that are advected by the flow, providing at the same time a sufficiently high optical contrast. The ideal operating condition for PIV is the one where the tracers appear as “point-like”, well resolved, highly contrasted objects, with a number density that is large enough to ensure a good spatial resolution. If these conditions are met, even approximately, and a sequence of images of a fixed region of interest (ROI) is acquired with a temporal separation Δt between consecutive frames, the local displacement of the tracers between two consecutive frames can be estimated by locating the center of mass of the peak of the cross-correlation function of the intensity distribution of the two frames. If this operation is performed on interrogation windows that are smaller than the ROI, a space-resolved map of the tracer displacement between the two frames is obtained, which becomes a velocity map after division for Δt.

The introduction of PIV to analyze cell biology imaging experiments is relatively recent.14 Its application for characterizing the motility of epithelial monolayers was pioneered about ten years ago in a seminal work of Petitjean et al.15 and it is now a widely adopted approach. In most of the cases where PIV has been used to monitor multicellular motion, this is done without using added tracers, which represents an important difference with PIV in fluids. Instead, inherent structures of the cells (the nucleus, organelles, components of the cytoskeleton, the plasma membrane…) play the role of “tracers” for the PIV algorithm, as long as they are capable of producing a textured, non-uniform intensity distribution in the image. Despite its popularity, the use of PIV with cells suffers from a number of intrinsic limitations. In the first place, the accuracy in the determination of the cross correlation peak is strongly affected by the finite size of the tracers and by their continuous shape changes and fluctuations. Moreover, nearby “tracers” inside the cell can have a significant relative motion, which is not necessarily related to the global migration of the collective. This makes the results of the PIV analysis dependent on the choice of the size of the interrogation window. Often generous smoothing over space and time and careful filtering of the estimated velocity field are required to obtain meaningful results and cogent velocity maps.16

Other velocimetry methods, such as Optical Flow Velocimetry, are less affected by some of the aforementioned problems, and have therefore become increasingly popular alternatives to PIV for the quantification of cellular motility.17 Nevertheless, the main limitation, shared by all velocimetry methods, is that they implicitly rely on the assumption of spatio-temporal “smoothness” of the probed displacement field. This is valid if a finite and experimentally accessible spatial correlation length in the particles displacements ensures meaningful spatial averages, and the temporal persistence of directed motion is larger than the sampling interval Δt.

A paradigmatic, non-biological example where velocimetry methods are doomed to fail is a suspension of Brownian particles. In this case, there are negligible correlations in the displacements of distinct particles and there is no temporal persistence in the particles trajectories, at least on timescales accessible with microscopy. As a consequence, the instantaneous velocity in an arbitrarily sized interrogation window is intrinsically ill-defined. This limiting case demonstrates that, in the presence of a random and poorly correlated motion, PIV does not give reliable results. While this may seem a case that is rarely applicable to cells, it has been recently observed that during the late stages of maturation of an epithelial monolayer the velocity correlation length decreases dramatically as cell collective jam.18 This can make difficult obtaining reliable velocimetry data. Single particle tracking (PT), which is based on tracking over time the position of the individual tracers, does not suffer from these limitations but represents a suitable option only in the presence of high quality, specific, fluorescently labelled, sub-cellular entities, typically the nucleus. Notwithstanding the impressive improvement in performance brought in the last few years by the implementation of machine learning-based approaches,19,20 PT in biological samples remains technically challenging with fluorescence images and even more so is with bright-field or phase-contrast images. Thus, very often, instead of PT, pseudo-trajectories by integration of PIV maps are generated9 and used to infer single-cell information, such as the mean square displacement or the persistence length.

During the last years, we have been developing tools to obtain single-particle information with tracking-free algorithms such as Differential Dynamic Microscopy (DDM)21–25 or, more generally, Digital Fourier Microscopy.26 In particular, DDM is based on the simple idea that the motion of tracer particles can be captured by calculating the difference between two image frames acquired at different times. Intuitively, the contrast of the signal in the different images increases if the temporal separation between the two considered frames is increased. This contrast can be quantify by calculating the variance of the image difference.21 This property was recently exploited to quantify the dynamics of dense, glassy two-dimensional (2D) colloidal systems by means of a so-called Differential Variance Analysis (DVA).27

Here, for the first time, we apply DVA to characterize the motility of cells belonging to two-dimensional (2D) and three-dimensional (3D) collectives. Our procedure decouples the global dynamics, obtained via image registration, from the local one, which is obtained with standard DVA, as well as with the here-introduced sr-DVA, to obtain space resolved activity maps, and 3D-DVA, to quantify dynamics in more complex 3D settings. The image registration step enables the robust determination of the global features of cellular motion, in terms of mean linear/angular velocity, tissue stretching, and so on. It also removes effectively all affine contributions and highlights local contributions to the dynamics, which are mainly due to mutual cell displacements and structural rearrangements, for the subsequent DVA. To demonstrate potential and versatility of our approach, we consider three paradigmatic model systems: (i) a confluent monolayer of epithelial cells initially seeded at relatively low density, evolving toward a kinetically arrested, jammed state,10,18,28 (ii) a fully jammed monolayer where the overexpression of RAB5A, a small GTPase that regulate endosomal biology, induces coordinated flocking-like motion and unjamming,10,29,30 (iii) a 3D spheroid of tightly packed epithelial cells from a clonal breast cancer line embedded in a thick collagen matrix. In this system, the expression of RAB5A is sufficient to induce persistent rotation and local fluidisation.12

2 Materials and methods

2.1 Cell culture

MCF-10A cells were a kind gift of J. S. Brugge (Department of Cell Biology, Harvard Medical School, Boston, USA) and were maintained in Dulbecco's Modified Eagle Medium: Nutrient Mixture F-12 (DMEM/F12) medium (Invitrogen) supplemented with 5% horse serum, 0.5 mg ml−1 hydrocortisone, 100 ng ml−1 cholera toxin, 10 μg ml−1 insulin and 20 ng ml−1 EGF. The cell line was authenticated by cell fingerprinting and tested for mycoplasma contamination. Cells were grown at 37 °C in humidified atmosphere with 5% CO2. MCF-10A cells were infected with pSLIK-neo-EV (empty vector control) or pSLIK-neo-RAB5A lentiviruses and selected with the appropriate antibiotic to obtain inducible, stable cell lines. Constitutive expression of EGFP-H2B was achieved by retroviral infection of MCF-10A cells with pBABE-puro-EGFP-H2B vector.

MCF10.DCIS.com obtained from Dr John F Marshall (Barts Cancer Institute, London, UK) were infected with pSLIK-neo-EV (empty vector, CTR) or pSLIK-neo-RAB5A lentiviruses and selected with the appropriate antibiotic to obtain stable inducible cell lines. Constitutive expression of EGFP-LifeAct, mCherry or EGFP-H2B was achieved by lentiviral and retroviral infection of MCF10A and MCF10DCIS.com cells with EGFP-LifeAct-puro or pBABE-puro-mCherry-H2B/pBABE-puro-EGFP-H2B vectors, respectively.

Transfections were performed using either calcium phosphate or FuGENE HD Transfection Reagent (cat. no. E2311, Promega), according to the manufacturers instructions. Lentiviral infection was performed as described previously.10

2.2 Microscopy observations

MCF10A monolayer at low confluency. Cells were seeded in six-well plates at 0.125 × 106 cells per well in complete medium and cultured for 24 hours to allow them to attach and reach subconfluency. At the time of recording, fresh media containing 20 ng ml−1 of EGF was added. Picture were taken every 10 min for a total 82 hours using an Olympus ScanR inverted microscope with 10× objective. Five fields of view (FOVs) were taken for each time point in the center of the six well and triplicate experiments were commonly run.
RAB5A and CTRL MCF10A confluent monolayer. Cells were seeded in six-well plates (1.5 × 106 cells per well) in complete medium and cultured until a uniform monolayer had formed. RAB5A expression was induced, where indicated, 16 h before performing the experiment by adding fresh complete media supplemented with 2.5 μg ml−1 doxycycline to cells. At the time of recording, fresh media containing EGF was added. An Olympus ScanR inverted microscope with 10× objective was used to acquire images with a frame rate of 10 min over a 24 h period. For each condition (RAB5A and control), five independent FOVs, much smaller than the entire culture plates and far from the boundaries, are captured. Each FOV is imaged both in fluorescence and in phase-contrast perfusion. After cell induction, doxycycline was maintained in the media for the total duration of the time-lapse experiment.
3D spheroids kinematic assay. MCF10DCIS.com cells were plated on ultra-low-attachment-surface six-well plates (cat. no. 3471 Corning) at a density of 5 × 103 cells per well. Cells were grown in serum-free conditions for 10 days by adding fresh culture media every 2 days. Then spheres from every single well were collected and resuspended in 150 μl of 6 mg ml−1 collagen type I (cat. no. 35429 Corning), and diluted in culture media, 50 mM Hepes, 0.12 NaHCO3 and 5 mM NaOH. The unpolymerized mix of spheres/collagen was placed in eight-well chamber slides and incubated at 37 overnight. The following day, before imaging, 2.5 μg ml−1 doxycycline hyclate was added to the polymerized collagen mix to induce RAB5A expression. Time-lapse imaging of the motility of spheroids was performed on a Leica TCS SP8 laser confocal scanner mounted on a Leica DMi8 microscope equipped with motorized stage; a HC PL FLUOTAR 20×/0.5 NA dry objective was used. A white-light laser was used as the illumination source. Leica Application Suite X was the software used for all the acquisitions. In all cases, the assays were performed using an environmental microscope incubator set to 37 °C and 5% CO2.

2.3 Image processing

The processing of a typical time-lapse movie entails two main steps, as schematically shown in Fig. 1. For each pair of consecutive frames (panel a), we first identify the affine transformation that better accounts for the changes occurred in the images. This step enables the time-resolved determination of the global features of the multicellular motion, in terms of directed migration (net displacement of the tissue in a given direction), rigid rotation, contraction or expansion of the aggregate along specific axes (panel d). In the second step (panel c), we quantify the non-affine part of the cellular motion, i.e. changes in the mutual positions of the cells. To this aim, we consider the movies registered in the first step and we adapt the standard DVA in ref. 27 to calculate the overlap parameter Q(τ) for our experiments. Q(τ) captures the progressive loss of correlation that occurs as the cells displace from their original mutual positions.31 In particular, the characteristic timescale τ* associated with the decay of Q(τ) from 1 to 0 corresponds to the time scale of the structural relaxation. Whether this decorrelation process occurs “smoothly” in space and time or, on the contrary, is characterized by spatial heterogeneity and temporal intermittency, can be assessed by considering the so-called dynamical susceptibility χ4(τ), that represents the mean squared amplitude of the fluctuations of Q(τ) around its mean value.32 Large values of χ4(τ), in correspondence of the structural relaxation time τ* are a hallmark of cooperative rearrangements. The peak value of χ4(τ) provides an estimate of the characteristic size (in terms of the number of cells) of a collectively rearranging region, provided that χ4(τ) is suitably normalized.10 To perform the entire image processing, we have developed custom MATLAB code. In particular, the function imregtform is used for identifying the linear transformation connecting two consecutive frames, while the warping operation involved in the calculation of the registered image I([T with combining macron](t,τ)x,t) (see eqn (5)) is performed using the function imwarp.
image file: d0sm01837f-f1.tif
Fig. 1 Schematic representation of the main steps of the image processing pipeline. (a) Snapshots of the same cell monolayer acquired at different times t0, t0 + Δt,…,t0 + τ1 − Δt,t0 + τ1; two subsequent frames I(x,t) and I(x,t + Δt), are connected by a matrix T(t) (eqn (1)) that maximizes the overlap between I(x,t) and I(T(t)x,t + Δt). (b) For an arbitrary time delay τ1 between two images, the registration step consists in identifying the affine transformation T(t0,τ1) connecting them, which is obtained as the composition of the affine transformations connecting pairs of consecutive images between them (eqn (4)). (c) Once the optimal registration is found, non-affine displacements of cells between two images separated by τ1 are quantified as q(τ1,t0) = |I(x,t0) − I(T(t)x,t0 + τ)|2. (d) An estimate of the instantaneous kinematic parameters (e.g. the migration velocity v and the angular velocity ω) is obtained from the parameters of the transformation found in (b). (e) After suitable normalization (see main text) the average and the variance of q(τ1,t0) over time provide estimates for the overlap parameter Q(τ1) and the dynamical susceptibility χ4(τ1), respectively.
Step I: image registration. The affine registration of a sequence of 2D or 3D images is a standard image processing operation, for which several algorithms are available.33 Here we consider a simple, intensity-based implementation, that proved to be robust and computationally effective for all the cases analyzed in the present work. For the sake of clarity and simplicity, the following description will be provided for the case of a time sequence of 2D images, even though we have used a similar algorithm also to register temporal sequences of 3D images, or volume stacks.

Let us consider a time sequence of 2D images I(x,t), where x = [x,y]. Let T(a) be an affine transformation of the 2D plane in itself, dependent upon a set of parameters a0. T can be represented by an extended 3 × 3 matrix

image file: d0sm01837f-t1.tif(1)
where 0 = [0,0], u = [ux,uy] is a translation vector and M is 2 × 2 matrix. Exploiting the so-called polar decomposition,34M can be written as the composition of a biaxial stretch (by factors λ1 and λ2, along two orthogonal directions forming and angle ϕ with the coordinate axes) and a rotation by an angle θ:
image file: d0sm01837f-t2.tif(2)
where R(θ) is a clockwise rotation by an angle θ
image file: d0sm01837f-t3.tif
Accordingly, T is univocally determined by the set a = {u,λ1,λ2,ϕ,θ}. A similar parametrization holds also in 3D. For each time point t we determine the set of parameters ā(t) = {ū(t),[small lambda, Greek, macron]1(t),[small lambda, Greek, macron]2(t),[small phi, Greek, macron](t),[small theta, Greek, macron](t)} that better reproduces the changes occurred in the image between time t and t + Δt, the criterion being the minimization of the squared difference
image file: d0sm01837f-t4.tif(3)
where Δt is the time delay between consecutive images.

In particular, image file: d0sm01837f-t5.tif provides an estimate for the velocity associated to directed migration, image file: d0sm01837f-t6.tif is the instantaneous angular velocity characterizing the in-plane rotation of the tissue, while image file: d0sm01837f-t7.tif (ν = x,y) can be interpreted as the strain rates associated with global deformations.

Given two non-consecutive frames, taken at times t and t + τ, the affine transformation connecting them can be obtained as a composition of the affine transformations connecting pairs of consecutive images between them

[T with combining macron](t,τ) = T[ā(t + τ)]○T[ā(t + τ − Δt)]○…○T[ā(t + Δt)]○T[ā(t)](4)
where the symbol ○ indicates a composition of operators that becomes a matrix product for the corresponding matrices.

We note that, while in discussing image registration, we consider the case of a general affine transformation, the class of transformations to be used (and so the span of the parameters space) can be suitably restricted to better fit specific experimental geometries. For example, in the case of a confluent tissue confined in a circular patch (see ref. 35 for example), one could consider only rotations around the patch center, while in a thin channel (as in ref. 36), the main features of the collective migration could be captured by considering only a translation and a stretch along the channel axis.

Step II: difference variance analysis. In the absence of mutual cell displacements, the image I(x,t + τ) and the registered image I([T with combining macron](t,τ)x,t) are expected to overlap perfectly, as far as detection noise and interpolations errors in the image registration process can be neglected. This observation suggests that the squared difference
d(τ|x,t) = [I(x,t + τ) − I([T with combining macron](t,τx,t)]2(5)
is thus a good proxy of the heterogeneous, non-affine component of the cellular motion. As schematically shown in Fig. 1c, eqn (5) is a differential map defined over a suitable ROI (consisting of, say, M pixels) included in the intersection between the region covered by the original images and its registered replica under the transformation [T with combining macron](τ,t) for all the considered values of τ and t.

According to ref. 27, we define the instantaneous overlap parameter as

image file: d0sm01837f-t8.tif(6)
where the symbol 〈fx indicates a spatial average image file: d0sm01837f-t9.tif. The quantity d can be either estimated as the plateau of 〈d(τ|x,t)〉 for large τ (as suggested in ref. 27) or using the formula d = 2[ 〈I2(x,t)〉x,t − 〈I(x,t)〉x,t2], where 〈ft indicates a time average. This formula holds for ergodic samples also when the system did not relax completely.

According to eqn (6), the presence of a significant amount of delta-correlated noise in the images can introduce an systematic offset value in the overlap parameter: image file: d0sm01837f-t10.tif, where A is mean squared amplitude of the noise. If a reliable estimate of the noise amplitude A can be obtained (see for example ref. 37), the noise contribution can effectively removed via the substitution dd − 2A. In any case, we observe that delta-correlated noise only affects the amplitude of overlap parameter and becomes irrelevant if a dynamical, τ-dependent information is extracted with a suitable fitting procedure.

The time-averaged overlap parameter Q(τ) is defined as

Q(τ) = 〈q(τ,t)〉t.(7)

A characteristic relaxation time τ* of the overlap parameter can be extracted via a fit (e.g. with an exponential or a stretched exponential function), or considering its half-life. As discussed in detail in ref. 38, the shape of Q(τ) is determined by a combination of both dynamic and static properties of the sample. Let us assume that the intensity distribution in each image can be modeled as a collection of moving “objects” or “features” with a typical linear size l, which could estimated for example as the width of the peak of the spatial self-correlation function the images CI(r) = 〈I(x0 + r,t)I(x0,t)〉x0,t. In this scenario, we can identify τ* as the characteristic time required for these objects to displace from their initial positions by an amount of the order of l. This implies that, in general, the results of a DVA on a given sample do depend on the specific imaging modality, as this influences which features contributes the most to the optical signal and determine how they appear (e.g. in terms of contrast and effective size) in the collected images. The dynamic susceptibility χ4(τ), that captures the intrinsic fluctuations associated with the overlap parameter, reads

χ4(τ) = Nq2(τ|t)〉t − 〈q(τ|t)〉t2.(8)
The normalization factor N, which coincides with the average number of cells within the considered ROI, enables interpreting the amplitude of χ4(τ) as the number of cells in a collectively rearranging region. In this work, N was determined by operator-supervised automatic cell counting.

Three-dimensional DVA (3D-DVA). With a few adaptations, the entire image processing procedure described in previous paragraphs can be readily extended to sequences of 3D images, or volume stacks, obtained, for example, from time-lapse, multi-z, confocal microscopy acquisitions. One obvious difference with the 2D case is the number of parameters needed to identify the spatial transformation involved in the registration step: for a generic 3D affine transformation the number of independent parameters is twelve, whereas for a rigid transformation (roto-translation) this number decreases to six.

The numerical implementation of the procedure should include a proper scaling step accounting for the fact that, usually, 3D stacks have non-cubic voxels, as a consequence of the axial optical resolution being typically worse than the lateral one. This imposes also an intrinsic limitation to the quality of any image registration procedure, in particular if rotations of large amplitude about a rotation axis lying on the xy plane are involved.33 Moreover, in 3D fluorescence imaging, even in the presence of a uniform spatial distribution of fluorescent molecules, a systematic variation of the image intensity along the optical axis is often present (see e.g.Fig. 5b). This is due to the optical absorbance of the sample that makes the intensity collected on a given plane increasingly dimmer as one goes deeper in the sample.

To minimize the impact of the above-mentioned effects on the analysis results one can restrict the calculation of the overlap parameter to suitably short delay times τ, such that the global transformation to be compensated is not too far from identity and remapping artifacts are less severe. Moreover, the 3D ROI over which spatial averages are performed has to be chosen well within the region of optimal illumination/collection.

Space-resolved DVA (sr-DVA). The above definitions implicitly assume the space-invariance of the system under consideration. While this assumption holds in some cases (for example for a confluent monolayer observed far from the boundaries), in general it can be worth considering a time-averaged but space-resolved indicator. To this end, in analogy with eqn (6), we define a local time-averaged overlap parameter as
image file: d0sm01837f-t11.tif(9)
where the normalization factor d(x) = 2[〈I2(x,t)〉t − 〈I(x,t)〉t2] compensates for potential optical inhomogeneity within the sample, due for example to uneven illumination or non-uniform concentration/expression of a fluorescent tag. A simple and effective way to extract a space-resolved activity map from q(τ|x) is to consider a fixed delay time τ0, shorter than the characteristic relaxation time, and look at the quantity q(τ0|x).

We note that, if the sample is non-ergodic or if its dynamics is so slow that a complete relaxation cannot be observed during the experimental window, the above definition of d(x) can lead to a systematic underestimate of q.

Particle tracking analysis. In order to obtain, at least for the case of 2D monolayers with fluorescently labelled nuclei, an independent estimate of the kinematic parameters, of the overlap parameter and of the corresponding susceptibility, we also perform PT on the image sequences analysed with DVA. Nuclei are identified frame by frame by finding local maxima in the Laplacian transform of the image and subsequently segmented through a watershed segmentation process. Nuclear trajectories xn(t) are then obtained by linking the positions of the same nucleus in subsequent frames. Our custom MATLAB algorithm is able to correctly capture complete trajectories for a fraction of the nuclei present in the time lapse (see ESI for details). An estimate of the instantaneous velocity of each nucleus is obtained as vn(t) = [xn(t + Δt) − xn(t)]/Δt, while the mean migration velocity (i.e. the velocity of the center-of-mass of the monolayer) is calculated as the average over all tracked nuclei image file: d0sm01837f-t12.tif. Let image file: d0sm01837f-t13.tif be the position on the n-th nucleus at time t in the center-of-mass reference frame. In order to calculate the overlap parameter, we define an “overlap function”31 as image file: d0sm01837f-t14.tif, where a0 is a parameter controlling the overlap length-scale. Qualitatively, wn(τ|t) ≃ 1 if the displacement of the n-the nucleus (evaluated in the monolayer center-of-mass reference frame) during the interval [t,t + τ] is much smaller than a0, while wn(τ|t) ≃ 0 when the displacement largely exceeds a0. The time-dependent overlap parameter qPT(t,τ), is then obtained by averaging the overlap function over all the N tracked nuclei: image file: d0sm01837f-t15.tif, while the time-averaged overlap parameter are defined in close analogy with eqn (7) and (8): QPT(τ) = 〈qPT(τ|t)〉t and χ4,PT(τ) = NqPT2(τ|t)〉t − 〈qPT(τ|t)〉t2. We note that according the above definitions, QPT(τ) can be interpreted as the average fraction of nuclei that, during a time interval of duration τ, displace from their initial positions by less than a0.

3 Results and discussion

In this section, we present the results obtained by applying the image processing pipeline described in Methods to a number of in vitro model systems: (i) a confluent monolayer of mammary epithelial cells undergoing a jamming transition,10 (ii) a mature, kinetically arrested monolayer, induced to acquire a large-scale coordinated migration and to unjam by the expression of RAB5A,10,30 (iii) a 3D spheroid of tightly packed epithelial cells from a clonal breast cancer line embedded in a thick collagen matrix that display persistent rotation and local fluidisation upon RAB5A expression.12

3.1 Cellular jamming in 2D

In Fig. 2a–c, we show some representative frames extracted from a 88 h long time-lapse video acquisition of a monolayer of MCF-10A epithelial cells initially seeded at a density of 0.125 × 106 cells per well (see also Movies SM01, ESI). At the beginning of the experiment (Fig. 2a), the cells display a fast and highly coordinated directed migration, characterized by a migration speed v0 ≃ 7 μm h−1. After an initial increase of v0 to values of the order of 9 μm h−1, the migration speed rapidly decreases after about 15 hours from the beginning of the experiment to reach a value of ≃0.2 μm h−1 after about 50 hours. The migration speed profile shown in Fig. 2d is the main output of the image registration step, which was performed by imposing a rigid translation.
image file: d0sm01837f-f2.tif
Fig. 2 Jamming transition in a confluent monolayer of MCF10A epithelial cells. (a–c) Representative phase-contrast images (taken 0, 32 and 64 hours from the beginning of the experiment, respectively) of the confluent monolayer. The size of each box corresponds to 0.5 mm. (d) Step I of analysis gives an estimate of the monolayer collective migration velocity v0, here shown for one representative field of view (FOV). (e) Step II provides the overlap parameter Q(τ) of the monolayer for different ages (16.7 h, 25 h, 33.3 h, 41.7 h, 50 h, 58.3 h, 66.7 h, 75 h, 83.3 h). Each curve is obtained by considering a time interval of 1000 minutes and averaging over 5 independent FOV. (f) Age-dependent relaxation rate γ* (blue circles) associated with the decay of Q(τ); rate γ4 (orange squares), obtained as the inverse of the time delay for which χ4(τ) has a peak. Inset: The height χ4* of the peak of χ4(τ) is plotted as a function of the relaxation rate γ*. (g) Dynamical susceptibility χ4(τ) obtained in the same conditions as (e).

Once the migration velocity is determined, the registered stack (see Movies SM01_reg, ESI) is analyzed with DVA to quantify the residual fluid-like motility of the cells. In order to obtain an age-resolved determination of the dynamics, the entire image sequence is divided into ten evenly spaced, partially overlapping sub-sequences. Each sub-sequence covers a time interval of duration 1000 min and is analysed separately. We observe a progressive slowing-down of the dynamics over time, which can be clearly appreciated from inspection of the overlap parameter Q(τ) curves obtained for the different monolayer ages (Fig. 2e). The characteristic decay time of Q(τ) can be estimated for instance by means of its half-life τ*, whose inverse, the relaxation rate γ*, displays a temporal evolution (Fig. 2f, blue circles) that appears well distinct from the one characterizing the global migration speed v0 reported in Fig. 2d.

The dynamical susceptibility χ4(τ) exhibits (Fig. 2g) a peak in correspondence of a characteristic time τ4. Such peak shifts toward larger lag times as the monolayer ages and is in very good agreement (Fig. 2f), orange squares with the relaxation time τ* of Q(τ). The height χ4* of the peak of χ4 also displays an overall decrease over time, suggesting that the size (in terms of number of cells) of the cooperatively rearranging regions decreases progressively during the monolayer maturation, reaching values that are indicative of a single cell in the latest stages of the process.

This finding, namely the evidence of a decreasing dynamical length scale accompanying the kinetic arrest is in stark contrast with the expected behaviour based on a naive analogy between a maturating monolayer and a passive systems undergoing the glass transition.39 It is worth mentioning that a similar phenomenology, related to a progressive “loss of organization” of the motility pattern in jamming monolayers, was previously observed in ref. 18 and, in the very latest stages of jamming, also in ref. 28. Interestingly, in ref. 18 this feature was found to be associated not so much to the gradual increase in cell density, but rather to the maturation of cell–cell and cell substrate contacts.

3.2 Cellular unjamming in 2D

The jammed state of a monolayer, in which mutual cell rearrangements are prevented, is thought to be representative of the physiological condition of a real tissue. Recently, however, it has been shown that densely packed, kinetically arrested jammed monolayers can reawaken their motility and unjam in response to a variety of external stimuli, notable example being mechanical compressive stresses,9 exposure to ionizing radiation,40 serum deprivation and re-stimulation,41 over-expression of the endocytic protein RAB5A.10,12 Of note, experiments performed with MCF-10A cells showed that RAB5A over-expression does not significantly alter the migration of individual cells, which behave identically to control cells in terms of mean velocity, directional persistence and extension of protrusions.10

To gain some insight into the unique cooperative features of the RAB5A-induced unjamming we applied our imaging analysis approach to time-lapse videos of RAB5A-expressing and control cells. We analyzed movies of H2B-GFP-expressing MCF10A cells, which display brightly fluorescent nuclei as the main source of contrast in the images. Movies were also simultaneously acquired in phase-contrast mode to determine whether the results are influenced by the modes of imaging. (see Movies SM02_CTRL, SM02_RAB5, SM02_CTRL_reg and SM02_RAB5_reg, ESI) The temporal profile of the collective migration speed obtained from the registration of the fluorescent movies by imposing a rigid translation (Fig. 3a) is in excellent agreement with the one found by analyzing the phase-contrast movie (Fig. 3d), both for the RAB5A cells (red line) and control cells (blue line). Conversely, the overlap parameter Q(τ) is more sensitive to the imaging mode (Fig. 3b and e) with relaxation times that are systematically shorter (by approximately a factor of 2) for the phase-contrast movies than for their counterparts obtained with fluorescence images. According to the discussion in Section 2.3, we ascribe this discrepancy to the different characteristic size l of the features contributing to the image intensity in the two imaging conditions. This hypothesis is supported by the spatial self-correlation functions of the image intensity obtained with fluorescence and phase-contrast microscopy (Fig. S3, ESI). Taking the half width at half maximum (HWHM) of the self-correlation functions (Fig. S4, ESI), we obtain lFLUO = 4.8 ± 0.2 μm for fluorescence images, while for phase contrast images, where a large number of small-scale features are present, we obtain and lPC = 2.0 ± 0.2 μm.

image file: d0sm01837f-f3.tif
Fig. 3 Unjamming transition in a confluent monolayer of MCF10A epithelial cells. Collective migration velocity (a), overlap parameters (b) and dynamic susceptibility (c) for RAB5A-overexpressing cells (red curves and symbols) and control cells (blue curves and symbols) obtained by analyzing fluorescence images of H2B-GFP cells. Results from PT analysis of the same fluorescence image sequences are also shown in panels (a), (b), (d) and (e) as continuous and dashed black lines for RAB5A-overexpressing cells and control cells, respectively. The overlap parameter is calculated by assuming the threshold value a0 = 0.75R in panel (b), while the value a0 = 0.3R is adopted for obtaining the curves in panel (e). Here R = 4.7 μm is the average nuclear radius. Similarly, (d–f) show the collective migration velocity, the overlap parameters and the dynamic susceptibility obtained from analysis of phase-contrast movies of the same cells. DVA to obtain the curves shown in panels (b), (c), (e) and (f) is performed over the time window4–18 hours outlined in gray in panels (a) and (d), respectively.

To further explore these discrepancy, we compared our results with the ones obtained from a PT analysis performed on the same fluorescent image sequences. We evaluate the overlap parameter QPT(τ) for different values of the overlapping lengthscale a0, ranging from 0.1R to 2R, where R = 4.7 μm is the average nuclear radius (see Fig. S1 in ESI). As it can be appreciated from Fig. 3b and e an excellent agreement is found between Q(τ) and QPT(τ), provided that a suitable value of a0 is assumed for each imaging conditions. In particular, DVA results on fluorescence images for both RAB5A-overexpressing and control cells agrees well with PT results if the value a0,FLUO = 0.75R is used, whereas a smaller value (a0,PC = 0.3R) is required to obtain the best superposition on the DVA results for the same monolayer observed in phase contrast microscopy. We note the ratio between the overlapping lengthscales a0,FLUO/a0,PC ≃ 0.25 is in excellent agreement with the ratio lFLUO/lPC ≃ 0.24, showing that the “spatial sensitivity” (i.e. the lengthscale over which the dynamics is probed) of DVA is ultimately set by spatial correlation length of the considered set of images.

A similar quantitative agreement between our analysis and PT is not found for the dynamic susceptibility χ4, whose estimate with PT is systematically lower than its DVA counterpart (see Fig. S2, ESI). We attribute this discrepancy to the fact that a significant fraction of the nuclei cannot be tracked properly and cannot thus be included in the analysis. This leads to a systematic underestimate of size of the collectively rearranging regions, when estimated from PT.

3.3 Cellular unjamming in 3D: global rotation and local fluidisation

We have proved that our approach works effectively with cell assemblies moving over a 2D flat substrate. A case that is rapidly gaining popularity among cell biologists and biophysicists is represented by cell collectives in 3D settings, which are generally thought to better reproduce the micro-environmental and physiological conditions of living tissues.42,43

We thus apply our approach to a system of 3D spheroids of control and RAB5A-MCF10DCIS.com (a clonal breast cancer cell line) expressing mCherry-H2B that are embedded into a matrix of native collagen type-I. Experimental details can be found in Section 2.2 and in ref. 12. The global kinematic features of each spheroid are determined via registration of the corresponding sequence of 3D images with a rigid transformation (roto-translation). With the exception of small (typically subpixel) uncorrelated displacements, likely due to errors in the positioning of the microscope stage between consecutive frames, no net translational motion is detected. Our analysis (step I) extracts the angular velocity ω = [ωx,ωy,ωz] associated with the spheroid rotation along each of the three axes (x,y,z). Explicitly, ω is obtained from the transformation matrix M (see eqn (2)) as

image file: d0sm01837f-t16.tif

As it can be appreciated in Movie SM03_RAB5 (ESI), a persistent global rotational motion can be identified in all RAB5A-expressing spheroids. The angular speed |ω(t)| reaches its maximum around 12 hours from the beginning of the experiment and for almost its whole duration it remains above 0.5 rad h−1, a value that enable to complete a full rotation in about 12 h. (Fig. 4a).

image file: d0sm01837f-f4.tif
Fig. 4 Unjamming transition for a MCF10DCIS.com spheroid embedded in a thick collagen matrix. Temporal evolution of the angular speed ω(t) of RAB5A (a) and control (b) spheroids. Time autocorrelation function of the unit vector n identyfing the direction of the rotation axis of RAB5A (c) and control (d) spheroids. Overlap parameter Q(τ) obtain form 3D-DVA on RAB5A (e) and control (f) spheroids. For each conditions, N = 7 spheroids are considered. In each panel, results for single spheroids are shown as thin lines, while the thick line represents the average over all spheroids.

We find that, while in each spheroid the rotation axis shows no preferred orientation, its direction tends to strongly persist over time. This can be appreciated from Fig. 4c, where we show the time autocorrelation function of the unit vector n(t) = ω(t)/|ω(t)| that identifies the rotation axis. Over the investigated time window [0,8] h, the autocorrelation function decreases by less than 10%, indicating that the persistence time of the orientation largely exceeds the total duration of the experiments. Control spheroids do not show any global rotational motion. The average detected angular speed is typically of the order 0.03 rad h−1, corresponding to a tangential displacement on the spheroid surface of about 0.5 μm between consecutive frames, while the orientation of the axis of instantaneous rotation is almost delta-correlated in time (Fig. 4d).

In striking analogy with observations on 2D MCF10A cells, 3D-DVA of MCF10DICS.com spheroids shows that the global, collective migration promoted by RAB5A expression is accompanied by a marked increase in local motility and in the rate of cellular rearrangement (Fig. 4e and f). This “fluidization” process is well captured by the fast decay of the overlap parameter Q(τ) observed in RAB5A-expressing spheroids (Fig. 4e). An exponential fit of Q(τ) provides correlation times (τ* = 0.7 ± 0.2) h and (τ* = 2.6 ± 0.5) h for RAB5A and control spheroids, respectively.

Finally, we use sr-DVA to investigate whether RAB5A-induced fludization occurs uniformly across the spheroid. As described in Methods, a space-resolved map of the activity associated with mutual cell displacements can be obtained by considering the space-resolved overlap parameter q(τ0,x), defined in eqn (9). We choose the value τ0 = 0.3 h, corresponding to the delay time between consecutive 3D stacks. In Fig. 5a and b instantaneous intensity distribution on two orthogonal planes passing through the center of the spheroid are shown.

image file: d0sm01837f-f5.tif
Fig. 5 Space-resolved activity maps of RAB5A-overexpressing unjammed MCF10DCIS.com spheroids. xy (a) and yz (b) sections of a RAB5A-expressing spheroid. xy (c) and yz (d) sections of the time-averaged, space resolved variance d(x) on the same spheroid. xy (e) and yz (f) sections of the time-averaged, space resolved squared difference 〈d(τ*,x)〉t, with τ* = 0.3 h, on the same spheroid. (g) Azimuthally averaged radial profiles q(τ*|r) of the local overlap parameter, plotted as a function of the scaled radial coordinate r/R0, where r is the distance from the center of the spheroid and R0 is its radius. Each thin line corresponds to one of the N = 7 analyzed spheroids, while the thick line is an average. Panels (h–p) shows the same content for control spheroids.

As it can be also appreciated from panels (c) and (d), where the corresponding sections of the time-averaged, space resolved variance dinfty(x) are shown, the optical contrast is not uniform across the spheroid, but displays a significant radial gradient (indicating a stronger expression of the fluorescent tag in peripheral cells), as well as a systematic decrease with z due to the optical absorbance of the sample. Both effects are suitably accounted for by the space-resolved normalization of the differential map 〈d(τ|x,t)〉t (shown in panels (e) and (f)) prescribed by eqn (9).

The radial profile of the obtained local overlap parameter q(τ0|x) (Fig. 5g) displays a well-defined peak in correspondence of the outer cellular layer, revealing the presence of “melted region” at the spheroid periphery, while mutual rearrangement in the core are markedly slower. Such spatial organization of cellular dynamics is lacking in control spheroids, which show a rather uniform activity profile (Fig. 5p).

4 Conclusion

In this work, we have developed a new approach for the quantification of multicellular motility, based on extending the recently introduced Differential Variance Analysis27 for the analysis of the dynamics of cell collectives. Our approach – consisting in a two step procedure that accurately disentangles the collective motility (via image registration) from the local one (via DVA) – provides accurate results both with 2D cell assemblies as well as in more complex 3D scenarios, in which computational complexity and accuracy requirements are very different. Beside computational simplicity and speed, the main strength of our approach is robustness. Unlike cell tracking, which is often challenging and operator dependent, it does not require high quality images or selective labeling and, unlike PIV, which is optimized on extracting velocity fields, it can be effectively used to measure also more disordered forms of motility. We thus expect a wide applicability for soft and biological matter systems.

Code and data availability

The data and codes that support the findings of this study are available from the corresponding authors.

Author contributions

RC, GS and FG designed research. AP and EF prepared the samples and performed the experiments. SV performed the single cell tracking analysis of the 2D unjamming experiment. FG developed the DVA, wrote the codes and analyzed all the other data. RC and FG wrote the manuscript. All Authors contributed discussing the results and the manuscript.

Conflicts of interest

There are no conflicts to declare.


We thank R. Pastore for stimulating discussions on the DVA method. This work has been supported by: the Associazione Italiana per la Ricerca sul Cancro (AIRC) to F. G. and S. V. (MFAG#22083), and to G. S. (IG#18621 and 5XMille #22759); the Italian Ministry of University and Scientific Research (MIUR) to G. S. (PRIN: Progetti di Ricerca di Rilevante Interesse Nazionale – Bando 2017#2017HWTP2K).


  1. X. Serra-Picamal, V. Conte, R. Sunyer, J. J. Muñoz and X. Trepat, Methods in cell biology, Elsevier, 2015, vol. 125, pp. 309–330 Search PubMed .
  2. D. T. Tambe, C. C. Hardin, T. E. Angelini, K. Rajendran, C. Y. Park, X. Serra-Picamal, E. H. Zhou, M. H. Zaman, J. P. Butler and D. A. Weitz, et al. , Nat. Mater., 2011, 10, 469–475 CrossRef CAS PubMed .
  3. F. Serwane, A. Mongera, P. Rowghanian, D. A. Kealhofer, A. A. Lucio, Z. M. Hockenbery and O. Campas, Nat. Methods, 2017, 14, 181 CrossRef CAS PubMed .
  4. N. Hino, L. Rossetti, A. Marn-Llauradó, K. Aoki, X. Trepat, M. Matsuda and T. Hirashima, Dev. Cell, 2020, 53(6), 646–660 CrossRef CAS PubMed .
  5. R. Sunyer, V. Conte, J. Escribano, A. Elosegui-Artola, A. Labernadie, L. Valon, D. Navajas, J. M. Garca-Aznar, J. J. Muñoz and P. Roca-Cusachs, et al. , Science, 2016, 353, 1157–1161 CrossRef CAS PubMed .
  6. A. Shellard, A. Szabó, X. Trepat and R. Mayor, Science, 2018, 362, 339–343 CrossRef CAS PubMed .
  7. A. J. Ewald, A. Brenot, M. Duong, B. S. Chan and Z. Werb, Dev. Cell, 2008, 14, 570–581 CrossRef CAS PubMed .
  8. A. Mongera, P. Rowghanian, H. J. Gustafson, E. Shelton, D. A. Kealhofer, E. K. Carn, F. Serwane, A. A. Lucio, J. Giammona and O. Campàs, Nature, 2018, 561, 401–405 CrossRef CAS PubMed .
  9. J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim and B. Gweon, et al. , Nat. Mater., 2015, 14, 1040–1048 CrossRef CAS PubMed .
  10. C. Malinverno, S. Corallino, F. Giavazzi, M. Bergert, Q. Li, M. Leoni, A. Disanza, E. Frittoli, A. Oldani and E. Martini, et al. , Nat. Mater., 2017, 16, 587 CrossRef CAS PubMed .
  11. K. J. Cheung and A. J. Ewald, Science, 2016, 352, 167–169 CrossRef CAS PubMed .
  12. A. Palamidessi, C. Malinverno, E. Frittoli, S. Corallino, E. Barbieri, S. Sigismund, G. V. Beznoussenko, E. Martini, M. Garre and I. Ferrara, et al. , Nat. Mater., 2019, 1–12 Search PubMed .
  13. M. Raffel, C. E. Willert, J. Kompenhans, et al., Particle image velocimetry: a practical guide, Springer Science & Business Media, 2007 Search PubMed .
  14. W. Supatto, D. Débarre, B. Moulia, E. Brouzés, J.-L. Martin, E. Farge and E. Beaurepaire, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 1047–1052 CrossRef CAS PubMed .
  15. L. Petitjean, M. Reffay, E. Grasland-Mongrain, M. Poujade, B. Ladoux, A. Buguin and P. Silberzan, Biophys. J., 2010, 98, 1790–1800 CrossRef CAS PubMed .
  16. C. D. Meinhart, S. T. Wereley and J. G. Santiago, J. Fluids Eng., 2000, 122, 285–289 CrossRef .
  17. D. K. Vig, A. E. Hamby and C. W. Wolgemuth, Biophys. J., 2016, 110, 1469–1475 CrossRef CAS PubMed .
  18. S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Silberzan and N. S. Gov, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15314–15319 CrossRef CAS PubMed .
  19. F. Cichos, K. Gustavsson, B. Mehlig and G. Volpe, Nat. Mach. Intell., 2020, 2, 94–103 CrossRef .
  20. M. S. Manak, J. S. Varsanik, B. J. Hogan, M. J. Whitfield, W. R. Su, N. Joshi, N. Steinke, A. Min, D. Berger and R. J. Saphirstein, et al. , Nat. Biomed. Eng., 2018, 2, 761–772 CrossRef CAS PubMed .
  21. R. Cerbino and V. Trappe, Phys. Rev. Lett., 2008, 100, 188102 CrossRef PubMed .
  22. P. Edera, D. Bergamini, V. Trappe, F. Giavazzi and R. Cerbino, Phys. Rev. Mater., 2017, 1, 073804 CrossRef .
  23. F. Giavazzi, C. Malinverno, G. Scita and R. Cerbino, Front. Phys., 2018, 6, 120 CrossRef .
  24. R. Cerbino and P. Cicuta, J. Chem. Phys., 2017, 147, 110901 CrossRef PubMed .
  25. R. Cerbino, Curr. Opin. Colloid Interface Sci., 2018, 34, 47–58 CrossRef CAS .
  26. F. Giavazzi and R. Cerbino, J. Opt., 2014, 16, 083001 CrossRef .
  27. R. Pastore, G. Pesce and M. Caggioni, Sci. Rep., 2017, 7, 43496 CrossRef CAS PubMed .
  28. T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J. Fredberg and D. A. Weitz, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 4714–4719 CrossRef CAS PubMed .
  29. F. Giavazzi, C. Malinverno, S. Corallino, F. Ginelli, G. Scita and R. Cerbino, J. Phys. D: Appl. Phys., 2017, 50, 384003 CrossRef .
  30. F. Giavazzi, M. Paoluzzi, M. Macchi, D. Bi, G. Scita, M. L. Manning, R. Cerbino and M. C. Marchetti, Soft Matter, 2018, 14, 3471–3477 RSC .
  31. N. Lačević, F. W. Starr, T. Schrøder and S. Glotzer, J. Chem. Phys., 2003, 119, 7372–7387 CrossRef .
  32. L. Berthier, G. Biroli, J.-P. Bouchaud and R. L. Jack, Dynamical Heterogeneities in Glasses, Colloids, and Granular Media, 2011, vol. 150, p. 68 Search PubMed .
  33. A. A. Goshtasby, 2-D and 3-D image registration: for medical, remote sensing, and industrial applications, John Wiley & Sons, 2005 Search PubMed .
  34. B. Hall, Lie groups, Lie algebras, and representations: an elementary introduction, Springer, 2015, vol. 222 Search PubMed .
  35. M. Deforet, V. Hakim, H. Yevick, G. Duclos and P. Silberzan, Nat. Commun., 2014, 5, 1–9 Search PubMed .
  36. S. R. K. Vedula, M. C. Leong, T. L. Lai, P. Hersen, A. J. Kabla, C. T. Lim and B. Ladoux, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 12974–12979 CrossRef CAS PubMed .
  37. R. Cerbino, D. Piotti, M. Buscaglia and F. Giavazzi, J. Phys.: Condens. Matter, 2017, 30, 025901 CrossRef PubMed .
  38. F. Giavazzi, D. Brogioli, V. Trappe, T. Bellini and R. Cerbino, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 031403 CrossRef PubMed .
  39. L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D. El Masri, D. L'Hôte, F. Ladieu and M. Pierno, Science, 2005, 310, 1797–1800 CrossRef CAS PubMed .
  40. M. J. O'Sullivan, J. A. Mitchel, A. Das, S. Koehler, H. Levine, D. Bi, Z. D. Nagel and J.-A. Park, Front. Cell Dev. Biol., 2020, 8, 21 CrossRef PubMed .
  41. E. Lång, A. Połeć, A. Lång, M. Valk, P. Blicher, A. D. Rowe, K. A. Tønseth, C. J. Jackson, T. P. Utheim and L. M. Janssen, et al. , Nat. Commun., 2018, 9, 1–15 CrossRef PubMed .
  42. M. Simian and M. J. Bissell, J. Cell Biol., 2017, 216, 31–40 CrossRef CAS PubMed .
  43. M. J. Bissell, J. Cell Sci., 2017, 130, 3–5 CrossRef CAS PubMed .


Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sm01837f

This journal is © The Royal Society of Chemistry 2021