Traction microscopy to identify force modulation in subresolution adhesions.

We present a reconstruction algorithm that resolves cellular tractions in diffraction-limited nascent adhesions (NAs). The enabling method is the introduction of sparsity regularization to the solution of the inverse problem, which suppresses noise without underestimating traction magnitude. We show that NAs transmit a distinguishable amount of traction and that NA maturation depends on traction growth rate. A software package implementing this numerical approach is provided.

Mechanical forces, transmitted via integrin-based adhesions to the extracellular matrix (ECM)[1], play an important role in integrin signaling, environmental sensing, and directed migration[2, 3]. These forces, or tractions, are not only the result of inside-out coupling of cytoskeletal activities (e.g. actomyosin contraction or propulsion of the growing actin meshwork) to the ECM but are also the effector for adhesion maturation[4], as studied extensively at the level of large focal adhesions (FAs)[5]. How much force adhesions transmit during their nascent state and whether the fate of nascent adhesions (NAs) also depends on force transduction has remained unknown because of the technical difficulties in resolving traction in adhesions with a size below optical diffraction.

Traction microscopy (TM) derives the traction exerted by cells onto their environment from the displacement of fluorescent beads embedded in, or coated on, deformable gel substrates[6, 7]. This requires a solution to the inherently ill-posed inversion of the deformation field into the generating traction field, i.e. noise in the deformation field can generate out-of-bound traction values. A well-established remedy is regularization, which stabilizes the reconstruction by constraining spatial variation in the traction field[8, 9].

To examine the effect of regularization in TM, we simulated bead displacements in a virtual gel substrate that is exposed to multiple traction impact regions of varying sizes and magnitudes (Fig. 1a, Supplementary Fig. 1). Accounting for physical disturbances in the substrate like gel swelling (Supplementary Fig. 2), we added to the designed traction field white noise up to 100 Pa (Fig. 1a, inset) before constructing bead images of the deformed substrate that were further subjected to pixel-by-pixel white noise. Traction reconstruction with regularization relies on minimizing the sum of the residual norm, i.e. the difference between predicted and measured deformation fields, and the solution norm of the traction field (Fig. 1b). Two parameters have to be tuned for this minimization problem: the regularization parameter λ controlling the relative importance of the two norms[10] and the type of solution norm (Online Methods). For the latter, TM algorithms have consistently used the Euclidian norm, or L2 norm, which is referred to as Tikhonov regularization or L2 regularization[9-13]. However, in many L2-regularization-based TM applications λ is selected subjectively and without explicit documentation of criteria[1, 14, 15], which prevents consistency in processed data. Our best guess is that λ is chosen by visual inspection to avoid reconstruction error in the traction impact regions (λ), or to avoid notable traction in the non-adhering background (λ). Importantly, neither λ nor λ are known in real experiments, making this choice arbitrary. Moreover, dependent on the choice, noise in the background is amplified or the traction is underestimated (Fig. 1b).

Figure 1

L-curve analysis for L2- and L1-regularization. (a) Simulated input traction field. Inset: rescaled traction map displaying small traction impact regions in the dashed window. The simulated traction field includes a maximum 100 Pa of white noise to reflect physical disturbances in the substrate like local gel swelling. The orientation of the traction field is downward, parallel to the y-axis. (b) (Top) Minimization problem in inverse problem; the regularization parameter λ determines the weight of the regularization semi norm against the residual norm; p denotes the type of norm (e.g. L2 norm if p=2, or L1 norm if p=1). (Bottom) traction maps resulting from qualitative choices λ and λ applied to L2-regularization (p=2). (c) L-curve for L2-regularization. λ is determined by the second derivative of the semi norm with respect to the residual norm. The regularization parameter λ in L2-regularization is chosen as the first inflection point in the L-curve for which λ < λ (Supplementary Fig. 3). (d,e) Traction maps reconstructed under L2, λ (d), L2, λ (e) regularization. Insets: rescaled traction maps of small traction impact regions in the dashed windows; red arrowheads indicate noise spikes in comparison to true traction signals (yellow arrowhead) (f) L-curve for L1-regularization. The regularization parameter λ in L1-regularization is chosen as the first inflection point for which λ > λ (Supplementary Fig. 3). (g,h) Traction maps reconstructed under L1, λ (g) and L1, λ (h) regularization. Scale bar: 10 µm. Scale bar in inset: 2 µm. (i) Traction recovery ratio, i.e. magnitude ratio of reconstructed to original tractions in small traction impact regions (NAs) and in large traction impact regions (FAs). n = 9 for NAs, n = 26 for FAs; error bars, s.d. (j) Average traction stress background outside traction impact regions. n = 229,255 pixels in non-adhesion area; error bars, s.e.m. Dotted line indicates the noise level in original traction field. ***P < 0.001; as determined by the student t-test. N.S., not significant. (k) Detectability of traction at NAs with a diameter of 6 pixel and load of transmit 600 Pa. Five replicate simulations; error bar shows s.e.m. (l) Regimes of traction detectability as a function of adhesion diameter (72 nm/pixel) and traction magnitude for reconstructions under L2, λ, L2, λ, L1, λ and L1, λ. Five replicate simulations for each pair of traction magnitude and adhesion diameter.

One criterion for objective selection of λ is the L-curve[10], which relates the residual norm to the solution norm. The best λ-value is determined by the position along the curve where the combined differential between the norms is maximal (referred to as the L-corner). In the present simulation, the value of λ was two-orders of magnitude larger than that of λ (Fig. 1c). Accordingly, the traction field determined with λ showed substantial underestimation of the traction magnitude compared to the simulated reference traction field (Fig. 1d). To improve on L-corner criterion, we defined λ as the inflection point in the L-curve smaller than λ (Supplementary Fig. 3). Similar to the reconstructions with λ, the reduction in regularization strength compared to λ caused noise spikes with a magnitude in the range of small impact regions (Fig. 1d). We thus concluded that L2-regularization makes it difficult to choose the right regularization parameter: parameters derived from either the L-corner or minimization of background tractions lead to substantial under-estimation of the stress field, as discussed in a previous study[9], whereas a parameter derived from minimization of the error in adhesion tractions increases the background level, which obscures tractions in small adhesions.

As an alternative to L2-regularization, one can use an L1 norm[16]. A key feature of L1-regularization is that it forces the solution to be sparse[17], which could be beneficial to TM as the majority of the traction field is at the background level with a few sparsely located traction impacts at discrete adhesions[16, 18]. To test this, we reconstructed the traction field using a range of regularization parameters and examined the L-curve and reconstruction accuracy (Fig. 1f–h). Using λ for L1-regularization, tractions at both small and large traction impact regions were restored to a level much closer to the level of the simulated traction field (Fig. 1g). However, due to the disturbances in the simulated traction field, the reconstructed traction field showed noise spikes whose level was comparable to traction at small traction impacts (Fig. 1g, inset, red vs. yellow arrowheads). These noise spikes were suppressed when λ was used without affecting the magnitudes at large traction impacts (Fig. 1h,i). This background suppression increased the likelihood of detecting weak tractions in small impact regions (Fig. 1h, inset, yellow arrowhead, Fig. 1j) and led to manifold improved traction detectability (Fig. 1k). This conclusion was corroborated by simulations of single impact regions with variable traction magnitude and size (Online Methods), showing that L1-regularization affords a substantially lower detection limit for stress magnitude and traction impact region size compared to that found using L2-regularization (Fig. 1l).

Using L1-regularization with the choice of λ as the regularization parameter yielded a best-matching traction field that preserved tractions at NAs (Supplementary Fig. 4) even compared to traction maps that minimized traction errors on either traction impact regions or background (λ or λ, respectively). This indicates that L1-regularization provides an objective way to choose the best regularization parameter in real-world experiments where λ or λ cannot be defined. Simulations with varying adhesion density and size also showed that λ is detected robustly regardless of adhesion density (Supplementary Fig. 5), whereas L2-regularization showed inconsistent L-curve shapes (Supplementary Fig. 6). Spatial resolution analysis also indicated that an L1-regularized traction reconstruction can distinguish two local traction maxima if they are separated by at least 8 pixels = 0.56 µm whereas L2-regularized traction reconstructions resolved local maxima only at distances greater than 11 pixels = 0.79 µm (Supplementary Figs. 7–9). Together, these analyses on simulated traction fields established L1-regularization as the far superior regularization scheme for TM and removed subjectivity in the selection of the regularization parameter.

To compare the performance of L1- versus L2-regularization in reconstructing tractions in live cell imaging experiments, we filmed migrating PtK1 epithelial cells expressing TMR-paxillin on 8 kPa TIRF-compatible silicone gels with beads covalently-bonded to the gel tops[7]. As in synthetic experiments, beads were tracked using a subpixel correlation via image interpolation (SCII), which found to be as accurate as the continuous window-shift method[19] but is not iterative (Online Methods). To test the spatial resolution of the reconstructed traction, we used a single frame that contained a partial view of the protruding cell edge and compared the localization of adhesions as visualized by paxillin fluorescence (Fig. 2a) to the traction field computed based on L2-λ (Fig. 2b), L2-λ (Fig. 2c) and L1-λ (Fig. 2d). The overall traction magnitude obtained from L2-λ was lower than the one obtained from L2-λ or L1-λ. At adhesions, the traction magnitudes obtained from L2-λ and L1-λ were similar. However, the L1-based reconstruction again displayed markedly less noise in the background and especially outside the cell footprint (Fig. 2c,d).

Figure 2

Colocalization of tractions with adhesions. (a) Snap-shot of HaloTag-TMR-paxillin at the leading edge of a protruding PtK1 cell. Scale bar: 2 µm. (b–d) Traction maps reconstructed under L2, λ (b), L2, λ (c), and L1, λ (d) regularization. (e–i) Colocalization analysis in FA-rich region (box 1 in a–d). (j–n) Colocalization analysis in NA-rich region (box 2 in a–d). (e,f,j,k) Inverted HaloTag-TMR-paxillin image without (e,j), and with (f,k) tracked NAs (red circles), FCs (orange outlines), and FAs (black outlines). Scale bar: 2 µm. See also Supplementary Video. (g–i, l–n) Traction maps (overlaid by tracked adhesion regions) reconstructed under L2-λ (g,l), L2-λ (h,m), and L1-λ (i,n) regularization. (o) Traction magnitude in NAs compared to the magnitude of local stress maxima outside the cell. n = 10, 27, 14 background peaks in b, c, d, respectively, n = 67 each for NAs in b, c, d; error bars, s.e.m. (p) Fraction of NAs colocalized with local traction stress maxima in L2-λ, L2-λ, and L1-λ reconstruction. n = 7,653 NAs from 100 frames; error bars, s.e.m. **P < 0.01, ***P < 0.001; as determined by the student t-test. N.S., not significant.

We examined in more detail a cell region rich in focal adhesions (FAs) (Fig. 2a, box 1 and 2e) and a region rich in NAs (Fig. 2a, box 2 and 2j). We developed an image segmentation method that used the paxillin channel to distinguish focal complexes (FCs), FAs, and diffraction-limited NAs, and then co-localized the adhesions with the traction field (Fig. 2f,k, see Online Methods). In the FA-rich region, the traction magnitude reconstructed by L2-λ was underestimated (~1200 Pa at maximum; Fig. 2g). Reconstructions with either L2-λ or L1-λ had stress magnitudes up to ~1600 Pa and showed tighter co-localization with the adhesion sites (Fig. 2h,i). In the NA-rich region, reconstructions based on L2-λ showed consistently low stress (~450 Pa at maximum) and lacked fine features in the stress distribution (Fig. 2l). The traction field from L2-λ showed a more prominent variation around the segmented NAs, with stress maxima reaching 650 Pa. However, the adjacent extracellular region displayed spikes in the background with magnitudes similar to those associated with NAs (Fig. 2m), making it difficult to distinguish bona fide traction values at NAs from random stress maxima induced by the regularization. The L1-λ reconstruction displayed a background stress level much lower than the tractions associated with NAs (Fig. 2n). A direct comparison of the traction magnitudes at the sites of NAs and the magnitudes of local traction maxima in the background outside the cell confirmed that only with a reconstruction using L1-λ NAs had traction values that were significantly higher than noise spikes (Fig. 2o). Moreover, L1-λ showed a higher fraction of NAs associated with a local maximum in the traction than the L2-λ reconstruction (~55% vs. ~30%; Fig. 2p). Simulations revealed that at the bead density used in these experiments, up to 13% of adhesions with small tractions could be missed due to the random bead distribution (Supplementary Fig. 10). These results suggest that more than one half of NAs transmitted independent tractions that remained largely undetected with L2-regularization, and that some NAs (~32%) might not be involved in traction transmission.

Exploiting the enhanced resolution of the L1 TM, we studied the traction modulation at individual adhesions by tracking each NA (Supplementary Video, Online Methods). We measured the tractions in NAs that transitioned into FCs and FAs (Fig. 3a) and in NAs that disassembled before maturation (average lifetime of ~2 minutes; Fig. 3b). In NAs that underwent the maturation process the traction stress increased over time (Fig. 3c), whereas it remained constant in NAs that did not mature (Fig. 3d). Irrespective of the fate of the NA, the stress magnitude at the time of NA appearance in the paxillin channel was significantly higher than tractions in the background or even tractions in the 1 µm vicinity outside the cell edge (Fig. 3e). Hence, tractions measured in emerging NAs are above the noise level. In fact, high tractions might be required for the assembly of NAs. A population analysis of failing vs. maturing NA tracks over the first two minutes of their lifetimes (Fig. 3f) indicated that maturing NAs undergo on average a two-fold higher rate in stress increase than failing NAs (Fig. 3g). They also started at higher stress levels (Fig. 3h), which suggested that sufficient tension is a pre-condition for adhesion maturation. It is worth noting that the higher traction growth rate of maturing NAs was also observed when traction was reconstructed under L2-λ (Supplementary Figure 11). However, these tractions had higher uncertainty compared to L1-based tractions (Fig. 2m,o). To investigate how spatially independent the maturation fates of NAs are, we performed a spatial clustering analysis. At each time point of an maturing or failing adhesion track (Fig. 3i), we considered the fate of neighboring NAs within a circular vicinity and calculated the fraction of maturing NAs within a population (Fig. 3j). This analysis revealed that the fate of maturing adhesions is coupled over a distance of 4 µm (Fig. 3k), implying that adhesions tend to mature in a collective, possibly by experiencing spatially clustered mechanical inputs e.g. from edge protrusion, actin retrograde flow or myosin activity at the edge. In contrast, the fate of failing adhesions was independent of their neighbors, regardless of distance. Taken together, the data provided by L1 TM provide solid statistical evidence that the maturation of NAs is directly coupled to the development of the transmitted traction.

Figure 3

Traction analysis in individual NAs of a migrating PtK1 cell. (a,b) Time-lapse images of paxillin (top) and traction stress map reconstructed under L1 (bottom) for representative examples of a maturing NA (a) and a failing NA (b); (see also Supplementary Video. Black solid lines in a show segmentation of FCs and FAs associated with tracked adhesion. Dotted red circles in the traction maps indicate the positions of detected adhesions. Scale bar: 1 µm. (c,d) Fluorescence intensity (black) and traction stress (orange) as a function of time in a,b, respectively. (e) Comparison of tractions in emerging NAs (n = 158), in a 1 µm-wide band outside the cell edge (n = 332), and in the entire area outside the cell edge (n = 735). Boxes extend from the 25th to 75th percentiles, with a line at the median. Whiskers extend to 1.5× IQR (interquartile range). Asterisks, outliers. (f) Time courses of traction in maturing NA tracks (green lines, n = 40) and failing NA tracks (orange lines, n = 242); thick lines represent average time course. (g,h) Comparison of rate of traction increase (g) and traction magnitude in the time point of initial appearance (h) between failing NA tracks and maturing NA tracks. Sample numbers are the same as in f. Error bars, s.e.m. (i–k) Spatial clustering analysis of failing and maturing NAs. (i) Failing and maturing adhesion tracks at a particular time point overlaid on a paxillin image. Scale bar: 2 µm. (j) The ratio of maturing adhesions to all NA tracks present at each frame was averaged over all frames and compared to the fraction of either maturing or failing adhesions present in the circular vicinity of each adhesion at the time point of appearance. (k) The ratios of maturing adhesions for the total population, for neighbors of failing NAs and for neighbors of maturing NAs in cases of varying vicinity radii. The ratio for neighbors of maturing adhesions (n = 779) becomes similar to the ratios for the total population (n = 787) and for neighbors of failing adhesions (n = 787) with circular vicinities of radius 4 µm and larger. Hence, the maturation of NAs is spatially coupled over ~4 µm. ***P < 0.001, **P < 0.01; as determined by ANOVA. The ratio for neighbors of failing adhesions is similar for all distances, as determined by Tukey’s HSD test.

Our work demonstrates a computational improvement in TM, which permits identification of tractions in nascent adhesions. The key ingredient of this method is L1- regularization, which relies assumes solution sparsity, i.e. instead of minimizing the magnitudes or derivatives of the solution the L1-norm tends to minimize the amount of non-zero coefficients in the solution[16]. L1-regularization goes back to the earliest numerical approaches to solving inverse problems[20]. However, it was only recently that L1-regularization was theoretically justified as equivalent to L0-regularization[16]. This lead to a burst of applications in MRI, electron tomography, network traffic and information technology[21]. Our results demonstrate that L1-regularization also has substantial advantages in TM (Figs. 1 and 2), including increased fidelity in the estimate of traction magnitude and impact location and the provision of an objective reference for selecting the regularization parameter (Fig. 1f). Hence, with the application of L1-based TM it should now be possible to obtain reproducible and mutually comparable traction fields.

Our result is consistent with reports based on traction reconstruction using point forces (TRPF)[22]. As discussed by Sabass et al.[9], while TRPF seems an attractive approach to overcome some of the resolution limitations of regular TM, it relies on the strong assumption that all adhesions generate forces. Our results show that only 55 % of NAs spatially coincide with local force maxima (Fig. 2p). Thus, constraining the locations of force transmission can lead to notable inaccuracy.

Although the application examples of this study focus on forces in NAs in order to showcase the gain in sensitivity and spatial resolution, L1 TM offers a more accurate resolution of traction distributions in general. This will now allow for the systematic study of molecular mechanisms that regulate cell-matrix adhesions, from the earliest lay-down of NAs to their maturation into and disassembly of larger adhesion structures.

Online Methods

Software

TM software featuring SCII bead tracking modules and traction reconstruction based on L1 or L2 norm regularization is available as Supplementary Software and at http://lccb.hms.harvard.edu/software.html; adhesion tracking software is also available at the same website.

Subpixel correlation by image interpolation (SCII) approach for quantifying substrate deformation

The displacement field was measured by tracking each bead via correlation-based image tracking between two images: an image reflecting the bead positions in the deformed substrate (f) and one reflecting the bead positions in the relaxed substrate (f) after release of cell adhesions by trypsinization. As used in previous TM studies[9, 11, 12], correlation-based tracking relies on the computation of a correlation score S between a static square template with side length 2L+1 and placed at position (x) of a particular bead in f and a moving template of the same size that is initially placed in f at the same position (x) and then shifted pixel-by-pixel to positions x + u and y + u, respectively: Here I(f,i,j)| denotes the image intensity at pixel position (x+j) within the template derived from image f, and Ī (f) is the mean intensity of the template. As shown previously[9], subtraction of the mean intensity and normalization of the correlation score relative to the intensity standard deviation eliminates effects from brightness variation between the images of the relaxed and deformed substrate.

The template size critically influences the position of the peak score. In an implementation of correlation-based tracking of fluorescent speckle flows, where the image contrast can be extremely heterogeneous, we optimized the template size for each speckle on the fly[23]. However, applied to TM this would mean that displacement vectors across the field of view originate from a different level of granularity, which can introduce significant error in the traction reconstruction. Thus, we maintained a fixed template size for all beads. Where the score function computed at the default template size showed multiple strong maxima, we increased the template size until the ambiguity disappeared. We then returned to the default template size and searched among candidate peaks in the score function for the one closest to the unique peak found at the larger template size. Subpixel localization of the maximal score value was accomplished by parabolic fitting of the 9 score values about the maximum of the pixel-based score function. Because the final displacement was defined with sub-pixel precision, we refer to this approach as pixel correlation with subpixel fitting (PCSF). This is in contrast to the subpixel correlation by image interpolation (SCII) approach, adapted from Gui et al. for particle image velocimetry[19], which relies on linearly pre-interpolated templates: where R is the refinement factor (e.g. 10) and I is an image interpolated bilinearly, where i and j are the largest integers not greater than i/R and j/R, and x and y are the fractional part of i/R and j/R, respectively, such that i/R=i+x, j/R=j+y. To minimize the increased computational cost associated with the image interpolation, we limited the interrogation range of u and u to and , respectively, where u, and u, denote the position of the correlation maximum identified by PCSF. In SCII, the template window shifts by 1/R. From the score, the maximum peak was found and better precision was achieved by interpolation through parabolic fitting with nine points around the new peak. Other types of image interpolation methods, including cubic and spline interpolations, were also tested.

Using synthetic bead images (Online Methods in Force simulation for the generation of synthetic displacement fields), we compared the performance of PCSF tracking with SCII (Supplementary Fig. 12, Supplementary Note 1). We found that SCII tracking was more accurate with enhanced detection of tractions in small adhesions and allowed for a smaller support area (template window) than PCSF. SCII is different than the correlation-based continuous window shift (CCWS) approach described in Gui et al.[19] in that the template window is interpolated for both the images of relaxed and deformed substrates, whereas Gui et al. interpolated only the template derived from the target image (equivalent to the image of the deformed substrate) and matched it with a non-interpolated template derived from an iteratively deformed source image (equivalent to the image of the relaxed substrate). The accuracy of SCII was found to be nearly the same as that of CCWS even though SCII is non-iterative (Supplementary Fig. 13).

Prior to the displacement measurement, bead images from deformed substrates were registered for translation against the bead image of trypsinized substrate : We first determined the coarser global translation between the two substrates by cross-correlation of a region that is at least 40 µm away from the nearest cell boundary followed by SCII tracking of individual beads in this region. Since cells typically occupied two corners of the image, we did not perform rotational registration, which requires all four corners of the field of view free of tractions and thus outside the cellular footprint. Nonetheless, the software release includes the option of rotational registration for data sets that fulfill this requirement.

Solution of inverse problem

2D TM assumes that tractions from cell adhesions on a gel are exerted in parallel to the gel surface, with negligible traction in the normal direction[24]. Under these conditions the displacement field u(x) is the product of a convolution between the traction field f(ξ) and a Green’s function g(x,ξ): where Ω is the domain within the cell boundary where traction is generated. Eq. 4 describes the substrate displacement in the i-th direction at the location x=(x1,x2) as a result of the superposition of all the tractions generated in Ω, where the contribution of a traction in location ξ=(ξ1,ξ2) is defined by the Green’s function. In 2D TM, the Green’s function is given by the Boussinesq solution[25]: where E is the Young’s modulus of the gel, ν is the Poisson ratio, δij is the Kronecker delta function, and |x−ξ| is the Euclidian distance between x and ξ. To infer f(ξ) from a measured displacement field u(x), we used the boundary element method (BEM), which approximates a traction field on discrete nodes[12]. For each node we defined a pyramidal-shaped basis function h(x) that is 1 at the node location itself and drops linearly to zero in all adjacent nodes. Hence the continuous traction field is approximated by discrete traction coefficients (β) as allowing us to rewrite Eq. 4 as where n is the number of nodes in Ω. In this notation we can define a discrete forward map in matrix form with matrix coefficients In this notation, u is a 2m×1 vector containing the displacement u1 and u2 of m beads and β is a 2n×1 vector containing the traction coefficients β1 and β2 for x- and y- directions. For a general mesh, calculating M is computationally expensive. We simplified this process by creating an evenly spaced square mesh so that equation (9) can be transformed and calculated in the frequency domain. This is reminiscent of the approach taken in Fourier transform traction cytometry (FTTC); however, we still solve the inverse problem in the spatial domain. Thus, we refer to this type of BEM as FastBEM.

The goal of traction reconstruction is to solve Eq. 8 for β, given a bead displacement field u(x). Due to measurement noise, the inversion of Eq. 8 is ill-posed (see main text). Therefore, it is necessary to introduce a regularization scheme and solve a minimization problem of the kind: The first term in Eq.10 defines the residual norm between the measured and predicted displacement fields. The second term defines a semi-norm, i.e. a norm of low-level spatial derivatives of the traction coefficients that are encoded by finite differences in the matrix R. The regularization parameter λ determines how much the semi-norm penalizes the minimization of the residual norm. For a semi-norm with 0th order derivatives and a regular mesh, R is the identity matrix. R can also accommodate 1st or 2nd order derivatives of the traction coefficients, which would implement the assumption that the spatial variation of the traction field is small at the length scale of the mesh size. However, tractions are transmitted at discrete spots, with inherently large spatial variation. Accordingly, it has been shown before[9] that the use of a semi-norm with 0th order derivatives better preserves punctate patterns of traction impact. Thus, our implementation of TM relies on 0th order derivatives.

In general, TM reconstructions minimize the semi-norm ‖Rβ‖2 (see Eq. 10), which is referred to as L2 regularization. In combination with 0th order derivatives (R=I) in the semi-norm, this type of regularization is known as Tikhonov regularization[10]. Tikhonov regularization suppresses the magnitude of the tractions, leading to an inherent underestimation of the traction level (see Fig. 1).

The regularization scheme involving minimization of ‖Rβ‖1 is called L1 regularization, or sparsity regularization, which can be formulated with an objective function, k, Ideal sparsity regularization would minimize ‖Rβ‖0, which is the number of nonzero entries in the traction coefficients vector β, while minimizing the residual norm. However, this kind of optimization problem is extremely difficult to solve. L1 regularization (Eq. 11) has been proposed as a numerically feasible alternative to L0-regularization[26] that provides solutions close to the ideal L0-regularization. To solve Eq. 11, we used an iterative reweighted least squares (IRLS) method[27] that solves a weighted least squares problems iteratively until solutions converge to an L1-norm minimizing solution. A brief summary of IRLS is described in Supplementary Note 2. Calculation of the forward matrix M is computationally expensive. For example, for one frame of 900 × 900 pixels (about 104 nodes in the force mesh) with 8000 displacement vectors, the forward matrix would need more than 1.2 GB of RAM, which requires access to a computer cluster.

One might think that L1-regularized problems with Boussinesq solutions could be solved in the frequency domain like in Fourier Transform Traction Cytometry (FTTC). However, doing so would not force the traction solution to be spatially sparse at each iteration. Instead, specific frequencies of tractions would be suppressed. That is, specific frequencies whose amplitudes are weakly mapped to the data in the forward problem would be impossible to recover in the inverse problem. Another fundamental difference between FTTC and BEM is that in FTTC displacement vectors are averaged and defined at grid locations to accommodate a 2D Fourier transform, whereas BEM uses raw displacement vectors defined at the locations of individual beads. Interpolation of displacement vectors onto a grid introduces uncertainties, mostly underestimation, into the reconstruction process.

Force simulation for the generation of synthetic displacement fields

To validate the traction reconstruction with synthetic data, we generated a circular area of radius r with traction impact only in the y-direction: Here, (x0,y0) denotes the center position of the circle, f denotes the maximal stress at the center, and α denotes the standard deviation of the Gaussian traction distribution controlling the traction gradient. For all simulations α was set to 1.1r. The purpose of introducing the Heaviside function H is to set the traction to zero outside the circle. To test the effect of the template side length on traction detectability, we chose an input stress f of 2 kPa and a traction impact region diameter of 6 pixels. To test the effect of the traction impact region diameter, an input stress f of 400 Pa was used. Eventually we used a range of f (200 Pa – 4 kPa) and d (2 – 12 pixel) for traction reconstructions under both L2-regularization and L1-regularization (Fig. 1l).

Synthetic bead images were created by randomly placing 2D Gaussian image signals with normalized amplitudes varying from 0.3 to 1. The standard deviation of the Gaussians was set to match a Gaussian approximation of the microscope’s point spread function[28]. The bead image of the deformed substrate was created with new bead positions calculated by the forward Boussinesq solution (Eq. 5) to Eq. 4 for a given traction distribution. With these two images, the displacement field was measured using correlation-based tracking, and traction fields were calculated based on L2 or L1 regularization. To validate the displacement measurements or traction reconstructions, we defined the traction detectability Ψ as where the numerator is the average traction in the circular traction impact region (A0), and the denominator is the maximal traction magnitude among all tractions in the background. The RMS error of the traction reconstruction was estimated as ζ=‖f−f‖2.

For validations with multiple traction impact regions, reflecting the presence of multiple adhesions, we first superimposed the traction fields of all hypothetical adhesion sites and then continued with the same procedure for displacement field simulation and traction reconstruction as discussed for the single traction impact regions. To better reflect the co-existence of NAs and FAs, we introduced traction impact regions of different sizes and different levels of shape anisotropy. For all validation tests, we used Young’s modulus E of 8 kPa and Poisson’s ratio ν of 0.5 for a virtual gel substrate.

Nascent adhesion (NA) detection, focal adhesion (FA) segmentation, and adhesion tracking

NAs in paxillin images were detected using an existing algorithm from single particle tracking[29, 30]. Briefly, TMR-paxillin images were filtered using the Laplacian of Gaussian filter and then local maxima were detected. Each local maximum was then fitted with an isotropic Gaussian function (standard deviation: 2.1 pixel) and outliers were removed using a goodness of fit test (p=0.05). The fluorescence intensity of the adhesion was considered to be the amplitude of the Gaussian, and the traction magnitude was found as the maximum value in the circular region around the NA of radius of one mesh length (0.5 µm). Likewise, NAs and local traction maxima were considered colocalized (Fig. 2p) if they were less than 0.5 µm apart, which is the length scale of the mesh used for traction reconstruction.

FA segmentation was performed based on a combination of Otsu and Rosin thresholding after image preprocessing with noise removal and background subtraction. Segmented areas were categorized as focal contacts (FCs) or FAs based on their areas following the criteria prescribed by Gardel et al.[4] (0.24 µm2 for FCs, 0.6 µm2 for FAs). NAs ware tracked with the uTrack software[30] using the gap closing option (maximum gap = 5 frames) and a Brownian search band of 2–5 pixels. Tracks that overlapped with an FC or FA were labeled as a maturing adhesion. Tracks that initiated and terminated as NAs were labeled as a failing adhesion..

Cell culture

PtK1 cells (ATCC, mycoplasma tested) were grown in Ham’s F-12 Nutrient Mix supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin. The cells were cultured in a humid atmosphere containing 5% CO2 at 37°C.

Gel substrate preparation and cell transfection

Silicone gel substrates were prepared in glass-bottomed 35 mm dishes (GWSB-3522 with #1.5 coverglass bottoms from WillCo Wells B.V.). Tracer particles, 40 nm carboxylated far-red fluorescent beads (excitation/emission 690/720 nm, from Invitrogen) were deposited on the coverglass bottoms at a low concentration to serve as fiduciaries of the glass surface. Pre-polymers for 2 and 8 kPa gels were prepared by mixing the components B and A of QGel 920 (Quantum Silicones) high refractive index (n = 1.49) silicone gel at ratios of 1.07:1 and 1.15:1, respectively, as described in Gutierrez et al.[7] The pre-polymers were spin-coated onto the coverglass bottoms at 2500 rpm for 30 sec and then baked for 2 hrs in a 100 °C oven, producing ~34 µm thick layers of cured gels. The gels were then treated with (3-Aminopropyl)triethoxysilane (APTES) to functionalize their surfaces. The same 40 nm carboxylated far-red fluorescent beads were covalently linked to the gel surfaces by incubating the gels under a suspension of the beads (1:10000 dilution from the 5% stock suspension) in 20 mM HEPES, pH = 8, with 0.01% EDC as a catalyst. To assess the each gel’s thickness, a fluorescence microscope with an oil-immersion objective was focused on beads on the surface of the gel and then on beads on the surface of the glass; the difference between the readings of the nosepiece knob of the microscope (z-axis positions of the objective) was calculated (and multiplied by 0.98 to correct for the difference between the refractive indices of the immersion oil, 1.515, and the gel, 1.49). The elastic moduli of the gels were measured with the microfluidic technique described in Gutierrez et al.[31] by assessing the deformation of a ~70 µm thick layer of the gel on a 35×50 mm #1.5 microscope cover glass under a controlled shear flow. To minimize the experimental uncertainty, the ~70 µm thick samples were prepared for each batch of the gel pre-polymer and cured together with the ~34 µm thick gels used in experiments on cells (2 hrs in a 100 °C oven). The Green’s function for a gel with this thickness was found to be minimally different from the Boussinesq solution that we used for this study (Supplementary Fig. 14).

To facilitate cell adhesion, the substrates were coated with fibronectin (FN) by incubation with 50 µg/ml of FN and 100 µg/ml 1-Ethyl-3-(3-dimethylaminopropyl) carbodiimide in PBS, (pH 7.4) for 15 min at room temperature. The coated dishes were washed 3 times with PBS and filled with cell media before cell seeding. Cells were transfected with HaloTag-paxillin using the Neon Transfection System (Invitrogen) right before being seeded. After 48 h in the incubator, cells were labeled by incubation with 5 µM HaloTag TMR ligand (Promega) for 15 min at 37°C with 5% CO2. Cells were then washed three times with warm growth medium to remove unbound ligands. After the final wash, cells were incubated in fresh growth medium for 30 min in the incubator. Right before imaging, in order to protect the cells from photo-damage during live imaging while allowing for subsequent trypsin treatment, the imaging chamber was fitted with inlet and outlet tubing, sealed with Valap[32], and filled with imaging media – Leibovitz’s L-15 medium (Life technology) without phenol red supplemented with 1% penicillin-streptomycin, 10% FBS, 0.45% glucose and 1% oxyrase – using a 10-mL syringe.

Live cell imaging

Cells were imaged in an enclosed 37°C, 5% CO2 incubation chamber. Images were taken on the 561 nm (paxillin TMR) and the 642 nm (beads) channels with a 500 ms exposure time every 5 seconds for 200 frames using a Nikon Ti Total Internal Reflection Fluorescence (TIRF) microscope using a 100× objective with 1.5× additional magnification factor with a Hamamatsu ORCA-D2 CCD camera (Hamamatsu Corporation, Bridgewater. NJ, USA, final resolution: 43 nm/pixel). To obtain a bead image from the undeformed substrate, cells were removed from the substrate by injecting a high dose (0.5%) of 5 mL Trypsin/EDTA (Invitrogen) for 30 min.