STRESS, an automated geometrical characterization of deformable particles for in vivo measurements of cell and tissue mechanical stresses

From cellular mechanotransduction to the formation of embryonic tissues and organs, mechanics has been shown to play an important role in the control of cell behavior and embryonic development. Most of our existing knowledge of how mechanics affects cell behavior comes from in vitro studies, mainly because measuring cell and tissue mechanics in 3D multicellular systems, and especially in vivo, remains challenging. Oil microdroplet sensors, and more recently gel microbeads, use surface deformations to directly quantify mechanical stresses within developing tissues, in vivo and in situ, as well as in 3D in vitro systems like organoids or multicellular spheroids. However, an automated analysis software able to quantify the spatiotemporal evolution of stresses and their characteristics from particle deformations is lacking. Here we develop STRESS (Surface Topography Reconstruction for Evaluation of Spatiotemporal Stresses), an analysis software to quantify the geometry of deformable particles of spherical topology, such as microdroplets or gel microbeads, that enables the automatic quantification of the temporal evolution of stresses in the system and the spatiotemporal features of stress inhomogeneities in the tissue. As a test case, we apply these new code to measure the temporal evolution of mechanical stresses using oil microdroplets in developing zebrafish tissues. Starting from a 3D timelapse of a droplet, the software automatically calculates the statistics of local anisotropic stresses, decouples the deformation modes associated with tissue- and cell-scale stresses, obtains their spatial features on the droplet surface and analyzes their spatiotemporal variations using spatial and temporal stress autocorrelations. The automated nature of the analysis will help users obtain quantitative information about mechanical stresses in a wide range of 3D multicellular systems, from developing embryos or tissue explants to organoids. Author summary The measurement of mechanical stresses in 3D multicellular systems, such as living tissues, has been very challenging because of a lack in technologies for this purpose. Novel microdroplet techniques enable direct, quantitative in situ measurements of mechanical stresses in these systems. However, computational tools to obtain mechanical stresses from 3D images of microdroplets in an automated and accurate manner are lacking. Here we develop STRESS, an automated analysis software to analyze the spatiotemporal characteristics of mechanical stresses from microdroplet deformations in a wide range of systems, from living embryonic tissues and tissue explants to organoids and multicellular spheroids.

Mechanical forces are known to play an essential role in the control of cell behavior as 2 well as in tissue and organ morphogenesis. During embryogenesis, cells apply differential 3 forces to guide tissue flows and build embryonic structures [1,2]. While mechanical 4 forces are generated at the cell and subcellular scales, it is the collective force generation 5 by many cells and its transmission over supracellular length scales that shapes 6 functional tissues [3]. Understanding the characteristics of mechanical stresses at 7 different length-and time-scales is necessary to bridge the gap between force generation 8 at the cell scale and organ formation at much larger scales [1][2][3]. Beyond their 9 important role in morphogenesis, mechanical forces are also known to play an important 10 role in the control of cell behavior. Mechanobiology studies have shown that important 11 cell behaviors, such as cell differentiation or proliferation, are affected by the mechanical 12 forces acting on them [4,5]. Therefore, having the ability to quantify the spatiotemporal 13 characteristics of mechanical stresses in 3D multicellular systems can help advance our 14 understanding of how mechanics affects both individual cell behavior and their 15 collective dynamics. 16 While several measurement techniques exist to quantify forces in 2D cell monolayers 17 on synthetic culture dishes [6], techniques to measure mechanical forces in 3D 18 multicellular environments are scarce [7,8]. The development of microdroplet techniques 19 enabled the quantification of forces in 3D multicellular environments, from 3D cell 20 culture to developing tissues, in vivo and in situ [9][10][11][12]. Mirroring oil microdroplets, gel 21 microbeads have been recently used to measure stresses in multicellular systems [13][14][15][16] 22 and also for single cell studies [17]. Both oil microdroplets and gel microbeads enable 23 measurements of mechanical stresses in various 3D systems, allowing quantitative 24 mechanobiology studies in a more physiological context. 25 The accurate measurement of mechanical stresses with deformable particles relies on 26 precise measurements of surface deformations of the probe used (oil microdroplets or gel 27 microbeads). Oil microdroplets enable precise measurements of anisotropic stresses 28 solely from knowledge of their surface geometry and the droplet interfacial tension 29 (which can be calibrated in vivo and in situ) [9,12], whereas measurements of stresses 30 with gel microbeads require the quantification of strain fields inside the probe [15]. 31 However, approximate methods have been recently developed to obtain mechanical 32 stresses using gel microbeads from their surface deformations and the bead stiffness [17]. 33 While different analysis methods exist to reconstruct the geometry of the probes in 34 3D [12,17,18], an automated and and reliable tool to obtain stresses in 3D and time, as 35 well as to analyze their spatiotemporal characteristics, is lacking. the test case used here, this is an oil microdroplet inserted in developing zebrafish 44 embryo) and automatically reconstructs the surface deformations of the particle, obtains 45 its geometrical characteristics and analyzes the spatiotemporal changes of mechanical 46 stresses. STRESS automatically decouples mechanical stresses occurring at 47 supracellular (tissue) scales and cell-scale by analyzing different deformation modes as 48 well as maximal and minimal values of stress on the droplet surface. It also obtains the 49 characteristics of spatial inhomogeneities of the stresses around the droplet, providing 50 insight on the spatial structure of the stresses in the tissue. Moreover, it automatically 51 calculates the persistence of total stresses, as well as of cell-and tissue-scale stresses, 52 from temporal autocorrelation functions. Finally, the analysis also provides accurate Imaging and sample preparation of surface-labeled microdroplets in tooth mesenchymal 60 cell aggregates and volume-labeled (or interior-labeled) droplets in developing zebrafish 61 tissues was done as previously described in references [12] and [9], respectively.

62
Surface-labeled droplets were coated with Cy5-streptavidin. Volume-labeled droplets 63 were prepared by dissolving FCy5 dye [19] in Novec 7300 fluorocarbon oil (3M) at a 25 64 µM concentration. A fluorinated Krytox-PEG(600) surfactant (008-fluorosurfactant, 65 RAN Biotechnologies) was also diluted in the fluorocarbon oil at a 2.5% (w/w) 66 concentration. This fluorescently-labeled oil was injected directly in the embryo, as 67 previously described [9], and the volume of the droplet was set by the injection pressure 68 and the duration of the injection.

69
Obtaining the point cloud 70 A point cloud representation of the droplet surface is constructed by tracing rays (or 71 lines) through the droplet surface in the fluorescence 3D image of the droplet and fitting 72 the sampled intensity profiles with an appropriate edge profile model, which depends on 73 what feature of the droplet is labelled fluorescently. Ray tracing paths were determined 74 as previously described [18]. For droplets labeled at their surface ( Fig. 1a), we employ a 75 Gaussian edge profile model [18] (Fig. 1d). However, for droplets prepared with 76 fluorescently labelled fluorocarbon oil [9,19,20] (Fig. 1b), which display fluorescence 77 signal throughout the droplet, we used an attenuated sigmoid profile model that 78 captures both the fluorescence intensity fast change at the droplet surface and the mild 79 attenuation of fluorescence inside the droplet (due to imaging a 3D object). Specifically, 80 we use the fit model F (x) = (d + c * (x − b))/(1 + exp(a * (x − b))) + e for the intensity 81 profile along the traced rays within the image volume at distance x from a starting Once a point cloud representing the droplet surface has been obtained, we use a 89 spherical harmonic basis to obtain a global approximation of the surface. A global 90 function approximating the object surface is constructed by finding the least-squares fit 91 of spherical harmonics to the point cloud [21], up to a given degree. The degree of the 92 fit is limited by the number of points in the point cloud, since we need at least one 93 degree of freedom in our data for each degree of freedom in our fit. However, in order to 94 have the fitting functions properly constrained, we restrict ourselves to two degrees of 95 freedom (data points) in our point cloud for each degree of freedom in our fit. This 96 means that we would ideally fit no more than 50 global basis functions to a point cloud 97 with 100 points. The spherical harmonic fits provide a smooth representation of the 98 surface with spherical topology, whose variations are limited by the maximum degree of 99 harmonics chosen. Importantly, the spherical harmonics converge spectrally to smooth 100 functions on the sphere [21], so it is possible to achieve a faithful representation of a 101 smooth droplet surface with a small to moderate number of modes. This approach, 102 adapted from recent work on surface flows on thin-shell membranes [22][23][24], has also 103 been used recently to analyze cellular shapes [25].

104
Mean curvature map 105 Once we have a smooth global representation of the droplet surface, we calculate its 106 mean curvature, which depends on the coordinate system chosen to parameterize our 107 global fit. In the simplest case, we define the surface points in a cartesian coordinate 108 system, namely r s = (x s (θ, φ), y s (θ, φ), z s (θ, φ)), where the coordinates of each point 109 are given in terms of the azimuthal and polar coordinates, θ and φ respectively, of our 110 spherical harmonic basis. By taking the first and second derivatives of our spherical 111 harmonic basis functions, we can calculate any geometric quantity associated to the 112 droplet surface, including the mean curvature H and the Gaussian curvature K, using 113 standard differential geometry [26].

114
All of the calculations are done using python's numpy library [27]. We validated the 115 convergence of the representations of these geometric quantities using the method of 116 manufactured solutions [23,24]. The symbolic expressions needed for these validation 117 tests, performed on various surface geometries, were generated by python's sympy 118 package [28]. 119 Surface integrals and particle volume 120 One major advantage of using spherical harmonics to obtain a global surface 121 representation is that, for a given degree, there is a set of lebedev quadrature points 122 that integrates these harmonics exactly on the surface of the unit sphere [29]. For our 123 calculations, we use a quadrature of up to 5810 lebedev points, which can integrate 124 inner-products of up to degree 65 spherical harmonics exactly. We specifically use these 125 lebedev points to represent the surface in our global fit and also to display the surfaces. 126 In order to calculate the object's (microdroplet or gel microparticle) volume V , we 127 adapt our surface quadrature. Parameterizing the droplet surface in spherical 128 coordinates (r, θ, φ) the object's volume is given by where R(θ, φ) is the droplet radius at each point of the surface and S 2 denotes the 130 surface of the unit-sphere. We tested the accuracy of volume measurements by 131 calculating the volume of a very eccentric ellipsoid (with semi-axes lengths 3 µm, 5 µm, 132 March 23, 2021 4/20 and 7 µm) using Eq. 1 and comparing it to its analytically known volume of 140 πµm 3 . 133 Using only 590 quadrature points (much less than the 5810 quadrature points used for 134 all our calculations), we find the relative error in the ellipsoid volume to be 10 −5 %.

135
Tissue-scale stresses 136 To obtain the ellipsoidal mode of the droplet shape, we perform a least-squares 137 ellipsoidal fit on the segmented point cloud. Specifically, we fit the equation: 138 Ax 2 + Bxy + Cy 2 + Dxz + Eyz + F z 2 + Gx + Hy + Iz = 1 to the point cloud and to 139 find the 9 coefficients that represent the least squares ellipsoid. To obtain the mean 140 curvature of this ellipsoid, we transform this ellipsoid into the coordinates system 141 (x 1 , x 2 , x 3 ) defined by its principal directions. These coordinates are along the major, 142 medial, and minor axes of the ellipsoid, respectively, with the same center, namely 143 x 1 = a cos(u) sin(v)  where γ is the droplet interfacial tension.

153
Cell-scale stresses 154 To quantify the stresses arising from higher order deformation modes (beyond the 155 ellipsoidal mode), which are associated with length scales closer to the cell size, we 156 calculate the deviations of the droplet deformations from the ellipsoidal mode. For any 157 given timepoint, we subtract the calculated total anisotropic stress value σ A from the 158 tissue-scale stresses associated with the ellipsoidal deformation, namely σ A T , to obtain 159 the cell-scale stresses (Eq. 6). Specifically, to compute the cell-scale stress anisotropy 160 between points p and q on the surface, we compute )], where we choose the points p and q on the 162 ellipsoid that correspond with the same ellipsoidal coordinates as p and q on the droplet. 163 In order to obtain the extrema on the deformed surface, we use the surface 164 triangulation of the 5810 lebedev nodes on the surface to relate them to each other.

165
Two points are classified as neighbors if they are connected by an edge. Local maxima 166 are points whose value of cellular stress is greater or equal to each of their neighbors' 167 values, and local minima are points whose value of cellular stress is less than or equal to 168 each of their neighbors' values. The extrema stresses (or adjacent extrema stresses) are 169 obtained by calculating the difference in cell-scale stresses between maxima and minima 170 (or adjacent maxima and minima).

Geodesic distances 172
To calculate distances between points on the curved droplet surface, we first triangulate 173 the surface using a Delauney Triangulation of the lebedev quadrature points on the 174 sphere using python's scipy package [31]. The geodesic distance between lebedev 175 quadrature points on the droplet surface is calculated from the surface triangulation 176 graph using using the gdist package [32]. To obtain maximal resolution and minimize 177 noise, we use the full set of 5810 lebedev nodes in our calculations.
where k indicates that this calculation can be performed for any stress field on the 186 droplet surface, including the total stress σ A and the tissue-and cell-scale stresses σ A T 187 and σ A C , respectively. We have checked that changing the value of does not change our 188 results. Very small values of (less than 1% of the geodesic diameter of the droplet) 189 lead to more noisy results for the autocorrelation, but do not change the features of the 190 obtained autocorrelation function.

191
For temporal correlations, we use a radial parameterization of the droplet (a system 192 of spherical coordinates with the same center as the droplet). This allows us to relate 193 points on the droplet surface with the same spherical coordinates as the droplet changes 194 over time. We then calculate the temporal correlations of the stress values at the same 195 spatial location (same angular coordinates) and average the results for all the points on 196 the droplet surface (Eq. 8).

198
To explain and exemplify the STRESS analysis pipeline, we focus on the use of In order to measure stresses from deformable particles and, in particular, from oil 207 microdroplets, it is necessary to image the particles in 3D. While it is also possible to 208 measure stresses from partial 3D microdroplet reconstructions and even from 2D 209 confocal sections [9,12], the analysis presented herein focuses on complete particle . Before imaging, the droplets were inserted in a 3D aggregate of tooth mesenchymal cells (A; cytoplasmic label, yellow) and in the presomitic mesoderm of developing zebrafish embryos (B; cell membranes, yellow). C, Sketch of a deformed particle labelled fluorescently at its surface. The inset shows a detail of the particle surface with its local normal vectorê n and the detected surface points (point cloud; red). D-E, Fluorescence intensity along rays traced along the normal direction to the surface for both surface labelled (D) and internally labelled droplets (E). Fits to the intensity profile along a ray (magenta) are used to detect the location of the surface (red dashed line). F, A local representation of the surface geometry can be done using Monge patches to estimate the local mean curvature at every point on the droplet surface. G, A global representation of the surface geometry can be obtained using a spherical harmonic basis, which progressively capture finer details as the number of modes is increased. The contribution of each mode depends on every point on the surface. The mean curvature of the deformed particle can be obtained at every point of the surface with high spatial resolution.

218
Previous raytracing algorithms have been shown to provide a faithful representation 219 of the particle's surface for surface-labeled deformable particles [18]. This method 220 identifies first the center of the droplet and traces rays in all spatial directions (Fig. 1C). 221 By measuring the intensity profile along the rays and detecting the intensity maximum, 222 it is possible to detect the location of the surface along the ray with high spatial 223 resolution (Fig. 1D). Repeating this procedure for all the rays provides a point cloud 224 that represents the surface of the object. A refined point cloud is then obtained by 225 repeating the procedure for rays retraced along normal directions to the surface [18] 226 (Fig. 1C). We have extended this algorithm for particles fluorescently labeled in their 227 interior (Fig. 1B). To do so, we employ the same raytracing procedure, but implement a 228 different fitting function for the intensity profile along the traced rays because particles 229 with interior label display a high fluorescence signal in their interior and a sharp 230 decrease at their surface (Fig. 1B,E; Methods). Fitting the intensity profile along a ray 231 allows the detection of the surface location on that same ray ( Fig. 1E; Methods). By 232 repeating this procedure along all traced rays, we obtain the point cloud that represents 233 March 23, 2021 7/20 the location of the particle's surface.

234
Once the point cloud that represents the deformed particle's surface has been 235 determined, there are several approaches to characterize its geometry. We previously 236 developed a methodology to obtain the surface mean curvature map from local 237 quadratic fits [18] (Fig. 1F). This strategy relies on the construction of local Monge 238 patches of adaptable size and the approximation of the local surface geometry within 239 the patch. Because this method employs only on local information, it allows the 240 calculation of curvature maps on both fully and partially reconstructed particles, but is 241 more prone to errors associated with limited spatial resolution or poor signal-to-noise 242 ratios, especially for very small patch sizes. To overcome these issues in full 3D particle 243 reconstructions, we mathematically characterize the surface using global fits of 244 controlled degree (Fig. 1G). Specifically, we use least-squares fits of spherical harmonics 245 (up to a desired degree) to the segmented point cloud (Methods). The obtained surface 246 representation using global fits provides a smooth, highly-resolved and accurate the stresses in the tissue. We first define the total anisotropic stresses, σ A ( x s , t i ), which, 266 for oil microdroplets, are obtained directly from the mean curvature map as previously 267 described [12], namely where γ is the droplet interfacial tension and H 0 is the average of the mean curvature Methods). Knowing the droplet interfacial tension, we use Eq. 2 to obtain the 3D map 276 of total anisotropic stresses on the droplet surface at each timepoint ( Fig. 2A,B), as the 277 droplet follows morphogenetic flows in the tissue. Since the global surface  (Fig. 2C). For this reason, all our calculations obtain the value of the average 283 mean curvature from its surface integral.

284
In order to evaluate the accuracy of the global surface representation, we make use 285 of the Gauss-Bonnet Theorem, which states that the integral of the Gaussian curvature 286 K over any smooth surface with spherical topology should equal 4π. We integrate the 287 Gaussian curvature K over the droplet surface and measure the relative deviation (or 288 error), , from the expected 4π, namely ≡ |1 − M KdA/4π|, at each timepoint 289 (Fig. 2D). Inaccurate droplet segmentations (i.e., those containing errant points or large 290 holes in the point cloud) will cause the error to increase as the degree of the surface 291 representation increases, allowing the automated detection of incorrectly segmented 292 images, which we do not include in our measurements (Methods). For the example 293 shown here, the relative errors are very small (Fig. 2D), indicating no issues with 294 droplet segmentation, and fluctuate over time as the droplet segmentation is slightly 295 different at each timepoint.

296
In addition to surface integrals, it is also possible to accurately calculate volume  (Fig. 2D). In contrast to microdroplets, changes in 307 the volume of gel microbeads provide a measure of the local isotropic stress or tissue 308 pressure. Therefore, the ability to accurately measure volume over time can be 309 employed to track temporal changes in tissue pressure when using gel microbeads.

310
While the full stress map on the droplet surface provides spatiotemporal information 311 about stresses ( Fig. 2A,B), in most applications it is necessary to define a single 312 numerical value of the total anisotropic stress for a droplet at a given timepoint, so that 313 it is possible to perform statistics and monitor the temporal evolution of stresses in the 314 tissue. To do so, we first calculate the statistical distribution of the magnitude of total 315 anisotropic stresses from all the values of the stresses at different points on the droplet 316 surface (probability density function; Fig. 2F). This distribution has important 317 information about the variations in the magnitude of local stresses in the tissue. We . If α = 0, the minimal and maximal values 324 correspond to the actual minimal and maximal values of the total anisotropic stress on 325 the droplet surface. However, since extreme values are prone to noise, it is preferable to 326 define a small value of α to cut off extreme fluctuations. Here we used α = 0.05 (or 5%; 327 removing the 5% largest and 5% smallest values of the total stress). By measuring the 328 magnitude of total anisotropic stresses at each timepoint, we can track the time 329 evolution of total anisotropic stresses ∆σ A (t) (Fig. 2H).  While high order deformation modes (localized higher/lower curvature changes) provide 335 information about specific length scales, the lowest order deformation mode, namely the 336 ellipsoidal mode, provides information about the mechanical stress anisotropy occurring 337 at the length scale of the droplet size (Fig. 3A,B). By making droplets bigger than the 338 cell size (typically ∼ 4 cell diameters; Fig. 3B), the ellipsoidal deformation mode enables 339 the measurement of the local value of the anisotropic stress at supracellular scales 340 (tissue-scale), σ A T . This tissue-scale stresses essentially average out smaller deformations 341 on the droplet surface and are thus associated with stress anisotropy in the tissue 342 propagating at length scales larger than the cell size.   Fig. 3D), reveal the principal directions of stress anisotropy in the tissue. By 347 construction, we choose a to be the longest semiaxis, andê 1 thus reveals the direction of 348 droplet elongation. The tissue-scale stress anisotropy between two principal directions i 349 and j (i, j = 1, 2, 3 and i = j) is obtained by calculating the stress anisotropy between 350 these directions, namely where H e corresponds to the mean curvature of the ellipsoid and x i and x j are the 352 surface locations where the principal axisê i andê j cross the ellipsoid surface (Fig. 3B). 353 In order to obtain a single measure of the tissue-scale stress anisotropy, we calculate the 354 difference between the maximal and minimal stresses on the ellipsoid (Methods), which 355 are associated with two principal directions, namely

378
Beyond the tissue-scale stresses associated with the ellipsoidal droplet deformation, 379 it is also possible to define the total stress anisotropy along the principal directions, which differ from the the tissue-scale stresses defined above (Eq. 3) in that it employs 382 the mean curvature H of the actual droplet surface at x i and x j rather than the mean 383 curvature of the ellipsoidal mode at these points. While the total stress anisotropy 384 along principal directions σ A T T provides the actual value of stress anisotropy, it is more 385 noisy than σ A T as it is affected by higher order deformations (Fig. 3H). In contrast, since 386 σ A T (Eq. 4) uses the mean curvature of the fitted ellipsoid, it removes the noise 387 associated with higher order deformation modes and provides a better defined measure 388 of tissue-scale stresses.

389
Cell-scale stresses 390 To quantify cell-scale stresses we study the deviations of the droplet deformations from 391 the ellipsoidal mode. We first express the droplet deformations on the ellipsoidal 392 reference frame and then calculate the deviations of the deformations from the 393 ellipsoidal mode (Fig. 4A), which, by definition, only contain deformation modes of 394 higher order than the ellipsoidal one. To calculate the deviations from the ellipsoidal 395 mode, we perform a least squares fit of the distances between the droplet surface and 396 the ellipsoid, using the coordinate system of the fitted ellipsoid to parameterize the 397 residuals (or deviations). For any timepoint t i , cell-scale stresses are given by where H e0 is the average mean curvature of the ellipsoidal mode. Since H 0 H e0 399 (Fig. 2c), the term H 0 − H e0 vanishes in Eq. 6.

400
As for the total anisotropic stresses and the tissue-scale stresses, it is important to 401 define a single numerical measure to quantify the time evolution of the magnitude of 402 cell-scale stresses. This can be done in an analogous way as done for the total 403 anisotropic stresses above, by obtaining the normalized distribution of the magnitudes 404 of cell-scale stresses (Fig. 4B), calculating its cumulative distribution (Fig. 4C) and 405 defining maximal and minimal stresses associated with a percentage cut-off, σ A C,M and 406 σ A C,N , respectively. The time evolution of the magnitude of anisotropic cell-scale stresses 407 is then given by ∆σ A C (t) ≡ σ A C,M (t) − σ A C,N (t) and can be monitored over time (Fig. 4D). 408 The magnitude of cell-scale stresses is larger than that of tissue-scale stresses and is very 409 similar to the measured magnitude of total anisotropic stresses (Fig. 2h), indicating the 410 cell-scale stresses dominate the total magnitude of the stresses.

411
To better understand the structure of deformation modes higher than the ellipsoidal 412 mode, we decompose the map of cell-scale stresses into spherical harmonics up to degree 413 20. We use the fact that the coefficients of spherical harmonics show a symmetric 414 structure when representing real quantities and combine the contributions of 415 off-diagonal modes of a given order with their opposite to evaluate the contribution of 416 each mode (Fig. 4E). These combined amplitudes (or coefficients) provide a measure of 417 the relative contribution of each mode to the total droplet deformation. We find that 418 some small number of modes contribute more to the droplet deformations than others 419 (Fig. 4F), but it requires about 50 modes to capture 90% of the droplet deformations 420 (Fig. 4G). The most represented modes (orders 3 and 4) are associated with periodic 421 deformations of the droplet surface at length scales of approximately 12 to 17µm, since 422 the droplet radius is approximately 16µm (Fig. 4H), suggesting the existence of a 423 characteristic length scale of droplet deformations. Beyond the specific example studied 424 here, it is possible to understand the contribution of stresses at different length scales by 425 using mode decomposition. extrema stresses σ A X , and obtain their distribution (Fig. 5E), which provides information 443 about the spread in the magnitude of anisotropic stresses between extrema on the 444 droplet surface. Our results show a well-defined peak at 80 Pa, but with a distribution 445 characterized by a long tail (average 230 Pa; median 150 Pa). Because all pairs of 446 maxima and minima were considered in this calculation, it is not possible to reveal any 447 local structure associated with these extrema. To characterize the length scales and 448 stresses associated with adjacent extrema (Fig. 5B), we calculate the anisotropic stresses 449 between an extrema (maximum or minimum) and its closest neighbor opposite extrema 450 (Fig. 5D). The distribution of geodesic distances between adjacent pairs of opposed 451 extrema shows a well-defined peak at approximately 3 µm, with the average distance 452 being 4 µm (Fig. 5D), indicating the existence of a well-defined length scale between 453 adjacent extrema at the droplet surface. While the cell diameter is approximately 12 454 µm in the studied tissues [20], the average length of cell-cell contacts in the tissue is 455 approximately 5 µm [20], indicating that extrema are likely associated with the spatial 456 inhomogeneities in stresses due to cell-cell junctions contacting the droplet surface. surface. Since extrema occur at distances much smaller than the droplet size, the results 462 obtained by considering only cell-scale stresses or the total stresses are very similar.

463
Spatial and temporal autocorrelations of anisotropic 464 stresses 465 The analysis described above enables the characterization of the total, tissue-scale and 466 cell-scale anisotropic stresses and shows the existence of spatial features of the stresses 467 on the droplet. To quantify the degree to which stresses on the droplet are spatially and 468 temporally correlated, we calculate their spatial and temporal autocorrelations, as these 469 provide information about the spatial structure and persistence of stresses, respectively. 470 Spatial autocorrelations are typically used in condensed matter physics to identify 471 the existence of repetitive spatial structures [33]. In this case, the spatial 472 autocorrelation of anisotropic stresses provides information about the spatial structure 473 of stress inhomogeneities around the droplet. The normalized spatial autocorrelation 474 function is given by where s is the geodesic distance between two given points on the droplet surface, x s is 476 the coordinate of a point on the surface and k indicates that this calculation can be 477 performed for any stress field on the droplet surface (k = total, tissue-scale or cell-scale 478 stresses). Since the spatial autocorrelation is calculated at each timepoint, we obtain the 479 time evolution of the spatial correlations in the system for the total, the cell-scale and 480 the tissue-scale stress anisotropies (Fig. 6A-C). The spatial autocorrelation of total and 481 cell-scale stresses is very similar because the total anisotropic stresses are dominated in 482 magnitude by the cell-scale stresses (Fig. 6A,B). The correlation decays rapidly and 483 reaches a minimum at a length scale of about the cell-size ( 12 µm), displaying a small 484 anti-correlation at this point. In contrast, the variations in the spatial autocorrelation of 485 tissue-scale stresses reflect the geometry of the ellipsoidal mode (Fig. 6C). Beyond the spatial structure of stress inhomogeneities around the droplet, it is 487 important to characterize how long different stresses persist in the tissue. To do so, we 488 calculate the temporal autocorrelation of anisotropic stresses, which is defined as approximately a few minutes (Fig. 6B). The total stresses show a similar, albeit slightly 494 longer, persistence, in agreement with our results that cell-scale stresses dominate the 495 total stresses (Fig. 6A). In contrast, tissue-scale stresses are much more persistent, with 496 stresses displaying substantial correlations even after 30 minutes, and only change 497 significantly (and sharply) after approximately 50 minutes (Fig. 6C). These results are 498 in agreement with previous analysis done with 2D confocal sections of microdroplets [9], 499 instructions, can be found at https://github.com/campaslab/STRESS.

507
The analysis implemented in the STRESS software allows automated surface 508 reconstructions of fluorescently-labeled deformable objects (such as oil microdroplets or 509 gel microbeads) and the analysis of its geometrical features. For the specific case of oil 510 microdroplets, the geometry of the droplet contains the necessary information to

521
While using a local representation of the geometry of deformable particles is useful 522 in some cases, especially when it is not possible to reconstruct the entire droplet or a 523 coordinate-free system is more adequate, the global surface representation described 524 herein provides a smoother and more accurate representation of surface deformations 525 that is less sensitive to noise. Moreover, the Gauss-Bonnet test allows a quantification 526 of the error in the reconstruction of the particle deformations, enabling the automated 527 detection of poorly resolved measurements. Finally, the global surface representation 528 allows the calculation of surface integrals, as well as the volume of the droplet with very 529 high accuracy. We believe that accurately calculating the particle volume will be useful 530 for measurements of the isotropic pressure with gel microbeads.

531
Performing instantaneous measurements of stresses at few timepoints could 532 potentially be done by separately analyzing each timepoint using previous analysis 533 methods. However, the analysis of the time evolution of stresses at high temporal 534 resolution generates a large number of 3D particle reconstructions to analyze, 535 precluding the separated analysis of each time frame. STRESS enables the analysis of 536 the time evolution of stresses from oil microdroplets, in vivo and in situ, at high 537 temporal resolution.

538
STRESS automatically decouples the contributions of different deformation modes, 539 enabling the study of stresses associated to different length scales. In the present form, 540 STRESS decouples the ellipsoidal deformation mode, which is associated to tissue-scale 541 stresses (stresses occurring at supracellular scales), and the cell-scale stresses, associated 542 with higher order modes, especially those characterized by length scales similar to the 543 cell size. Moreover, we have developed a new method to characterize stresses associated 544 with extrema (maxima and minima) in mean curvature of the droplet surface, as these 545 extrema are linked to the spatial inhomogeneities in the structure of the tissue (e.g., 546 cellular structure).

547
Finally, to better understand any spatial features of the stress inhomogeneities in the 548 tissue around the droplet, STRESS calculates the spatial correlations of the measured 549 March 23, 2021 17/20 stresses at each timepoint. This spatial correlation quantifies the length scale over 550 which stresses are correlated on the droplet surface and can be obtained for the stresses 551 associated to any deformation mode. Moreover, the temporal evolution of the spatial 552 autocorrelation provides information about the temporal changes in the spatial 553 inhomogeneities of stresses around the droplet. Finally, for a single droplet, we calculate 554 the temporal autocorrelation of total stresses, as well as tissue-and cell-scale stresses, 555 which provides a measure of how long each of these stresses persist.

556
Altogether, the STRESS software enables the accurate quantification of multiple 557 types of stresses, as well as their spatial and temporal characteristics, using deformable 558 microdroplets in 3D multicellular systems, from living tissues in developing embryos, to 559 3D cell culture and organoids.