Avoiding tensional equilibrium in cells migrating on a matrix with cell-scale stiffness-heterogeneity

Intracellular stresses affect various cell functions, including proliferation, differentiation and movement, which are dynamically modulated in migrating cells through continuous cell-shaping and remodeling of the cytoskeletal architecture induced by spatiotemporal interactions with extracellular matrix stiffness. When cells migrate on a matrix with cell-scale stiffness-heterogeneity, which is a common situation in living tissues, what intracellular stress dynamics (ISD) emerge? In this study, to explore this issue, finite element method-based traction force microscopy was applied to cells migrating on microelastically patterned gels. Two model systems of microelastically patterned gels (stiff/soft stripe and stiff triangular patterns) were designed to characterize the effects of a spatial constraint on cell-shaping and of the presence of different types of cues to induce competing cellular taxis (usual and reverse durotaxis) on the ISD, respectively. As the main result, the prolonged fluctuation of traction stress on a whole-cell scale was markedly enhanced on single cell-size triangular stiff patterns compared with homogeneous gels. Such ISD enhancement was found to be derived from the interplay between the nomadic migration of cells to regions with different degrees of stiffness and domain shape-dependent traction force dynamics, which should be an essential factor for keeping cells far from tensional equilibrium.


Introduction
Changes in intracellular stress brought about by cell deformation and shaping have been shown to affect various cell functions. For example, cell contractility sensitively regulates DNA synthesis/cell proliferation 1,2,3,4 , lineage-specific gene expression/cell differentiation 5,6,7,8 , and the generation of cell-polarity/directed cell movement 9,10,11 .
The strength and spatial distribution of intracellular stresses are determined by the interaction between mechanical properties of the extracellular matrix (ECM) and intrinsic cellular prestress which the cell actively generates inside itself as a result of the activities of adhesion machinery and cytoskeletal architecture 12,13,14,15 . In principle, this interaction is dynamic because of continuous cell-shaping and remodeling of the cytoskeletal architecture during cell movement. Since cells in vivo generally migrate in tissues with characteristic mechanical properties under a microscopic spatial distribution, intracellular stress dynamics (ISD) provide fingerprint-like information on the mechanical interactions between cells and ECM. Especially when the surrounding mechanical field is heterogeneous comparable to a cell-size scale 16,17,18,19 , ISD in an entire cell should exhibit temporal fluctuations that are markedly different from those observed on a homogeneous mechanical field. Although this situation is rather common in vivo, the detailed characteristics of ISD on a matrix with cell-scale heterogeneity remain unclear.
To better understand the actual nature of ISD, traction force microscopy (TFM) is a powerful tool for determining the traction behavior at cell adhesion interfaces which reflect cellular contractility and intracellular stresses. Various types of TFM have been proposed for 2D and 3D analyses 20,21,22,23 . When a cell adheres to a flat homogeneous gel surface, analytical solutions of linear elastic theory, the Boussinesq solution, are often used to reconstruct the traction force 20,21 . However, for a matrix with cell-scale stiffnessheterogeneity, such as a microelastically patterned substrate, the elastic modulus sharply changes within the adhered area of a single cell, and an analytical solution for the equation of equilibrium cannot be obtained in general. Thus, for micromechanically-heterogeneous systems, TFM based on a numerical calculation such as with the finite element method (FEM) must be applied 22,24 . In addition, the time-scale for analysis is essential for investigating the effect of cell-scale stiffness-heterogeneity on ISD. In particular, the long-period fluctuation of traction force in an entire cell is important since it can cause field size-dependent periodic perturbation of the tensional steady-state of the cells. To analyze this phenomenon, TFM should be performed for more than several tens of hours.
In this study, to reveal the feature of ISD affected by cell-scale stiffness-heterogeneity, we sought to determine the long-term and whole cell-scale dynamics of traction force for cells migrating on a heterogeneous matrix. By employing FEM-based TFM 22,24 , we measured the long-term dynamics of traction force of mesenchymal stem cells (MSCs) cultured on microelastically patterned hydrogels 25,26 . Several tens of microns-scaled patterns were applied with stiff/soft stripes and periodic triangular stiff domains. In contrast to previous typical TFM results on a single linear elasticity gradient with pillar arrays 27,28 , the size and curvature of the stiff domain were found to crucially affect the magnitude of the traction stress. We clarified that cell-scale variability of stiffness in an adhesion area determines the magnitude of the traction force as well as the local stiffness of the substrate. Reflecting these behaviors, the long-term computation of traction force dynamics revealed that cell-scale stiffness-heterogeneity strongly enhances the strength and frequency of fluctuation of the ISD in migrating MSCs in comparison to those on homogeneous gels.

Preparation of photocurable sol solution
Photocurable styrenated gelatin (StG) was used for photolithographic microelasticity patterning of the gel, which was prepared as described previously 25,26,29  A soft base gel was prepared by irradiation of the entire sample with visible light for 120s -210s (45 -50 mW / cm 2 at 488 nm; light source: MME-250; Moritex Saitama, Japan). Next, stiff regions were prepared by local irradiation of the base gel using a set irradiation pattern and a custom-built, mask-free, reduction-projection-type photolithographic system. Equilateral-triangular irradiated regions with side lengths L of 150 µm and 250 µm were projected on the gel using an EB-1770W liquid crystal display projector (SEIKO EPSON, Nagano, Japan). The center of the triangle was periodically placed on the honeycomb lattice. The distance W between the corners of the triangles was 120 µm (Fig. 1 (inset)). Finally, the gels were detached from the NIPAAm-coated normal glass substrate and washed thoroughly with PBS at 28 °C to remove the adsorbed PNIPAAm and Tween20. Stripe patterns with a unit length L of 100 and 600 µm ( Fig. 1 (inset)) and homogeneous gels were also prepared.

Measurement of the surface elasticity distribution around the elasticity boundary
The surface elasticity of the StG gel was determined by nano-indentation analysis using an atomic force microscope (JPK NanoWizard 4, JPK Instruments, Bruker Nano GmbH, Germany). A commercial silicone-nitride cantilever with a nominal spring constant of 0.03 -0.09 N/m was used (qp-BioAC-CI CB3, Nanosensors). Young's moduli of the surface were evaluated from force-indentation curves by nonlinear least-square fitting to the Hertz model in the case of a parabolic indenter (tip radius of curvature 30 nm; Poisson ratio µ: 0.5). To make a spatial distribution of elasticity in the triangular pattern gel, we measured about 100 tiles of a 40 x 40 µm elasticity map with a distance of 2 µm between each measured point. Next, we reconstructed the spatial distribution of elasticity from the tiles.
The topographic condition of the sample surface was also measured by using the atomic force microscope. Although the soft regions were 15 -25 µm higher than the stiff regions, the connection between the two regions was smooth enough to not disturb natural durotaxis. In addition, we did not observe the accumulation of MSCs at the basin of the topography of triangular-pattern gel. Thus, in our experiment, curvotaxis is expected to have little effect 30 .

Time-lapse observation of cell movement
The migratory motion of cells on the microelastically patterned gel was monitored using an automated all-in-one microscope with a temperature-and humidity-controlled cell

Analysis of cell-shaping dynamics
Movement trajectories and the shapes of the cells were determined and analyzed using MATLAB software. Based on edge-detection of a cell, we extracted the shape of the cell from phase-contrast images. The details of image-processing have been explained previously 32 . We traced each cell and measured the time evolution of the trajectory and shape. If the cells collided with each other or if a cell was divided, we stopped the trace.
When the cells separated again, we renumbered the cells and restarted the trace. Through the image analysis, cell trajectories x(t) were calculated.

Traction force microscopy
To calculate the traction force of the cells, we used finite-element-method (FEM) based traction force microscopy 22,24 . The displacement field of fluorescent beads was measured by comparing images with and without cells. The displacement of the beads was calculated using commercial PIV software (Flownizer 2D; Detect Corporation, Tokyo, Japan). In the calculation of FEM, a hexahedron was used for discretization (the size of the unit is about 8 x 8 x 8 µm). The nodes of a hexahedron at the substrate surface were chosen to be identical to the nodes of PIV. We assume that the patterned gel is a linear elastic body with spatial modulation according to Young's modulus E(x): where i, j, and k are x, y, z. σik and uik are the stress tensor and strain tensor, respectively. ν is the Poisson ratio.
where Kijkl is a stiffness matrix, Uik is displacement of the k component at node i, and Fjl is force in the l direction applied at node j. From Eq.
where T denotes the transposition. N is the number of the node where the force is reconstructed. J is the number of the node where the displacement is measured experimentally. In the case that the forces are distributed on surface nodes, displacement of the surface U can be expressed by using M: U MF = By solving the inverse problem of Eq. (7), we can reconstruct the distribution of the force.
The flow chart of FEM-based TFM is shown in Fig. S1. Because of the ill posed nature of force reconstruction, some filtering schemes are required to reduce the noise. In this paper, F is reconstructed by minimizing function S, which is a combination of the leastsquare displacement error and a weighted norm of the forces: where Ur is displacement of the substrate measured experimentally. The matrix A is a diagonal matrix containing the local penalty weights. Next, we calculate F, which gives the minimum of S, dS/dF = 0: We then iteratively change A depending on the estimated force. We use a local penalty matrix that was proposed previously 33 .
We repeatedly calculate Eqs. (9) and (10) until convergence is reached. When the regulation parameter α is sufficiently small compared to a typical value of M T M, the reconstructed traction field becomes very noisy. If α is too large, traction stress tends to be underestimated. Since M T M is inversely proportional to the square of the elasticity of the substrate, we vary α depending on the magnitude of M T M. We set regulation parameter α as The validity and accuracy of FEM-based TFM are discussed in the SI (Figs. S2 and S3).

Calculation of the force spot
To identify the force spot, we calculated the threshold of the traction stress. As reported previously 34 , the cumulative distribution of the traction stress around the cell had exponential tails (Fig. S4): where the P0 is the typical traction stress of the cell. Since P0 does not depend on segmentation of the cell shape, we used P0 as the threshold for the force spot.

Definitions of mean traction stress and mean elasticity
The mean traction stress Trm(t) of a cell and mean elasticity Em(t) just below a cell are defined as

Tr t Tr x y t dxdy dxdy E t E x y dxdy dxdy
Ω Ω Ω Ω = = ∫ ∫ ∫ ∫ (13) where Tr(x,y,t) and E(x,y) are the traction stress and elasticity at point (x,y). Ω(t) is the area of the cell that was calculated from the segmentation of the phase-contrast image.

Definitions of magnitude of fluctuation of traction stress and other variables
Similar to the approach described in the literature 35  where σ(Trm) and

Microelastically patterned substrates with cell-scale stiffness-heterogeneity
To investigate the effect of cell-scale stiffness-heterogeneity on cell migration and ISD, we prepared two types of microelastically patterned substrates by photolithography with photocurable styrenated gelatin (StG) (Fig. 1) Table 1.
For the striped gels, the unit widths of 100 µm and 600 µm were designed to model single cell-size and multicellular-size of MSCs, respectively. Young's modulus was ca.
10 kPa for the soft regions and ca. 40 kPa for stiff regions over a transition region of ~50 µm in width (Figs. 1e and f). A stiffness gradient of 30 kPa / 50 µm with a 10 kPa soft region has been previously clarified to have differential effects on the induction of durotaxis between these two unit widths, i.e., durotaxis is induced in the 100 µm-wide stripes, but not in the 600 µm-wide stripes depending on the cell-scale stiffnessheterogeneity 29 . To better understand the reason for these differential behaviors in durotaxis from the perspective of a traction force analysis, the present stiffness gradient was adopted.
For the triangular patterned gels, the side lengths of a stiff equilateral triangle L were designed to be 150 µm and 250 µm to model single cell-size and multiple cell-size, and the unit sizes of the as-prepared soft triangular base including stiff triangles were 254 µm and 354 µm, respectively. (Figs. 1c, d, and g). We changed L with a fixed width of soft regions W (see Fig. 1c and d inset). Young's modulus varied from around 4 kPa to 20 kPa over a transition region of ca. 80 µm in width (Figs. 1h and i). The corner of a stiff triangle has a convex elasticity boundary toward the soft region with a curvature radius of less than 50 µm, where reverse durotaxis from stiff to soft regions can be induced, as we previously reported 26 . Thus, the triangular stiff domains are expected to induce both usual soft-to-stiff durotaxis 36 and reverse stiff-to-soft durotaxis at the side and at the corner of the triangle, respectively. Such a system can be applied as a model of the microscopic mechanical environment of ECM in which different types of elasticity boundary in terms of curvature closely coexist. In living tissue, the presence of curved architecture and the coexistence of stiffness gradients with different nonlinearities are rather common 19,37 .
What features of ISD emerge in a system with different types of cues of competing cellular taxis? To investigate this issue, triangular-patterned gels with a similar range of stiffness gradient as with striped-patterned gels were used for comparisons.

Stripe width-dependent modulation of traction force
We previously reported that, generally, the stiffness gradient required to induce durotaxis is altered by the size of the stiff/soft regions 29   analyzed for the differential responses of durotaxis on different stripe widths with the same stiffness gradient.
As shown in Fig. 2a 29 . Even though the strength of a stiffness gradient is set to be the same for these two conditions, only Stripe-100 induced durotaxis affected by the size of the stiff region. The mechanobiological basis for this differential induction of durotaxis was explored by using FEM-based TFM. Figures 2b and c show a representative deformation field and reconstructed traction field generated by a MSC on Stripe-100. The background black/white shading represents the elasticity field of the gel. As the shading becomes dark, the gel becomes stiffer. Although the gel was highly deformed on soft regions (Fig. 2b), the magnitude of the traction force tended to be large on stiff regions (Fig. 2c). On the other hand, Figs. 2d and e show these fields on Stripe-600. Compared with the results for Stripe-100, large traction forces were exerted on both soft and stiff regions (Fig. 2e). These observations suggested that the striped pattern with a single-cell size could enhance the magnitude of traction force on a stiff domain.
This effect of stripe width on the traction force seemed to dominantly arise at the force spots where large traction stress was localized 34 . Therefore, we next analyzed the dependence of the force magnitude in the force spot on the spatial heterogeneity of stiffness. From the cumulative distribution of the traction stress of an entire cell, we calculated the threshold of traction stress to identify the force spots ( Fig. 2f and Method).
The left side of Fig. 2f shows an example of the spatial distribution of traction stress of a representative cell. Yellow domains on the right side of A narrow stiff region enhanced the strength of traction stress enough to induce durotaxis, while a wide stiff region did not, which was found to be the cause of the emergence of differential responses of durotaxis for the same stiffness gradient.

Traction forces around the triangular stiff domains
In the above case of mechanically striped-patterned gels, only the linear stiffness gradient across an elasticity boundary existed on the gels, and the substrate mechanicsdependent taxis cue could induce one-way directed cell movement, i.e., soft-to-stiff durotaxis. However, stiffness-heterogeneity in a natural living body should not exhibit such a simple situation, and rather includes a variety of nonlinear stiffness gradients, such as a curved elasticity boundary. As mentioned above, microscopic stiffness-heterogeneity of ECM may provide different effects for cellular taxis which may compete with each other. To model this situation, a triangular stiff domain can be applied because it has dual Considering the spatial symmetry, the smallest stiff unit was surrounded by white dashed lines. We calculated the probability distribution on the smallest unit. Next, we replicated the smallest unit to construct the unit triangle. Thus, the probability distribution is highly symmetric. (e) Deformation fields measured on Triangle-150. (f) Distribution of the traction forces estimated from Fig. 4 (e). (e, f) A vector field of different time points was superimposed on a figure. A longer vector represents greater deformation and force. As the color of the vector changes from blue to red, deformation and force increase. (g, h) Color map of the traction stress of the force spot as a function of the position x and y. The Mann-Whitney U test was used to calculate the P value. *** P < 0.001. n =46 for Triangle-150. n =23 for Triangle-250. tactic cues of usual durotaxis 36 and reverse durotaxis 26 . What behaviors of dynamics in ISD appear in the system? To investigate this issue, we analyzed traction forces on triangular patterned gels.
First, we measured cell trajectories to characterize migration (Figs 3a and b). Here, the cell trajectory is defined as the path of the geometric center of the cell shape projected on an x-y plane, which was superimposed on the elasticity distribution estimated from AFM measurement; background black/white shading denotes stiff and soft regions, respectively.  (Fig. 4).
For patterned gels, the time series of mean traction stress was found to largely fluctuate during cell movement between soft and stiff regions (Figs. 4a and b, blue lines) compared to the control (Fig. 4c, blue lines). This tendency was confirmed for 10 cells in each condition (Figs. 4d-f), suggesting that fluctuation of the mean traction stress was greatest on triangular gels and least on the homogeneous gel. In the case of striped gel (Fig. 4d), they seemed to be intermediate. To characterize the typical time-scale of the observed fluctuations of traction stress, we calculated the population-averaged autocorrelation function of the time series of traction stress (Fig. 5a). The autocorrelation functions exponentially decayed to zero within 2 -4 hours without a clear peak, indicating that the fluctuation of the traction force was rather noisy. Through fitting of the autocorrelation function by an exponential exp(-κt), the correlation time Tf = ln2/κ for each cell on a triangular patterned gel was found to be about 50 % smaller than that on homogeneous gels (Fig. 5 b). This result quantitatively confirmed that the temporal fluctuation of traction stress is strongly enhanced on triangular patterned gel.
Regarding the magnitude of the temporal fluctuation of traction stress, the standard deviation of the time series of traction stress divided by the average traction stress of each cell (see Method) was about two times larger on triangular patterned gels than on homogeneous gels (Fig. 5 c). The fluctuation on Stripe-100 was intermediate between those on triangular and homogeneous gels. No significant difference was noted between Stripe-600 and homogeneous gel. Similar to the fluctuation of the traction stress, the magnitude of the deformation strongly fluctuated on the triangular pattern gels (Fig. 5 d).

Discussion
In this study, we focused on the features of intracellular stress dynamics (ISD) in cells migrating on a matrix with cell-scale stiffness-heterogeneity. To address this issue, two model systems of microelastically patterned gels were used: stiff/soft stripe and stiff triangular patterns. These stripe and triangular patterns were designed to characterize the effects of a spatial constraint for cell-shaping and of different coexistent cues to induce competing cellular taxis (usual and reverse durotaxis) on the ISD, respectively. One of the most striking ISD responses in migrating cells was that the strength of long-term fluctuation of traction force on a whole-cell scale was markedly enhanced in the latter system of triangular patterns compared with a stiffness-homogeneous matrix, which made the cells far from tensional equilibrium. Based on this finding, we discuss the following Does nomadic cell migration on a heterogeneous mechanical field exist in a living organism? Some types of cells should migrate on micromechanically heterogeneous milieu, such as typically seen in wound healing of disordered tissues 50,51 , in the biological development of immature organs 52,53,54 , and in stem cell localization inside a specific niche 55,56 . In particular, MSCs are delocalized in niches composed of hard bone, soft vessels, and very soft blood cell surroundings in bone marrow, and are nomadically delocalized in the niche with such stiffness-heterogeneity 56 . This in vivo situation regarding nomadic cell migration, as inferred from the model system in the present study, should also induce the avoidance of tensional equilibrium.

Conclusion
The  . For FTTC, we adopted a Gaussian filter with a cut-off frequency of 0.066 µm -1 for the deformation field. The two methods give a very similar estimated traction field. However, for FEM-TFM, a large traction force tends to be localized in a narrow region. This phenomenon could be due to the regulation method we applied. Figure S2 (c) shows the time course of the mean traction stress under a cell. Blue and red lines represent the traction stress estimated by FEM-TFM and FTTC, respectively. The figure shows complete agreement between these two methods. For FTTC, the mean traction stress increases as the cut-off frequency of the Gaussian filter increases.
Next, through a numerical calculation, we checked the accuracy of the estimation method 1 of the traction force on the elastically patterned substrate. First, we provided an artificial force distribution that mimicked traction forces exerted by a polarized MSC (Fig.   S3 (a)). By using FEM, we calculated the deformation field at the surface of the substrate. Next, we added Gaussian noise to the deformation field with a standard deviation of 0.1 µm, which was the typical value of the background noise observed in the experiment ( Fig.  S3 (b)). Finally, we reconstructed the force distribution from the disturbed deformation field. As shown in Fig. S3 (c), we successfully estimated the force distribution. Compared to the forces in the soft region, estimated forces in the stiff region tended to be noisy, because the signal-to-noise ratio was lower in the stiff region. Plain-stiff. Cells on Plain-soft had a rather circular shape, while cells on Plain-stiff had an elongated shape. As shown in Fig. S5 (c), traction stress in the force spot on Plain-stiff was significantly greater than that on Plain-soft.

Residence time in soft and stiff regions.
We calculated the residence time, which is the amount of time a cell centroid remained in the soft and stiff regions (Fig. S6). Here, we defined the irradiated region of the photomask as a stiff region and the masked region as a soft region (Fig. 1 (c) and (d) inset). Residence times in soft and stiff regions were almost equal for wide stripes and a small triangular pattern. Otherwise, MSCs stayed in stiff regions longer. Compared to the residence time on stripe patterns, the residence time on triangular patterns was short, which indicates that the cell frequently moved between soft and stiff regions.