Measuring cell displacements in opaque tissues: dynamic light scattering in the multiple scattering regime

: Coherent light scattered by tissues brings structural and dynamic information, at depth, that standard imaging techniques cannot reach. Dynamics of cells or sub-cellular elements can be measured thanks to dynamic light scattering in thin samples (single scattering regime) or thanks to diﬀusive wave spectroscopy in thick samples (diﬀusion regime). Here, we address the intermediate regime and provide an analytical relationship between scattered light ﬂuctuations and the distribution of cell displacements as a function of time. We illustrate our method by characterizing cell motility inside half millimeter thick multicellular aggregates.


Introduction
In the past two decades, we have witnessed an unprecedented development of optical microscopy for biology. On the one hand, the advent of super-resolved microscopy [1][2][3] allowed to break the diffraction barrier. On the other hand, new techniques were introduced to allow optical sectioning of thick samples, including confocal and two photon microscopy [4,5], light sheet microscopy [6][7][8], optical coherence tomography [9][10][11][12] and phase conjugation [13,14]. Apart from the latter, these techniques are based on the signal carried by ballistic photons, those that have not been scattered in their path through the tissue. Unfortunately, cells have a very large scattering cross section that severely limits the mean free path of photons (called s hereinafter) to few tens of microns, in biological tissues [15,16]. As the non-scattered photons decrease exponentially with the penetration depth, the signal becomes unexploitable beyond a thickness corresponding to few cells [17].
To go beyond this limit, several techniques based on the exploitation of the scattered photons have been developed. Even for very thick tissues, fluctuations of coherent scattered light are informative of the dynamics of the scatterers population, though individual behaviors cannot be measured. In laser Doppler flowmetry [18,19] (LDF), scatterers are red blood cells and blood flow is calculated from the fluctuations of backscattered light. Laser speckle flowmetry [20,21] (LSF) works similarly, but the blood flow is evaluated from the speckles contrast relatively to the exposure time. Dynamic light scattering (DLS) [22,23] is generally used to measure the diffusion coefficient of colloidal suspensions and was recently applied on sub-cellular components [24][25][26]. All three techniques rely on the same physical principles, although the formalism varies as they were developed in different communities relatively to their field of application. Light fluctuations are analyzed either by mean of the power spectrum of the intensity over time (LDF), or by mean of the speckle contrast (LSF), related to the power spectrum as explained in [27], or by mean of the autocorrelation function of the intensity over time (DLS). As LDF and LSF were developed to measure in vivo blood flow, a backscattering configuration (scattering angle of 180 • ) was adopted. On the contrary, in DLS an additional information can be retrieved from the angular dependence of the signal. To obtain similar information with backscattering, the acquisition must be performed on several axes [28].
Relating light fluctuations to a quantitative measurement of the scatterers dynamics is made difficult by multiple scattering (MS). Several methods have been proposed to reject multiply scattered photons if they represent a small fraction of the scattered light [29][30][31][32]. Low-coherence interferometry proved to be efficient in selecting photons depending on their path-length, allowing single scattering conditions either by setting a path-length smaller than s or by use of confocal optics at a depth of the selected path-length (DLS-optical coherence tomography) [33]. Apart from technical solutions, MS can be integrated into the model used to relate light fluctuations and scatterers movements. For the specific case of red blood cells moving behind a highly diffusive tissue (LDF,LSF), a model was proposed by Bonner et al. [34]. For a backscattered configuration (LDF), a correction of the power spectrum under MS was reported by Wax et al. [35]. As for DLS, diffusive wave spectroscopy [36,37] (DWS) was developed as an extension to optically thick samples. More precisely, DLS is only valid if the sample size L< s because photons are assumed to be scattered only once, while DWS applies to samples thicker than the transport mean free path * [38]. * is the average distance after which photons are isotropically oriented [see Eq. (10)]. In biological tissues, * is roughly 10 times larger than s [39] because the scattering is mainly forward. DWS was applied to biological samples [40,41] and was also used to perform dynamical tomography by acquiring scattered light in numerous directions and numerically solve an inverse problem (Diffuse Correlation Tomography [42], Speckle Contrast Optical Tomography [43]).
In aforementioned models, the angular dependence of the signal was not studied, either because only backscattered light was collected or because the signal was isotropic in the L> * regime. However, the intermediate case ( s <L< * ) is rather frequent and the forward-scattering configuration (or small-angle scattering) is more beneficial than backscattering for two reasons. First, multiply scattered photons are fewer because in tissues light is mainly scattered in the forward direction [39]. Secondly, the signal can be measured at different angles, bringing more information than backscattering. Using numerical resolution, some approaches address the lower limit of DWS [44,45] or even bridge DLS and DWS regimes [46].
In this paper, we derive the relationship between scatterers displacements and scattered light fluctuations in the regime L< * , thanks to an analytical model of the multiply scattered photons propagation. Our method, denoted Multiple Scattering DLS (MS-DLS) in the following, is suitable for small-angle scattering in the intermediate regime and does not require time-consuming simulations. The angular dependence allows the measurement of the displacements distribution over time, without having to assume any type of motion (ballistic or Brownian). The scatterers that interested us were cells rather than sub-cellular elements studied in [24][25][26]. This is consistent with a measurement at small angles where the scattered light is dominated by the cells as a whole rather than its sub-cellular elements [47]. Using DLS on cell monolayers, we previously showed that cell migration speed could be measured without the need for usual cell tracking procedures [48]. Thanks to MS-DLS, we extend here the measurement of cell displacements over time to spherical multicellular aggregate (Multi Cellular Spheroid -MCS) where MS cannot be neglected. Our experiments revealed a decrease in the average cell speed from the periphery toward the center, as well as preference in the tangential direction for cells motion. We also found that an isotropic compression of the MCS inhibits cell motility, probably by limiting the access to the intercellular space, filled by extracellular matrix and interstitial fluid.

Sample preparation
Colon carcinoma cells CT26 (ATCC CRL-2638 TM ; ATCC, Manassas, VA) were cultured in 25 cm 2 sterile flask (VWR, Radnor, PA) containing 6 to 10 mL of DMEM with GlutaMAX (THERMO FISHER, Waltham, MA) supplemented with 10% of fetal bovine serum and 1% Penicillin-Streptomycin (SIGMA-ALDRICH, St. Louis, MO). Incubated at 37 • and 5% CO 2 , cells were passaged using trypsin every 2 to 3 days to keep the cell surface concentration between 30 and 1000 cells/mm 2 . Spheroid were prepared in 96 wells plates coated with 1% agarose gel [49]. 300 cells suspended in a 150 µL culture medium solution were seeded in each well, leading to the formation of a L = 250 ± 30 µm (standard deviation) diameter spheroid after 3 days of incubation ("small" spheroids with ∼3000 cells) or a L = 470 ± 50 µm diameter spheroid after 7 days of incubation ("big" spheroids with ∼20000 cells). To apply an isotropic pressure on the spheroid, we used the osmotic stress exerted by Dextran molecules (SIGMA-ALDRICH) in solution (r>15 nm, MW>200 kDa) as previously described in [50,51]. Dextran molecules are excluded from the spheroid volume and exert a pressure at its surface, following the empirical relationship indicated in Table 1. For each experiment, a single spheroid was put into a 35 mm VWR round petri dish containing 2 mL of culture medium with or without Dextran. To prevent translational motion of the spheroid, the petri dish was treated beforehand with Pol-L-lysine (SIGMA-ALDRICH) for 10 minutes and rinsed with phosphate-buffered saline (PBS), improving the adhesion of the spheroid to its surface.

Experimental setup
The light scattering setup was implemented on a commercial microscope (Axiover 100; ZEISS, Jena, Germany) in order to switch easily between the phase contrast image of the sample and scattered light image [see Fig. 1]. Monitoring both images simultaneously is possible, as described in [48], but it would lower the quality of the light scattering signal. Implementing a scattering detection over a microscope have been previously done in numerous works, such as [52], [24] and [53]. The microscope was equipped with a phase-contrast condenser, a 10× objective (EC Plan-Neofluar; ZEISS, NA = 0.3) and a Charge-Coupled Device camera (CCD1 -Pike; ALLIED VISION, Stadtroda, Germany). The sample was illuminated using a 850 nm laser (CPS850V, 0.85 mW; THORLABS, Newton, NJ). The beam was attenuated (2.0 ND; THORLABS) and narrowed with two lenses (L1: +150 mm, L2: +50 mm) to obtain a Gaussian beam with a waist of 192 ± 12 µm. The front focal plane of L3 lens (+60 mm, LA1401-B; THORLABS) [see Fig. 1] contains the scattered light image formed at infinity (Fraunhofer scattering) and is imaged on the CCD2 camera (Guppy F 046B; ALLIED VISION) using L4 lens (same as L3). L3 was positioned as close as possible to the sample to maximize the range of collected scattering vectors q: from 0.01 µm −1 (laser divergence) to 1.7 µm −1 . L4 was positioned to adjust this q range to the CCD detector size, leading to an angular resolution of 0.0014 µm −1 per pixel. Angular calibration was performed using a diffraction grating. To maintain temperature, humidity and pH, the sample was put inside a top-stage incubator system (CU-501; Live Cell Instrument, Seoul, South Korea). We choose to keep the scattering setup simple but one could use low-coherence interferometry and confocal optics to have a z-resolve measurements [33]. Light scattering dectection was added to a commercial microscope. The illumination part was composed of a laser (850 nm) attenuated (neutral density) and narrowed (L1 and L2). The detection part was composed of two lenses L3 and L4 positioned so that the front focal plane of L3 is imaged on the CCD2 camera.

Data acquisition and analysis
Data acquisition and analysis were performed using MATLAB (THE MATHWORKS, Natick, MA). Scattered light was recorded using CCD2 at one frame every 20 seconds, over a period of at least 2.5h. To capture the wide range of scattered intensity, each frame was composed of 5 concentric rings taken at different time exposures. From largest scattering angles, where the intensity was lower, to smallest angles, exposure times were: 512 ms, 128 ms, 32 ms, 8 ms and 2 ms. To measure the phase function p [see Eq. (2)], the background light was subtracted and vignetting effects were corrected [54]. Then, the intensity was averaged over the azimuthal angle ϕ. More specifically, the average over ϕ was performed on a set of rings of scattering vectors q satisfying k∆q<|q|<(k + 1)∆q, k being an integer and ∆q the discretization step [see Fig. 2 As the speckle grain area was about 900 pixels, rings contained from 3 (inner ring) to 150 (outer ring) speckle grains.
When cells move inside the spheroid, each spot of the diffraction pattern fluctuates. Figure 2 (a) displays the temporal evolution of the intensity I(q, t), recorded from a given pixel at a scattering vector q. To quantify light fluctuations, we computed the intensity autocorrelation function C [see Fig. 2 Qualitatively, at a given |q|, the time needed for the correlation function C to reach zero is related to the time needed for one cell to travel a distance of 1/|q|. The more times photons are scattered, the shorter the time to zero is because photons probe a combination of several cells displacements. Different scattering patterns corresponds to different cell configurations inside spheroids. We assume that all different orientation configurations of spheroids can be explored while looking at a single direction over time (ergodic hypothesis) : I(q, t) ϕ = I(q, t) t . This hypothesis is reasonable because of the rotational symmetry of spheroids. The spheroid azimuthal symmetry also led to the same symmetry in the correlation function and C was averaged over ϕ [see Fig. 2 (b) and (d)]: C(q, τ) = C(q, τ) ϕ . In order to measure cells dynamics at steady-state, we waited until the correlation function became stationary. In practice we discarded the first tens of minutes of acquisition after MCS manipulation.

DLS extension to multiple scattering
Dynamic light scattering relies on the hypothesis of single scattering (Born approximation). However, this approximation is limited to thin samples as the average distance between two scattering events, the so-called mean free path s , range from 15 to 100 µm inside tissues [15,16].
To address the case of a small number of scattering events, intermediate between classical DLS and DWS, we modelled the propagation of light as corpuscular. While crossing a spheroid, a photon can be scattered several times by different cells [see Fig. 3 (a)]. Each scattering event is described as a change of direction, from u to u (unit vectors), following a probability proportional to the angular intensity I s of a single scattering event. This description ignores interference with other photons during the propagation following the scattering event. This is the so-called ladder or independent scattering approximation valid for a dilute medium such that the wavenumber satisfies k 0 s 1. The probability function I s corresponds to the phase function of a single cell and is defined as where k 0 is the wavenumber k 0 = 2π/λ. The phase function is usually defined as a dimensionless quantity. Here we chose to define it with the dimension of a surface (1/k 2 0 ) in order to have consistent dimensions in the following Eq. (3). The phase difference between each photon is calculated geometrically, according to its trajectory. Finally, the electrical field on the detector is taken as the addition of all electrical fields associated to photons. Based on this model, we deduced the expression of the intensity correlation function C(q, τ), defined in Eq. (1), with the delay τ and the scattering vector q = k 0 (u s − u i ) [see Fig. 3 (a)]. The calculus is detailed in the appendix. Briefly, all possible photon trajectories are integrated, with weights given by the angular distribution p s and the proportion, noted α k , of k-times scattered photons in the scattered light. The result gives the relationship between the correlation function C and the cell displacement distribution at a given delay, noted D(∆r, τ) : with where p s k p s denotes p s being convoluted with itself k times, with the convention p s 0 p s = p s ; D(q, τ) is the 3D Fourier Transform of D(∆r, τ) along q. p corresponds to the phase function of the whole spheroid, defined as in Eq. (2) by replacing I s by the average intensity scattered by the spheroid I(q, t) t . In the single-scattering approximation, Eq. (3) simplifies into: (3) is valid under the following hypotheses: i) light fluctuations are mainly due to cell displacements rather than cell shape changes, ii) cell displacements are uncorrelated, iii) displacements are small as compared to the mean free path, iv) the sample is homogeneous in terms of cell displacement behavior. Regarding hypothesis i), the smaller the scattering angle is, the less cell shape details impact the scattered pattern. Hypothesis ii) is a strong hypothesis in biological samples. However, measured speeds may still be interpreted as speeds of independent groups of cells. Finally, q only takes values on a surface, approximated at small angles as a plane perpendicular to the optical axis. As a result, D(∆r, τ) can only be measured on this plane and movements parallel to the optical axis cannot be detected.
To evaluate α k , we start from the probability for a photon to be scattered k-times while traveling a distance s inside the sample, which follows a Poisson distribution [55] (we neglect absorption over scattering): Then, the proportion α * k of photons scattered k-times by a spheroid of radius R can be calculated by integrating Eq. (5) over the cross section of the spheroid (polar coordinates r, θ), weighted by the incoming laser intensity I l : The traveling distance s(r) was approximated as the thickness of the spheroid (s(r) = 2 √ R 2 − r 2 ) because photons trajectories are nearly straight as the forward direction is lost after a distance of * = 950 µm, as measured hereafter. The integration over r goes beyond R to take into account the light that misses the spheroid (s(r) = 0 in this case).
Finally, to obtain the proportion α k of photons k-times scattered within the scattered light, the unscattered light was taken out: The mean free path s was deduced from the measurement of α * 0 , the ratio of unscattered light (central spot) with and without the sample. We found s = 76 ± 4 µm at a λ = 850 nm wavelength. Figure 3 (b) and 3 (c) shows α k computed for two spheroids, the diameters of which are respectively 250 µm and 470 µm. It shows that the Born approximation is not valid for our samples, hence the need of Multiple Scattering DLS.

Displacement distribution measurement
To illustrate our approach, we deduce the displacement distribution D(∆r, τ) of cells from the correlation function C(q, τ), inside an opaque multicellular aggregate. The correlation function is calculated from the evolution of the scattered intensity pattern, as explained in the method section. Then, to deduce the cell displacement distribution from Eqs. (3) and 4, we empirically chose a parametric form for the phase function p s , combining a Gaussian function for small angles and a Henyey-Greenstein function for larger angles [56]: where the dot product u · u is the cosine of the scattering angle. The parameters a, q 0 and b were adjusted so that p calculated using Eq. (4) fitted (least squares) the data points. The calculation is easier in Fourier space as convolution products become products. Fourier transforms were performed analytically when possible (Gaussian, Henyey-Greenstein). Otherwise, analytical expressions were discretized with a sampling rate and a range making the numerical error negligible as compared to data noise. Figure 4 (c) shows the adjustment of the measured phase function of a spheroid of 146 µm in diameter. The best fit was obtained for a = 0.23, q 0 = 0.14 µm −1 and b = 0.90. Incidentally, we deduced * from this measure of p s . First, we calculated the anisotropy factor g, defined as Then, * was obtained as * = Similarly, we propose an empirical parametric function for D(∆r, τ). A Gaussian distribution showed a relative agreement with the data [ Fig. 4 (a)], while preserving the azimuthal symmetry showed by the correlation functions C: where ∆r = |∆r|. σ τ is an adjustable parameter that scales with τ in a manner that depends on the nature of cell motion (subdiffusive, diffusive, ballistic). Because this approach is based on curve fitting, it can be used even when certain angles are not measurable. This is our case as, due to the direct beam, C is measurable only for q>0.09µm −1 .

Full multiple scattering model
To consolidate our model, we compared the prediction given by Eq. (3) to a full theoretical model of light transport in the MS regime. This model is presented in detail in Ref. [57]. Starting from first principles (i.e.generalized scattering matrix for a moving scatterer), a transport equation is derived directly for the temporal correlation function of the field g (1) . This equation, similar to the well-known Radiative Transfer Equation (RTE) reads in the steady-state regime: where h is the 1D Fourier transform of h(v) is the velocity distribution of the scatterers. It is related to the displacement distribution by the relation In Eq. (12), g (1) (r, u, τ) is the Wigner transform of the field time correlation function [57]. It gives access to the time decorrelation of the field at a given point r and for a given direction u inside or outside the medium. It is strictly equivalent to the specific intensity in the case of a static disordered medium. This RTE-type equation states that light propagates following a random walk process where the average step-size is s and the angular redistribution at each scattering event is given by p s . A time decorrelation occurs also at each scattering event because of the Doppler shift taken into account in h in Eq. (12). This model is valid in a dilute inhomogeneous medium such as k 0 s 1 and for all possible optical thicknesses b = L/ s (from single-scattering to diffusion). To the best of our knowledge, analytical solutions of a RTE-type equation are limited to very simple geometries and the case of a spheroid, even with a homogeneous density of scatterers, requires a numerical treatment. The simplest numerical algorithm to solve this problem is the well-known Monte Carlo method. This is particularly adapted to obtain directly the field temporal correlation function for the beam radiated in far-field. To reproduce experimental conditions, parameters of the simulation were fixed at their experimental values. The phase function was taken from Eq. (8) with the subsequent coefficients values, the wavelength was λ = 0.85µm and the scattering mean free path was s = 76µm. The dilute assumption was then largely fulfilled. The spheroid was shined with a Gaussian beam of waist w = 192µm centered in the spheroid center. Finally, the Siegert relation was used to deduce the intensity temporal correlation function assuming that we had a Gaussian statistical field [58]: Fig. 3 (d) and (e) show the comparison, on virtual samples, between C obtained by simulation of the full multiple scattering model (dash black line) and C calculated with Eqs. (4) and 3 (plain orange line). In silico samples were homogeneous spheres of L = 250 µm (d) or L = 470 µm (e), containing scatterers whose displacement distribution follows Eq. (11) with σ τ = 5 µm. As an indication, the C that would have been obtained in the single scattering regime is represented (plain blue line). In all cases, we found a very good agreement between our model prediction and the simulations. This confirms the validity of the hypotheses we used, in particular the approximation of straight line trajectories in the calculation of α k .  11), using a unique function σ(τ) =vτ/(2 2/π), wherev is the average speed of the cells, v = 11.4 µm/h. In every measurement performed on spheroids, ballistic motion was observed for delays smaller than 30 min. However, a progressive deviation toward a Brownian motion (σ(τ) ∝ τ 1/2 ) was observed for longer delays, suggesting cells perform a random walk with steps lasting more than 30 min. Finally, the shape of the speed distribution gives information about the cell population diversity, with a standard deviation of 12.4 µm/h.

Cell motility under pressure
There are evidences that a moderate compressive stress (500 Pa to 5000 Pa) inhibits cell proliferation inside MCS [51]. With MS-DLS we can now easily determine whether the pressureinduced proliferation arrest correlates with a reduction in cell motility. To do so, we measured the mean cell velocity inside spheroids submitted to pressures up to 25kPa, which are the typical stresses that tumors receive in vivo from its surrounding environment [59]. Applying a pressure to spheroids did not affect the type of motion nor the shape of the distribution. For this reason, the impact of pressure on cells motility was only characterized by the evolution of the average speed of cellsv. Blue squares in Fig. 4 (d) display the mean cell speedv(P) as a function of pressure, for spheroids with a diameter of 250µm. The velocity exponentially decreases (blue lines) fromv(0kPa) = 22 ± 2 µmm/h tov(25kPa) = 7 ± 3 µm/h. Interestingly, cells inside larger spheroids (L = 470 µm) are systematically slower, even without pressure, and the plateau is reached at lower pressures (red triangles). In large spheroids, we measuredv(0) = 9 ± 1 µm/h andv(∞) = 4.1 ± 0.3 µm/h.

Radial profile of cell speed inside spheroids
We made the approximation of a homogeneous sample (hypothesis iv), but spheroids are known to be heterogeneous, with a radial structure [60]. To test whether a heterogeneity also existed in term of cell speed, and measure if it does, we probed the local cell velocity. To do so, the laser waist was reduced to 120 ± 40 µm: smaller than the spheroid but large enough to have a reasonable statistic (∼500 cells). The center of the laser was positioned at different distances from the center of the spheroid to measure cells speed radial profile [see Fig. 5 (a)]. To compare spheroids of different diameters (240 µm to 480 µm), we introduced the normalized coordinates r/R, where R is the spheroid radius and r the distance from the center, and v(r/R)/v 0 , where v 0 is the mean velocity of the cells measured at r/R = 0. When the laser is off-center, the rotational symmetry of the spheroid is lost. This is clearly visible in the speckle pattern of the scattered light [see Fig. 5 (b)]. As a consequence, the correlation function was measured for each axis. For a given direction and a given q norm, pixels used for C calculation are comprised inside two π/8 sections centered on the axis [see Fig. 5 (b)]. Radial and tangential C were then fitted together as described above, with an empirical parametric function for the displacement distribution D(∆r, τ). As previously, we chose a Gaussian distribution, but with two independent axes: where ∆r r , ∆r t and ∆r z are respectively the radial, the tangential and the optical axis components of ∆r. The average speed along the optical axis was chosen to be equal to the one tangential one, but it actually has no impact on the result of the calculation. Radial (v r ) and tangential (v t ) velocities are fitted for each laser position with respect to the spheroid center [see Fig. 5 (c)]. At the center of the spheroid, both directions are equivalent and therefore speeds are identical (v t /v r = 0.98 ± 0.04). Moving the laser toward the surface, two features appear: both speeds increase toward the surface and tangential speed becomes greater than radial speed. Note that, due to the finite beam waist, r/R may be greater than 1, while the laser beam still intersects the spheroid.

DLS extension to multiple scattering
In the diffusion regime (DWS), the large number of scattering events effectively mixes up the signal from different scattering angles, with a loss of information. As a result, the autocorrelation function is flattened, as illustrated in simulations showed in Fig. 3 (d) and (e). If we consider C as a function of the delay for a given scattering vector, the flattening implies a faster loss in intensity auto-correlation as compared to the single-scattering regime, meaning intensity fluctuates more rapidly. As expected, this flattening is exacerbated in thicker samples [ Fig. 3 (e) as compared to (d)]. Using MS-DLS, this effect can be corrected, by adjusting the parametric model of the MS signal to experimental measurements. An analytical solution allows for tremendously faster curve fitting as compared to Monte Carlo simulations. Moreover, our model helps understanding how MS modifies the signal, making it easier to choose a parametric model for the dynamics of the sample. Two aspects limits MS-DLS to samples thinner than * . First, because it relies on angular dependence of the correlation function, which is lost for L> * . Secondly, in the calculation of α k (Eq. (6)), we assumed that photons follow straight line trajectories. This is not true for L> * and proper calculation requires random walk Monte Carlo simulations.
Results from Fig. 5 show that the hypothesis of homogeneity was not valid, which might induce a bias in the average of cells speed. This assumption was used to transit from Eq. (39) to Eq. (41) in the appendix, but a more restrictive hypothesis is enough: that the dynamics of the population of cells involved in k-times scattered photons does not depend on k. This is not true for spheroids because photons scattered a small amount of time mainly come from the sides of the spheroid where cells are faster. To evaluate the bias induced in spheroid-averaged speed, we fitted data in Fig. 4 considering a simpler, but exaggerated, model where fast cells (standard deviation of 2σ τ in Eq. (11)) are involved in slightly-scattered photons (k ≤ 4, 56% of scattered light) while slow cells (σ τ /2) are involved in highly scattered photons (k>4, 44% of scattered light). Instead of previous average distances 1.3, 2.2, 3.4 and 4.6 µm, we found 0.8, 1.3, 2.2 and 2.5 µm. Under homogeneity hypothesis, speed are overestimated by a factor 1.7 in average because phase dynamics is dominated by slightly scattered photons, involving fast cells in our case. In reality, this factor should be smaller as highly scattered photons also involve cells from the surface. In any case, relative measurement of speed is still valid and absolute measurement can be achieved by reducing the width of the probing laser beam.

Cell motility under pressure
Pressure inside multi-cellular system result from two components: a growth induced internal pressure [61,62] and an external pressure coming from the surrounding tissues. Both are present in our experiments, the second one being mimicked by the pressure induced by Dextran molecules. The resulting pressure increases from the periphery to the core [63] which correlates with the decrease of speed observed in Fig. 5 (c). To check whether pressure causes a cell displacements reduction, we varied its intensity by adding an external pressure or by making larger (older) spheroids with an increased growth pressure [64]. Figure 4 (d) shows that in both cases, the average cell speed is decreased when pressure is increased. The speed is more than halved at 7kPa and 15kPa, respectively for large and small spheroids ; Then a plateau is reached. It is noticeable that the same pressure slows cell proliferation in spheroids made from cells of the same line [50]. Interestingly, this pressure also corresponds to the maximal stress that a growing tumor can exerts on its surroundings [59], defined as the homeostatic pressure in [65]. We conclude that aggregates are sensitive to stresses up to the homeostatic pressure, with no additional effect above this threshold. The effect of pressure on cell displacements depends on the cell line. Janet et al. [66] showed on cells sheets that pressure could increase the migration of aggressive cancerous cells while decreasing the migration of less invasive cells. It might also depend on the direct surroundings of a cell. Using the same cell line as ours, Alessandri et al. [67] measured an increased migration under pressure at the surface of encapsulated spheroids.

Conclusion
MS-DLS is an analytical method that can measure cell motion in opaque samples thinner than the transport mean free path * . Based on small-angle light scattering, it can be implemented on any microscope devoted to biological imaging with the addition of few optical elements. We applied it to the measurement of cell motion inside multicellular aggregates, quantifying how much external mechanical constraints impair cell displacements. This experiment could be extended to the study of cell migration regulation under numerous conditions: mechanical constraints, genetic mutations, temperature, drugs. More generally, our technique would prove useful in studying the migration of specific cells (cancer cells, macrophages) inside complex environments (gels with inclusions, ex vivo tissues).

Appendix : MS-DLS analysis formula demonstration
In this appendix, we demonstrate the formula used to obtain cell displacement distribution from the intensity correlation function (Eqs. (3) and (4) in Results section).

Light propagation
Light propagation is modeled as photons going from a scatterer to another inside the spheroid to finally emit a wave at the last scattering. First, we calculate the relative phase shift between photons, depending on the changes of propagation direction at each scattering event. Figure 6 shows an example of a photon being scattered twice: at positions r 1 and r 2 . As phase reference in the calculation of the phase shift, we take a photon scattered once, at the center of the spheroid. u i , u 1 and u s are the directions of the photon incoming, scattered once and outcoming respectively. The scattering angle of observation is q = k 0 (u s − u i ) and we note q 1 = k 0 (u 1 − u i ) and δ = r 2 − r 1 . As in DLS or DWS, the phase shift is calculated from the difference in the optical path of the photon. In this example, the phase shift φ is: = −q 1 .r 1 + k 0 (u 1 − u s ).r 2 = −q 1 .r 1 + (q 1 − q).r 2 (17) This result is equivalent to what we would obtain from the transport equation (Eq. (12)) which was derived from first principles by controlling rigorously the approximations made [57].
In the general case of k scattering events at positions r m for m in [1.
.k], with directions u m associated to q m = k 0 (u m − u i ), the phase shift is: With q 0 = 0 and q k = q. The electrical field on the detector is taken as the addition of all waves emitted by N scattering cells: where E j is the amplitude of the field scattered by the j th cell and φ j its phase shift. E j is determined by the cell shape. We assume that light fluctuations due to cell shape changes are negligible as compared to those due to cell displacements (hypothesis i). The smaller the scattering vector is, the less details of the cell shape impact E j (q) and the more valid this hypothesis is. Following hypothesis i, we assume E j to be constant over time.

Phase function
Let us start with the phase function in Eq. (4). The temporal average intensity is: We make the hypothesis that cell positions are not correlated (hypothesis ii). Because different phases φ rely on different cells, they are not correlated and we have: (20), we have: As described in the main section, cells considered here correspond to the last event of scattering in the scattering chain of photons. We group them depending on the number of scattering in the chain.
Each group includes an average of Nα k cells (see the definition of α k in the Results section DLS extension to multiple scattering), so that: where the average is over cells of indexes j(k) that are last scatterers of k times scattered photons. E 2 j can be decomposed into an amplitude factor I 0 and a phase function. The phase function is the one of a single scatterer, noted p s in the main section, with an orientation determined by the direction before the last scattering: The average over j(k) is then taken as an average over possible trajectories. For k = 1, q k−1 = 0 and: For k = 2, the photon is first scattered in a direction q 1 with a probability p s (q 1 ) and the average calculated as: where [p s p s ] is the convolution product of p s with itself. In the general case, the photon undergoes k-1 scattering in the successive directions q 1 ,. . . , q k−1 , associated with a probability p s (q 1 )p s (q 2 − q 1 ) · · · p s (q k−1 − q k−2 ) and the final phase function is p s (q − q k−1 ). Following the same reasoning, we obtain: where [p s k p s ] is p s convoluted to itself k-times, with the convention [p s 0 p s ] = p s . Replacing in Eq. (22), we get: To obtain p, Eq. (27) needs to be divided by the total intensity scattered: Finally, we have:

Intensity correlation function
This section details the calculation leading to Eq. (3) on C.
For our sample, I(q, t) ϕ = I(q, t) t . The average intensity is thus written I(q) hereafter and is time and ϕ independent. Under this condition, the intensity correlation function C becomes: C(q, τ) = I(q, t)I(q, t + τ) t I(q) 2 − 1 We introduce cell displacements over the period τ as: ∆r(t, τ) = r(t + τ) − r(t) (31) Phase modifications are then expressed as: ∆φ(q, t, τ) = φ(q, t + τ) − φ(q, t) Where we made the approximation that q m are constant over time (hypothesis iii). This approximation is valid if the distance between scattering events of a single photon trajectory (≈ mean free path) is large as compared to cell displacements. First, we calculate I(q, t)I(q, t + τ) t using Eq. (19): E j (q)E l (q)E m (q)E p (q)e i(φ j (q,t)−φ l (q,t)) e i(φ m (q,t+τ)−φ p (q,t+τ)) t = N j,l,m,p E j (q)E l (q)E m (q)E p (q) e i(φ j (q,t) e −iφ l (q,t) e iφ m (q,t+τ) e −iφ p (q,t+τ) t (33) As in the calculation of the phase function, A = e i(φ j (q,t)−φ l (q,t)) e i(φ m (q,t+τ)−φ p (q,t+τ)) t is non null only if: • j = l and m = p, leading to A = 1.
• j = p and l = m with m p, leading to A = e i(φ j (q,t)−φ j (q,t+τ)) t e i(φ m (q,t+τ)−φ m (q,t)) t = e −i∆φ j (q,t,τ) t e i∆φ m (q,t,τ) t • j = m and l = p with m p, leading to A = e i(φ j (q,t)+φ j (q,t+τ)) t e −i(φ p (q,t)+φ p (q,t+τ)) t = 0  (36) and the correlation function becomes: Now, we take separately the sum, and rearrange the terms depending on the number k of scattering events, as we did previously in Eq. (22): Nα k E 2 j (q)e i∆φ j (q,t,τ) t,j(k) (38) where the average is now also over cells of indexes j(k) that are last scatterers of k times scattered photons. Finally, using Eqs. (23) and (32) we calculate: We assume that all cells follow the same displacement probability D τ (∆r), i.e. the sample is homogeneous in terms of cell displacements (hypothesis iv). As a consequence, photon trajectories, given by the q j , are independent on specific cells and the average over trajectories can be separated from the average over time.