Persistent random deformation model of cells crawling on a gel surface

In general, cells move on a substrate through extension and contraction of the cell body. Though cell movement should be explained by taking into account the effect of such shape fluctuations, past approaches to formulate cell-crawling have not sufficiently quantified the relationship between cell movement (velocity and trajectory) and shape fluctuations based on experimental data regarding actual shaping dynamics. To clarify this relationship, we experimentally characterized cell-crawling in terms of shape fluctuations, especially extension and contraction, by using an elasticity-tunable gel substrate to modulate cell shape. As a result, an amoeboid swimmer-like relation was found to arise between the cell velocity and cell-shape dynamics. To formulate this experimentally-obtained relationship between cell movement and shaping dynamics, we established a persistent random deformation (PRD) model based on equations of a deformable self-propelled particle adopting an amoeboid swimmer-like velocity-shape relationship. The PRD model successfully explains the statistical properties of velocity, trajectory and shaping dynamics of the cells including back-and-forth motion, because the velocity equation exhibits time-reverse symmetry, which is essentially different from previous models. We discuss the possible application of this model to classify the phenotype of cell migration based on the characteristic relation between movement and shaping dynamics.

Cell migration plays important roles in various physiological and pathological processes in living organisms such as embryogenesis, morphogenesis, immunological response 1 , wound healing 2 , cancer metastasis 3 , etc. The ability to characterize and predict the migration behaviors of various kinds of cells is an important issue not only from a biomedical viewpoint, but also from the perspective of basic science in molecular cell biology. In general, cells dynamically change their shape as a result of contraction by actomyosin and extension through protrusion of the plasma membrane driven by actin polymerization 4 . In a time-scale of from minutes to hours, an entire cell can move based on the sum of such local fluctuations in shape. For example, in the case of keratocytes, extension of the front part and retraction of the rear part occur simultaneously at a constant speed. As a result, the cell experiences ballistic motion with a constant shape 5 . In the case of Dictyostelium cells, local extension and contraction fluctuate spatiotemporally 6 . As a result, cell movement consists of an alternating series of directed motion and random turning, which is called persistent random motion 7 .
With regard to such persistent random motion, random walk-based models, such as the persistent random walk (PRW) model, have been proposed to reproduce the migration patterns, but only if the trajectory does not have strong spatiotemporal correlations [8][9][10][11][12][13] . However, the PRW model does not adequately explain ordered patterns of migration, such as rotation, oscillation, and zig-zag trajectories, because this model assumes Brownian motion. These ordered motions have been reported to derive from the spatiotemporal dynamics of pseudopodia 6,14-17 , i.e., cell-shape dynamics. Thus, to explain spatiotemporally correlated motion, we should consider the effect of the shaping dynamics. However, previous approaches to formulate cell-crawling have not adequately quantified the relationship between cell movement and shape fluctuations based on experimental data regarding actual shaping dynamics. Recently, as a model for the migration of keratocytes and Dictyostelium cells, a phenomenological cell-crawling model was proposed based on the assumption that cell velocity is determined by the cell shape 18 . However, such a shape-based formulation of movement has not been experimentally examined for persistent random motion.
In this study, we aimed to elucidate and formulate the relationship between movement and shape fluctuations through the quantitative analysis of cell-shaping dynamics. First, to clarify the quantitative relationship between velocity and shape, we experimentally characterized the crawling of fibroblast cells in terms of shape fluctuations, especially extension and contraction, by using an elasticity-tunable gel substrate to modulate cell shape. Through a Fourier-mode analysis of the shape, the cell velocity was found to follow the cell-shape dynamics, where the obtained velocity-shape relationship was equivalent to that of an amoeboid swimmer 19 . Next, to formulate such shape fluctuation-based cell movement, we proposed a persistent random deformation (PRD) model by incorporating the amoeboid swimmer-like velocity equation 19 into model equations for a deformable self-propelled particle 18 . The PRD model fully explains the statistics and dynamics of not only movement but also cell shape, including the characteristic back-and-forth motion of fibroblasts. This reciprocating motion is due to the time-reverse symmetry of the amoeboid swimmer-like velocity equation 19 , which is essentially different from previous migration models. Through fitting of experimental data with the model, we quantitatively evaluated fitting parameters, such as mobility, relaxation time of shaping, and magnitude of the internal force. The dependence of the fitting parameters on elasticity revealed that cells showed strong adhesion and large internal force on stiffer gels, as previously reported 20 . Finally, we discuss the possible application of this PRD model to classify the phenotype of the migration of different kinds of cells based on their characteristic relations between movement and shaping dynamics.

Results
Movement and deformation of cells on a gel surface. First, to elucidate the phenomenological relationship between cell movement and deformation, we studied the movement and deformation of NIH 3T3-fibroblast cells on a hydrogel surface with different degrees of stiffness. Photocurable styrenated gelatin (StG) was used because the elasticity of the gel can be adjusted from 1 to 1000 kPa by changing the duration of light irradiation in a photocrosslinking procedure 21 (see Method). To determine the cell response to the elasticity of the substrate, we used 35, 120, and 410 kPa gels. Because cell movement is restricted to the flat gel surface, we analyzed the shape and trajectory of the geometric center of the cell shape projected on an x-y plane. In this paper, we define cell movement as the translocation of a cell body from one site to another. The word "cell migration" is used when we focus on the properties of movement at a cell-population level. Figure 1A,B show phase-contrast images of cell trajectories on 35 and 410 kPa gels, respectively. As the elasticity of the gel increases, the motility and persistence of motion decrease, and the cell body extends 20,22,23 . In all cases, the cells tend to show back-and-forth motion (Fig. 1C). Figure 1D shows time evolutions of movement and cell shape for cells on 35 kPa gel: Fig. 1D1-D3 show slow, middle, and fast cells, respectively. In all cases, the cells migrate along the long-axis of elongation 24,25 . When the cells change their direction of movement, they contract pseudopodia and then extend new ones. Thus, at the turning point of cell movement, the cells change their shape. As long as the cells remain elongated, they persistently migrate. A comparison of the slow and fast cells (Fig. 1D1,D3) shows that the fast cell repeatedly extends and contracts during migration. On the other hand, cells that maintain constant elongation do not move fast (Fig. 1D1). Figure 1D implies that cell movement is correlated with the cell-shape dynamics. Here, we seek to quantify the relation between cell movement and shape. To evaluate cell shape quantitatively, we calculated the complex Fourier coefficient C n of cell-shape R (see Method). The inset in Fig. 1D shows a phase-contrast image of a cell and its periphery detected by image analysis. The distance R(θ) from the centroid to the rim is calculated as a function of the angle θ, where θ is measured from the x axis (Fig. 1D, Inset). We averaged R over θ ± 6° for smoothing. Next, complex Fourier coefficient C n (t) is calculated from R(θ, t). Here, we focus on the elongation mode C 2 (t) and triangular mode C 3 (t). Figure 1E shows the time series of the magnitude of cell velocity |v(t)|, elongation |C 2 (t)|, time derivative of elongation C 2  | |, and time derivative of triangular deformation  | | C 3 . |v(t)| tends to be large when elongation |C 2 (t)| varies rapidly. Thus, the time series of C 2  | | closely resembles that of |v|. |v| is also moderately correlated with C 3  | |. This result indicates that a cell moves when the cell body extends or contracts. Figure 1F shows the relation between |v| and | |  C 2 . Each symbol denotes the time average data of a cell. As expected, there is a strong positive correlation between these values (correlation coefficient r is 0.90). |v| is also positively correlated with | |  C 3 (r = 0.77, Fig. 1G). These results indicate that cells which exhibit frequent extension and contraction migrate faster. In a previous cell-crawling model 18 , |v| is conjectured to be proportional to |C 2 | and |C 3 |. However, Fig. 1D-G indicate that velocity v is positively correlated with C 2  and  C 3 , instead of C 2 and C 3 .

Migration law based on deformation.
The results in the previous section suggest that cell movement is determined by the cell-shape dynamics. Since an equation that describes the relation between cell shape and velocity has not been elucidated for persistent random motion, we explored the velocity equation based on deformation, C n and  C n , experimentally. The previous cell-crawling model 18 conjectured that velocity is determined by and C-2 is a complex conjugate of C 2 . Here, we examine the velocity equation in terms of correlation among phases (arguments) of variables, v 1 , C n , and C n  . The phases correspond to the direction of movement and deformation (see Method). When the previous model β =  Figure 2A shows the joint probability distribution of arg(v 1 ) and − C C arg( ) 2 3 . As shown, arg(v 1 ) is not correlated with arg(C −2 C 3 ). Thus, v C C c 1 2 3 β = − does not hold for our experiment. Instead, we found that arg(v 1 ) and C C arg( ) 2 3 −  are closely correlated. As shown in Fig. 2B, the joint probability distribution of arg(v 1 ) and  − C C arg( ) 2 3 has a region of high probability around the diagonal line. This result indicates that  = − v C C arg( ) arg( ) 1 23 . By analogy to the previous model, also satisfies the requirement that |v| is positively correlated with C 2  | | (Fig. 1E,F). Similarly, we found a moderate correlation SCIEntIfIC REPORTS | (2018) 8:5153 | DOI:10.1038/s41598-018-23540-x between arg(v 1 ) and C C arg( ) 2 3  − − (Fig. 2C). Figure 2C suggests that the velocity equation also includes the term . Based on the above consideration, we propose that the velocity equation where β 1 and β 2 are fitting parameters. As β 1 and β 2 increase, cell velocity increases. Thus, β 1 and β 2 represent the mobility of the cell. This equation appeared to be equivalent to the migration law for a swimming amoeboid in a 3D fluid 19 . In general, terms with a higher mode, − +  C C n n 1 and − +  C C n n 1 with n > 3, could appear in the velocity equation 19 . However, in our experiment,  C C 2 3 − and C C 2 3 −  are the largest modes, and higher modes make only a small contribution. Thus, for simplicity, we do not include higher modes in Eq. (1). The equations for higher mode, C n (n > 3), are discussed in the Supplementary Information (SI).
Here, we examine Eq. (1) experimentally, in terms of the measured velocity and deformation. We calculated the mobility β 1 and β 2 by least-squares fitting of Eq. (1). Estimated β 1 and β 2 for the three gels are listed in Table 1. As expected, β 2 is smaller than β 1 in all cases. As the elasticity of the gels increases, the mobility decreases. This result suggests that the resistance to motion increases as the elasticity of the substrate increases. By using estimated β 1 and β 2 , we can compare the velocity v 1 predicted from the model to the actual velocity V 1 = V x + iV y  Fig. 2D implies that there is a positive correlation with large fluctuation. The correlation coefficient r between v x and V x is 0.39. To clarify the detailed relation, we average V j in sections, where the boundaries of the sections are shown as dashed lines in Fig. 2D with γ = 12 μm h −1 . In Fig. 2E, the symbols and error bars indicate mean values and standard deviations of V x and V y . Blue and red bars at the lower side of the figures denote the region where 99% of the data points are included. For both V x and V y , the symbols are well aligned on the black line that represents v x = V x and v y = V y . Next, autocorrelation functions of | | v 1 and V 1 | | were compared (Fig. 2F). For comparison, autocorrelation functions of the velocity | | v c in a previous cell-crawling model 18 were also calculated, While the autocorrelation functions of | | v 1 and | | V 1 are very close, that of v c | | differs significantly. As shown in Fig. 2D, it is still hard to precisely predict the instantaneous velocity. However, on average, Eq. (1) describes the velocity dynamics well (Fig. 2E,F). Similar results were obtained for fibroblasts on 120 kPa and 410 kPa gels (see SI).
What can be predicted for cell movement based on Eq. (1)? Next, we discuss examples of the predicted dynamics in cell movement and deformation ( Fig. 3A to C). Here, we consider the case of increasing |C n | with constant phases. In this case, a cell extends its body toward a certain direction. When the first term C C (1) is dominant, arg(C 2 ) = arg(C 3 ) = 0 gives the dynamics that one of three pseudopodia extends (Fig. 3A). As a result, the cell body elongates and the centroid moves in the direction of the extended pseudopod. When the second term  C C  the case of decreasing |C n |, movement with contraction is described. Since Eq. (1) exhibits time-reverse symmetry, contraction is described as the reverse of extension. Figure 3D to I show the corresponding dynamics observed in the experiment, including both extension and contraction. As expected from Eq. (1), it is hard to distinguish between extension and contraction based solely on a snapshot of the cell shape ( Fig. 3D-I). However, there are differences in the cell-shape dynamics. For example, contraction is faster than extension. Consequently, we need to know the cell-shape dynamics to predict movement. This migration property is quite different from those in a previous cell-crawling model 18 and keratocytes 5 , where velocity is determined by the instantaneous shape.
Migration model and fitting of the experimental data. Equation (1) describes the dynamics of the velocity. However, we still don't have the governing equations for C 2 and C 3 to describe the actual movement of a cell, i.e., the dynamics of the cell trajectory and shaping, because the velocity is determined by the time evolution of C 2 and C 3 . Here, we propose a migration model by introducing time-evolution equations for C 2 and C 3 that are based on a previous cell-crawling model 18 . In the cell-crawling model, elongation and triangular deformations are induced by the force dipole F 2 and force quadrupole F 3 of the internal force acting on the periphery of the cell. We derive time-evolution equations for C 2 and C 3 by considering the force multipole and possible coupling terms (see Method): In these equations, κ 2 C 2 and K 3 C 3 cause relaxation to a circular shape, which corresponds to a restoring force. The terms α 2 v −1 C 3 and α 3 v 1 C 2 represent changes in shape due to movement. As we explain later, β 3 C −3 F 6 is introduced to fit the long tail of the probability distribution function of C 3 . γ | | C 3 3 2 is a nonlinear damping term to suppress the divergence of C 3 due to β 3 C −3 F 6 . Similar to the noise term in the persistent random walk model 8-10 , we assume that the force multipole fluctuates randomly. Here, we use red noise with cut-off frequency κ f for the force multipole, because the spatial distribution of the internal force may not change so frequently: where ξ i = ξ ix + iξ iy . ξ ix and ξ iy are white Gaussian noise; 0, ij ij kl ik jl ξ ξξ δ δ 〈 〉 = 〈 〉 = , i, k = 2, 3, … and j, l = x, y. In the experiment, cells move along the long-axis of elongation (Fig. 1D). In our model, Eqs (1) to (5) are not sufficient to reproduce such movement, because Eq. (1) does not include a coupling term between v 1 and C 2 . To resolve this discrepancy, we modify the velocity equation, Eq. (1), by adding a coupling term between v 1 and C 2 : The additional term α v v −1 C 2 causes movement along the long-axis of elongation for α v > 0 25,26 . In the previous model, velocity regulates the direction of elongation 18 . In our model, elongation regulates the direction of movement of the cell. We also add the nonlinear damping term γ | | C v 2 2 to suppress the divergence of v 1 due to α v v −1 C 2 . We confirmed that Eq. (6) is a better model than Eq. (1) in terms of Akaike's Information Criterion (AIC) 27 . For 35 kPa gel, Eq. (1) gives AIC = .
× 5 352 10 4 , and Eq. (6) gives AIC = 5 263 10 4 . × (see SI). Eventually, velocity v 1 is fully determined by cell deformations C 2 and C 3 that are randomly activated by internal forces F 2 and F 3 . Since C 2 and C 3 should have persistence through terms κ 2 C 2 and K 3 C 3 , the migration model, Eqs (2) to (7), can be referred to as a persistent random deformation (PRD) model. Figure 4A and B show the trajectories obtained by experiment and numerical simulation with Eqs (2) to (7) (see also SI and movies). Although the model does not have an inertia term in the velocity equation, the trajectory consists of persistent motion and rapid turning (Fig. 4B). Figure 4C show the cell-shape dynamics calculated in the simulation. In the simulation, we calculate C 4 , C 5 , and C 6 (see SI). Since inertial forces F 2 and F 3 fluctuate, the cell shape repeatedly extends and contracts. At the turning point, cells change the direction of the deformation. Cells that experience frequent extension and contraction move faster (Fig. 4C). Thus, the PRD model well explains the cell movement and shaping in the experiment (Fig. 1D). In the PRD model, the cell-shape dynamics and force gradually follow the fluctuation due to the time-retardation term, 1/κ i . As a result, elongation and triangular deformation persist for a few hours. Since the direction of movement is approximately parallel to the  (Fig. 4C), persistent random motion in the PRD model arises from persistent fluctuation of the cell shape, as we expected.
The detailed properties of migration dynamics are shown in Fig. 5, along with a comparison of the experimental and simulated results. In this figure, red and blue symbols represent the data obtained on 35 and 410 kPa gels, respectively. For clarity, symbols for the 410 kPa gel are shifted downward. Black and green lines are fitting curves for Eqs (2) to (7) for cells on 35 and 410 kPa gels, respectively (see Method). Note that the statistical properties of the results in Fig. 4 are shown as data for the 35 kPa gel in Fig. 5. The experimental and simulated data for 120 kPa gel are shown in the SI. The properties of migration and deformation are well-fitted by the model. For comparison, we also fitted the same experimental data by the conventional PRW model (gray lines in Fig. 5A,H,I, see Method). The PRW model reproduces the MSD and autocorrelation function of velocity (see SI). However, the model cannot explain the PDFs of velocity, persistent length, or rotation angle (Fig. 5A,H,I). Fitting by the PRW model is described in detail in the SI. Statistical properties of cell velocity, deformation and trajectory. We evaluated the statistical properties of v 1 , C 2 , and C 3 through a probability distribution function (PDF) for the migration of a cell population. Due to spatial symmetry, the real and imaginary parts of v 1 , C 2 , and C 3 obey the same distribution function. Thus, in Fig. 5A,C, we averaged the PDFs of the real and imaginary parts of the variables. The PDF of the velocity has a non-Gaussian tail (Fig. 5A). We found that the PDF of velocity approximately obeys an exponential distribution, as previously reported 28,29 . In the PRD model, the term α v v −1 C 2 in Eq. (6) causes the exponential tail of the velocity distribution (Fig. 5A). Without the nonlinear damping term, C v 2 2 γ | | in Eq. (7), the velocity would show a power law distribution 30,31 . The tail of the PDF of elongation resembles a Gaussian distribution, but the peak of the PDF is sharper than that of a Gaussian distribution (Fig. 5B). The PDF of triangular deformation has an exponential tail (Fig. 5C). In the PRD model, the term β 3 C −3 F 6 in Eq. (3) gives the exponential tail of the PDF of triangular deformation (Fig. 5C), because the term behaves as multiplicative noise 30,31 . In addition to the PDFs of v 1 , C 2 , and C 3 , those of phase differences among velocity and deformations are also well fitted (see SI). Figure 5D,E show autocorrelation functions of v 1 , C 2 , and C 3 for 35 and 410 kPa gels. The dashed lines represent fitting curves by the model. The relaxation times of C 2 and C 3 are long compared to v 1 , which denotes a long persistence of the cell shape. The non-Gaussian distribution of C 2 and C 3 (Fig. 5B,C) and the long persistence (Fig. 5D,E) indicate the persistent non-Gaussian fluctuation of cell deformation. Because the cells migrate in the direction of elongation (Fig. 1D) and the relaxation of elongation is slow (Fig. 5D,E), the trajectories of the cells tend to follow a straight line. Figure 5F shows the mean square displacement (MSD) of the cells. Similar to previous results 13,28,32 , the exponents of MSD were 1.6-1.7 at a short time interval. At a long time interval, the exponents came close to 1 and the dynamics became diffusive. The cell trajectory consisted of a series of persistent motion and rapid turning 7,13 . Thus, we investigated the persistence and turning angle of the trajectory (Fig. 5G), which are defined in the Method section. Figure 5H show the complementary cumulative distributions (CDF) of the persistent length. The CDF of the persistent length has an exponential tail, as reported previously 13,32 . The PDF of the rotation angle Δθ has a large peak at 180° (Fig. 5I), which is quite different from previous studies with persistent random walk-based models 7,13 . The angle 180° represents reciprocal motion (Fig. 1C), but the back-and-forth motion is not periodic (Fig. 5D,E). The peak decreases as the elasticity of the substrate increases. Thus, back-and-forth motion is suppressed on the stiff gel. In contrast to PRW models that violate time-reverse symmetry, the time-reverse symmetry of the velocity equation, Eq. (6), causes this strong anti-correlation of turning angles.
Dependence of the cell properties on the elasticity of the substrate. We now summarize the cell response to the elasticity of the substrate. Some of the fitting parameters for the three gels are shown in Table 2. All sets of fitting parameters are listed in SI. Table 2 shows that the mobility β 1 and β 2 significantly decreased as the elasticity of the substrate increased. This result suggests that, for a stiff substrate, cells strongly stick to the substrate 20 . The relation κ 2 < κ 3 indicates that elongation is the most unstable and 'active' Fourier mode. As the substrate become stiffer, κ 3 significantly decreases while κ 2 is constant. This implies that triangular deformation is enhanced on a stiff gel.
Finally, we performed a dimensional analysis to estimate the magnitude of internal forces. In the PRD model, we assume that the coefficients of  C 2 and C 3  are unity. As a result, the dimensions of σ 2 and σ 3 become μm h −1 . Here, we consider C n /R 0 to be the strain for the cell body, where R 0 is the mean radius of the cell. We also define E as the elasticity of the cell. If we consider that the terms κ 2 C 2 and κ 3 C 3 in Eqs (2) and (3) cause relaxation to a circular shape, these terms should correspond to the restoring force EC n /R 0 . Thus, by multiplying Eqs (2) to (4) by E/(R 0 κ n ), we can estimate the magnitude of the internal forces to be Eσ n /(R 0 κ n ). In our analysis, we can calculate the non-dimensional forces, σ n /(R 0 κ n ). Table 2 shows a constant value in the non-dimensional force dipole, σ 2 / (R 0 κ 2 ), while the non-dimensional force quadrupole, σ 3 /(R 0 κ 3 ), increases significantly on a stiffer substrate. It has been reported that the elasticity of fibroblast cells increases as the substrate becomes stiffer 23 . Eventually, internal forces should increase as the elasticity of the substrate increases 20 .

Discussion
In this study, we applied a Fourier mode analysis of the cell shape to quantitatively examine the relation between movement and shaping dynamics. The Fourier mode analysis revealed that the velocity of fibroblast cells is a function of the time derivative of the elongation and triangular deformation of the cell. We sought to apply our analytical method to other cells. For example, keratocytes are well-known migratory cells 5 . Keratocytes migrate in a manner that is quite different from that of fibroblasts. Fibroblasts migrate along the long axis of the cell body through extension and contraction. On the other hand, keratocytes migrate along the short axis of the cell body while retaining an almost constant shape 5 . Equation (1) cannot explain the migration dynamics of keratocytes because Eq. (1) with a constant shape gives a motionless state. Thus, the previous model of cell-crawling is Designation β 1 (μm −1 ) β 2 (μm −1 ) κ 2 (h −1 ) κ 3 (h −1 ) σ 2 /(R 0 κ 2 ) σ 3 /(R 0 κ 3 )  adequate for a migration model of keratocytes 18 , which assumes that the velocity is proportional to C −2 C 3 . If we combine Eq. (1) with the previous model, the general form of the migration law can be written as where β c = 0 gives a migration model of fibroblast-like cells and β 1 = β 2 = 0 gives a migration model of keratocytes. Dictyostelium cells seem to have intermediate properties. When we focus on the phases of velocity and deformation, arg(v 1 ) is correlated with arg(C −2 C 3 ) 24 . This implies that the velocity equation includes the term β c C −2 C 3 . However, in contrast to keratocytes, the shape of Dictyostelium cells varies considerably over time 6 . Thus, extension and contraction, Eq. (1), could also be included in the velocity equation. Although the above discussion is merely a conjecture, a systematic investigation of the relation between cell velocity and shape may be useful for classifying the migration type of cells.
In this work, we found that the movement of fibroblast cells on a gel surface is not explained by previous migration models. We quantitatively show the correlation between cell movement and extension/contraction of the cell body. We propose a persistent random deformation (PRD) model that is based on extension/ contraction-based movement. The PRD model shows good agreement with the statistical properties of the trajectory, velocity and shape of the cell, including persistent non-Gaussian fluctuation of deformation. By fitting experimental data to the model, we quantitatively evaluate the coefficients in Eqs (2) to (7), such as motility parameters, the relaxation time of shaping, and the magnitude of internal force. With regard to a theoretical viewpoint, it is important to clarify the physical meaning of the coefficients by solving the continuous migration model 33,34 based on the dynamics of focal adhesion and actin 35 . The movement and shaping of cells should be regulated by the dynamics of focal adhesion, stress fibers and an actin network 35 . Thus, with regard to an experiment, the simultaneous measurement of such internal structures could provide a better understanding of migration law Eq. (1) and the PRD model. Since the PRD model can be applied regardless of the details of the cell and its environment, we can evaluate the dependence of cell properties on the culture environment through estimation of the fitting parameters. In summary, the proposed model and analytical method provide a new tool for investigating the migration dynamics of cells.

Methods
Preparation of StG gel substrate. Photocurable styrenated gelatin (StG) was used as a culture substrate.
The preparation method has been described previously 21,36,37 . StG (30 wt%) and sulfonyl camphorquinone (2.5 wt% of gelatin; Toronto Research Chemicals, ON, Canada) were dissolved in phosphate-buffered saline (PBS). The mixed solution was centrifuged and aspirated to exclude deposits and dissolved oxygen. The sol solution was then conditioned using an AR-100 deforming agitator. 30 μl of the StG sol solution was spread between vinyl-silanized glass substrates (vinyl-glass) and a normal glass substrate coated with poly(N-isopropylacrylamide) (PNIPAAm, Sigma Aldrich, St. Louis, MO). The gel was then prepared by irradiation of the entire sample with visible light (62 or 63 mW/cm 2 at 488 nm; light source: MME-250; Moritex Saitama, Japan). Finally, the gels were detached from the PNIPAAm-coated normal glass substrate and washed thoroughly with PBS at 28 °C to remove adsorbed PNIPAAm. The elasticity of the gel was varied by changing the duration of irradiation from 300 s to 660 s. As we reported previously, the surface biochemical conditions were ensured to be the same for all the gel samples with different elasticities 21 . Therefore, the present experiments enabled us to investigate the pure biomechanical aspects of cell migration. The surface elasticity of the StG gel was determined by nano-indentation analysis using an atomic force microscope (JPK NanoWizard 4, JPK Instruments). A commercial silicon-nitride cantilever with a tetrahedral tip and a nominal spring constant of 0.1 N/m was used (BioLever mini, Olympus). We made three hydrogels with Young's moduli of 34 ± 18 kPa, 121 ± 46 kPa, and 412 ± 69 kPa, respectively. Cell culture. NIH 3T3 fibroblast cells were cultured in Dulbecco's Modified Eagle Medium (DMEM, Nacalai Tesque) at 37 °C in a humidified atmosphere containing 5% CO 2 .

Time-lapse observation of cell migration. The migratory motion of cells on gels was monitored using
an automated all-in-one microscope with a temperature-and humidity-controlled cell chamber (BIO REVO BZ-9000; Keyence Corporation, Osaka, Japan). Prior to the time-lapse observations, cells were seeded onto the gel surface at a density 1.5 × 10 3 cells/cm 2 and cultured for 6 hours under 5% CO 2 . Phase-contrast images of cells were captured every 5 min for 15-20 h.
Analysis of cell trajectories and cell shape. Movement trajectories and the shape of the cells were determined and analyzed using MATLAB software. Based on the edge detection of the cell, we extracted the shape of the cell from the phase-contrast images (Fig. 1D). The details of image-processing are explained in SI. We traced each cell and measured the time evolution of the trajectory and shape. If the cells collided or if a cell replicated, we stopped the trace. When the cells separated again, we renumbered the cells and restarted the trace. Thus, the cells have different durations of data. We only analyzed data that accumulated for longer than 8.3 h (35 kPa and 120 kPa gels) or 4.2 h (410 kPa gel). Through the image analysis, cell trajectories x(t) and shape R(θ, t) were calculated (Fig. 1D inset). The complex Fourier coefficient C n (t) of the spatiotemporal shape R(θ, t) is defined as where R 0 is the mean radius and m is the number of data points. The mode n = 1, C 1 (t), almost always vanishes because C 1 (t) corresponds to displacement of the centroid. Thus, C 1 (t) is included in the velocity of the centroid 25,26 . The amplitude |C n (t)| corresponds to the magnitude of deformation, and the phase φ n represents the n . Before calculating velocity v t ( ) and  C t ( ) n , we take the moving average of x t ( ) and C n (t) over 3 consecutive data points. β 1 and β 2 in Eq. (1) are estimated by minimizing the distance S between the velocity V t ( ) obtained experimentally and the velocity v t ( ) in the model; Coupling terms in the model. Based on the symmetry argument, we explain how coupling terms in the model are determined. A necessary condition for the coupling term is to satisfy the fundamental symmetries, uniformity and isotopy of the space 25,26,38 . The model equations should be invariant under a rotational transformation of coordinates. For simplicity, we define = C v 1 1 and φ φ = v 1 . When we rotate the cell θ 0 radians, Similarly, the 2 nd -order nonlinear term C C l m is transformed into θ + C C i l m exp ( ) l m 0 . In the equation for C n , we require that C n and C C l m undergo the same transformation, otherwise the equation changes under rotational transformation. Thus, n = l + m is required. For example, C −1 C 2 and C −2 C 3 are acceptable nonlinear terms in the equation of C 1 . In the same way, the 3 rd -order nonlinear term C C C k l m must satisfy the condition n = k + l + m. In the case of a coupling term that includes  C l and F m , the same condition holds.
Numerical simulation of the model. For numerical calculation of the force terms, Eq. (6), we use the standard Euler-Maruyama scheme for time discretization: To calculate Eqs (3) and (4), we use the Euler method with Δt = 0.5 min. The trajectories are calculated as In the procedure for fitting, we numerically sampled x(t), y(t), and C n (t) with a time interval of 5 min, which is identical to the time interval of the image sequence in the experiment. To include observation error 13 , we add Gaussian white noise to x(t), y(t), and C n (t); ξ ξ ξ ξ Re( ( )) Re( ( )) , Im( ( )) Im( ( )) , (11) ox oy n n rn n n in where ξ i are white Gaussian noises with standard deviation σ 0 . We estimated that the magnitude of the observation error of image processing was around 1 μm (see SI). In the fitting procedure, similar to the previous work 13 , the magnitude of the observation error σ 0 was treated as a fitting parameter. We then calculate v t ( ),  C t ( ) n after a moving average of x t ( ) and C t ( ) n n n with τ = 5 min. The properties of the trajectory are analyzed with the same procedure as that used for the experimental data.
Fitting of the experimental data by the PRD and PRW models. In this section, we explain the method for fitting of the PRD and PRW models. Both fittings are performed at the cell-population level. For the PRD model, the number of fitting parameters is large. Thus, it is difficult to find the best fit either analytically or numerically. Thus, we manually search for a set of fitting parameters that well reproduces PDFs of velocity, deformation, phase differences (see SI), autocorrelation function of velocity and deformations, persistent length and time, rotational angle, and mean square displacement. The parameter range of manual search is shown in SI. Next, we numerically search the set of fitting parameters that locally minimize an error function ERR around the manually searched fitting parameters.
where R i 2 are coefficients of determination. Here we use R 2 of quantile-quantile plots of experimentally-measured and numerically-calculated PDFs of velocity, deformations, and phase differences. We also include R 2 of the autocorrelation functions of velocity and deformations in the error function. Since the numerical results slightly fluctuate due to the finite size of the data points, the calculated local minimum of ERR slightly fluctuates from time to time. Thus, we searched for the local minimum 20 times. In Table 2, we show the average and standard deviation of fitting parameters that give a local minimum of ERR.
For fitting by the PRW model, we fit the mean square displacement of experimental data to that of the analytical solution of the PRW model 13 (see SI). In the fitting procedure, we numerically minimize the weighted residual sum of squares by using a nonlinear programming solver in MATLAB software.
Calculation of persistence and turning angle. To extract the persistent motion from the trajectory, we used a time series of velocity correlation CV(t) with a short time interval; x x y y where Δt = 10 min. CV(t) has a small value for two cases; the cell shows a rapid change in the direction of migration, or the cell almost stops. Therefore, we define that motion is persistent when CV(t) has a large value. For SCIEntIfIC REPORTS | (2018) 8:5153 | DOI:10.1038/s41598-018-23540-x the threshold of persistent motion, a value that was 2/3 of the median of CV(t) was used, where the median was taken for each trajectory. In Fig. 5G, the trajectory is drawn with red lines when the cell moves persistently. We define the persistent length as the length of the red lines in Fig. 5G. We also define the rotation angle as the angle between two successive persistent trajectories 7 , shown as Δθ in Fig. 5G.
Data and materials availability. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.