Mechanochemical subcellular-element model of crawling cells

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 de-adhesion processes on intracellular biochemical signals. By introducing reaction-diffusion equations that reproduce traveling waves of local chemical concentrations, we clarify that the chemical dependence of the cell-substrate 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. Our present work sheds light on how intracellular chemical reactions are converted to a directional cell migration under the force-free condition. 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.


Introduction
Autonomous migration is a fundamental function of biological cells, and it is of essential importance in many biological processes during development and homeostasis. A number of studies have been conducted to reveal how intracellular biochemical reactions break the symmetry of a cell so that it migrates directionally Swaney et al. (2010); Beta and Kruse (2017); Devreotes et al. (2017). In order to bridge the gap between intracellular chemical signals and spatial migration of cells, we have to consider how the OPEN ACCESS EDITED BY force that a cell produces by itself is controlled by the intracellular chemical reactions and how it drives the cell.
From the viewpoint of mechanics, 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, i.e., vanishing simple sum of the force. This is because of the action-reaction law, which requires the existence of the counter force of internal force acting on another part of the active object with the same magnitude but in the opposite direction. Note that the internal force is generated inside the active object. Under this force-free condition, it is necessary to break symmetry to achieve spontaneous motion, such as directional motion. The same applies to the torque that active matter generates. That is, spontaneous rotation requires torque-free condition. For microorganisms that swim in a fluidic environment, the scallop theorem Purcell (1977) describes the importance of breaking reciprocality to achieve a net migration via internal cyclic motions. In particular, studies on a simple model of Purcell's three-bead swimmer revealed that the phase shift between the two periodically stretching bonds can break time reversal symmetry to realize a net migration instead of reciprocating motion Najafi and Golestanian (2004); Günther and Kruse (2008).
In contrast to the locomotion of microswimmers that stir surrounding fluid to propel, adhesion to the substrate such as the extracellular matrix and other cells plays an important role in the locomotion of cells crawling on the substrates, on which the cells exert traction force. Typical mechanism of such cell crawling includes the following four processes, 1) protrusion, 2) adhesion to the substrate, 3) de-adhesion from the substrate, and 4) contraction, as sketched in Figure 1. This simple mechanism of crawling cycle is currently accepted in the biology community Ananthakrishnan and Ehrlicher, (2007). In this case, the forcefree condition requires the existence of the counter force of the extensile and contractile force of cytoskeleton during the protrusion and contraction, that act on another part of the cell. See also the sketch in Figure 1. Note that the protrusion and contraction are induced by the dynamics of intracellular cytoskeleton Ananthakrishnan and Ehrlicher, (2007); Blanchoin et al. (2014); Tarama and Shibata (2022). In our previous study Tarama and Yamamoto (2018), we clarified the role of substrate adhesion by using a simple mechanical model in which a cell is described by two subcellular elements, that can switch their substrate friction periodically between the adhered stick state and the de-adhered slip state, and that are connected by a viscoelastic bond including an actuator that elongates and shrinks cyclically. By tuning the phase shifts between the actuator elongation and the substrate friction change of each subcellular element, we demonstrated that coupling between the cyclic intracellular force and the dynamic asymmetry of the substrate friction, has a great impact on the crawling distance and efficiency, as well as the crawling direction. By analogy to the Purcell's three-bead swimmer model, this model cell can achieve a net crawling motion by adjusting the phase shift among the periodic stretching of cell body and the adhesion/de-adhesion dynamics at the front and rear of the cell. Several similar studies Sketch of the crawling cycle of a cell on substrate, consisting of the four processes: ① protrusion, ② adhesion, ③ de-adhesion, and ④ contraction. The cell depicted in cyan enclosed by membrane (black line) is crawling on the substrate (gray). The two magenta circles schematically represent the front and rear parts of the cell, and the orange bars symbolize the interaction with the substrate underneath. Protrusion is led by the extensile force f e on the leading edge generated by intracellular cytoskeleton, and thus, its counter force − f e acts on another part (i.e., the rear in the sketch) of the cell because of the law of action and reaction. In the same way, the contraction process is caused by the contractile force f c at the trailing edge due to the force generation of intracellular cytoskeleton, and its counter force − f c act on the front of the cell. The local interaction of the cell with the substrate underneath (orange bars) can change its characteristics between the adhered state and the de-adhered state. The processes ② and ③ in the box on the right represent the recovery of the substrate adhesion to close the crawling cycle.
In contrast to such mechanistic tuning of the phase shift, in real cells, these extension and contraction, as well as the adhesion/de-adhesion processes, are controlled by intracellular chemical reactions. Therefore, in this paper, we address the question how intracellular chemical reactions can drive the spatial migration of a crawling cell by controling periodic extension and contraction of cytoskeleton under the force-and torque-free conditions as well as the substrate adhesion dynamics. To this end, we develop a basic particlebased mechanical model for cell crawling based on the typical mechanism of cell crawling sketched in Figure 1, which is coupled to intracellular chemical reactions. By considering a cell crawling on a flat substrate, we describe a cell by a set of many subcellular elements connected by viscoelastic bonds Newman (2007) as a simple extension of our previous mechanical model to two dimensions Tarama and Yamamoto (2018). In addition, intracellular chemical reactions are represented by simple reaction-diffusion (RD) equations Murray (2002Murray ( , 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.

Materials and methods
First, we introduce our mechanical model of a crawling cell and the RD equations representing intracellular chemical reactions. As for RD equations, we employ a previously introduced model for an intracellular chemical reaction observed experimentally. 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.

Subcellular-element model
We describe a single cell by a set of subcellular elements Newman (2007) connected by Kelvin-Voigt type viscoelastic bonds, as schematically depicted in Figure 2. The elastic spring and damper of the Kelvin-Voigt bonds represent intracellular elasticity and dissipation. In addition, a linear actuator is included to the Kelvin-Voigt bond, which models expansion and contraction due to intracellular force generation by cytoskeleton. Each subcellular element interacts with the substrate underneath via the substrate friction proportional to the local velocity, where the substrate friction coefficient can change its characteristics between adhered stick state and deadhered slip state. Since the typical size of a cell is on the order of 10 μm, the effect of inertia is negligible. Then, the force balance equation of element i is given by where v i is the velocity of the element i located at the position r i . Here, the abbreviations r = |r| andr r/r are used for the relative position r ij = r j − r i . The summation is over the set of the elements Ω i connected directly to the element i by the viscoelastic bonds. Note that in this paper we consider the case that the connection of the bonds between subcellular elements are permanent without breakup nor reconnection nor creation of new bonds for simplicity. The first term on the left-hand side of Eq. 1 represents the substrate friction with coefficient ζ i (t), which changes over time due to intracellular activity. The second term represents intracellular dissipation with the rate ξ. The first term on the right-hand side represents intracellular elasticity with the elastic modulus κ and free length ℓ ij . Intracellular activity is also Sketch of the subcellular element model of a cell crawling on a substrate. (A) The cell is described by a set of subcellular elements (magenta circles) connected by viscoelastic bonds (blue lines). The shape of a cell at rest is assumed to be a perfect hexagonal lattice. The element indicated by the star is the activator element. (B) Details of the subcellular elements and the connecting viscoelastic bond. Each element possesses the chemical concentrations c i . The actuator length ℓ act ij (t) and the substrate friction coefficient ζ i (t) change over time.
Frontiers in Cell and Developmental Biology frontiersin.org 03 included in the actuator, which tends to elongate the connecting bond by changing the free length over time as ℓ ij + ℓ act ij (t). Here, ℓ act ij (t) represents the actuator elongation, from which the force generated by the actuator is calculated as f act ij (t) −κℓ act ij (t)r ij /ℓ ij . We emphasize that the model Eq. 1 satisfies the force-free 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 ( where is the intracellular force acting on the element i. The last term on the right-hand side of Eq. 1 prevents the collapse of the subcellular element network. It is given by f area i −zU area /zr i , where U area 〈i,j,k〉 σ/S 2 ijk with σ = 10 -6 . This potential U area penalizes shrinking of the area of each triangle 〈i, j, k〉 formed by connected subcellular elements i, j, and k, which is defined by S ijk (r ij × r ik ) · e z /2 withê z as the unit vector perpendicular to the 2D substrate. We scale the system by L 0 = 10 μm for length and T 0 = 1 min 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 F 0 = 10 nN, 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 Table 1. We note that the dimensionless parameter constructed by the typical time scale, cell elasticity, and dissipation rate, that is related to the fluid-solid transition of cellular mechanics, is consistent with the typical value for a living cell summarized in Table 2.

Intracellular chemical reaction
In the model Eq. 1, the effects of the intracellular activities are included in the actuator elongation ℓ act ij (t) and the change in the substrate friction coefficient ζ i (t). The former represents the protrusion and contraction processes. The latter corresponds to the adhesion and de-adhesion 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.
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), PI3-kinase (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): where the reaction terms are defined as The global couplings are given by where the summation is over all subcellular elements. Here, U i and V i represent the PIP2 and PIP3 concentrations for subcellular element i, respectively. The first and second terms in G U represent the phosphorylation of PIP2 and the dephosphorylation of PIP3, respectively. The counterparts appear in the first two terms in G V . The third and fourth terms in G U are constant production and degradation of PIP2. The last term of G V represents a constant degradation of PIP3 Taniguchi et al. (2013). The global coupling terms originate from the conservation of the PI3K and PTEN concentrations.
To integrate Eq. 4, we calculate the Laplacian by using the moving particle semi-implicit (MPS) method Koshizuka et al. (1998). That is, for the chemical component where n ij = (n i + n j )/2, and n i = j≠i w (r ij ) is the number of neighboring elements around element i. The weight function that we employed is defined by where r e is the cutoff length, which is set to 4ℓ 0 . The normalization factor λ is given by r 2 e /6. By using this method,

Mechanical parameters Model Living cells References
Migration speed ≈ 10 μm(min) −1 Maeda et al. (2008) Traction stress ≈ 100 Pa Tanimoto and Sano (2014) Traction force 10 nN (estimated) Cell-substrate friction coefficient Elastic modulus κ cell κ/L 2 cell 10-100 Pa Bausch et al. (1999); Micoulet et al. (2005) Elastic constant Dissipation rate ξ cell = ξL bond Fluid/solid transition time ξ cell /k cell tens of sec Kasza et al. (2007)   Frontiers in Cell and Developmental Biology frontiersin.org 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 (δU, δV) = (−I excite , + I excite ) to Eq. 4 on the activator element, assuming that the phosphorylation of PIP2 to PIP3 is enhanced for this element. The parameters of Eq. 4 that we used in this paper are summarized in Table 3. Here, the diffusion coefficient D U = 0.48 corresponds to 0.8 μm 2 (sec) −1 , which is within the range of the experimentally-observed 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 U-V space, as plotted in Figure 3.

Model parameters Values used in simulations
The important property of the RD equations, Eq. 4, is that they are of the Grey-Scott type Scott (1983, 1984). One of the advantages of the Grey-Scott model is that it can show either an excitable or a bistable nature depending on the parameters. Series of studies showed that the signaling pathway around PIP2-PIP3 reactions is excitable Nishikawa et al. (2014) 2014), which is also claimed in Taniguchi et al. (2013). Interestingly, similar RD equations were studied by Shao et al. Shao et al. (2010) in the context of cell crawling, where they assumed a bistable regime to reproduce the steady migration of keratocyte cells.

Mechanochemical coupling
To combine the cell mechanics, Eq. 1, and the RD equations, Eq. 4, we consider the coupling of the chemical concentrations to the actuator elongation ℓ act ij (t) and the substrate friction coefficient ζ i (t) individually.

Actuator elongation
First, we introduce the coupling between the RD equations and the actuator elongation. Higher concentration of PIP3 enhances actin intensity Taniguchi et al. (2013), which we assume leads to actuator elongation. Then, we introduce the dependence of the actuator elongation on the PIP3 concentration, such that it elongates with PIP3 concentration as where V ij (t) = (V i (t) + V j (t))/2 is the mean PIP3 concentration for the bond connecting the elements i and j. Although V i (t) is a positive quantity, its maximal value depends on the strength of the initial fluctuation because of the excitable nature of Eq. 4. Therefore, tanh is introduced on the right hand side of Eq. 9 to prevent an extremely large elongation. a is a constant denoting sensitivity, and ℓ V is the magnitude of the elongation. Here, we set a = π and ℓ V = ℓ ij .

Substrate adhesion
Next, we consider the adhesion to the substrate underneath and the de-adhesion from it. We model the adhesion/deadhesion processes by the transition of the substrate friction coefficient between the adhered stick state and the de-adhered slip state. Here, we consider the dependence of the substrate friction coefficient on both mechanical and chemical signals (see Figure 4A): where τ ζ is the time delay. A v takes a value between 0 and 1, representing the ratio of the mechanical and chemical dependence of the stick-slip transition of the substrate friction.
Here, we assume a rapid transition between the adhered stick state ζ stick and de-adhered slip state ζ slip , where the transition  Table 1.
Frontiers in Cell and Developmental Biology frontiersin.org sharpness is given by ϵ ζ . Here, we set as ϵ ζ = ζ slip /2. Note that the friction coefficient is smaller at the slip state than at the stick state, ζ slip < ζ stick . 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 h v (v) in Eq. 10 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 where v* is the threshold value. The subcellular element tends to adhere to the substrate (ζ = ζ stick ) if the speed is smaller than the threshold value, i.e., v < v*, while the element slips on the We set the threshold value to v* = 1. The formula of the stick-slip transition of the cell-substrate friction depending on the local velocity, i.e., Eq. 10 with A v = 1, was introduced in Ref. Barnhart et al. (2015). As a result of the balance of the two functions, h ζ (ζ) and h v (v), 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: In Eq. 13, σ V stands for the sensitivity, and V* is the threshold concentration. Due to this term, a large value of V prevents strong adhesion if σ V > 0, while large V enhances the adhesion if σ V < 0. However, we are not sure whether PIP3 enhances or diminishes adhesion. Therefore, in the next section, we numerically solve the time-evolution equations for both cases and see what will happen. In Figures 4B,C, we show two example trajectories of ζ that changes depending only on the local velocity v ( Figure 4B for A v = 1) and on the PIP3 concentration V ( Figure 4C for We note that the mechanical dependence of the stick-slip transition of the substrate friction is consistent with that of usual materials such as adhesive plastic tapes, which tends to de-adhere under strong mechanical force or velocity. In contrast, the chemical cue seems rather characteristic to biological cells, although to our knowledge the signaling pathway that regulates the substrate adhesion mediated by the focal adhesion is yet to be clarified.

Numerical analysis
To integrate the time-evolution Equations 1, 4 with Eqs 9, 10 numerically, we employ the Euler method with time increment δt = 10 -5 . We add the stimulus (δU, δV) = (−I excite , + I excite ) to Eq. 4 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 Figure 2A, except in the case of random motion, for which an activator element is chosen randomly every time t = 0.15. Since the traveling wave is generated by the stimuli on the activator element, we excite the activator element every t = 1.5 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 I excite = 0.75.
We note that most of the analyses in the following sections are performed by using a hexagonal cell shape, depicted in Figure 2. This is because it enables us to prepare a regular lattice of subcellular elements and, thus, is free from the complexity that the inhomogeneity of the position and connection of the subcellular elements may raise. In order to see the impact of this choice, we compare the results with the case of a circular cell shape in Sect. 3.5, which looks more relevant to realistic cells.

Inhomogeneity of substrate adhesion
We start with considering the role of the inhomogeneity of the substrate adhesion on cell crawling. From Eq. 1, the centreof-mass velocity is calculated as where f int i is defined by Eq. 3. If the substrate adhesion is homogeneous, ζ i (t) = ζ 0 (t) ≠ 0, Eq. 14 becomes Here, the force-free condition, Eq. 2, is used at the second equality. Note that the time dependence of the homogeneous substrate adhesion does not lead to a finite centre-of-mass velocity. Note also that the summation in Eq. 15 can be replaced by the integral over the cell in the continuum limit without loosing generality. This demonstrates that spatial translational motion is not achieved without symmetry breaking of substrate adhesion under the forcefree condition.
Frontiers in Cell and Developmental Biology frontiersin.org 07

Sinusoidal traveling chemical wave
The analysis of the two-element case, i.e., dumbbell model in our previous paper Tarama and Yamamoto (2018) showed that the coupling between the asymmetric substrate friction and the actuator elongation affects the crawling efficiency. In fact, there exists a reciprocating motion where a cell propels in one direction and moves in the opposite direction for the same distance for the rest of the period, resulting in no net migration. In this section, we test if this coupling also changes the crawling motion of a cell composed of many subcellular elements.
As a simple realization of inhomogeneous substrate adhesion, we consider a traveling harmonic wave of a single intracellular chemical component c(t): where the phase changes as Here c 0 is the maximum concentration, which is set to c 0 = 0.5 so that 0 ≤ c i (t) ≤ 1. ω and q w are the frequency and the wavenumber of the traveling wave, and x* is the x position of the activator element, from which the traveling wave occurs. The dependence of the actuator elongation on the intracellular chemical signal Eq. 9 is replaced by where ℓ c represents the magnitude of the elongation. The dynamics of the substrate friction coefficient ζ i (t) is also replaced by a simple two-state function that switches between the adhered stick state and the de-adhered slip state: where m i is an integer that satisfies 2m i π < ϕ i (t) − ψ ζ ≤ 2 (m i + 1)π. ψ ζ is a phase shift between the asymmetric substrate friction and the actuator elongation. Since these internal motions, i.e., the actuator elongation and the stick-slip transition of substrate friction, are perfectly cyclic, this phase shift controls the coupling between the asymmetric substrate friction and the actuator elongation. We carried out numerical simulations of Eq. 1 together with Eqs 16-19 with the time increment δt = 10 -4 . In the simulation, we choose the element indicated by the star in Figure 2A as the activator element, and we set q w > 0 so that the wave travels from Frontiers in Cell and Developmental Biology frontiersin.org right to left with the frequency ω = 2π. We let the actuator elongate twice as long as its equivalent length: ℓ c = ℓ ij . For these parameters, we vary the wavenumber q w and the phase shift ψ ζ and measure the migration distance for one cycle. Figure 5A shows that the migration distance for one cycle ΔR alters depending on the phase shift ψ ζ . The sign of ΔR distinguishes the migration direction. Positive ΔR represents rightward motion, whereas negative ΔR corresponds to leftward motion. Since the wave of the intracellular chemical concentration c i travels from rightmost subcellular element toward the left of the cell, the direction of the rightward cell migration is opposite to that of the traveling wave, as depicted in Figure 5B. Note that the cell shape at time t = 3T/4 looks like a flipped version of the cell at t = 0 in Figure 5B. This is because the phase of the intracellular chemical concentration in Eq. 17 is defined as a function of the distance with respect to the rightmost subcellular element and the length of the bonds connecting subcellular elements changes due to the actuator elongation. On the other hand, in the case of the leftward migration, the cell migrates in the same direction as that of the intracellular traveling wave of c i as shown in Figure 5C. At the phase shift ψ ζ where the migration direction switches from the rightward migration to the leftward motion, reciprocating motion appears, as depicted in Figure 5D. In reciprocating motion, the cell migrates in one direction (to the left in the case of Figure 5D) during some time of the period, and migrates in the other direction (to the right in Figure 5D) for the same distance during the rest of the period. Therefore, the reciprocating motion results in no net migration, and thus, ΔR = 0. Note that the asymmetry of substrate friction exists even in the cell undergoing the reciprocating motion. The migration distance ΔR also depends on the wavenumber q w . There is an optimal q w around q w = 4π/5 for which the cell can achieve the largest migration distance.
These results highlight that in addition to the substrate friction asymmetry, the spatial motion of a cell changes depending on the phase shift between the force generation and the formation of asymmetric adhesion, which corresponds to the time delay in Eq. 10 where the intracellular chemical reaction is controlled by reaction-diffusion equations.

Crawling by direct and retrograde waves
Now we consider the situation where the cell mechanics Eq. 1 is driven by intracellular chemical reaction Eq. 4. First, we study the effect of the sign of σ V , by numerically integrating Eqs 1, 4

FIGURE 6
Cell crawling obtained from Eqs 1, 4 with Eqs 9, 10 for different signs of σ V : (A) σ V = 2π and (B) σ V = −2π. The cell for positive σ V crawls in the opposite direction against the traveling chemical wave as shown in panel (A), whereas, for negative σ V , it moves in the same direction as the wave as displayed in panel (B). The other parameters are set to A v = 0 and τ ζ = 0.01. In both panels, the first column depicts two-dimensional snapshots, and the second column shows the corresponding values of V i (by dots with color that also represents the value of V i , connected by black line) and ζ i (in gray) along the cross-section at y = 0. The position of each subcellular element is plotted by a circle whose size and color indicate the value of ζ i and V i , respectively. The color of the connecting bonds corresponds to V ij . The number in the bottom left corner of each subplot represents the time.
Frontiers in Cell and Developmental Biology frontiersin.org with Eqs 9, 10 for both positive and negative σ V . Figure 6 depicts a time series of snapshots of a crawling cell for σ V = 2π in Figure 6A and σ V = −2π in Figure 6B, respectively. In both Figures 6A,B, the first column depicts two-dimensional snapshots, and the second column shows the corresponding values of V i (by dots with color that also represents the value of V i , connected by black line) and ζ i (in gray) along the cross-section at y = 0. We set the threshold in Eq. 13 to V* = 0.5 throughout this paper. If σ V is positive, the cell moves in the opposite direction to the PIP3 traveling wave, as shown in Figure 6A. Interestingly, however, if the sign of σ V is negative, a qualitatively different result appears: namely, the cell starts to move in the same direction as the traveling wave, as displayed in Figure 6B. This difference in the direction of migration can be understood from the time-series of the values of V i and ζ i along the cross-section at y = 0 in Figure 6. In the case of positive σ V , the subcellular elements tend to be at the de-adhered slip state ζ slip for high values of V i , when the actuator tries to elongate. This leads to the protrusion of cell edge around the origin of the traveling PIP3 wave, i.e., on the right side of the cell at the earlier time period (0 ≲ t ≲ 0.24), which is then adhered strongly to the substrate for the later time period (0.36 ≲ t ≲ 0.6) when V i in this region has already returned close to the resting state (V i~0 ). Therefore, in this case, protrusion around the origin of the wave leads to the cell migration in the opposite direction to the traveling wave. On the other hand, for negative σ V , the subcellular elements tend to adhere strongly to the substrate for higher V i . This makes the actuator elongation at the earlier time t~0.12 result in an almost symmetric deformation of the cell along the x axis at y = 0. At the later time 0.36 ≲ t ≲ 0.48, however the recovery of the PIP3 concentration V i to the resting state makes the subcellular elements de-adhere from the substrate around the back of the PIP3 wave, which causes the actuator contraction to bring the right half of the cell towards the PIP3 wave traveling at the left half of the cell where the subcellular elements adheres strongly to the substrate. This effective contraction drives the cell in the same direction as the traveling wave in the case of negative σ V . Note that the actuator elongates for high values of the PIP3 concentration V i .
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. In this sense, the above crawling motion for positive σ V in Figure 6A corresponds to the motion with the retrograde wave, and the one for σ V < 0 in Figure 6B corresponds to the motion with the direct wave. We note that the crawling by the direct and retrograde waves are also reported in the model for adhesive locomotion of gastropods Iwamoto et al. (2014). 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 σ V = 2π for the rest of this paper.

Mechanical vs chemical control of adhesion
Now, we study the effect of the mechanical and chemical dependence of substrate friction coefficient on cell crawling. In particular, we focus on the following two factors. One is the time delay τ ζ controlling the phase shift of the change in the substrate friction characteristics with respect to the intracellular chemical wave and, thus, to the intracellular force generation due to the actuator elongation, which is also induced by the chemical wave. The other is the parameter A v , which gives the dependence ratio of the substrate friction ζ on the local velocity and on the intracellular chemical signal. A v is varied between A v = 0, corresponding to the case where the substrate friction is fully controlled by the intracellular chemical signals, and A v = 1, where the characteristics of the substrate adhesion changes depending only on the mechanical load. To characterize the cell migration, we measure the migration distance ΔR of the cell in one cycle, i.e., with one traveling wave.
In Figure 7A, we show the dependence on the time delay for different values of A v as indicated in the legend. In all cases of A v , the dependence on the time delay appears strongest for large τ ζ . In this case, the response of the adhesion characteristics to the mechanical and chemical signals is so slow that the actuator elongation induced by the intracellular chemical wave of V is not well converted to a spatial migration of the cell. For a small τ ζ , the impact of the time delay becomes small. However, in the case of purely mechanically controlled substrate friction A v = 1, the migration distance decreases for small τ ζ , and the migration is most efficient for the intermediate time delay around τ ζ = 0.2. In addition, the migration distance ΔR is always much worse in the case of A V = 1 than in the case of A v = 0. When A v = 1, the substrate friction characteristics changes depending only on the local speed and tends to be at the de-adhered slip state for high speed. Since this dependence does not distinguish the expanding and contracting processes, the substrate friction becomes the de-adhered slip state in both processes. Therefore, although the actuator elongation drives the cell front forward, the following actuator contraction moves it backward, resulting in a worse migration of the cell.
Interestingly, the mixing of the mechanical and chemical dependences of the substrate friction may result in larger values of ΔR than purely mechanical or chemical control. In order to highlight this point, we plot ΔR against A v for τ ζ = 0.01 and 0.1 in Figure 7D. The plot in Figure 7D clearly shows the maximum ΔR at the intermediate A v~0 .6. In Figures 7A,D, ΔR/L cell reaches approximately 0.36 for A v = 0.6 and τ ζ = 0.01. This value is comparable to the measurement of a crawling Dictyostelium cell, which moves its body length in approximately two cycles Tanimoto and Sano (2014).
Frontiers in Cell and Developmental Biology frontiersin.org

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, where the subcellular elements are connected as depicted in Figure 8A. The length of the cell was set to L cell = 1, and the total number of subcellular elements was set to N = 100. The total number of bonds connecting the subcellular elements was then N bond = 267. We set ℓ ij of each bond of the circular cell to the initial element distance, as shown in Figure 8A. Then, the mean bond length was calculated to be ℓ 0 = 0.105724. The other parameters were kept the same as for hexagonal cells. Figure 8B shows an example of a time series of circular cell crawling obtained numerically from Eqs 1, 4 with Eqs 9, 10. Here, the activator element is set as the one denoted by the star in Figure 8A, which is then stimulated by (δU, δV) = (−I excite , + I excite ) with I excite = 0.75 every t = 1.5. We then measure the migration distance for different values of A v and τ ζ . The results are plotted in Figures 7B,E, which are qualitatively the same as those of the hexagonal cell in Figures 7A,D. This means that the rest shape of the model cell has less impact on the crawling dynamics.  Frontiers in Cell and Developmental Biology frontiersin.org 3.6 Impact of cell size Now, we study the impact of the size of the cell. We prepare a hexagonal cell of size L cell = 1.4, which is approximately twice the area of the previous hexagonal cells. We again measure the migration distance ΔR, which is normalized by the cell length L cell to facilitate comparison with the previous results for L cell = 1.
The results are summarized in Figures 7C,F. Qualitatively, the tendency is the same as that of the results in Figures 7A,D for L cell = 1. 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.

Random excitation
In reality, cells change their migration direction over time. In our model, we can reproduce such motion by introducing intracellular stochasticity, which may originate from, e.g., the complexity of intracellular chemical reaction processes. Here, we randomly choose one element in every t = 0.15 and add to that element the stimulus (δU, δV) = (−I excite , + I excite ) with intensity I excite = 0.75.
The results are summarized in Figure 9. In Figures 9A,B, the parameter values are kept the same as in Figure 6A. Due to the stochasticity, however, the cell changes its migration direction frequently, as shown in the trajectory of the center-of-mass position in Figure 9A. From the snapshots of the cell in Figure 9B, we 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, in principle, a new wave tends to originate from the elements that are near the origin of the previous wave.
To quantify the persistence of the migration direction, we measure the mean-square displacement (MSD), which is defined by 〈(X cm (t 0 + t) − X cm (t 0 )) 2 〉, where X cm is the center-ofmass position of the cell and 〈x (t 0 )〉 represents the time average of x (t 0 ), i.e., the average of x (t 0 ) over time t 0 . The results are displayed in Figure 9E. It shows a crossover from ballistic (∝ t 2 ) to diffusive regimes (∝ t 1 ) at around t~1. In fact, the velocity autocorrelation function (VAC) defined by 〈Ṽ cm (t 0 + t) ·Ṽ cm (t 0 )〉 for the direction of the center-of-mass velocityṼ cm V cm /|V cm |, almost vanishes at this time scale, as shown in Figure 9F. If we compare this result with the experimental measurement in Ref. Li et al. (2008), which  Table 3.
Frontiers in Cell and Developmental Biology frontiersin.org reported a persistent time of t p = 8.8 ± 0.1 min, the value that our model reproduced is slightly smaller. We think this originates from the fact that our model cell lacks explicit polarity that memorizes the direction of migration. If the parameter K K in the reaction term of Eq. 5 is slightly increased from 5 to 5.5, the cell changes its migration direction more frequently, as depicted in Figure 9C. Depending on the random stimuli, the cell switches from directional motion to spinning motion as a spiral wave appears, as shown in Figure 9D. 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. In this case, the MSD shows a subdiffusive feature, where the MSD increases with the slop slightly smaller than t 1 for long time intervals, as shown in Figure 9E. This is probably because of the existence of the spinning motion that occurs from time to time. The spinning motion also leads to an oscillatory damping behaviour of the VAC, as plotted in Figure 9E.

Traction force multipoles
In many experiments, traction force that crawling cells exert on the substrate is measured Style et al. (2014) because of its fundamental importance for cell motility. The traction force results from the interaction between the cell and the substrate: f traction i ζ i (t)v i . Since we are interested in cellular scale dynamics, we calculate the traction force multipoles, which was also measured experimentally Tanimoto and Sano (2014).
First, the traction force monopole is defined by  Figure 9A, whereas those in panels (B,D) are for the cell in Figure 9B. Note that the axis of the M (2) 11 is inverted, to match the plot in Ref. Tanimoto and Sano (2014) in panels (C,D).
Frontiers in Cell and Developmental Biology frontiersin.org where i represents the subcellular element i and the summation runs over the entire cell. The subscript α indicates spatial component α = 1, 2, corresponding to x and y coordinates. Here, to compare with the experimental results in Ref. Tanimoto and Sano (2014), the traction force multipoles are calculated in the comoving coordinate on the cell, and thus, α = 1 and 2 represent the components parallel and perpendicular to the instantaneous center-of-mass velocity, respectively. In the numerical simulation of the model crawling cell, the traction force monopole is equal to 0, as shown in Figure 10, which is consistent with the experimental result Tanimoto and Sano (2014). This is readily understood from the force balance equation, Eq. 1, and the force-free condition, Eq. 2. That is, the traction force monopole vanishes because of the fact that the inertia is negligibly small for crawling cells on top of the force-free condition. Given that the force-free condition is a requirement for the intracellular force generation, we can understand that the vanishing force monopole obtained from the analysis of the experiment Tanimoto and Sano (2014) indicates the fact that the inertia is negligibly small for the crawling cell. The next lowest mode is the traction force dipole, the element of which is defined by On the one hand, from Figures 10A,B, 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 dipoles, respectively. 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 off-diagonal components M (2) 12 and M (2) 21 take the same values, indicating that the traction force dipole is symmetric, although such symmetry is not presumed in its definition, Eq. 21. 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 torque-free nature of the cell. Here, the torque is defined by which, in a two dimensional space, becomes by using Eq. 21. Note that, if one describes the force dipole tensor with its invariants as in the fourth panels in Figures 10A,B, the torque appears as the imaginary part of the eigenvalues. Here, the trajectories of the two eigenvalues λ 1 and λ 2 of the traction force dipole look identical to those of M (2) 11 and M (2) 22 , in particular in Figure 10A. This is because the direction θ 1 of the eigenvector of M (2) is basically in the same direction as the center-of-mass velocity (φ v ) as shown in the last panel of Figure 10A. Here the eigenvalue λ 1 is selected such that the direction of the corresponding eigenvector (θ 1 ) is closer to φ v than that of the other eigenvalue λ 2 . We emphasize that the torque-free property of the traction force is satisfied even by the spinning cell, as plotted in Figure 10B. In this case, the deviation of the angle θ 1 of the eigenvector corresponding to λ 1 from the angle φ v of the center-of-mass velocity is larger as plotted in Figure 10B, resulting in the larger deviation of the trajectories of λ 1 and λ 2 from those of M (2) 11 and M (2) 22 . Finally, in Figure 10C, we plot the time evolution of the traction force dipole M (2) 11 and the traction force quadrupole M (3) 111 , which is compared with the measurement in the experiment Tanimoto and Sano (2014). Here the traction force quadrupole is defined by Interestingly, the trajectory in M (2) 11 -M (3) 111 space shows a counterclockwise rotation, which is qualitatively consistent with the experimental results Tanimoto and Sano (2014). For the spinning cell, this counter-clockwise rotation is not clear in Figure 10D because it shows more rapid change in the magnitude of the traction force multipoles. Note that in the corresponding plot in Ref. Tanimoto and Sano (2014) the axis of the force dipole is inverted to highlight the magnitude of contraction. We followed this in Figures 10C,D where the axis of M (2) 11 is inverted, to make the comparison easy.

Summary and discussion
To summarize, we investigate how intracellular chemical reactions can drive the spatial migration of a crawling cell by controling periodic extension and contraction of cytoskeleton under the force-and torque-free conditions as well as the substrate adhesion dynamics. To this end, we constructed a mechanochemical model of a cell crawling on a substrate based on the typical mechanism of cell crawling sketched in Figure 1. The mechanical part is described by a subcellularelement model, where we extend our previous model Tarama and Yamamoto (2018), and the chemical part is described by RD equations proposed by Taniguchi et al. (2013). To combine them, Frontiers in Cell and Developmental Biology frontiersin.org we introduce two mechanical activities. One is the actuator elongation, which depends on the intracellular chemical concentration. The other is the substrate friction, which shows a sharp transition between the adhered stick state and the deadhered slip state depending on the intracellular chemical concentration in addition to the local velocity. We also introduce a time delay of the substrate friction change.
Although we assumed a simple function form of the substrate friction dynamics and its dependence on the local velocity and intracellular chemical concentration, we believe that this time delay can be measured experimentally. By using this model, we demonstrated that the substrate adhesion dynamics affect how the intracellular force leads to crawling motion under the force-and torque-free conditions. Both the substrate adhesion and intracellular force generation are controlled by the traveling wave of intracellular chemical reactions on the cellular scale. Depending on the sign of the sensitivity of the substrate friction coefficient to the intracellular chemical concentration, the model cell exhibited crawling with the retrograde flow or with the direct flow. In the former case, 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, by introducing stochasticity in the RD equations, the cell changes its migration direction and even switches 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.
The persistence time of the cell crawling in our current model is slightly smaller than the experimental observation Li et al. (2008). The origin of the persistence in our current model is the excitable nature of the RD equations. That is, the intracellular chemical wave is triggered by a stimulus on the elements that is close to the resting state. Therefore, a new wave tends to occur from the region close to the origin of the previous wave, where the chemical signals have more chance to relax back to the resting state than the other part of the cell. This explains why our model cell cannot be persistent over too many wave cycles. Real cells, however, often establish a rather stable polarity, which allows the cells to migrate more persistently. This polarity can be modelled by bistable chemical reactions. In fact, recent experimental studies shed light on the coupling of the excitable chemical signaling pathways and the bistable chemical reactions, representing polarity Fukushima et al. (2019);Flemming et al. (2020). Therefore, it is important to include polarity to the model and study how it changes the relation among the intracellular chemical reactions and the force-and torque-free intracellular mechanics and the spatial migration of cells to realize a persistent crawling migration. We will come to this point in near future. Finally, we discuss other possible extensions of our current model.
Contraction process: In our current model, the protrusion and contraction processes are both modeled by the actuator elongation, ℓ act ij (t), of the bond connecting two subcellular elements. The two processes are distinguished by the sign of ℓ act ij (t). In this paper, however, we consider only the protrusion process, i.e., ℓ act ij (t) ≥ 0, 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 should extend our current model to include more detailed dynamics of the concentration of adhesion molecules and their diffusion to other subcellular elements. Then, we can discuss realistic structures such as the footstep-like 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 Figure 6A. 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 V i . Therefore, this inconsistency may be resolved by, e.g., introducing dependence on the gradient of V i .In addition, the network of the subcellular elements is permanent in our current model. However, it is interesting to see what will happen by including the reconnection of the subcellular elements that evolves with the dynamics of the subcellular elements. With this extension, we expect more complicated deformation of the cell shape, which may look more realistic.
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.
In conclusion, the modeling of crawling cells is still a challenging task because of the complexity in intracellular processes and of the variety of the cell crawling mechanism. By focusing on the mechanism how intracellular chemical reactions control the periodic extension and contraction of intracellular cytoskeleton to achieve a net directional migration under the force-and torque-free conditions (Figure 1), we demonstrated that the symmetry breaking Frontiers in Cell and Developmental Biology frontiersin.org due to inhomogeneous substrate adhesion and its time delay are of fundamental importance. In real cells, there are a huge variety of intracellular chemical reactions, and often they are related to each other. Nevertheless, we believe that this study can provide a basic framework to understand the synergetic mechanism of intracellular chemical reactions to realize effective migration of complex real cells by controlling the extension and contraction of cytoskeleton under the forceand torque-free conditions.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
MT and RY designed the research; MT derived the model; MT and KM conducted analyses; MT and RY wrote the manuscript.

Funding
This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (17H01083, 19K14673, 21KK0127, 22K14017, 22K06214) grant, as well as the JSPS bilateral joint research projects.