Hindrances to precise recovery of cellular forces in fibrous biopolymer networks

How cells move through the three-dimensional extracellular matrix (ECM) is of increasing interest in attempts to understand important biological processes such as cancer metastasis. Just as in motion on flat surfaces, it is expected that experimental measurements of cell-generated forces will provide valuable information for uncovering the mechanisms of cell migration. However, the recovery of forces in fibrous biopolymer networks may suffer from large errors. Here, within the framework of lattice-based models, we explore possible issues in force recovery by solving the inverse problem: how can one determine the forces cells exert to their surroundings from the deformation of the ECM? Our results indicate that irregular cell traction patterns, the uncertainty of local fiber stiffness, the non-affine nature of ECM deformations and inadequate knowledge of network topology will all prevent the precise force determination. At the end, we discuss possible ways of overcoming these difficulties.


I. INTRODUCTION
The migration of eukaryotic cells in complex environments plays a significant role in many biological processes, such as embryonic morphogenesis, immune defense, and tumor invasion [1]. One widely encountered biomechanical environment for migrating eukaryotic cells in vivo is the three-dimensional (3D) extracellular matrix (ECM), composed of a dense network of biopolymers such as collagen and fibrin [2]. To make their way through ECM, cells apply a variety of different strategies, involving mechanisms of cytoskeleton force generation, protease production and cell adhesions. Whatever strategy cells use, cell-generated forces acting on the ECM are valuable clues to infer what is happening within a migrating 3D cell.
Recently, experimental advances have been made in quantifying the ECM's response to migrating cells [3][4][5][6]. In these experimental setups, cells are often cultured in artificially synthesized extracellular matrix (ECM), such as type I collagen gel, which efficiently mimics the environment in living tissues [7]. As the cells migrate, they deform the surrounding environment; this deformation is trackable by, for example, placing marker beads in the gel or imaging collagen itself. However, from a theoretical perspective, a gap still exists between knowing the deformation of the ECM and determining what forces cells have exerted on that ECM. The inversion from the former to the latter remains elusive, because the ECM displays very complex properties such as strain-stiffening [8] and non-affine deformation [9].
In recent work, Steinwachs et al attempted a reconstruction scheme based on a continuous elasticity model, which phenomenologically captures the strain-stiffening property of collagen gels [6]. This effort goes beyond previous approaches which used linear elastic assumptions and hence represents a significant step forward [3,4,10]. However it remains unclear how accurately this method would capture the mechanics of a real biopolymer network. In particular real networks are expected to exhibit micromechanical fluctuations [11] in its properties on the scale of the network elements which are not very different than that of the embedded cell. Therefore, it is possible that a more detailed understanding of ECM networks is essential for a quantitatively successful reconstruction of cellular forces in ECM.
Here, we use a lattice-based mechanical model to study this force reconstruction problem. In particular, we make use of recent progress in the soft-matter physics community towards an understanding of the ECM system [12]. It has been shown that these systems can be modeled as a disordered network of semiflexible polymers with an interplay between bond stretching and bending. On the basis of this idea, computational models have been built to capture the critical properties such as strain stiffening, negative normal stress, and non-affine deformations. These models roughly fall into two categories: latticebased models [13][14][15][16][17][18][19][20][21] and off-lattice models [21][22][23][24][25][26][27][28][29] . The former places straight fibers on a regular lattice; these fibers are determined by straight segments of bonds on a diluted network. The other approach consists of placing stochastically positioned fibers, intersecting with each other and forming crosslinks; this is usually referred to as a Mikado model. As in both of these cases the mechanics is controlled by critical behavior [13,19,21] around the Maxwell point (at which bending becomes the dominant response mechanism at small strain), the results from these different approaches are extremely consistent with each other [30,31].
In this paper, we study the general reconstruction problem based on a two-dimensional (2D) diluted lattice model. We choose the lattice model for its computational efficiency, The fact that lattice-based models in 2D and 3D capture similar non-affine behavior and nonlinear elastic response of ECM [13], allows us to quantitatively explore the feasibility of the force reconstruction without a need to carry out full 3D simulations. Our approach enables us to study the feasibility of doing this reconstruction even if we do not know the exact microstructure of the material. We find several factors that prevent precise recoveries of cellular forces, including the irregular cell traction pattern, the uncertainty of local fiber stiffness, the non-affine nature of ECM and inadequate knowledge of network topology. Our findings sug-gest continuum theories adopted in most existing traction force microscopy (TFM) methods may generate large inaccuracies in force estimation. At the end of our paper, we discuss how future studies may be able to correctly characterize micromechanical information and give more accurate force predictions.

II. MODELS AND SIMULATIONS
We model the ECM as a diluted triangular lattice (Fig  1a). In this lattice, each bond, with stretching stiffness k and bending stiffness κ , exists stochastically with a probability p. The Hamiltonian is: in which a is the rest length of each bond. g ij = 1 if the bond between node i and j is present and g ij = 0 if the bond is removed. The first term refers to the stretching energy: < i, j > sums over all neighboring lattice sites and |R ij | is the length of bond in the deformed state. The second term represents the bending energy: < i, j, k > sums over all groups of three co-linear consecutive lattice sites in the reference state and ∆θ ijk is the change of angle in the deformed state. In lattice models, p satisfies pZ = z , Z is the coordination number, which is 6 for a triangular lattice, and z is the average connectivity of the biopolymer network. Since experiments have shown z ≈ 3.4 [32], we study p in the range [0.5, 0.65] and use p = 0.57 to illustrate our findings. The same p value was adopted in our previous study [11], and it predicts micromechanical properties consistent with experimental observation.
We insert a round cell into the network by cutting a circular hole in the middle (Fig 1b). The cell and the network are connected by the attached nodes, which are the intersection points between the cell boundary and network bonds in the reference state. Then we let the cell stretch by radially displacing attached nodes; for the i th attached node, its displacement towards the cell center satisfies: in which θ i is the angular position of the i th node, the vector P = (A 0 , A 1 , ...A N , B 0 , ..., B N ) are the determining parameters of cell stretching and the length of P determines the number of degrees of freedom. Note that we use cell contraction to generate motion of the points where the cell is attached to the surface; exerting this displacement is equivalent to putting traction forces on the attached nodes, since the Hamiltonian of our simple model is known. In real biological systems, the Hamiltonian may not be known in detail, which of course causes difficulties in modeling the system. In this paper, we address the other aspect of the problem, namely what effects might prevent the successful inference of the traction forces from marker beads displacements even given that we have a models that accurately describes the system.
Once the cell is deformed, we relax the network into its energy minimum state by the conjugate gradient method. Each parameter setting P leads to a particular cell stretching pattern: which results in ECM deformation. In our scheme, we assume that ECM deformation is measured through the displacements of a set of M marker beads, as following: We stochastically set the parameters P = P set and generate the corresponding "observed" displacements of marker beads, D observed . By minimizing we can realize the force reconstruction by approximating the P set with P guess . To minimize the function f , we apply Particle Swarm Optimization (PSO) [33].

Intrinsic Limitations in Reconstruction
We carry out our simulations in a 60 x 60 lattice, with p ∈ [0.50, 0.65]. To check the possible boundary effects of the lattice, we have verified that a 100 x 100 lattice gives similar results. We set the lattice spacing to be a and our round cell has an initial radius R 0 = 8a. Physically, a corresponds to the persistent length of collagen fibers ( ≈ 1µm). Suppose the cell contraction can be characterized by the N longest wavelength modes: N ), and ξ i are stochastic variables uniformly distributed between 0 and 1. A constant a 0 ( usually around 2a to 3a in our simulations) is added to make sure that the cell mostly contracts, as most cells studied in such ECM systems are contractile. To track the ECM deformation, we initially let all lattice sites be inhabited by marker beads. As mentioned above, we solve for the stretching pattern that leads to the observed bead displacements with Particle Swarm Optimization algorithm.
Since the Particle Swarm Optimization algorithm is a heuristic algorithm, instead of giving a deterministic answer, we measure how likely is it that the reconstruction will succeed. Here, we define successful reconstructions as predictions for all traction forces of all attached nodes deviating by no more than 5% in magnitude and 5 degrees in direction from the input data . Of course, the chances to succeed for a PSO procedure depends on the number of searching rounds and ideally one can always get to the right answer with infinite computer time. Due to limited resources, we assign at most three rounds of PSO searches, each with a maximum of 100 steps.
We ran groups of 50 simulations for a variety of N values. The results indicate the existence of a limit of resolution in force reconstructions for spatial frequencies. The reconstructions are mostly successful if the cell traction input is restricted to the 3 or 4 longest wavelength modes. Conversely, the algorithm experiences a sharp drop in its performance when higher frequencies are present (Fig  2a). That is to say, a highly non-uniform cell contraction can reduce the resolution of both the magnitude and the direction of the reconstructed forces (Fig 2b-g). It may be surprising that the minimization of an error function in as few as 10 dimensions cannot converge to the global minimum, but the problem here is highly nonlinear, i.e. the cell contraction will be transformed into marker beads displacements in a highly nonlinear way. When high-frequency modes are present in the cell contraction, the landscape of the error function f becomes highly rugged (Fig 3), which traps our optimization algorithm in local minima. Sometimes (as is the case in Fig. 2) the reconstruction can still recover information regarding the low-order modes even while failing to accurately determine the higher-order ones. Other times, as for example in Fig. 4, even the low-order modes are estimated incorrectly. Our results indicate, therefore, that only relatively smooth force profiles can be reliably recovered, Next, we consider the effect of errors in the experimental measurements. We stochastically disturb the positions of all marker beads after the network has relaxed to its energy-minimized state by adding Gaussian noise. We then redo the inverse problem. We found a reasonable degree of robustness in the reconstruction, with accurate forces found even with the variance of the added noise as large as 30% bond length (see Fig 5a). Another possible source of error is the inaccurate estimate of the parameters in the Hamiltonian. Accordingly, we vary the stretching and bending stiffnesses of our system and see how it influences the traction force calculation. As shown in Fig 5b, our scheme fails as soon as stretching stiffness varies by more than 10%, indicating the traction force recovery is very sensitive to fiber stiffness measurements.
We have discussed the role of cell contraction pattern and measurement errors in our general scheme of the traction force reconstruction. In these results, we have assumed that we can measure the displacements of all nodes in our network. In other words, we are considering the high marker beads concentration limit, where we have complete information of the network deformation. A straightforward question that arises here is how the bead marker distribution affects the force reconstruction pre- cision. To answer this question, we start by "turning off" the marker beads far away from the cell center ( > R max ). Noise-free simulations with R max = 10a, 15a, 20a, 25a all lead to perfect performance, indicating that the full information of network deformation is redundant. To further investigate the effects of the distribution of marker beads near the cell, we only keep the marker beads within a range of R min to R max from the center of the cell. As Table 1 shows, the reconstruction performance drops sharply as the marker beads are located further than 15 lattice spacings from the cell, which suggests information loss in long-distance force transmission. That is to say, the beads closer to the cell contains more information regarding cellular forces. In addition, inversion attempts with both randomly and regularly diluted (up to 50%) distributions of beads in the neighborhood of the cell show similar performance (data not shown), showing the insensitivity of the reconstruction to details of the local bead distributions. So far, all our studies focus on the contraction of round cells. In reality, however, cells often display long protrusive structures. As a simple extension of our approach, we have also studied elliptical cells (Fig6 (a,b)). We allow an elliptical cell to translate and to shrink along its long axis. In this case, our simulations can quickly converge to the correct reconstruction. (Fig 6(c,d,e)). Noticing that slight tangential movements are also possible, we also verify that the reconstruction works well for a uniformly rotated round cell ( Fig 7(a, b) ).

Roles of Non-affinity and Micromechanics
As we have discussed above, the response of the ECM to cell-induced deformation is highly non-linear and nonaffine. This is ultimately due to the fact that the network lies well below the Maxwell point and hence exhibits a transition from bending-dominated to stretchingdominated as strain increases and the fibers become more aligned. The other significant feature is that the ECM, both in our model and in experiments, is highly heterogeneous on the cellular scale. These features argue against being able to accurately reconstruct cell contractions by spatially-uniform linear elastic models. To investigate the seriousness of the problem, we first compare our results with calculations on a full triangular lattice, with all bonds present except for the ones attached to the cell, representing an approximately affine model. To make a fair comparison, we have to relate the parameters of this system with our diluted lattice model. Both the bulk modulus and shear modulus of the full lattice are solely determined by stretching stiffness k, since the deformation of the full lattice is affine and the bending energy is always zero. We tune the value of k to match the bulk modulus of the full lattice with that of the dilute lattice.
(Alternatively, we can match the shear modulus, which leads to similar conclusions.) In a forward comparison, we contract our cell with exactly the same pattern in a diluted lattice and in a full lattice and compare the displacements of marker beads in the range between ten and fifteen lattice spacings from the center of the cell. As seen in Fig 8(a-c), significant differences in the network responses are observed between the non-affine case (p ∼ 0.6) and the affine case (p = 1).
In the reconstruction comparison, we relax our system on a diluted lattice so as to locate the marker beads and then attempt to recover the cell deformation based on a full lattice (Fig 9(a,b) ); this would be analogous to doing the inverse problem with a linear elastic model. In all of our 50 attempts, reconstructions fail with large deviations (Fig 9 (c,d,e) ). Both the forward and backward calculations imply that not surprisingly, the non-affinity of ECM is critical in the force transmission. Simple affine models mischaracterize the system and lead to highly inaccurate predictions.
Finally, we turn to the role of micromechanics, namely the local variation of elastic response properties that of course correlate with the precise local network structure. Even if an ECM model captures the non-affinity aspect properly, the precise network topology can induce notable differences in the network deformation. To further elucidate the idea, we apply similar forward and backward comparison between two networks of the same p value but different network topology. In the forward comparison, deformations may differ by a magnitude about twice the lattice spacing in these two cases for the same cell contractions (Fig 10a). As for the backward comparison, we relax the network within one network, while reconstructing cell contraction based on a differently connected network with the same macrscopic properties. All simulations, without exception, lead to remarkably inaccurate predictions in both force magnitude and direction (Fig 10(b,c) ). An immediate consequence is that a nonlinear continuum model would be incapable of correctly describing the deformation field of cellularized networks, since network topology detail is missing in any continuum model. Therefore, we need much more information about network structure to successfully reconstruct traction forces.

DISCUSSION
In this study, we explore how the mechanical complexity of the ECM material limits our ability to invert ECM deformations to obtain traction forces. To be more specific, we have explored the feasibility of reconstructing cell-generated traction forces from ECM deformation based on a general discrete ECM model. Our model is based on a 2D diluted triangular lattice, which has been shown to mimic both macroscopic properties of rheology such as strain-stiffening and alignment, and local heterogeneities on the scale of individual cells of real 3D ECMs [11]. Moreover, the lattice model in this study is able to capture the buckling effects of single fibers, as discussed in Ref [31]. We limit our investigation in 2D in this study for the sake of computational efficiency. The extension to a 3D model is straightforward and it should exhibit similar physics.
Within our model, we find that there is an intrinsic limitation in the spatial resolution of cell contraction. This limitation is not too surprising. We can imagine a pair of closely neighboring adhesion sites on the cell membrane which undergo very different amounts of stretching; if they exchange their stretching, almost the same displace- ments of marker s will be observed, thus it is hard to distinguish the different configurations in the reconstruction. It explains why a relatively smooth stretching configuration, where closely neighboring bonds stretch similarly, will suffer much less from this degeneracy. In addition, we find that the inaccurate estimates of stretching stiffness or bending stiffness of fibers can lead to imprecise predictions of traction forces. Note that we assume constant fiber stiffness in our network; In reality, however, each fiber has different thickness and length, and fiber stiffness varies from fiber to fiber, which further hinders the accuracy of traction force reconstruction.
By comparing affine networks (p ∼ 0.6) and non-affine networks using a forward and a backward approach, we find that the non-affinity of networks can fundamentally change the way cell deforms ECMs. Furthermore, even two networks with a same p value but different topology can exhibit remarkable differences in both forward and backward comparisons. These results strongly suggest that continuum models, which leave out microscopic structures and heterogeneity of ECMs, are not appropriate tools in inverting traction forces from ECM deformations. It also suggests that the inversion must make use of direct measurements of local network structure as op-  To overcome all the hindrances to precise recoveries of traction forces, we suggest that the only feasible option is to couple local active local measurements to the inversion algorithm. These measurements can be implemented by actively perturbing the beads located in a variety of regions of the ECM, as was done for example in a recent study of ECM micromechanics [11]. With observations of the mechanical responses for each of these perturbations observed, we believe that the micromechanical information can be sampled well enough to reconstruct a network model for that specific piece of ECM and thereby recover precisely the forces the cells are exerting on that ECM. We will leave an attempt to verify this approach to future work.