Roberto
Cerbino
*^{a},
Stefano
Villa
^{a},
Andrea
Palamidessi
^{b},
Emanuela
Frittoli
^{b},
Giorgio
Scita
^{bc} and
Fabio
Giavazzi
*^{a}
^{a}Università degli Studi di Milano, Dipartimento di Biotecnologie Mediche e Medicina Traslazionale, 20090 Segrate, Italy. E-mail: cerbino@gmail.com; fabio.giavazzi@unimi.it
^{b}IFOM-FIRC Institute of Molecular Oncology, Milano, Italy
^{c}Università 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.

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 generated^{9} 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}

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}

MCF10A monolayer at low confluency.
Cells were seeded in six-well plates at 0.125 × 10^{6} 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 × 10^{6} 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 × 10^{3} 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 NaHCO_{3} 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% CO_{2}.

Fig. 1 Schematic representation of the main steps of the image processing pipeline. (a) Snapshots of the same cell monolayer acquired at different times t_{0}, t_{0} + Δt,…,t_{0} + τ_{1} − Δt,t_{0} + τ_{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(t_{0},τ_{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},t_{0}) = |I(x,t_{0}) − I(T(t)x,t_{0} + τ)|^{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},t_{0}) 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.

where 0 = [0,0], u = [u_{x},u_{y}] is a translation vector and M is 2 × 2 matrix. Exploiting the so-called polar decomposition,^{34}M 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 θ:

where R(θ) is a clockwise rotation by an angle θ

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),_{1}(t),_{2}(t),(t),(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

where Δt is the time delay between consecutive images.

where the symbol ○ indicates a composition of operators that becomes a matrix product for the corresponding matrices.

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 a_{0}. T can be represented by an extended 3 × 3 matrix

(1) |

(2) |

Accordingly, T is univocally determined by the set a = {u,λ

(3) |

In particular, provides an estimate for the velocity associated to directed migration, is the instantaneous angular velocity characterizing the in-plane rotation of the tissue, while (ν = 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,τ) = T[ā(t + τ)]○T[ā(t + τ − Δt)]○…○T[ā(t + Δt)]○T[ā(t)] | (4) |

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,τ)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

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) for all the considered values of τ and t.

where the symbol 〈f〉_{x} indicates a spatial average . 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[ 〈I^{2}(x,t)〉_{x,t} − 〈I(x,t)〉_{x,t}^{2}], where 〈f〉_{t} indicates a time average. This formula holds for ergodic samples also when the system did not relax completely.

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.

d(τ|x,t) = [I(x,t + τ) − I((t,τ)·x,t)]^{2} | (5) |

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

(6) |

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: , 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 d → d − 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 C_{I}(r) = 〈I(x_{0} + r,t)I(x_{0},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}(τ) = N〈q^{2}(τ|t)〉_{t} − 〈q(τ|t)〉_{t}^{2}. | (8) |

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

where the normalization factor d_{∞}(x) = 2[〈I^{2}(x,t)〉_{t} − 〈I(x,t)〉_{t}^{2}] 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).

(9) |

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 x_{n}(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 v_{n}(t) = [x_{n}(t + Δt) − x_{n}(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 . Let 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 , where a_{0} is a parameter controlling the overlap length-scale. Qualitatively, w_{n}(τ|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 a_{0}, while w_{n}(τ|t) ≃ 0 when the displacement largely exceeds a_{0}. The time-dependent overlap parameter q_{PT}(t,τ), is then obtained by averaging the overlap function over all the N tracked nuclei: , while the time-averaged overlap parameter are defined in close analogy with eqn (7) and (8): Q_{PT}(τ) = 〈q_{PT}(τ|t)〉_{t} and χ_{4,PT}(τ) = N〈q_{PT}^{2}(τ|t)〉_{t} − 〈q_{PT}(τ|t)〉_{t}^{2}. We note that according the above definitions, Q_{PT}(τ) can be interpreted as the average fraction of nuclei that, during a time interval of duration τ, displace from their initial positions by less than a_{0}.

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 v_{0} 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.

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 l_{FLUO} = 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 l_{PC} = 2.0 ± 0.2 μm.

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 a_{0} = 0.75R in panel (b), while the value a_{0} = 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 window^{4–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 Q_{PT}(τ) for different values of the overlapping lengthscale a_{0}, 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 Q_{PT}(τ), provided that a suitable value of a_{0} 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 a_{0,FLUO} = 0.75R is used, whereas a smaller value (a_{0,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 a_{0,FLUO}/a_{0,PC} ≃ 0.25 is in excellent agreement with the ratio l_{FLUO}/l_{PC} ≃ 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.

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

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).

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.

As it can be also appreciated from panels (c) and (d), where the corresponding sections of the time-averaged, space resolved variance d_{infty}(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).

- 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 .
- 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 .
- 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 .
- 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 .
- 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 .
- A. Shellard, A. Szabó, X. Trepat and R. Mayor, Science, 2018, 362, 339–343 CrossRef CAS PubMed .
- A. J. Ewald, A. Brenot, M. Duong, B. S. Chan and Z. Werb, Dev. Cell, 2008, 14, 570–581 CrossRef CAS PubMed .
- 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 .
- 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 .
- 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 .
- K. J. Cheung and A. J. Ewald, Science, 2016, 352, 167–169 CrossRef CAS PubMed .
- 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 .
- M. Raffel, C. E. Willert, J. Kompenhans, et al., Particle image velocimetry: a practical guide, Springer Science & Business Media, 2007 Search PubMed .
- 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 .
- 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 .
- C. D. Meinhart, S. T. Wereley and J. G. Santiago, J. Fluids Eng., 2000, 122, 285–289 CrossRef .
- D. K. Vig, A. E. Hamby and C. W. Wolgemuth, Biophys. J., 2016, 110, 1469–1475 CrossRef CAS PubMed .
- 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 .
- F. Cichos, K. Gustavsson, B. Mehlig and G. Volpe, Nat. Mach. Intell., 2020, 2, 94–103 CrossRef .
- 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 .
- R. Cerbino and V. Trappe, Phys. Rev. Lett., 2008, 100, 188102 CrossRef PubMed .
- P. Edera, D. Bergamini, V. Trappe, F. Giavazzi and R. Cerbino, Phys. Rev. Mater., 2017, 1, 073804 CrossRef .
- F. Giavazzi, C. Malinverno, G. Scita and R. Cerbino, Front. Phys., 2018, 6, 120 CrossRef .
- R. Cerbino and P. Cicuta, J. Chem. Phys., 2017, 147, 110901 CrossRef PubMed .
- R. Cerbino, Curr. Opin. Colloid Interface Sci., 2018, 34, 47–58 CrossRef CAS .
- F. Giavazzi and R. Cerbino, J. Opt., 2014, 16, 083001 CrossRef .
- R. Pastore, G. Pesce and M. Caggioni, Sci. Rep., 2017, 7, 43496 CrossRef CAS PubMed .
- 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 .
- F. Giavazzi, C. Malinverno, S. Corallino, F. Ginelli, G. Scita and R. Cerbino, J. Phys. D: Appl. Phys., 2017, 50, 384003 CrossRef .
- 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 .
- N. Lačević, F. W. Starr, T. Schrøder and S. Glotzer, J. Chem. Phys., 2003, 119, 7372–7387 CrossRef .
- 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 .
- A. A. Goshtasby, 2-D and 3-D image registration: for medical, remote sensing, and industrial applications, John Wiley & Sons, 2005 Search PubMed .
- B. Hall, Lie groups, Lie algebras, and representations: an elementary introduction, Springer, 2015, vol. 222 Search PubMed .
- M. Deforet, V. Hakim, H. Yevick, G. Duclos and P. Silberzan, Nat. Commun., 2014, 5, 1–9 Search PubMed .
- 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 .
- R. Cerbino, D. Piotti, M. Buscaglia and F. Giavazzi, J. Phys.: Condens. Matter, 2017, 30, 025901 CrossRef PubMed .
- F. Giavazzi, D. Brogioli, V. Trappe, T. Bellini and R. Cerbino, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 031403 CrossRef PubMed .
- 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 .
- 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 .
- 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 .
- M. Simian and M. J. Bissell, J. Cell Biol., 2017, 216, 31–40 CrossRef CAS PubMed .
- M. J. Bissell, J. Cell Sci., 2017, 130, 3–5 CrossRef CAS PubMed .

## Footnote |

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

This journal is © The Royal Society of Chemistry 2021 |