2.2.1 Deformation fields and TFM.

Deformation fields are calculated from the images of the substrate in a tensed and relaxed state using the cross correlation-based Particle Image Velocimetry (PIV) algorithm implemented in the openPIV Python package [9]. PIV is performed by selecting for example a 50x50 pixel tile around a given pixel from the tensed image and shifting the tile by pixel increments in all directions across the corresponding tile in the relaxed image. This yields a correlation matrix of in this case 99x99 pixels. The deformation vector is then obtained by calculating the vector between the position of the highest correlation and the center of the matrix. The initial deformation vector is further refined to sub-pixel accurate values by fitting a 2D Gauss curve to the directly neighbouring correlation values. To reduce noise, deformation vectors with a signal-to-noise ratio smaller than 1.03 are exclude and replaced by the local mean of the surrounding deformations at distances < = 2 pixel. The signal-to-noise ratio of each deformation vector is defined as the ratio of the correlation of the highest peak and the correlation of the second-highest peak outside of a neighborhood of 2 pixels around the highest peak.

A common issue when calculating the deformation field is a global drift of the images, which needs to be corrected prior to PIV. This drift can reach several μm and is caused by the combined effect of positioning inaccuracies of the motorized x, y-stage, mechanical handling or shaking during the addition of trypsin-EDTA, and slow temporal or temperature-induced mechanical drift. pyTFM offers a global drift correction that works in three steps: First, the drift is estimated with sub-pixel accuracy by cross-correlating the entire first image with the entire second image. This is done with the “phase_cross_correlation” function from the scikit-image [25] python package. Next, the first image is shifted by the drift, and finally both images are cropped to the overlapping field of view.

Tractions are calculated with the Fourier Transform Traction Cytometry (FTTC) method [13]. Deformations () and tractions () are related by the convolution of the traction vector field and a Greens tensor K: (1)

In the case of a linearly elastic semi-infinite substrate, K is given by the Boussinesq equations [27]. Inverting Eq 1 and solving for the tractions is difficult in real space. However, by exploiting the convolution theorem, the equation simplifies to a multiplication in Fourier space: (2) where , and are the Fourier transforms of the deformation field, the traction field and the Greens tensor. The latter can be found in [13].

Eq 2 can be analytically solved and thus allows for the direct calculation of tractions in Fourier space. Tractions in real space are then obtained by applying the inverse Fourier transform.

One particular challenge of Traction Force Microscopy is that noise in the deformation field can lead to large errors in the traction field. This can be remedied by regularization of the reconstructed forces, e.g. by adding the L1 or L2 norm to the cost function of the inverse minimization problem [28]. By contrast, pyTFM does not use explicit regularization but instead smooths the calculated traction field with a user-defined Gaussian low-pass filter, with a sigma of typically 3 μm. This effectively suppresses all tractions with high spatial frequencies but inherently limits the spatial resolution of the pyTFM algorithm. The appropriate degree of smoothing depends on the spatial resolution of the deformation field, which in turn depends e.g. on the density of fiducial markers, the window size for the PIV algorithm or image noise. The user is encouraged to test different values for sigma and to select the smallest value for which the noise in the cell-free areas is still tolerable in comparison to the magnitude of cell tractions.

The original TFM algorithm assumes that the underlying substrate is infinitely thick, which is justified in the case of single cells with dimensions that are smaller than the thickness of the elastic substrate. In the case of cell patches, however, this assumption is inadequate. We have therefore included a correction term for finite substrate thickness [29].

2.2.2 Monolayer stress microscopy.

Stresses in a cell sheet are calculated with an implementation of Monolayer Stress Microscopy as described in [16, 17]. For computing stresses in small cell patches, we implemented a method that corrects for the limited spatial resolution of unconstrained TFM, which otherwise would lead to a substantial underestimation of stresses [20]. Details of this correction are described below.

In the absence of inertial forces, tractions and stresses are balanced according to the relation: (3) where σ_xx, σ_yy are the normal stresses in x- and y- direction, σ_yx is the shear stress, and t_x and t_y are the x- and y-components of the traction vector. This differential equation is solved using a Finite Element method (FEM) where the cell patch is modeled as a 2-dimensional network of nodes arranged in a grid of quadrilateral elements. Each node in the FEM-grid is loaded with a force of the same magnitude but opposing direction as the local tractions. In the standard FE method, the nodal displacements of the cell patch are calculated by solving the equation (4) where are the vector of nodal forces, and K⁻¹ is the inverse stiffness matrix. The nodal displacements are converted to strains by taking the derivative in x- and y-direction. Then, the strain is used to calculate the stress from the stress-strain relationship of a linearly elastic 2-dimensional material: (5) where E and v are Young’s modulus and Poisson’s ratio of the material, and ϵ₁₁, ϵ₂₂ and ϵ₁₂ are the components of the strain tensor. Most of the FEM calculation is performed using the solidspy Python package [30].

The stiffness matrix K in Eq 4 depends on the Young’s modulus in such a way that the Young’s modulus in Eq 5 cancels out. The traction-stress relation is therefore independent of the Young’s modulus of the cell patch [17]. Furthermore, the Poisson’s ratio has only a negligible influence on the stress prediction [17]. In the pyTFM algorithm, the Young’s modulus is set to 1 Pa, and the Poisson’s ratio is set to 0.5.

The FEM algorithm assumes that there are no torques or net forces acting on the cell patch. This must also be true in reality as the cell patch would not be stationary otherwise. However, the TFM algorithm only ensures that forces and torques are globally balanced (across the entire image), but not necessarily across a cell patch. These unbalanced net forces and torques acting on a cell patch must be corrected prior to performing the FEM algorithm to accurately compute the cellular stresses. pyTFM corrects unbalanced net forces by subtracting the sum of all force vectors of the FEM system from the force vector at each node. The unbalanced net torque is corrected by rotating the direction of all force vectors by a small angle (typically below 5°) until the torque of the entire system is zero.

By constraining the FEM system to zero rigid rotational or translational movement, Eq 4 is uniquely solvable [31]. These constraints can be applied in two ways: The first option is to apply a boundary condition of zero displacement in x- and y- direction to an arbitrarily chosen node of the FEM grid, and a boundary condition of zero rotation between this fixed node and another arbitrarily chosen second node [17]. In practise, this is implemented by selecting a second node with the same y-coordinate as the fixed node and applying a zero x-displacement boundary condition. The second option, which is implemented in pyTFM, is as follows: Instead of subjecting individual nodes to displacement boundary conditions, we formulate zero rigid displacement and rotation conditions on the whole system in three separate equations (Eqs 6–8), add them to the system of equations in Eq 4 and finally solve the combined system numerically using a standard least-squares minimization. Eqs 6 and 7 ensure that the sum of all nodal displacements in x- and y-direction is zero, and Eq 8 ensures that the rotation of all nodes around the center of mass of the FEM system is zero. (6) (7) (8)

r_x and r_y are the components of the distance vector of the corresponding node to the center of mass of the FEM-grid. Note that the FEM algorithm described above models the cell sheet as a linear elastic material. However, in the future, other FEM algorithms suitable for non-linear elastic materials could also be applied.

The analysis of stresses in small cell patches poses a second challenge: The FEM-grid should be of the same size and shape as the cell patch, as outside nodes add additional stiffness, leading to an underestimation of the stress field. However, the limited spatial resolution of both PIV and TFM implies that some forces generated close to the edge of the cell patch are predicted to originate from outside the cell patch. Neglecting these forces would lead to an underestimation of the stress field. This can be avoided by extending the FEM-grid by a small margin so that all cell-generated forces are included in the analysis. In practice, the user outlines the area with clearly visible tractions (red outline in Fig 2), over which pyTFM then spans the FEM-grid. We explain further details of this approach in Section 3.1.1.

2.2.3 Limits of applicability of monolayer stress microscopy and TFM.

The TFM and Monolayer Stress Microscopy algorithms can only be applied if a number of conditions are met. 2-dimensional TFM relies on the assumption that tractions in z-direction generate only small deformations in the x- and y-plane. This is valid if z-tractions are small, or if the substrate is almost incompressible (Poison’s ratio close to 0.5) [11] Additionally, TFM assumes that the matrix is a linearly elastic material. Both assumptions are valid for polyacrylamide and PDMS, two popular substrates for TFM [32–35].

For Monolayer Stress Microscopy, cells are modeled as a linearly elastic material with homogeneous and isotropic elastic properties. Although many cell types have been shown to exhibit stress stiffening (the Young’s modulus increases with the cellular stress [36]), this has only a small effect on the stresses recovered by Monolayer Stress Microscopy [17].

Furthermore, Monolayer Stress Microscopy models the cells as a 2-dimensional flat sheet. Deviations from this assumption can introduce an error to the stress calculation on the order of (h/l)² [17], where h is the cell height and l is the wavelength of the tractions in Fourier-space (in the case of a single cell that generates tractions in the form of two opposing force monopoles, l corresponds to the distance between the force monopoles). This error can become relevant for isolated round cells but not for larger flat cell colonies.

3.1.1 Effect of FEM-grid area on stress recovery.

pyTFM requires the user to select an area of the traction field over which pyTFM then computes contractility and strain energy and draws the FEM-grid for computing monolayer stresses. The size of this area influences the accuracy of the stress and force measurements. Selecting an area that is too small leads to an underestimation of stresses and contractility. Selecting an area that is too large also leads to an underestimation of stresses. To systematically analyze which effect the size of the user-selected area has on the traction and stress reconstruction, we expand the traction area and analyze the average normal stress and the contractility for the synthetic test data described above (Fig 4A) and for an MDCK cell patch grown on a polyacrylamide substrate (Young’s modulus 49 kPa, Fig 4B). In the case of the synthetic data, we normalize the computed average normal stress and contractility to the known input stress (1 N/m) and to the known contractility of the input traction field (600 N), respectively. In the case of the experimental data, we normalize the computed average normal stress and contractility to their respective maximum values as the true stress and contractility is unknown.

The influence of FEM-grid area on the stress calculation may also depend on the size of the cell patch. To analyze this effect, we reduce the cell patch size in the synthetic test data (Fig 4A) and calculate the average mean normal stress (Fig 4C). The margin of expansion for each cell patch is chosen so that we recover the maximal stress (the optimal margin of expansion Δx is between 4-5 μm, independent of patch size, as it is mainly determined by the spatial accuracy of the TFM algorithm).

We find that the normalized stress rapidly increases (by approximately 40%) with increasing area until it reaches a maximum, after which it declines at a slower rate. The contractility displays a similar initial increase but then remains approximately constant. The maximum of the normalized stress occurs when the traction area just covers all cell-generated tractions, including those that appear outside the cell patch. In the cases of the synthetic data, the maximum is reached at a traction area expansion distance of 5 μm beyond the cell patch outline, whereas in the case of the MDCK cell patch, it is reached at at an expansion distance of 20 μm. The reason for this larger distance in the MDCK data is the additional blurring of tractions introduced by the PIV algorithm (whereas no PIV was needed for analyzing the synthetic data). The traction area corresponding to the maximal normal stress can be regarded as the optimum, as approximately 93% of the input stress is recovered. Expanding the traction area and thus the FEM-grid beyond the optimum distance adds elastic material to the monolayer and thereby reduces the average stress. This stress reduction, however, occurs only gradually (Fig 4B), which implies that in practice it is best to choose the traction area rather generously to include all cell-generated tractions. The contractility reaches its maximum values at almost the same expansion distance as the stress. Thus, it is possible to use the same area to accurately compute both stress and contractility.

The stress recovery also depends on the cell patch sizes (Fig 4C). The recovered mean normal stress increases rapidly with patch size but reaches a plateau above 50 μm so that more that 90% of the input stress is recovered. Expanding the grid area by 4-5 μm, as noted above, always results in a higher stress recovery compared to an FEM-grid area that matches the patch area (Fig 4c red vs. turquoise curves).