## Abstract

Adherent cells exert tractions on their surroundings. These tractions can be measured by observing the displacements of beads embedded on a flexible gel substrate on which the cells are cultured. This paper presents an exact solution to the problem of computing the traction field from the observed displacement field. The solution rests on recasting the relationship between displacements and tractions into Fourier space, where the recovery of the traction field is especially simple. We present two subcases of the solution, depending on whether or not tractions outside the observed cell boundaries are set to be zero. The implementation is computationally efficient. We also give the solution for the traction field in a representative human airway smooth muscle cell contracted by treatment with histamine. Finally, we give explicit formulas for reducing the traction and displacement fields to contraction moments, the orientation of the principal axes of traction, and the strain energy imparted by the cell to the substrate.

- displacement fields
- fluorescent beads

cells exert tractions on their surroundings in the course of a variety of cell functions including contraction, spreading, crawling, and invasion. These functions are associated with complex mechanical interactions between the substrate, adhesion molecules, cytoskeletal elements, and molecular motors. Dembo and Wang (3) recently have shown that the traction field that a cell exerts on its surroundings can be mapped from knowledge of the displacement field in a flexible substrate on which cells are adherent. Measurement of the displacement field is accomplished by tracking small beads, typically 0.2 μm in diameter, embedded near the surface of the substrate gel. Computing the tractions from the measured displacement field is difficult and computationally intensive, however, so much so that Pelham and Wang (5), for example, suggested using the raw displacements themselves as a qualitative map of the local tractions. More recently, Balaban et al. (1) implemented a simplified version of the Dembo and Wang (hereafter denoted DW) approach, which is restricted to the case of isolated point sources of traction.

In this article, we describe the foundations of Fourier transform traction cytometry (FTTC), a new and computationally efficient method for computing the traction field given the displacement field. FTTC divides into two subcases (unconstrained and constrained), described below, and in several respects is fundamentally different from the DW method. The first key difference is that the DW method requires the cell boundary to be drawn by hand and constrains the tractions exterior to the boundary to be zero while still retaining approximate matching of the exterior displacements. By contrast, unconstrained FTTC exactly matches all the displacement data, independent of the perceived cell boundary; constrained FTTC forces the tractions exterior to the cell boundary (drawn by hand) to be zero, as in the DW method, but ignores the exterior displacement data. The second key difference is that the DW approach utilizes Tikhonov regularization (8, 9) with a particular choice and intensity of smoothing functional (6,10). By contrast, FTTC utilizes no smoothing and is exact in the sense that it yields a traction map for which the induced displacements exactly match the given displacement field. The discussionelaborates these issues.

We present the explicit formulation and solution to the traction recovery problem and introduce the concept of contraction moments and strain energy, which are particularly robust measures that characterize the cell-substrate interaction. Finally, we include one example of the traction field associated with a human airway smooth muscle cell and an illustration of the displacement and traction fields associated with a particularly simple simulation. In two companion studies (7, 11), we have utilized FTTC to address the biologically important questions of the relationship between cell prestress and cell rigidity, and the extent to which microtubules play a role in this connection.

## THEORY

### The Boussinesq Solution

The traction field is defined as the stress, i.e., local force per unit area, imposed on the gel surface by an adherent cell. The traction field, in turn, determines the displacement field of the gel surface. If the gel can be approximated as a semi-infinite solid (see*Limits of applicability*), the displacements can be computed from the distribution of surface tractions as follows. First, one finds the displacement field, or Green's function, associated with a point traction on the surface. This is a classic problem, the solution to which was found by Boussinesq (see Ref. 4). Second, integration of this function over the given traction field then yields the corresponding displacement field; this is the so-called forward problem. The problem we address in this article is the inverse problem, namely, inferring the traction field from measured displacements.

The FTTC method is based on Fourier analysis and arises from the observation that the displacement at any point *r⃗* on the surface due to a point traction source at another point*r⃗*′ is (apart from the direction of the displacement and the direction of the traction) a function only of the difference ‖*r⃗* − *r⃗′*‖. We denote the displacement vector at *r⃗*′ as*u⃗*(*r⃗*) and the traction vector at*r⃗*′ as *T⃗*(*r⃗*′). We denote the Green's function, or kernel, mapping traction to displacement by the tensor **K** =**K**(‖*r⃗* − *r⃗′*‖). The displacements are then given by the convolution *u⃗*= **K** ⊗ *T⃗*. In this representation,*u⃗*(*r⃗*) and*T⃗*(*r⃗*′) are 2-vectors with elements labeled *x* and *y* (we ignore displacements in the*z* direction, and the traction in the *z* direction, or normal stress, is taken to be 0), the kernel **K** is a 2 × 2 matrix, and ⊗ denotes integration over*r⃗*′. The major difficulty in the inversion of this equation is that **K** is not diagonal in *r⃗* and*r⃗*′ (if it were, then the solution would only involve inverting a 2 × 2 matrix). The fact that **K** is not diagonal in real space (i.e., tractions at one point are coupled to displacements at different points) is the origin of why its inversion in real space necessarily requires the construction and inversion of very large matrices (as in the DW approach). In Fourier space, these difficulties do not arise.

#### The Boussinesq solution is diagonal in Fourier space.

The key to the FTTC method is the exploitation of the Faltung or convolution theorem, which states that the Fourier transform of a convolution is the simple product of the Fourier transforms of the functions convolved. The forward problem then becomes
=**K̃**(*k⃗*)
, where the tilde overbar denotes the (two dimensional) Fourier transform with wave vector *k⃗*. In this form it is clear that **K̃** is diagonal in that there is no coupling between different wave vectors *k⃗*. Of course,**K̃** remains a nondiagonal 2 × 2 matrix insofar as tractions in the *x* or *y* direction separately induce displacements in both the *x* and *y*directions, but this presents no difficulty because **K̃**remains strictly diagonal in *k⃗* space. It follows that**K̃**
^{−1} is trivial to compute if**K̃** is known. The solution to the inverse problem is then given by
Equation 1where FT
denotes the (two dimensional) inverse Fourier transform.

#### Explicit evaluation of the kernel in Fourier space.

Implementation of *Eq. 1 *requires an explicit formula for**K̃**(*k⃗*), the Fourier transform of the Boussinesq solution **K**(*r*), where *r*= ‖*r⃗*‖. The forward kernel, written in matrix form, for a point source at the origin is given by Landau and Lifshitz (4), which when reduced to *x* and *y*displacements with zero normal traction is given by
Equation 2where *A* = (1 + ς)/π*E*, in which ς is Poisson's ratio and *E* is Young's modulus. The components of this matrix are denoted *K _{i,j}
*, where the indices

*i*and

*j*run through (

*x, y*); this notation will be used consistently with all matrices and vectors. Thus the two-dimensional Fourier transforms of

*r*

^{−1},

*x*

^{2}/

*r*

^{3},

*xy*/

*r*

^{3}, and

*y*

^{2}/

*r*

^{3}are required. These can be derived in several different ways. One especially clear method can be sketched as follows; it relies on a generalization of the problem to three dimensions and a subsequent reduction to two dimensions. We begin with the singular solution to Laplace's equation in three dimensions, ∇

^{2}

*r*

^{−1}= −4πδ

^{3}(

*r⃗*), where δ is the Dirac delta function. We denote the three-dimensional Fourier transform by FT

_{3}, which when applied to the Laplace equation yields FT

_{3}(

*r*) = 4π/(

^{−1}*k*+

*k*

^{2}), where we separated out the wave vector

*k*from the wave vectors in the

_{z}*x, y*plane (here

*k*

^{2}=

*k*+

*k*, and in what follows,

*k*is the nonnegative square root of

*k*

^{2}). Now note that the two-dimensional Fourier transform is the single inverse three-dimensional transform in

*z*, evaluated at

*z*= 0, i.e. The desired answer reduces to an elementary integral The other three transforms can be obtained by simple manipulations. Consider the expression ∂/∂

*k*FT

_{x}_{2}[(∂/∂

*x*)

*r*

^{−1}]. This can be evaluated in two distinct ways. First, directly differentiating with respect to

*x*and

*k*yields Second, integrating by parts with respect to

_{x}*x*yields where we use the transform of

*r*

^{−1}obtained above. Performing the indicated derivative with respect to

*k*and setting the two independent evaluations equal, we obtain FT

_{x}_{2}(

*x*

^{2}/

*r*

^{3}) = (2π/

*k*

^{3})(

*k*

^{2}−

*k*) = 2π

*k*/

*k*

^{2}. By symmetry, the transform of

*y*

^{2}/

*r*

^{3}can be written down by inspection, and the transform of

*xy*/

*r*

^{3}can be obtained by evaluating ∂/∂

*k*FT

_{y}_{2}[(∂/∂

*x*)

*r*

^{−1}] in the same two distinct ways as above. In summary, the desired transformed matrix is given by Equation 3

#### Contraction moments.

The Fourier approach also gives robust measures of certain low order moments of the tractions. The zeroth order moment of the tractions is given by and is equal to the net force applied by the cell to the substrate. For isolated adherent cells, this is known a priori to be zero. However, registration shifts from one image to another of a given pair (say before and after contractile activation) will induce a spurious traction field corresponding to a non-zero net force. This artifact can be trivially accounted for in Fourier space simply by setting (0) = 0, which guarantees no net force by the cell on the substrate.

The first order moments are associated with contraction/dilation tractions (radially oriented tractions) and tractions corresponding to torques (circumferentially oriented tractions). These correspond to the four combinations of the tractions in the *x* and *y*directions weighted by their *x* and *y* coordinates. For example, a positive traction in the *y* direction (*T _{y}
*) at a location with a positive

*x*coordinate (or a negative

*T*at some

_{y}*x*< 0) will contribute to a counterclockwise rotational torque through a term proportional to

*xT*. As a visual aid, the four terms are shown schematically in the following diagrammatic representation

_{y}When symmetrized (since the net torque conferred by the cell on the substrate must be zero) and integrated over the surface, this is the contraction/dilation and shear moment matrix **M**. Written in component notation, this matrix is explicitly given by
Equation 4
We approximate the derivatives by discrete differences in*k* space. Because
(0) = 0 by construction (to eliminate registration artifact), this expression then involves only the Fourier transforms of the tractions at the lowest non-zero wave numbers, Δ*k _{x}
* and Δ

*k*. Explicitly,

_{y}*Eq. 4*reduces to Equation 5The interpretation of the elements of the moment matrix,

*M*, is as follows. The total contribution of the cell to contracting the substrate in the

_{ij}*x*and

*y*directions is given by

*M*and

_{xx}*M*, respectively. M

_{yy}_{xy}(or

*M*) is the contribution of the cell to deformation of the substrate arising from variations in

_{yx}*x*tractions with

*y*and variations in

*y*tractions with

*x*. There are additional anisotropic contributions arising from unequal variations of

*x*tractions with

*x*and from

*y*tractions with

*y*, i.e., when

*M*≠

_{xx}*M*. A simple way to characterize this is to apply a rotation operator

_{yy}**R**to

**M**such that

**M**

^{rot}=

**R**

^{−1}

**MR**is diagonal, i.e.,

*M*=

*M*= 0. This form puts all the tractions of the cell into their principal axes. The orientation of the principal tractions can be obtained from the

*x*,

*y*coordinate axes of the original images and the angle of rotation of

**R**. This orientation of the cell is then clearly independent of the coordinate system, and may be an important measure of directionality in cell motility assays, including chemotaxis. In this context, the ratio of the principal tractions is a direct measure of the traction polarity of the cell.

The net moment tending to dilate or contract the substrate is given by the trace of the moment matrix; we thus define the net contractile moment μ of the cell by
Equation 6(Here either **M** or **M**
^{rot} can be used because the trace is invariant under coordinate rotations.) The net contractile moment μ is a coordinate invariant scalar measure of the cell's contractile “strength.”

#### Strain energy.

The total energy *U* transferred from the cell to the elastic distortion of the substrate is given by
Equation 7and is another measure of contractile strength. The use of Fourier analysis may involve some artifactual behavior at the field boundaries, which is particularly pertinent in computing the strain energy. We therefore evaluate this integral over the strict interior of the field, i.e., without the boundary points. Because of this, the value obtained is different from the evaluation of *Eq. 7 *in Fourier space with Parseval's theorem. The source of this discrepancy lies in the use of Fourier analysis over a finite domain, wherein periodic boundary conditions are imposed. In general, the displacements at the edges of the field of view will not approximate continuous periodic functions, and therefore there will be artifactual tractions present along the bounding nodes of the lattice grid. These in turn will contribute to the estimated strain energy if they are included in the integral above. Such artifacts also are present in the Fourier domain, although they are in general spread over all Fourier components. It is therefore simpler to avoid this problem by direct integration of the strain energy density over the strict interior of the domain.

#### Limits of applicability.

Note that both the FTTC method presented here as well as the DW method and that of Balaban et al. (1) explicitly approximate the elastic gel as a semi-infinite medium. It has been stated that this is valid if the displacements are small compared with the gel thickness (3), but this is not correct (Wilson TA, personal communication). In fact, the ratio of displacements to gel thickness being small is a necessary (but not sufficient) condition for the applicability of linear elasticity theory. By contrast, the use of a semi-infinite elastic continuum to approximate the finite thickness gel is valid to the extent that the lateral dimensions of the cell and the lateral distances over which displacements are measured are both small compared with the gel thickness. A simple example will illustrate this. Consider the case of uniform traction *T* over a circular region of radius *R*, on top of an incompressible slab of thickness *h*, shear modulus *G*, and with a fixed bottom. If *R *≪ *h*, then the gel is effectively infinitely thick, the Boussinesq solution applies, and the displacement of the disk is approximately *RT*/*G*. By contrast, if *R *≫ *h*, then the medium is approximately in simple shear, and the displacement of the disk is approximately*hT*/*G*. This implies that, in the latter case, the displacements will be lower than would be observed in the semi-infinite medium, which in turn leads to an underestimate of the tractions.

## IMPLEMENTATION

The implementation of all methods of computing tractions from bead displacements naturally divides into two parts. The first is the estimation of the displacement field itself on some appropriate lattice or mesh, and the second is the computation of the traction field for that given displacement field. This section describes these two processes.

### Estimating the Displacement Field From the Images

To estimate the displacement field of the substrate, we compared digital images of the same region of the gel, taken at different times. Images showed fluorescent microbeads (0.2 μm) embedded in the gel, before the cells were plated. In our experiments, there were usually 1,000–2,000 beads in an image. Images were of the size 1,024 × 1,280 pixels. Each bead image occupied an area of 4–6 pixels in diameter, and neighboring beads were about 5–30 pixels apart.

The processing of these images began with the correction of the pair of images for relative translational shifts. Here we used the correlation theorem, which says that the Fourier transform of a correlation of two functions is a product of the Fourier transform of one function and the complex conjugate of the Fourier transform of the other function. We formed the normalized two-dimensional cross-correlation function between the two images (normalized by the square root of the product of the maximal values of the autocorrelation functions of these 2 images). We identified the coordinates of the peak of the correlation function and translated one of the images with respect to the other by that uniform displacement. For the calculations of the correlation functions, we utilized the two-dimensional fast Fourier transform (FFT) algorithm in MATLAB. Even though the images were represented by relatively large (1,024 × 1,280) matrices, use of the correlation theorem and FFT algorithm made the actual computations reasonably fast.

Having the corrected images, we divided them into a number of small window areas. For these images, we chose a window size of typically 64 × 64 pixels. The window areas overlapped; the distance between the centers of successive windows was chosen to be 16 pixels. The displacement of each window area was then calculated by correlating a window in one image with the window at the same coordinates in the other image, in the way described above for the correlation between the whole images. The coordinates of the peak of the cross-correlation function between two windows were assigned to the center of the window as the window's displacement vector. Repeating the same procedure over all windows yielded the uniform discretized displacement field between two entire images. (The fact that this technique yields a lattice with uniform spacing means that simple FFT algorithms can be used in the subsequent analysis.)

The window size of 64 × 64 pixels was chosen to guarantee that at least one fluorescent marker was located within the window, regardless of the window position. To evaluate the potential smoothing effect of widow size on the recovered traction field, we also experimented with smaller windows down to 16 × 16 pixels. This required us to manually eliminate windows that did not include fluorescent markers and to substitute the missing displacements with those obtained from runs that employed larger windows. In the cell that we chose as an example here (see Figs. 3-5), smaller window sizes did not appreciably alter the recovered traction field.

Images of different gels differed in bead number and distribution. For images with a relatively low number of beads or less uniform bead distribution, there were window areas for which the computed cross-correlation was below a threshold that we set at 0.95. For the displacement of such window areas, we used values obtained by fitting the rest of the displacement field by a third-order polynomial and calculating the values of the polynomial in the points of the lattice for which the value was missing.

Our method of obtaining the displacement field is different from the DW approach, which relies on the measurement of the*x*, *y* coordinates of beads in each of the two images. We encountered two difficulties with that method. First, there is the issue of determining bead identity; the displacement requires the initial and final positions of the same bead. This is not a problem if bead density is low, but in attempting to achieve higher resolution with higher bead densities, some bead identities can become ambiguous; errors here can lead to spurious displacements and, hence, artifactual tractions. Second, there may be areas on some images where beads are sparse or even absent. In this case, no estimates are possible save by interpolation from neighboring regions. Our method of maximizing the cross-correlation of small windows between the two images is less sensitive to these problems. In the first place, if there are beads with unambiguous identity, their contrast dominates the cross-correlation function and the estimated displacements are similar to those computed by direct position measurements. However, ambiguities in any one bead identity contribute less to the displacement field to the extent that other beads are present in the same window. This is important in areas where the bead density is high and where there may be clusters of beads with ambiguous identities. Areas with no beads also can show reasonable correlation between the two images, depending on the displacements of other features that still carry sufficient contrast to be measured. Such features may include heterogeneities in the gel and embedded beads that are out of focus. In summary, the advantages of the cross-correlation approach are that it permits semiautomated estimates of the displacement field and is insensitive to ambiguities regarding bead identification between images.

### Computing the Traction Field From the Displacement Field

We have implemented the solution for computing cell tractions in two distinct ways. The first method, unconstrained FTTC, uses all displacement data from an image pair obtained as described in*Estimating the displacement field from the images*, does not use any constraints on the recovered tractions, and is a direct application of the methods described in theory. The second method, constrained FTTC, is the solution to the mixed boundary value problem, which ignores the displacement field outside the boundary of the cell and constrains the tractions outside the cell boundary to be zero. It is important to note that this method requires additional information beyond the displacement field, namely, an independent estimate of the location of the cell boundary, drawn by hand. As described below, there are advantages and disadvantages to both methods; they should be viewed as complementary approaches.

#### Unconstrained FTTC.

Here we use the direct solution given by *Eqs. 1
* and *
3
*. The specific procedure is as follows. *1*) Calculate the Fourier transform of the discrete displacement field (and set the Fourier component at *k⃗* = 0 to zero to eliminate translation artifact). *2*) Multiply the transformed displacements by**K̃**(*k⃗*)^{−1} to map the transformed displacements to transformed tractions. *3*) Take the inverse Fourier transform of the result to obtain the tractions.

#### Constrained FTTC.

This is a mixed boundary value problem, with the displacements under the cell being specified (by measurement) and with the tractions outside the cell boundary specified (by assumption) to be zero. The Fourier approach above also can be used iteratively to solve this problem. It consists of the following procedures. *1*) Calculate the traction field in the way described in*Unconstrained FTTC. 2*) Define a new traction field by setting the tractions outside of the cell boundary to zero.*3*) Calculate the displacement field induced by this traction field. This is done by using the Fourier approach in a forward direction: calculate the FT of the traction field, multiply by**K̃**(*k⃗*) to obtain the transformed displacements; the inverse FT is the new displacement field.*4*) Define a new displacement field by replacing the displacements of the calculated displacement field within the cell boundary by the experimentally observed displacements. *5*) Repeat *steps 1–4 *until convergence is reached at some level of tolerance. There are a variety of criteria that can be used. In our case we chose to terminate the iterative procedure when the variation in the maximum magnitude of the tractions within the cell was less than 1 part in 10^{6} on succeeding steps.

Note that in both constrained and unconstrained FTTC, there is an ambiguity in the off-diagonal elements of**K̃**(*k⃗*)^{−1} at the Nyquist frequency, because positive and negative Nyquist frequency components are indistinguishable but *k _{x}k_{y}
* ≠ −

*k*. This problem is avoided by setting the off-diagonal elements of

_{x}k_{y}**K̃**(

*k⃗*)

^{−1}to be zero when either

*k*or

_{x}*k*is a Nyquist frequency.

_{y}## EXPERIMENTAL METHODS

A technique for preparation of polyacrylamide gel sheets (3) was modified and used to make flexible gel disks. A mixture of acrylamide (2%), bis-acrylamide (0.25%), and fluorescent latex beads (diameter 0.2 μm, 1:125 dilution by volume) was added to activated glass coverslips. The droplet of the solution was covered by a small circular coverslip. After polymerization (45 min), the circular coverslip was removed. Type I collagen was attached to the surface of the gel. Gel disks were typically 50–70 μm thick and had a diameter of 12 mm. The elastic modulus (Young's modulus) of the gel was determined to be 1,200 Pa; Poisson's ratio was taken to be 0.48.

HASM cells were cultured in plastic dishes and serum deprived for 2 days before the experiments. Cells at *passage 3–6 *were plated on the gel disks in a serum-free medium and allowed to spread and stabilize for 6 h. Cells were then stimulated with histamine (0.01 mM) for 5 min. Photomicrographs were taken of the cells both with phase-contrast optics to visualize the cells and with 470-nm ultraviolet illumination to excite the beads, which fluoresce at 515 nm. To assess the distribution of beads (and other features with contrast) in the unstressed gel, the cells were detached from the substrate with trypsin (∼2%); this therefore leaves the flexible gel with no surface tractions. [Note that the image of the traction-free gel (i.e., “pretreatment”) is taken after the images of the posttreatment distribution of beads.]

## RESULTS

Figure 1 shows a phase-contrast image of a representative HASM cell, cultured on the flexible polyacrylamide gel covered with collagen type I, prepared as described in experimental methods, after histamine treatment.

Figure 2 shows a fluorescence image of the same field of view as in Fig. 1; the 0.2-μm beads embedded in the gel are easily visualized. Superposed on this image is the outline of the cell boundary, drawn by hand from Fig. 1. (This outline is drawn somewhat larger than the appearance of the cell in Fig. 1. This is to ensure the inclusion of potential stress bearing interactions between the cell and the substrate that may not be visible and to avoid interactions between the cell boundary and the discretized lattice on which the solution is defined. See discussion.) The corresponding fluorescent image of the beads after cell detachment with trypsin (pretreatment condition) looks virtually indistinguishable from this picture because the actual bead displacements are very small.

Figure 3 shows the displacement field computed from the two fluorescent images of the beads pre- and posttreatment, as described in *Computing the displacement field from the images*. The arrows in Fig. 3 show the relative magnitude and direction of the displacement field of the gel under the adherent smooth muscle cell. Figure 3 also is color coded by the absolute magnitude of the displacements.

Figure 4 shows the traction field as computed by unconstrained FTTC, the direct computation of tractions from the Fourier decomposition of the displacements. Also shown is the boundary of cell, although it is important to note that this information was not used in computing the tractions. The arrows in Fig.4 show the relative magnitude and direction of the tractions, and the colors show the absolute magnitude of the traction vectors.

Figure 5 shows the traction field calculated by constrained FTTC, iterating the Fourier approach until convergence was reached. Note that the tractions are zero outside the cell boundary by construction. As in Fig. 4, the arrows show the relative magnitude and direction, and the colors show the absolute magnitude of the traction vectors. Notice first that the maps in Figs.4 and 5 are similar but that there is a traction concentration near the boundary of the cell computed by constrained FTTC (Fig. 5). This results from the requirement, by construction, that all tractions exterior to the cell must be zero.

Table 1 displays the moment matrices**M** and **M**
^{rot}, the net contractile moment μ, the orientation of the principal tractions, and the strain energy *U* exerted by the cell as computed by unconstrained FTTC. Note that the moments are given in units of picoNewton meters (pNm) and energy is given in units of picoJoules (pJ). We use these units to distinguish clearly between moments (forces multiplied by distances from the origin) and energy (forces multiplied by displacements), despite the fact that pNm and pJ are formally equivalent. Table 2 displays the same quantities as computed by constrained FTTC, with the cell boundary drawn by hand as shown in Fig. 2. All non-zero tractions are constrained to lie within this boundary.

Note that there are two methods for obtaining the moment matrices, given by *Eqs. 4
* and *
5
*; the former is obtained by integration of the recovered tractions over the (real) *x-y*space, whereas the latter is obtained by the formal equivalent of the vector derivative of the tractions in Fourier space evaluated at the origin of *k* space. These two expressions in practice can be quite different. In principle, given that the moment matrix is defined by an integration over real space, it might be thought that using*Eq. 4
* directly would be preferable, but the existence of boundary tractions associated with the Fourier decomposition in the first place implies a potentially substantial error from the field boundary, especially in unconstrained FTTC. By contrast, the use of the first non-zero Fourier coefficient (*Eq. 5
*) as a measurement of the moment matrix has the advantage of multiplying the entire field by the lowest frequency sine component, thus exactly canceling the major dipole artifact arising from the boundary of the field of view. This is the method we use.

As an aid to interpreting the relationship between the traction map and moments, we show in Fig. 6, *A*and *B*, the displacement field and traction field for an artificially constructed example. The simulation consists of two pairs of point traction sources of different magnitudes, scaled to be similar to those seen in real cells. Note that the traction map has non-zero tractions only at the four source points, whereas the displacement field is non-zero over broader regions. (Because this is a simulation with no artificially added noise, the recovery of the traction map was identical between the unconstrained and constrained FTTC methods.) For this simulation, the moment matrices, orientation of principal tractions, and strain energy are listed in Table3. This shows how the off-diagonal elements of **M** arise when the laboratory coordinate system does not coincide with the principal axes and how when suitably rotated (here by ∼30° from the *x*-axis), the moment matrix becomes diagonal (**M**
^{rot}). The fact that*M*
and*M*
are both negative means that, as expected from Fig. 6, the principal tractions are contractile. The magnitude of *M*
is larger than that of *M*
, corresponding to stronger contraction along that axis.

### Effect of Noise in the Displacement Data

With respect to questions of resolution and accuracy in the presence of noisy data, it is important to recognize that in the FTTC approach, noise arises only in the estimation of the displacement field from the image pairs; the calculation of the traction field is, apart from round-off errors in finite word-length arithmetic, an exact procedure. There are many different methods by which to estimate the displacement field from image pairs; the one we have chosen (locating peaks in the cross-correlation function on windowed regions of the two images) is convenient for our purposes, but it is not the focus of this paper. We have, however, characterized the effect of displacement noise on the recovered moment matrices, net contractile moment, and strain energy, as described below.

We performed simulations where the displacement field consisted of pure noise and examined the departures of the recovered moments and strain energy from zero. We performed 100 simulations using the same grid and gel characteristics as in the real images in Figs. 2 and 3. The displacement noise was Gaussian with a mean of zero and a standard deviation of 1 μm. The traction field was recovered by using both the unconstrained and constrained Fourier methods, where in the constrained case we used the same cell outline as in Fig. 3. The results of these simulations on the moment matrices, the net contractile moment, and the strain energy are as follows. The elements of the moment matrices were not significantly different from zero. This was to be expected, because the tractions are linear functions of the displacements and the expectation values of the noise in the displacements are zero. Because this is true of all elements of the moment matrices, the effect of the noise on the net contractile moment also is not significantly different from zero. In any given single realization, however, it is important to know the expected magnitude of the departure from zero. The standard deviation of the net contractile moment arising from pure noise, from these simulations, is 0.29 and 0.21 pNm per micrometer standard deviation of displacements in the unconstrained and constrained case, respectively. The effect of noise on strain energy is quite different; unlike the tractions (and therefore the moments), which are linear functions of displacement, the strain energy is quadratic. This implies that the expectation value of the strain energy due to pure noise is non-zero. In our simulations, we found that the energy associated with displacement noise is 11.1 and 0.88 pJ per square micrometer of displacement variance in the unconstrained and constrained case, respectively. This substantial difference in strain energy in the constrained and unconstrained cases (roughly a factor of 12) is precisely what was expected, because the total field area is roughly 12 times the area bounded by the cell (Fig. 3), and there is no strain energy in the gel conferred by surface tractions exterior to the cell in the constrained case.

## DISCUSSION

### Advantages and Disadvantages in Unconstrained and Constrained FTTC

There are a number of advantages (+) and disadvantages (−) associated with unconstrained and constrained FTTC.

#### Unconstrained FTTC.

(+/−) All of the observed data are used, including the displacements of beads exterior to the perceived cell boundary. The falloff in displacements exterior to the cell constitutes additional information regarding the overall traction field, especially when the distribution of beads is sparse within the cell boundary. This may seem to be an obvious advantage; it is in large measure true. However, the falloff in displacements from any particular traction source point is like 1/*r*, so the information “content” of displacements exterior to the cell is correspondingly low, and important information may in fact not be lost.

(+) The cell boundary need not be identified, and, as such, no investigator judgment is necessary to identify this boundary. This is an advantage over all constrained methods insofar as force generation associated with small filopodia, flat lamellipodia, or fibrous connective tissue elements may be missed in the original micrograph images.

(−) As with all discrete Fourier problems implemented over a finite space, the inherent periodicity introduces artifactual tractions at the boundary of the field because the measured displacements are not strictly periodic. To the extent that these are remote from the cell, they pose no difficulty. Moreover, such an artifact is equivalent to a dipole field on the boundary, which does not strongly influence the computed cell tractions. These boundary artifacts are easy to recognize and can be safely ignored.

(+) Errors in the recovered tractions exterior to the real cell boundary, secondary to noise in the displacement field, will have zero mean if the noise has zero mean. This follows from the linearity of relationship between tractions and displacements.

(−) By contrast, the strain energy (which is quadratic in the displacements) will be artifactually high due to the contribution of noise in the traction field exterior to the cell.

#### Constrained FTTC.

(−) If the cell is exerting non-zero tractions on the substrate at locations exterior to the perceived cell boundary, the interior tractions will necessarily be in error, because they must compensate for these tractions that were incorrectly constrained to be zero. This results in artifactually large tractions, especially in the vicinity of the imposed cell boundary. The danger here is that the large tractions at the perceived cell boundary may be interpreted mistakenly as reflecting real traction concentrations.

(+/−) Constrained FTTC requires an iterative approach, and so computational efficiency is not guaranteed. In our experience, however, we have found that the iterative scheme described above typically converges quite quickly (typically <10 iterations) so that this approach is also not computationally intensive.

(−) The nature of the mixed boundary condition implies that an assessment of the noise on the final traction field is difficult. This is because the induced level of noise depends on the boundary of the cell.

(+) The strain energy in constrained FTTC is confined, by construction, to the interior of the perceived cell boundary, so there is no contamination of the net strain energy from noise in the exterior displacement data.

There are several advantages to FTTC common to both the constrained and unconstrained implementations. These include the following: (+) The moment matrix is an especially simple formulation in Fourier space. No additional calculations are needed.

(+) The use of the FFT implies that the entire problem of traction measurement is no longer computationally intensive; large image pairs can be analyzed in seconds or minutes.

### Effect of Noise on Traction Recovery in the Example Smooth Muscle Cell

Here we describe the effect of noise on our actual computations of moments and strain energy associated with the smooth muscle cell shown in Figs. 1-5. A comparison of the root mean square tractions exterior to the cell boundary with that associated with pure noise, as described in results, gives a rough estimate of the noise in the displacement field and is a conservative estimate insofar as there appear to be patches of non-zero correlated tractions in Fig. 4, possibly secondary to other cells exterior to the field of view. In our case, this results in an estimate of ∼0.05 μm root mean square noise level in the displacement field. From the pure noise simulations quoted in results, this amounts to an uncertainty in the net contractile moments of 0.015 and 0.010 pNm in the unconstrained and constrained cases, respectively. These numbers should be compared with those in Tables 1 and 2, where the net contractile moment of the cell is ∼3 pNm in magnitude and shows that noise is negligible in its contribution to these estimates. By contrast, as remarked above, the strain energy is quadratic in the noise level, so displacement noise of 0.05 μm corresponds to roughly 0.025 and 0.0025 pJ in the unconstrained and constrained cases, respectively. The difference in the estimated strain energies for the real cell in the unconstrained and constrained methods, from Table 1and Table 2, is roughly 0.1 pJ, and we conclude that noise may account for some quarter of this difference. Inspection of Fig. 4, however, shows that this is also to be expected, because the patches of correlated tractions will certainly contribute to these differing estimates.

### Remarks on the DW Method

The DW method specifies both tractions (exactly) and displacements (approximately) exterior to the perceived cell boundary. This particular specification of the problem requires special techniques because these are approximate Cauchy conditions, for which the elliptic Navier equations of elasticity in general have no solution. This is an ill-conditioned problem that necessitates smoothing or regularization to obtain stable solutions. The DW method utilizes the regularization method introduced by Tikhonov in the 1940s (summarized in Refs. 8 and 9) with the choice of smoothing functional and level of smoothing introduced by Phillips (6) and Twomey (10). In brief, the residuals of the displacement field plus a certain amount of the L^{2}norm of the gradient of the traction field (the regularizing functional) are minimized, and the displacement residuals are examined. The level of smoothing is then varied until there is an appropriate level of variation in the predicted displacements, given a priori information about the noise in the displacement field (6). For problems that are highly ill conditioned, and when there are unambiguous a priori constraints, this kind of approach is often useful and appropriate (2). By contrast, the displacement kernel in the Boussinesq solution decays like *r*
^{−1}, which is sufficiently rapid that in FTTC we do not find unacceptably large oscillations in the recovered traction field that would necessitate a Tikhonov-type approach.

### Recommendations and Conclusion

On the basis of the results presented, the best strategy for measuring traction field using FTTC may be summarized as follows. An initial examination of the traction field recovered with the use of unconstrained FTTC will reveal the extent of significant tractions exterior to the perceived cell boundary. If these can be determined to result from contractile cells exterior to the field of view, the traction maps obtained with constrained FTTC may be more accurate. On the other hand, such tractions might arise from real structural elements that are not seen in phase-contrast microscopy and that are preserved in unconstrained FTTC. Whichever method is used, the traction moments are easily computed in Fourier space, whereas the strain energy is best computed with constrained FTTC, integrating over real space. If desired, the noise level in the displacement field may be estimated by examination of its Fourier spectrum at high frequencies, where the white noise is manifest as a constant level of intensity.

In conclusion, Fourier transform traction cytometry is a new solution to the problem of mapping of the traction field between a cell and its substrate, given the displacement field between two micrograph images. This method has the advantages of being exact, computationally efficient, and not subject to certain artifacts that can lead to misleading conclusions. This approach also yields simple measures of the net contractile moment of the cell, the strain energy imparted to the substratum, the orientation of the principal tractions, and a quantitative index of cell polarity. FTTC may represent a new and important tool for studying the mechanical properties and function of adherent cells.

## Acknowledgments

We thank S. M. Mijailovich for performing finite element simulations of test cases early in the course of this work. We especially thank E. J. Millet and N. Wang for helpful discussions. This work was stimulated in large measure by collaboration with N. Wang on his studies of cell mechanics. We thank J. Chen for technical assistance in cell and gel preparations. Cells were kindly supplied by R. Panettieri.

## Footnotes

This work was supported by National Heart, Lung, and Blood Institute Grant HL-P01–33009.

Address for reprint requests and other correspondence: J. P. Butler, Physiology Program, Harvard School of Public Health, 665 Huntington Ave., Boston, MA 02115 (E-mail:jbutler{at}hsph.harvard.edu).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.10.1152/ajpcell.00270.2001

- Copyright © 2002 the American Physiological Society