Superresolved and reference-free microparticle traction force microscopy (MP-TFM) reveals the complexity of the mechanical interaction in phagocytosis

Force exertion is an integral part of cellular behavior. Traction force microscopy (TFM) has been instrumental for studying such forces, providing both spatial and directional force measurements at subcellular resolution. However, the applications of classical TFM are restricted by the typical planar geometry. Here, we develop a particle-based force sensing strategy, specifically designed for studying ligand-dependent cellular interactions. We establish an accessible batch approach for synthesizing highly uniform, deformable and tunable hydrogel particles. In addition, they can be easily derivatized to trigger specific cellular behavior. The 3D shape of such particles can be resolved with superresolution (< 50 nm) accuracy using conventional confocal microscopy. We introduce a computational method that allows us to infer surface traction forces with high sensitivity (~ 10 Pa), directly from the particle shape without embedding tracers in the particle. We illustrate the potential and flexibility of this approach by measuring forces throughout phagocytic engulfment. This strategy can readily be adapted for studying cellular forces in a wide range of applications.


INTRODUCTION
Cells continuously exert forces on their environment, allowing them to probe and respond to their mechanical surroundings 1 , while also engaging in interactions with neighbouring cells [2][3][4] . Traction force microscopy (TFM) provides the ability to study such forces in detail 5,6 . However, classical TFM is largely limited to in vitro applications 7 , and a range of cellular processes warrant different (i.e. non-planar) geometries, such as cell-cell interactions like phagocytosis or the formation of the immunological synapse, in which mechanical forces have a widespread role 8,9 . Moreover, with a few notable exceptions [10][11][12] , TFM focuses on the study of shear forces, often neglecting forces normal to the surface.
Force sensing within tissueand in vivois a new area of exploration and first became possible by the introduction of oil droplet-based technologies 13 . Recently, deformable hydrogel microparticles (MPs) were also used as force reporters 14 . The use of microgels for MP based traction force microscopy (MP-TFM) shows promise for broad applicability because of the mechanical tunability over orders of magnitude and the potential to measure both shear and normal stresses. However, synthesis of suitable hydrogel MPs for such applications can be rather complex, often requiring specialized expertise in microfluidics 14,15 . Moreover, particle properties, especially their uniformity in mechanical properties, must be optimized and thoroughly characterized for accurate traction measurements 15 . Like in TFM, multiple strategies are possible for measuring deformations and for calculation of the stress fields 16 , including reference-free methods 17 , which have not yet been explored in MP-TFM. Finally, deformable microparticles have so far only been used as passive force reporters without biologically specific ligands, despite the fact that MPs can be chemically modified to trigger and elicit specific cellular behaviors.
Here, we describe a novel particle-based force sensing strategy, specifically designed for studying ligand-dependent cellular interactions ( Supplementary Fig. 1). We established a batch production approach for synthesizing deformable tunable hydrogel particles that is simple and reproducible. The resulting particles are homogeneous throughout their volume and from particle to particle. The incorporation of carboxyl groups allows easy and efficient derivitization. Using confocal microscopy, we can

Batch production of homogeneous, readily conjugatable, hydrogel microparticles with cell-like rigidity
We used membrane emulsification with Shirasu Porous Glass (SPG) membranes with uniform pore size to produce monodisperse spherical hydrogel microparticles 19,20 .
Acrylamide mixtures containing monomeric acrylamide (AAm), acrylic acid (AAc) and the crosslinker N,N'-methylenebisacrylamide (BIS) were extruded through a hydrophobic SPG membrane under nitrogen pressure, generating droplets in an oil phase ( Fig. 1a). After the entire aqueous mixture was dispersed, acrylamide polymerization was initiated within the droplets by addition of the free radical initiator Azobisisobutyronitrile (AIBN). The resulting deformable poly-AAm-co-AAc microparticles (DAAM-particles) were collected and resuspended in PBS. As shown previously 21,22 , this strategy yields highly monodisperse particles (Fig 1b). The particle size could be tuned between 4 -15 m (each with coefficient of variation (CV) ≲ 0.1) using SPG membranes with various pore sizes, revealing the expected linear relation between pore and particle size 22 (Fig.   1c). Furthermore, since membrane emulsification is a bulk method, we were able to synthesize an estimated 10 10 -10 11 particles per batch ( Supplementary Fig. 2). To further tune important particle properties, such as their rigidity, we changed the composition of the acrylamide mixture. In particular, the total mass concentration of acrylic compounds (CT = CAAm + CAAc + CBIS = 100 mg/mL), as well as the relative concentration of acrylic acid (mAAc/(mAAm + mAAc + mBIS) = 10%) were kept constant, while the crosslinker concentration Cc = mBIS/(mAAm + mAAc + mBIS) was varied between 0.3 and 2.3 % (Supplementary Fig. 2). As expected, Cc affects the hydrogel swelling properties, leading to an inverse relation between Cc and particle size (Fig. 1d).
Coating MPs with ligands is critical for triggering specific cellular behavior and conjugation with fluorescent molecules is required for visualization of the particles and their deformation in microscopy applications. We incorporated acrylic acid into the gel mixture to allow straightforward and efficient covalent coupling of primary amine groups to the DAAM-particles using carbodiimide crosslinking chemistry 23 . DAAM-particles were functionalized with both fluorescently conjugated bovine serum albumin (FITC-BSA) and small fluorescent molecules (TRITC-Cadaverine (Cad)) ( Fig. 1e). Confocal not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint imaging revealed highly homogeneous coupling of both protein and reactive dye to the microporous DAAM-particles. The ideal force-reporting microparticle has a homogeneous and isotropic network structure, which underlies equally homogeneous elastic properties. However, initiating polymerization at the droplet boundaries, or subsequent gel swelling, could potentially lead to a radially inhomogeneous network structure. Derivatization of the particles could also affect particle mechanical properties and lead to further inhomogeneities. Therefore, we evaluated the radial fluorescent intensity profile (n > 100 particles) of both FITC-BSA and TRITC-Cad conjugated to DAAM-particles. (Fig. 1f, Supplementary Fig. 2) Simulation of the fluorescent signal of a perfectly homogeneous sphere, by convolution with an experimentally determined point spread function (PSF) ( Supplementary Fig. 3), accurately reproduced the experimental radial intensity profile, independent of crosslinker density (Fig. 1f, Supplementary Fig. 2). The lack of radial deviations from a homogeneous sphere, combined with the homogeneous appearance of individual particles, strongly suggests that these particles have uniform and isotropic polymer networks.
In addition to intraparticle homogeneity, strong interparticle uniformity is critical for application of microparticles as force sensors. To evaluate particle-to-particle variation, we measured quantitative 3D refractive index (RI) maps of individual DAAMparticles (Fig. 1g). These measurements revealed RI values of 1.3349 ± 0.0001 (standard deviation (s.d.), n = 8) to 1.3381 ± 0.0001 (s.d., n = 29) for particles with Cc = 0.32% and Cc = 2.3%, respectively (Fig. 1g). The coefficient of variation (CV: calculated as the standard deviation divided by the relative RI difference between particle and medium) within each condition were markedly low at 0.03 -0.1. Moreover, the similarity of the standard deviations between samples suggests that the observed variation may be dominated by the measurement error, and the particle-to-particle variation in polymer network density is potentially even lower.
The close match of RI between these particles (1.3349 -1.3381) and the medium (RIPBS = 1.334) also represents a significant benefit over glass (RIGLASS = 1.45) or polystyrene (PS, RIPS = 1.59) beads for many microscopy applications, including their use as target models for phagocytosis. The optical distortions, i.e. the reflection and not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint refraction of light, caused by a particle depends on RImedium -RIMP. For DAAM-particles this value is 0.001 -0.005, whereas for glass or polystyrene these values are typically 100-fold higher (Fig. 1h). We compared confocal images of DAAM and PS particles with volumetric and surface labels (Fig. 1i). Z-slices of both particles revealed apparently equivalent signals. Axial slices of the volume of PS-particles, however, showed a distorted shape (the upper hemisphere appears axially elongated by ~RImedium /RIMP) and an inhomogeneous signal. Furthermore, the fluorescent signal from the upper hemisphere's surface is reflected off the PS-particle surface and is largely not captured on the camera. Such lensing effects are barely noticeable for DAAM-particles. Finally, due to the porous nature of DAAM-particles, the advantageous optical properties are rather independent of the medium, as evidenced from repeating these measurements in vectashield with RI 1.45 (Fig. 1h, Supplementary Fig. 4).
To allow quantitative force sensing with deformable DAAM-particles, we characterized their mechanical properties by atomic force microscopy (AFM). We used pyramidal tips with a large end radius (~600 nm), which allow imaging of spherical particles, and indentations in the Hertzian regime because of the large end radius (Fig.   1j). Imaging of single particles before indentation allows precise localization of the centroid of the particle for subsequent indentations (Fig. 1j). Moreover, it provides an estimate of the particle shape, revealing that DAAM-particles are significantly deformed due to surface adhesion, forming spherical caps with contact radius > 4 um ( Supplementary Fig. 5). This renders the indentation-induced deformation of the bottom of the DAAM-particles, which is in contact with the surface, negligible compared to the deformation of the upper part, which is in contact with the AFM tip. Nano-indentations at the particle's centroids showed typical Hertzian force responses and revealed cell-like rigidities, as quantified by the Young's moduli (Ey) between 0.3 and 3.9 kPa, for Cc = 0.32 -2.3 % (Fig. 1k). For the more rigid particles, we found up to 1.5-fold higher moduli after conjugation of the particles with BSA. Importantly, the Young's moduli are 5-10 fold lower than bulk hydrogels of similar composition ( Supplementary Fig. 5). Such discrepancy may be caused by the differences in polymerization conditions (e.g. different initiator and temperature) or gel swelling, and illustrate the importance of direct mechanical characterization of MPs. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Superresolved shape analysis of DAAM-particles
To determine the accuracy with which the particle shape can be reconstructed and to generate a reference shape measurement of undeformed particles, we imaged adherent particles in Vectashield mounting medium (RI ~1.45). The 3D shape of individual particles was reconstructed by processing confocal imaging stacks, encompassing deconvolution, and edge detection and localization steps ( Supplementary Fig. 6). We varied the laser power to obtain particle images with various signal-to-noise ratios (SNR) ( Fig. 2a) and estimated our localization precision as the local surface roughness Rrms observed over a small area (1 m 2 ). Only the DAAM-particles' upper hemispheres were analyzed to exclude effects of deformation due to adhesion to the glass substrate. Due to the non-isotropic resolution in confocal microscopy, the local Rrms depended on the position on the particle. Therefore, we separately quantified the Rrms at the equator (latitude < 15°) and at the particle apex (latitude > 75°), which are dominated by the resolution in z and xy respectively (Fig. 2b). Even at modest SNR (~15) the localization precision in z (~200 nm) and xy (~40 nm) exceeded the nominal resolution (full width at half maximum of the PSF) in the corresponding direction. The principle enabling us to superresolve the particle boundary is similar to single molecule localization techniques 24,25 (Supplementary Discussion). With increasing SNR (>100), the localization further improved (z: ~40 nm, xy ~20 nm), also leading to a 2.5-fold reduction in the resolution anisotropy (Fig. 2b).
Aside from edge localization precision, the accuracy with which deformations can be measured (and forces recovered) also depends on the particle shape in the absence of externally applied forces. Deviations from a perfect spherical shape can, for example, be caused by deformation of the emulsion droplets during gelation. In addition, apparent shape deviations can arise from imaging artifacts like drift during recording of confocal stacks. To measure such deviations we quantified 1) the global Rrms, calculated over the entire upper particle hemisphere and 2) the sphericity, defined as the surface area of a sphere with equal volume to the particle divided by the observed particle surface area (1 for a perfect sphere). Global Rrms measurements revealed small (50  25 nm) overall surface roughness, corresponding to deviations of less than 1% of the particle radius ( Fig. 2c). Sphericity calculations corroborated that DAAM-particles are very spherical not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Macrophage-like cells induce localized target deformations during phagocytosis
Cellular forces have an important role during target engulfment in phagocytosis.
Such phagocytic forces have previously only been measured in the frustrated state on flat substrates 26 , and our conjugatable and deformable DAAM-particles present an ideal tool for measuring them on particles with more physiological curvature. We exposed J774 murine macrophage-like cells to soft DAAM-particles (0.3 kPa, Cc = 0.32%) that were labeled with TRITC-Cad and functionalized with BSA and anti-BSA immunoglobulin G (IgG). J774s showed strongly ligand-dependent attachment to, and internalization of, antibody-coated DAAM-particles ( Supplementary Fig. 7). After fixation of the cells, a secondary antibody was added, which is sterically unable to diffuse into the cell-target contact area 27 . Because the secondary antibody only binds the exposed DAAM-particle surface, it provides a measurement of the degree to which the phagocytic process has progressed (Fig. 3). Confocal images of DAAM-particles undergoing phagocytosis revealed local cell-induced deformations at all stages of the process.
To comprehend the full extent of phagocytic target deformation by macrophagelike J774s, the 3D shape of particles in various stages of phagocytic engulfment was reconstructed ( Supplementary Fig. 6, Fig. 4). In addition, determination of the fluorescent intensity of the secondary antibody over the particle surface allowed determination of the contact area between the phagocyte and the DAAM-particle, and the location of the base of the phagocytic cup ( Supplementary Fig. 6, Fig. 4). This analysis revealed complex and distinct patterns of deformation, suggesting an intricate series of mechanical interactions throughout phagocytosis.

Reference-free estimation of normal and shear stresses from particle shape
To reveal the cellular forces exerted on phagocytic targets, we calculated the traction forces from the observed particle shapes. In classical TFM, the displacement not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint field is measured directly, and hence estimating the traction forces on the surface is relatively straightforward. In our method, the 3D shape of DAAM-particles is measured with high resolution instead. The surface displacement field is not uniquely determined by the surface shape, since multiple displacement fields can lead to the same shape. To derive traction forces, we thus solved the inverse problem of inferring the traction forces given the observed particle shapes and the traction-free regions from the fluorescent immunostaining ( Supplementary Fig. 6). This is accomplished by an iterative optimization procedure described below (Fig. 5a).
We start with a trial surface displacement field ( , ), where , are co-latitude and longitude respectively, from which we can rapidly obtain the full elasticity solution, including the stress field everywhere in the particle, the traction force ( , ) on the surface, and the elastic energy using the spherical harmonics-based method 18 . Given the surface displacement field, we also can compute the corresponding shape as described The minimization of the cost function requires the solution of the elasticity problem for millions of trial displacement fields ( , ) . Therefore, each elasticity problem needs to be solved very fast. Using the spherical harmonics-based method, we can solve each elasticity problem within 0.007 second, so that the cost function can be minimized within a day (Methods). This would not have been possible if the elasticity not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint problem were solved using general-purpose elasticity solvers such as the finite element method (FEM).

High resolution and 3D analysis of phagocytic forces and target deformation
The method presented here allows studying cellular forces throughout phagocytosis and resulting target deformation in great detail. Observations in the initial stage of phagocytosis (particle 1, 9% engulfed) provided evidence of the cell compressing its target in discrete spots (Fig. 4a, Supplementary Video 1). Computation of the mean curvature H over the particle surface clearly showed that the deformation induced by the cell consists of three distinct spots of approximately 1 m in diameter.
The compressive stresses causing this deformation were ~50 Pa (Fig. 5b). Approximately 9% of the DAAM-particle was covered by the cell and the compressive forces were located centrally in this area (Fig. 5b).
Significantly larger target deformations of up to 1 m in the radial direction were observed during the stage of pseudopod extension (particle 2, 34% engulfed) (Fig. 4b,   Supplementary Video 2). The majority of compressive deformation by the cell was localized in a ring, and was accompanied by a ~6% elongation in the direction normal to the ring. The peak of the inward deformations coincided closely with the extension of the pseudopod as estimated from the exposed surface labeling, while the indentation induced within the ring is markedly heterogeneous with clear dents and bumps (Rrms = 120 nm) Supplementary Fig. 9). This heterogeneity appears to be caused by equally varied exertion of normal stresses from 20 to 100 Pa (Fig. 5b). Shear stresses up to 50 Pa were also concentrated in the ring (Fig. 5c). Perhaps surprisingly, the cell also appeared to pull from the cup base with tensile stresses having similar magnitude to the compressive stresses. The entire engulfed surface exposed greater heterogeneity in curvature than the particle surface not in contact with the cell (Fig. 4b). In particular, two distinct pits were visible close to the base of the phagocytic cup (indentation depths: ~250 and ~700 nm). Such deformations were caused by localized compressive stresses up to 100 Pa, and were not expected since in this phase of phagocytosis the cytoskeletal actin is being cleared from the cup base 28 . not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint Observations in the stage of cup closure (particle 3, 84% engulfed) indicated compression by the cell in a rather smooth ring, and an approximately 7% particle elongation normal to the ring (Fig. 4c, Supplementary Video 3). Although deformations in the ring were less strong than observed in particle 2, the underlying normal stresses were mostly higher (50 -100 Pa), which can be explained by compression of the target throughout the cell-target contact area (Fig. 5b). The smaller variation in force also resulted in a more homogeneous appearance of the ring (Rrms = 60 nm) than in particle 2 ( Fig. 4c, Supplementary Fig. 9). Interestingly, the peak of the inward deformation and force was not localized at the cell boundary but trailing approximately ~1 m behind.
Curvature and traction patterns were visible throughout the phagocytic cup, and seemed somewhat concentrated at the phagocytic cup base with 2 notable dents (~90 Pa). Such pits, which were also visible in particle 2, resemble the size and shape of podosomes 29,30 , and could indicate a role for these mechanosensitive structures in phagocytosis.
Finally, a fully internalized particle (particle 4), which lacked fluorescent secondary immunostaining (Fig. 4d, Supplementary Video 4), was subjected to some of the strongest indentations (~ 1 m, ~600 Pa) by the cell. These strong normal stresses were also accompanied by large shear stresses up to 100 Pa (Fig. 5c). The deformations consisted of multiple (at least 8) dents, ranging between 1-2 m in diameter. The brightfield image (Fig. 3) was used to estimate the location of the cup base, revealing that these compressive deformations arise from the site of phagocytic cup closure. Perhaps this indicates that actin that accumulates during cup closure 31 contributes to the repositioning of the maturing phagosome towards the center of the cell. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

DISCUSSION
Here we have introduced batch-produced, deformable, and readily derivatizable microparticles as cellular stress sensors, as well as a novel numerical method to calculate forces exerted on them by cells. Aside from force-sensing applications, such particles can offer several key advantages as model targets for studying cellular processes, such as phagocytosis and immunological synapse formation. They allow the unique possibility to mimic the size, rigidity, and relevant chemical characteristics of cells. The tunability of the particles' Young's moduli between 0.3 and 10 kPa covers a significant part of the range exhibited by various celltypes 32 , and can likely be expanded further. In contrast, polystyrene or silica, which are often used as model targets, have Young's moduli of 2 GPa and 70 GPa, respectively. This is 10 5 -10 7 fold larger than cells, and most other biological materials. Moreover, for microscopy applications, the small refractive index difference between DAAM-particles and the surrounding medium ensures little optical distortion in imaging of particles themselves, and will also result in improved visualization of surrounding (cellular) structures. Two recent reports have described similar deformable hydrogel particles consisting of AAm or alginate that can be used as passive force sensors in multicellular aggregates and tissues 14,15 . Compared to the single microfluidic channel methods chosen previously, our batch particle synthesis method is simpler, and requires little specialized expertise or facilities. In terms of resulting particle properties, our production method leads to monodispersity (CV ≲ 0.1) approaching that reached with a single microfluidic channel (CV ~ 0.05). Importantly, our observed particle-to-particle variation in optical and mechanical properties is markedly low. For MPs with refractive index ~1.4, previously a spread on the order of 0.001 was reported 15 , while we observe a significantly lower spread of ~0.0001. Our particles also appear mechanically homogeneous, with CVs in Young's moduli between 0.08 -0.30 for stiff to soft particles. Such particle uniformity may be related to our choice for an oil-soluble polymerizing agent (Supplementary Discussion). We also do not observe any increase in variation of particle material properties after functionalization, presumably because of the efficient and homogeneous conjugation to the acrylic acid incorporated in DAAM-particles. Such efficient not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint conjugation is also crucial for robust triggering of specific cellular responses and, by allowing homogeneous staining, particle edge localization. Currently, the dominant strategy in traction force microscopy is (random) embedding of fluorescent tracer particles in hydrogels, which, when compared to a reference state, allow measurement of the displacement field and relatively straightforward calculation of the stress field. Recently, this approach was extended to deformable microparticles (~25 m) by Mohagheghian et al. 14 . Our contrasting approach, which focuses on accurate measurement of the particle boundary, has several key benefits: 1) The magnitude of deformation within elastic materials typically scales with 1/l, where l is the distance to the site of application of the stress. Thus, measuring stresses directly at the particle boundary enhances the sensitivity to small forces. Indeed, we are able to resolve small tractions (~10 Pa). 2) For evaluation of traction forces, the stress field in the elastic gel is multiplied with its surface normals. Without precise knowledge of the location of the surface, certain assumptions need to be made. For microparticles, it was previously assumed that the deformed shape was convex 10 . However, throughout phagocytic engulfment we observed many concave surface regions (i.e. areas with mean curvature H < 0, Fig. 4), which are hence critical for a complete description of exerted forces by the cell. 3) In contrast to embedding randomly distributed tracer particles, our method provides a reference-free estimate of cellular stresses. This simplifies experimental procedures, is less disruptive to the sample, and will allow larger throughput studies. Finally, our approach is best adapted for relatively small particles (~10 m), where the advantages mentioned above become increasingly important, since relatively more information is contained near the surface.
Although we have focused on phagocytosis and have illustrated the mechanical complexity of this process, the method presented here is broadly applicable in the study of inter-and intracellular forces in vitro and in vivo.
not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Synthesis of deformable acrylamide acrylic-acid microparticles
All reagents were reagent grade and purchased from major chemical suppliers unless  Table 1). placing the dialysis bag in solid sucrose, and desalted over a G25 column equilibrated with PBS.

Microparticle size measurements
For microparticle size measurements, particles were imaged with a 20x NA 0.4 objective using phase-contrast. For the smallest particles, which were extruded through 0.5 m pore size SPG filters, a 40x NA 1.3 objective was used for imaging. Images were analyzed in ImageJ, by smoothing (2x), performing edge detection, binarizing images by thresholding and filling holes. Finally, adjacent particles were separated using the watershed algorithm, and particle radii were estimated from particle areas.

Microparticle refractive index measurements
3D quantitative refractive index (RI) maps were measured on a 3D Cell Explorer (NanoLive) and analyzed using ImageJ. Z-slices of the 3D tomogram through the center of individual DAAM-particles were selected for analysis. Histograms of the RI in a small region of interest around the DAAM-particle revealed a clear bimodal distribution. RIMP -RImedium was then measured as the distance between the two peaks of the RI histogram.
RImedium was measured for PBS (1.334), or used as specified by the manufacturer as 1.45 for vectashield mounting medium (Vectorlabs, H-1000). not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Atomic force microscopy measurements
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint Before indentation, particles were imaged in peak force tapping mode (peak force setpoint 0.5 -1 nN). Subsequently, 2-5 nanoindentations were performed on the center of particles. AFM data was analyzed using custom Matlab software. DAAM-particle shape was determined as described previously for small vesicles. 34 Briefly, a spherical arc was fit to upper part of the DAAM-particle (exceeding half its height) to obtain the radius of curvature RMP, which was corrected for tip convolution by subtracting the tip radius (Rtip).
Assuming that adherent particles form spherical caps, the shape was then approximated from RMP and the particle height ( Supplementary Fig. 5). Indentation curves (or force deformation curves, FDCs) were first processed by correction for drift in the baseline and subtraction of a fit to the cantilever response, which was obtained by performing indentations on bare glass surfaces. The resulting DAAM-particle indentation curves were analyzed by fitting a Hertzian model up to 0.5 Rtip, with the contact point as a fitting parameter. For these fits, the effective radius R * , where (1/R * = 1/Rtip + 1/RMP), was used.
FDCs for the stiffer DAAM-particles (Cc = 0.65 % and Cc = 2.3 %) were analyzed above 50 pN force. Due to the low rigidity of the softest DAAM-particles (Cc = 0.32%) a force threshold of 10 pN was used. The Young's moduli (Ey) obtained from each of the 2-5 FDCs for individual particle were averaged to reduce the measurement error of the DAAM-particles. For the Young's moduli calculations, DAAM-particles were assumed to be incompressible (Poisson ratio = 0.5), which was previously shown to be a good approximation for bulk 35 and microparticle 15 acrylamide gels. For bulk gel measurements, hydrogels, covalently bound to glass substrates, were prepared as described under "PSF measurements", but with compositions, and buffer conditions, identical to the DAAM-particles (for all three crosslinker concentrations). Measurements were performed with spherical tips (SiO2) with ~300 nm radius (NovaScan Technologies). Indentation curves were analyzed as described for the DAAM-particles.

Cell culture
Murine macrophage-like cells J774A.1 (ATCC, TIB-67) were maintained and subcultured according to the methods recommended by ATCC. Briefly, cells were maintained in DMEM medium (Gibco, 11965-092) supplemented with 10% FBS not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Phagocytic assays
For phagocytic assays, cells were transferred to 12-well glass bottom plates (Cellvis, P12-1.5H-N) (1.5x10 5 cells/well). 1h before addition of DAAM-particles, the medium was replaced by L-15 medium without serum. 15 minutes before addition of the DAAMparticles Hoechst 33258 (ThermoFisher, H3569) was added for a final concentration of 1 g/mL in each well. The medium was then replaced with 200 L L-15 per well including the DAAM-particles (7.5x10 5 DAAM-particles/well), the plate was spun for 1 min. at 300 g and incubated at 37 °C. Cells where fixed by direct addition of 1 ml of 4% formaldehyde (J.T. Baker, 2106-01) in PBS after 5 min for the phagocytic deformation assay (Fig. 3-5) and after 30 min for the uptake efficiency assay (Supplementary Fig. 7). Cells were rinsed vigorously by pipetting up and down with PBS to remove particles not in contact with cells. Exposed surface of the particles was immunostained with 4 g/mL Alexa Fluor 647 donkey anti-rabbit IgG (ThermoFisher Sci., A31573) in PBS. After 30 minutes the wells were rinsed with PBS (3x) and imaged on a confocal microscope.

Confocal Microscopy
3D confocal imaging was performed using a spinning disk CSU-W1 scanner unit (Yokogawa) and a DMi6 inverted microscope, equipped with a 100x (NA 1.40) objective (Leica Microsystems). Excitation was done with solid state lasers (100 mW 488 nm, 50 mW 560 nm, 100 mW 642 nm) and images were captured on an Ixon Ultra EMCCD camera (Andor). Z-stacks with 100 nm step size were recorded using a piezo stage insert (LUDL Electronic Products, 96A900). The setup was operated using MicroManager 36 .

Microparticle 3D shape reconstruction
Image analysis was performed with custom software in Matlab. Individual microparticles were segmented based on intensity thresholding confocal z-stacks, using the dip in the distribution of the logarithm of voxel intensities as the threshold. The centroid of DAAMparticles was estimated, and they were isolated in ~20 × 20 × 30 m regions of interest not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint (ROIs) that were used for further processing. First, images were deconvolved with a point spread function (PSF), measured ~8 m from the sample surface ( Supplementary Fig. 3).
Deconvolution was performed with 25 iteration steps using a Landweber algorithm, and was implemented using DecovolutionLab 2 37 . Background subtraction was performed by averaging the Z-slices outside the ROI along the axial direction and subtraction of the resulting 2D image from each slice of the ROI. At this point, cubic interpolation along the z-axis was used to equalize the voxel dimensions. Edge detection was performed using the 3D Sobel operator 38 (Supplementary Fig. 6), after which the total edge magnitude was calculated as the root mean squared of the three resulting images. For superlocalization of edge coordinates, first a volume estimate of each DAAM-particle was made, based on the number of thresholded voxels, which was converted into a particle surface area estimate assuming a perfect sphere. The number of edge coordinates N was determined as the particle surface area divided by the area of a circle with radius 250 nm (for the DAAM-particles subjected to phagocytosis) or 500 nm (for undeformed spheres, in which case we wanted to prevent correlation between neighbouring points).
Cubic interpolation was used to estimate the intensity values along lines originating from the particle centroid in the directions of the angles determined above. Then a Gaussian function was fit to each line profile to determine the location of the edge with subpixel precision.

Smoothing of edge coordinates
Edge coordinates were smoothed (for sphericity, curvature and force calculations) using an equivalent of a 2D moving average, operating on the radial component of edge coordinates. Great circle distances between edge coordinates with indexes i and j were calculated along a perfect sphere: = arccos(sin sin + cos cos cos( − )) , not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. where R is the equivalent diameter of a sphere to the particle. The radial component of the edge coordinates was then averaged within the given window size (1 m 2 ).

Calculation of particle and surface properties
First a triangulation was generated between points with the same angles (longitude and latitude) as the particle edge coordinates, but on a unit sphere. This allows using a convex hull approach to connect all data points and not leave any gaps. The particle surface area (S) was then calculated as the sum of the surface areas of all triangles. The particle volume (V) was calculated as the sum of the (signed) volumes of tetrahedra, which were formed by connecting all triangles with an arbitrary point. The particle centroid was determined as the volume-weighted sum of the tetrahedra centroids. Before calculation of sphericity and surface curvature, the radial component of edge coordinates was first smoothed within a 1 m 2 window (see above). This was done because of the high sensitivity of these measures to high frequency noise. Sphericity was calculated as  = ( 1/3 *6*V 2/3 /S)/2 1/3 , where the division by 2 1/3 was used when only the upper hemisphere was analyzed, as it is a correction that guarantees that  cannot exceed 1 for non-closed shapes. For surface curvature calculations, first principal curvatures (k1 and k2) of the triangulated mesh were determined as described previously 39 . Finally, mean curvature H = (k1 + k2)/2 and Gaussian curvature K = k1k2 were calculated. Local root-mean-squared surface roughness (Rrms) measurements were performed by using great circle distances to find points within a 1 m 2 area (as above, typically 6-8 points), and then taking the standard deviation of the radial component of these points. Global Rrms was calculated similarly, but over the entire upper hemisphere.

Stress-free boundary analysis
Analysis of the immunostaining of the free particle surface was performed in custom Matlab software (Supplementary Fig. 6). First, images were deconvolved as described above. Then, cubic interpolation was used to estimate the intensity values along lines originating from the particle centroid (as determined during particle shape analysis). A regular grid in spherical coordinates (as opposed to an approximately equidistant grid used in particle shape analysis) was used, to simplify later image processing steps (in not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint particular the active contour algorithm described below). The number of points in the regular grid was approximately equal to the number of points used for shape analysis.
After generating line profiles, the integrated intensity Itot under the Gaussian was determined. Itot was approximated as 1.065 * Imax * FWHM, where Imax is the maximum fluorescent intensity, FWHM is the full width at half maximum of the Gaussian peak in the intensity profile and 1.065 is the prefactor that appears during integration. Itot appeared more uniform over the DAAM-particle surface than the Imax, because of the non-isotropic resolution in confocal microscopy. After generating a regular grid with total intensity values, this signal was binarized, initially by using a dip in the distribution of the logarithm of pixel intensities as the threshold. Then, a regionbased active contour algorithm 40 ("snake") was used to optimize the mask. Finally, interpolation of the mask was used to determine for each particle edge coordinate if it was part of the freely exposed surface.

Reference-free estimation of both normal and shear stresses from particle shape
Force calculations were performed with custom Python software. The particle shapes were smoothed within a 1 m 2 window, and the edge of the traction-free boundary was dilated and smoothed (Supplementary Note 2). Resolving traction forces from the particle shape can be defined as a linear elasticity problem on the sphere Ω with mixed boundary conditions: the surface shape expressed by the radial function 0 ( , ), and the region Ω on which the traction should be zero. The equilibrium condition of the elastic continuum Ω in terms of the displacement field is: where and is the shear modulus and the Poisson's ratio. Given the displacement on the spherical surface, ( , ), the corresponding surface shape, i.e. the radial function ( , ), can be constructed. Hence the mixed boundary condition can be written as ( , ) = 0 ( , ), and | Ω t = 0.
To solve this problem, we use an iterative approach which minimizes a cost function of trail surface displacement ( , ) . For each trial ( , ) , the elasticity equation is solved 18 , and the corresponding stress fields, surface traction forces, and not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint elastic energy el are obtained. The trial surface displacement field is expanded into spherical harmonic space with a truncation of the degree number < max : where ̂ are coefficients, and is the spherical harmonic function defined in Supplementary Note 3. In this work lmax = 20, which was chosen as a good balance between lateral spatial resolution (~800 nm) and computational time (~0.007s per traction evaluation). The traction is also expressed in the spherical harmonic space: where ̂ are coefficients that are completely determined by ̂.
Therefore, our goal is to determine the coefficients ̂ that satisfy the mixed boundary conditions described above. Given that multiple displacement fields can give rise to the same surface shape, this problem does not have a unique solution. To infer the true traction forces, we solve the following minimization problem: The shape difference ( ; 0 ) is calculated as follows: We first generate a collection of sampling points using GLQ meshing (Supplementary Note 4) on the not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint undeformed spherical surface, whose spherical coordinates are ( , , ), = 1,2, … , , where is the total number of points. Given the displacement field, we obtain the coordinates of these points after deformation ( , , ). The deviation between the deformed surface shape from the measured shape 0 ( , ) is calculated as: where ̂ is the unit vector in radial direction, the weight function and ( ) are introduced to avoid the noise introduced by GLQ meshing and measuring (Supplementary Note 4); The residual traction magnitude on the traction-free region is calculated as: where is the total number of sampling points that is on the traction-free region.

Displacement solution based on spherical harmonics
As we are using an iterative method for solving the elasticity problem, it is essential to have fast computation of the shape ( , ) from the surface displacement field ( , ).
The spherical harmonics method was previously developed as a fast solver for linear elasticity problems with spherical interfaces 18

Confocal voxel size calibration
Differences in axial (z) versus lateral (x,y) length measurements can arise from imperfect calibration of the microscope stage, and from refractive index mismatch between the sample and the immersion oil/cover glass. Hence, accurate calibration is necessary for proper reconstruction of the particle shape. First, confocal image stacks of 10 m yellowgreen fluorescent polystyrene beads (Polysciences, 18142) were recorded. Then, maximum intensity axial slices through individual particles were selected, and smoothed with a Gaussian filter (with a standard deviation of 3 pixels). Subsequently, edges of the homogeneously stained microspheres were detected using the 2D Sobel operator. The zlocation of the equator, and the diameter of the particle, was determined by finding the horizontal line were the distance between the edges (at opposite sides of the particle) was maximal. Gaussian fits to the edge profiles were made to accurately determine the particle diameter. Next, the edge at the apex of the lower hemisphere of the particle (the "south pole") was localized using a similar approach. This difference in z-position of the equatorial plane with the south pole was used to estimate the apparent particle radius in zdirection. This approach, were only the lower hemisphere is used, prevents these measurements to be affected by the lensing effect caused by the particle (Fig. 1i,   Supplementary Fig. 4). Finally, the ratio of the particle radius measured in xy and the not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint apparent particle radius in z was determined for multiple particles and average. In PBS, we found an axial elongation of 1.100  0.003 (standard error of the mean (s.e.m.), n = 34). In vectashield (VS) we did not perform such a correction, since particles did not appear elongated in the axial direction (Fig. 2).

Point spread function measurements
For measurements of point spread functions (PSFs), red or yellow-green 100 nm fluorescent beads (ThermoFisher Sci., F8801/F8803) were embedded in flat adherent pAAm gels. Such hydrogels were prepared similarly to described previously 42,43   not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint murine cells were exposed to DAAM-particles for 5 minutes and subsequently fixed. After fixation, the freely exposed particle surface was immunostained with Alexa-647 secondary antibody. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint the surface from the volumetric centroid. Middle, mean curvature of the particle surface. The colormap is centered around the median of the mean curvature (~1/R, with R the particle radius). Right, intensity signal (integrated along the radial direction) from immunostaining of the particle surface not in contact with the not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint cell. Top rows, visualization of particle surfaces from 3 different viewpoints. Bottom rows, Equirectangular map projections (standard parallel taken at latitude 0) showing the full particle surface, although necessarily distorted (most strongly around the polar regions). In these maps, the base of the phagocytic cup is visualized at latitude 0 and longitude -/2. Stars (white or black) mark the base of the phagocytic cups (as determined from the binarized secondary antibody signal), and dashed lines (white or black) mark the outlines of the phagocytic cups. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted September 30, 2018. . https://doi.org/10.1101/431221 doi: bioRxiv preprint not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.