Polyacrylamide Bead Sensors for in vivo Quantification of Cell-Scale Stress in Zebrafish Development

Mechanical stress exerted and experienced by cells during tissue morphogenesis and organ formation plays an important role in embryonic development. While techniques to quantify mechanical stresses in vitro are available, few methods exist for studying stresses in living organisms. Here, we describe and characterize cell-like polyacrylamide (PAAm) bead sensors with well-defined elastic properties and size for in vivo quantification of cell-scale stresses. The beads were injected into developing zebrafish embryos and their deformations were computationally analyzed to delineate spatio-temporal local acting stresses. With this computational analysis-based cell-scale stress sensing (COMPAX) we are able to detect pulsatile pressure propagation in the developing neural rod potentially originating from polarized midline cell divisions and continuous tissue flow. COMPAX is expected to provide novel spatio-temporal insight into developmental processes at the local tissue level and to facilitate quantitative investigation and a better understanding of morphogenetic processes.

contrast to oil microdroplets, hydrogel probes are compressible and therefore allow quantification of local pressure changes by tracking the bead volume change. However, as the shear or elastic moduli of these PAAm beads remained undetermined, their application was limited to the identification of isotropic stresses. More recently, Lee et al., 2019, confirmed the applicability of PAAm beads for the determination of mechanical stress in multicellular spheroids 26 . In contrast to the PAAm beads introduced by Dolega et al., ellipsoidal shape deformations could be evaluated. However, a monodisperse bead production was not realizable precluding a reliable in vivo application of these bead sensors.
In another study, elastic round microgels (ERMGs) made of alginate, loaded with fluorescent nanoparticles, were used to quantify isotropic and anisotropic compressive stresses in living tissues by tracking nanoparticle displacement 27 . This represents a promising method for identifying stresses, both in vitro and in vivo. The quantitative analysis of microsphere deformations rests on the definition of a stress-free spherical shape as reference configuration. Thus, the accuracy of this general approach critically depends on the availability of compressible, elastic spheres whose shape, size, and mechanical properties are well established prior to injection.
Here, we present a method that circumvents the limitations of existing cell-scale stress sensors. We validate PAAm as suitable material to create cell-like, spherical beads with a narrow size distribution using droplet microfluidics. We extensively characterize both the compressibility and elastic modulus of the beads to enable the quantification of isotropic and anisotropic stresses. This detailed determination of size and mechanical properties before injection provides the methodical basis for the in vivo use of the computational analysis-based cell-scale stress sensing (COMPAX). Importantly, it removes the necessity to recover the beads observed after the experiment in order to acquire knowledge of their individual stress-free reference configuration. We demonstrate the utility of COMPAX using PAAm beads by spatially and temporally quantifying for the first time the stresses acting in the developing zebrafish neural rod. We show the presence of oscillatory pressure propagation, which may potentially be generated by polarized cell divisions, and quantify compressive stress distributions within the tissue during neural rod development. COMPAX will not only provide novel spatio-temporal insight into morphogenetic processes in embryos, but can also be used for quantifying stresses in adult tissues and organoids, or when beads encounter individual cells during phagocytosis or migration.

Results
Fabrication and mechanical characterization of PAAm beads. The production of standardized celllike elastic PAAm beads has recently been described by our group elsewhere 28 . In order to use these beads as force sensors, we adjusted their size and elastic properties to be almost uniform and similar to those of cells. Spherical PAAm beads were prepared by controlled polymerization of acrylamide and N,N´-methylenebisacrylamide in a flow-focusing microfluidic device (left inset Fig. 1a) to obtain uniform particle size with a diameter of 17.0 µm ± 0.5 µm (mean ± SD, Fig. 1a). By altering the total monomer concentrations, the elastic modulus of the beads was adjusted to mimic the typical stiffness range of eukaryotic cells 29 (assuming that cells can be represented as a homogeneous isotropic elastic material at the lowest order). To render the inert PAAm beads bioadhesive, and to allow a visualization of the stress response in real time, they were covalently modified with Poly-L-Lysine (PLL) conjugated with Cy3 fluorophores (PLL-Cy3) via NHS-ester, after production (right inset Fig. 1a) 28 .
Atomic force microscopy (AFM)-based colloidal probe nanoindentation was applied to determine the elastic modulus of the PAAm beads in an initial (small) strain regime and to analyze the effect of different processing steps and environmental parameters on bead stiffness. Figure 1b shows the effect of PLL-Cy3 functionalization on the elastic modulus of the PAAm beads. Unmodified PAAm beads exhibited a Young's modulus of 1.4 ± 0.6 kPa (mean ± SD). PAAm beads which were modified with NHS-ester did not significantly increase bead elasticity (1.6 ± 0.6 kPa). The addition of PLL-Cy3 to NHS modified beads slightly increased bead elasticity to 1.8 ± 0.7 kPa. This value was then used for further analysis of the bead deformation in vivo. The insensitivity of the mechanical bead properties to different physiological temperatures was confirmed by AFM measurements of the same PAAm beads at temperatures from 24-50 °C (Fig. 1c). At all conditions tested, the elastic modulus of the beads remained constant. Young's moduli determined using AFM were validated by numerical reconstruction of the indentation process and revealed that the reconstructed radial displacement matched the increase in bead diameter obtained by confocal microscopy during colloidal probe indentation ( Supplementary Fig. 1).
Measuring the volume variation of the PLL-Cy3 PAAm beads as function of osmotic pressure using dextran solutions 25,30 was used to determine the bulk modulus (Fig. 1d). Linear fitting of the curve at small strains in the region of elastic deformation resulted in a bulk modulus of ĸ = 6.2 ± 0.2 kPa (mean ± SD). Considering the Poisson ratio ν, which was previously determined to be 0.443 28 , the correlation between bulk modulus ĸ and elastic modulus E can be described using the equation E = 3ĸ (1-2ν), which corresponds to a Young's modulus of 2.1 ± 0.3 kPa. This independently confirmed the value obtained by AFM indentation.
In addition, plate compression of PAAm bulk gels also verified bead elasticity at small compressions ( Supplementary Fig. 2). Simultaneously, these measurements displayed a non-linear material behavior of PAAm at large strains, which precluded the application of the small strain (geometrically linear) framework. However, for moderate strains up to 15% the material response was only slightly nonlinear, which indicated the possibility to use a rather simple finite strain formulation (e.g. Neo-Hookean model) for the analysis of bead deformations.
Compressibility of the PAAm beads was further confirmed by incorporation of the beads into human mesenchymal stromal cell (MSC) aggregates ( Supplementary Fig. 3). At 24 h and 48 h of culture, PLL-Cy3 functionalized PAAm beads showed significant volume reduction. After dissociation of the cells after 14 days of culture, the bead volume returned back to its initial value (shape recovery), illustrating their elastic material properties even over extended periods of time. Together, this detailed material characterization of the PAAm beads' size, compressibility and elastic modulus provide the basis for subsequent computational analysis of in vivo bead deformations.
deformations are analyzed and converted to stresses using a newly developed method called computational analysis-based cell-scale stress sensing (COMPAX). COMPAX consists of three basic steps. In step one (Fig. 2a), we registered two-dimensional confocal images of the deformed PAAm beads embedded in the tissue of a zebrafish embryo. Since we were aiming to capture the detailed shape of the beads while avoiding as much noise as possible in the confocal images, the distance between the individual two-dimensional sections was chosen as 1 μm. Based on this, we constructed a three-dimensional geometry of the deformed PAAm bead and a suitable finite element (FE) discretization using the 3D visualization and processing software Avizo. The first step in geometry generation, using the confocal images, involves the segmentation of voxels belonging to the microbead. Then the geometry of the bead is interpolated between confocal images according to certain criteria. Since the geometry creation process is a fixed order, the quality of COMPAX results remains unaffected.
In step two (Fig. 2b) the undeformed configuration and approximate boundary conditions were estimated. Since each deformed bead cannot (yet) be assigned to its undeformed configuration due to technical reasons, the undeformed state is generally not known. However, encouraged by the almost uniform fabrication of the PAAm beads resulting in narrow distributions in shape and diameter, we considered the reference geometry of the beads to be a sphere with a mean diameter of 17.0 μm. To construct the FE mesh for the respective reference configuration, we computed radial distance vectors pointing out from the surface nodes of the deformed mesh to the assumed spherical reference geometry and utilized those distance vectors as Dirichlet boundary conditions during a preprocessing FE analysis. The resulting deformed configuration was then used as the mesh for the assumed undeformed geometry of the bead. At this point, an approximation of the mean pressure is already feasible by computing the volume change of the bead based on the meshes produced (see section 'Numerical simulations'). However, since this approximation involves the simplifications of the theory for small strains, the accuracy is limited.
In step three ( Fig. 2c) we applied the inverse radial distance vectors as surface displacements to the mesh constructed in step 2. Then, the main FE simulation was performed to obtain the stress distribution that was assumed to mimic the one in the real PAAm bead. For this simulation we considered the Neo-Hooke material model in a geometrically nonlinear continuum mechanics setting that allowed for large displacements (see also Supplementary Fig. 2). The two material parameters E and ν were estimated from the linearized situation obtained previously, i.e. they coincided with the Young's modulus and the Poisson ratio from the linear elasticity framework that describes the initial stress-strain response under small strains. Considering large strains, the model yielded a slightly nonlinear increase in stresses with increase in strains. This was motivated by the fact that large displacements occurred in the bead; however, the resulting stress-strain response was only slightly nonlinear in the regime of large strains up to approximately 10% ( Supplementary Fig. 2). Note that strains in the beads analyzed were mostly of this level, while local maximal values reached to approximately 20%. Numerical validation. COMPAX constitutes a new approach to quantify cellular stresses inside soft tissues, which permits the determination of isotropic stresses as well as distinct shape changes of the beads. Therefore, a numerical validation was of major importance to evaluate the accuracy of the computed stress state from a qualitative and quantitative point of view. For this purpose, we analyzed several artificially created load scenarios applied to a spherical bead with a diameter of 17.0 μm, a Young's modulus of 1.8 kPa, and a Poisson ratio of 0.443. These artificially deformed beads served as virtual experiments, that compared stresses obtained using COMPAX by analyzing the shapes resulting from the reference simulation with known data used for the simulation (Fig. 3). Note that the same diameter and material parameters were used for COMPAX as for the reference simulation such that the differences in the results should be on the order of computer accuracy in a perfect scenario.
The quality of our approach was assessed via four criteria: the local differences in surface pressure Pres σ , the magnitude and direction of the volumetric mean of the principal stresses 1 Ø σ , 2 Ø σ and σ 3 Ø , and the volumetric mean of pressure σ Pres Ø . Here, pressure refers to negative hydrostatic stress, defined as the arithmetic mean of normal stresses. Positive values of the pressure indicate compressive isotropic stresses (for more information about definitions, see section 'Numerical simulations'). Note that by "normal" stresses, contrary to shear stresses, we refer to those components of the stress tensor which are associated with the normal directions of the particular cross-sections in the body. As a plausibility check, we first compared the results for a homogeneous surface pressure of 1000 Pa, where the resulting displacements are indeed radial. The resulting stresses of COMPAX differed from the virtual data by values in the order of computer accuracy, implying that the method had been correctly implemented.
In a first relevant test scenario, we defined a periodic surface pressure whose intensity varied sinusoidally with circle coordinates angle ϕ in the x-y-plane from a minimum of 600 Pa to a maximum of 1000 Pa with a wavelength of π 4 . In this case, the comparison revealed an excellent agreement between virtual data and the COMPAX-method (Fig. 3a). In particular, the almost identical contour plots indicate that the detection of large fluctuations in the surface pressure represents no difficulty for the COMPAX-method.
The second test scenario considered a quadratic distribution of surface pressure along one rotation axis of the bead with a maximum pressure of 1000 Pa at the equator and a minimum pressure of 600 Pa at the poles. The contour plots exhibited some local quantitative differences, while also the values of the volumetric means of the principal stresses varied (Fig. 3b). Nevertheless, the direction of principal stresses as well as the mean pressure showed good conformity. In a final load scenario, we combined a homogeneous surface pressure of 800 Pa with www.nature.com/scientificreports www.nature.com/scientificreports/ two different shear loads of 120 Pa. Each of the two loads was applied tangentially on one of the two halves of the bead. The first and second loads were oriented from the north to the south pole and from the east to the west pole, respectively. Comparing the contour plot in Fig. 3c, significant quantitative and qualitative differences in the local pressure distributions could be observed, even though the directions of the principal stresses and volumetric mean of the pressure remained comparable. In addition, Supplementary Fig. 5 shows contour plots of the shear stresses for the final load scenario. As expected, the shear stresses generated by the COMPAX-method differ from the shear stresses in the reference simulation. This is due to our assumption of radial displacements at the surface of the bead which results only in moderate shear deformations.
These artificial situations did not yet include potential errors resulting from measurement uncertainties e.g., associated with the shape of the undeformed bead and the material parameters. Therefore, an uncertainty analysis was performed for diameter and Young's modulus, which showed an influence of up to 30% standard deviation in volumetric mean pressure as a result of measured differences in these two parameters ( Supplementary Fig. 4).
Taken together, COMPAX was less reliable when extracting surface shear and larger non-radial deformations, which led to local deviations in stress computation and simultaneously affected the values of principal stresses. However, we obtained very satisfying results for the determination of mean pressure and directions of the principle stresses, which were thus primarily analyzed in the experiments.
Quantifying cell-scale stresses during zebrafish development. Cell-generated forces that emerge during neural rod formation in zebrafish embryos were quantified by computationally reconstructing the deformations of embedded PLL-Cy3 PAAm beads. The beads were microinjected into the developing neural plate at the tailbud stage (10 hours post fertilization (hpf)) and time-lapse confocal microscopy was used to visualize three-dimensional shape changes of the beads in the neuro-epithelium (Fig. 4).
First, we analyzed cell-induced deformations of a bead embedded in the developing neural rod of a zebrafish embryo. Bead and embryo were imaged for different periods of time at 14 hpf, 15 hpf, and 20 hpf (Fig. 5a). Throughout the imaging period, the bead remained constantly positioned in the basal region of neural rod adjacent to the developing otic placodes of the embryo. At the indicated time points, COMPAX showed that the surrounding cells exerted forces on the bead that induced its deformation (see the deformed shapes in Fig. 5b). www.nature.com/scientificreports www.nature.com/scientificreports/ Here, we wanted to distinguish between local stresses on distinct areas of a bead and the global volumetric mean of pressure.
While local tensile isotropic stresses (negative pressure values) were initially detected at 14 hpf, volumetric compression (positive pressure values) dominated at later time points. The temporal evolution of volumetric mean of the pressure for each time period, as well as the corresponding volume of the PAAm bead, is shown in Supplementary Fig. 6a. Additionally, COMPAX-based rendering of local pressure could be related to the activity of adjacent cells, where cell spreading/division caused greater local compressions. An exemplary representation of such pressure distributions can be found in Fig. 5c, in which contour plots of the distributed pressure in selected cross-sections of the bead are depicted within the corresponding confocal images of the z-stack. To quantify prevailing stress directions during neural rod formation, we computed the normal stress components (acting in the x-y-plane) as a function of the direction in space described by the angle ϕ (Fig. 5d). For better visualization, we amplified the curviness of the resulting stress curves by amplifying the fluctuation around the mean. This was achieved by first calculating the mean values of the normal stresses acting in the x-y-plane and secondly scaling the differences of the normal stresses to the mean by a factor of five. This allowed us to detect the directions of the main normal stresses in the x-y-plane for serial time points (see red marks in Fig. 5d). In addition, we were interested in the development of these normal stresses over time. Although the mean value at each time point was not modified by the amplification, different intensities of fluctuations led to absolute values of the normal stresses which are not comparable for different time points. Therefore, we normalized the normal stresses by the maximal absolute value of normal stress.
We observed a continuous pulsatile pressure propagation pattern in the basal part of the neural rod during the evaluated periods. The intensity of pressure oscillations decreased with time but were clearly detectable for the entire period of evaluation. Decreased changes between stress distributions of serial time points implied the dominance of the global stress state over stress changes caused by local developments. Furthermore, deformation of the beads indicated a change in stress orientation during the first imaging period at 14 hpf, as evidenced by www.nature.com/scientificreports www.nature.com/scientificreports/ the change in the direction of the maximum value of normal stress. Stress orientation remained unaltered during subsequent periods.
This occurrence of periodic pressure fluctuations during zebrafish development was confirmed by microinjecting additional beads into the neural plate and tracking their deformations at identical developmental stages (biological replicates). A PAAm bead at a position comparable to the bead shown in Fig. 5 -in the developing structure of the otic placode (basally localized) -also exhibited oscillatory normal stresses at 14 hpf and 15 hpf (Fig. 6a, Supplementary Fig. 6b). In contrast, a bead positioned close to the midline (apically localized) of the neural rod within the same embryo (Fig. 6b, Supplementary Fig. 6c) was simultaneously exposed to pulsatile tensile stresses. Deformation patterns of another bead embedded near the midline in another embryo (Fig. 6c,  Supplementary Fig. 6d) also confirmed the overall trend of oscillatory stress propagation, but showed compressive stresses only at 14 hpf, which might be explained by the continuous movement of this bead from the midline towards the basal border of the neural rod (inset Fig. 6c).
To visualize prevailing forces after the opening of ventricles in the neural rod, deformations of a bead located at the ventricle at 19 hpf were imaged (Fig. 6d, Supplementary Fig. 6e). The bead was surrounded by cells on two sides but unattached on the other two sides. COMPAX analysis revealed oscillating compressive stresses within the developing cerebral rod. Overall, these results suggest the presence of spatially varying oscillatory stresses during development of the zebrafish neural rod.

Discussion
The present work has demonstrated that PLL-Cy3 functionalized PAAm beads can be efficiently injected into living tissue/embryos to quantitatively determine in vivo cellular stresses. For the first time, a reliable reference configuration has been proposed to numerically analyze bead deformations that deviate from the spherical shape in vivo. As a result, the COMPAX-method reconstructed credible data on pressure changes as well as the direction of main shape changes. Detailed mechanical characterization of the beads and elaborate computational analysis of in vivo bead deformation revealed pulsatile stresses during zebrafish embryo development.
The preparation of PAAm beads using standardized microfluidic flow focusing 28 resulted in an extremely narrow size distribution which, to the best of our knowledge, has not been achieved in other studies. In general, computation of bead shape changes is based on the knowledge of the stress-free state as reference. Mohagheghian et al. 27 have shown that this reference configuration could be accurately determined by cell removal during in vitro application of alginate microspheres. However, as this is not realizable in vivo, the initial state after injection was used to define the stress-free state, which introduced uncertainties in stress calculation. In contrast, the size uniformity of the PAAm beads used here was essential to render COMPAX insensitive to variations in initial size and permitted the application of this approach to living developing tissues/animals over extended periods of time.
The elastic modulus of the PAAm beads (~2 kPa) was similar to that of cells or tissues which, once injected into a living tissue/embryo, allowed the surrounding cells to deform the beads, and thus enabled sensitive cell force studies at different stages of zebrafish embryonic development. Since knowledge of elastic bead properties is an essential prerequisite for a reliable stress analysis, we performed extensive material characterization measurements directly on the PAAm beads. AFM-based colloidal probe indentation was used to determine the Young's modulus of the beads at various stages of the PLL-Cy3 functionalization process. The average Young's modulus of NHS-modified beads slightly increased after PLL-Cy3 modification due to electrostatic interactions (additional crosslinking points). The PLL-Cy3 enhances interactions with the cell membranes (negatively charged ions) and its fluorescent label enables the detection of the beads with confocal microscopy. Modification of the gel with PLL introduced positive charges that may have altered bead elasticity by introducing additional crosslinking points into the polymer network 31 . Once the beads were modified with PLL-Cy3, the whole beads were homogenously fluorescent and displayed stable stiffness at temperatures ranging from room temperature to 50 °C.
In a majority of studies on the use of colloidal microgels as stress sensors, mechanical characterization of the test system was performed on bulk gels 27,32 . Here, we have directly compared bead properties at the macro-and microscale. Pressure tests on PAAm bulk gels (identical composition as the beads) were performed and revealed similar initial elastic moduli as determined by AFM-based indentation. Nevertheless, bulk gel tests exhibited minor dispersions in curve shape compared to AFM measurements, which could be caused by variations in sensitivity to the porosity of the gels during macro-and microscale responses. Moreover, a comparison of axial material behavior of PAAm gels with the Neo-Hookean material law showed conformity for strains up to 10%. Within that range, deformations of the PAAm beads were observed. Simultaneously, a Poisson ratio of 0.443 was determined 28 . This justifies the application of the compressible Neo-Hookean material law, which allows the detection of volumetric changes, which is excluded for incompressible beads. However, a new non-linear elastic material law needs to be developed to capture with equal accuracy larger deformations that were observed locally in the beads. This conclusion is contrary to the conventional approach of previous studies, which describe hydrogel stress sensors as linear elastic materials. The resulting numerical simulations were restricted to being geometrically linear -a scenario that did not allow for an appropriate incorporation of large strains. Additionally, bulk tests have demonstrated that further investigations of the material at the microscale could facilitate development of a more detailed understanding of mechanical behavior. As AFM-based colloidal probe indentation only analyzed the Young's modulus for deformations less than 10%, the rising slope of the strain-stress curve from the bulk tests could not be captured. Both, macroscopic and microscopic material characterization were prerequisite for a detailed understanding of the material behavior which enabled us to constrain uncertainties during computations of isotropic and anisotropic stresses to a minimum.
Numerical validation of COMPAX exhibited good results for all load scenarios with four different criteria, along with significant accuracy during determination of principal stress orientation and the mean pressure. Only shear stress distribution could not be captured quantitatively using COMPAX, which can be explained by discrepancies of the real displacement direction from the assumed radial orientation. To increase the accuracy of the www.nature.com/scientificreports www.nature.com/scientificreports/ method, fluorescent and highly dispersed markers could be embedded within the beads to enable a point-wise reconstruction of displacement vectors, at least at the marker positions, and to obtain volume displacement information in addition to surface displacements. A similar approach for reconstructing shape deformations using fluorescent markers has already been used for the ERMG method 27 . However, because of the lack of a reference configuration corresponding to the shape of the stress-free bead, the marker could not be utilized to obtain detailed information on displacements. Further, another approach will be required to map the markers in the deformed configuration to corresponding markers in the undeformed configuration. For this purpose, a few different markers with a predefined distance between each other in order to re-identify the orientation of the bead in the deformed state and to map the highly dispersed markers in the undeformed configuration. This method modification could serve to identify the exact reference configuration of the undeformed initial bead and would constitute a substantial improvement. Furthermore, deconvolution, or multi-angle light-sheet microscopy offers the possibility of improving the axial resolution of COMPAX. Nonetheless, since the average bead diameter is used to construct an undeformed sphere, which was assumed to be the reference geometry for the main FE analysis, uncertainty in the bead diameter had to be considered. The narrow size distribution of the fabricated PAAm beads enabled us to constrain the standard deviation of the related uncertainty analysis to 30%. The implantation of fluorescent markers could reduce this variance even further.
Initial observations on tissue stress distributions in zebrafish embryos were possible by ERMGs made of alginate injected into the blastula and at early gastrulation stages of embryonic development 27 . In that study, the authors provided evidence for spatial differences in tensile and compressive stresses during the early stages of tissue development. As an extension of the previous study, we now show that PAAm beads can be used to quantify changes in stresses generated by cells during zebrafish neurulation, a dynamic morphogenetic process that transforms the neural plate into the neural rod. We also show that tensile and compressive stresses occur simultaneously at different positions within the developing neural rod and that they coincide with observations derived from early developmental studies that used ERMGs. Importantly, both studies indicate that spatial variations can arise in prevailing stresses during tissue formation. In the present study, we also provide evidence that compressive stresses might dominate neural rod formation between 14 hpf and 20 hpf. Such compressive stresses are attributable to persistent and dense tissue packing and associated feedback control of neurogenesis at the apical-basal polarity axis 33 .
Using PAAm bead stress sensors we show, for the first time, differences in spatial and temporal stress distributions not only during gastrulation, but also directly within the tissue during zebrafish neural rod formation. The observed pulsatile variations in relative stress fields during neurulation may possibly result from polarized cell divisions, especially for beads located in the neural rod at 14 hpf and 15 hpf. Such oriented cell divisions occur during the neural keel-rod period of zebrafish neurulation and induce a uniform distribution of polarity in the neural rod through synchronized cell divisions at the midline [34][35][36] . As neurulation progresses, tissue flow in combination with general cell shape changes 37 , cell movements along the apical-basal axis 33 and axis elongation 24 may also be a cause of the observed stress distribution. The observed oscillations are well resolved by the sampling frequency in Fig. 5d at 20 hpf and Fig. 6. Further, the amplitude of oscillations gradually changes from large oscillations at 14 hpf to small oscillations (and lower frequency) at 20 hpf (Fig. 5d) or vice versa (Fig. 6a,b), which would not occur if the effect was generated by random noise. Oscillations in morphogenetic processes are not uncommon and to be expected 38 . A good example are the well-documented and understood oscillations of the amnioserosa cells during dorsal closure in Drosophila morphogenesis 39,40 . Along these lines, due to the biological and morphogenetic processes that occur at this time, we believe that the observed variations over time are real oscillations. Since the mean values per time point are not changed by the amplification, the indication of the presence of significant oscillations is additionally supported. From Fig. 5d we observe a 20-80% deviation of maximal normal stress from one time point to the next time point. Such large deviations were not even observed for the virtual numerical analysis where exaggerated scenarios were analyzed. Therefore, an insufficient accuracy of COMPAX itself cannot explain the oscillations. Variations in the elastic modulus or radius from bead to bead (as analyzed in Supplementary Fig. 4) can also not be associated with noise in this regard because the same bead is analyzed over time and, thus, the influence from incorrect modulus or radius would be consistent at each time point. The remaining potential source for noise is the image data and the associated segmentation process. However, considering the significantly deformed beads shown in Fig. 5, the displacements at the boundary can be expected to be several orders of magnitude larger than the variation of shape in the image resulting from noise. Therefore, the significant oscillations observed for the normal stresses can also not be explained by noise in the image data. Furthermore, our observations also imply that the influence of single cell divisions related to the entire development of the fish probably decline with time as we show dominance of global stress state over stress changes due to local developments and only minor changes in stress distributions during serial time points.
In conclusion, we have used PAAm beads and COMPAX to quantify compressional forces in real-time and estimate the direction of main shape changes during neural rod formation in vivo in zebrafish embryos. The same principle can be applied to other fields of application in order to investigate prevailing stresses in vivo and in vitro 2,41 . From a developmental biology perspective, further investigation of stresses acting in tissues/organisms during early or later development as well in adult animals can be addressed using this system. Moreover, exploring stress changes in cultured cells, organoids, or monolayers will provide further knowledge on mechanical aspects of cell-cell interactions. Also at a single cell level, processes such as migration or phagocytosis, can be detected when the beads encounter individual cells. The versatility of this approach has been recently exemplified by Vorselen et al., who have used the same general approach to quantify the forces generated by macrophages during phagocytosis of PAAm beads 42 . We predict that, many other applications can be envisioned where the quantification of cell-scale stresses will improve our understanding of the role of mechanics in biology and medicine.
www.nature.com/scientificreports www.nature.com/scientificreports/ Methods Fabrication and modification of PAAm beads. The fabrication of PAAm beads using a microfluidic device as well as their subsequent functionalization with PLL has been described elsewhere in detail by Girado et al. 28 . Briefly, the generation of PAAm beads, which were sufficiently compliant to allow cell force induced deformations and enabled a convenient handling during microinjection at the same time, was achieved by using a total monomer concentration of acrylamide (AAm) and bis-acrylamide (BIS) of 7.9%. During the production process PAAm droplets were functionalized with NHS-ester to enable their modification with Poly-L-Lysine (PLL) conjugated with Cy3 fluorophores (NANOCS). A PLL-Cy3 concentration of 28 pg/bead was used to enable a homogeneous peptide functionalization.
Atomic force microscopy (AFM). The AFM indentation measurements were performed using a Nanowizard I AFM (JPK Instruments) mounted on an inverted optical microscope (Axiovert 200, Zeiss). A tipless cantilever (Arrow TL-1, nominal spring constant k = 0.035-0.045 N/m, NanoAndMore GmbH), equipped with a polystyrene microsphere (diameter: 5 µm, Microparticles GmbH), was used. For gluing the microsphere to the end of the cantilever, a two-component epoxy glue (Araldite) was used. The cantilever was calibrated by the thermal noise method before the experiment. For indentation, the cantilever tip was aligned over the center of the bead and individual force-distance curves were acquired with an approach velocity of 5 µm/s and with a contact force of 2 nN (typical indentation depth: 1 µm). Extracting the Young's modulus of the beads was realized by fitting the approach force-distance curve with the Hertz model for an spherical indenter and applying a double-contact correction considering an additional deformation derived from the counter pressure at the bottom side of the bead during indentation 28,43,44 . A cell-adhesive protein solution (CellTak, Cell and Tissue Adhesive, Corning) was used to immobilize the PAAm beads on the bottom of a petri dish to prevent bead motion during the experiment. The Young's modulus was determined using the JPK data processing software (JPK Instruments). The measurements were executed in PBS and at room temperature, except the temperature stability tests where the temperature was set to 24 °C, 30 °C, 37 °C, 44 °C and 50 °C, respectively. For AFM measurements combined with confocal fluorescent microscopy the identical AFM set up and parameters were used in combination with a confocal microscope (Zeiss 510 Meta).

Bulk modulus measurements of PAAm beads.
Compressive osmotic stress for the determination of the bulk modulus of the PAAm beads was induced by PBS supplemented with fluorescein labeled dextran (molecular weight of 2 MDa, Sigma Aldrich 52471-1 G). Well-defined amounts of dextran were utilized to ensure a controlled application of osmotic stress 25,30 . The three-dimensional shape of the PAAm beads was captured using an inverted confocal microscope (Zeiss LSM700) before and 30 min after exposing to the dextran solution, and the volume of the PAAm beads was analyzed using the open source software Fiji (3D image counter) 45 .

MSC aggregate formation.
Human telomerase reverse transcriptase (hTERT) immortalized mesenchymal stromal cells (MSCs) were transduced according to Girardo et al. 28 . and maintained under humidified 5% CO 2 atmosphere in low-glucose Dulbecco´s modified Eagle medium (DMEM, Gibco) enriched with 10% fetal bovine serum (Gibco, Invitrogen). For multicellular aggregate formation, MSCs in suspension were mixed with PLL-Cy3 functionalized PAAm beads and cultured in form of 70 µl droplets on inverted petri dish lids for 24 h, 48 h and 14 days according to the classical hanging drop method. numerical simulations. Confocal images were imported to the 3D visualization and processing software Avizo. Therein, the bead was segmented in every image utilizing segmentation tools e. g., thresholding. Based on the segmentation the 3D geometry was interpolated and, subsequently, the FE mesh was constructed. All simulations were conducted in ABAQUS, a commercial finite element software, which provides the Neo-Hookean material law as internal material subroutine.
All stress computations were based on the Cauchy stress tensor σ. Contour plots of PAAm microbeads show the distribution of surface pressure, which is defined as tr( ) where σ e is the Cauchy stress tensor at a single integration point within the finite elements and v e is the corresponding volume. We label the principal stress of the volumetric mean of the Cauchy stress tensor to be σ i Ø with = i [1,3] and the volumetric mean of the pressure Pres Ø σ . An approximation of the volumetric mean of the pressure σ Pres Ø may be computed by where κ represents the bulk modulus, V c is the volume of the current configuration and V r is the volume of the reference configuration. However, due to the nonlinearities included in the Neo-Hooke model, the volume change locally varying inside the bead will induce a locally varying, nonlinear pressure. Its volume average is not necessarily (2019) 9:17031 | https://doi.org/10.1038/s41598-019-53425-6 www.nature.com/scientificreports www.nature.com/scientificreports/ identical to the evaluation of the constitutive law at the total bead scale, i.e.  Zebrafish strain and maintenance. Zebrafish (Danio rerio) embryos and adults were obtained, raised, and maintained as described previously 46,47 . Embryos were staged as hours post fertilization (hpf) 48 . To obtain membrane-tagged eGFP embryos, adults from the transgenic strain Tg(Bactin:HRas-eGFP) (ZDB-ALT-061107-2) 49 were outcrossed with those from wild type AB strain. Furthermore, a transgenic line which expressed a nucleus-targeted venus fluorescent protein at the otx2 locus (knock-in) was established using the CRISPR knock-in strategy as previously described 50 ; this line was used to identify the midbrain hindbrain boundary.
All animal experiments were carried out in accordance with animal welfare laws and local authority requirements (Landesdirektion Sachsen, Germany) and performed in accordance with appropriate guidelines and regulations. Protocols for the generation (24-9168.11-1/2013-14) and maintenance (DD24-5131/346/11 and DD24-5131/346/12), and experimentation with transgenic animals (24-9168.24-1/2014-4) were also appropriately approved and performed in accordance with appropriate guidelines and regulations (Landesdirektion Sachsen, Germany). All experiments utilized zebrafish embryos only until 120 hours post fertilization, which do not come under animal experimentation (Directive 2010/63/ EU and in accordance with German animal protection law and Landesdirektion Sachsen, Germany). Adult transgenic animals were used only for maintaining transgenic lines and breeding purposes.
Zebrafish microinjection. The embryos were grown in a 28.5 °C incubator, dechorionated at 9 hpf using pronase (Sigma-Aldrich, P8811) and incubated in calcium-free ringer solution (recipe from the zebrafish information network (ZFIN) database; http://zfin.org) for 10 minutes. The bead implantation was carried out in the same solution. Standard borosilicate capillaries without filament (World Precision Instruments) were pulled using a micropipette puller (Sutter Instruments). The tips were cut to obtain an opening of about 15 µm and back-loaded with the injection solution containing hydrogel beads. The beads were implanted into the developing neural plate close to the prospective midbrain hindbrain boundary using a microinjector system (PV820, pneumatic picopump, World Precision Instruments). The injected embryos were incubated in E2 buffer (recipe from ZFIN protocol) until imaging. The implantation of the beads did not result in any morphological or developmental abnormalities. Fish were pre-screened for successful implantation in the desired location. The selected embryos were mounted as described below and imaged. embryo mounting and imaging. Embryos with beads in the desired location were anesthetized using MS-222 (Sigma Aldrich, A5040) in E2 solution and dorsally mounted on a glass-bottomed petridish (MatTek) in 1% low melting agarose prepared in E2 solution (Fig. 4a). Time-lapse images were acquired using an inverted confocal microscope (Zeiss LSM780) using a 40x water immersion objective (NA 1.2) with laser lines 488 nm for eGFP and 561 nm for Cy3. Images were analyzed using the open source software FIJI 43 .
Statistical analysis. The bin widths of the histogram (Fig. 1a) was chosen according to the Freedman-Diaconis rule 51 . In the box plots (Fig. 1b,c; Supplementary Fig. 4b) the mean is shown as filled square symbol, the median as straight line, and the boxes are determined by the 25 th and 75 th percentiles. The whiskers represent the standard deviation and the 1 st and 99 th percentiles are indicated by crosses in Fig. 1b,c. In Fig. 1d data points with error bars indicate mean ± standard deviation. The number of measurements (n) is given at the respective data point boxes ( Fig. 1d; Supplementary Fig. 4b) and in the figure description (Fig. 1b,c), respectively. The data in Fig. 1b and Supplementary Fig. 4b was analyzed using OriginPro Software. Both data sets rejected normally (Shapiro-Wilk test). To evaluate the statistical difference, the non-parametric Wilcoxon-Mann-Whitney test was used. The asterisks define the statistical difference as follows: ***P < 0.001.