Mechanochemical subcellularelement model of crawling cells
Abstract
Constructing physical models of living cells and tissues is an extremely challenging task because of the high complexities of both intra and intercellular processes. In addition, the force that a single cell generates vanishes in total due to the law of action and reaction. The typical mechanics of cell crawling involve periodic changes in the cell shape and in the adhesion characteristics of the cell to the substrate. However, the basic physical mechanisms by which a single cell coordinates these processes cooperatively to achieve autonomous migration are not yet well understood. To obtain a clearer grasp of how the intracellular force is converted to directional motion, we develop a basic mechanochemical model of a crawling cell based on subcellular elements with the focus on the dependence of the protrusion and contraction as well as the adhesion and deadhesion processes on intracellular biochemical signals. By introducing reactiondiffusion equations that reproduce traveling waves of local chemical concentrations, we clarify that the chemical dependence of the cellsubstrate adhesion dynamics determines the crawling direction and distance with one chemical wave. Finally, we also perform multipole analysis of the traction force to compare it with the experimental results. To our knowledge, our present work is the first study that accomplishes fully forcefree migration utilizing intracellular chemical reactions. Although the detailed mechanisms of actual cells are far more complicated than our simple model, we believe that this mechanochemical model is a good prototype for more realistic models.
I Introduction
Cells are the basic units of living organisms. Despite the diversity of cellular types, they share similar structures. Biological cells are composed of a number of polymers, proteins, and lipids, which form stable structures such as the lipid membrane. At the same time, cells exist in a state far from equilibrium. Inside a cell, complex chemical reactions take place and are converted to mechanical forces by molecular motors. Consequently, biological cells spontaneously exhibit various dynamics, such as locomotion and proliferation.
In general, objects that exhibit spontaneous motion are called active matter Lauga and Powers (2009); Ramaswamy (2010); Vicsek and Zafeiris (2012); Cates (2012); Marchetti et al. (2013); Bechinger et al. (2016). In contrast to objects passively driven by external forcing, active matter, including living cells, generates force in itself, which is characterized by a vanishing force monopole due to the actionreaction law. Under this forcefree condition, it is necessary to break symmetry to achieve spontaneous motion, such as directional motion. For microorganisms that swim in a fluidic environment, the scallop theorem Purcell (1977) describes the importance of breaking reciprocality to achieve net migration via internal cyclic motions.
In nature, there exist a number of microorganisms that crawl on substrates, such as the extracellular matrix and other cells. In contrast to the locomotion of microswimmers, adhesion to the substrate plays an important role in the locomotion of crawling microorganisms. Such crawling motion is often observed in eukaryotic cells, such as Dictyostelium cells and keratocytes.
Many models have been introduced to explain various aspects of the dynamics of crawling cells. Examples include the vertex model Honda (1978); Honda et al. (2004, 2008); Bi et al. (2016), the cellular Potts model Nishimura et al. (2009); Niculescu et al. (2015), the continuous model Marchetti et al. (2013); Nier et al. (2016), the phase field model Ziebert and Aranson (2016); Shao et al. (2010); Taniguchi et al. (2013); Shi et al. (2013); Ziebert et al. (2012); Tjhung et al. (2015), and the particlebased model Newman (2007); Sandersius and Newman (2008); Basan et al. (2011); Zimmermann et al. (2016); Smeets et al. (2016). Here, we construct a model based on the particlebased model.
In our previous study Tarama and Yamamoto (2018), we systematically investigated the cycle of the typical crawling mechanism Ananthakrishnan and Ehrlicher (2007): 1) protrusion at the leading edge, 2) adhesion to the substrate at the leading edge, 3) deadhesion from the substrate at the trailing edge, and 4) contraction at the trailing edge. To clarify the role of this cycle in efficient crawling motion, we introduced a simple mechanical model in which a cell is described by two subcellular elements connected by a viscoelastic bond, which includes an actuator that elongates and shrinks cyclically. The substrate friction of each element switches cyclically between the adhered stick state and the deadhered slip state. By tuning the phase shifts between the actuator elongation and the substrate friction of each element, we demonstrated that the order of the four basic processes of the typical crawling mechanism has a great impact on the crawling distance and efficiency, as well as the crawling direction.
If we consider the extension of the mechanical model to a model cell consisting of far more than two elements or to two or threedimensional (2D or 3D) space, adjusting the phase shift of each element “by hand” becomes less feasible. Instead, we need to consider the underlying processes that regulate the actuator elongation and the stickslip transition of the substrate friction.
The aim of this paper is to develop a basic particlebased model for cell crawling. To this end, we consider a cell crawling on a flat substrate and extend our previous mechanical model to two dimensions. We describe a cell by a set of many subcellular elements connected by viscoelastic bonds Newman (2007). In addition, intracellular chemical reactions are represented by simple reactiondiffusion (RD) equations Murray (2002, 2003); Kuramoto (1984); Epstein and Pojman (1998); Pismen (2006), which trigger mechanical activities. We then couple the RD equations and mechanical models to achieve efficient migration. In particular, we focus on the time delay between the intracellular chemical reactions and cell mechanics, which corresponds to the ordering of the basic crawling processes.
This paper is organized as follows. In the next section, we introduce the model equations that couples cell mechanics and intracellular chemistry. Then, we show that the dependence of the substrate adhesion on the intracellular chemistry determines the direction of the cell crawling in Sect. III. In Sect. IV, we investigate how the cell crawling changes depending on the mechanical and chemical signals that control the substrate adhesion. The impact of the cell shape and size are studied in Sets. V and VI, respectively. In Sect. VII, random crawling motion is realized by random excitation of intracellular chemistry, for which we analyze the traction force multipoles in Sect. VIII. Sect. IX is devoted to the summary and discussion, and this paper concludes with Sect. X.
Ii Model
First, we introduce our mechanical model of a crawling cell and the RD equations representing intracellular chemical reactions. The choice of RD equations is arbitrary, and we employ a previously introduced model. We then couple the mechanical model and the RD equations, which regulate the intracellular mechanical activities. In particular, we confine ourselves to studying possible couplings between the intracellular chemical and mechanical models.
ii.1 Subcellularelement model
We describe a single cell by a set of subcellular elements Newman (2007) connected by KelvinVoigt type viscoelastic bonds, as schematically depicted in Fig. 1. Since the typical size of a cell is on the order of ten micrometers, the effect of inertia is negligible. Then, the force balance equation of element is given by
(1) 
where is the velocity of the element located at the position . Here, the abbreviations and are used for the relative position . The summation is over the set of the elements connected to the element by the viscoelastic bonds.
The first term on the lefthand side of Eq. (1) represents the substrate friction with coefficient , which changes over time due to intracellular activity. The second term represents intracellular dissipation with the rate . The first term on the righthand side represents intracellular elasticity with the elastic modulus and free length . Intracellular activity is also included in the actuator, which tends to elongate the connecting bond by changing the free length over time as . Here, represents the actuator elongation, from which the force generated by the actuator is calculated as . We emphasize that the model Eq. (1) satisfies the forcefree condition since the intracellular force acts symmetrically on the pair of subcellular elements. Namely, the sum of the intracellular force in Eq. (1) vanishes as
(2) 
where
(3) 
is the intracellular force acting on the element .
The last term on the righthand side of Eq. (1) prevents the collapse of the subcellular element network, which is given by , where . This potential penalizes shrinking of the area of each triangle formed by connected subcellular elements , , and , which is defined by with as the unit vector perpendicular to the 2D substrate.
We scale the system by for length and for time, which are physiologically relevant values for typical living cells Maeda et al. (2008); Tanimoto and Sano (2014). In addition, the scale of the force is set to , which is on the order of the traction force that cells exert on the substrate. The typical values of the mechanical parameters of the model Eq. (1) are summarized in Appendix B.
ii.2 Chemical reaction
In the model Eq. (1), the effects of the intracellular activities are included in the actuator elongation and the change in the substrate friction coefficient . The former represents the protrusion and contraction processes. The latter corresponds to the adhesion and deadhesion of the cell to the underlying substrate. In actual cells, such cellular activities are caused by various intracellular chemical signals. However, it is not realistic to include all chemical components and their signaling pathways. Therefore, we model the intracellular chemical reactions by simple RD equations.
Here, we employ the RD equations proposed by Taniguchi et al. Taniguchi et al. (2013), which are twocomponent activatorinhibitor equations:
(4)  
(5) 
where and represent the inhibitor and activator concentrations for subcellular element , respectively. In the original paper Taniguchi et al. (2013), Eq. (5) were introduced to model the phosphoinositide signaling pathway of Dictyostelium cells, and the activator and the inhibitor corresponded to phosphatidylinositol (3,4,5)trisphosphate (PIP3) and phosphatidylinositol (4,5)bisphosphate (PIP2) concentrations, respectively. The details of these RD equations are given in Appendix C.
The important property of the RD equations, Eq. (5), is that they are of the GreyScott type Gray and Scott (1983, 1984). One of the advantages of the GreyScott model is that it can show either an excitable or a bistable nature depending on the parameters. Taniguchi et al. claimed that the signaling pathway that they were modeling was excitable, and thus, they considered the parameter region of the excitable case to successfully reproduce the experimental results Taniguchi et al. (2013). Interestingly, similar RD equations were studied by Shao et al. Shao et al. (2010) in the context of cell crawling. However, they assumed a bistable regime to reproduce the steady migration of keratocyte cells.
In this paper, we consider the excitable case with the parameters summarized in Table 3, following the study by Taniguchi et al. Taniguchi et al. (2013). The Laplacian terms in Eq. (5) are calculated by using the moving particle semiimplicit (MPS) method Koshizuka et al. (1998). See also Appendix C for further information.
ii.3 Mechanochemical coupling
To combine the cell mechanics, Eq. (1), and the RD equations, Eq. (5), we consider the coupling of the chemical concentrations to the actuator elongation and the substrate friction coefficient individually.
Before introducing the coupling, we tested a specified traveling wave that couples to the actuator elongation and the substrate friction change. We introduced a time delay between the elongation and substrate friction change. The results showed that there exist an optimum time delay and optimum wavelength for which the cell exhibits the largest migration distance, as summarized in Appendix D.
ii.3.1 Actuator elongation
First, we introduce the coupling between the RD equations and the actuator elongation. In Ref. Taniguchi et al. (2013), actin polymerization was found to be enhanced with increasing PIP3 concentration. Therefore, we presume that the actuator elongation depends on the PIP3 concentration as
(6) 
where is the mean PIP3 concentration for the bond connecting the elements and . Although is a positive quantity, its maximal value depends on the strength of the initial fluctuation because of the excitable nature of Eq. (5). Therefore, is introduced on the right hand side of Eq. (6) to prevent an extremely large elongation. is a constant denoting sensitivity, and is the magnitude of the elongation. Here, we set and .
ii.3.2 Substrate adhesion
Next, we consider the adhesion to the substrate underneath and the deadhesion from it. We model the adhesion/deadhesion processes by the transition of the substrate friction coefficient between the adhered stick state and the deadhered slip state. Here, we consider the dependence of the substrate friction coefficient on both mechanical and chemical signals:
(7) 
where is the time delay. takes a value between 0 and 1, representing the ratio of the mechanical and chemical dependence of the stickslip transition of the substrate friction. See also Fig. 2 for the plot of these functions.
The function is defined by
(8) 
where and are the substrate friction coefficients in the adhered stick state and the deadhered slip state, respectively. The small parameter indicates the sharpness of the adhesiondeadhesion transition. Here, we set the transition sharpness to . We note that, if there are no change in the substrate friction coefficient, the cell does not exhibit any translational motion.
If we consider an artificial vesicle or droplet sitting on a substrate, its adhesion strength changes depending on the force acting on it Schwarz and Safran (2013). The term in Eq. (7) represents this dependence of the cell adhesion to the substrate. Here, instead of the force acting on each subcellular element, we presume that the local velocity changes the adhesion strength through
(9) 
where is the threshold value. The subcellular element tends to adhere to the substrate () if the speed is smaller than the threshold value, i.e., , while the element slips on the substrate () if . We set the threshold value to .
The formula of the stickslip transition of the cellsubstrate friction depending on the local velocity, i.e., Eq. (7) with , was introduced in Ref. Barnhart et al. (2015). As a result of the balance of the two functions, and , the substrate friction switches between the stick and slip states with a sharp transition.
In addition to the mechanical dependence, the adhesion strength of a cell can change depending on its internal chemical conditions Schwarz and Safran (2013). Since the molecular details of cell adhesion are complicated, we assume here that it changes depending on the PIP3 concentration as the actuator elongation:
(10) 
In Eq. (10), stands for the sensitivity, and is the threshold concentration. Due to this term, a large value of prevents strong adhesion if , while large enhances the adhesion if . However, we are not sure whether PIP3 enhances or diminishes adhesion.
Iii Crawling by direct and retrograde waves
To determine the sign of , we numerically integrated Eqs. (1)–(10) for both positive and negative . See Appendix A for the details of the simulation. Fig. 3 depicts a time series of snapshots of a crawling cell for in Fig. 3(a) and in Fig. 3(b), respectively. Here, we set the threshold in Eq. (10) to throughout this paper.
If is positive, the cell moves in the opposite direction to the PIP3 traveling wave, as shown in Fig. 3(a). Interestingly, however, if the sign of is negative, a qualitatively different result appears: namely, the cell starts to move in the same direction as the traveling wave, as displayed in Fig. 3(b).
With respect to the migration direction, the traveling wave in the same direction is called the direct wave, while the one in the opposite direction is referred to as the retrograde wave Iwamoto et al. (2014). In this sense, the above crawling motion for positive in Fig. 3(a) corresponds to the motion with the retrograde wave, and the one for in Fig. 3(b) corresponds to the motion with the direct wave. Since the experiments in Ref. Taniguchi et al. (2013) show that cells move in the direction in which PIP3 concentration is increased and thus actin polymerization is enhanced, we set for the rest of this paper.
Iv Mechanical vs. chemical control of adhesion
Now, we study the effect of the mechanical and chemical dependence of substrate friction coefficient on cell crawling. To characterize the cell migration, we measure the migration distance of the cell in one cycle, i.e., with one traveling wave, varying the time delay . Varying between 0 and 1 tunes the ratio of the mechanical and chemical dependences.
The results are summarized in Fig. 4(a), where different lines correspond to different values of as indicated in the legend. Interestingly, the mixing of the mechanical and chemical dependences of the substrate friction may result in larger values of than purely mechanical or chemical control. In Fig. 4(a), reaches approximately 0.36 for and . This value, however, is smaller than the crawling of a Dictyostelium cell, which moves its body length in approximately two cycles Tanimoto and Sano (2014).
V Impact of cellular shape
Thus far, we have assumed a hexagonal cell shape, where the structure of the subcellular elements is given by a perfect hexagonal lattice when the cell is at rest. However, this structure does not describe real cells, which are instead circular or often of more complicated shapes. To elucidate the impact of the cell shape on the crawling motion, we prepare a cell of circular shape as described in Appendix E. We then perform the simulation and measure the migration distance for different values of . The results are plotted in Fig. 4(b) and are qualitatively the same as those of the hexagonal cell in Fig. 4(a).
Vi Impact of cell size
Now, we study the impact of the size of the cell. We prepare a hexagonal cell of size , which is approximately twice the area of the previous cells. We again measure the migration distance , which is normalized by the cell length to facilitate comparison with the previous results for . The results are summarized in Fig. 4(c). Qualitatively, the tendency is the same as that of the results in Fig. 4(a) for . Namely, the migration distance can be larger if the substrate friction depends both on the mechanical and chemical signals than if it depends on only either one of them.
Vii Random excitation
In reality, cells change their migration direction over time. In our model, we can reproduce such motion by introducing stochasticity, which may originate from, e.g., the complexity of intracellular processes. Here, we randomly choose one element in every and add to that element the stimulus with intensity .
The results are summarized in Fig. 5. Due to the stochasticity, the cell changes its migration direction frequently, as shown in the trajectory of the centerofmass position in Fig. 5(a). From the snapshots of the cell in Fig. 5(b), we can see that the migration direction depends on the position at which the chemical wave occurs. Because of the excitable nature of the RD equations, a stimulus on the element that has relaxed back to the resting state is more likely to be the origin of the next wave. Therefore, a new wave tends to originate from elements that are near the origin of the previous wave. As a result, the cell tends to maintain the same migration direction for a time, as shown in Fig. 5(a). In Figs. 5(a) and (b), the parameter values are kept the same as in Fig. 3(a).
If the parameter is slightly increased from 5 to 5.5, the cell changes its migration direction more frequently, as depicted in Fig. 5(c). Depending on the random stimuli, the cell switches from directional motion to spinning motion as a spiral wave appears, as shown in Fig. 5(d). The spinning motion is rather stable, but the cell can also switch back to directional motion in response to a stimulus. Note that the RD equations maintain their excitable nature at this parameter value.
Viii Traction force multipoles
In many experiments of cell crawling, traction force is measured Style et al. (2014) since it is of fundamental importance for cell motility. Here, we measure the traction force of our model cell. Since traction force is the force that a cell exerts on the substrate underneath, it is described by the interaction between the cell and the substrate: . We calculate the traction force force multipoles, which are defined as follows, to compare with the experimental result Tanimoto and Sano (2014).
First, the traction force monopole is defined by
(11) 
where represents the subcellular element and the summation runs over the entire cell. The subscript indicates spatial component . Here, the traction force multipoles are calculated in the comoving coordinate on the cell, and thus, and 2 represent the components parallel and perpendicular to the centreofmass velocity, respectively, to compare with the experimental results in Ref. Tanimoto and Sano (2014). In the numerical simulation of the model crawling cell, the traction force monopole is equal to 0, as shown in Fig. 6, which is consistent with the experimental result Tanimoto and Sano (2014). This is readily understood from the force balance equation, Eq. (1), and the forcefree condition, Eq. (2)(a). That is, the traction force monopole vanishes because of the fact that the inertia is negligibly small for crawling cells on top of the forcefree condition.
The next lowest mode is the traction force dipole, the element of which is defined by
(12) 
On the one hand, from Fig. 6(a), the diagonal components of the traction force dipole are oscillating around 0. Here, note that the positive and negative force dipoles represent the extensile and contractile force. In the experiment Tanimoto and Sano (2014), only the contractile traction force dipole was observed, which our results fail to reproduce. The reason why our model does not reproduce the contractile force dipole is not clear yet. One possibility is that the friction coefficient of the protrusion process may be different from that of the contraction process, which are set to the same in our current model.
On the other hand, the offdiagonal components and take the same values, indicating that the traction force dipole is symmetric, although such symmetry is not presumed in its definition, Eq. (12). Such symmetric property of the traction force dipole was also obtained in the experiment Tanimoto and Sano (2014).
Now, we consider the meaning of the symmetric property of the traction force dipole that is obtained in our simulation as well as in the experiment. Actually, this symmetric property of the traction force dipole indicates the torquefree nature of the cell. Here, the torque is defined by
(13) 
which, in a two dimensional space, becomes
(14) 
by using Eq. (12). Note that, if one describes the force dipole tensor with its invariant, the torque appears as the imaginary part of the eigenvalues.
Finally, in Fig. 6(b), we plot the time evolution of the traction force dipole and the traction force quadrupole , which was also measured in the experiment Tanimoto and Sano (2014). Here the traction force quadrupole is defined by
(15) 
Interestingly, the trajectory in  space shows a counterclockwise rotation, which is qualitatively consistent with the experimental results Tanimoto and Sano (2014).
Ix Summary and Discussion
To summarize, we have constructed a mechanochemical model of a cell crawling on a substrate. The mechanical part is described by a subcellularelement model, and the chemical part is described by RD equations. To combine them, we introduce two mechanical activities. One is the actuator of the bond connecting each pair of the subcellular elements, which elongates depending on the intracellular chemical concentration. The other is the substrate friction coefficient of each subcellular element, which shows a sharp transition between the adhered stick state and the deadhered slip state. We consider the dependence of the substrate friction coefficient on both the local velocity and the intracellular chemical concentration. We also introduce a time delay of the substrate friction change.
By using this model, we clarified that the substrate adhesion dynamics affect how the intracellular force is converted cell crawling motion. Depending on the sign of the sensitivity of the substrate friction coefficient to the PIP3 concentration, the model cell exhibited crawling with the retrograde flow or with the direct flow. For the former case, which is consistent with experimental observations, our model showed that there is an optimum time delay and that the combined effect of the mechanical and chemical signals on the substrate friction coefficient can increase the migration distance. We also investigated the impact of the cell shape and the cell size, which led to qualitatively the same results. In addition, we included stochasticity in the RD equations, enabling the cell to change its migration direction and to change its dynamical mode from translational motion to spinning motion. Further, we performed multipole analysis of the substrate traction force, which was qualitatively consistent with the experimental results except the contractile nature of the traction force dipole.
Finally, we discuss some possible extensions of our current model.
 Contraction process

In our current model, the protrusion and contraction processes are both modeled by the actuator elongation, , of the bond connecting two subcellular elements. The two processes are distinguished by the sign of . In this paper, however, we consider only the protrusion process, i.e., , which is related to the PIP3 concentration. One reason of this is that the chemical reactions that regulate contraction process are not well understood yet. Therefore, in principle, we can also consider the contraction process by introducing dependence on relevant chemical concentrations.
 Adhesion dynamics

The cell adhesion is simply modeled by the switching of the substrate friction coefficient of each subcellular element in our model. However, in real cells, adhesion is mediated by adhesion molecules, which can diffuse and form focal adhesions. To represent these processes of adhesion molecules, we include detailed dynamics of the concentration of adhesion molecules and their diffusion to other subcellular elements. Then, we can discuss detailed structures such as the footsteplike focal adhesion observed for Dictyostelium cells Tanimoto and Sano (2014).
 Shape deformation

Our model cell shows a lateral expansion with respect to the crawling direction, as shown in Fig. 3(a). However, real cells, e.g., Dictyostelium cells, tend to elongate in the direction of motion Maeda et al. (2008); Bosgraaf and Van Haastert (2009). One possible reason that our current model fails to reproduce this elongated shape is that actuator elongation depends only on the absolute value of . Therefore, this inconsistency may be resolved by, e.g., introducing dependence on the gradient of .
 Three dimensions

In this paper, we modeled a cell as a twodimensional network of viscoelastic springs by assuming crawling on a flat substrate. In reality, however, cells are three dimensional object. The extension of our current model to three dimensions is straightforward.
X Conclusion
In conclusion, the modeling of crawling cells is still a challenging task due to the complexity of intra and intercellular processes. The force that a cell generates should satisfy the forcefree condition, where the total force vanishes. To achieve net migration under the forcefree condition, the cell needs to break symmetry.
In our mechanochemical subcellularelement model, the intracellular force acts symmetrically on each pair of subcellular elements; therefore, it naturally satisfies the forcefree condition. Symmetry breaking occurs due to the switching of the substrate friction coefficient between the adhered stick state and the deadhered slip state. Therefore, our model clearly distinguishes intracellular force and external force. To control those mechanical activities, we included RD equations representing intracellular chemical reactions.
The RD equations that we employed in this study were introduced to explain the chemical traveling wave observed within Dictyostelium cells Taniguchi et al. (2013). However, a number of chemical reactions occur inside a cell, and which chemical reactions are relevant may depend on the phenomena of interest. Nevertheless, we believe that our model can provide a basic framework for the future construction of mechanochemical models of crawling cells by replacing the RD equations with suitable ones for each specific phenomenon.
Acknowledgements.
This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (17H01083, 19H14673) grant, as well as the JSPS bilateral joint research projects.Appendix A Simulation method
We integrate Eqs. (1)–(7) numerically by using the Euler method with time increment . We add the stimulus to Eq. (5) on one of the subcellular elements, which we call the activator element. We set the activator element to the rightmost element denoted by the star in Fig. 1(a), except in the case of random motion, for which an activator element is chosen randomly every time . Since the traveling wave is generated by the stimuli on the activator element, we excite the activator element every so that the wave travels to the other edge of the cell within that period, and all the subcellular elements relax back to the resting state when the next wave is generated. The intensity of the initial stimulus is set to .
Appendix B Parameter values
The typical values of the mechanical parameters that appear in the model Eq. (1) are summarized in Table 1 and Table 2. The length of the cell at rest is set to , and the total number of subcellular elements is set to , which gives the total number of bonds as , unless otherwise stated. For a larger cell, we set , , and . For a cell consisting of a perfect hexagonal lattice, as in Fig. 1(a), the element distance at rest is set to . Then, the mean bond length is , where is the total bond length. The elastic modulus of each bond is set to and the intracellular dissipation rate to , which gives the dimensionless parameter . From the mechanical parameter values of typical living cells in Table 2, the value of this dimensionless parameter is estimated about 0.1–1. The substrate friction coefficient at the deadhered slip state is set to , so that the dissipation by the substrate friction is on the same order as the intracellular dissipation, . The substrate friction coefficient at the adhered stick state is set to be ten times larger than at the deadhered slip state, i.e., . The value of in is set to .
Model parameters  Simulation values 

Diameter of cell at rest,  1 (1.4 for a larger cell) 
Number of subcellular elements,  91 (169 for a larger cell) 
Number of bonds,  240 (462 for a larger cell) 
Bond length at rest,  
Total bond length at rest,  
Mean bond length,  
Substrate friction coefficient  
in the slip state,  
in the stick state,  
Elastic modulus of each bond,  10 
Intracellular dissipation rate,  
Dimensionless parameter, 
Mechanical parameters  Model  Living cells  References 

Migration speed  Maeda et al. (2008)  
Traction stress  Tanimoto and Sano (2014)  
Traction force  (estimated)  
Cellsubstrate friction coefficient  (estimated)  
Elastic modulus  10–  Bausch et al. (1999); Micoulet et al. (2005)  
Elastic constant  –  (estimated)  
Dissipation rate  
Fluid/solid transition time  Kasza et al. (2007) 
Appendix C Intracellular RD equations
Here, we explain the RD equations proposed by Taniguchi et al. Taniguchi et al. (2013). Taniguchi et al. gathered from their experimental observations of Dictyostelium cells that phosphatidylinositol (3,4,5)trisphosphate (PIP3) promotes actin polymerization and protrusion of the cellular membrane, and therefore, they considered a signaling pathway around PIP3 including phosphatidylinositol (4,5)bisphosphate (PIP2), PI3kinase (PI3K), and phosphatase and tension homolog (PTEN). By eliminating the dynamics of PI3K and PTEN, they obtained the following set of RD equations Taniguchi et al. (2013):
(16)  
(17) 
where the reaction terms are defined as
(18)  
(19) 
The global couplings are given by
(20) 
where the summation is over all subcellular elements. Here, and represent the PIP2 and PIP3 concentrations for subcellular element , respectively.
Some experiments on Dictyostelium cells Taniguchi et al. (2013); Arai et al. (2010); Shibata et al. (2012) revealed that the actin filaments accumulate around the region where the PIP3 concentration increases, and thus, PIP3 promotes actin polymerization, which induces membrane protrusion. Therefore, we assumed that the actuator elongation depends on the PIP3 concentration as in Eq. (6).
The parameters we used in this paper are summarized in Table 3. Here, the diffusion coefficient corresponds to , which is within the range of the experimentallyobserved diffusion coefficient of proteins inside a cell near the membrane Golebiewska et al. (2008). We note that the choice of the other parameters is arbitrary, but what is of more importance than their absolute values is the nature of the RD equations, namely, the excitability, which is apparent from the nullclines and the flow field in the – space, as plotted in Fig. 7.
To integrate Eq. (17), we calculate the Laplacian by using the moving particle semiimplicit (MPS) method Koshizuka et al. (1998). That is, for the chemical component , the diffusion term is modeled by
(21) 
where , and is the number of neighboring elements around element . The weight function that we employed is defined by
(22) 
where is the cutoff length, which is set to . The normalization factor is given by . By using this method, we checked whether the traveling and spiral waves are formed in the absence of the mechanical changes. To generate the wave, we add the stimulus to Eq. (17) on the activator element, assuming that the phosphorylation of PIP2 to PIP3 is enhanced for this element.
Model Parameters  Values used in simulations 

0.48  
0.48  
240  
90  
5  
5  
30  
6  
30 
Appendix D Sinusoidal traveling wave
In this appendix, we show the cell crawling motion induced by a specified intracellular chemical traveling wave.
For simplicity, we consider only a single intracellular chemical component , which we assume travels periodically inside the cell. The chemical component for element changes as follows:
(23) 
where is the maximum concentration, which is set to so that .
For simplicity, we assume a harmonic plane wave that travels in the direction. Then, the phase of the chemical concentration changes as follows:
(24) 
where and are the frequency and the wavenumber of the traveling wave, and is the position of the activator element, from which the traveling wave occurs.
We assume that the cellular activity depends on the value of the intracellular chemical concentration . The actuator elongation depends on the mean chemical density:
(25) 
where represents the magnitude of the elongation. In addition, the substrate friction coefficient of element changes depending on . Here, we assume the substrate friction to be described by a twostate function that switches sharply between the adhered stick state and the deadhered slip state as follows:
(26) 
where is a phase shift, is an integer that satisfies , and and are the values of the friction coefficient in the stick and slip states, respectively. Here, since the internal motion, i.e., the elongation and the stickslip transition, is perfectly cyclic, we consider the phase shift instead of a time delay. The two are actually equivalent.
We carried out numerical simulations of Eq. (1) together with Eqs. (23)–(26) with the time increment . In the simulation, we choose the element indicated by the star in Fig. 1(a) as the activator element, and we set so that the wave travels from right to left. The period of the traveling wave is set to the unit time scale, , and thus, . We let the actuator elongate twice as long as its equivalent length: . For these parameters, we vary the wavenumber and the phase shift and measure the migration distance for one cycle.
The results are summarized in Fig. 8. In Fig. 8(a), the migration distance for one cycle is plotted against the phase shift for different values of , and time series of a typical crawling cell are depicted in Fig. 8(b) and Fig. 8(c).
Fig. 8(a) shows clearly that the migration direction changes depending on the phase shift . The positive represents forward motion, i.e., migration towards the activator element. That is, the cell migrates against the intracellular traveling wave, as depicted in Fig. 8(b). In contrast, the negative corresponds to backward motion, where the cell migrates in the same direction as the intracellular traveling wave, as shown in in Fig. 8(c). Therefore, there are maximum and minimum values of , where the cell achieves maximum migration distance in the forward and backward directions, respectively.
The migration distance also depends on the wavenumber . There is an optimal around for which the cell can achieve the largest migration distance. That is, in this case, the cell is in a fully stretched state.
Appendix E Circular cell
For cells with a circular shape at rest, we prepared the subcellular elements and the connecting bonds as depicted in Fig. 9(a). The length of the cell was set to , and the total number of subcellular elements was set to . The total number of bonds connecting the subcellular elements was then . We set of each bond of the circular cell to the initial element distance, as shown in Fig. 9(a). Then, the mean bond length was calculated to be . The other parameters were kept the same as for hexagonal cells.
References
 Lauga and Powers (2009) E. Lauga and T. R. Powers, Reports on Progress in Physics 72, 096601 (2009).
 Ramaswamy (2010) S. Ramaswamy, Annual Review of Condensed Matter Physics 1, 323 (2010) .
 Vicsek and Zafeiris (2012) T. Vicsek and A. Zafeiris, Physics Reports 517, 71 (2012), collective motion.
 Cates (2012) M. E. Cates, Reports on Progress in Physics 75, 042601 (2012).
 Marchetti et al. (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013).
 Bechinger et al. (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
 Purcell (1977) E. M. Purcell, Am. J. Phys 45, 3 (1977).
 Honda (1978) H. Honda, Journal of Theoretical Biology 72, 523 (1978).
 Honda et al. (2004) H. Honda, M. Tanemura, and T. Nagai, Journal of Theoretical Biology 226, 439 (2004).
 Honda et al. (2008) H. Honda, N. Motosugi, T. Nagai, M. Tanemura, and T. Hiiragi, Development 135, 1407 (2008) .
 Bi et al. (2016) D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Phys. Rev. X 6, 021011 (2016).
 Nishimura et al. (2009) S. I. Nishimura, M. Ueda, and M. Sasai, PLoS Comput Biol 5, e1000310 (2009).
 Niculescu et al. (2015) I. Niculescu, J. Textor, and R. J. de Boer, PLOS Computational Biology 11, 1 (2015).
 Nier et al. (2016) V. Nier, S. Jain, C. T. Lim, S. Ishihara, B. Ladoux, and P. Marcq, Biophysical Journal 110, 1625 (2016).
 Ziebert and Aranson (2016) F. Ziebert and I. S. Aranson, npj Computational Materials 2, 16019 (2016).
 Shao et al. (2010) D. Shao, W.J. Rappel, and H. Levine, Phys. Rev. Lett. 105, 108104 (2010).
 Taniguchi et al. (2013) D. Taniguchi, S. Ishihara, T. Oonuki, M. HondaKitahara, K. Kaneko, and S. Sawai, Proceedings of the National Academy of Sciences 110, 5016 (2013) .
 Shi et al. (2013) C. Shi, C.H. Huang, P. N. Devreotes, and P. A. Iglesias, PLOS Computational Biology 9, 1 (2013).
 Ziebert et al. (2012) F. Ziebert, S. Swaminathan, and I. S. Aranson, Journal of The Royal Society Interface 9, 1084 (2012) .
 Tjhung et al. (2015) E. Tjhung, A. Tiribocchi, D. Marenduzzo, and M. Cates, Nature communications 6 (2015).
 Newman (2007) T. J. Newman, Modeling multicellular structures using the subcellular element model, in SingleCellBased Models in Biology and Medicine, edited by A. R. A. Anderson, M. A. J. Chaplain, and K. A. Rejniak (Birkhäuser Basel, Basel, 2007) pp. 221–239.
 Sandersius and Newman (2008) S. A. Sandersius and T. J. Newman, Physical Biology 5, 015002 (2008).
 Basan et al. (2011) M. Basan, J. Prost, J.F. Joanny, and J. Elgeti, Physical Biology 8, 026014 (2011).
 Zimmermann et al. (2016) J. Zimmermann, B. A. Camley, W.J. Rappel, and H. Levine, Proceedings of the National Academy of Sciences 113, 2660 (2016) .
 Smeets et al. (2016) B. Smeets, R. Alert, J. Pešek, I. Pagonabarraga, H. Ramon, and R. Vincent, Proceedings of the National Academy of Sciences 113, 14621 (2016) .
 Tarama and Yamamoto (2018) M. Tarama and R. Yamamoto, Journal of the Physical Society of Japan 87, 044803 (2018) .
 Ananthakrishnan and Ehrlicher (2007) R. Ananthakrishnan and A. Ehrlicher, International journal of biological sciences 3, 303 (2007).
 Murray (2002) J. D. Murray, Mathematical Biology I. An Introduction, 3rd ed., Interdisciplinary Applied Mathematics, Vol. 17 (SpringerVerlag New York, 2002).
 Murray (2003) J. D. Murray, Mathematical Biology II. Spatial models and biomedical applications, 3rd ed., Interdisciplinary Applied Mathematics, Vol. 18 (SpringerVerlag New York, 2003).
 Kuramoto (1984) Y. Kuramoto, Chemical oscillations, waves, and turbulence, Vol. 19 (SpringerVerlag Berlin Heidelberg, 1984).
 Epstein and Pojman (1998) I. R. Epstein and J. A. Pojman, An Introduction to Nonlinear Chemical Dynamics : Oscillations, Waves, Patterns, and Chaos., Topics in Physical Chemistry (Oxford University Press, 1998).
 Pismen (2006) L. M. Pismen, Patterns and Interfaces in Dissipative Dynamics (SpringerVerlag Berlin Heidelberg, 2006).
 Maeda et al. (2008) Y. T. Maeda, J. Inose, M. Y. Matsuo, S. Iwaya, and M. Sano, PLoS ONE 3, e3734 (2008).
 Tanimoto and Sano (2014) H. Tanimoto and M. Sano, Biophysical journal 106, 16 (2014).
 Gray and Scott (1983) P. Gray and S. Scott, Chemical Engineering Science 38, 29 (1983).
 Gray and Scott (1984) P. Gray and S. Scott, Chemical Engineering Science 39, 1087 (1984).
 Koshizuka et al. (1998) S. Koshizuka, A. Nobe, and Y. Oka, International Journal for Numerical Methods in Fluids 26, 751 (1998).
 Schwarz and Safran (2013) U. S. Schwarz and S. A. Safran, Rev. Mod. Phys. 85, 1327 (2013).
 Barnhart et al. (2015) E. Barnhart, K.C. Lee, G. M. Allen, J. A. Theriot, and A. Mogilner, Proceedings of the National Academy of Sciences 112, 5045 (2015) .
 Iwamoto et al. (2014) M. Iwamoto, D. Ueyama, and R. Kobayashi, Journal of Theoretical Biology 353, 133 (2014).
 Style et al. (2014) R. W. Style, R. Boltyanskiy, G. K. German, C. Hyland, C. W. MacMinn, A. F. Mertz, L. A. Wilen, Y. Xu, and E. R. Dufresne, Soft Matter 10, 4047 (2014).
 Bosgraaf and Van Haastert (2009) L. Bosgraaf and P. J. M. Van Haastert, PLoS ONE 4, e5253 (2009).
 Bausch et al. (1999) A. R. Bausch, W. Möller, and E. Sackmann, Biophysical journal 76, 573 (1999).
 Micoulet et al. (2005) A. Micoulet, J. P. Spatz, and A. Ott, ChemPhysChem 6, 663 (2005).
 Kasza et al. (2007) K. E. Kasza, A. C. Rowat, J. Liu, T. E. Angelini, C. P. Brangwynne, G. H. Koenderink, and D. A. Weitz, Current Opinion in Cell Biology 19, 101 (2007) .
 Arai et al. (2010) Y. Arai, T. Shibata, S. Matsuoka, M. J. Sato, T. Yanagida, and M. Ueda, Proceedings of the National Academy of Sciences 107, 12399 (2010) .
 Shibata et al. (2012) T. Shibata, M. Nishikawa, S. Matsuoka, and M. Ueda, Journal of Cell Science 125, 5138 (2012) .
 Golebiewska et al. (2008) U. Golebiewska, M. Nyako, W. Woturski, I. Zaitseva, S. McLaughlin, and T. U. Martin, Molecular Biology of the Cell 19, 1663 (2008) .