Validation tool for traction force microscopy

Traction force microscopy (TFM) is commonly used to estimate cells' traction forces from the deformation that they cause on their substrate. The accuracy of TFM highly depends on the computational methods used to measure the deformation of the substrate and estimate the forces, and also on the specifics of the experimental set-up. Computer simulations can be used to evaluate the effect of both the computational methods and the experimental set-up without the need to perform numerous experiments. Here, we present one such TFM simulator that addresses several limitations of the existing ones. As a proof of principle, we recreate a TFM experimental set-up, and apply a classic 2D TFM algorithm to recover the forces. In summary, our simulator provides a valuable tool to study the performance, refine experimentally, and guide the extraction of biological conclusions from TFM experiments.


Introduction
Tissue cells are anchorage dependent, i.e. must adhere to the extracellular matrix (ECM) to grow and survive (Galbraith & Sheetz 1998;Discher et al. 2005). Through their focal adhesions, cells exert forces on their surrounding environment, remodeling the ECM and promoting tissue morphogenesis (Ingber 2006). These traction forces are crucial in a number of biological processes including cell signaling and cell migration (Chicurel et al. 1998;Galbraith & Sheetz 1998). Therefore, accurate quantification of cell forces will help understanding the mechanisms underlying normal tissue morphogenesis and tissue remodeling both normal and diseased.
Traction force microscopy (TFM)  is widely used to quantify traction forces exerted by cells laying or embedded in a biomimetic hydrogel. This is done by measuring the displacement of fluorescent beads embedded in the gel, which serve as tracking points of the deformation caused by the forces. Any TFM experiment is affected by a number of experimental and computational variables that must be carefully and separately evaluated before extracting biological conclusions. This can be easily done in silico using computer simulations, thus avoiding expensive and time-consuming experimental repetitions. In addition, using simulations one can decouple the intrinsic computational errors of the method from those introduced by the non-ideal characteristics of the experimental setup. In this context, several quantitative studies (Schwarz et al. 2002;Sabass et al. 2008;Stricker et al. 2010) have revealed that in 2D, force recovery algorithms are very sensitive to errors in the calculation of the displacement of the beads, the density of measurements taken at focal adhesions, the distance between neighboring focal adhesions, and the amount of regularization used by the force recovery algorithms. Interesting as they are, these studies present two significant limitations. On the one hand, these simulations consider that cells only exert shear forces (parallel to the plane of the substrate), thus ignoring forces exerted in the third, vertical dimension. The importance of these forces on the interactions between the cell and its substrate is being increasingly acknowledged (Delanoë-Ayari et al. 2010;Hersen & Ladoux 2011;Legant et al. 2013). On the other hand, these simulations only take into account errors originated during the force recovery step, while the errors derived from the calculation of the displacements of the beads are simply approximated as a loss of spatial resolution plus an additive noise. However, since the quality of the calculated force field depends on the accuracy of the displacement estimation, it seems reasonable to separately analyze both sources of error.
To address these two relevant issues, we have developed a novel TFM simulator. This simulator generates 3D images of elastic hydrogels containing fluorescent beads, as acquired by a non-ideal optical microscope. Then, it creates an artificial cell, located on top of the hydrogel, exerting traction forces that in turn produce 3D displacements of the fluorescent beads. Using this set-up, any combination of bead displacement calculation and force recovery algorithms can be applied and analyzed under multiple experimental set-ups (i.e. density of beads, intensity, and distribution of forces). In particular, unlike other simulators reported in the literature, the accuracy of the displacement calculation can be analyzed independently of the precision of the force retrieval process. Specifically, we can study the effects of the size and separation of focal adhesions or the range of exerted forces, and decouple them from those of the size and density of the fluorescent beads.
As a proof principle of how our simulator can be used to analyze the results of a particular TFM algorithm and experimental set-up, we have simulated the classic 2D TFM algorithm that uses particle image velocimetry (PIV) to calculate the displacement field  and Fourier transform traction cytometry (FTTC) to recover the shear forces . The results were analyzed using a battery of error metrics, and relevant conclusions were extracted about the magnitude, orientation, and spatial distribution of the estimated forces.

Simulator
Our TFM simulator consists of two modules that can be used separately or combined sequentially. The first module simulates 3D forces generated by a synthetic cell, along with the displacement field caused by the forces in the surrounding medium; the second module generates images of hydrogel volumes containing embedded fluorescent beads, as acquired by an optical microscope. When both components are combined, the forces generated by the first module are fed into the second module, thus producing two images of hydrogel volumes, one in a relaxed state (no cell forces being exerted) and the other in an stressed state (under the deforming effect of the cell forces) as in a real TFM experiment (see Figure 1). We describe these two modules next.

3D force generator
To relate cell forces with the displacement field that they cause in the substrate, we implement the direct solution of the elasticity problem under the assumption of homogeneous, isotropic and linear materials where F i is the ith component (i [ fx; y; z}) of the force field F, D k is the kth component (k [ fx; y; z}) of the displacement field D, G is the Green function, and x ¼ ðx 1 ; x 2 ; x 3 Þ is a vector representing the ðx; y; zÞ spatial components.

Force field distribution
The simulator uses a synthetic mask of a cell that can be manually drawn or imported from the segmentation of a real cell. The location of the cell focal adhesions can be manually specified by the user as well. Alternatively, the system can pseudo-randomly locate them, after the user has defined their number and the minimum distance between them. The synthetic focal adhesions are of elliptical shape with their principal axis oriented toward the cell centroid. The magnitude of the forces is made proportional to the area of the focal adhesion as proposed by Balaban et al. (2001). The azimuth angle of the 3D forces points toward the centroid of the cell. Regarding the Z-component of the force and its inclination angle, the user can select between two different push-pull configurations (see Supplementary Note 1): pivot-like behaving cells with focal adhesions exerting rotational moments (Legant et al. 2013); or Dictyostelium-like cells (Delanoë-Ayari et al. 2010), where focal adhesions at the cell periphery pull upwards, and an area surrounding the cell centroid pushes downwards (see Figure 2).

Mechanical response of the hydrogel
We approximate the hydrogel as an elastic half-space, thus the Green function G can be described by the Boussinesq solution (Landau & Lifshitz 1970). It is well known that this solution presents a singularity at the origin. To avoid it, we use an approximation based on the integration of the Boussinesq solution over a small rectangular area. The equations for the 2D force and displacement fields were presented in Huang et al. (2009), and we have extended them to 3D (see Supplementary Note 2).

Displacement field calculation
We consider the hydrogel as a linear and shift-invariant system, being the Green function G the impulse response of the hydrogel. Therefore, when using the integral Boussinesq solution, the convolution in Equation (1) becomes a product in the Fourier domain and the displacement field can be easily computed componentwise where F and F 21 represent the forward and inverse Fourier transforms, respectively.

Hydrogel simulator
Hydrogel volumes with embedded fluorescent beads are simulated by uniformly distributing synthetic spheres in a 3D space. The beads are distorted as described in this section, simulating the effect of the ensemble microscope -CCD camera. Figure 3 shows an example of real and synthetic images of fluorescent beads embedded in the hydrogel, as acquired with a widefield microscope.

Effect of the microscope optics
We simulated the acquisition of 3D hydrogel volumes using two microscopy set-ups: a conventional microscope or a spinning disk confocal microscope. Conventional microscope. A microscope can be characterized by its impulse response, also known as point spread function (PSF) (Young 1996). The shape of the PSF is affected by increasing spherical aberrations as we move from the objective lens, deeper into the sample (Sibarita 2005). This is aggravated by the mismatch between the refractive indexes of the different media traversed by light (i.e. the objective, the immersion medium, the coverslip, and the ECM-like hydrogel). As a result, the PSF is not symmetric with respect to the focal plane and is non-stationary along the optical Z-axis. To account for this, our simulator implements the analytical PSF model proposed by Gibson and Lanni (1992). Following this model, the depth-variant PSF of the widefield microscope PSF z;WF is expressed by the Kirchhoff's diffraction integral formula as where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , J 0 ½: is the zero order Bessel function of the first kind, NA the objective numerical aperture, l em the emission wavelength of the fluorescent beads, OPD the optical path difference associated with the refractive index mismatch, and r is the normalized radius of the objective back focal plane.
Spinning disk confocal microscope: this microscope provides high contrast confocal quality images at reasonable speed. If the ratio between the excitation (l ex ) and emission wavelengths (l em ) is smaller than 2, the difference between the wavelengths has little effect on the PSF (Conchello & Lichtman 1994), which is the case of typical bead fluorochromes. Then, the PSF z;SD can be approximated as where^denotes the convolution operator and P is a function describing the pinhole aperture pattern of the spinning disk. These apertures are modeled by a Dirac delta function and distributed in a hexagonal pattern. A comparison between the PSFs of a widefield and spinning-disk systems is shown in Supplementary Figure 1 for a sample located at 0, 5, and 90 mms depth. Summarizing, our simulator models the acquisition of a volume containing the spherical particles I p as where PSF z is given by Equations (3) or (4) depending on the type of microscope.

Effect of the image acquisition system
First, the image obtained from the simulated microscope is resampled to fit the resolution of the CCD camera (XYplane) and the stage (Z-direction). This is done by downsampling I blurred over the Z dimension by a factor D z and integrating every XY plane over squared regions of size D xy . Second, three sources of noise are added to the signal: (1) photon noise inherent to the light signal of the fluorescent beads; (2) dark current noise caused by the thermal energy in the CCD camera; (3) readout noise associated with the on-chip camera amplification. The first two sources are characterized as Poisson processes (Aguet et al. 2005), while the third one follows an additive Gaussian distribution (Van Vliet et al. 1998). Therefore, the output of the bead simulator can be written as where R D xy ;D z denotes the resampling function in both XY and Z, and N p and N G are the Poisson (shot and thermal) and Gaussian (readout) noise distributions, respectively.

Implementation
The TFM simulator was built using custom-made code written in Matlab (The Mathworks Inc., Natick, MA, USA). All the simulator parameters can be tuned by the user. Given the size of the 3D data-sets and the subpixel resolution required for the calculation of bead displacements and recovered forces, we divided the data volumes into smaller, tractable blocks. This partition optimizes the use of memory, at the expense of higher time consumption.

Proof of principle
To show the use of our simulator, along with its potential benefits, we have simulated several artificial experimental set-ups, and we have applied a well-known 2D TFM algorithm based on PIV  and FTTC , to retrieve the shear forces at the surface of the 3D synthetic gel. It must be emphasized that the use of this proof of principle is not to present or analyze the algorithm itself, which is well known and has been analyzed elsewhere (Sabass et al. 2008;Stricker et al. 2010), but to illustrate how our simulator can be used to analyze the results of any TFM algorithm. Introducing novel, more advanced algorithms that solve the limitations of the 2D PIV -FTTC algorithm is far beyond the scope of this article, which focuses on the simulator itself. In the next paragraphs, we describe the algorithm that we used, the simulations performed, the metrics of error used, and the obtained results.

TFM algorithm
The 2D PIV -FTTC uses a PIV-like algorithm to calculate the displacement of the beads and recovers the forces using Fourier domain deconvolution (FTTC). The displacement field is calculated using local correlation windows, following a multiscale resolution strategy. In this study, 64 £ 64 and 16 £ 16 pixel-sized windows were used as the coarsest and finest resolution, respectively. Then, the displacement field is converted to the Fourier space and the force field is retrieved by inverting the Green function G in the Fourier domain. Since this function behaves as a low-pass filter, its inverse amplifies high frequency noise in the displacement field.
To overcome this problem, the computed displacement field was low-pass filtered before deconvolution.

Simulated data-sets
We performed 48 TFM experiments consisting in 12 different simulated cells showing a Dictyosteliumlike push-pull behavior, under four experimental conditions: widefield/spinning-disk microscopy and medium (0.33 beads/mm 3 ) or high (1 bead/mm 3 ) bead density. The 12 cells had different number and distribution of synthetic focal adhesions. The magnitude of the forces varied randomly between 5 and 65 nN. For each experiment, we generated images of two 90 mm thick hydrogels, one under stress and one relaxed. Only the first 6 mm from the hydrogel surface was simulated. The Young modulus, the Poisson ratio, and refractive index were set to 9 kPa, 0.5 and 1.39, respectively. The hydrogels contained 0.5 mm fluorescent beads (l em ¼ 605 ; nm) and were generated as imaged through a 40 £ dry objective lens (0.95 NA), captured using a CCD (8 mm £ 8 mm pixel size). The spatial resolution was set to 0.22 mm in the XY-plane and 0.3 mm along the Z-axis. The images were generated using five times the aforementioned values. Later, they are converted back to the system resolution to account for the signal degradation that occurs at the acquisition step, and to test the sub-pixel resolution capabilities of TFM algorithms. We handled the non-invariance generated by the depthvarying model of image formation by dividing the image into sub-volumes along the Z-axis, inside of which the PSF variation can be considered negligible (Preza & Conchello 2004) (details are in Supplementary Note 3). To emulate the methodology of a real 2D TFM experiment, we only used the top section of the simulated fluorescence bead volume, which represents the surface of the hydrogel, to calculate the bead displacement and recover the cell forces.

Error metrics
We quantified the error in the recovery of the magnitude and the direction of both the displacement and force fields. Furthermore, we have also measured the force spread beyond the focal adhesion.
Being V retrieved the retrieved vector field, V model the original vector field generated by the model, and N T the total number of points within the domain where both vector fields are compared, the error metrics are defined as described next.

Magnitude error
The error in the estimated magnitude was quantified using both the absolute error and its bias.
The mean absolute normalized error in magnitude (MANE m ) quantifies the relative absolute difference between V retrieved and V model where MANE m is a symmetric score, thus equally penalizes the errors due to overestimation and underestimation. The lower bound MANE m ¼ 0 indicates no errors, whereas the upper bound MANE m ¼ 1 implies an infinite error. We also calculate the mean bias normalized error in magnitude (MBNE m ) A negative value of MBNE m reveals underestimation whereas a positive value implies overestimation of the displacement or the force fields. We combined MANE m and MBNE m into a single error metric to facilitate the use and interpretation of the results. In this context, the mean percentage error in magnitude (MPE m ) is defined as where SF m being the scale factor between the recovered and original magnitudes, computed as

Direction error
The errors in the direction of V retrieved relative to V model were measured as the mean shortest rotation absolute angle (MSRAA) between both fields where a i is the shortest rotation absolute angle given by and f i is the direction angle of the i th component of V in degrees.

Sparsity of the recovered force field
The forces applied by the cell are confined to the focal adhesions. Consequently, the retrieved force field should ideally be sparse (i.e. a small number of non-zero entries contain a large proportion of the energy). The Gini index G is a measure of sparsity for which a value of 0 means uniformly distributed forces and a value of 1 expresses maximal inequality (i.e. all the force is exerted at one pixel). The sparsity ratio (SR ¼ G retrieved =G model ) measures how much the recovered forces are spread compared to the original. The Gini index is calculated as follows. First, the magnitude of the force field ( F k k) is arranged as a columnarray f ¼ ½f 1 ; f 2 . . . f N T and reordered from the smallest to largest, f ð1Þ # f ð2Þ # · · · # f ðN T Þ , where ð1Þ; ð2Þ; . . . ðN T Þ are the new indices after sorting. Then, the Gini index (0 # G # 1) is measured (Hurley & Rickard 2009)

Discussion
In this paper, we present a tool that simulates the 3D force fields generated by an adherent cell lying on top of an elastic hydrogel, along with the corresponding displacements suffered by the substrate. Any TFM algorithm can then be used to recover the forces from the displacement field. This way, different TFM algorithms can be compared and evaluated in silico against changes in the experimental set-up.
Our simulator produces cell forces of variable magnitude and distribution. It simulates the displacement of fluorescent beads as they would be observed under a real microscope, including the resolution loss caused by the CCD acquisition and the noise and the blurring introduced by the optical system.
Our TFM simulator improves the existing solutions (Schwarz et al. 2002;Sabass et al. 2008) in two ways: on one hand, our TFM model allows creating 3D force fields, which have been reported to be essential in the study of the cell -substrate interactions (Delanoë-Ayari et al. 2010;Legant et al. 2013). On the other hand, the fluorescence bead simulator allows evaluating the performance of TFM algorithms under a variety of conditions, including the size and distribution of focal adhesions, the magnitude of the bead displacements, the influence of the bead density, or the effects caused by the optical system. Nevertheless, our simulator presents some limitations that have to be considered. First, we have characterized the Green function of the hydrogel by a solution based on Boussinesq equations which assumes linear elasticity. Although this has no effect on the evaluation of the capability of the TFM algorithms to recover the displacement and force fields, the assumption of linearity has been recently questioned, and a solution based on FEM has been proposed for more complex hydrogel models (Palacio et al. 2013). In addition, due to practical reasons, the simulator assumes a correlation between the force magnitude and the size of the focal adhesion. However, it is still not clear whether this correspondence holds in real cells for all focal adhesions (Stricker et al. 2011).
As a proof of principle, the classic 2D PIV -FTTC has been evaluated with 3D synthetic data generated by our simulation tool under different conditions, including two bead densities and widefield/spinning-disk microscopy set-ups. Three complementary error metrics have been used to evaluate the TFM results with regard to the orientation and the magnitude of the computed vector fields, and the spread of the focal adhesions. As clearly shown in the representative Figure 4 and confirmed quantitatively (Table 1), 2D PIV -FTTC tends to underestimate retrieved forces and spread them out in larger areas, merging the individual contributions of closely located focal adhesions. Both errors are caused by the PIV algorithm, where the window correlation method used introduces a resolution loss as large as the window size and an averaging of the displacements inside each window. Specifically, the resolution loss has a considerable negative impact on the results by preventing force recovery algorithms from retrieving small forces and from confining the force into closely located individual focal adhesions. In addition, it should be emphasized that 2D PIV -FTTC assumes negligible vertical forces, which contribute as well to the underestimation of recovered total force magnitude. Finally, it is important to note that when using spinning disk confocal microscope, the results improve with increasing density of fluorescent beads. This is not the case when using a widefield microscope, maybe due to the background signal generated by the out-of-focus light, which increases with the number of beads and complicates the calculation of the displacement field.

Conclusion
We have developed a computer simulator that mimics all the steps involved in a real TFM experiment and can be used to study and compare in silico the performance of TFM algorithms. This is a valuable tool to interpret the results obtained using existing and novel algorithms on different experimental set-ups. Ultimately, this information will be helpful to support the consistency of biological findings extracted from real experiments.

Supplemental data
Supplemental data for this article are available.