A Novel Cell Traction Force Microscopy to Study Multi-Cellular System

Traction forces exerted by adherent cells on their microenvironment can mediate many critical cellular functions. Accurate quantification of these forces is essential for mechanistic understanding of mechanotransduction. However, most existing methods of quantifying cellular forces are limited to single cells in isolation, whereas most physiological processes are inherently multi-cellular in nature where cell-cell and cell-microenvironment interactions determine the emergent properties of cell clusters. In the present study, a robust finite-element-method-based cell traction force microscopy technique is developed to estimate the traction forces produced by multiple isolated cells as well as cell clusters on soft substrates. The method accounts for the finite thickness of the substrate. Hence, cell cluster size can be larger than substrate thickness. The method allows computing the traction field from the substrate displacements within the cells' and clusters' boundaries. The displacement data outside these boundaries are not necessary. The utility of the method is demonstrated by computing the traction generated by multiple monkey kidney fibroblasts (MKF) and human colon cancerous (HCT-8) cells in close proximity, as well as by large clusters. It is found that cells act as individual contractile groups within clusters for generating traction. There may be multiple of such groups in the cluster, or the entire cluster may behave a single group. Individual cells do not form dipoles, but serve as a conduit of force (transmission lines) over long distances in the cluster. The cell-cell force can be either tensile or compressive depending on the cell-microenvironment interactions.

Efforts at visualizing cellular traction forces may be traced back to 1980s when Harris et al. used thin polymeric silicone substrates for cell culture, and observed the wrinkling phenomena caused by the traction of migrating cells [45]. However, quantitative estimation of the traction from the wrinkling of silicone substrates is challenging due to the inherent non-linearity of the problem. From 1995 on, Lee, Jacobsen and Dembo et al., as well as other groups, developed several traction force microscopy techniques (TFM) to quantify the cellular traction produced by migrating or stationary cells on soft substrates [46,47,48,49,50,51,52,53,54]. TFM computes the cell traction forces from the deformation of a soft substrate with known elastic properties, such as polyacrylamide (PA) gel, on which cells are cultured. The deformation is measured from the displacements of micro-fluorescent markers embedded in the substrate. The motion is measured from two images. First image is taken with the cells adhered to the substrate. Here, the cells have generated traction force on the substrate, and the image gives the deformed configuration of the soft substrate. Then cells are removed from the substrate through trypsinzation, and a second image is taken. Subsequently, the substrate is relieved of cell traction, and the image shows the un-deformed configuration of the substrate. A comparison of the two images gives the displacement field of the substrate's top surface due to cell tractions. Digital image correlation method (DICM) is used to quantify the displacement field. The traction field is estimated from the displacement field. Several methods have been proposed for force estimation ranging from analytical methods, i.e. the Boussinesq formulation (either using Bayesian likelihood regularization method [51,55] or Fourier transformed approach [49]), to computational methods like finite element analysis (FEA) [56]. The Boussinesq formulation approach, which assumes the substrate as a semi-infinite elastic half space [57], was first adopted by Dembo and Wang, et al., to compute the traction forces from the displacement fields followed by regularization [51,55,58,59]. Since the Boussinesq formulation involves solving an inverse problem, the solution demands computational regularization schemes to predict the approximate traction solutions. Importantly, Butler, Trepat and Fredberg, et al. [49,60,61,62,63] made significant progress in mitigating some pitfalls of the regularization scheme by solving the Boussinesq equation using Fourier transform. Later Schwarz et al. introduced a new method to compute traction forces only at the focal adhesion site of the cell by assuming that the cell force transfer occurs only through these sites [50]. Some novel platforms, such as the photobleachingactivated monolayer with adhesive micro-patterns developed by Scrimgeour et al. [64] and the elastic substrates with micro-contact printing demonstrated by Stricker et al. [65], were also used to characterize the cell traction force. Furthermore, a FEA-based technique was also developed by Yang et al. to greatly improve the accuracy of traction force calculations [56]. The FEA method no longer depends on the Boussinesq formulation and thus is not limited by the semi-infinite elastic half space assumption [66,67]. Recently additional contribution has been made in traction force computation in three dimensions [19,68,69,70,71,72,73]. 3D TFM techniques compute the 3D traction force fields from the cell induced 3D displacement and strain fields obtained using laser scanning confocal microscopy (LSCM) and digital volume correlation (DVC). However, it is challenging to obtain the Zdimension displacement field and the technique can only be applied to single cell cases, rather than multiple cells or cell clusters.
The above studies focused on traction force computation for single cells far from their neighbors, i.e. cells that do not interact mechanically with each other. However, live cells do interact with their neighbors chemo-mechanically and form cell clusters [7,29,37,74,75,76]. In this paper we present a novel finiteelement-based TFM technique to compute the traction fields generated by multiple cells and clusters. We first present a theoretical proof showing that the 3D traction field computed from prescribed displacement field of the substrate is unique. We verify the uniqueness by considering a 2-cell case. We test the accuracy of the computational technique by applying a known force on PA gel substrate using a micro-needle, and by comparing the experimental force with the computed one. Finally, we compute the traction fields generated by multiple cancerous and fibroblast cell clusters, and reveal that cells might be under compression in such 2D clusters. We believe that the present technique may enable better examination and understanding of a variety of biological phenomena involving homotypic and heterotypic cells and cell cluster interactions [77,78,79].

Uniqueness of traction field computed from displacement field in 3D linear elastic solids
Consider a 3D linear elastic solid with volume V in static equilibrium. Its boundary, S, consists of S u and S s (S = S u +S s ) where displacements u B i and traction t B i are prescribed respectively.
Proposition: Given displacement field u B i at S u and traction t B i at S s , the corresponding traction t i at S u is unique. (Note: indices i, j = 1,2,3 correspond to x,y,z Cartesian coordinates respectively; all equations follow standard tensor notation and summation convention). Supporting material Text S1 presents the proof of the proposition.

Simple 1D examples of uniqueness
Displacement boundary condition. To gain an intuitive insight on the uniqueness of traction solution, we present a simple 1D model. The stiffness and compliance matrices of a uniform elastic bar have been derived (see Supplementary Materials Text S2). In Fig. 1a, a uniform bar is subjeted to three concentrated force F 1 , F 2 , F 3 at nodes 1, 2, 3 with corresponding displacements u 1 , u 2 , u 3 and linear stiffness k 1 , k 2 , k 3 respectively. For simplicity, let k 1 = 1, k 2 = a, k 3 = b, then the displacement-force relationship is given by: Thus, if displacements (u 1 , u 2 , u 3 ) are given, nodal forces (F 1 , F 2 , F 3 ) can be obtained from Eqn (1). For a given displacement (u 1 , u 2 , u 3 ), there is only one possible value of (F 1 , F 2 , F 3 ), i.e., solution for the nodal forces is unique, since Young modulus is positive, compliance and stiffness matrices are positive definite and hence are non-singular and invertible. In other words, there is always a unique relationship between displacement and forces on the nodes.

Author Summary
Adherent cells sense, transduce and respond to their microenvironment by generating traction forces on their surroundings. To accurately understand these mechanotransduction processes, it is critical to have a robust and reliable method for traction force visualization and quantification. However, most cell traction force microscopy methods are limited to only single cell traction force analysis. Considering that most physiological processes are essentially collective multi-cellular events, there is a need for traction force microscopy methods capable of analyzing traction forces resulting from multiple cells. We have developed a novel and robust multi-cellular traction force microscopy method for computing cell traction on soft substrates, and applied it to compute traction field generated by both multiple cells and cell clusters. We verified the accuracy, robustness, and efficiency of the method by theoretical, numerical and experimental approaches. Our method provides a powerful toolset to pursue the mechanistic understanding of collective biological activities, such as cancer metastasis and neuromuscular interactions.
Mixed boundary value problem. In Eq (1), given all displacement data at the nodes, the force vectors can be directly calculated, and vice versa. However, there are cases where a combination of displacement and forces on boundaries are given, called mixed boundary value problems (MBVP). For example, suppose a cell is adhered at nodes 1 and 2 and creates contractile forces F 1 and F 2 . The corresponding displacements u 1 and u 2 are measured but there is no anchorage at node 3 (i.e zero traction) and hence F 3 = 0 (Fig. 1b). Given: u 1 = 0, u 2 = 21, F 3 = 0, let us find contractile forces F 1 , F 2 and displacement u 3 . Our unknowns are a combination of forces and displacements. Applying the given boundary conditions into Eqn (1) and solving for the unknowns, ,u 3~{ 1: Note that F 1 ?0 although u 1 = 0 and similarly u 3 ?0 while F 3 = 0. Therefore, zero displacement does not necessarily result a zero force (traction) at a node and vice versa. The displacement at zerotraction node 3 is due to the displacements at other nodes. This example illustrates the counter intuitive possibility of non-zero displacements at points on the body where traction is zero. Also note that, F 1 zF 2~0 which confirms the traction field under the cell is self-equilibrated.

Finite element approach for solving cell traction on 2D substrates
We illustrate our computational scheme as follows. Consider two separate cells on a soft elastic substrate. The substrate is adhered to a rigid surface (such as glass) at the bottom. The lateral boundary of the substrate is far from the cells. In the finite element scheme, the substrate is modeled as a rectangular pyramidal solid body. It is discretized as a collection of small cubes with common nodes. We need to prescribe three boundary conditions, namely any combination of forces (F x , F y , F z ) and displacements (u x , u y , u z ), at each of the surface nodes. For example, (F x, u y , u z ) can be a boundary condition at a surface node. To ensure that the body is at rest (no rigid body translation or rotation), at least two of the nodes are prescribed with u x , = u y , = u z = 0. Given the boundary conditions, finite element scheme calculates the deformation of the solid body such that the total energy is minimized. Thus the displacements at each node within the body, and at the surface nodes where forces are prescribed are evaluated. This leads to the evaluation of strains and stresses using the elastic properties of the solid (Young's modulus and Poisson's ratio for the isotropic gel). Surface traction is calculated from the stress near the surface and normal vector to the surface (t i~sij n j ), as shown in Supplementary Materials S1. Surface nodal forces are calculated from an area integral of traction at the vicinity of the node. Thus, the analysis provides the forces at nodes where displacement is prescribed, and displacements where forces are prescribed. If (F x , u y , u z ) is prescribed at a surface node for example, one gets (u x , F y , F z ) at that node. Even though the solution is unique in principle, errors are introduced if the discretization is coarse. With finer discretization, the solution converges to the correct one. This convergence test is often employed to gage the accuracy of the solution.
In our problem with two cells, we prescribe zero displacement boundary conditions at the bottom surface and at the four vertical sides of the body (Fig. 2). Thus all the nodes on the bottom and the vertical sides are fixed. For simplicity of illustration, consider that there are a few nodes on the top free surface outside the cell boundary, and a few nodes within (Fig. 2). Our objective is to calculate the traction on these nodes. We can experimentally Figure 1. Modeling of cell contraction on 1D elastic substrate with mixed boundary conditions. (a) An elastic bar is discretized with 3 nodes with concentrated forces applying on each node along with their respective displacements. Note that this general loading is used for deriving stiffness matrix which uniquely relates nodal forces to the nodal displacement subject to any boundary condition. (b) A cell applies contractile forces on nodes 1 and 2 (i.e with known, measured displacements) while node 3 is free. This set of inputs constitutes a Mixed Boundary Condition, in that a combination of nodal displacements and forces are given (u 1~0 ,u 2~{ 1,F 3~0 ) and their respective unknowns (F 1 ,F 2 ,u 3 ) are computed by the model. doi:10.1371/journal.pcbi.1003631.g001 measure displacements (u x , u y , u z ) at all the nodes on the surface. They are generated by cell forces, although we do not know the precise locations of these forces. We also know that the surface outside the cells has no traction, and that each cell or cell cluster produces a traction field that is self-equilibrated, i.e., the sum of forces applied by the cell or the cell cluster on the substrate is zero. Cell traction can be evaluated by prescribing either of the two boundary conditions: i) Whole-field displacement boundary conditions (BC1): (u x , u y , u z ) are prescribed at all the nodes on the surface, and their traction is analyzed by FEM. ii) Mixed boundary condition (BC2): zero traction is prescribed at all the nodes outside the cells (F x , = F y , = F z = 0), and displacements (u x , u y , u z ) are prescribed for nodes within the cell boundaries (Fig. 2).
Remarks. (1) The mixed boundary scheme applies exact boundary condition (zero force) at nodes outside the cells. Hence none of the displacements (u x , u y , u z ) need to be prescribed at these nodes. Thus, it is not necessary to measure the displacements of the beads outside the cells. Due to the exact boundary conditions outside the cells, the traction solution is expected to be more accurate. However, errors will be introduced if the cell boundary is incorrectly defined and there are nodes that fall outside the cell boundary where cells apply traction. In cases where the cell boundaries cannot be identified due to imaging conditions (Fig. 3), displacements should be prescribed for regions nearby the cells.
(2) Displacement u z and Poisson's ratio: It is shown in the supplementary material (Supplementary materials text S3, Fig.  S1b and c), that if the Poisson's ratio of the gel approaches 0.5, then the in-plane displacements, (u y , u z ), on the surface of the gel are independent of the out-of-plane component of traction (F z ). That is, (u x , u y ) are determined by (F x , F y ) on the surface. Similarly, u z is determined by F z on the surface only. Thus, in order to evaluate the in-plane traction only, one needs to measure and prescribe in-plane displacements only at the surface nodes, and prescribe arbitrary boundary condition in z direction (i.e F z = 0 or u z = 0) at all surface nodes, when Poisson's ratio is close to 0.5. We experimentally measured the Poisson's ratio of our gel as 0.4760.02 (Fig. S3b, n = 5). In order to estimate the in-plane traction only, we have prescribed F z = 0 for all nodes within the cells in the rest of the paper. This results in an error of less than 2% in the calculation of in-plane forces F x and F y (Supplementary materials text S3 and Fig. S3b). If F z is desired, one needs to measure and prescribe (u x , u y , u z ) at the surface nodes. Also, if Poisson's ratio is much less than 0.5 (e.g., 0.35), (u x , u y , u z ) must be prescribed at the nodes within the cells even when only in-plane traction is desired.

Validation of uniqueness of solution in finite-element models
In this section, we demonstrate computationally that the traction solution from finite element simulation is unique as long as the full 3D boundary conditions are prescribed. We define two circular boundaries representing two cells with half-cell distance apart on a soft gel surface. The diameter of each boundary is chosen as 20 mm, close to real cell size. A three-dimensional finiteelement (FEM) block model is generated (ANSYS 12.0 Workbench Package) to represent the PA gel substrate [79][80][81][82][83][84][85][86][87][88][89][90][91][92][93][94][95][96][97][98]. The gel is presumed linear elastic, isotropic, and homogeneous in their mechanical properties for a wide range of deformations [78,99]. The Elastic modulus, E, of the gel is 1KPa (our experimental value is 1.0560.17 kPa, measured by AFM indentation (n = 15; Fig.   S3a), [99][100][101]]). The model height is 70 mm, same as the thickness of PA gel used in experiments. We first apply an in-plane force field (Fig. 4a) within each boundary, and compute the corresponding displacement field, u x , u y , u z (Fig. 4b). Second, we use the computed u x , u y and u z within the cell boundary on the surface (Fig. 4c), and zero-traction conditions outside the boundaries to calculate the traction within the cells (Fig. 5d). A comparison between the prescribed and the calculated forces from the two steps shows close quantitative agreement (within 1%) (Fig. 4e-f). Note that individual cells or cell clusters generate selfequilibrated traction on the substrate. Hence, we use a measure of accuracy of the traction solution by defining the error ratio, where F xi and F yi , are the nodal force components within the individual cells, and i = 1, n, the number of nodes within the cell or cell cluster boundary. For exact solution, e = 0.

Demonstration of 2-cell experiments using mixedboundary condition method
In this section, we demonstrate the applicability of the method by evaluating the traction induced by two neighboring cells. Here, two monkey kidney fibroblasts were plated on PA gel (1 kPa) with Poisson's ratio of 0.47 (Fig. 5a). Two different regions (two sets of S u and S s ) were selected to prescribe the displacement boundary conditions: (1) displacement field underneath the two cells were prescribed in the model (the white parts in Fig. 5b), whereas the traction-free condition was applied outside the cells (the black part in Fig. 5b); (2) the displacement field within a region enclosing both cells was prescribed (the white part in Fig. 5c), whereas the traction-free condition was applied outside this region (the black part in Fig. 5c). The out-of-plane force, Fz, was prescribed as zero within the cellular regions in (1) and (2). The traction fields were calculated for both cases (Fig. 5d,e,g,h), and compared ( Fig. 5f and i). The RMS~ffi

Demonstration of whole-field displacement boundary conditions method and comparison
In this section, we compare our mixed-boundary condition method with traditional whole-field displacement boundary condition method, which requires iterative calculation and has been successfully used by Fredberg, et al [49,102]. Briefly, the iteration calculation proceeded as follows: (a) we assigned the complete 2D DICM (digital image correlation method) displacement data (u x ,u y ) for all nodes of the top surface of the gel (both intracellular and extracellular regions ; Fig 6a-b). We prescribe F z = 0 within the cluster for both the mixed boundary condition and iterative methods. (b) The traction field was solved using FEM. Then all the forces in the extracellular region were replaced by F x = F y = F z = 0 to satisfy the traction-free condition, while the forces in the intracellular region were retained intact. (c) The new traction field was used to generate a new displacement field using FEM. Thus a new displacement field was computed within the intracellular region. (d) The computed intracellular displacement field was replaced with the DICM displacement field (u x and u y ), while the computed extracellular u x , u y , and u z from previous step were retained intact. (e) The steps (b), (c), (d) were repeated until  the solution converged, i.e., the difference between the root mean square (RMS) of surface nodal forces in two consecutive cycles became less than 5% (Fig. 6c-e).
Our computational results showed that the solutions from mixed-boundary and iterative methods converge (Fig. 6c-e). We found, the difference between the root mean square (RMS) value of traction from the two methods was 1.6610 21 kPa (Fig. 6f), less than 3.8% of the maximum computed cell traction. The difference between the RMS of the nodal forces was 0.2 nN, which is 0.25% of the maximum nodal force at cell cluster -substrate interface (Fig. 6g). The distribution of traction |t| and forces DF D at nodes (Fig. 6h-6j) shows good agreement between the two methods. We used e (Eqn 2) as a measure of accuracy of the traction solution.

Mesh size effect
In FEM, convergence test is required to determine the optimal mesh size needed to obtain the accurate solution. Three mesh sizes, 3.23 mm, 4.84 mm, and 6.45 mm were tested, as shown in Fig. 7a-c, and used to calculate the traction field of the same cell cluster by mixed-boundary condition method. The distribution of nodal traction and forces showed minor difference between the three mesh sizes (Fig. 7a-c and 7e). The values of e were 4.74%, 6.69%, and 6.12% for mesh size of 3.23 mm, 4.84 mm, and 6.45 mm respectively (Fig. 7d). Therefore, in the following computations, mesh size of 4.84 mm was used for analysis. The upper limit of mesh size is dependent on the specific cell size and the gradient of the traction field produced by the cell. A starting point on mesh size can be ,20% of cells size.

Traction calculation for multiple cell clusters
A key attribute of the present method is the computation of traction fields generated by multiple cell clusters interacting with each other. Each cluster may consist of multiple cells, and the cluster size might be similar to or larger than the thickness of the soft substrate. Hence the effect of the glass-gel interface needs to be considered, and the gel may not be treated as half space. In the following, we study several cell clusters (Figs. 8-10) and outline the main biological findings. The mixed-boundary condition method was used to compute the traction fields.
Cells in clusters exert cell-cell forces. Fig. 8a shows a cluster of colon cancer cells (HCT-8) on 2 kPa substrate. These cancer cells are epithelial in nature, and have low metastatic potential. Figs. 8b and c show the traction generated by the cluster and the corresponding nodal forces on the substrate. Here the grid size in the analysis is about 5 mm. We note that the cells within the cluster do not form individual force dipoles. The entire cluster behaves as a single contractile unit, and generates forces along the periphery. The concave peripheries generate larger inward forces. These forces are balanced by far away opposing forces. Thus the cells behave as generators and transmitters of forces from one edge to the other of the cluster. The mechanics of this force transmission can be visualized from a free body diagram (Fig. 8d). Here the cluster exerts contractile forces on the substrate through the adhesion sites of the outer cells. The cells inside the cluster contribute and transmit the force possibly through cadherin junctions and force-bearing cytoskeleton. Thus the cell-substrate traction is transduced to cell-cell contractile forces. As if, the peripheral cells pull the interior cells outward. This is consistent with the observation that advancing edge of a cell cluster pulls the interior cells [49,102,103,104,105,106,107]. Previously we reported that HCT-8 clusters are sensitive to their mechanical microenvironment and display a metastasis-like phenotype transition (MLP) when cultured on appropriately soft substrates [7,8,28,29,30,31]. This MLP transition always initiates from the periphery of clusters. It remains to be seen whether the difference in forcing on the peripheral cells compared to those in the interior contributes to the MLP transition.
Cell clusters may generate traction interior to the periphery. Fig. 9a shows two clusters of pre-MLP HCT-8 cells close to each other, cultured on 2 kPa substrate for 5 days. The traction and the force maps (Fig. 9b, c) show that there are two regions well within the boundaries where the traction is high, unlike the cluster of Fig. 8 where the traction was mostly peripheral. Here the interior forces are balanced locally, i.e., these forces form local dipoles leaving the rest of the cluster nearly traction free. This may explain the spherical morphology of the clusters, which minimizes their surface areas. The cells of the clusters might be under compression due to line tension of the peripheral cells.
Traction domains of cell clusters are dynamic. The cell clusters of Fig. 9a merge on the 6th day of culture, as shown in Fig. 9d. Soon after merging, the traction and force pattern changes  (c-e) The traction field calculated by mixed-boundary method, and full-field displacement boundary method (with iterative calculation 1 time and 2 times, respectively), were shown respectively. The difference of RMS of the traction between mixed-boundary method and full-field displacement boundary method with 1 time iteration was 1.6610 21 kPa, less than 3.8% of the maximum computed traction at cell cluster and substrate interface. The difference of RMS of their nodal force was 0.2 nN, which was 0.25% of the maximum nodal force at cell cluster and substrate interface. Dashed lines in orange outline the cells boundaries. (f-g) Histograms of nodal traction and force obtained by the two methods demonstrated good agreement between each other. (h) Sum of net forces and absolute forces calculated by the above three conditions. The force equilibrium was best satisfied in mixed boundary condition method, which is 6.69% of total force. (i) Sum of surface nodal force distribution calculated by above three conditions. The RMS results of nodal force calculated by mixed BC method and 1-time iteration method agreed within 4.96%. (j) Sum of surface nodal traction distribution calculated by the above three conditions. The RMS results of nodal force calculated by mixed BC method and 1-time iteration method agreed within 9.27%. doi:10.1371/journal.pcbi.1003631.g006 dramatically (Fig. 9e-f). The net force increases by about 20 folds, although the direction of the net force does not change. Many more cells in both the clusters now participate in generating traction. The new traction regions are also mostly interior to the periphery, and the merged cluster takes a smooth elliptical shape. The merger between two cell clusters mimics that between two droplets, as they both tend to minimize the surface energy. It is known that cells may interact with each other through substrate strain fields [4,11,108,109,110]. In case of the two neighboring clusters, the displacement fields were localized well within the clusters. It is thus unlikely that their merger was induced by strain fields.
Evidence of cell-cell compression. It is generally understood that cells generate contractile forces produced by actomyosin machinery [11,76,111,112]. However, in a 2D cluster, cells may be subjected to compression as shown in Fig. 10. Here, monkey kidney fibroblasts (MKF) form several large and small clusters on 1 kPa substrate. Each cluster is sufficiently far away from the others so that there is no mechanical coupling between them. The displacement field between them is negligible (Fig. 10b). The traction within each cluster is shown in Fig. 10c. Unlike the previous two examples, here many more cells in the clusters participate in traction generation. Fig. 10d shows the nodal forces of the larger cluster. Here, several regions generate dipole type forces within the cluster. However, there are interior boundaries where opposing forces appear on the substrate, i.e., cells ''push'' against each other. This can be explained by the schematic of Fig. 10e where neighboring cells have adhesion sites with the substrate. Due to the low stiffness of the substrate, the cells have less likelihood of spreading or wetting the substrate, though they may adhere to the substrate due to the fibronectin functionalization. Now, if growth occurs in any of the cells next to a neighbor, it would push out the neighbor generating an outward force on the substrate. This results in a compressive stress between the cells. There are three regions of such cell-cell compression in the cluster of Fig. 10d (enclosed by dashed lines).

Discussion
The majority of fundamental physiological processes in tissue development, health, and disease are coordinated by the collective activities of multiple cells [60,62,76,102], rather than single cells [10,103]. To understand how mechanical traction applied by neighboring cell cluster groups could specify or mediate the tissue functionalities [7,8,11,75,104,105,106], robust cellular traction evaluation method is indispensable. In the present study, we developed a finite element element-based traction force microscopy (TFM) to accurately compute and visualize the traction maps resulting from multiple cell clusters. The uniqueness, convergence, and correctness of traction solutions are substantiated. We showed that as the gel Poisson's ratio .0.4, the in-plane traction can be obtained with minimal error from the in-plane displacement field alone. For Poisson's ratio ,0.4, both in and out of plane traction depend on both in and out of plane displacement boundary conditions, and it is essential to measure these displacements to compute any of the traction components. The method presented is applicable to substrates with any value of the Poisson's ratio. It calculates the full 3D traction field given the 3D displacement boundary condition within cells or cell clusters. Moreover, unlike the classical TFM methods that are based on Boussinesq solutions [39,40,48,49], the FEM takes into account the effect of substrate thickness and nearby environment. It is now known that cells can sense the substrate depth within the cellular length scales by showing distinct morphological variation on the gel substrate with same Elastic modulus but with varying thickness [22,107]. We applied the method to compute the traction generated by multiple cell clusters. Some of the clusters were more than 100 mm in size consisting of many cells, while others were in close proximity to each other. The computational scheme presented here is ideal for studying such clusters, since the domain of traction field is much larger than the thickness of the gel, and one needs to account for the finite thickness of the substrate. A few interesting biological insights emerge from these analyses. First, the cluster may behave as a single contractile unit where the peripheral cells serve as anchorage sites. Force is transmitted between distant peripheries by the cells inside the cluster. Thus the cells are subjected to tensile intercellular forces, as if the peripheral cells are pulling the interior cells outward. It needs to be seen whether there are specific cells within the cluster that generate the force, or all the cells behave as contractile actuators. In any case, the cells probably use cell-cell junctions and cytoskeleton to transmit the force through the cluster.
We also found instances where traction is limited to small regions well within the clusters. These regions can have locally balanced traction (forming dipoles), leaving the rest of the clusters nearly traction free and weakly adhered to the substrate. These clusters are spherical in morphology, as expected. The traction free regions tend to minimize the surface area by being circular, just as a free-standing cell cluster takes a spherical shape. It is plausible that the cells within the circular clusters are under compression due to the surface tension of the peripheral cells. In any case, the interior traction maps can be highly dynamic. When cell clusters merge, the traction map can change their orientations, and the net force can increase by an order of magnitude over short times.
It is known that cells generate contractile forces. Thus, it is expected that the cells in a 2D cluster will be under intercellular tension. We found evidence to the contrary. If the cells are on soft substrates where they do not spread much, but they adhere to the substrate, then some of the cells in the cluster may be subjected to compression. We found regions within such clusters where the neighboring cells apply repulsive forces on the substrate, i.e., the cells are pushing against each other while being adhered to the substrate. One possible explanation might be that the neighboring cells are growing, but their adhesion sites are stationary.
In conclusion, we developed a robust FEM-based cell traction force microscopy technique to estimate the traction forces produced by multiple cells and clusters. The utility of the technique is exemplified by computing the traction force fields generated by multiple monkey kidney fibroblast (MKF) and pre-MLP human colon cell (HCT-8) clusters in close proximity. The developed technique is user-friendly and computationally inexpensive. Our FEM-based traction force microscopy provides a powerful tool to probe multi-cell questions involving assembly/ disassembly dynamics of cell ensembles, tissue network formation, and wound healing. Future work is needed to determine the subcellular processes involved in mechano-sensing and regulation, and their respective timescales.

PA gel substrate preparation and ECM conjugation
Polyacrylamide (PA) gel substrates with 1 kPa stiffness used in present study were made by mixing 12.83% (v/v) of acrylamide (Sigma-Aldrich, Inc.), 1.54% (v/v) of N, N-methylene-bisacrylamide (Sigma-Aldrich, Inc.), 2% (v/v) of 1 mm diameter fluorescent micro-beads (Invitrogen, Inc.) and 10 mM Hepes (Gibco., Inc.) [7,11]. Solution was vortexed thoroughly for 5 min to obtain uniform distribution of beads. TEMED and ammonium persulfate (Fisher Scientific, Inc.) were used to initiate PA gel crosslinking. Chemical modification of glass slides and preparation of PA gels were carried out following the procedures described previously [50,80,81,82,83,84,85]. Briefly, a circular glass coverslip (Fisher Scientific, Inc.) of 1.2 cm in diameter was placed on an acrylamide solution drop on activated coverslip and placed on the bottom of a petri dish. Capillarity spreads the drop and fills the space between the circular coverslip and the activated coverslip. The gel was cured at room temperature and reached to the stabilized thickness of 70 mm [82,85,86]. The circular glass coverslip was peeled off from the gel that remained on the activated cover slip. The surfaces of the air dried PA gels were activated by incubating in 97% hydrazine hydrate (Acros Organics.) for 12 h followed by a complete rinsing with DI water and 30 minutes incubation along with gentle shaking in 5% acetic acid (Avantor Performance Material, Inc.) [7,8,11,13,81]. Solution of human fibronectin (25 mg/ml, BD Biosciences) was prepared by dissolving in phosphate buffer saline (PBS) and the carbohydrate groups of fibronectin were oxidized by sodium periodate (Sigma-Aldrich, Inc.). To minimize the displacement noise and rigid body motion during imaging, the glass slides was firmly adhered to the bottom of 30 mm petri dish using adhesive glue (Henker Consumer Adhesive, Inc.). Full experiment procedures and sample characterization are provided in Supporting Materials Text S4-S9 and Figures S1-S5. Figure S1 Confocal microscopy images of monkey kidney fibroblasts (MKF) cells on gel substrate with immunofluorescent stained F-actin cytoskeleton (green) and focal adhesion protein, Vinculin (red). The x-y plane shows the horizontal view of spread MKF cells. The z-y and z-x cross-sectional planes, A and B, show the vertical structures of spread MKF. They display that the height-to-length ratio of spread MKF cells is in the range of 1/40,1/50. Vinculin staining (red) indicates the basal surface of MKF cells. This low height-tolength ratio implies that the cells exert their traction forces mostly along the x-y plane through their contractile filaments. The cartoon of spread cells on top right of (a) shows that the angle Q between contractile cytoskeleton (green) and substrate is within the range of 1,10 o . (b) To estimate the error due to out-of-plane forces on the evaluation of in-plane traction, a general 3D forcedisplacement model for the cell is developed. In the model, the cell applies both in-plane and out-of-plane forces on the substrate, Q and P, with corresponding deformation u and w. (c) The error index plotted for all three cases, P = 0, w = 0, and general loading P = kQ v.s Poisson ratio n. For P = 0, there is no error in planar force calculation for all position x and displacement u. For other cases, however, there are deviations due to the presence of out-ofplane force, P, at different boundary conditions. It is evident from the plot that for small values of Poisson's ratio, the z-component of deformation w will influence the in-plane force Q and thus create varying results depending on loading modes and the value of Poisson's ratio. Therefore, excluding out-of-plane deformation w will introduce error in calculating the in-plane force, Q. However, as Poisson's ratio approaches K, most of the discrepancies in planar force calculations becomes negligible, and all set of curves converge to a unified value corresponding to P = 0, regardless of value and direction of the out-of-plane force P and deformation w. (TIF) Figure S2 A representative elastic body subject to the most general form of mixed boundary condition. Displacement field and traction field are given on separate surfaces S u and S s respectively. The body has total surface S = S u +S s and total volume V. The general state of stress tensor s ij and respective displacement vector u i are shown at an arbitrary point P within the body. Cauchy traction vector t i applies on an arbitrary, infinitesimal surface denoted by the unit normal vector n j . (TIF)  Figure S5 (a) A Tungsten probe with known stiffness of 10.74 nN/mm (calibrated with weight) was vertically held by a high-resolution x-y-z piezo-stage to apply horizontal force on the flexible hydrogel surface. (b) The deflections of probe tip with respect to reference base, as well as the resultant displacement fields of beads on gel's top surface, were recorded. The displacement fields were assigned to FEM model to compute the resulting force. The double-headed arrows indicated the gap between micro-needle and reference base. Multiplying this gap with spring constant of the micro-needle provided the force applied on the substrate. (c) The sum of nodal reaction forces on PA gel was calculated using present traction force microscopy and compared with the needle force. The relative error in force estimation is within 6.5%. (TIF)

Supporting Information
Text S1 Proof of uniqueness of traction field computed from displacement field in 3D linear elastic solids.