Non-invasive measurement of nuclear relative stiffness from quantitative analysis of microscopy data

Abstract The connection between the properties of a cell tissue and those of the single constituent cells remains to be elucidated. At the purely mechanical level, the degree of rigidity of different cellular components, such as the nucleus and the cytoplasm, modulates the interplay between the cell inner processes and the external environment, while simultaneously mediating the mechanical interactions between neighboring cells. Being able to quantify the correlation between single-cell and tissue properties would improve our mechanobiological understanding of cell tissues. Here we develop a methodology to quantitatively extract a set of structural and motility parameters from the analysis of time-lapse movies of nuclei belonging to jammed and flocking cell monolayers. We then study in detail the correlation between the dynamical state of the tissue and the deformation of the nuclei. We observe that the nuclear deformation rate linearly correlates with the local divergence of the velocity field, which leads to a non-invasive estimate of the elastic modulus of the nucleus relative to the one of the cytoplasm. We also find that nuclei belonging to flocking monolayers, subjected to larger mechanical perturbations, are about two time stiffer than nuclei belonging to dynamically arrested monolayers, in agreement with atomic force microscopy results. Our results demonstrate a non-invasive route to the determination of nuclear relative stiffness for cells in a monolayer. Graphic abstract


Introduction
In all the phases of their life cycle, biological tissues are dynamic entities: during morphogenesis, cells move and find their optimal location and shape in the organized spatial distribution of the tissue; in homeostatic equilibrium, dying cells are continuously extruded and replaced by new ones, but a variety of processes-such as tissue repair, inflammatory response, carcinogenesis, and tumor progression-can drive the tissue outof-equilibrium. In all these cases, the dynamic state of the tissue changes as a result of the appropriate physical and chemical stimuli being exchanged at different scales, from the one of the single cells to the one involving the entire tissue. However, the details of this process remains elusive.
With the eyes of a physicist, these changes of state of the tissue can be thought of as phase transitions that bring the tissue from one state to a different one. Focusing on cancer, it is possible to consider cancer progression toward metastasis and invasion as a sequence of a e-mail: roberto.cerbino@univie.ac.at b e-mail: fabio.giavazzi@unimi.it (corresponding author) steps involving a transition of the tissue from a solidlike to a fluid-like one; this fluidization transition is also referred to as unjamming [1][2][3][4][5][6][7][8]. The fluid tissue phenotype is typically accompanied at the single cell level by a softening of the cell, which may impact the infiltration of tumor cells [9,10], and at the same time, by an increased cell motility, which makes it easier for sick cells to invade surrounding healthy tissues. Since forces external to the cells propagate through the cytoplasm directly to the nucleus [11], larger stresses within the tissue resulting from unjamming can lead to nuclear envelope rupture and consequently lead to DNA leakage and damage [12][13][14][15]. In this framework, nuclear mechanical properties and shape are important regulators of the state of a tissue [16]; for example, nontrivial correlations between nuclear stiffness and tissue dynamics have been recently observed [17,18].
As a consequence of the complexity of considering cellular tissues as materials, an exhaustive characterization of their mechanical and rheological properties is extremely challenging. One of the main reasons is the intricacy of the cell, which results from the composition of several elements: cell membrane, cytoplasm, 50 Page 2 of 18 Eur. Phys. J. E (2022) 45 :50 cytoskeleton, nuclear envelope, organelles, nucleoplasm, chromatine, etc.. Moreover, cell mechanical properties depend on a variety of factors, including the cell status, the substrate, the treatments and the applied stress [19][20][21][22][23][24]. Modelization of the cell as an ensemble of few elastic and viscous elements is often incomplete and unrealistic [9]. To further complicate the picture, different techniques for the measurement of the cell viscoelastic properties have been deployed (e.g., AFM [19,20,25], micropipette aspiration [26,27], optical tweezers [28], magnetic twist cytometry [21]), and they generally obtain different values for these properties, as they often probe different cellular components [23,29].
Specifically for the cell nucleus, when considering its stiffness within the force field of a tissue one expects values that differ substantially from the ones obtained experimentally in controlled settings and, moreover, different cell lines can have very different mechanical properties, with cells in similar situations exhibiting nuclei either stiffer or softener than cytoplasm [19]. Recently, Parreira et al. [18] have shown how the presence of a single stiffer nucleus alters the dynamics of a whole monolayer, impairing collective migration. Through a refined image segmentation of nuclei, allowing the collection of several dynamical parameters, the Authors show that the decrease in net migration velocity is counterbalanced by an increase in the rotational dynamics of cell nuclei, accompanied by an increment of nuclear deformation. Interestingly, nuclear deformation affects all the cells and not only the first neighbors of the stiffer nucleus, thus suggesting that the presence of the obstacle influences nuclei far away from it. These results provide an important motivation for the development of non-invasive in situ tools for quantifying the nuclear stiffness as a proxy of mechanobiological alterations of the tissue.
In this work, we propose an automated procedure enabling the extraction of time-resolved information on morphology and dynamics of the nuclei in a monolayer, and the investigation of the interplay between nuclear shape fluctuations and collective motility within the cell aggregate. A custom algorithm is developed to identify, segment and track fluorescent nuclei over time in timelapse microscopy image sequences. The data are subsequently analyzed to investigate correlations between nuclear deformations and kinematics parameters characterizing the monolayer. Our estimators allow extracting non-invasively a robust estimate of the relative stiffness between the nucleus and the cytoplasm. We test our approach with different cell lines approaching cell jamming, as well as on tissues undergoing an unjamming transition via flocking induced by the overexpression of the RAB5A protein [4].

Sample preparation and imaging
MCF10A 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 (Gibco) 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% CO 2 . MCF10A 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 MCF10A cells with pBABE-puro-EGFP-H2B vector.
Transfections were performed using either calcium phosphate or FuGENE HD Transfection Reagent (cat. no. E2311, Promega), according to the manufacturer's instructions. Lentiviral and retroviral infections were performed as described if Ref. [4]. Cells were seeded in six-well plates in complete medium and cultured until a uniform monolayer had formed. For all the experiments, the seeding density was 1.5×10 6 cells per well, with the exception of the experiments indicated as LC (low confluency), for which the seeding density was 0.75 × 10 5 cells per well.
RAB5A expression was induced, where indicated, 16 h before performing the experiment by adding fresh complete media supplemented with 25 μgml −1 doxycycline to cells. At the time of recording, fresh media containing EGF and doxycycline was added. After cell induction, doxycycline was maintained in the media for the total duration of the time-lapse experiment.
An Olympus ScanR inverted microscope with a 5 × objective (experiments on jamming monolayers of MCF10.DCIS.com, MCF10AneoT, and MCF10CA1a cells) or a 10×objective (all the other experiments) was used to acquire images with a frame rate of 0.5 frames/min (experiments on jamming monolay-ers), 0.1 frames/min (experiment on MCF10A flocking monolayers) and 0.4 frames/min (experiment on MCF10.DCIS.com flocking monolayers) over a 24 h period. For each sample, 4 or 5 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 microscopy.

Nuclear tracking and segmentation
In order to extract both static and dynamic features of the monolayer at the single nucleus level, a Matlab algorithm was developed for the automatic nuclear segmentation and tracking. The algorithm operates on timelapse microscopy images of monolayers of cells with fluorescently tagged nuclei. The main steps of the algorithm are detailed below.

Image background subtraction and registration
Images are first corrected for background intensity inhomogeneities and microscope stage positioning errors.
Background subtraction is obtained through the implementation of the iterative procedure described in Ref. [30]. We first perform a cubic interpolation of the intensity in each image, by randomly choosing a set of points within a mask which initially coincides with the whole image. After the interpolation, the mask is redefined by excluding the points corresponding to high intensity values, which are likely to belong to a nucleus rather than to the background, and a second interpolation is performed based on the redefined mask. After an adequate number of iterations (in our case ten), the mask should include only regions belonging to background, evenly distributed across the image. The final interpolation step provides thus a reliable estimate of the background intensity distribution. This process is repeated over different frames equally spaced in time, covering the whole duration of the experiment. The average of all the so-obtained background images represents our best estimate for the background intensity distribution, which is thus subtracted to each image in the sequence.
Microscope stage translations between consecutive frames, due to positioning errors while imaging different locations at each time point, can introduce a bias in the reconstructed monolayer dynamics, by adding a spurious contribution to the velocity within the field of view (FOV). In order to address this issue, the Matlab intensity-based registration function imregister is used to recover the displacement Δx reg (t i ) of the FOV center of mass between each pair of consecutive frames (i−1) and i, considering pure translations. We then calculate the trajectory x reg (t i ) = i i =2 Δx reg (t i ) of the FOV center of mass. In the presence of both directed cell migration and stage movement, x reg can be written as x reg (t i ) = x cm (t i ) + x noise (t i ), where x noise (t i ) is a random, delta-correlated, noise associated with stage movement, and x cm (t i ) is the genuine displacement of the center of mass of the cell monolayer, which we esti-

Nuclear segmentation procedure
Once a registered stack of background-subtracted images I (t i ) is obtained, we process each frame i to identify single nucleus. Typical nuclei segmentation methods identify nuclei either localizing their centers (for example, looking for intensity maxima or centers of symmetry) or through the identification of nuclear edges [31]. In our case, the identification of the nuclei from their centers is made extremely difficult by the presence of intensity heterogeneity mirroring the chromatin nuclear distribution (Fig. 1a); we thus prefer to rely on nuclear edges rather than on nuclear centers. To identify nuclear edges, we apply a watershed transform to the image spatial gradient ∇I (t i ), which locates nuclear edges of I (t i ) as ridge lines. However, direct watershed segmentation is not very efficient in properly segmenting the nuclei, in particular in the case of jamming monolayers, where the signal-to-noise ratio is relatively low and partial superposition of two or more nuclei is relatively frequent. In order to improve the quality of the result, we introduce, before the application of watershed transform, a pre-processing where the image gradient ∇I (t i ) is set to zero in correspondence of suitable "seed points" located both within and outside the boundary of each nucleus. The procedure leading to the identification of the seed points is described in detail in the following.
Noise in each frame I (t i ) is first reduced by applying a Wiener filter, an adaptive noise-removal filtering that preserves nuclei edges [32]. A Laplacian-of-the-Gaussian (LoG) filter is then applied to enhance nuclear edges. In the obtained map L G (see Fig. 1b), higher values (toward yellow), corresponding to nuclear edges of the original image, surround deep intensity wells (blue), within which the nuclei centers of mass locate. Differences in the fluorescent intensity of different nuclei are reduced by dividing L G by an intensity map obtained via bicubic interpolation of the local minima of L G . The corrected L G is then binarized through the application of a suitable threshold k th , setting to 0 (1) the pixels whose intensity is larger (lower) than k th . The threshold value is determined as the one maximizing the number of domains of connected pixels equal to 1 (Fig.  2d). Indeed, as it can be appreciated by comparing Fig.  2a, b, when, starting from the minimum value of L G , the value of k th is increased, the number of connected domains initially increases as well. This is due to the fact that the correction based on bicubic interpolation of the local minima of L G does not perfectly level out the intensity differences between nuclei. On the other hand, when k th reaches the typical intra-nuclear value of L G , domains corresponding to different nuclei starts to merge and their number tends to decrease (Fig. 2bc).
Once the proper threshold k th is imposed, the resulting binarized image L BW (Fig. 1c) is an ensemble of domains of connected pixels, each one of which is a first approximation of the inner part of each nucleus.
On L BW , additional minor operations are then performed to correct false counting and spurious merging of neighboring nuclei. False counting are in most cases due to background noise and results in domains of few connected pixels. A morphological"open"operation (an erosion followed by a dilation) has been implemented to remove these artifacts.
Three examples of mistakenly merged nuclei are marked in red in Fig. 3. In order to correct this problem, the possible cases in which it may be occurred are automatically identified considering the domains with an aspect ratio larger than 3 (against a typical distribution of domains aspect ratios with median around 1.4 and standard deviation of 0.5) 1 . In order to roughly localize in them the centers of the merged nuclei, for each selected domain the two deeper local minima of L G are identified and two disks of radius 2 pixels are cen-1 The domains' aspect ratio limit has been determined through experience, has no real nuclei has been found in the analyzed data set with larger aspect ratios. It can be eventually adjusted for samples with different domains' aspect ratio distribution. The geometrical centers of mass of the domains of connected pixels are chosen as inner seeds for the watershed segmentation (red points in Fig. 1c).
The determination of the seeds external to the nuclei is simpler. The connected domains in L BW do not only indicate the position of the corresponding nuclei, but have also an area which is typically larger for nuclei with wider projected area. Consequently, the space between two nuclei in I (t i ) is typically equidistant from the edges of the corresponding connected domains in L BW . For the identification of the external seeds, we first determine the Euclidean distance transform L ED of L BW [33], a map where each pixel value is the distance from the closest domain (Fig. 1d). On L ED a watershed segmentation is finally applied to determine the equidistance lines between noighboring domains in L BW (yellow points in Fig. 1d), which are then used as external seeds.
In Fig. 1e, the recovered internal (red) and external (yellow) seeds are superimposed on the original image I (t i ). By setting ∇I (t i ) to zero in correspondence of all seeds and evaluating the watershed transform, the nuclei segmentation is finally obtained (Fig. 1d, green pixels). The recovered segmentation is pixel resolved; subpixel resolution can be achieved by running the Matlab code subpixelEdges 2.13 [34] on the obtained profiles ( Fig. 1e).
After nuclei segmentation, several static quantities are easily accessible. In this work, we consider in particular the center of mass x j , the projected area A j and the aspect ratio A j of the j-th nucleus. The latter is defined as the ratio between the nuclear major and minor axis, evaluated as the square root of the eigenvalues of the covariance matrix of the segmented objects [35].

Dynamic parameters
Time-resolved information on nuclear shape and mobility is recovered by linking the positions x j of each nucleus in subsequent frames. To this aim, we employed the publicly available Matlab implementation by D Blair and E Dufresne [36] of the linking algorithm of JC Grier and DG Crocker [37]. Once nuclear trajectories are built, the time evolution of the single-nucleus parameters x j (t i ), A j (t i ), p j (t i ) and A j (t i ) can be also obtained. The instantaneous velocity on the j-th nucleus across frames i and Our definition of the instantaneous velocity in terms of the average velocity over a time interval equal to δt relies on the fact that, on this timescale, the cellular motion can be assumed to take place with approximately constant velocity. We checked this assumption by measuring the nuclear mean square displacement and estimating the characteristic persistence time of the ballistic-like motion observed for short time delays. As discussed in detail in "Appendix 1," for all the considered data sets, the persistence time was found to be larger than δt.
The velocity of the center of mass of the monolayer is evaluated as the instantaneous mean velocity of the where . . . j denotes the average over all the nuclei in the field of view. The amplitude of velocity fluctuations in the monolayer is evaluated as the root mean square velocity of the nuclei in the center of mass reference frame: Comparison between nuclei tracked automatically by the algorithm and manually by an operator reveals that the described segmentation procedure is effective in identifying about 80-90% and 90-95% of the nuclei present in each FOV for the jamming and the flocking monolayers, respectively. This difference in tracking efficiency reflects differences in noise level and spatial resolution between different experiments (see previous section for details). Despite the effort to minimize segmentation artifacts, such as multiple segmentation of the same nucleus and spurious merging of multiple nuclei, such segmentation errors occur, especially in those cases where the signal-to-noise ratio is low or partial superposition of different nuclei is frequent. To minimize the impact of segmentation errors on the analysis of nuclear features, a "quality check" has been implemented to select a subset of reliably segmented nuclei. To this end, we evaluate the total instantaneous intensity J j (t i ), which is obtained as the sum of the intensity of all pixel within the segmented area of the j-th nucleus at frame i. The instantaneous value J j (t i ) is then compared with its median evaluated over the previous 10 frames. If the difference between J j (t i ) and the median is larger than 10%, the segmentation of the j-th nucleus at frame i is considered unreliable, and the corresponding parameters are not included in the statistics. Trajectories which, after the application of this quality filter, lose more than 20% of frames are entirely excluded. The number of nuclei that "pass" the quality check varies from sample to sample between 600 and 5000, depending on the FOVs size and image quality.

Particle image velocimetry
Particle image velocimetry (PIV) of fluorescent microscopy images of confluent monolayers is performed by using the Matlab PIVLab software [38]. We choose an interrogation area with size slightly larger than the average inter-nuclear distance, corresponding to approximately 14 µm. Outliers in the reconstructed velocity field, whose components v x and v y exceed fixed threshold values, are identified and replaced with the median value of the velocity over neighboring grid points 3 . From the velocity field obtained from PIV, v CM (t) and v RMS (t) are recovered as for PT averaging over grid points coordinates x in place of j. Additionally, the local divergence of the velocity field ∇ · v (x, t) is computed from v x (x, t) and v y (x, t) using the Matlab divergence function.
For the same samples, the PIV algorithm has been tested on both fluorescent and phase contrast time lapses. As shown in Appendix 1, both the PIV analyses return velocity mean values and distributions in agreement with the ones obtained from PT, thus validating the obtained dynamics.

Jamming monolayers
A first set of experiments was performed on mature highly confluent monolayers undergoing a jamming transition, a progressive slowing down toward a dynamically arrested state [39][40][41].
We considered three distinct MCF10A-derived cell lines (MCF10.DCIS.com, MCF10AneoT, and MCF10 CA1a cells), each one of those seeded at two different densities (see Sect. 2.1 for details).

Confluent monolayers evolve toward a high density kinetically arrested state
In Fig. 4, we report the time evolution of the RMS velocity v RMS (Fig. 4a) and of the center of mass velocity (Fig. 4b) for the jamming monolayers, as obtained from PIV. Each color refers to a different cell line (thick and thin lines correspond to high and low seeding density, respectively). For each cell line and seeding density, 5 independent FOVs are considered.
The general decreasing trend of both the considered indicators toward a plateau follows the expected dynamic arrest characterizing jamming transition and is accompanied by a progressive reduction of the mean cell area, evaluated of the order of 30% 4 , as the number of cells grow.
The expected increasing density within the cell monolayer can be observed in   Besides these trends, we are here interested in understanding the interplay between the dynamics of the monolayer and changes in nuclear shape at the singlecell level.
On top of the monotonic trends described above, each nucleus undergoes shape fluctuations, as it can be seen in Fig. 6 where the contours of the same nucleus are reported for different times.
Two parameters we have direct access to in order to characterize nuclear shape fluctuations are the area and the aspect ratio. They appear to be strongly correlated. As it can be appreciated from Fig. 5f, where we report the time evolution of the average nuclear aspect ratio, the aspect ratio displays a monotonically increasing trend with time. This indicates that, while nuclear projected area reduces as the density increases, nuclear shape anisotropy tends to increase. In order to deepen this observation, we consider the association between the instantaneous aspect ratio A j and the instantaneous area A j at the single nucleus level. An example of the resulting scatter plot, binned in equally spaced intervals of area, is reported in Fig. 5g. Different colors refer to 5 subsequent and overlapping time intervals of 8 hours covering the whole duration (24h) of the experiment. A striking negative correlation is observed, proving that in nuclei smaller projected areas are systematically more anisotropic. This observation is in agreement with what observed in Ref. [42]. No great differences emerge for different times, suggesting that, at a first glance, the correlation between A and A j is independent on age and monolayer dynamics.
Since the two parameters are strongly correlated, in the following we only focus on the projected area, whose value can be determined with lower uncertainty.

Nuclear deformation rates
Due to the non stationarity of the mean projected area on the experimental time scales, deformation amplitudes are difficult to be uniquely retrieved, as the obtained value are strongly dependent on the way the mean dynamics of the projected area is subtracted. In analogy with the mean squared displacement, we therefore introduce the mean square strain as where the nuclear strain Δa j (τ |t) of the jth nucleus between time t and t + τ is calculated as Representative MSS (τ ) curves obtained for the DCIS HC monolayer are reported in Fig. 7a. Each curve refers to one of the five, partially overlapping time intervals the experimental window was divided into. The linearity of MSS (τ ) at low τ points to a diffusive-like evolution of the area at short time scales. For larger time delays, MSS (τ ) becomes sublinear. The limited experimental time window does not allow establishing whether an asymptotic plateau value ∼ 1 is eventually attained. In order to extract the key parameters characterizing nuclear deformation, we use an exponential model function MSS (τ ) = σ w +ȧ 0 τ c 1 − e −τ/τc to fit the data. The diffusive-like growth of the area fluctuations is captured by the model in the limit of small τ : MSS (τ ) ∼ȧ 0 τ , with a characteristic strain rateȧ 0 . The short-time regime is followed by an exponentiallike relaxation trend toward a plateau valueȧ 0 τ c . The model also includes an offset σ w accounting for ran- dom, delta-correlated noise associated with the determination of projected area. In Fig. 7d, we report the values ofȧ 0 obtained from the fitting procedure for all the jamming experiments, at different time points. A common feature that can be observed in all datasets is the progressive decrease in the nuclear strain rate over time, indicating that mechanical deformation of nuclei decreases as kinetic arrest is approached (see also Fig. 4).

Nuclear deformation rates correlate with the local dynamics of the monolayer
In order to investigate in more detail the interplay between nuclear deformation and cell motility, we consider the divergence of the velocity field as a suitable parameter to capture local density fluctuations (i.e. compressions or dilations) within the monolayer, as schematically illustrated in Fig. 8. We evaluate the divergence field by computing for each time t the divergence ∇·v (t, x) of the velocity field obtained from PIV analysis (see "Appendix 1" for details). We then esti-(a) (b) Fig. 8 a, b Overplot of the divergence field over two FOVs of an MCF10A cell monolayer. Divergence field is represented as a shaded colormap where the passage from negative to positive values of divergence is marked by the color progression from red to white. Arrows represent the velocity field. Blue shades highlight two segmented nuclei subjected to a negative (a) and to a positive (b) divergence mated the value ∇ · v j (t) of the divergence in the position x j (t) corresponding to the center of mass of the tracked nuclei via cubic interpolation. In Fig. 7b, we report the experimentally determined probability distribution functions (PDFs) of ∇ · v j (t) for the DCIS HC monolayer over the same time intervals considered also in panel a. We observe that the PDF of the divergence tends to become narrower and narrower as the monolayer ages. This can be also appreciated from Fig. 7e, where the standard deviation σ ∇ of the PDF of the divergence is reported as a function of time. The overall trend is very similar to the one found forȧ 0 .
The origin of the similarity between the average values reported in Figs. 7d, e can be investigated by considering the direct association between the projected area of the j-th nucleus at a given time and the corresponding instantaneous local value of the divergence ∇ · v j (t). To this end, we computed the instantaneous single-nucleus area strain rate, defined aṡ In order to probe the correlation betweenȧ j (t) and ∇ · v j (t), we evaluate the corresponding time crosscorrelation function Kȧ ,∇v (τ ) averaged over the nuclei in the FOVs. For the precise definition of Kȧ ,∇v (τ ) and for consideration and comparison with corresponding self-correlations, see appendix 1. In Fig. 18c of appendix 1, we report the cross-correlation function for both negative and positive values of time delay τ . As it can be appreciated from the figure, a marked cross-correlation peak, centered in τ = 0, is present. The characteristic width of this peak varies between 5 and 20 minutes, depending on the cell line and on the age. The fairly symmetric shape of the cross-correlation peak with respect to the vertical axis indicates that there is no systematic delay between divergence and nuclear deformation, at least within the experimental temporal resolution (2 minutes).
The evidence of an instantaneous correlation betweeṅ a j (t) and ∇ · v j (t) prompts us to further investigate it to understand its character and origin. For each one of the 8 hours time intervals, we therefore generate a scatter plot ofȧ j (t) versus ∇ · v j (t). To make the plot more readable, we introduced a uniform binning of the horizontal axis. For each bin, containing on average hundreds or thousands of points of the scatter plot, we report in Fig. 7c the average valueȧ n of the strain rate as a function of the average value ∇·v of the divergence. Remarkably, at least in the vicinity of the origin, the obtained curves display a fairly linear behavior with zero intercept, thus denoting a direct proportionality between the nucleus strain rate and the corresponding divergence of the velocity field the nucleus is subjected to. In order to rationalize the physical meaning of the proportionality constant m between ∇ · v anḋ a n , we note that, under the hypothesis of weak spatial dependence of the monolayer cell number density ρ c , the continuity equation for the monolayer density can be written as where a c is instantaneous cell area divided by its mean value andȧ c is the cell strain rate. The coefficient m can be therefore interpreted as the ratio between the nuclear area strain rate and the cell area strain rate. The reciprocal value e n = m −1 is related to the relative stiffness of the nucleus compared to the one of the entire cell: the larger is e n , the lower is the deformation per unit time of the nucleus compared to the one of the cell. Estimates of m for the different cell lines and ages are obtained via a linear fit of the correspondingȧ n versus ∇ · v in the neighborhood of the origin. The corresponding values of e n are reported in Fig. 7f. Interestingly enough, the relative stiffness has a time dependence which is rather different from the one displayed by both the divergence and the strain rate. While in all jamming monolayers both σ ∇ anḋ a 0 reach a stable plateau value after about 12 hours, the relative stiffness typically keeps on decreasing over time, with a roughly constant rate. As it can be appreciated from Fig. 7c, for large values of ∇·v a significant deviation from linearity is observed, withȧ n seemingly approaching a saturation value for both negative and positive strain rates. The ratio between nuclear strain rate and cell strain rate becomes thus lower and lower as the imposed strain rates increase.

In flocking monolayers, nuclei experience stronger deformations and are stiffer
To further investigate the impact of cell motility on nuclear deformation, we repeated the experiments using two different cell lines (MCF10A and MCDF10A. DCIS.com) whose dynamical state is perturbed by inducing the overexpression of the RAB5A protein, a master regulator of endocytosis [43]. Upon RAB5A overexpression, mature, almost completely kinetically arrested monolayers experience a dramatic reawakening of motility, characterized by highly coordinated directed migration (flocking) and by the presence of local cell rearrangements (fluidization). [4,44,45]. For both considered cell lines, we compare RAB5A overexpressing monolayers with the corresponding controls in order to investigate the impact of the flocking transition on nuclear deformation. As detailed in "Appendix 1," PIV analysis confirms the striking motility phenotype of RAB5A overexpressing monolayer previously reported [4,5]. In Fig. 9, we show the average nuclear strain rate (Fig. 9a) and the standard deviation of the divergence of the velocity field (9b) for the two cell lines.
As it can be seen, bothȧ 0 and σ ∇ are systematically larger in flocking monolayers, where the relative motion of cells is enhanced (leading to larger values of divergence) and, as a consequence, nuclei are subjected to larger stresses and deform more. Despite the significant difference inȧ 0 and σ ∇ between the two cell lines, a striking correspondence between the two series of experiments can be observed in the slopes of the linear regime of the scatter plots reported in Fig. 9d, e: as it can be clearly seen in Fig. 9c, values of e n are very similar for the corresponding conditions, with a systematic difference by a factor ∼ 2 between the RAB5A overexpressing monolayers and the control ones. This suggests that the nuclei of RAB5A overexpressing cells deform less, compared to the control ones, when subjected to the same local density variation. Theȧ n versus ∇ · v curves reported in panels 9d,e display some interesting features even beyond the linear regime. As far as MCF10A cells are concerned, control monolayers display a clear deviation from the linear regime when |ȧ n | is above ∼ 0.05 h -1 , while the corresponding curve for RAB5A overexpressing cells remains linear for |ȧ n | < 0.1 h -1 . On the other hand, for MCF10.DCIS.com cells, both control and RAB5A overexpressing monolayers display a nonlinear behavior at nuclear strain rates larger than ∼ 0.1 h -1 and, apparently, the nuclear strain rate seems to attain a plateau at |ȧ n | ∼ 0.2 h -1 . The occurrence of such a plateau suggests the presence of mechanisms preventing the nuclei to be deformed above a certain threshold. Moreover, the marked difference between the two cell lines is even more interesting considering that MCF10.DCIS.com is a tumoral cell line deriving from MCF10A [46]. Further investigation is required to clarify whether the expanded linear regime (and the wider dynamic range ofȧ n ) displayed by MCF10.DCIS.com control monolayers compared to the corresponding MCF10A ones could be related to the tumorogenic potential of the cells.

Model
Our in-plane observations do not allow discriminating whether the observed fluctuations in nuclear and cellular projected area can be assumed to occur at constant volume, with each change in the projected area corresponding to a height variation of opposite sign. In principle, our observations are compatible with an actual variation of cell and nuclear volume, mediated by fluid exchange between nucleus and cytoplasm and extracellular medium. Large volume fluctuations in multicellular aggregates, due to active water transport in and out of the cell, have been indeed reported. The complex interplay between single-cell volume fluctuations and collective monolayer dynamics has been investigated via a combination of velocimetry and singlecell segmentation, similar to one exploited also in the present work [47,48]. Although the reported timescale of these volume fluctuations (a few hours) [49] is typically much larger than the correlation time of the projected area fluctuations observed in our experiments (a few minutes), we cannot rule out the possibility that the observed behavior is actually due to a poroelastic-like response of the cell [50]. Even without making assumptions on the details of the underlying process, we can provide a description of the mechanical properties of cells and nuclei in terms of effective quantities, building on the experimentally accessible observables described in the previous sections. Given the limited temporal resolution of our experiments (2-10 min), we do not expect to be able to capture the full rheological response of the cells. In particular, we can hardly probe the process of viscoelastic relaxation upon stress application, which typically occurs on the scale of few tens of seconds [51]. We are then most likely sensitive to the elastic response of cell components. Indeed, the lack of delay betweeṅ a n and ∇ · v observed from their cross-correlation (see Fig. 18) is consistent with this hypothesis.
We introduce a simple phenomenological model (sketched in Fig. 10) for the mechanical response of the cell to in-plane homogeneous compressive or tensile stresses, where an elastic nucleus and an elastic cytoplasm, with constant elastic moduli E n and E cy , respectively, are connected in series. Within this simple model, we can derive an explicit relationship between the ratio E n /E cy of the elastic moduli and relative stiffness e n introduced in Sect. 3.1.4. For two elastic components in series, the total area deformation ΔA c is the sum of the deformations ΔA n and ΔA cy of the components, while the external stress σ applied on the system coincides with the one applied on each element σ ∝ (E n a n ) = (E cy a cy ) .
By introducing the ratio between the average projected nuclear area and the average projected cell area β .
= A n / A c , it is possible to write the total strain associated to the cell projected area a c = ΔA c /A c in terms of the nuclear a n and the cytoplasmic strain a cy , as a c = βa n + (1 − β) a cy .
Combining this equality with the time derivative of equation 7, we get A time-resolved estimate of β can be obtained directly from the images, by combining the average nuclear area (Fig. 5) with the average cell area, which can be simply estimated as the total area of the FOV divided by the number of identified cells. In Fig. 11a, we report the resulting time evolution of β for the jamming monolayers, evaluated over the same time intervals considered in panels 7d-f. We note that β takes slightly different values for different samples, and it is approximately time independent.
Once β is known, we can use Eq. 9 to obtain an estimate of the ratio between the elastic moduli. As it can be seen from Fig. 11b, DCIS moduli ratio uniformly decreases from values above unity to values less than 1. Their nuclear modulus is therefore initially larger than the one of cytoplasm, but the relation reverses as far as cell packing increases. Similarly, in MCF10AT initially both nucleus and cytoplasm have similar moduli but with time cytoplasm becomes stiffer than nucleus. Finally, in MCF10CA E cy is systematically larger than E n . Possibility to observe nuclei either with lower or In the inset of (d), it is reported a zoom of the scatter plot over the area where both samples exhibit a linear behavior Fig. 10 Schematic representation of a cell deforming under the action of an in-plane compressive stress. On the right, 1D model of the cell as a series of two elastic elements with different elastic moduli En and Ecy, corresponding to the nucleus and the cytoplasm, respectively. Corresponding areas An and Acy are represented as lengths of the springs in the 1D model. The total area Ac given by the sum of An and Acy higher elastic modulus is not surprising as already observed in the literature [19] Robustness of the presented results is granted by the small relative errors, obtained from the standard deviation of the mean evaluated over 5 different FOVs for each sample, and by the Through the same procedure, we estimate β (Fig.  12a) and E n /E cy (Fig. 12b) for the flocking monolayers. While in control samples the elastic modulus of the nuclei is always close to the one of the cytoplasm, in RAB5A overexpressing cells it is almost twice as large as the cytoplasmic one, for both MCF10A and MCF10.DCIS.com. In principle, these results would be compatible with both nuclear stiffening and cytoplasmic softening, as a consequence of RAB5A overexpression. However, AFM measurements of monolayer rigidity (which are known to probe mainly nuclear stiffening) have already documented in MCF10A cell line an increase by a factor two in the Young's modulus in RAB5A overexpressing confluent cells [4].

Conclusions
In this work, we have introduced and demonstrated an experimental non-invasive procedure aimed at extracting information on the relative stiffness of nuclei compared to the cytoplasm, for epithelial cell monolayers whose state is monitored during time-lapse microscopy experiments. The procedure involves segmentation and tracking of the fluorescently tagged nuclei, combined with PIV and PT analyses of the nuclear motion.
More specifically, space-resolved and time-resolved studies of the monolayer dynamics and the nuclear deformation allow the study of their interplay, from which we find that for small strain rates, a linear correlation holds between the local density changes, estimated via the divergence of the velocity field, and the nuclear strain rate. In this linear regime, consistent with an elastic response we extract information on the effective relative stiffness of the nuclei, as we show with experiments performed with both jamming and flocking monolayers. For jamming monolayers of three different cell lines, we find nuclear moduli close to the ones of cytoplasm; we also find that, as the monolayers jam over time and the dynamics lowers and internal agitation within the monolayer decreases, the nucleus soft-ens compared to the cytoplasm. In agreement with this observation, a set of experiments on flocking monolayers obtained by overexpression of the RAB5A protein in two cell lines shows that the nuclear relative stiffness is higher than in control samples, despite the fact that nuclei deform more. The difference in the relative nuclear moduli is significant (about a factor of two) and is surprisingly similar for the two cell lines.
For large strain rates, we observe nuclear stiffening, with strain rates eventually approaching a plateau as a function of the divergence of the velocity field. We interpret this result as a mechanoprotective response of the nuclei to oppose large stresses that arise from the increased internal agitation of the nuclei in a highly dynamic monolayer. Interestingly, our results are compatible with the view that the value of the nuclear elastic modulus for small strains is mainly determined by chromatin, whereas for larger strains the strainstiffening contribution of lamins dominates [52,53].
The possibility of non-invasive monitoring of the relationship between mechanical stimuli, nuclear deformations and nuclear relative stiffness in cell monolayers may have immediate biological applications. An early version of the methodology that we describe here was recently used in Ref. [54] to assess in real time the nuclear mechanical stress and response and the associated role in causing a long-term transcriptionaldependent phenotype as a consequence of the shortterm adaptive response to stress. We thus believe that our methodology represents a useful addition to the portfolio of non-invasive tools to characterize the mechanobiological response of tissues.
the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

A. Mean square displacement analysis
In this "Appendix," we discuss the mean square displacements (MSDs) associated with the trajectories obtained from nuclear tracking, as described in Sect. 2.2. For each set of experiments, the MSD is calculated within each of the five time windows considered in the main text. In Fig.  13a-c, we display averages over all five time windows as a function of the lag time τ . There, continuous and dashed black lines indicate the slope corresponding to a ballistic scaling (quadratic in τ ), and a diffusive one (linear in τ ), respectively. As it can be appreciated also visually, in all experiments a crossover from a ballistic-like scaling at short time delays to a diffusive-like scaling at larger time scales is observed. The characteristic time τp marking the transition between the two regimes, which corresponds to the persistence time during which the velocity of a given cell is constant, can be estimated by fitting to the data the function where v b is the ballistic velocity and σMSD is an offset to account for experimental noise [55]. The persistence times recovered from the fit are reported in Fig. 13d-f. As it can be seen, they are systematically larger than the time step between two acquired frames in the corresponding experiments, represented by horizontal dotted lines. This observation implies that, on the timescale δt corresponding to the delay between two consecutive images, cells are moving with good approximation with constant velocity. This confirms the validity of the assumption made in Sect. 2.2.3 that the instantaneous nuclear velocity can be reliably estimated as the average velocity between two consecutive frames.

B. Comparison of PIV results from fluorescent and phase contrast time lapses
In this "Appendix," we compare the dynamics parameters obtained from PT and from PIV made over fluorescent and phase contrast time lapses. We then deepen the discrepancies between the velocities and velocity divergences recovered from PT and from PIV analysis of fluorescent and phase contrast time lapses, in order to motivate the choice of the fluorescent PIV analysis to evaluate the divergence of the velocity field as an indicator of the local density fluctuations within the monolayer. In Fig. 14 are reported the velocity distributions along the two main axis for a flocking MCF10A RAB5A overexpressing monolayer. The reported dynamics refers to the 8-9 hours time interval, corresponding to the cell motility peak [4]. Blue, orange and green lines refer to PIV computed over phase contrast and fluorescent time lapses and to PT, respectively. The recovered distributions are very similar to each other, thus validating the obtained dynamics.
Mean values obtained from PT and PIV are also very close, as shown in Fig. 15 where is reported for the six jamming experiments the time evolution of vRMS as recovered from PIV made over fluorescent time lapses (Fig. 15a) and from PT (Fig. 15b).
Given the correspondence between the results obtained from PT and from PIV, it is possible to choose which of the two methods to use to compare velocity divergence with the nuclear deformations shown in Sect. 3.1.4. A first possibility would be to start from the single-nucleus velocities vj (ti) to build up frame by frame a velocity field through a proper interpolation of the known velocities in the plane. This choice has been discarded in order to avoid spurious correlations deriving from the common origin from the localisation and segmentation of nuclei of the correlating data (divergence of the velocity field and time derivative of the area). We therefore chose to use a divergence obtained from PIV to ensure the authenticity of any correlation. Additionally, even a 10% missed nuclei in PT while allowing a fair enough agreement with PIV results in the velocity distribution (Fig. 14), at the same time prevents a proper interpolation of the velocity field. Test made revealed a discrepancy even by a factor 2 in the width of the ∇ · v distributions A deeper analysis revealed that even between the velocity fields recovered applying PIV analysis on fluorescent and phase contrast images have some discrepancies. We evaluate the normalized space cross-correlation between the velocity fields moduli v Fluo (x) and v PhC (x) obtained from the PIV analysis of the two data images, respectively: where • x stands for the spatial averages over all the positions where the velocity fields are evaluated. In Fig. 16b, it is reported the azimuthal average of Kv Fluo ,v PhC together with the space self-correlations of v Fluo (x) and v PhC (x). The spatial scale of decorrelation is similar in the three cases, but the cross-correlation value in r = 0 is lower by 20% than the self-correlations, thus testifying a weak but existent difference in the velocity fields.
Similarly, in Fig. 16d we report the azimuthal average of the cross-correlation relative to the divergences of the velocity fields defined as: As in the previous case, self-correlations of ∇ · v Fluo (x) and ∇ · v PhC (x) are reported in the same plot. For the divergences, the cross-correlation in r = 0 is much lower than the correspondent self-correlations, thus pointing out much larger differences between the divergence fields evaluated with the two methods. Discrepancy between the fields of ∇ · v much more pronounced than the ones in the corresponding v can also be appreciated comparing Fig. 16a, c, where the relative difference between the fluorescent and the phase contrast analysis is represented for the modulus of the velocity field and the divergence, respectively. While for v the discrepancies only emerges in the regions where there is a low density of nuclei (dark blue regions in Fig.  16a), very different appear to be ∇ · v fluo from ∇ · v phc . We  note that such differences are only present when spatial resolution is considered. Distributions of the divergence, frame by frame, are very close in the two cases. Given the differences, a choice has to be made between the two methods for the comparison with the nuclear strain rates. To this aim, we report in Fig. 17 scatter plots of the single-nucleus velocity components vx,j (Fig. 17a, c) and vy,j (Fig. 17b, d) obtained from PT as a function of the corresponding components of the cubic interpolation of the velocity field v fluo (Fig. 17a, b) and v phc (Fig. 17c, d) in the center of mass of the same nuclei. Looking at dispersions, one can find out that correspondence is better between PT and PIV analyses of fluorescent time lapses. This can be more easily sensed by looking at the width of the distributions along their minor dimension. This can be made by computing the minor eigenvalue of the covariance matrix between the velocity component obtained from PT and from PIV. The resulting widths are σx = 0.012 and σy = 0.012 for the fluorescent time lapses case and σx = 0.016 and σy = 0.017 for the phase contrast one.
Larger discrepancies from PT in the latter case can be partially rationalized considering that in phase contrast inner dynamics of the cells is also visible. Moreover, visual inspection of the videos revealed that several bright spots appear in phase contrast localized at the edges between adjacent cells. They move along cell edges faster than the cells themselves, thus affecting the velocity obtained from PIV.
In conclusion, because of this larger discrepancy between PT and phase contrast, we choose to compare area strain rate with the divergence obtained from PIV analysis of fluorescent time lapses.

C. Self-and Cross-correlations between local divergence and nuclear strain rate
In the main text, we report the normalized cross-correlation between nuclear strain rateȧj (t) and the corresponding divergence of the velocity field ∇ · vj. Normalized crosscorrelation for generic variables qj (t) and pj (t) is defined as In Fig. 18 are reported the nuclear strain rate (Fig. 18a) and divergence (Fig. 18b) self-correlations Kȧ,ȧ and K∇,∇ together with the cross-correlation Kȧ,∇ (Fig. 18c). In Fig.  18d the average of positive and negative τ of Fig. 18c is made. A strong decorrelation at the first delay can be observed in Kȧ,ȧ. This probably is a consequence of the high level of noise present on the determination ofȧ, and not of a real lack of time correlation of the strain rate. White noise origin of the strong time decorrelation is confirmed by the fact thatȧ and ∇v instead correlates in time Fig. 17 a, b Correlation between the velocities components recovered from PT and the ones obtained interpolating at the nuclei center coordinates the velocity field obtained from PIV on fluorescent time lapses. The minor widths of the distributions, evaluated as the lower eigenvalue of the covariance matrix, are σx = 0.012 and σy = 0.012 for vx and vy, respectively. c, d Correlation between the velocities components recovered from PT and the ones obtained interpolating at the nuclei center coordinates the velocity field obtained from PIV on phase contrast time lapses. The minor widths of the distributions are now σx = 0.016 and σy = 0.017 for vx and vy, respectively. Scatter plots are made considering the frames in the time interval 1-8 hours relative to one of the field of view of DCIS HC jamming experiments c Normalized cross-correlation between the instantaneous nuclear strain rateȧn and the divergence of the velocity field evaluated in the center of mass of the same nucleus. Data relative to jamming experiment DCIS HC in the time interval 0-24h. d Same data of c. but averaging over positive and negative τ . Black dashed line represents the best exponential fit detailed in the text. e Decay rate of Kȧ,∇ as a function of the considered time interval for the different experiments on jamming monolayers. Red, blue and green refer to DCIS, MCF10AT and MCF10CA monolayers, respectively. Thin and thick lines correspond to monolayers seeded at low and high concentration, respectively (18c), with a characteristic time comparable with the selfcorrelation characteristic time τK of ∇v (t, x) (Fig. 18b). The fact that an apparently delta-correlated signal-as the nuclear strain rate seems to be-correlates in time with the divergence might seem weird. As we show in the following, it is explained assuming that also nuclear strain rate has a finite time correlation hidden behind a strong deltacorrelated white noise. For a generic quantity q (t) whose real measurement q (t) = q (t) + σq is affected by the delta-correlated white noise σ q , the non-normalized self-correlation is: For the nuclear strain rate, it is possible to recover an indicative value of σp from the MSS as the square root of the offset σw. For the reported DCIS HC monolayer, √ σw = 0.07. If we compare this value with the square root of the non-normalized self-correlation evaluated in τ = 0, equal to 0.08, we perceive how small is the hidden genuine correlation compared to white noise. If the large white noise impedes the observation of time correlation of nuclear strain rate, it does not hinder information on the correlation with independently measured variables like divergence. Indeed, in cross-correlation between noisy datap andq with uncorrelated errors σp and σq it results: p (t) ·q (t + τ ) t = p (t) · q (t + τ ) t (15) as other terms vanish because of the delta correlation of noises σp and σq.
Oscillations at negative values of correlation at finite τ in Kȧ,ȧ and K∇,∇, and consequently in Kȧ,∇, can be rationalized considering that, as [56]: in order to reach a plateau of MSS for τ → ∞ the integral of the time correlations must be zero and thus correlations must assume negative values.
From the cross-correlation, we recovered the characteristic time τ k through an exponential fit excluding point in τ = 0. Inclusion of τ = 0 results in a failure of both exponential and stretched exponential fits. Apparently, a first very fast decorrelation arrived at time scales not accessible in our experiments. An example of the fit is superposed to data in Fig. 18d, and the resulting decorrelation rates 1/τ k are reported in Fig. 18e as a function of time for the experiments on jamming monolayers.

D. Dynamics of flocking monolayers
In Fig. 19 are reported the time evolution of vRMS (Fig. 19a) and the modulus of vCM (Fig. 19b). As expected [4,44], both MCF10A and MCF10.DCIS.com cell lines exhibit a reawakening of motility as a consequence of RAB5A overexpression, which is characterized by the emergence of a collective migration as a consequence of mutual velocity aligment of neighboring cells-testified by the motility peak in vCMand by the acquisition of a fluid-like dynamical state in the center of mass reference system, as evident from the dramatic difference in vRMS between RAB5A overexpressing monolayers and control ones.