3D Traction Forces in Cancer Cell Invasion

Cell invasion through a dense three-dimensional (3D) matrix is believed to depend on the ability of cells to generate traction forces. To quantify the role of cell tractions during invasion in 3D, we present a technique to measure the elastic strain energy stored in the matrix due to traction-induced deformations. The matrix deformations around a cell were measured by tracking the 3D positions of fluorescent beads tightly embedded in the matrix. The bead positions served as nodes for a finite element tessellation. From the strain in each element and the known matrix elasticity, we computed the local strain energy in the matrix surrounding the cell. We applied the technique to measure the strain energy of highly invasive MDA-MB-231 breast carcinoma and A-125 lung carcinoma cells in collagen gels. The results were compared to the strain energy generated by non-invasive MCF-7 breast and A-549 lung carcinoma cells. In all cases, cells locally contracted the matrix. Invasive breast and lung carcinoma cells showed a significantly higher contractility compared to non-invasive cells. Higher contractility, however, was not universally associated with higher invasiveness. For instance, non-invasive A-431 vulva carcinoma cells were the most contractile cells among all cell lines tested. As a universal feature, however, we found that invasive cells assumed an elongated spindle-like morphology as opposed to a more spherical shape of non-invasive cells. Accordingly, the distribution of strain energy density around invasive cells followed patterns of increased complexity and anisotropy. These results suggest that not so much the magnitude of traction generation but their directionality is important for cancer cell invasion.

1. Characterizing the mechanical properties of the ECM; 2. Measuring the cell-induced deformation in the ECM; 3. Computing the strain energy stored in the deformed ECM.
The following sections give details for each step.

Mechanical properties of collagen gels
Reconstituted collagen gels show a fibrous networks structure (Fig. S1a). At the collagen concentration used here (2.4 mg/ml collagen), the average pore size of the network is 1.3 µm [1]. Fluorescent microbeads of 1 µm diameter were dispersed in the gels prior to collagen polymerization and served as fiducial markers. These beads are trapped in the gel by tightly binding to collagen fibers. The average bead-bead distance is (24 ± 8) µm (Fig. S1b). At this length scale, we can neglect the porous or fibrous microstructure of the collagen and treat the collagen gels as a continuum with isotropic properties [1].
The elastic properties of the collagen gels were measured with a rate-controlled cone-plate rheometer (AR-G2, TA Instruments, New Castle, DE, USA) by applying oscillatory shear strain with an amplitude of 5%. Gels were prepared as described in Methods but allowed to polymerize in situ (between the rheometer plates) to ensure good contact with the probing cone-plate. Measurements were carried out at a temperature of 37°C. The frequency response of the complex modulus shows predominantly elastic behavior at all frequencies measured (Fig. S1d). Collagen gels display a nonlinear stress-strain relationship ( Fig. S1c) with prominent strain stiffening at strains above 5%. The method used to extract the strain energy described below ignores non-linear behavior and treats the collagen gels as isotropic linear elastic materials with a shear modulus of G = 118 Pa.

Figure S1: Properties of collagen gels. a) Confocal section through a collagen gel. b) Distribution of bead-bead distances. c) Stress-strain relationship of the collagen network measured in a cone-plate rheometer during a strain ramp (speed: 1%/s). The behavior is approximately linear for strains below 5%; d) Frequency response measured in a cone-plate
rheometer. The amplitude of sinusoidal oscillations was 5%. The storage modulus G' is 118 Pa at a frequency f of 1 Hz. Storage modulus G' and loss modulus G" increased weakly with frequency according to a power-law with exponent 0.08, indicative of predominantly elastic behavior. Figure S2: Intensity distribution in a subvolume around a bead at two distinct time points. To track the bead position, the subvolume at t=2 is shifted until a best match with the subvolume at t=1 is achieved.

Measurement of the displacement field
Matrix displacments are inferred from measurements of the 3D positions of fluorescent beads that are scattered throughout the gel. Optical sections (z-stacks) are taken throughout the depth of the gel at 2 µm intervals at two or more time points. The 3D positions of the beads and their displacements over time can then be determined as follows: In a first step, the positions of the beads during the first time point are determined from the positions of all pixels in the image stack with a local intensity maximum above a specific intensity threshold. Second, all pixels in a 7x7x7 pixel neighborhood ((4.5 x 4.5 x 14) µm) around these local maximum pixels are selected. The displacement vector of each bead is determined using a 3D difference-with-interpolation algorithm that shifts the 7x7x7 neighborhood pixels around each bead from the second image stack by sub- This procedure yields subpixel resolution. The accuracy of the method was estimated by repeatedly measuring bead positions in a cell-free collagen gel and found to be 22 nm in the x-and y-direction, and 130 nm in the z-direction (20x 0.4 NA objective, 0.5x optical coupler, z-section distance of 2 µm). This section outlines the computation of the strain energy. More details can be found elsewhere [2].

Calculation of strain energy
In a continuum material, the elastic strain ε at any spatial position can be calculated from the displacements u according to ( If we assume isotropic linear elastic properties, the stress σ is related to strain ε by Hook's law where D is the material stiffness matrix. D is fully defined by two material constants, e.g. Young's modulus and Poisson's ratio. The strain energy U due to elastic deformation is then the volume integral of the inner product of stress and strain, which according to Eq. (3) takes the form To solve the above integral, the deformation of the material everywhere in the measurement volume needs to be known. However, we only know the material displacements at the discrete locations of the beads. Displacements between bead positions are therefore obtained by tri-linear interpolation through a Delaunay tessellation [3] with linear tetrahedral finite elements. Thus, four bead positions serve as nodes for a tetrahedral element as illustrated in Fig. S3. For convenience, we write the 12 displacement components (4 nodes with 3 displacement components in x-, y-, and z-direction) as a These strain components are then calculated from the displacements in Eq. (5) as with the strain-displacement matrix B containing the spatial derivatives of the shape functions N.
Inserting Eq. (9) in Eq. (4) gives the discretized form for the strain energy where the volume integral defines the stiffness matrix K for the elements. The strain energy of one element can therefore be calculated from the measured displacements of four beads by The method outlined in this section was implemented in Matlab (MathWorks, Natick, MA, USA) [4].

Accuracy of strain energy measurements
Sources of error A main source for error in the calculation of strain energy is the measurement noise in the bead displacements. How displacement noise is translated into erroneous strain energy is governed by two element properties, its volume and its shape, as explained below. It can be shown that the strain energy density of a linear tetrahedral element U due to given displacements (e.g. noisy bead displacements) scales inversely with the volume of the element by a power law with exponent 2/3: Consequently, increasing the number of beads in the measurement volume with the aim of improving the resolution of the strain energy also leads to elements with smaller volumes and therefore increased levels of measurement noise.
The optimal shape of a tetrahedron in terms of insensitivity to displacement noise is a regular tetrahedron with all six edges of equal length. After tessellation of randomly distributed beads in a given volume, however, the resulting tetrahedra can have all possible shapes including degenerate forms such as an almost flat tetrahedron. Flat tetrahedral elements are exceptionally prone to erroneous strain energy as small displacement errors can translate into large strains. To characterize the shape of tetrahedral elements, we determine the ratio of the volume of the element to the sum of its cubed edge lengths, normalized to range between 0 and 1: Thus, a regular tetrahedron has a shape of q=1. The value of q decreases with increasingly degenerate forms (Fig. S4a-c). In a typical collagen gel, tetrahedra of all shapes can be found, with shape factors around 0.5 being the most frequent (Fig. S4d). The dependence of erroneous strain energy on the shape factor was estimated from the measured positions of ~40,000 beads in a typical collagen gel, superimposed by random Gaussian displacement noise according to our experimentally determined values of 22 nm (x,y-direction) and 130 nm (zdirection). We find that Gaussian displacement noise causes an erroneous strain energy density with approximately log-normal distribution (Fig. S5a). The width of the distribution and the expected value of the erroneous strain energy increase dramatically as q falls below 0.5 ( Fig. S5a-b). U of a tetrahedron with a shape factor of 0.1) with a log-normal probability density distribution. The expected value and the width of the distribution depend on the shape of the tetrahedra. As the shape quality factor gradually decreases, the distribution widens and the expected value increases. The inset shows the geometric mean of the erroneous strain energy as a function of the shape quality factor q. For values of q<0.5, the tetrahedral elements become excessively sensitive to noisy displacements.

Error correction for individual elements
In the following, we describe how the erroneous strain energy from measurement noise contributes to the total strain energy, and how this error can be corrected for.
The measured displacements m u are normally distributed with mean µ and standard deviations (rms measurement noise) of x σ , y σ and z σ in the x-, y-, and z-direction, respectively. The expected value of the strain energy ) ( m U u is the sum of the true strain energy and a term U ∆ that represents the erroneous strain energy due to measurement noise in the bead positions and that is independent of the mean displacements.

Strain energy error correction
The following strategy was employed to correct for errors in the total strain energy: 1. Elements with a shape factor of 5 . 0 < q are discarded. 2. A strain energy error evaluated with Eq. (15) and Eq. (16) is subtracted from the strain energy values of an element. 3. Gaps in the strain energy density distribution resulting from step 1 above are filled by nearest neighbor interpolation.
The impact of this error correction on the strain energy density distribution around a cell is illustrated in Fig. S6. The strain energy density distribution is plotted for each of the various correction steps described above. While the strain energy density seems little affected by the correction steps, the importance of error correction becomes evident when the total strain energy is computed by integration of the strain energy density over increasingly larger volumes, starting from the cell center (Fig. S6d). In all cases, i.e. without or with error correction, the strain energy rises rapidly at first, followed by an increase with smaller slope that is sensitive to noise. Because normally-distributed displacement noise results in a log-normal distribution of the erroneous strain energy in tetrahedral elements with similar shape factors (Fig. S5), the error subtraction approach (Eqs. (15) and (16)) is necessarily incomplete.

Strain energy in a fixed volume around two representative cells evaluated on a subset of beads, i.e. the number of beads in an experimental data set was systematically reduced (a, b, and c show resulting finite element meshes) to obtain the curves in d.
In a linear continuum material and under noise-free conditions, the accuracy for estimating the strain energy is limited only by the interpolation accuracy within the tetrahedral elements. Using an increasing number of beads would enable mapping out increasingly finer details of the strain energy distribution until eventually convergent results are obtained. In practice, there exists an upper bound for the number of beads that can be dispersed in the gel, because of the fibrous nature of the material and due to optical effects, e.g. blur from out-of-focus beads in a conventional epifluorescence microscope.
To estimate the sensitivity of the strain energy to the number of beads used, we performed our analysis on a subset of beads with varying densities as illustrated in Fig. S7. A number of beads were randomly drawn from the measured data set, and the strain energy was computed on this subset. As expected, the strain energy around a cell increases with the number of beads (solid curves). With a 10-fold higher bead density, the total strain energy approximately doubled. The strain energy did not saturate even at the highest bead density, however, suggesting that the method underestimates the true strain energy, although part of the raise in strain energy with increasing bead density can be explained by an incomplete noise correction.

Experimental verification of a continuum mechanical description of collagen gels
To test whether the continuum assumption of the collagen gel holds at the length scale of typical bead-to-bead distances (24 µm), we deformed the gel with a superparamagnetic bead that was embedded in the collagen gel and laterally forced with magnetic tweezers. The collagen gel was prepared according to our protocol but with the addition of superparamagnetic beads of 4.5 µm diameter. A 20 µl droplet of the mixture was placed on a cell culture dish and allowed to polymerize upside down so that the magnetic beads drifted to the surface of the droplet. After polymerization, we recorded z-stacks through the gel before and after the application of a 58 nN force acting on an isolated magnetic bead. Force calibration of the magnetic tweezer setup was performed as explained previously [5]. The displacement vectors of the fluorescent beads surrounding the magnetic bead between the force-on and force-off states are shown in Fig. S8a-b. The displacements in a continuum linear elastic half-space due to a point force acting on the surface are given by the Boussinesq solution [6] and were fitted to the measured displacements of the fluorescent beads. The measured bead displacements closely follow the Boussinesq displacement field (Fig. S8a), confirming that the deformation field of the collagen gel can be described by a continuum approach. The Poisson's ratio enters the fit as an independent parameter and was found to be 0.35; this value was used for computing the strain energy in all cell experiments.

Experimental verification of strain energy computation
To quantitatively assess the accuracy of the measured strain energy, we computed the strain energy in the collagen droplet surrounding the magnetic bead between the force-on and force-off states. The strain energy density around the magnetic bead is shown in Fig. S8c and falls off rapidly with increasing distance from the magnetic bead, as expected. The integrated total strain energy stored in the gel was measured to be 0.08 pJ. This value agrees well with the expected strain energy of 0.1 pJ, calculated from the work done on the magnetic bead as it moved towards the magnetic needle (force acting on the bead times bead displacement times 0.5 according to Clapeyron's theorem [7]).
Furthermore, we carried out a simple indentation experiment. A spherical steel ball of 100 µm diameter was placed on the surface of a collagen gel. The ball indented the gel by 36.1 µm due to a gravitational force of 35.4 nN. The displacements of the measured bead positions between the undeformed and deformed positions are plotted in Fig. S9b. Bead displacements are projected on a cylindrical coordinate system with the axis going through the center of the ball. Because bead displacements larger than 10 µm could not be reliably tracked, we determined the strain energy in a donut shaped subvolume of the gel as indicated by the red outline in Fig. S9b and Fig. S9c. We have furthermore calculated the expected strain energy of the same volume from a finite element simulation of the experiment as illustrated in Fig. S9a, with the measured indentation depth as a prescribed condition for the displacement of the ball. Due to symmetry, the finite element model geometry was constructed with a planar section through the center of the ball together with cylindrically symmetric boundary conditions and fixed walls. Contact between the ball and the gel was modeled with a no-slip condition. The simulation was carried out with Abaqus (SIMULIA, Providence, RI, USA). The displacement results of the simulation were then transferred to the 3D geometry of the experiment and the strain energy was calculated with our algorithm. The result of the simulation was 1.5 pJ and agrees well with that of the experimentally obtained strain energy of 1.6 pJ.

Total strain energy, cell shape anisotropy and strain energy anisotropy follow a log-normal distribution
Cumulative probabilities of all measured cell data were evaluated and are shown in Fig. S10 for strain energy values, in Fig. S11 for cell shape, and in Fig. S12 for strain energy anisotropy indices. The latter two measures are bounded by 1 and this bias was taken into account before computation. The cumulative probability assuming an underlying log-normal distribution is plotted for comparison showing that the distribution of total strain energy among individual cells closely follows a lognormal distribution.

Time-lapse measurements of strain energy in 3D collagen gels
The technique for measuring the 3-D strain energy can be readily combined with time-lapse recordings by repeatedly taking image stacks around cells embedded in 3D collagen gels. A last stack imaged after the release of cell tension (e.g. by treatment with cytochalasin D) serves as the reference stack for obtaining the bead positions at the force-free state. With this approach, further insight can be gained into the dynamics of cellular traction generation during cell migration.
In the following we present two examples, one of a highly contractile yet non-invasive vulva carcinoma cell (A-431), and one of an invading breast carcinoma cell (MDA-MB-231). In both cases, cells have been mixed to the collagen solution prior to polymerization and allowed to spread for 24h as described previously. Image stacks were then taken at time intervals of 15 min and 20 min, respectively. To reduce phototoxicity, light intensities were minimized as far as possible, and the exposure shutter was closed between imaging.
During the measurement period of 165 min before force release, the vulva carcinoma cell kept the surrounding collagen matrix under high tension at all times but the total strain energy fluctuated by up to 25% around the time average (Fig. 4a).
The fluctuations of total strain energy are caused by fluctuations of strain energy density that varied both temporally and spatially around the cell. An example of the temporal and spatial distribution of strain energy fluctuations is shown in Fig. S13-S14 and Movie S3. The cell moved during the measurement period in a random manner but remained within a few µm of its starting position. In concert with this random cell movement, the strain energy density distribution around the cell shifted as well. This is illustrated by contour lines of the local change in strain energy density between two time steps t 1 and t 2 , computed as ) (  Negative values with corresponding contour lines towards the blue end of the spectrum point to regions where tension is relaxed between successive time steps. Contour lines towards the red end of the spectrum point to regions where tension is increased between successive time steps. In the case of the non-invasive vulva carcinoma cell, we observe periods of comparatively small changes (Fig.S13a,b and Fig.S14b,e) and periods of comparatively large changes in the strain energy density distribution (Fig.S13c-f and Fig.S14a,c,d). It is interesting to note that areas of increasing tension at one time step generally turn into areas of decreasing tension in the following time step. In summary, this non-invasive carcinoma cell keeps the collagen gel under high tension all time steps, with temporally fluctuating regions of high and low strain energy density that shift the cell around in a random, non-persistent fashion during seemingly uncoordinated cycles of tension generation and release.
In the second example, an invasive breast carcinoma cell also keeps the collagen gel under high tension during the measurement period of 360 min (Fig.4b), but with larger fluctuations (more than 50% around the time average) compared to the non-invasive cell. Moreover, the cell migrated through the collagen in all three spatial directions by a distance spanning more than 100 µm (Fig. S15-S17 and Movie S4).
Compared to the non-invasive cell, the temporal and spatial pattern of strain energy density fluctuations shows two striking differences. First, the changes in strain energy density are considerably more persistent in time, meaning that a local increase during one time step is usually not followed by a corresponding decrease over the following time step. Second, local changes in strain energy density often preceded shape changes of the cell and are followed by a cell movement in the same direction for positive strain energy changes, or in the opposite direction for negative changes.

Software
The algorithms presented in this paper implemented in MATLAB (MathWorks, Natick, MA, USA) are available as Supporting Information File S1.