A Bayesian traction force microscopy method with automated denoising in a user-friendly software package
aa r X i v : . [ q - b i o . CB ] M a y A Bayesian traction force microscopy method withautomated denoising in a user-friendly software package
Yunfei Huang a , Gerhard Gompper a , Benedikt Sabass a, ∗ a Theoretical Soft Matter and Biophysics, Institute of Complex Systems and Institute forAdvanced Simulation, Forschungszentrum Juelich, 52425 Juelich, Germany
Abstract
Adherent biological cells generate traction forces on a substrate that playa central role for migration, mechanosensing, differentiation, and collectivebehavior. The established method for quantifying this cell-substrate inter-action is traction force microscopy (TFM). In spite of recent advancements,inference of the traction forces from measurements remains very sensitiveto noise. However, suppression of the noise reduces the measurement accu-racy and the spatial resolution, which makes it crucial to select an optimallevel of noise reduction. Here, we present a fully automated method fornoise reduction and robust, standardized traction-force reconstruction. Themethod, termed Bayesian Fourier transform traction cytometry, combinesthe robustness of Bayesian L2 regularization with the computation speed ofFourier transform traction cytometry. We validate the performance of themethod with synthetic and real data. The method is made freely availableas a software package with a graphical user-interface for intuitive usage.
Keywords:
Traction force microscopy, Bayesian inference.
1. Introduction
Traction force microscopy (TFM) is a technique for measuring surfacetraction forces on an elastic substrate. Due to the unique possibilities offeredby the technique, TFM enjoys wide popularity among biologists, materialsscientists, and experimental physicists, see Refs. [1, 2, 34, 3, 4, 5, 6] for a non- ∗ Corresponding author.
E-mail address: [email protected] U describing thedeformation of the substrates is usually done by comparing images of thefluorescent markers before and after removal of the adherent cells. Estab-lished methods for calculating the displacement fields include particle imagevelocimetry, single particle tracking, or optical flow tracking [34, 2, 36]. Hereand in the following, positions on the two-dimensional surface are describedin a Cartesian system with coordinates x = ( x , x ). Using the planar defor-mations ( U , U ) as input, a spatial map of traction ( F , F ) on the gel surfacecan be mathematically reconstructed if out-of-plane forces are assumed to benegligible. For the reconstruction, the substrate is typically assumed to bea homogeneous, isotropic and linearly elastic half-space. Thus, the relationbetween the continuous displacement field U i ( x ) and the traction force field2 j ( x ′ ) on the surface of substrates can be expressed as [27] U i ( x ) = Z Ω 2 X j =1 G ij ( x − x ′ ) F j ( x ′ ) d x ′ , (1)where G ij ( x ) is the Green’s function and Ω covers the whole substrate sur-face. The traction forces can be calculated in real space with finite elementmethods [28, 29] or boundary element methods [26, 30]. To calculate thetractions based on Eq. (1) numerically, the integral equation needs to be dis-cretized. In real space, one can employ linear shape functions [26, 30, 31] towrite Eq. (1) as a linear matrix equation u = Mf , where the lower case letters u and f denote the discretized displacements and tractions. The coefficientmatrix M results from integration of the shape functions. Such real-spacemethods are very flexible since they permit the study of various linear ma-terial responses encoded in the Green’s function and spatial constraints areeasily incorporated. However, accurate construction of the matrix M requiressignificant computation time on desktop machines. Alternatively, Eq. (1) canbe solved in Fourier space by making use of the convolution theorem. Thisapproach is called Fourier transform traction cytometry (FTTC) [8, 26]. Weemploy a spatial wave vector k = ( k , k ) with absolute value k = | k | . Instandard FTTC, the integral Eq. (1) is written as ˜ u i k = { P j ˜ G ij ˜ f j } k , wherethe tilde denotes the Fourier-transformed quantity. Using a matrix formula-tion analogous to the real-space expression, we have ˜ u = ˜ M ˜ f with ˜ M havinga tri-diagonal structure. For conceptual clarity, in the following we will writethe measurement noise in the recorded displacement explicitly as s in thereal-space domain and as ˜ s in Fourier space. This noise can be estimated inthe experiment by quantifying the variance of the measured displacementsin absence of traction. The discretized equations then read ( u = Mf + s in real space , ˜ u = ˜ M ˜ f + ˜ s in Fourier space . (2)For traction force microscopy, either of these equations is employed to deter-mine the tractions f . The removal of noise is critical in most TFM methods.In real-space TFM calculations, the condition number of M , defined as theratio of the largest singular value to the smallest, is almost always muchlarger than unity, typically above 10 . M is therefore ill-conditioned whichimplies that small noise produces drastic changes in the calculated traction3orces. For FTTC, spatially varying random noise occurs mainly at high spa-tial wave numbers. Hence, noise suppression can be achieved by suppressinghigh frequency data. In Ref. [31], we systematically tested a variety of trac-tion reconstruction approaches based on Eq. (2). The standard approach forsolving the equation in real space is L2 regularization [32, 26, 33, 2], which in-vokes a penalty on the traction magnitude to robustly suppress the effects ofnoise. With Fourier space methods, a low-pass filter is frequently employed tosuppress noise in the displacement field before direct inversion of Eq. (2) [8].Alternatively, Fourier-space traction reconstruction can also be combinedwith L2 regularization, which conveys additional robustness [26, 34, 35, 36].Virtually all standard methods for traction calculation require the implicitor explicit choice of a parameter that suppresses noise and leaves as muchof the true signal conserved as possible. Within a Bayesian framework, thisparameter choice can be rationalized by relating filter- or regularization pa-rameters to prior distributions that represent prior knowledge about the data.Maximizing the likelihood of the prior distributions yields the correspondingoptimal trade-off between noise suppression and faithful data reconstruction.Bayesian regularization has been used for example in astrophysics [37, 38]and mechanical structure monitoring [39]. For inference of internal stressin a cell monolayer, an iterative maximum a posteriori estimation has beenemployed [14]. TFM with Bayesian L2 regularization (BL2) was introducedin Ref. [31] and is based on an established framework of Bayesian fitting [40].Bayesian L2 regularization was first employed for real-space TFM methodssince this variant allows comparison of a broad variety of approaches. Forpractical applications, however, calculations in Fourier domain have signifi-cant advantages in terms of robustness and speed. In this work, we presentthe corresponding method that we term Bayesian Fourier transform trac-tion cytometry (BFTTC). We compare BFTTC with other methods such asclassical L2 regularization, Bayesian L2 regularization in real-space (BL2),and regularized Fourier transform traction cytometry (FTTC). We find thatBTTC is a computationally fast method that provides robust traction calcu-lations without requiring manual adjustment of the noise-suppression level.We also provide a Matlab software package for BFTTC that is freely avail-able for download. This software package is intended to provide a simple androbust solution for data analysis in the hands of experimentalists. A graph-ical user interface allows intuitive use of the program and little theoreticalbackground knowledge is required. 4 . Methods and implementation Assuming that the substrate is a semi-infinite half-space, the Green’sfunction in Eq. (1) is given by the standard expression G ij ( x ) = (1 + ν ) πE (cid:20) (1 − ν ) δ ij r + ν x i x j r (cid:21) , (3)where E and ν represent the Young modulus and Poisson ratio, respectively.Here, r = | x | and δ ij is the Kronecker delta function. Denoting the wavevector by k = ( k , k ) with absolute value k = | k | , the Fourier-transformedGreen function is given by˜ G ij k = 2(1 + ν ) E (cid:20) δ ij k − νk i k j k (cid:21) . (4)The continuous traction and displacement fields are discretized by rectan-gular meshes where m and n are the number of discretization nodes fortractions and displacements respectively. In the discretized equations (2),the size of the displacement vector u is 2 m × f is 2 n × m = n . The classical approach to solve Eq. (2) is L2 regularization, which is alsocalled Tikhonov regularization or ridge regression. L2 regularization is arobust procedure that suppresses noise and produces a smoothed tractionfield [31]. Here, the residual k u − Mf k = ( Mf − u ) T ( Mf − u ) is minimizedtogether with the solution norm λ k f k = λ f T f . The factor λ is calledregularization parameter. The reconstructed traction ˆ f satisfiesˆ f = argmin f (cid:2) k Mf − u k + λ k f k (cid:3) . (5)This approach can be employed for real-space TFM and in Fourier space,where the square norms can be calculated conveniently with Parseval’s theo-rem. The proper choice of the regularization parameter λ is critical in caseaccurate traction calculations are required. A popular heuristic for choosing5he regularization parameter is based on a double-logarithmic plot of the so-lution residual vs. the traction norm for varying λ . Often, the plotted lineresembles an “L” shape and the regularization parameter is chosen to lie inthe corner of this curve, thus providing a trade-off between faithful recon-struction and smoothing [41]. However, this “L-curve criterion” is often oflittle use, in particular when the corner is either absent or cannot be localizedprecisely on the double-logarithmic scale. Moreover, it has also been shownthat the L-curve criterion can fail systematically [43, 42]. Therefore, theL-curve criterion is often complemented with other methods for finding theregularization parameter, such as cross-validation [31]. In any case, a manualparameter variation is mandatory to check the validity of the solution. Bayesian methods can be used to regularize data in a systematic andautomated way. Our approach is based on an established iterative inferenceprocedure [40]. In the first step, a model is fitted to the data. In the secondstep, the evidence for the chosen model is calculated. Traction computationswith Bayesian L2 regularization (BL2) were first introduced as a real-spaceapproach in Ref. [31]. Here we describe the adaptation of this method toFourier-space traction calculation. It is assumed that the noise s in Eq. (2)has a Gaussian distribution with vanishing mean and a variance of 1 /β .Therefore, given a traction vector f , the likelihood of measuring a particular2 m × u is p ( u | f , β ) = exp[ − βE u ] Z u = exp[ − β ( Mf − u ) T ( Mf − u ) / Z u , (6)where Z u = (2 π/β ) m . As a prior distribution for the 2 n × f we choose a Gaussian distribution with variance 1 /α as p ( f | α ) = exp[ − αE f ] Z f = exp[ − α f T f / Z f , (7)where Z f = (2 π/α ) n . According to Bayes’ rule, the posterior distribution of f is given by p ( f | u , α, β ) = p ( u | f , β ) p ( f | α ) p ( u | α, β ) = exp[ − K ( f )] p ( u | α, β ) Z u Z f , (8)where K ( f ) = βE u + αE f and p ( u | α, β ) = R d n f exp[ − K ( f )] / ( Z u Z f ). Tofind the traction vector with the highest posterior probability, we maximize6 ( f | u , α, β ) with respect to f . The calculation yields f MP = argmin f (cid:2) β k Mf − u k / α k f k / (cid:3) , which is equivalent to our formula for L2 regularization,Eq. 5, when λ = α/β [34].Next, the values of the hyperparameters α and β have to be deter-mined. In principle, both values can be found by maximizing the evi-dence p ( α, β | u ) that depends on the measured displacements u . However,the noise variance 1 /β can also be estimated directly from the measure-ment uncertainty. Thereby, the maximization of p ( α, β | u ) can be reducedto a robust one-dimensional search for the optimal value of α . Bayes’ lawyields p ( α, β | u ) = p ( u | α, β ) p ( α, β ) /p ( u ). We next assume a uniform prior p ( α, β ) ≃ const . and note that p ( u ) does not play a role for the optimization.Thus, we only need to maximize p ( u | α, β ) ∼ R d n f exp[ − K ( f )] with respectto α . The integral can be analytically calculated by completing the square.On defining A = α I + β M T M one finds p ( u | α, β ) = R d n f exp[ − K ( f )] Z u Z f = (2 π ) n (det A ) − / Z u Z f exp[ − K ( f MP )] . (9)Since f MP and A both depend on α , the maximization of Eq. (9) with re-spect to α needs to be done iteratively. This iteration can be sped upby performing the calculations in Fourier space. For notational clarity,we will write Fourier-space variables and derived quantities with a tilde.The Fourier-transformation of f MP yields ˜f MP = ( ˜M † ˜M + α/β I ) − ˜M † ˜u [26],where the complex transpose is indicated by a † . Parseval’s theorem allowsconvenient expression of Eq. (9) through Fourier-space variables. We have˜ E u = ( ˜ M ˜ f − ˜ u ) † ( ˜ M ˜ f − ˜ u ) / (2 m ), ˜ E f = ˜ f † ˜ f / (2 n ), and ˜ A = α I /n + β ˜ M † ˜ M /m .Using these expressions, the logarithm of the evidence, cf. Eq. (9), can bewritten aslog p ( ˜u | α, β ) = − β ˜ E u ( ˜f MP ) − α ˜ E f ( ˜f MP ) −
12 log(det ˜ A )+ n log α + m log β − m log(2 π ) . (10)This expression is evaluated numerically. The calculation of log(det ˜ A ) isdone by a Cholesky decomposition of the positive matrix ˜ A = LL T aslog(det( LL T )) = 2 log Π i L ii = 2Σ i log( L ii ) [31]. To determine the valueof α = ˆ α that maximizes log p ( ˜u | α, β ) we employ a golden-section search.Finally, the L2 regularization parameter follows as ˆ λ = ˆ α/β .The calculation of the parameter value ˆ λ requires a well-defined maxi-mum of the logarithmic evidence as a function of α . To assess whether this7aximum exists, we investigate the condition dd α log p ( ˜u | α, β ) = 0. For eval-uation of the derivatives of ˜ E u ( ˜f MP ) and ˜ E f ( ˜f MP ) we use that n = m and that ˜M commutes with ( ˜M † ˜M + α/β I ) − since the Fourier-transformed Green’sfunction is a real, symmetric matrix. A straightforward calculation yields dd α ˜ E u ( ˜f MP ) = − λ dd α ˜ E f ( ˜f MP ). Therefore, the condition determining the max-imum becomes 0 = dd α log p ( ˜u | α, β ) = − ˜ E f ( ˜f MP ) − n Tr[ ˜ A − ] + nα . We nextperform a symbolic eigenvalue decomposition of ˜M and denote the eigenval-ues by { m i } , the matrix of eigenvectors by V T , and define ˆ u i = V ij ˜ u j . Thecondition determining the maximum of the logarithmic evidence then reads12 n n X i =1 β ˆ u † i ˆ u i m i ( m i + λ ) = 12 n X i =1 m i λ ( m i + λ ) . (11)Solutions exist if the functions of λ on the left hand side and on the righthand side of Eq. (11) cross each other. Both functions decrease monotonouslywith λ . However, for λ → λ ≥
0. Inthe limit of λ → ∞ , the condition for the occurence of a maxium becomes n P ni =1 β ˆ u † i ˆ u i m i / (cid:16)P nj =1 m j (cid:17) ≥
1. For the TFM data, we find that the val-ues of ˆ u † i ˆ u i roughly decrease with decreasing squared eigenvalues m i sincethe displacement magnitudes typically decrease with higher Fourier modes,as do the entries of ˜ M . Assuming that the approximate ordering of m i and ˆ u † i ˆ u i holds strictly, we can invoke Chebyshev’s sum inequality to ob-tain n P ni =1 β ˆ u † i ˆ u i m i / (cid:16)P nj =1 m j (cid:17) ≥ n P ni =1 β ˆ u † i ˆ u i . Since for all reasonableTFM datasets the mean squared displacement is larger than the noise vari-ance, we expect that n P ni =1 β ˆ u † i ˆ u i = βn P ni =1 u i >
1. Therefore, the condi-tion for the occurrence of a maximum in log p (˜ u | α, β ) should be fulfilled forsome λ >
0. The resulting maximum is unique due to defined signs of the λ derivatives of Eq. (11). In summary, a semi-quantitative argument supportsthe existence of a unique maximum of the logarithmic evidence log p (˜ u | α, β )when appropriate TFM data is used. In our tests, a maximum was found forall datasets. To confirm that the Bayesian approach yields a correct estimate for theregularization parameter we employ synthetic data sets with known proper-8ies. In our first test series, we generate random traction fields by drawingindividual traction vectors from Gaussian distributions with fixed variances,as illustrated in Fig. 2 (a-i) and (a-ii). The traction field is produced on a50 ×
50 grid with a Young modulus of E = 10 kPa and a Poisson ratio of ν = 0 .
3. For example, we employ a Gaussian traction distribution with avariance of 10 Pa and therefore α = 10 − Pa − . After calculation of the dis-placements from the traction, Gaussian noise with a variance of 10 − Pix isadded, thus β = 10 Pix − . In the second test series, we construct syntheticdata to study the reconstruction quality for localized traction patterns. Asin previous work [26, 31], we assume that the traction is localized in circu-lar spots, each having a constant traction magnitude. For every individualspot, the step-like traction profile can be integrated analytically to producea displacement field. Due to the linearity of the problem, displacements fromdifferent spots can be added to produce the final result. Explicit formulasfor the displacement field are provided in the supplementary of Ref. [31]. Forgeneration of this data, we fix the Young modulus E = 10 kPa and the Pois-son ratio ν = 0 .
3. The traction patterns consist of 10 −
20 circular tractionspots, as illustrated in Fig. 3(a). The diameter of the spots is 2 µ m and themesh size of the reconstructed traction is 0 . µ m. The traction force magni-tude in the spots is randomly chosen in the range [0 − For the synthetic test data with circular spots the traction force is ex-actly known. Therefore, we can qualitatively calculate the reconstructionerrors. Here, we use four different error measures introduced in our previ-ous work [26, 31]. To provide simple definitions of the error measures, werewrite the 2 m × f as a m × t = { t x , t y } at every grid node. Real traction and reconstructed traction aredenoted by t true and t recon , respectively. • The Deviation of Traction Magnitude at Adhesions (DTMA) is definedas DTMA = 1 N i X i mean j (cid:0) k t recon j,i k − k t true j,i k (cid:1) mean j (cid:0) k t true j,i k (cid:1) , (12)9here N i is the number of circular traction patches and the index i runs over all patches. The index j runs over all traction vectors in onepatch. The DTMA lies between − • The Deviation of Traction Magnitude in the Background (DTMB) isthe normalized difference between the reconstructed and real tractionmagnitude outside the circular patchesDTMB = mean k ( k t recon k k − k t true k k ) N i P i mean j (cid:0) k t true j,i k (cid:1) , (13)where the index k runs over all traction vectors outside the patches.The DTMB lies in the range [0 ,
1] and a value close to 0 indicates lowbackground noise in the reconstructed traction. • The Signal to Noise Ratio (SNR) is defined in this context asSNR = N i P i mean j ( k t recon j,i k )std k ( t recon k ) . (14)The index k runs over all traction vectors outside the patches while j is the index of each traction vector in the patch i . The SNR measuresthe detectability of a real signal within a noisy background. Its valueranges from 0 to infinity where a SNR that is much larger than unityindicates a good separation between traction and noise. • The Deviation of the traction Maximum at Adhesions (DMA) measureshow peak values of the traction over- or underestimate the true value.The quantity is defined asDMA = 1 N A X i (cid:2) max j ( k t recon j,i k ) − max j ( k t true j,i k ) (cid:3) max j ( k t true j,i k ) , (15)where the maximum traction vector with index j is calculated for eachtraction patch separately. A DMA of 0 indicates that the local trac-tion maxima in the reconstruction and in the original data are equal.Positive or negative values of the DMA imply that the maximum oftraction is overestimated or underestimated.10 .6. Software for traction force calculation We provide a Matlab software package containing the presented Fourier-space methods for calculating traction forces. Note that the program requiresthe input of substrate-deformation data. Usually, substrate deformations arequantified by measuring the lateral displacements of fluorescent marker beadsin a stressed substrate with respect to the marker positions recorded in astress-free state. The standard computational image analysis method for thistask is called particle image velocimetry (PIV) and various well-establishedsoftware packages are available [45, 46, 47]. Once the displacement datahas been extracted, our program can be used to calculate the traction forceswith standard L2 regularization or with Bayesian L2 regularization in Fourierspace. The software is split into a routine for loading data and two routinesfor TFM. The routine “get input data” allows the user to select folders con-taining the data for the measured displacements, the noise, and for images.The required data structure in the file with the displacement data is illus-trated in Fig. 5 (a). Parameters of the experimental setup, including theYoung modulus and the Poisson ratio, also need to be provided. Next, theuser can choose between “Regularization” and “Bayesian regularization”, asshown in Fig. 5 (b) and (c). Selecting “Regularization” allows the choice ofa regularization parameter, which is then held fixed for the whole sequenceof images that are analyzed in the data set. For “Bayesian regularization”,an optimal regularization parameter is selected automatically from the dataset and the noise variance. The standard deviation of the noise can either beprovided as an input or can be determined by manually selecting an image re-gion that is far away from the cell, as illustrated in Fig. 5 (c). Once selected,the region used for determining the noise remains the same throughout thewhole data set of multiple images. After pressing “Analyze sequence” theresults are calculated and saved in automatically named files, see Fig. 5 (c).Since the regularization parameter λ depends in our framework on the noiseand the traction magnitudes, it should be adapted if the signal-to-noise levelchanges significantly. However, note that a change of the parameter withinone image sequence is not always necessary, which reduces the computationaleffort and may be advantageous for data postprocessing.11 . Results To check whether the proposed method actually finds the correct reg-ularization parameter, synthetic data sets with exactly known underlyingdistributions are required. Therefore, we create random traction patterswith traction vectors at each grid point drawn from a Gaussian distribu-tion. Exemplary data is shown in Fig. 2 (a-i). The calculated displace-ment field is then corrupted with a controlled level of noise, see Fig. 2(a-ii). For the reconstruction, we search for the hyperparameter α that max-imizes the log-likelihood function, Eq. (10). As illustrated in Fig. 2 (a-iii),log p ( u | α, β ) has a unique, clear maximum. The regularization parameterdetermined from the optimization compares favorably with the true optimalparameter resulting from the distributions used for simulating the data, hereˆ λ = 9 . × − Pix / Pa ≃ α/β = 10 − Pix / Pa . Visual comparison of thetraction patterns in Figs 2 (a-i,a-iv), as well as a comparison of the tractiondistributions in Fig. 2 (b), confirm that the Bayesian traction reconstructionyields correct results. Note that the measured (posterior) traction distribu-tion does not agree with the original traction distribution when the noisemagnitude is large. This fact is a result of the deviation of the posteriorprobability distribution, Eq. (8), from the prior probability distribution. InFig. 2 (c), we illustrate the difference between the measured traction dis-tribution and the original traction distribution for the synthetic data. Therelative difference of the traction standard deviations is plotted against thevariance of the noise-free displacement field divided by the noise variance, σ /σ . The relative difference of the standard deviation of the measuredposterior and the original traction distribution scales with the relative noisevariance. Figure 2 (d) illustrates how the measurement uncertainty affectsthe mean traction error. For the experimentally relevant regime of measure-ment uncertainties, 0 . & σ noise /σ u & .
1, the relative mean traction error isalmost proportional to the relative measurement uncertainty σ noise /σ u . Forvery low measurement noise, the mean traction error is dominated by nu-merical inaccuracy and aliasing effects. Note that the Bayesian estimate forthe regularization parameters ˆ λ produces errors that are close to the opti-mal errors resulting from regularization with the known parameters α/β forsynthetic data. 12 .2. Quality assessment of traction reconstruction with BFTTC To quantify the reconstruction quality for localized traction patterns, weconstruct synthetic data consisting of circular spots of constant traction asshown in Fig. 3 (a). We employ two classical methods where the regular-ization parameter value is selected by the L-curve criterion, namely a realspace calculation with L2 regularization and regularized Fourier transformtraction cytometry (FTTC). The results are compared with the correspond-ing parameter-free approaches, namely Bayesian L2 regularization (BL2) inthe real-space domain and Bayesian Fourier transform traction cytometry(BFTTC), see Fig. 3 (b). For the real-space TFM results shown exemplarilyin Fig. 3 (c), the L-curve can have a visible corner. Note that the calcu-lations in real space are done with standardized data [31], which rendersthe regularization parameter dimensionless. For the Bayesian real-space ap-proach, illustrated in Fig. 3 (e), the logarithmic evidence always exhibits aclear maximum in our experience. The resulting optimal regularization pa-rameter is usually close to the value from the L-curve criterion. However, inthe Fourier-space approach, illustrated with the example in Fig. 3 (d), the L-curve often does not show a clear corner and it becomes challenging to selectan appropriate regularization parameter. This weakness of the Fourier-spaceapproach is overcome with BFTTC. As illustrated in Fig. 3 (f), the loga-rithmic evidence calculated in BFTTC has a pronounced maximum, whichprovides a clear criterion for the automated choice of the optimal regular-ization parameter. To generate statistics on the performance of the differentmethods, we next record the traction reconstruction quality in 8 separatetests with different traction magnitudes and patterns. The resulting errornorms show that all four methods offer similar traction reconstruction ac-curacies, see Figs. 3(e)(i-iv). The most noticeable reconstruction errors arean underestimation of mean traction (negative DTMA) and a pronouncedtraction background (positive DTMB) [31]. The similarity in reconstructionaccuracy is expected because all methods are based on L2 regularizationand also make use of the same spatial grid for discretization. However, thenumerical effort required for the four methods is very different. Table 1 sum-marizes the computation time required for building the coefficient matrices M or ˜ M and for reconstructing the traction forces. While ˜ M is rapidly builtin Fourier space, the assembly of a large coefficient matrix M in real spacerequires can require many hours. Inferring the optimal regularization pa-rameter requires additional computation time. Overall, real-space methodsare not prohibitively slow but quite impractical for every-day use by experi-13ental scientists. BFTTC, however, requires acceptable computation timesranging from seconds to a few minutes. Table 1:
Computation time for different methods
Reconstruction method L2 BL2 FTTC BFTTCBuilding of M or ˜ M To provide an application example for BFTTC, we quantified the tractionforces generated by NIH 3T3 fibroblasts and mouse podocytes on polyacry-lamide gel substrates. The experiments were done precisely as described inRefs. [34, 44]. The gel substrates had a Young’s modulus of E = 32 kPa anda Poisson’s ratio of roughly ν = 0 .
48. Figure 4(a) shows a cell outline andthe measured displacement data. After recording images of the cell and thenanobeads, the cell was removed from the substrate to provide a stress-freereference for tracking the motion of the nanobeads. We estimate the vari-ance of the noise in the displacement data by quantifying the displacementvariance in a small region that is very far away from the cell and contains nosystematic displacement. Plotting the logarithmic evidence as a function of α yields a curve with a clearly defined maximum, see Fig. 4(b–i,c–i), whichresults in an unambiguous selection of the regularization parameter. Furthertests were performed where, for visualization of the force-generating struc-tures, the cell–substrate adhesions were labeled with GFP-paxillin. It wasfound that BFTTC produces traction maps with defined foci that co-localizewith the GFP-labeled sites of focal adhesion, see the example in Fig. 5.
4. Summary
Traction force microscopy is a popular technique for studying minuteforces generated by biological cells, as well as wetting or frictional forces,on soft substrates. The technique is based on the measurement of substrate14isplacements below the specimen, which allows calculation of the tractionforces. Usually, this calculation is done by solving an inverse linear probleminvolving elastic Greens functions. The procedure requires methods for noisesuppression. Dealing with noise appropriately is an essential issue since thelinear system can be ill-conditioned, which means that the noise can becomeamplified to an extent that the true solution is entirely degraded. A simpleway to remove the effects of noise is to filter the displacement field prior totraction reconstruction. This strategy usually works if the linear problemis solved in Fourier space because the resulting linear system is sparse. Analternative strategy for dealing with noise is regularization, most popular isL2 regularization. With L2 regularization, spatial high-frequency variationsin the data are suppressed, which leads to a robust solution of the inverseproblem of calculating the traction. Regularization is more versatile thandata filtering since it can deal with higher levels of noise, works both inreal-space and Fourier-space approaches, and ensures robust reconstructionif non-standard Green’s functions are employed, for example to take intoaccount three-dimensional substrate topography and tractions. Regardlessof the method, suppression of noise always reduces the spatial resolution.Optimal resolution of the fine details of the traction field can only be gainedif the level of noise suppression is adapted for each sample. For L2 regular-ization, this adaptation is done by changing the regularization parameter,which is usually a manual process based on heuristics, which introduces aconsiderable degree of subjectivity in the resulting traction.Here, we have introduced a Bayesian method for automatic inference of theL2 regularization parameter for traction reconstruction in Fourier space. Us-ing synthetic data of different type, we demonstrate that Bayesian Fouriertransform traction cytometry (BFTTC) is a fast and reliable method. Ourtests show that BFTTC can handle large measurement noise. However, thenoise- and displacement variances ideally satisfy σ /σ u . .
01 for accu-rate traction reconstruction. While the quality of traction reconstructionwith BFTTC is comparable to other methods based on L2 regularization,the choice of the regularization parameter is now automated. Heuristics likethe L-curve criterion, which is particularly ambiguous in Fourier space, areno longer required. The additional computation time required for determin-ing the optimal regularization parameter in BFTTC is only a few secondsto minutes for large data sets. In our experience, the logarithmic evidencealways exhibits a maximum that is sufficiently pronounced to yield a reg-ularization parameter estimate. However, it is important to keep in mind15hat the algorithm is based on the assumption of a Gaussian prior distribu-tion that is symmetric around the origin. Thus, the use of BFTTC is notrecommended if the traction forces in the field of view do not balance eachother. Moreover, if complex, non-Gaussian traction distributions, e.g., multi-modal distributions, are expected, it may be preferable to resort to Bayesianmethods with prior distributions tailored to the specific problem in order tomaximize the reconstruction quality.To provide users from biology, physics, and materials sciences with an easy-to-use tool to analyze their TFM data, we implemented BFTTC as well asregularized FTTC as a Matlab package. The package comes with a user-friendly graphical interface, requires minimal knowledge of the algorithmicdetails, and is freely available [48].
5. Acknowledgments
We thank S. V. Plotnikov (University of Toronto) and C. Schell (Albert-Ludwigs-University Freiburg) for providing the experimental test data.
References [1] William J Polacheck and Christopher S Chen. Measuring cell-generatedforces: a guide to the available tools.
Nat. meth. , 13(5):415, 2016.[2] Ulrich S Schwarz and J´erˆome RD Soin´e. Traction force microscopyon soft elastic substrates: A guide to recent computational advances.
Biochim. Biophys. Acta , 1853(11):3095–3104, 2015.[3] Pere Roca-Cusachs, Vito Conte, and Xavier Trepat. Quantifying forcesin cell biology.
Nature cell biology , 19(7):742–751, 2017.[4] Robert W Style, Rostislav Boltyanskiy, Guy K German, Callen Hyland,Christopher W MacMinn, Aaron F Mertz, Larry A Wilen, Ye Xu, andEric R Dufresne. Traction force microscopy in physics and biology.
Softmatter , 10(23):4047–4055, 2014.[5] Michael J Siedlik, Victor D Varner, and Celeste M Nelson. Pushing,pulling, and squeezing our way to understanding mechanotransduction.
Methods , 94:4–12, 2016. 166] Jeffrey A Mulligan, Fran¸cois Bordeleau, Cynthia A Reinhart-King, andSteven G Adie. Traction force microscopy for noninvasive imaging of cellforces. In
Biomechanics in Oncology , pages 319–349. Springer, 2018.[7] Micah Dembo and Yu-Li Wang. Stresses at the cell-to-substrate interfaceduring locomotion of fibroblasts.
Biophys. J. , 76(4):2307–2316, 1999.[8] James P Butler, Iva Marija Toli´c-Nørrelykke, Ben Fabry, and Jeffrey JFredberg. Traction fields, moments, and strain energy that cells exert ontheir surroundings.
Am. J. Physiol., Cell Physiol. , 282(3):C595–C605,2002.[9] Gwendolen C Reilly and Adam J Engler. Intrinsic extracellular matrixproperties regulate stem cell differentiation.
Journal of biomechanics ,43(1):55–62, 2010.[10] Timo Betz, Daniel Koch, Yun-Bi Lu, Kristian Franze, and Josef A K¨as.Growth cones as soft and weak force generators.
Proceedings of theNational Academy of Sciences , 108(33):13420–13425, 2011.[11] Wesley R Legant, Colin K Choi, Jordan S Miller, Lin Shao, Liang Gao,Eric Betzig, and Christopher S Chen. Multidimensional traction forcemicroscopy reveals out-of-plane rotational moments about focal adhe-sions.
Proc. Natl. Acad. Sci. U.S.A. , 110(3):881–886, 2013.[12] Sangyoon J Han, Kevin M Dean, Austin J Whitewood, Alexia Bachir,Edgar Guttierrez, Alex Groisman, Alan R Horwitz, Benjamin T Goult,and Gaudenz Danuser. Formation of talin-vinculin pre-complexes dic-tates maturation of nascent adhesions by accelerated force transmissionand vinculin recruitment.
BioRxiv , page 735183, 2019.[13] Christian Franck, Stacey A Maskarinec, David A Tirrell, and Gu-ruswami Ravichandran. Three-dimensional traction force microscopy: anew tool for quantifying cell-matrix interactions.
PloS one , 6(3):e17833,2011.[14] Vincent Nier, Shreyansh Jain, Teck Lim Chwee, Shuji Ishihara, BenoitLadoux, and Philippe Marcq, Inference of internal stress in a cell mono-layer.
Biophys. J. , 110(7):1625–1635, 2016.1715] Julia Gerber, Tobias Lendenmann, Hadi Eghlidi, Thomas M Schutz-ius, Dimos Poulikakos, Wetting transitions in droplet drying on softmaterials.
Nat. commun. , 10(1):1-10, 2019.[16] Ye Xu, Wilfried C Engl, Elizabeth R Jerison, Kevin J Wallenstein,Callen Hyland, Larry A Wilen, Eric R Dufresne, Imaging in-plane andnormal stresses near an interface crack using traction force microscopy.
Proc. Nat. Acad. Sci. USA , 107(34):14964–14967, 2010.[17] Joshua Bush, Venkat Maruthamuthu, In situ determination of exertedforces in magnetic pulling cytometry.
AIP advances , 9:035221, 2019.[18] Alexandra S Piotrowski, Victor D Varner, Nikolce Gjorevski, and Ce-leste M Nelson. Three-dimensional traction force microscopy of en-gineered epithelial tissues. In
Tissue Morphogenesis , pages 191–206.Springer, 2015.[19] Aereas Aung, Young N Seo, Shaoying Lu, Yingxiao Wang, Colin Jamora,Juan C del ´Alamo, and Shyni Varghese. 3d traction stresses acti-vate protease-dependent invasion of cancer cells.
Biophysical journal ,107(11):2528–2537, 2014.[20] Julian Steinwachs, Claus Metzner, Kai Skodzek, Nadine Lang, IngoThievessen, Christoph Mark, Stefan M¨unster, Katerina E Aifantis, andBen Fabry. Three-dimensional force microscopy of cells in biopolymernetworks.
Nature methods , 13(2):171, 2016.[21] Xavier Trepat, Michael R Wasserman, Thomas E Angelini, Emil Millet,David A Weitz, James P Butler, and Jeffrey J Fredberg. Physical forcesduring collective cell migration.
Nat. Phys. , 5(6):426, 2009.[22] Agust´ı Brugu´es, Ester Anon, Vito Conte, Jim H Veldhuis, MukundGupta, Julien Colombelli, Jos´e J Mu˜noz, G Wayne Brodland, BenoitLadoux, and Xavier Trepat. Forces driving epithelial wound healing.
Nature physics , 10(9):683, 2014.[23] Benedikt Sabass, Matthias D Koch, Guannan Liu, Howard A Stone, andJoshua W Shaevitz. Force generation by groups of migrating bacteria.
Proc. Natl. Acad. Sci. U.S.A. , 114(28):7266–7271, 2017.1824] Marie-C´ecilia Duvernoy, Thierry Mora, Maxime Ardr´e, Vincent Cro-quette, David Bensimon, Catherine Quilliet, Jean-Marc Ghigo, MartialBalland, Christophe Beloin, Sigol`ene Lecuyer, et al. Asymmetric adhe-sion of rod-shaped bacteria controls microcolony morphogenesis.
Naturecommunications , 9(1):1120, 2018.[25] Janice H Lai, Juan C del Alamo, Javier Rodr´ıguez-Rodr´ıguez, andJuan C Lasheras. The mechanics of the adhesive locomotion of terres-trial gastropods.
Journal of Experimental Biology , 213(22):3920–3933,2010.[26] Benedikt Sabass, Margaret L Gardel, Clare M Waterman, and Ulrich SSchwarz. High resolution traction force microscopy based on experimen-tal and computational advances.
Biophys. J. , 94(1):207–220, 2008.[27] Lev D Landau and EM Lifshitz. Theory of elasticity, vol. 7.
Course ofTheoretical Physics , 3:109, 1986.[28] Zhaochun Yang, Jeen-Shang Lin, Jianxin Chen, and James HC Wang.Determining substrate displacement and cell traction fieldsa new ap-proach.
J. Theo. Biol. , 242(3):607–616, 2006.[29] Xin Tang, Alireza Tofangchi, Sandeep V Anand, and Taher A Saif. Anovel cell traction force microscopy to study multi-cellular system.
PLoSComput. Biol. , 10(6):e1003631, 2014.[30] Sangyoon J Han, Youbean Oak, Alex Groisman, and Gaudenz Danuser.Traction microscopy to identify force modulation in subresolution adhe-sions.
Nat. Methods , 12(7):653–656, 2015.[31] Yunfei Huang, Christoph Schell, Tobias B Huber, Ahmet Nihat S¸im¸sek,Nils Hersch, Rudolf Merkel, Gerhard Gompper, and Benedikt Sabass.Traction force microscopy with optimized regularization and automatedbayesian parameter selection for comparing cells.
Scientific reports ,9(1):539, 2019.[32] Micah Dembo, Tim Oliver, A Ishihara, and K Jacobson. Imaging thetraction stresses exerted by locomoting cells with the elastic substratummethod.
Biophys. J. , 70(4):2008–2022, 1996.1933] Sebastian Houben, Norbert Kirchgeßner, and Rudolf Merkel. Estimatingforce fields of living cells–comparison of several regularization schemescombined with automatic parameter choice. In
Joint Pattern Recogni-tion Symposium , pages 71–80. Springer, 2010.[34] Sergey V Plotnikov, Benedikt Sabass, Ulrich S Schwarz, and Clare MWaterman. High-resolution traction force microscopy.
Methods in cellbiology , 123:367, 2014.[35] Jean-Louis Martiel, Aldo Leal, Laetitia Kurzawa, Martial Balland, IreneWang, Timoth´ee Vignaud, Qingzong Tseng, and Manuel Th´ery. Mea-surement of cell traction forces with imagej. In
Methods in cell biology ,volume 125, pages 269–287. Elsevier, 2015.[36] Claude N Holenstein, Unai Silvan, and Jess G Snedeker. High-resolutiontraction force microscopy on small focal adhesions-improved accuracythrough optimal marker distribution and optical flow tracking.
Sci.Rep. , 7:41633, 2017.[37] Sherry H Suyu, PJ Marshall, MP Hobson, and RD Blandford. Abayesian analysis of regularized source inversions in gravitational lens-ing.
Mon. Not. R. Astron. Soc. , 371(2):983–998, 2006.[38] Abhik Ghosh, L´eon VE Koopmans, Emma Chapman, and Vibor Jeli´c.A bayesian analysis of redshifted 21-cm h i signal and foregrounds: sim-ulations for lofar.
Mon. Not. R. Astron. Soc. , 452(2):1587–1600, 2015.[39] Yong Huang, James L Beck, Stephen Wu, and Hui Li. Robustbayesian compressive sensing for signals in structural health monitor-ing.
Comput.-Aided Civ. Infrastruct. Eng. , 29(3):160–179, 2014.[40] David JC MacKay. Bayesian interpolation.
Neural Comput. , 4(3):415–447, 1992.[41] Per Christian Hansen. Regularization tools version 4.0 for matlab 7.3.
Numerical algorithms , 46(2):189–194, 2007.[42] Curtis R. Vogel Non-convergence of the L-curve regularization parame-ter selection method.
BIT Inverse problems , 12(4):535–547, 1996.2043] Martin Hanke. Limitations of the l-curve method in ill-posed problems.
BIT Numerical Mathematics , 36(2):287–301, 1996.[44] Christoph Schell, Benedikt Sabass, Martin Helmstaedter et al. ARP3controls the podocyte architecture at the kidney filtration barrier.
Dev.Cell u= Mf Gel f Measured displacements Reconstructed traction
Figure 1:
Schematic of traction force microscopy (TFM) to measure cellulartraction on flat elastic substrates.
Adherent cells deform the substrate and the dis-placement field u is obtained by tracking markers within the gel. The traction force field f generated by the cell is calculated by inverting a linear equation system. f =1e4 mean of force=0 Gaussian traction [ P a ] u . [ P i x ] displacement mean of noise=0 Gaussian noise f maximize evidence obtainLog evidence: of regularization parameter [Pix /Pa ] [ P a ] Reconstructed ^ Automatic selection
Log evidenceSimulated
Add noise [Pa ] -2 [Pix ] -2 noise =1e-4 [Pix ] [Pa ] Arti (cid:1) cial (true) traction traction ^ i ii iii [Pix /Pa ] =1e-4 =1e4 ^^ [Pa] cb ^ ^ d Figure 2:
Validation of the Bayesian method for regularization parameterchoice. (a,i) Traction force vectors, discretized on a quadratic mesh, are randomly chosenfrom a Gaussian distribution with fixed variance σ f . Space bar: 100 Pix = 10 grid spacingson a 50 ×
50 mesh. (a,ii) Using the prescribed traction as input, a displacement field iscalculated and Gaussian noise with a variance σ is added. (a-iii) The regularizationparameter is determined by localizing the maximum in the log evidence curve and tractionforces are subsequently calculated. (b) Histogram of the tractions for the sample shownin (a). In the limit of weak noise, the histogram of the reconstructed traction matchesthe true traction distribution. (c) Relative difference between the standard deviation ofthe measured traction distribution σ BFTTC f and the width of true traction distribution σ true f . The grid mesh sizes are denoted by dx. σ u is the variance of the synthetic displace-ment data prior to corruption with noise. Increasing the noise level produces a measured(posterior) traction distribution that no longer agrees with the true traction distribution.(d) Mean error of the reconstructed traction as a function of the relative measurementuncertainty σ noise /σ u . The Bayesian estimate for the regularization parameter ˆ λ and theoptimal regularization parameter α/β produce comparable errors for all noise levels. Arti fi cial traction i [ P a ] [Pix] [Pix ] [Pix ] [Pix ] [ P a ] g L2 L2 L-curve
B(cid:3)(cid:4) FT(cid:5)C(cid:6)(cid:7)(cid:8)(cid:9)(cid:10)
Spatial domain (cid:11)r equency domain
L2 ReconstructionBL2 Log(Evidence) BL2 Reconstruction a bc d
FTTC L-curve FTTC Reconstruction [ P a ] BFTTC Log(Evidence) BFTTC Reconstruction [ P a ] RegularizationRegularizationBayesian Manual selection of Automaticselection of e f [ P a ] [ P a ] i λ Manual selection of λ λ ^ Automaticselection of λ ^ [-] ii [-] [Pix] [ P i x ] iiii iii S i g n a l t o N o i s e R a t i o iviii M a x i m u m a t A d h e s i o n s D e v i a t i o n o f T r a c t i o n ii D e v i a t i o n o f T r a c t i o n M a g n i t u d e a t B a c k g r o un d D e v i a t i o n o f T r a c t i o n i M a g n i t u d e a t A d h e s i o n [Pix ] i no corner detectable Figure 3:
Reconstruction quality of BFTTC compared to other regularizationmethods. (a) Synthetic traction force pattern that is used for testing the reconstruc-tion. Space bar: 5 µ m. (b) Tabulated overview of the compared traction reconstructionmethods. (c) Classical traction reconstruction in real space with L2 regularization. TheL-curve shows a slight “corner”, which is used to determine the value of the regularizationparameter. Note that calculations in real space are done with standardized data [31],which renders the regularization parameter dimensionless. (d) Bayesian L2 regularization(BL2) in real space determines the regularization parameter value automatically. The au-tomatically determined regularization parameter is close to the one predicted in (c) fromthe L-curve. (e) Classical, regularized Fourier transform traction cytometry (FTTC). TheL-curve does not show a “corner”, which makes it difficult to determine an appropriateregularization parameter. (f) Bayesian Fourier transform traction cytometry (BFTTC)determines an optimal regularization parameter automatically. (f) Comparison of thereconstruction quality measures in 8 synthetic data sets; error bars are the standard de-viations of the samples. The reconstruction accuracy of all four methods is found to besimilar. [Pa] [Pix ] (cid:15) m BFTTC (cid:16)(cid:17)(cid:18)(cid:19) Pa [-] (cid:20)(cid:21)(cid:22)(cid:23) Pa Frequency domain (cid:24)(cid:25)(cid:26)(cid:27) [Pa] (cid:28) m (cid:29) m BL2
S(cid:30)(cid:31) !" * m + m bc iiii iiiiii a , m iii - m ./01 3 457 89:;<=>?@AD t Figure 4:
Test of Bayesian Fourier transform traction cytometry (BFTTC) us-ing experimental data. (a,i-ii) Image of an adherent cell and the measured gel dis-placements. Only every 7-th displacement is shown for better visibility. The cell edge isoutlined in white. (b) Results from traction calculation with BFTTC. (b,i) A plot of thelogarithmic evidence reveals a clear maximum, which serves to determine the regulariza-tion parameter. (b,i-ii) Calculated traction forces. (c) Results from traction calculationwith the real-space method BL2 for comparison with BFTTC. While the two methodsproduce similar results, traction fields calculated with BFTTC are slightly smoother thanthe fields calculated with the real-space method due to the different discretizations. nput_data.mat.noise(i) .displacement(i).vec.pos i frame .vec.posResults Manual selection Two methods for reconstruction Automatic selection
Regularization cb Bayesian regularizationGet input data a Manual noise selectionStructure of input data
Figure 5:
Graphical user interface of the provided software for regularizedFTTC and BFTTC. (a) The “get data” interface allows users to input data locationsand parameters of experimental setup. The data structure of the input files can handle awhole video sequence or individual traction recordings. (b) If the “Regularization” optionis chosen, a regularization parameter in units of Pix must be provided by the user. (c) Ifthe option “Bayesian regularization” is chosen, the regularization parameter is automati-cally determined from the measured displacement data and its noise variance. A samplewith displacement noise can either be provided with the input file or it can be determinedfrom a manually selected region that is far away from the cell. A “Preview” button offersthe possibility to visually inspect the solution before one presses “Analyze sequence” tocalculate and save the results.must be provided by the user. (c) Ifthe option “Bayesian regularization” is chosen, the regularization parameter is automati-cally determined from the measured displacement data and its noise variance. A samplewith displacement noise can either be provided with the input file or it can be determinedfrom a manually selected region that is far away from the cell. A “Preview” button offersthe possibility to visually inspect the solution before one presses “Analyze sequence” tocalculate and save the results.