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.

1 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 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 force-free 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 of simple inchwarm-style model have also shown that the regulation of substrate adhesion plays an important role in directional motion of cells on a substrate Kumar et al. (2008); Lopez et al. (2014); Mai and Camley (2020). We also note that many elaborate models have been introduced to explain various aspects of the cell crawling dynamics, including 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 particle-based model Newman (2007); Sandersius and Newman (2008); Basan et al. (2011); Zimmermann et al. (2016); Smeets et al. (2016).

FIGURE 1

FIGURE 1. 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 particle-based 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 (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.

2 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.

2.1 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 de-adhered 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

$\begin{array}{ll}\hfill & {\zeta }_i\left(t\right){\mathit{v}}_i+\sum _{j\in {\mathrm{\Omega }}_i}\xi {\ell }_{ij}\left({\mathit{v}}_i-{\mathit{v}}_j\right)\hfill \\ \hfill & =\sum _{j\in {\mathrm{\Omega }}_i}\frac{\kappa }{{\ell }_{ij}}{\widehat{\mathit{r}}}_{ij}\left\{r_{ij}-\left({\ell }_{ij}+{\ell }_{ij}^{\text{act}}\left(t\right)\right)\right\}+{\mathit{f}}_i^{\text{area}},\hfill \end{array}(1)$

where v_i is the velocity of the element i located at the position r_i. Here, the abbreviations r = |r| and $\widehat{\mathit{r}}=\mathit{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.

FIGURE 2

FIGURE 2. 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 ${\ell }_{ij}^{\text{act}}(t)$ and the substrate friction coefficient ζ_i(t) change over time.

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 included in the actuator, which tends to elongate the connecting bond by changing the free length over time as ${\ell }_{ij}+{\ell }_{ij}^{\text{act}}(t)$. Here, ${\ell }_{ij}^{\text{act}}(t)$ represents the actuator elongation, from which the force generated by the actuator is calculated as ${\mathit{f}}_{ij}^{\text{act}}(t)=-\kappa {\ell }_{ij}^{\text{act}}(t){\widehat{\mathit{r}}}_{ij}/{\ell }_{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

$\sum _i{\mathit{f}}_i^{\text{int}}=0.(2)$

where

$\begin{array}{ll}\hfill {\mathit{f}}_i^{\text{int}}& =-\sum _{j\in {\mathrm{\Omega }}_i}\xi {\ell }_{ij}\left({\mathit{v}}_i-{\mathit{v}}_j\right)\hfill \\ \hfill & +\sum _{j\in {\mathrm{\Omega }}_i}\frac{\kappa }{{\ell }_{ij}}{\widehat{\mathit{r}}}_{ij}\left\{r_{ij}-\left({\ell }_{ij}+{\ell }_{ij}^{\text{act}}\left(t\right)\right)\right\}+{\mathit{f}}_i^{\text{area}}\hfill \end{array}(3)$

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 ${\mathit{f}}_i^{\text{area}}=-\partial U^{\text{area}}/\partial {\mathit{r}}_i$, where $U^{\text{area}}={\sum }_{⟨i,j,k⟩}\sigma /S_{\mathit{ijk}}^2$ 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_{\mathit{ijk}}=({\mathit{r}}_{ij}\times {\mathit{r}}_{ik})\cdot {\widehat{\mathit{e}}}_z/2$ with ${\widehat{\mathit{e}}}_z$ as the unit vector perpendicular to the 2D substrate.

We scale the system by L₀ = 10 μm for length and T₀ = 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₀ = 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.

TABLE 1

TABLE 1. Mechanical parameters in the model and their values for the simulation of the hexagonal cell and the corresponding values in physical units. See also Table 2. The values in parentheses are for a large cell.

TABLE 2

TABLE 2. Mechanical parameter values for typical living cells and corresponding model parameters.

2.2 Intracellular chemical reaction

In the model Eq. 1, the effects of the intracellular activities are included in the actuator elongation ${\ell }_{ij}^{\text{act}}(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):

$\begin{array}{l}\hfill \frac{\partial U_i}{\partial t}=D_U{\nabla }^2{U}_i+G_U\left(U_i,V_i\right),\\ \hfill \frac{\partial V_i}{\partial t}=D_V{\nabla }^2{V}_i+G_V\left(U_i,V_i\right),\end{array}(4)$

where the reaction terms are defined as

$\begin{array}{ll}\hfill & G_U\left(U,V\right)=-\frac{\alpha U{V}^2}{K_K+\langle V^2\rangle }+\frac{\beta UV}{K_P+\langle U\rangle }+S-\gamma U,\hfill \\ \hfill & G_V\left(U,V\right)=+\frac{\alpha U{V}^2}{K_K+\langle V^2\rangle }-\frac{\beta UV}{K_P+\langle U\rangle }-\mu V.\hfill \end{array}(5)$

The global couplings are given by

$\langle U\rangle =\frac{1}{N}\sum _i{U}_i,\langle V^2\rangle =\frac{1}{N}\sum _i{V}_i^2,(6)$

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 c_i = {U_i, V_i}, the diffusion term is modeled by

$D_c{\nabla }^2{c}_i=\frac{4D_c}{\lambda }\sum _{j\ne i}\frac{1}{n_{ij}}\left(c_j-c_i\right)w\left(r_{ij}\right),(7)$

where n_ij = (n_i + n_j)/2, and n_i = ∑_j≠iw (r_ij) is the number of neighboring elements around element i. The weight function that we employed is defined by

$w\left(r\right)=\left\{\begin{array}{cc}\frac{r_e}{r}-1\hfill & \mathrm{f}\mathrm{o}\mathrm{r}0\le r<r_e\hfill \\ 0\hfill & \mathrm{f}\mathrm{o}\mathrm{r}{r}_e\le r\hfill \end{array}\right.(8)$

where r_e is the cutoff length, which is set to 4ℓ₀. The normalization factor λ is given by $r_e^2/6$. 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 (δ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² (sec)⁻¹, 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.

TABLE 3

TABLE 3. The parameter values of the RD equations in the simulation.

FIGURE 3

FIGURE 3. The nullcline and the flow field in the U-V space of Eqs 4, 5 under the assumption of uniform U(r) = U and V(r) = V. The orange solid line and the red dashed line correspond to the nullclines of G_U(U, V) = 0 and G_V(U, V) = 0, respectively, which correspond to dU/dV = 0 and dV/dt = 0. The gray lines show example trajectories starting from various (U, V), whereas the arrow heads and their color show the direction and magnitude of the flow (G_V(U, V), G_U(U, V)) at each phase point. The black dot represents the stable fixed point located at (V, U) = (0, S/γ).

The important property of the RD equations, Eq. 4, is that they are of the Grey-Scott type Gray and 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); Xiong et al. (2010); Huang et al. (2013); Devreotes et al. (2017); Gerisch et al. (2011, 2012); Gerhardt et al. (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.

2.3 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 ${\ell }_{ij}^{\text{act}}(t)$ and the substrate friction coefficient ζ_i(t) individually.

2.3.1 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

${\ell }_{ij}^{\text{act}}\left(t\right)={\ell }_V\tanh \left[a{V}_{ij}\left(t\right)\right],(9)$

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.

2.3.2 Substrate adhesion

Next, we consider the adhesion to the substrate underneath and the de-adhesion from it. We model the adhesion/de-adhesion 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):

${\tau }_{\zeta }\frac{d{\zeta }_i}{dt}=h_{\zeta }\left({\zeta }_i\right)-A_v{h}_v\left(v_i\right)-\left(1-A_v\right)h_V\left(V_i\right),(10)$

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.

FIGURE 4

FIGURE 4. (A) Functional forms of h_ζ(ζ), h_v(v), and h_V(V), which determine the substrate friction. Example trajectory of ζ that changes in time depending (B) only on the velocity v (A_v = 1) and (C) only on the PIP3 concentration V (A_v = 0). The parameters are set to the values given in Table 1.

The function h_ζ(ζ) is defined by

$h_{\zeta }\left(\zeta \right)=-\frac{1}{2}\tanh \left(\frac{\zeta -{\zeta }_{\text{stick}}}{{ϵ}_{\zeta }}\right)-\frac{1}{2}\tanh \left(\frac{\zeta -{\zeta }_{\text{slip}}}{{ϵ}_{\zeta }}\right).(11)$

Here, we assume a rapid transition between the adhered stick state ζ_stick and de-adhered slip state ζ_slip, where the transition 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

$h_v\left(v\right)=\frac{{\left(v/v^{*}\right)}^2}{1+{\left(v/v^{*}\right)}^2}-\frac{1}{2},(12)$

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 substrate (ζ = ζ_slip) if v > v*. 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:

$h_V\left(V\right)=\frac{1}{2}\tanh \left({\sigma }_V\left(V-V^{*}\right)\right).(13)$

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 A_v = 0). Since both h_v(v) and h_V(V) take values between [ − 0.5, 0.5] and since 0 ≤ A_v ≤ 1, −0.5 ≤ A_vh_v(v) + (1 − A_v)h_V(V) ≤ 0.5. Then, from Eq. 10, ζ takes a value between [ζ_slip, ζ_stick].

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.

2.4 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.

3 Results

3.1 Inhomogeneity of substrate adhesion

We start with considering the role of the inhomogeneity of the substrate adhesion on cell crawling. From Eq. 1, the centre-of-mass velocity is calculated as

${\mathcal{V}}_{\text{cm}}=\frac{1}{N}\sum _i{\mathit{v}}_i=\frac{1}{N}\sum _i\frac{1}{{\zeta }_i\left(t\right)}{\mathit{f}}_i^{\mathrm{i}\mathrm{n}\mathrm{t}},(14)$

where ${\mathit{f}}_i^{\mathrm{i}\mathrm{n}\mathrm{t}}$ is defined by Eq. 3. If the substrate adhesion is homogeneous, ζ_i(t) = ζ₀(t) ≠ 0, Eq. 14 becomes

${\mathcal{V}}_{\text{cm}}=\frac{1}{N{\zeta }_0\left(t\right)}\sum _i{\mathit{f}}_i^{\text{int}}=0.(15)$

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 force-free condition.

3.2 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):

$c_i\left(t\right)=c_0\left(1-\cos {\varphi }_i\left(t\right)\right),(16)$

where the phase changes as

${\varphi }_i\left(t\right)=\omega t-q_w\left(x_i-x^{*}\right).(17)$

Here c₀ is the maximum concentration, which is set to c₀ = 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

${\ell }_{ij}^{\text{act}}\left(t\right)={\ell }_c\frac{c_i\left(t\right)+c_j\left(t\right)}{2},(18)$

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:

${\zeta }_i\left(t\right)=\left\{\begin{array}{cc}{\zeta }_{\text{slip}}\hfill & \mathrm{i}\mathrm{f}2m_i\pi <{\varphi }_i\left(t\right)-{\psi }_{\zeta }\le \left(2m_i+1\right)\pi \hfill \\ {\zeta }_{\text{stick}}\hfill & \mathrm{i}\mathrm{f}\left(2m_i+1\right)\pi <{\varphi }_i\left(t\right)-{\psi }_{\zeta }\le 2\left(m_i+1\right)\pi \hfill \end{array}\right.(19)$

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 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.

FIGURE 5

FIGURE 5. Cell crawling driven by a prescribed simple wave of intracellular chemical concentration, that is given in the form of a harmonic wave with the frequency ω = 2π and the wavelength L_w that travels from the rightmost subcellular elements to the left. (A) Migration distance for different values of the wavelength L_w and the phase shift ψ_ζ. The displacement in one cycle is measured after several cycles of relaxation. The sign of ΔR distinguishes the direction of motion. The positive and negative signs of ΔR correspond to the rightward and leftward motion, respectively. The frequency of the traveling wave is set to ω = 2π. [(B)–(D)] Time series of crawling cells for q_w = π. (B) Rightward motion for ψ_ζ = 5π/9, (C) leftward motion for ψ_ζ = 14π/9, and (D) reciprocating motion for ψ_ζ = 4π/90. In panels (B)–(D), the numbers show the time, and the color indicates the distribution of the intracellular chemical component c(t). The size of the subcellular elements represents the substrate friction; large and small elements correspond to the adhered stick state and the de-adhered slip state, respectively.

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.

3.3 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 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.

FIGURE 6

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.

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.

3.4 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.

FIGURE 7

FIGURE 7. Normalized migration distance ΔR/L_cell (A–C) against the time delay τ_ζ and (D–F) against A_v for (A,D) a hexagonal cell (L_cell = 1), (B,E) a circular cell (L_cell = 1), and (C,F) a large hexagonal cell (L_cell = 1.4). In the panels (A–C), the results for several values of A_v are plotted by different lines, whereas the dependence on A_v is shown for τ_ζ = 0.01 and 0.1 in the panels (D–F).

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).

3.5 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.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.

FIGURE 8

FIGURE 8. Cell with a circular shape. (A) Circular cell at rest. Subcellular elements plotted as magenta circles are connected by viscoelastic bonds indicated by blue lines. The details of the bond are the same as in Figure 2B. (B) Time series of crawling circular cells for A_v = 0.6, v* = 1, V* = 0.5, and τ_ζ = 0.01. In each subplot in panel (B), the number in the bottom left corner shows the time, and the color indicates the value of V(t). The size of each subcellular element represents the value of the substrate friction coefficient.

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.

3.7 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.

FIGURE 9

FIGURE 9. Cell crawling with randomly-chosen activator element for (A,B) K_K = 5 and (C,D) K_K = 5.5. (E) Mean-square displacement (MSD) and (F) velocity-autocorrelation function (VAC) for K_K = 5 and 5.5. In panels (A) and (C), the trajectory of the center-of-mass position is plotted, where the gray circles indicate the initial position at t = 0, and the red triangles represent the position at every time interval of 10. The orange squares are the final position at time 40. (B,D) Time series of snapshots, (B) where the crawling direction of the cell changes and (D) where the cell undergoes a spinning motion. The numbers at the left bottom corners of the subplots in panels (B) and (D) show the time. The other parameters are the same as in Table 3.

To quantify the persistence of the migration direction, we measure the mean-square displacement (MSD), which is defined by $⟨{({\mathcal{X}}_{\text{cm}}(t_0+t)-{\mathcal{X}}_{\text{cm}}(t_0))}^2⟩$, where ${\mathcal{X}}_{\text{cm}}$ is the center-of-mass position of the cell and ⟨x (t₀)⟩ represents the time average of x (t₀), i.e., the average of x (t₀) over time t₀. The results are displayed in Figure 9E. It shows a crossover from ballistic (∝ t²) to diffusive regimes (∝ t¹) at around t ∼ 1. In fact, the velocity autocorrelation function (VAC) defined by $⟨{\tilde{\mathcal{V}}}_{\text{cm}}(t_0+t)\cdot {\tilde{\mathcal{V}}}_{\text{cm}}(t_0)⟩$ for the direction of the center-of-mass velocity ${\tilde{\mathcal{V}}}_{\text{cm}}{=\mathcal{V}}_{\text{cm}}/{|\mathcal{V}}_{\text{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 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¹ 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.

3.8 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: ${\mathit{f}}_i^{\text{traction}}={\zeta }_i(t){\mathit{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

$M_{\alpha }^{\left(1\right)}=\sum _i{\mathit{f}}_{i,\alpha }^{\text{traction}},(20)$

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.

FIGURE 10

FIGURE 10. Traction force multipole of the randomly crawling cells in Figure 9. (A,B) Times series of each component of traction force monopole ($M_1^{(1)}$ and $M_2^{(1)}$) and the diagonal ($M_{11}^{(2)}$ and $M_{22}^{(2)}$) and off-diagonal components ($M_{12}^{(2)}$ and $M_{21}^{(2)}$) of the traction force dipole, and the two eigenvalues of the force dipole tensor (${\lambda }_1^{(2)}$ and ${\lambda }_2^{(2)}$) and the torque (T), as well as the angle of the velocity (φ_v) and the two eigenvectors of the traction force dipole (θ₁ and θ₂). (C,D) Time evolution of traction force dipole $M_{11}^{(2)}$ and quadrupole $M_{111}^{(3)}$. The color indicates the time. The data in panels (A,C) are for the cell in Figure 9A, whereas those in panels (B,D) are for the cell in Figure 9B. Note that the axis of the $M_{11}^{(2)}$ is inverted, to match the plot in Ref. Tanimoto and Sano (2014) in panels (C,D).

The next lowest mode is the traction force dipole, the element of which is defined by

$M_{\alpha \beta }^{\left(2\right)}=\sum _i{r}_{i,\alpha }{f}_{i,\beta }^{\text{traction}}.(21)$

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_{12}^{(2)}$ and $M_{21}^{(2)}$ 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

$\mathit{T}=\sum _i{\mathit{r}}_i\times {\mathit{f}}_i^{\text{traction}},(22)$

which, in a two dimensional space, becomes

$T=\sum _i\left(r_{i,1}{f}_{i,2}^{\text{traction}}-r_{i,2}{f}_{i,1}^{\text{traction}}\right)=M_{12}^{\left(2\right)}-M_{21}^{\left(2\right)}(23)$

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 λ₁ and λ₂ of the traction force dipole look identical to those of $M_{11}^{(2)}$ and $M_{22}^{(2)}$, in particular in Figure 10A. This is because the direction θ₁ of the eigenvector of M⁽²⁾ 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 λ₁ is selected such that the direction of the corresponding eigenvector (θ₁) is closer to φ_v than that of the other eigenvalue λ₂. 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 θ₁ of the eigenvector corresponding to λ₁ 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 λ₁ and λ₂ from those of $M_{11}^{(2)}$ and $M_{22}^{(2)}$.

Finally, in Figure 10C, we plot the time evolution of the traction force dipole $M_{11}^{(2)}$ and the traction force quadrupole $M_{111}^{(3)}$, which is compared with the measurement in the experiment Tanimoto and Sano (2014). Here the traction force quadrupole is defined by

$M_{\alpha \beta \gamma }^{\left(3\right)}=\sum _i{r}_{i,\alpha }{r}_{i,\beta }{f}_{i,\gamma }^{\text{traction}}.(24)$

Interestingly, the trajectory in $M_{11}^{(2)}$-$M_{111}^{(3)}$ 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_{11}^{(2)}$ is inverted, to make the comparison easy.

4 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 subcellular-element 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, 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 de-adhered 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, ${\ell }_{ij}^{\text{act}}(t)$, of the bond connecting two subcellular elements. The two processes are distinguished by the sign of ${\ell }_{ij}^{\text{act}}(t)$. In this paper, however, we consider only the protrusion process, i.e., ${\ell }_{ij}^{\text{act}}(t)\ge 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 two-dimensional 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 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 force- and 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.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ananthakrishnan, R., and Ehrlicher, A. (2007). The forces behind cell movement. Int. J. Biol. Sci. 3, 303–317. doi:10.7150/ijbs.3.303

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnhart, E., Lee, K.-C., Allen, G. M., Theriot, J. A., and Mogilner, A. (2015). Balance between cell?substrate adhesion and myosin contraction determines the frequency of motility initiation in fish keratocytes. Proc. Natl. Acad. Sci. U. S. A. 112, 5045–5050. doi:10.1073/pnas.1417257112

PubMed Abstract | CrossRef Full Text | Google Scholar

Basan, M., Prost, J., Joanny, J.-F., and Elgeti, J. (2011). Dissipative particle dynamics simulations for biological tissues: Rheology and competition. Phys. Biol. 8, 026014. doi:10.1088/1478-3975/8/2/026014

PubMed Abstract | CrossRef Full Text | Google Scholar

Bausch, A. R., Möller, W., and Sackmann, E. (1999). Measurement of local viscoelasticity and forces in living cells by magnetic tweezers. Biophys. J. 76, 573–579. doi:10.1016/S0006-3495(99)77225-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bechinger, C., Di Leonardo, R., Löwen, H., Reichhardt, C., Volpe, G., and Volpe, G. (2016). Active particles in complex and crowded environments. Rev. Mod. Phys. 88, 045006. doi:10.1103/RevModPhys.88.045006

CrossRef Full Text | Google Scholar

Beta, C., and Kruse, K. (2017). Intracellular oscillations and waves. Annu. Rev. Condens. Matter Phys. 8, 239–264. doi:10.1146/annurev-conmatphys-031016-025210

CrossRef Full Text | Google Scholar

Blanchoin, L., Boujemaa-Paterski, R., Sykes, C., and Plastino, J. (2014). Actin dynamics, architecture, and mechanics in cell motility. Physiol. Rev. 94, 235–263. doi:10.1152/physrev.00018.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosgraaf, L., and Van Haastert, P. J. M. (2009). The ordered extension of pseudopodia by amoeboid cells in the absence of external cues. PLoS ONE 4, e5253. doi:10.1371/journal.pone.0005253

PubMed Abstract | CrossRef Full Text | Google Scholar

Cates, M. E. (2012). Diffusive transport without detailed balance in motile bacteria: Does microbiology need statistical physics? Rep. Prog. Phys. 75, 042601. doi:10.1088/0034-4885/75/4/042601

PubMed Abstract | CrossRef Full Text | Google Scholar

Devreotes, P. N., Bhattacharya, S., Edwards, M., Iglesias, P. A., Lampert, T., and Miao, Y. (2017). Excitable signal transduction networks in directed cell migration. Annu. Rev. Cell. Dev. Biol. 33, 103–125. PMID: 28793794. doi:10.1146/annurev-cellbio-100616-060739

PubMed Abstract | CrossRef Full Text | Google Scholar

Epstein, I. R., and Pojman, J. A. (1998). “An introduction to nonlinear chemical dynamics : Oscillations, waves, patterns, and chaos,” in Topics in physical chemistry (Oxford, United Kingdom: Oxford University Press).

Google Scholar

Flemming, S., Font, F., Alonso, S., and Beta, C. (2020). How cortical waves drive fission of motile cells. Proc. Natl. Acad. Sci. U. S. A. 117, 6330–6338. doi:10.1073/pnas.1912428117

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukushima, S., Matsuoka, S., and Ueda, M. (2019). Excitable dynamics of ras triggers spontaneous symmetry breaking of pip3 signaling in motile cells. J. Cell. Sci. 132, jcs224121. doi:10.1242/jcs.224121

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerhardt, M., Ecke, M., Walz, M., Stengl, A., Beta, C., and Gerisch, G. (2014). Actin and pip3 waves in giant cells reveal the inherent length scale of an excited state. J. Cell. Sci. 127, 4507–4517. doi:10.1242/jcs.156000

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerisch, G., Ecke, M., Wischnewski, D., and Schroth-Diez, B. (2011). Different modes of state transitions determine pattern in the phosphatidylinositide-actin system. BMC Cell. Biol. 12, 42–16. doi:10.1186/1471-2121-12-42

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerisch, G., Schroth-Diez, B., Müller-Taubenberger, A., and Ecke, M. (2012). Pip3 waves and pten dynamics in the emergence of cell polarity. Biophys. J. 103, 1170–1178. doi:10.1016/j.bpj.2012.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Golebiewska, U., Nyako, M., Woturski, W., Zaitseva, I., McLaughlin, S., and Martin, T. U. (2008). Diffusion coefficient of fluorescent phosphatidylinositol 4, 5-bisphosphate in the plasma membrane of cells. Mol. Biol. Cell. 19, 1663–1669. doi:10.1091/mbc.e07-12-1208

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, P., and Scott, S. (1983). Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Isolas and other forms of multistability. Chem. Eng. Sci. 38, 29–43. doi:10.1016/0009-2509(83)80132-8

CrossRef Full Text | Google Scholar

Gray, P., and Scott, S. (1984). Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Oscillations and instabilities in the system a + 2b → 3b; b → c. Chem. Eng. Sci. 39, 1087–1097. doi:10.1016/0009-2509(84)87017-7

CrossRef Full Text | Google Scholar

Günther, S., and Kruse, K. (2008). A simple self-organized swimmer driven by molecular motors. Europhys. Lett. 84, 68002. doi:10.1209/0295-5075/84/68002

CrossRef Full Text | Google Scholar

Huang, C.-H., Tang, M., Shi, C., Iglesias, P. A., and Devreotes, P. N. (2013). An excitable signal integrator couples to an idling cytoskeletal oscillator to drive cell migration. Nat. Cell. Biol. 15, 1307–1316. doi:10.1038/ncb2859

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwamoto, M., Ueyama, D., and Kobayashi, R. (2014). The advantage of mucus for adhesive locomotion in gastropods. J. Theor. Biol. 353, 133–141. doi:10.1016/j.jtbi.2014.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Kasza, K. E., Rowat, A. C., Liu, J., Angelini, T. E., Brangwynne, C. P., Koenderink, G. H., et al. (2007). The cell as a material. Curr. Opin. Cell. Biol.Cell Struct. Dyn. 19, 101–107. doi:10.1016/j.ceb.2006.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Koshizuka, S., Nobe, A., and Oka, Y. (1998). Numerical analysis of breaking waves using the moving particle semi-implicit method. Int. J. Numer. Methods Fluids 26, 751–769. doi:10.1002/(SICI)1097-0363(19980415)26:7¡751:AID-FLD671¿3.0.CO;2-C

CrossRef Full Text | Google Scholar

Kumar, K. V., Ramaswamy, S., and Rao, M. (2008). Active elastic dimers: Self-propulsion and current reversal on a featureless track. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 77, 020102. doi:10.1103/PhysRevE.77.020102

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuramoto, Y. (1984). Chemical oscillations, waves, and turbulence. Berlin, Germany: Springer-Verlag Berlin Heidelberg.

Google Scholar

Lauga, E., and Powers, T. R. (2009). The hydrodynamics of swimming microorganisms. Rep. Prog. Phys. 72, 096601. doi:10.1088/0034-4885/72/9/096601

CrossRef Full Text | Google Scholar

Li, L., Nørrelykke, S. F., and Cox, E. C. (2008). Persistent cell motion in the absence of external signals: A search strategy for eukaryotic cells. PLoS ONE 3, e2093. doi:10.1371/journal.pone.0002093

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopez, J. H., Das, M., and Schwarz, J. M. (2014). Active elastic dimers: Cells moving on rigid tracks. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 90, 032707. doi:10.1103/PhysRevE.90.032707

PubMed Abstract | CrossRef Full Text | Google Scholar

Maeda, Y. T., Inose, J., Matsuo, M. Y., Iwaya, S., and Sano, M. (2008). Ordered patterns of cell shape and orientational correlation during spontaneous cell migration. PLoS ONE 3, e3734. doi:10.1371/journal.pone.0003734

PubMed Abstract | CrossRef Full Text | Google Scholar

Mai, M. H., and Camley, B. A. (2020). Hydrodynamic effects on the motility of crawling eukaryotic cells. Soft Matter 16, 1349–1358. doi:10.1039/C9SM01797F

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchetti, M. C., Joanny, J. F., Ramaswamy, S., Liverpool, T. B., Prost, J., Rao, M., et al. (2013). Hydrodynamics of soft active matter. Rev. Mod. Phys. 85, 1143–1189. doi:10.1103/RevModPhys.85.1143

CrossRef Full Text | Google Scholar

Micoulet, A., Spatz, J. P., and Ott, A. (2005). Mechanical response analysis and power generation by single-cell stretching. ChemPhysChem 6, 663–670. doi:10.1002/cphc.200400417

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, J. D. (2002). “Mathematical biology I. An introduction,” in Interdisciplinary applied mathematics. 3 edn (New York: Springer-Verlag).

Google Scholar

Murray, J. D. (2003). “Mathematical Biology II. Spatial models and biomedical applications,” in Interdisciplinary applied mathematics. 3 edn (New York: Springer-Verlag).

Google Scholar

Najafi, A., and Golestanian, R. (2004). Simple swimmer at low Reynolds number: Three linked spheres. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 69, 062901. doi:10.1103/PhysRevE.69.062901

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, T. J. (2007). Modeling multicellular structures using the subcellular element model. Basel: Birkhäuser Basel, 221–239. doi:10.1007/978-3-7643-8123-3_10

CrossRef Full Text | Google Scholar

Niculescu, I., Textor, J., and de Boer, R. J. (2015). Crawling and gliding: A computational model for shape-driven cell migration. PLoS Comput. Biol. 11, 10042800–e1004322. doi:10.1371/journal.pcbi.1004280

PubMed Abstract | CrossRef Full Text | Google Scholar

Nier, V., Jain, S., Lim, C. T., Ishihara, S., Ladoux, B., and Marcq, P. (2016). Inference of internal stress in a cell monolayer. Biophys. J. 110, 1625–1635. doi:10.1016/j.bpj.2016.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishikawa, M., Hörning, M., Ueda, M., and Shibata, T. (2014). Excitable signal transduction induces both spontaneous and directional cell asymmetries in the phosphatidylinositol lipid signaling system for eukaryotic chemotaxis. Biophys. J. 106, 723–734. doi:10.1016/j.bpj.2013.12.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishimura, S. I., Ueda, M., and Sasai, M. (2009). Cortical factor feedback model for cellular locomotion and cytofission. PLoS Comput. Biol. 5, e1000310. doi:10.1371/journal.pcbi.1000310

PubMed Abstract | CrossRef Full Text | Google Scholar

Pismen, L. M. (2006). Patterns and interfaces in dissipative dynamics. Berlin, Germany: Springer-Verlag Berlin Heidelberg.

Google Scholar

Purcell, E. M. (1977). Life at low Reynolds number. Am. J. Phys. 45, 3–11. doi:10.1119/1.10903

CrossRef Full Text | Google Scholar

Ramaswamy, S. (2010). The mechanics and statistics of active matter. Annu. Rev. Condens. Matter Phys. 1, 323–345. doi:10.1146/annurev-conmatphys-070909-104101

CrossRef Full Text | Google Scholar

Sandersius, S. A., and Newman, T. J. (2008). Modeling cell rheology with the subcellular element model. Phys. Biol. 5, 015002. doi:10.1088/1478-3975/5/1/015002

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwarz, U. S., and Safran, S. A. (2013). Physics of adherent cells. Rev. Mod. Phys. 85, 1327–1381. doi:10.1103/RevModPhys.85.1327

CrossRef Full Text | Google Scholar

Shao, D., Rappel, W.-J., and Levine, H. (2010). Computational model for cell morphodynamics. Phys. Rev. Lett. 105, 108104. doi:10.1103/PhysRevLett.105.108104

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, C., Huang, C.-H., Devreotes, P. N., and Iglesias, P. A. (2013). Interaction of motility, directional sensing, and polarity modules recreates the behaviors of chemotaxing cells. PLoS Comput. Biol. 9, 10031222–e1003217. doi:10.1371/journal.pcbi.1003122

PubMed Abstract | CrossRef Full Text | Google Scholar

Smeets, B., Alert, R., Pešek, J., Pagonabarraga, I., Ramon, H., and Vincent, R. (2016). Emergent structures and dynamics of cell colonies by contact inhibition of locomotion. Proc. Natl. Acad. Sci. U. S. A. 113, 14621–14626. doi:10.1073/pnas.1521151113

PubMed Abstract | CrossRef Full Text | Google Scholar

Style, R. W., Boltyanskiy, R., German, G. K., Hyland, C., MacMinn, C. W., Mertz, A. F., et al. (2014). Traction force microscopy in physics and biology. Soft Matter 10, 4047–4055. doi:10.1039/C4SM00264D

PubMed Abstract | CrossRef Full Text | Google Scholar

Swaney, K. F., Huang, C.-H., and Devreotes, P. N. (2010). Eukaryotic chemotaxis: A network of signaling pathways controls motility, directional sensing, and polarity. Annu. Rev. Biophys. 39, 265–289. doi:10.1146/annurev.biophys.093008.131228

PubMed Abstract | CrossRef Full Text | Google Scholar

Taniguchi, D., Ishihara, S., Oonuki, T., Honda-Kitahara, M., Kaneko, K., and Sawai, S. (2013). Phase geometries of two-dimensional excitable waves govern self-organized morphodynamics of amoeboid cells. Proc. Natl. Acad. Sci. U. S. A. 110, 5016–5021. doi:10.1073/pnas.1218025110

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanimoto, H., and Sano, M. (2014). A simple force-motion relation for migrating cells revealed by multipole analysis of traction stress. Biophys. J. 106, 16–25. doi:10.1016/j.bpj.2013.10.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarama, M., and Shibata, T. (2022). Pattern formation and the mechanics of a motor-driven filamentous system confined by rigid membranes. Phys. Rev. Res. 4, 043071. doi:10.1103/PhysRevResearch.4.043071

CrossRef Full Text | Google Scholar

Tarama, M., and Yamamoto, R. (2018). Mechanics of cell crawling by means of force-free cyclic motion. J. Phys. Soc. Jpn. 87, 044803. doi:10.7566/JPSJ.87.044803

CrossRef Full Text | Google Scholar

Tjhung, E., Tiribocchi, A., Marenduzzo, D., and Cates, M. (2015). A minimal physical model captures the shapes of crawling cells. Nat. Commun. 6, 5420. doi:10.1038/ncomms6420

PubMed Abstract | CrossRef Full Text | Google Scholar

Vicsek, T., and Zafeiris, A. (2012). Collective motion. Phys. Rep.Collective motion 517, 71–140. doi:10.1016/j.physrep.2012.03.004

CrossRef Full Text | Google Scholar

Xiong, Y., Huang, C.-H., Iglesias, P. A., and Devreotes, P. N. (2010). Cells navigate with a local-excitation, global-inhibition-biased excitable network. Proc. Natl. Acad. Sci. U. S. A. 107, 17079–17086. doi:10.1073/pnas.1011271107

PubMed Abstract | CrossRef Full Text | Google Scholar

Ziebert, F., and Aranson, I. S. (2016). Computational approaches to substrate-based cell motility. npj Comput. Mat. 2, 16019. doi:10.1038/npjcompumats.2016.19

CrossRef Full Text | Google Scholar

Ziebert, F., Swaminathan, S., and Aranson, I. S. (2012). Model for self-polarization and motility of keratocyte fragments. J. R. Soc. Interface 9, 1084–1092. doi:10.1098/rsif.2011.0433

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimmermann, J., Camley, B. A., Rappel, W.-J., and Levine, H. (2016). Contact inhibition of locomotion determines cell–cell and cell–substrate forces in tissues. Proc. Natl. Acad. Sci. U. S. A. 113, 2660–2665. doi:10.1073/pnas.1522330113

PubMed Abstract | CrossRef Full Text | Google Scholar