Dynamic traction force measurements of migrating immune cells in 3D matrices

Immune cells such as natural killer (NK) cells migrate with high speeds of several µm/min through dense tissue, but the traction forces are unknown. We present a method to measure dynamic traction forces of fast migrating cells in non-linear biopolymer matrices. The method accounts for the mechanical non-linearity of the 3D tissue matrix and can be applied to time series of confocal or bright-field image stacks. The method is highly sensitive over a large range of forces and object sizes, from ∼1 nN for axon growth cones up to ∼10 µN for mouse intestinal organoids. We find that NK cells display bursts of large traction forces that increase with matrix stiffness and facilitate migration through tight constrictions.

cell behavior or cause photo-damage. All current methods 39 require knowledge of the force-free reference configuration 40 of the matrix, which is obtained either by detaching the cell 41 from the matrix, by the addition of drugs that suppress cell 42 forces (e.g. cytochalasin D), or by killing the cell, e.g., with 43 high-intensity laser light (Fig. 1b top). Such methods there-44 fore pose limitations for measuring traction forces of rapidly 45 moving cells, such as immune cells, or for investigating mul-46 tiple cells within a culture dish.

48
Here, we overcome these limitations by a locally adjusted 49 regularization approach, and by obtaining the force-free ref-50 erence configuration of the matrix from time-lapse images. 51 Our method builds on a previously developed non-linear 52 semi-affine network optimizer (Saeno) (12). This method it-53 eratively minimizes the discrepancy between measured and 54 reconstructed matrix deformations (Fig. 1a middle) by ad-55 justing point forces that can appear anywhere in the volume 56 but are regularized to favor few, large forces that are typically 57 located near the cell body. Cell forces are then summed over 58 a user-defined sphere around the cell (12,19). One advantage 59 of this approach is that the exact position of the cell surface 60 does not need to be known. A disadvantage is that balanc-61 ing forces from sources outside the image stack can occur 62 throughout the reconstructed force field. 63 In our new Python-based implementation (Saenopy), we reg-64 ularize the reconstructed forces only within the reconstruc-65 tion volume (which can be of the same size as the image 66 stack, or larger) but not at its boundaries. This ensures that 67 balancing forces from sources outside the image stack appear 68 only at the surface of the reconstruction volume. Therefore, 69 we do not need to limit the reconstructed cell forces to be 70 within a defined distance around the cell (12, 19) ( Fig. 1a 71 middle) and instead can integrate them over the entire recon-72 struction volume, excluding the surface. This allows the size 73 of the imaged volume to be greatly reduced without losing 74 accuracy (Fig. 1c). It is also no longer required that force-75 bioRχiv | August 23, 2023 | 1-35 10 0 10 1 10 2 Simulated Contractility (nN) Body force microscopy with global regularization reconstructs forces without knowledge of the cell surface. Forces within a defined distance (rmax) around the force epicenter are considered cell-generated forces (purple), forces beyond rmax (dotted gray line) are considered balancing forces (blue) that also account for noise, drift, or contractile cells outside the field of view. Bottom: Saenopy is a body force microscopy method where forces are regularized throughout the image stack except for the boundary regions. This ensures that balancing forces appear only at the boundaries. b, Schematic of two methods for measuring the force-free matrix configuration. Top: Cell forces are relaxed using actin-depolymerizing or myosin-inhibiting drugs. 3D matrix deformations are measured from two image stacks recorded before and after drug treatment. Bottom: The force-free matrix configuration is computed from the cumulative matrix displacements between consecutive image stacks of a time series, median-averaged over the observation period. Dynamic 3D matrix deformations are obtained from the cumulative matrix displacements relative to the median average. c, Influence of image stack size and reconstruction volume on force reconstruction accuracy. The cell is simulated by a force dipole with 20 nN and a distance of 30 µm. The image stack size (cube length) is increased between (20 µm-150 µm) 3 . The reconstruction volume (mesh) has either the same dimensions as the image stack (yellow) or has a constant size of (150 µm) 3 (blue). Force reconstruction is performed with Saenopy regularization (yellow and blue lines) or with global regularization (rmax = 50 µm, magenta line). d, Accuracy of force reconstruction for simulated dipoles (upper inset) and quadrupoles (lower inset) for different contractilities. Dipoles or quadrupoles are simulated as pairs of monopoles either 15 µm or 30 µm apart. Contractility is the sum of the force monopoles. Dashed gray line shows line of identity. e, Influence of noise on force reconstruction (top) and strain energy (bottom) for simulated dipoles (contractility of 20 nN) or quadrupoles (contractility of 40 nN). Maximum matrix deformations are around 3 µm in both cases. Noise is added as Gaussian noise with specified standard deviation σ to the deformation field (upper right images for three selected noise values) and is expressed as absolute sd noise σ or percentage of sd noise relative to the maximum matrix deformations σ dmax . Lower right images show reconstructed deformation fields. Reconstructed contractility and strain energy are expressed relative to the input values.
The ability to obtain the force-free reference state without 94 killing or relaxing the cells, combined with the ability to 95 recover the traction field from an incomplete deformation 96 field, allows us to measure the dynamic force generation 97 of natural killer cells (NK92 cell line) during migration in 98 1.2 mg/ml collagen gels (shear modulus in the linear range 99 of 100 Pa (SI Fig. 1-3), average pore size 4.4 µm (SI Fig. 100  4)). We record one image stack (123×123×123 µm, voxel-101 size 0.24×0.24×1 µm, acquisition time of 10 seconds using a 102 resonance scanner (8000 Hz) and a galvo stage) every minute 103 for 23 min. From these 23 image stack pairs, we compute 104 the differential deformations using particle image velocime-105 try (PIV) (20), integrate them over time, and estimate the 106 force-free matrix configuration (undeformed reference state) 107 from the median deformation (Fig. 1b). Cell forces at each 108 time point are then computed from the current matrix defor-109 Contractile phases during natural killer cell migration. a, Matrix deformations (relative to the dynamic reference state) around a natural killer cell migrating in a collagen type I hydrogel at 3 time points selected from a 23 min time series. b, Corresponding maximum-intensity projected confocal reflection image overlaid with the projected maximum matrix deformation field. c, Cell speed (blue), cell aspect ratio (1: round, >1: elongated, orange) and cell contractility (magenta) over a 23 min period. Dashed gray lines indicate the time points shown in a,b. d, Cells treated with Y-27632 (10 µM) or Blebbistatin (3 µM) show a decrease in maximum contractility (top left, mean±se), decrease in speed (bottom left, mean±se), and a decrease in the fraction of contractile cells (defined as contractility > 5 nN). Measurements of individual cells are shown as points, colors indicate independent experiments from different days. ** indicates p<0.01 for two-sided t-test. e, Cross-correlation between time courses of speed and force (left), speed and aspect ratio (middle), and aspect ratio and force (right). Signals were shifted relative to each other by ±5 min. The physical meaning of a positive correlation for different time points is indicated in the diagrams. Data points are mean values from 52 cells, shaded areas indicate one standard deviation determined by bootstrapping. f, Differential interference contrast (DIC) images show morphological changes of a NK cell during migration through collagen pores. Arrows indicate the local matrix deformation between consecutive frames. g, Cell contractility (maximum value within 23 min measurement period) and matrix deformation (99% percentile value of the absolute matrix deformation vectors is calculated per timestep. Then we choose the maximum value within 23 min measurement period) of NK92 cells embedded in collagen gels of different stiffness and collagen batches (A,C, see SI Fig. 1) (mean±se for 44-56 cells from three independent experiments; ** indicates p<0.01 for two-sided t-test, SI Fig. 6). The storage modulus of different collagen gels in the linear regime is derived from shear rheological measurements with 1% strain amplitude and and a frequency of 0.02 Hz (SI Fig. 3). mations relative to the undeformed reference state. Thus, 110 this method allows for dynamic deformation measurements 111 of highly motile cells with rapidly fluctuating traction forces.

112
From the cell positions and cell shapes, we also extract cell 113 velocity and aspect ratio (Fig. 2a-c).

114
Around 78 % of NK cells that are elongated at t=0 migrate 115 by more than one cell diameter (13 µm) during the obser- show bursts of substantial contractile forces between 5-50 nN 120 (Fig. 2d, SI Video 20-26), lasting typically 3-5 min, fol-121 lowed by non-contractile phases (< 5 nN). Cells migrate both 122 during contractile and non-contractile phases, but contrac-123 tile phases are significantly (p < 0.01) correlated with higher 124 speed and higher cell aspect ratio (Fig. 2c, e). Accordingly, 125 speed and aspect ratio are also positively correlated (Fig. 2e). 126 Myosin inhibition with blebbistatin and Rho-kinase inhibitor 127 Y-27632 markedly reduce the magnitude and frequency of 128 force bursts (Fig. 2d) and also reduce the migration speed 129 (Fig. 2d)  suggest that strong integrin-mediated adhesions may also be 157 involved in immune cell migration through tight constric-158 tions. This finding may have broad biological implications 159 and adds to our current understanding of how immune cells 160 migrate through tissue for their homing towards target sites 161 in inflammation or cancer. 162 We next apply our method to measure traction forces of 163 growth cones from retinal ganglion cell axons of African 164 clawed frogs (Xenopus laevis). From 2D TFM measure-165 ments, it is known that forces of axonal growth cones are 166 small (1-10 nN) (29-31) and highly dynamic (31, 32). In 167 3D collagen and PEG gels, axons growth cones have been 168 reported to generate no measureable traction forces (33). Be-169 cause axons are typically much larger than the 145 µm size of 170 our image stacks, traction forces can be partially unbalanced. 171 Despite these challenges, our method resolves traction forces 172 of axonal growth cones in 3D of around ∼3 nN. As seen in 173 immune cells, these forces act only intermittently, and longer 174 periods of axon growth can occur without detectable traction 175 forces (SI Fig. 7, SI Video 27).

176
In contrast to weakly contractile axon growth cones, fibrob-177 lasts are highly contractile and can induce extensive matrix 178 deformations that play an important role in the progression 179 of fibrosis or cancer (34-36). Even larger matrix deformations can be observed around spheroids or organoids consist- This description gives rise to a semi-affine behavior: Macro-232 scopic material deformations cause non-affine deformations 233 of the individual fibers, e.g. some fibers may buckle while 234 others are stretched, depending on their orientation. Using a 235 mean-field approach, we approximate the material behavior 236 by averaging over many isotropically oriented fibers so that 237 all spatial directions contribute equally.

238
-1<λ<0 0<λ<λ s λ s <λ Specifically, the energy function w(λ) is used to calculate the 239 energy density W stored in the material for a given deforma-240 tion gradient F, averaged over the full solid angle Ω (12, 44): 241 The strain λ of a unit vector ⃗ e r (Ω) (representing a fiber), λ = 242 |F · ⃗ e r (Ω)| − 1, is used to calculate the energy ω(λ) stored 243 in this spatial direction. The strain is averaged over the full 244 solid angle Ω, approximated by a finite set (typically 300) 245 of homogeneously and isotropically distributed unit vectors 246 ⃗ e r (Ω). 247 Finite element model. To describe the material behavior in 248 the case of an inhomogenous deformation field, we tessellate 249 the material volume into a mesh of linear tetrahedral elements 250 T . The deformation field is then modeled as deformations ⃗ U 251 of the nodes of these elements. The total energy E( ⃗ U ) of the 252 mesh is the sum over the energy density W of all tetrahedral 253 elements T multiplied by their volume V T : From the total energy E( ⃗ U ), we calculate the stiffness matrix 255 K( ⃗ U ) and the nodal forces ⃗ f ( ⃗ U ) (12, 44): Index i represents the nodes of the finite element mesh, with 258 index j representing the adjacent nodes.

259
Mechanical characterization of biopolymer networks. 260 We investigate the influence of three different collagen

299
For each stretch increment, we measure the gel thickness by 300 focusing on the top and bottom surface of the gel (highest and 301 lowest layer of beads) and correct for the refractive index of 302 water.

303
For collagen batch D (used only for axon growth cone mea-304 surements), the stress-strain relationship is determined using 305 a rheometer equipped with a 25 mm diameter cone-plate ge-306 ometry (MCR302e, AntonPaar, Graz). The polymerisation 307 of 50 µl of hydrogel between the rheometer plates is moni-308 tored for 90 minutes at 20°C using an oscillatory strain of 309 0.5%. After polymerization, a frequency sweep (0.1-100 310 rad/s) with 0.5% strain is performed, followed by a strain 311 ramp at a constant strain rate of 0.5%/s (SI Fig. 2).

312
Because a cone-plate rheometer or uniaxial stretching device 313 may not be available in every laboratory, and hence the cor-314 rect material parameter for a given collagen batch may be un-315 known, we also quantify the effects of using erroneous ma-316 terial parameters on cellular force reconstruction. (SI Fig. 317 16-17). Furthermore, we explore the possibility to derive the 318 linear stiffness parameter k 0 of the material model, which 319 most sensitively affects force reconstruction, from oscillatory 320 cone-plate rheology measurements in the linear range (SI Fig. 321 3).

322
Estimating material parameters from rheological mea-323 surements. As the material model is based on the mechan-324 ical behavior of microscopic biopolymer fibers, the model 325 parameters cannot be directly extracted from macroscopic 326 rheological measurements. Instead, we simulate the macro-327 rheological experiments using different model parameters 328 until a best-fit is achieved (SI Fig. 1). Data from uniaxial 329 stretching are required to fit the buckling coefficient d 0 of the 330 material. For the other parameters, either data from a shear 331 rheometery or an extensional rheometry experiment are suf-332 ficient as long as the strain range exceeds the linear range of 333 the material (43).

334
Collagen gel preparation. Collagen type 1 hydrogels are 335 prepared from acid-dissolved rat tail (R) and bovine skin (G1) 336 collagen (Matrix Bioscience, Mörlenbach). Both collagen 337 types are mixed ice-cold at a mass ratio of 1:2 and are dis-338 solved in a dilution medium containing 1 vol part NaHCO 3 , 339 1 vol part 10× DMEM and 8 vol parts H 2 O, adjusted to pH 9 340 using NaOH. The final collagen concentrations are adjusted 341 to 0.6 mg/ml, 1.2 mg/ml or 2.4 mg/ml using dilution medium. 342 The the final mixture is then adjusted to pH 9 using NaOH. 343 Collagen gels are polymerized for 1 h at 5% CO 2 and 37°C.

344
For axon growth cone measurements, collagen type 1 hy-345 drogels are prepared from acid-dissolved bovine collagen 346 (TeloCol-10, Advanced Biomatrix, Carlsbad    culated between a deformed state and a force-free reference 527 state (Fig. 1b top). The reference image stack is typically 528 recorded after drug-induced force relaxation using e.g. high 529 concentrations of cytochalasin D (12,19,21). For fast mov-530 ing cells, Saenopy offers a second option for obtaining the 531 reference state from a sufficiently long time series of image 532 stacks (Fig. 1b bottom): First, the change in matrix deforma-533 tions between consecutive image stacks, which correspond to 534 the time derivative of the matrix deformations, are calculated 535 at each voxel position using OpenPIV (20). We then calcu-536 late the cumulative sum of the time derivative at each voxel, 537 which corresponds to the absolute matrix deformation apart 538 from an offset. To remove the offset, we subtract the median 539 deformation at each voxel position, assuming that the median 540 matrix deformations around a fast moving cell tend towards 541 zero, when observed over a prolonged time period. For the 542 NK92 cells in this study, we employ the second option. By 543 comparing a region of the collagen matrix before and after an 544 NK cell has passed through, we confirm that these cells do 545 not permanently remodel the matrix.

546
For the OpenPIV algorithm, we use a window-size of 12 µm 547 with an overlap of 60% (corresponding to an element size 548 of 4.8 µm). Deformation vectors with low confidence (de-549 fined by OpenPIV for a ratio < 1.3 between the global maxi-550 mum of the image cross correlation divided by the second-551 highest local maximum value) are ignored. The accuracy 552 with which we can measure matrix deformations with our 553 setup is estimated by imaging collagen gels without cells. We 554 record image stacks at 28 non-overlapping positions in the gel 555 (each 123×123×123 µm), each over a period of 23 min with 556 a time interval of 1 min and a voxel-size of 0.24×0.24×1 µm 557 as in the cell experiments. We then compute the deforma-558 tion field at every voxel position as described above, and 559 the accuracy (noise level) is calculated as the standard devia-560 tion σ of all deformation values, separately evaluated in x,y, 561 and z-direction. We obtain values of σ x =41 nm, σ y =42 nm, 562 σ z =99 nm.

563
Cell tracking. Cell shape and position at each timepoint are 564 extracted from intensity-projected bright-field image stacks. 565 We compute the local entropy of the band-pass filtered im-566 ages, subtract the background using regional maxima filter-567 ing, binarize the image using Otsu-thresholding (49), and ap-568 ply binary morphological operations to remove small objects 569 and holes. From the binarized object (which is the segmented 570 cell), we calculate the area, the x,y position (center of mass), 571 and the aspect ratio defined here as the major axis divided by 572 the minor axis of an ellipse fitted to the cell shape. If a cell 573 reaches the edges of the image, this cell is not tracked further. 574 Correlation analysis. We cross-correlate the time develop-575 ment of cellular contractility, cell speed, and cell aspect ra-576 tio for all NK92 cells under control conditions (Fig. 2d,f). 577 Since cell velocity and cell forces are calculated from the 578 difference in cell position and matrix deformations between 579 two consecutive image stacks, whereas the aspect ratio is ob-580 tained for each single image stack, we linearly interpolate 581 the aspect ratio between two consecutive image stacks before 582 cross-correlation. We then calculate the Spearman rank cor-583 relation coefficient between speed and contractility, between speed and aspect ratio, and between contractility and aspect 585 ratio, for time shifts from -7 to +7 min.
The force field ⃗ f is computed from the simulated deformation

607
As f ( ⃗ U ) is non-linear, L( ⃗ U ) cannot be minimized easily. 608 We apply a first-order Taylor series expansion of the dis-609 placement ⃗ U ( ⃗ U + ∆ ⃗ U ), using the stiffness matrix K (Eq. 4)
The ∆ ⃗ U that minimizes this equation satisfies the following 612 normal equation (51).

This linear equation (of the form
ing the conjugate gradient method to obtain a value for ∆ ⃗ U .

615
Because of the pronounced non-linearity of the problem, for 616 the next iteration cycle we update the simulated deformation 617 field only by a fraction (stepper) of ∆ ⃗ U (typically 0.33): With the new displacement ⃗ U ′ , the stiffness matrix K( ⃗ U ′ )

619
(Eq. 4), the nodal forces ⃗ f ( ⃗ U ′ ) (Eq. 5), and the weight ma- . 7) are updated, and the linear Taylor ex-621 pansion (Eq. 9) is solved again. This procedure is iterated 622 until the total energy E( ⃗ U ) (Eq. 3) of the system converges. 623 Convergence is reached when the standard deviation of E( ⃗ U ) 624 divided by the mean of E( ⃗ U ) (coefficient of variation) for the 625 last 6 iterations is below a convergence criterion τ (typically 626 0.01).

627
From the resulting force vectors ⃗ f i at all nodes i, we compute 628 the coordinates of the force epicenter ⃗ c. To find it, we mini-629 mize the magnitude Q of the cross product of the force field 630 ⃗ f i with the vectors from the nodes ⃗ r i pointing to the epicenter 631 coordinates ⃗ c (12, 44): We determine the cellular contractility as the sum of the pro-633 jected forces in the direction of the force epicenter (12, 44): In this study, we perform the force reconstruction of migrat-635 ing NK92 cells using the regularization parameter α = 10 10 636 (Eq. 7), which maximizes contractility and is stable against 637 noise both for the experimental setup and for simulated cells 638 (SI Fig. 18). We interpolate the 3D matrix deformations onto 639 a cubic grid with a grid-size of 4 µm (in the case of single 640 cells, or larger as specified in Methods), and tesselate cubes 641 into 6 tetrahedra (SI Fig. 19).

642
Boundary forces. Because the boundary of the simulated, 643 finite-sized volume is fixed, the reconstructed force field con-644 tains cell-generated forces from the bulk, as well as dis-645 tributed counter-forces from the surrounding matrix that ap-646 pear at the volume boundary. The individual force vectors of 647 the counter-forces are small compared to cell forces. Hence, 648 the small-force penalization scheme (Eq. 7, first and second 649 case) would bundle them into fewer, larger force vectors lo-650 cated in the bulk where they could interfere with cell forces. 651 To avoid this, we set the regularization weights A ii to zero 652 for all tetrahedral mesh nodes at the boundary (surface) of 653 the simulated volume (Eq. 7, third case). Hence, we exclude 654 the surface of the mesh from the regularization cost function 655 so that virtually all counter-forces only appear at the volume 656 boundary and not in the bulk (Fig. 1a). Saenopy can therefore 657 directly sum all the forces in the bulk of the mesh as they are 658 contractile cell forces, without the need to arbitrarily define a 659 maximum radius around the cell center for the summation of 660 forces. This improves robustness and accuracy especially for 661 small stack sizes (Fig. 1c,e).

662
Simulations of force-generating dipoles and 663 quadrupoles. We simulate the 3D deformation fields 664 around force-generating dipoles and quadrupoles in a 665 1.2 mg/ml collagen matrix (Batch A; SI Fig. 1) with the 666 iterative scheme described above, with the following mod-667 ifications: The displacements at the mesh boundary are set 668 to zero, the forces at the dipole or quadrupole points are 669 given, and the forces in the bulk are set to zero. instead of the Eq. 9, we solve the following equation: and update the deformation field ⃗ U after each iteration ac- that an elastic material model is appropriate. Potential prob-723 lems arise only if the force-free matrix configuration is cal-724 culated from images taken at the beginning of the experiment 725 (or from the median deformations) so that non-elastic matrix 726 deformations, if present, can contribute to the total deforma-727 tions. In this case, a purely elastic model overestimates the 728 forces. For highly motile immune cells, however, we do not 729 observe permanent deformations of the collagen matrix after 730 the cell has migrated through it.    The collagen fiber structure of two different collagen batches (a-c: Batch A; e-g: Batch C) is imaged using confocal 998 reflection microscopy for three different collagen concentrations (imaged volume of 160x160x200 µm with voxel-sizes of 999 0.314x0.314x0.642 µm). Images show a single slice. 3D pore diameters are computed from the covering radius transform  a, Benchmark tests comparing the accuracy of force reconstruction using linear Saenopy (Boundary Method, linear elastic 1033 material model without parameters λ s , d 0 and d s , left), Saeno (Maximum Distance, r max =50 µm, middle), and non-linear 1034 Saenopy (Boundary Method, right), for artificial dipoles (blue) and quadrupoles (orange). Force reconstruction with linear 1035 Saenopy tends to overestimate total contractility. Force reconstruction with Saeno (maximum distance) tends to underestimate 1036 total contractility for high forces. Force reconstruction with non-linear Saenopy gives accurate contractility estimates over the 1037 entire simulated range. Images (right) show simulated deformation fields around a dipole (200 nN) and a quadrupole (400 nN). 1038 b, Reconstructed contractility for dipoles (circles) and quadrupoles (crosses) with low (6 and 12 nN), medium (60 and 120 nN) 1039 and high contractility (200 and 400 nN) versus maximum distance parameter r max for Saeno (blue). For comparison, results 1040 from the Saenopy Boundary Method (orange) and ground-truth values (black) are shown as lines. Compared to the boundary 1041 method, Saeno shows larger deviations for all r max values especially at high contractilities. c, Reconstruction of unbalanced 1042 forces from monopoles and tripoles (a single point force is removed from a dipole or quadrupole) using Saenopy (boundary 1043 method) and Saeno (maximum distance, r max =50 µm). For larger forces, Saenopy (boundary method) outperforms Saeno.

1048
A representative hepatic stellate cell embedded in collagen (1.2 mg/ml, Batch C) is imaged with confocal reflection microscopy. 1049 Image stacks are cropped in x,y,z-direction by different amounts. From the cropped image stacks, we calculate matrix deforma-1050 tions and reconstruct forces. We compare Saenopy (Boundary Method), linear Saenopy, Saeno (Maximum Distance method, 1051 r max set to 75 µm), and TFM Lab software (9). For Saenopy and Saeno, we reconstruct the cell forces in a volume of (370 µm)³. 1052 For reconstruction with TFM Lab, we chose an elastic modulus of 1010 Pa (Batch C, 1) and a Poisson's ratio of 0.25, and ad-1053 ditionally provide information of the cell surface based on fluorescence images (calcein). We choose a comparable number of 1054 finite elements for all four methods. From the reconstructed force vectors, we calculate the sum of absolute forces (top left). Lin-1055 ear elastic material assumption (linear Saeno, and TFM Lab) results in large forces (to compensate for long-range-propagated 1056 matrix deformations that are not explained by linear material properties). Normed total forces (top right) are normalized to the 1057 respective values obtained from uncropped stacks. Reconstructed forces decrease with cropping in the case of TFM Lab and 1058 linear Saeno, and are more stable in the case of Saenopy.     is more pronounced. In general, we find that the contractility and strain energy for differently contractile quadrupoles change 1115 proportionally, which implies that relative changes in contractility can be correctly inferred even with the incorrect material 1116 model, except for material parameters that cause extremely non-linear material behavior (very low d 0 and d s ). a, Cubic meshes are tessellated into either 5 tetrahedra (not of equal volume) or 6 tetrahedra (of equal volume). b, Reconstructed 1152 contractility (for 5 tetrahedra sub-tessellation) for differently contractile dipoles (two force monopoles 15 µm or 30 µm apart) 1153 and quadrupoles (force monopoles at the vertices of a regular tetrahedron with an edge length of 15 µm or 30 µm). For total 1154 contractilities above 50 nN, the reconstructed forces become increasingly underestimated. c, same as in (b) for 6 tetrahedra 1155 sub-tessellation. In contrast to a 5 tetrahedra sub-tessellation, forces are correctly reconstructed also for higher contractilities 1156 (c).
Supplementary Information 20: Video of an NK92 cell migrating through collagen (bright-field