 Research Article
 Open Access
 Published:
Full L_{1}regularized Traction Force Microscopy over whole cells
BMC Bioinformaticsvolume 18, Article number: 365 (2017)
Abstract
Background
Traction Force Microscopy (TFM) is a widespread technique to estimate the tractions that cells exert on the surrounding substrate. To recover the tractions, it is necessary to solve an inverse problem, which is illposed and needs regularization to make the solution stable. The typical regularization scheme is given by the minimization of a cost functional, which is divided in two terms: the error present in the data or data fidelity term; and the regularization or penalty term. The classical approach is to use zeroorder Tikhonov or L_{2}regularization, which uses the L_{2}norm for both terms in the cost function. Recently, some studies have demonstrated an improved performance using L_{1}regularization (L_{1}norm in the penalty term) related to an increase in the spatial resolution and sensitivity of the recovered traction field. In this manuscript, we present a comparison between the previous two regularization schemes (relying in the L_{2}norm for the data fidelity term) and the full L_{1}regularization (using the L_{1}norm for both terms in the cost function) for synthetic and real data.
Results
Our results reveal that L_{1}regularizations give an improved spatial resolution (more important for full L_{1}regularization) and a reduction in the background noise with respect to the classical zeroorder Tikhonov regularization. In addition, we present an approximation, which makes feasible the recovery of cellular tractions over whole cells on typical fullsize microscope images when working in the spatial domain.
Conclusions
The proposed full L_{1}regularization improves the sensitivity to recover small stress footprints. Moreover, the proposed method has been validated to work on fullfield microscopy images of real cells, what certainly demonstrates it is a promising tool for biological applications.
Background
Tissue remodeling implies the reorganization of the extracellular matrix (ECM), which is driven by the conversion of intracellulargenerated mechanical forces into extracellular traction, which reorganizes the ECM fibers. This is a crucial process during regeneration (e.g., in wound healing) but it is equally important in pathologic scenarios (e.g., in inflammation and/or cancer). In fact, some of the pathological changes associated to these diseases are due to the altered ability of cancer, or inflammatory cells, to exert abnormal traction on their microenvironment [1]. Thus, the precise quantification of the traction exerted by cells (and groups of cells) is key to understand the remodeling processes that occur during these biologically relevant events, including events on the cellular scale (e.g., cell migration and division) as well as at a subcellular scale (e.g., signal transduction events). Cells convert intracellulargenerated forces into extracellularapplied tractions and transmit them to the microenvironment through micronsized protein accumulations termed focal adhesions (FA) [2,3,4,5,6,7,8].
Traction Force Microscopy (TFM) is a technique that estimates the tractions exerted by adherent cells placed on top of a flexible hydrogel that mimics the mechanical properties of the ECM [9].
In a typical TFM experiment, a large number of fluorescent beads are randomly mixed inside the hydrogel. The fluorescent beads are displaced as the cells exert tractions on the substratum and act as fiduciary markers enabling the transformation of motion into force (see Fig. 1). Experimentally, fluorescence microscopy is used to acquire images of the beads before (stressed hydrogel) and after (relaxed hydrogel) suppressing cell tractions, either by druginduced relaxation of cells, or directly via cell detachment or lysis. The images of the stressed and relaxed hydrogel are then compared to track the motion of the embedded beads and obtain the matrix displacement field. This, together with a mechanical model of the hydrogel based on its height and its elasticity (Young’s modulus) allows estimating the cellular tractions using computational methods.
The displacement field of the beads can be obtained using different image processing methods. Two classical approaches used in TFM experiments include blockwise image correlation computation (Particle Image Velocimetry) [10]; and individual motion tracking of each bead (Particle Tracking Velocimetry) [11]. One of the novel aspects of the present study is that we use a nonrigid image registration approach, which provides very good accuracy with a reasonable computation time and without a loss in spatial resolution compared to previous methods [12].
Recovering cellular tractions (causes) from matrix displacements (observable effects) implies solving an illposed inverse problem. Consequently, the solution is very sensitive to the errors committed in the calculation of the displacement field.
To constraint the solution and prevent the estimated tractions from overfitting the noise present in the matrix displacements, it is common to regularize the traction reconstruction process, including a penalty term in the norm of the tractions [13]. Historically, L_{2}norm regularization, also known as Tikhonov regularization [14], is the most widely used scheme in TFM, as it allows an efficient traction recovery in the frequency domain [15]. However, recent studies showed that the use of the L_{1}norm in the imposed penalty term yields improved results in terms of spatial resolution, estimation of the traction magnitude and background noise reduction, giving sparser traction fields as is expected in TFM experiments [16,17,18]. For example, L_{1}regularization was introduced to study the focal adhesions maturation, which was possible by the huge improvement obtained in the spatial resolution [16]. Using a similar approach Bohr et al. in [17], demonstrated that traction fields with high spatial resolution could be computed from displacement fields with low resolution.
A major issue with L_{1}regularization is that it requires solving the problem in the spatial domain instead of the frequency domain, which greatly increases the computational demand. Therefore, its application has been limited to either reduced regions of interest with fullresolution displacements, or full field of views with lowresolution displacements.
To overcome this limitation, in this study, we introduce an approximation that reduces the computational cost of recovering the tractions in the spatial domain, thus, providing the means to work with both fullfield microscopy images (covering whole cells) and fullresolution displacements. In addition, we propose the use of full L_{1}regularization in TFM, which minimizes the L_{1}norm of both the penalty (i.e., recovered tractions) and the data fidelity terms (i.e., difference between the estimated displacements and the ones obtained from the given recovered tractions). Here, the performance of the proposed scheme satisfactorily compares to the classical approach (regularization minimizing the L_{2}norm for both the penalty and the data fidelity terms) and the simple L_{1}regularization (constraining the L_{1}norm of the penalty term) with synthetic and real data.
In the following sections, we present the foundation of 2D TFM traction recovery in the spatial domain and the method implemented to compute the traction forces for fullsize microscopy images, followed by the evaluation of synthetic data and experiments with live cells, respectively. Results section contains the experimental results, including a qualitative and quantitative comparison between the different regularization schemes, followed by a discussion on the limitations of the present approach and future perspectives on the application of this technique.
Methods
Regularized Traction Force Microscopy
In this work, we consider a bidimensional (2D) TFM scenario in which cells are cultured on the top of a controlled thickness hydrogel. One assumption is that cells only exert inplane (shear) forces on the surface of the hydrogel. In this case, hydrogel deformation is obtained as the solution of the elasticity problem:
where the jth component (j ∈ {x, y}) of the displacement field u at spatial location x is related to the tractions t exerted along the ¡th Cartesian direction (¡ ∈ {x, y}) at the location x′ by the Green’s function of the hydrogel g(x,x′).
The Green’s function models the elastic deformation of the hydrogel in response to a point force. It is usually given by the analytical Boussinesq solution in 2D TFM, where the hydrogel can be approximated by a homogeneous, isotropic, elastic and semiinfinite medium. Under these assumptions, the Green function can be considered linear shiftinvariant [19], which allows transforming Eq. 1 in a summation of convolutions:
In practice, displacement and traction fields are only estimated at discrete locations of the space, and Eq. (2) can be rewritten in matrix form as:
where u = [u _{ x }(x _{1}); … ; u _{ x }(x _{ N }); u _{ y }(x _{1}); … ; u _{ y }(x _{ N })] is a 2N × 1 column vector with N being the number of locations where displacements are measured, \( \mathbf{t}=\left[{t}_x\left(\boldsymbol{x}{'}_1\right);\dots; \left({\boldsymbol{x}}_M^{'}\right);{t}_y\left({\boldsymbol{x}}_1^{'}\right);\dots; {t}_y\left(\boldsymbol{x}{'}_M\right)\right] \) is a 2M × 1 column vector with M being the number of locations where tractions will be estimated, and K is the 2N × 2M stiffness matrix given by:
with
A trivial solution for the traction recovery is obtained by the direct inversion of Eq. 3. However, it implies the inversion of the stiffness matrix, which is illconditioned, and thus, it would amplify the errors present in the computed displacement field leading to unstable results. For this reason, a regularization scheme is commonly used to reach a stable solution. Then, the traction field recovery in the spatial domain is given by the minimization of the following cost functional:
where ‖·‖_{ q } denotes the L_{q}norm, ∥Kt − u∥_{ q } is the data fidelity term, ∥t∥_{ p } is the imposed penalty and λ is a parameter that controls the amount of regularization applied.
L_{2}Regularization in the Fourier Domain
Following Parseval’s theorem, L_{2}regularization has an efficient, and widely used, implementation in Fourier domain [15], where the summation of convolutions in the spatial domain (Eq. 2) is transformed in a summation of products:
where \( \overset{\sim }{u},\overset{\sim }{g},\overset{\sim }{t} \) are the Fourier transforms of u , g , t in Eq. 2 and k are the frequencies where these variables are evaluated.
Analogously to the spatial domain, the previous equation can be written in a matrix form:
where \( \overset{\sim }{\boldsymbol{A}} \) contains all the values of \( {\overset{\sim }{g}}_{ji} \) in Eq. 7.
Then, the L_{2}regularization can be formulated as the minimization of a cost functional, taking the L_{2}norm for both penalty and data fidelity terms:
being ℱ^{−1} the inverse Fourier transform.
The solution for the traction field is, thus, given by:
where \( {\overset{\sim }{\boldsymbol{A}}}^{\ast } \) is the complex conjugate of \( \overset{\sim }{\boldsymbol{A}} \) .
L_{2}Regularization in the Spatial Domain
Zero order Tikhonov regularization can be also formulated in spatial domain, using the L_{2}norm for both terms in Eq. 6 (p = q = 2). Then, the minimization problem has an analytical solution given by [20]:
with K ^{T} being the transpose of the stiffness matrix.
L_{1}Regularization in the Spatial Domain
A better alternative is to use the L_{1}norm for the penalty term (q = 2 and p = 1 in Eq. 6), which promotes sparser traction fields [21], as the ones expected for TFM experiments. Unfortunately, it is not possible to express the solution in a closed form, thus requiring an optimization algorithm. Here, the Iterative Reweighted Least Squares (IRLS) algorithm [22] is used to solve the convex optimization problem and estimate \( {\widehat{\mathbf{t}}}_{{\mathrm{L}}_1}, \) as was used in [16, 18]. This algorithm approximates the L_{1}norm by a weighted L_{2}norm problem, updating those weights iteratively until convergence is reached. Then, at each iteration s:
where W ^{s} is the weight matrix at the s th iteration. We considered that convergence was reached when the absolute value of the difference between the tractions estimated in two consecutive iterations was smaller than 10^{−6}.
Full L_{1}Regularization in the Spatial Domain
The alternative, we propose in this manuscript, is to use the L_{1}norm for both data fidelity and penalty terms (p = q = 1 in Eq. 6). Analogously to the simple L_{1}regularization, there is no a closed form solution for the convex optimization problem and the IRLS algorithm can be used to estimate the cellular tractions. Then, the solution at each iteration s is given by:
where \( {\mathbf{W}}_d^s \) is the weight matrix of the data fidelity term at the sth iteration.
Stiffness matrix reduction
The main issue when working in the spatial domain with fullsize real images and fullresolution displacements is the high computational requirements. For instance, for a typical field of view of 800 ∙ 800 pixels and storing each element of the stiffness matrix with double precision, its size in memory would be as large as is 2N ∙ 2M ∙ 8Bytes ~ 13 Terabytes, where N = M = 800 ∙ 800, being N the positions of the displacements used for the tractions estimation and M, the positions where the tractions would be recovered.
Consequently, solving the regularization problem in the spatial domain is not feasible on a standard desktop computer. To avoid this limitation and be able to work with both fullsize real images and fullresolution displacements, we perform a reduction of the stiffness matrix by including a priori information inherent to TFM experiments (see Fig. 2). In particular, it is known that adherent cells exert tractions at discrete locations; namely, the focal adhesions. Thus, it is reasonable to assume that cellular tractions should be only estimated in certain regions, outside of which they should be zero. Note that given the conditions of the experiments we performed, we are not able to detect individual focal adhesions, but clusters of them as determined by cytoskeletal fibers. We refer to this approximation as stress footprints throughout the manuscript.
First, we consider that the regions presenting large displacements correspond to the sites of traction application, which is a priori a sufficient approximation, as the support of the traction field is much reduced than the one of displacement field. Hence, we can limit the number of locations where the tractions will be recovered (M) to the points inside these regions.
Furthermore, the signaltonoise ratio of the displacements decreases with its magnitude and their information content quickly vanishes with the distance to the locations where tractions are exerted. Therefore, it is also reasonable to limit the locations of the measured displacements that will be used in the traction recovery (N) to those ones inside the segmented regions, as those outside would mainly contribute with noise in the system of equations.
To segment these regions, we apply an Otsu’s thresholding to the magnitude of the estimated displacement field. On our hands, the segmentation of the displacements extract disconnected regions of reduced size. Thus, the tractions in each of the disjoint regions can be recovered using only the information from the displacements located inside it. Consequently, the global stiffness matrix can be split in as many local stiffness matrices (with a drastically reduced size) as detected regions of interest, making the traction estimation feasible on a standard PC. Once all the local TFM problems have been solved, the contributions of each segmented region are combined together to maintain the original field size.
All the methods described in this section have been implemented in Matlab (MathWorks Inc.) code.
Synthetic data
Simulated data generation
To assess the performance of the different regularization schemes, we used our TFM simulator [23] to generate images of relaxed and stressed hydrogel containing embedded fluorescent beads as acquired by an optical microscope (see Fig. 3). Briefly, a simulated traction field was used to generate the displacements (Eq. 3) that move a large number of randomly distributed beads in the relaxed hydrogel to their new locations in the stressed hydrogel. Furthermore, the images containing the beads were generated taking into account the noise and the point spread function associated blurring that is introduced by the optical system and thus, represent a realistic synthetic dataset.
Then, the Bspline based nonrigid registration algorithm [12] was used to calculate a fullresolution displacement field from these images, as done in real TFM experiments. Lastly, the estimated displacements were fed to the regularized traction reconstruction algorithms under evaluation and the recovered tractions were compared to the groundtruth given by the simulated fields.
To create a realistic map of tractions that can be used with our TFM simulator, we used the tractions recovered from TFM experiments with real cells. Specifically, real traction fields computed with full L_{1}regularization were segmented using an extension of the classical isodata algorithm [24] that automatically selects five thresholds using a clustering approach. Then, we perform a binary segmentation based on the most restrictive of them to be used as the synthetic stress footprints mask.
Then, following the pipeline explained before, ten different realizations have been created for each simulated case, changing in each realization the position of the fluorescent beads. The selection of the optimal parameter λ for each regularization method has been performed by an exhaustive search. Namely, the inverse problem has been sequentially solved on a synthetic image randomly selected using a regularization parameter λ on the range [0,1] with a step size Δλ of 0.01. The regularization parameter λ used on the simulations for the evaluation of the algorithms was the one minimizing the cost functional defined on Eq. 6.
Error metrics
We have evaluated the error in the recovered traction field within a region of interest defined by the corresponding stress footprint. These footprints were obtained by segmenting the magnitude of the recovered tractions using an Isodata thresholding of two levels. These metrics are evaluated only in the recovered stress footprint areas.
Being t _{ gt } the simulated groundtruth traction field, \( \widehat{\mathbf{t}} \) the retrieved traction field, and P and Q the total number of points within the respective stress footprints, we have defined the following error metrics:
Error in magnitude
Absolute error in the recovered traction magnitude (in percentage) defined as
where \( \widehat{n} \) is the number of recovered stress footprints, and P _{ n } the points of the simulated stress footprint associated to each recovered footprint.
Error in angle
Average angular error (in percentage) in the recovered stress footprint, weighted at each point by the traction magnitude. It is normalized by the maximum error in degrees (180°) to give the result in percentage as given by
Error in area
Absolute error of the recovered stress footprint area (in percentage) as given by
with \( \widehat{A} \) and A _{ gt } being the area of the recovered and groundtruth footprints, respectively.
Loss ratio of stress footprint areas
It is the ratio (in percentage) between stress footprints that the regularized TFM methods cannot recover (the number of false negatives) and the total number of simulated stress footprints synthetized by cell. It is defined as
where \( \widehat{n} \) and n _{ gt } are the number of recovered and groundtruth stress footprints.
Smallest detectable stress footprint areas
It is the smallest product between the area (in μm ^{2}) and the peak magnitude (in kPa) of the recovered stress footprint areas. This metric reflects the dependence between the size and the magnitude of the stress footprints and it is defined as
where \( {A}_{gt}^i \) and \( {t}_{max}^i \) are the ideal area and maximum traction amplitude of a recovered stress footprint i.
Statistical analysis
The statistical significance among different realizations and regularization schemes was evaluated for each error metric using the tStudent test with signrank Matlab function. A pvalue smaller than 0.01 was considered to give statistical significance.
Real data
We have used the presented algorithms to quantify the tractions exerted by Chinese Hamster Ovary (CHOK1, American Type Culture Collection (ATCC; Rockville, MD), CCL 61) cells expressing LifeactGFP. These cells were cultured on a ~ 90 μm thick polyacrylamide hydrogel with a Young’s modulus of 5 kPa and a Poisson ratio of 0.45, containing embedded fluorescent polystyrene beads (0.2 μm in diameter, 647 nm of emission wavelength, ~1 bead per μm^{2}). Images of the cells were acquired at multiple locations of the hydrogel surface with a 60× glycerol objective (NA = 1.3) mounted in a Leica SP5 Laser Scanning Confocal microscope with 0.32 μm of resolution in each direction of the plane. Timelapse images of the stressed hydrogel were taken every 20 s during 25 min. An image of the relaxed hydrogel was taken after removing the cells.
Results
Synthetic data
To build a realistic synthetic dataset we used our TFM simulator with data from ten different real TFM experiments (i.e., cells) as explained in the previous section. A wide variety of test conditions were considered, specifically, the number of stress footprints per simulated case ranged between 10 and 120, their peak magnitude between 0.6 kPa and 4 kPa and their area between 0.25 μm^{2} and 20 μm^{2}. The parameters to simulate the mechanics of the hydrogel and the acquisition of the bead images were chosen to be the same than those of the real experiments and are given below. Furthermore, multiple (10) realizations were simulated for each synthetic case to take into account the variability introduced by the random distribution of the beads inside the hydrogels. Finally, for the sake of comparison, the parameters controlling the amount of regularization were fixed for all the setups to \( {\lambda}_{L_2}=0.2 \), \( {\lambda}_{L_1}=0.08 \), \( {\lambda}_{\mathrm{full}{\mathrm{L}}_1}=0.1 \). These values were independently determined for each regularization scheme after an exhaustive search minimizing the least square error between recovered and simulated traction fields. Additional file 1 shows an illustration of the changes in the recovered traction fields with the regularization parameter for each method.
Figure 4 shows a representative example of simulated traction field (Fig. 4a) and ideal displacement field (Fig. 4b). The noisy displacements computed from the synthetic images of relaxed and stressed hydrogel are shown in Fig. 4c. Their corresponding tractions recovered by each regularization scheme are shown in Fig. 4df. The mask defining the regions were tractions are estimated is overlaid. As expected, a visual inspection of Fig. 4 reveals that L_{1}regularizations give sparser traction fields with less background signal and better estimation of traction magnitudes than the classical L_{2}norm based Tikhonov regularization. Furthermore, full L_{1}regularization is able to recover more stress footprints than the other methods.
Those results are quantitatively confirmed in Fig. 5 and Additional file 2. We observe how L_{1}regularization presents smaller errors in magnitude and angle of the recovered stress footprint areas than Tikhonov regularization and full L_{1}regularization, giving around three times smaller magnitude error (12%) than the other methods. Regarding the area error of the recovered stress footprints, full L_{1}regularization shows the best performance with an average error less than 4%, while L_{1}regularization has about 7% of error and Tikhonov present the worst performance with an average area error about 16%. Looking at the ratio of the lost stress footprints, full L_{1}regularization is able to recover more than the double of stress footprint areas (around 75%) than the other regularization methods (around 30% of them). In addition, in Fig. 6, it can be seen how full L_{1}regularization has the smallest detectable stress footprint area values. Namely, this method detects smaller stress footprint areas with reduced traction magnitude than the other methods. In particular, the smallest stress footprint areas detected in all experiments for each method are (Peak Magnitude (kPa), Area (μm ^{2}), mFA(nN)): L_{2}regularization (0.89, 0.38, 0.33); L_{1}regularization (0.94, 0.31, 0.29); full L_{1}regularization (0.68, 0.25, 0.17).
Moreover, most of the differences between the error measures given for the different regularization schemes are statistically significant (p < 0.001). The exceptions are: The error in magnitude between Tikhonov and full L_{1}regularization and the smallest footprint detected between Tikhonov and L_{1}regularization.
Real data
In Fig. 7, an example of traction recovery on one frame of a real cell is shown. In Additional file 3, ten samples frames are presented. It can be seen how the results are qualitatively similar to those obtained for synthetic data.
Image of the cell (Fig. 7a), images of the stressed and relaxed polyacrylamide hydrogels (Fig. 7b) and the displacement estimation (Fig. 7c) are shown with the mask used for the tractions calculation, overlaid. In Fig. 7df, the recovered traction fields for each of the regularization schemes are shown. The Tikhonov regularization clearly smoothens the solution and introduces more background noise than L_{1}regularizations. Moreover, fullL_{1}regularization enables to distinguish between closer stress footprints as smaller stress footprints are recovered.
Discussion
The use of the full L_{1}norm was previously introduced [25] and has been applied in different research areas, including biomedical applications [26,27,28], image processing [29, 30] and computer vision [31]. Whereas L_{2}data fidelity is useful for fields contaminated with Gaussian noise, microscopy images often contain more complex noise (mixed Poisson and Gaussian noise) as well as outliers. In these cases, applying the L_{1}norm data fidelity constraint is more appropriate.
Moreover, there has been previous theoretical work exploring the consequences of the use of the L_{1} norm in the data fidelity for image reconstruction [30]. In particular, it has been shown that some undesirable characteristics derived from the minimization of the absolute error instead of the leastsquares error such as the lack of uniqueness of solutions, and the lack of continuous dependence on data, can be real assets. A major advantage of the L_{1} fidelity based model over the standard one is that the regularization imposed on solutions is more geometric, in the sense that the regularization process has less dependence on the contrast of image features than on their shapes; As distinct from the standard model, small features in the images maintain their contrast. In addition, L_{1} data fidelity works better reducing some artifacts like ringing and the presence of outliers [25, 29].
In Additional file 1, a comparison between the data fidelity and penalty term recovered by the regularization methods for a range of lambda values (from 0 to 0.22 with a step of 0.02) is shown. It is clear that full L_{1}regularization reduces the data fidelity term. This means that the matching between the estimated displacements and the inferred ones (the product of the stiffness matrix and calculated traction field) is greatly improved with respect to the other regularization schemes. This could be responsible for the strong decrease in the loss stress footprints ratio and it is the source of differences between the L_{1}regularizations.
Our results reveal that L_{1}regularizations give an improved spatial resolution (more important for full L_{1}regularization) in terms of sensitivity to detect smaller stress footprints in both synthetic and real data with respect to the classical zeroorder Tikhonov regularization.
In our hands, L_{1}regularization of the penalty term presents smaller errors in magnitude and angle of the recovered stress footprint areas than the other regularization methods. It was expected that L_{1}regularization of the penalty term overperformed the results obtained by Tikhonov regularization as was already demonstrated in previous publications [16, 18]. Nevertheless, it was a surprising result for full L_{1}regularization and it is maybe due to our error metric evaluation procedure. Namely, we compute the error metrics over all the recovered stress footprints being the number recovered by full L_{1}regularization much larger. Even more, the ones recovered only by full L_{1}regularization tend to be the smaller in magnitude and area. Those risk to be the ones over which the reconstruction errors could be larger penalizing the computation of the error metrics.
As detailed in the previous sections, working in the spatial domain implies a high computational demand in terms of RAM and CPU usage due to the size of the stiffness matrix and the need to use an iterative optimization algorithm to recover the tractions, respectively. To overcome this severe problem, we have presented a method to reduce the size of the stiffness matrix that allows working with fullfield microscopy images and full resolution displacement fields in the spatial domain. It is based on the calculation of forces in a few confined regions of interest defined with an Otsu thresholding over the magnitude of the estimated displacements. This approach results in a generous segmentation of the regions, guaranteeing no loss of information. Namely, more than 85% of the estimated displacements magnitudes are kept in all cases. This can be qualitatively observed in Fig. 8 for one of the synthetic cases. The stress footprints recovered using L_{2}regularization in the Fourier domain (using all data from the estimated displacement field) and in the spatial domain (using just the data covered by the mask) are visually similar, and no loss of information is apparent as shown in Fig. 8 (top) and (middle). Indeed, this approximation reduces the background signal but it does not cause loss of information of the recovered stress footprints (see Fig. 8 (bottom)). Moreover, there are not statistically significant differences between the quantitative error measures for both methods (ttest, p < 0.001) when computed on our synthetic dataset. As this solution is largely conservative, we believe that further refinements could reduce the computational demand even further. Finally, there exists an alternative TFM strategy called Traction Reconstruction with Point Forces (TRPF) [32] that also restricts the reconstruction of cellular tractions in the space to the focal adhesion sites at the cost of labeling them. Furthermore, in contrast to the method presented in this work, it would miss the tractions generated under the cellular nucleus in TFM setups considering transversal forces (known as 2.5D TFM setups) [33].
Recent publications have implemented novel methods incorporating relevant biomechanical constraints (i.e., the imposition of no forces outside the cell area and the assumption that if the cell is in equilibrium the sum of forces over the whole cell should be zero) [34, 35]. In our work, we have not imposed those constraints in the recovery of the tractions due, mainly, to two reasons: The first one is that in our real data, due to experimental limitations, the cell membrane might not be completely labeled. Consequently, it is not possible to accurately delineate the cell contour and, therefore, to impose that no forces are applied outside the cell. Regarding to the constraint of having the distribution of forces over the whole cell to sum up to zero, it has not been imposed in our experiments because in the stiffness matrix reduction approach the whole image is split in subregions to solve the traction recovery problem in small areas and, then, reduce the computational demand. This reduction implies solving the problem not in the whole cell and thus, it is not possible to add the force balance constraint. Nevertheless, for our real data, it has been verified that this sum is close to zero (less than 0.1% of the recovered maximum traction magnitude for all cases, see Additional file 4) (mean (Pa), standard deviation (Pa)): L_{2}regularization in Fourier Domain (−0.06, 0.24); L_{2}regularization in spatial domain (−1.26, 18.07); L_{1}regularization in spatial domain (−2.98, 18.59); full L_{1}regularization in spatial domain (−2.04, 17.88).
To reduce the computation time of the L_{1}regularization methods (full L_{1}regularization is about five times slower than L_{1}regularization and 20 times slower than L_{2}regularization in spatial domain), a more efficient implementation could be performed in a more suitable platform (i.e., using CUDA running over the FPGA), which would take advantage of every available resource of a desktop computer.
It is worth noting that we have simulated very small stress footprints in some of the synthetic traction maps. Those are very difficult to recover due to a variety of experimental constraints, namely, the random position of the beads, the limited microscope resolution and the possible proximity of other stress footprints with higher size or magnitude. A limitation of our experimental setup is the fact that we use relatively large beads (200 nm) restraining the density of fiducial markers embedded in the hydrogel to avoid clustering. This issue could be partly solved by using smaller beads (~40 nm) that could be packed with a higher density and closely to the surface [16]. A recently proposed and elegant alternative is to use a polydimethylsiloxane (PDMS) silicone hydrogel with printed static quantum dots of 200 nm of diameter distributed regularly (with 1,5 μm of separation) [36].
As per the experiments with real data we have only performed a qualitative comparison of the different regularization methods. Our current work focuses on designing experiments in which the cellular focal adhesions will be stained in vivo. For this purpose, we are constructing cells stably expressing bona fide markers of focal adhesion coupled to fluorescent markers (e.g. EGFP), which are commonly used in these types of experiments [6]. Our motivation is to build ground truth data in which would be possible to establish quantitatively the resolution reached by the different regularization schemes.
We believe this improvement in sensitivity and spatial resolution could be relevant when combined with advanced microscopy methods, for cell biomechanics studies willing to drive the focus toward the characterization of individual focal adhesions instead of clusters of them as historically performed.
Conclusion
In this manuscript, we have implemented full L_{1}regularization to recover cellular tractions in TFM experiments, showing its good performance compared with classical Tikhonov regularization and simple L_{1}regularization both in simulated and real data.
The main characteristic of full L_{1}regularization is the huge improvement in the sensitivity and the spatial resolution, which allows recovering smaller stress footprints compared with the stateoftheart regularization methods. In addition, we have demonstrated the improved performance of L_{1}regularizations, giving smaller errors in the recovered traction field and resulting in less background noise than Tikhonov regularization.
In this work, we have also presented a method to make feasible the recovery of the tractions exerted by whole cells on fullfield microscope images when working in the spatial domain. We have further demonstrated the suitability of the approach for the analysis of both realistic simulations and real data experiments.
Abbreviations
 CHOK1:

Chinese Hamster Ovary
 CUDA:

Compute Unified Device Architecture
 ECM:

Extracellular Matrix
 FA:

Focal Adhesion
 FPGA:

Field Programmable Gate Array
 GFP:

Green Fluorescent Protein
 IRLS:

Iterative Reweighted Least Squares
 LFA:

Lost ratio of stress Footprint Areas
 PDMS:

PolyDiMethylSiloxane
 PSF:

Point Spread Function
 SF:

Stress Footprint
 SFA:

Smallest detectable stress Footprint Areas
 TFM:

Traction Force Microscopy
 TRPF:

Traction Reconstruction with Point Forces
References
 1.
Butcher DT, Alliston T, Weaver VM. A tense situation: forcing tumour progression. Nat Rev Cancer. 2009;9(2):108–22.
 2.
Heisenberg CP, Bellaïche Y. Forces in tissue morphogenesis and patterning. Cell. 2013;153(5):948–62.
 3.
Friedl P, Gilmour D. Collective cell migration in morphogenesis, regeneration and cancer. Nat Rev Mol Cell Biol. 2009;10(7):445–57.
 4.
Brugués A, Anon E, Conte V, Veldhuis JH, Gupta M, Colombelli J, et al. Forces driving epithelial wound healing. Nat Phys. 2014;10(9):683–90.
 5.
Discher DE, Janmey P, Wang YL. Tissue cells feel and respond to the stiffness of their substrate. Science. 2005;310(5751):1139–43.
 6.
Gardel ML, Sabass B, Ji L, Danuser G, Schwarz US, Waterman CM. Traction stress in focal adhesions correlates biphasically with actin retrograde flow speed. J Cell Biol. 2008;183(6):999–1005.
 7.
VicenteManzanares M, Ma X, Adelstein RS, Horwitz AR. Nonmuscle myosin II takes centre stage in cell adhesion and migration. Nat Rev Mol Cell Biol. 2009;10(11):778–90.
 8.
Balaban NQ, Schwarz US, Riveline D, Goichberg P, Tzur G, Sabanay I, et al. Force and focal adhesion assembly: A close relationship studied using elastic micropatterned substrates. Nat Cell Biol. 2001;3(5):466–72.
 9.
Butler JP, Tolic’Nørrelykke IM, Fabry B, Fredberg JJ. Traction fields, moments, and strain energy that cells exert on their surroundings. Am J Phys Cell Phys. 2002;282(3):C595–605.
 10.
TolićNørrelykke IM, Butler JP, Chen J, Wang N. Spatial and temporal traction response in human airway smooth muscle cells. Am J Phys Cell Phys. 2002;283:C1254–66.
 11.
Legant WR, Miller JS, Blakely BL, Cohen DM, Genin GM, Chen CS. Measurement of mechanical tractions exerted by cells in threedimensional matrices. Nat Methods. 2010;7(12):969–71.
 12.
JorgePeñas A, IzquierdoAlvarez A, AguilarCuenca R, VicenteManzanares M, GarciaAznar JM, Van Oosterwyck H, et al. Free Form Deformation–Based Image Registration Improves Accuracy of Traction Force Microscopy. PLoS One. 2015;10(12).
 13.
Dembo M, Wang YL. Stresses at the celltosubstrate interface during locomotion of fibroblasts. Biophys J. 1999;76(4):2307–16.
 14.
Tikhonov A. Solution of incorrectly formulated problems and the regularization method. Soviet Math Doklady. 1963;5:1035–8.
 15.
Sabass B, Gardel ML, Waterman CM, Schwarz US. High Resolution Traction Force Microscopy Based on Experimental and Computational Advances. Biophys J. 2008;94(1):207–20.
 16.
Han SJ, Oak Y, Grisman A, Danuser G. Traction Microscopy to identify force modulation in subresolution adhesions. Nat Methods. 2015;12(7):653–6.
 17.
Bohr J, SinglaBuxarrais G, Vincent R, Trepat X. Compressed sensing traction force microscopy. Acta Biomater. 2015;26:286–94.
 18.
A. SuñéAuñón, A. JorgePeñas, H. Van Oosterwyck, A. MuñozBarrutia, “L_{1}regularized reconstruction for traction force microscopy”. 13th IEEE International Symposium In Biomedical Imaging (ISBI); 2016. pp. 140–144.
 19.
Landau L, Lifshitz E. Theory of elasticity: course of theoretical physics, vol. 7. Oxford: Pergamon Press; 1970.
 20.
Candes EJ, Romberg J, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory. 2006;52(2):489–509.
 21.
M. Fornasier, “Theoretical foundations and numerical methods for sparse recovery”. vol. 9. Ed. Walter de Gruyter; 2010.
 22.
Scales JA, Gerztenkorn A, Treitel SJ. Fast L_{p} solution of large, sparse, linear systems: Application to seismic travel time tomography. J Comput Phys. 1988;75:314–33.
 23.
A. JorgePeñas, A. MuñozBarrutia, E. M. deJuanPardo, C. OrtizdeSolorzano, “Validation tool for traction force microscopy”, Comput Methods Biomech Biomed Eng. 2015;18(13):1377–85.
 24.
Ball GH, Hall DJ. ISODATA, a novel method of data analysis and pattern classification: Stanford research institute Menlo Park CA; 1965.
 25.
Alliney S. Digital filters as absolute norm regularizers. IEEE Trans Signal Process. 1992;40(6):1548–62.
 26.
Yin W, Chen T, Zhou SX, Chakraborty A. Background correction for cDNA microarray images using the TV+ L_{1} model. Bioinformatics. 2005;21(10):2410–6.
 27.
Gao H, Zhao H. Multilevel bioluminescence tomography based on radiative transfer equation Part 2: total variation and l_{1} data fidelity. Opt Express. 2010;18(3):2894–912.
 28.
K. G. Chan, S. J. Streichan, L. A. Trinh, M. Liebling, “Simultaneous temporal superresolution and denoising for cardiac fluorescence microscopy”, IEEE Trans Comput Imaging. 2016;2(3):348–58.
 29.
MuñozBarrutia A, Blu T, Unser M. L_{p}multiresolution analysis: how to reduce ringing and sparsify the error. IEEE Trans Image Process. 2002;11(6):656–69.
 30.
Chan TF, Esedoglu S. Aspects of total variation regularized L_{1} function approximation. SIAM J Appl Math. 2005;65(5):1817–37.
 31.
Chen T, Yin W, Zhou XS, Comaniciu D, Huang TS. Total variation models for variable lighting face recognition. IEEE Trans Pattern Anal Mach Intell. 2006;28(9):1519–24.
 32.
Schwarz US, Balaban NQ, Riveline D, Bershadsky A, Geiger B, Safran SA. Calculation of forces at focal adhesions from elastic substrate data: the effect of localized force and the need for regularization. Biophys J. 2002;83(3):1380–94.
 33.
DelanoëAyari H, Rieu JP, Sano M. 4D traction force microscopy reveals asymmetric cortical forces in migrating Dictyostelium cells. Phys Rev Lett. 2010;105(24):248103.
 34.
Peschetola V, Laurent VM, Duperray A, Michel R, Ambrosi D, Preziosi L, Verdier C. Timedependent traction force microscopy for cancer cells as a measure of invasiveness. Cytoskeleton. 2013;70(4):201–14.
 35.
R. Michel, V. Peschetola, G. Vitale, J. Etienne, A. Duperray, D. Ambrosi, L. Preziosi, C. Verdier, “Mathematical framework for traction force microscopy”. vol. 42. In ESAIM: Proceedings EDP Sciences; 2013. pp. 61–83.
 36.
Bergert M, Lendenmann T, Zündel M, Ehret AE, Panozzo D, Richner P, et al. Confocal reference free traction force microscopy. Nat Commun. 2016;7:12814.
Acknowledgements
Not applicable.
Funding
This work was partially supported by the Spanish Ministry of Economy and Competitiveness (TEC2013–48552C2–1R, TEC2015–73064EXP and TEC2016–78052R) (AMB, ASA) and (SAF2014–54705R) (MVM, RAC), the European Research Council (ERC) under the EUFP7/2007–2013 through ERC Grant Agreement n° 308,223 (HVO, AJP). ASA is supported by an FPI grant of the Spanish Ministry of Economy and Competitiveness. MVM is supported by a Marie Curie Grant (CIG293719) and a Ramon y Cajal fellowship (RYC2010–06094) from the Spanish Ministry of Economy and Competitiveness.
The authors declare that the funding bodies have not had influence in the design of the study and the datasets, the analysis and interpretation of the data and in writing the manuscript.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Author information
Affiliations
Contributions
ASA contributed in the conception of the work, the search of bibliography, the methods’ implementation, the design and implementation of synthetic data, the analysis and presentation of the results and was the major contributor in writing the manuscript. AJP contributed in the conception of the work, the search of bibliography, the methods’ implementation, the analysis of the results and writing and revising the manuscript. RAC and MVM contributed in the design and acquisition of the real data and revising the manuscript. HVO contributed in the conception of the work and revising the manuscript. AMB contributed substantially in the conception of the work, the search of bibliography, the design of synthetic data, the analysis and presentation of the results and writing and revising the manuscript. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Arrate MuñozBarrutia.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1:
Comparison of the different terms involved in the cost functional for the different regularization methods and a swept of the regularization parameter (λ) values. The first row is for Tikhonov regularization, the second one for L_{1}regularization and the third row is for full L_{1}regularization. The columns from left to right show: traction field (in Pa), displacement field (in μm), K ∙ t in Eq. 6 (in μm) and data fidelity term (Kt − u_{ q } in Eq. 6) (in μm). The outline of the mask used for traction recovery is shown in white. The scale bar represents 30 μm. (GIF 18689 kb)
Additional file 2:
Quantitative results from synthetic data. Table with the different error metrics (mean±standard deviation) obtained by the different regularization methods on the synthetic data. Ten different traction maps have been considered and ten realizations for each one of them. The best results for each metric are highlighted in red. For all cases, the differences between the realizations are statistically significant (p < 0.001) as computed by a Student’s test. (TIFF 431 kb)
Additional file 3:
Ten samples frames illustrating the whole TFM experiment and comparing the traction recovery with the different regularization methods. For each frame: (Top row, left) CHO cell expressing LifeactGFP; (Top row, center) Pseudocolor image showing the fluorescent beads at the hydrogel surface. The beads of the unstressed and stressed hydrogels have been superposed and pseudocolored in red and green, respectively; therefore, beads are colored in yellow when not displaced. The contrast of the pseudocolor images has been modified to highlight the areas with bead displacements; (Top row, right) Magnitude (in μm) and direction (arrows) of inplane displacements estimated from the bead images. (Bottom row) Recovered traction magnitude (in Pa) and direction (arrows) using: (Left) Tikhonov regularization; (Center) L_{1}norm regularization; (Right) full L_{1}norm regularization. The outline of the mask used for traction recovery is shown in red (Top row, right) and white (for the rest). The scale bar represents 30 μm. Frame #5 corresponds to Fig. 7 in the main manuscript. (GIF 10461 kb)
Additional file 4:
Force balance over real cells. Table with the mean (in Pa and in percentage of the maximum traction magnitude) and the standard deviation (in Pa) of the sum of forces over the whole cell for each regularization scheme and for all real dataset. (TIFF 110 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Traction Force Microscopy
 Spatial domain
 Regularization
 Spatial resolution