Mechanochemical subcellular-element model of crawling cells

TL;DR

This paper introduces a mechanochemical model of cell crawling that links intracellular chemical signals to movement direction and distance, providing insights into force-free migration mechanisms.

Contribution

It presents the first fully force-free cell migration model driven by intracellular chemical reactions, integrating reaction-diffusion dynamics with mechanical processes.

Findings

Chemical dependence of adhesion controls crawling direction.

Traveling chemical waves influence cell movement.

Traction force analysis aligns with experimental data.

Abstract

Constructing physical models of living cells and tissues is an extremely challenging task because of the high complexities of both intra- and intercellular processes. In addition, the force that a single cell generates vanishes in total due to the law of action and reaction. The typical mechanics of cell crawling involve periodic changes in the cell shape and in the adhesion characteristics of the cell to the substrate. However, the basic physical mechanisms by which a single cell coordinates these processes cooperatively to achieve autonomous migration are not yet well understood. To obtain a clearer grasp of how the intracellular force is converted to directional motion, we develop a basic mechanochemical model of a crawling cell based on subcellular elements with the focus on the dependence of the protrusion and contraction as well as the adhesion and deadhesion processes on intracellular biochemical signals. By introducing 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. To our knowledge, our present work is the first study that accomplishes fully force-free migration utilizing intracellular chemical reactions. Although the detailed mechanisms of actual cells are far more complicated than our simple model, we believe that this mechanochemical model is a good prototype for more realistic models.

Figures

Tables

Table 1. Table 1: Mechanical parameters in the model and their values for the simulation of the hexagonal cell.

Model parameters	Simulation values
Diameter of cell at rest, $L_{\mathrm{cell}}$	1 (1.4 for a larger cell)
Number of subcellular elements, $N$	91 (169 for a larger cell)
Number of bonds, $N_{\mathrm{bond}}$	240 (462 for a larger cell)
Bond length at rest, ${\ell }_{ij}$	$L_{\mathrm{cell}}/\sqrt{N}$
Total bond length at rest, $L_{\mathrm{bond}}$	${\sum }_{⟨i,j⟩}{\ell }_{ij}=N_{\mathrm{bond}}{L}_{\mathrm{cell}}/\sqrt{N}$
Mean bond length, ${\ell }_0$	$L_{\mathrm{bond}}/N_{\mathrm{bond}}=L_{\mathrm{cell}}/\sqrt{N}$
Substrate friction coefficient
in the slip state, ${\zeta }_{\mathrm{slip}}$	$L_{\mathrm{cell}}^2/N$
in the stick state, ${\zeta }_{\mathrm{stick}}$	$10{\zeta }_{\mathrm{slip}}$
Elastic modulus of each bond, $\kappa$	10
Intracellular dissipation rate, $\xi$	$L_{\mathrm{cell}}^2/L_{\mathrm{bond}}$
Dimensionless parameter, $\omega {\xi }_{\mathrm{cell}}/k_{\mathrm{cell}}$	$2\pi {\xi }_{\mathrm{cell}}/k_{\mathrm{cell}}{T}_0=0.2\pi L_{\mathrm{cell}}^3$

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

Mechanical parameters	Model	Living cells	References
Migration speed		$\approx 10\mu m{\left(\min \right)}^{-1}$	Maeda et al. (2008)
Traction stress		$\approx 100\mathrm{Pa}$	Tanimoto and Sano (2014)
Traction force		$10\mathrm{nN}$	(estimated)
Cell-substrate friction coefficient	${\zeta }_{\mathrm{cell}}={\zeta }_{\mathrm{slip}}N$	$1\mathrm{nN}{\left(\mu m\right)}^{-1}\min$	(estimated)
Elastic modulus	${\kappa }_{\mathrm{cell}}=\kappa /L_{\mathrm{cell}}$	10–$100\mathrm{Pa}$	Bausch et al. (1999); Micoulet et al. (2005)
Elastic constant	$k_{\mathrm{cell}}=\kappa /L_{\mathrm{cell}}$	$0.1$–$1\mathrm{nN}{\left(\mu m\right)}^{-1}$	(estimated)
Dissipation rate	${\xi }_{\mathrm{cell}}=\xi L_{\mathrm{bond}}$
Fluid/solid transition time	${\xi }_{\mathrm{cell}}/k_{\mathrm{cell}}$	$\approx 1\sec$	Kasza et al. (2007)

Table 3. Table 3: The parameter values of the RD equations in the simulation.

Model Parameters	Values used in simulations
$D_U$	0.48
$D_V$	0.48
$\alpha$	240
$\beta$	90
$K_K$	5
$K_P$	5
$S$	30
$\gamma$	6
$\mu$	30

Equations

$${\zeta }_i(t){\mathbf{v}}_i+\sum _{j\in {\Omega }_i}\xi {\ell }_{ij}({\mathbf{v}}_i-{\mathbf{v}}_j)$$

$$=\sum _{j\in {\Omega }_i}\frac{\kappa }{{\ell }_{ij}}{\hat{\mathbf{r}}}_{ij}\left\{r_{ij}-\left({\ell }_{ij}+{\ell }_{ij}^{\mathrm{act}}(t)\right)\right\}+{\mathbf{f}}_i^{\mathrm{area}},$$

$$\sum _i{\mathbf{f}}_i^{\mathrm{int}}=0.$$

$${\mathbf{f}}_i^{\mathrm{int}}=$$

$$+\sum _{j\in {\Omega }_i}\frac{\kappa }{{\ell }_{ij}}{\hat{\mathbf{r}}}_{ij}\left\{r_{ij}-\left({\ell }_{ij}+{\ell }_{ij}^{\mathrm{act}}(t)\right)\right\}+{\mathbf{f}}_i^{\mathrm{area}}$$

$$\frac{\partial U_i}{\partial t}=D_U{\nabla }^2{U}_i+G_U(U_i,V_i),$$

$$\frac{\partial V_i}{\partial t}=D_V{\nabla }^2{V}_i+G_V(U_i,V_i),$$

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

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

$$h_{\zeta}(\zeta)=-\frac{1}{2}\tanh\Big{(}\frac{\zeta-\zeta_{\rm stick}}{\epsilon_{\zeta}}\Big{)}-\frac{1}{2}\tanh\Big{(}\frac{\zeta-\zeta_{\rm slip}}{\epsilon_{\zeta}}\Big{)},$$

$$h_v(v)=\frac{(v/v^{\ast }{)}^2}{1+(v/v^{\ast }{)}^2}-\frac{1}{2},$$

$$h_V(V)=\frac{1}{2}\tanh ({\sigma }_V(V-V^{\ast })).$$

$$M_{\alpha }^{(1)}=\sum _i{\mathbf{f}}_{i,\alpha }^{\mathrm{traction}},$$

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

$$\mathbf{T}=\sum _i{\mathbf{r}}_i\times {\mathbf{f}}_i^{\mathrm{traction}},$$

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

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

$$\frac{\partial U}{\partial t}=D_U{\nabla }^{2}U+G_U(U,V),$$

$$\frac{\partial V}{\partial t}=D_V{\nabla }^{2}V+G_V(U,V),$$

$$G_U(U,V)=-\frac{\alpha U{V}^2}{K_K+⟨V^2⟩}+\frac{\beta UV}{K_P+⟨U⟩}+S-\gamma U,$$

$$G_V(U,V)=+\frac{\alpha U{V}^2}{K_K+⟨V^2⟩}-\frac{\beta UV}{K_P+⟨U⟩}-\mu V.$$

$$⟨U⟩=\frac{1}{N}\sum _i{U}_i,⟨V^2⟩=\frac{1}{N}\sum _i{V}_i^2,$$

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

$$w(r)=\left\{\begin{array}[]{ll}\frac{r_{e}}{r}-1&\textrm{for~{}}0\leq r<r_{e}\\ 0&\textrm{for~{}}r_{e}\leq r\end{array}\right.$$

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

$${\varphi }_i(t)=\omega t-q_w(x_i-x^{\ast }),$$

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

$$\zeta_{i}(t)=\left\{\begin{array}[]{ll}\zeta_{\rm slip}&\textrm{if~{}}2m_{i}\pi<\phi_{i}(t)-\psi_{\zeta}\leq(2m_{i}+1)\pi\\ \zeta_{\rm stick}&\textrm{if~{}}(2m_{i}+1)\pi<\phi_{i}(t)-\psi_{\zeta}\leq 2(m_{i}+1)\pi\end{array}\right.$$

Full text

I Introduction

Cells are the basic units of living organisms. Despite the diversity of cellular types, they share similar structures. Biological cells are composed of a number of polymers, proteins, and lipids, which form stable structures such as the lipid membrane. At the same time, cells exist in a state far from equilibrium. Inside a cell, complex chemical reactions take place and are converted to mechanical forces by molecular motors. Consequently, biological cells spontaneously exhibit various dynamics, such as locomotion and proliferation.

In general, objects that exhibit spontaneous motion are called active matter Lauga and Powers (2009); Ramaswamy (2010); Vicsek and Zafeiris (2012); Cates (2012); Marchetti et al. (2013); Bechinger et al. (2016). In contrast to objects passively driven by external forcing, active matter, including living cells, generates force in itself, which is characterized by a vanishing force monopole due to the action-reaction law. Under this force-free condition, it is necessary to break symmetry to achieve spontaneous motion, such as directional motion. For microorganisms that swim in a fluidic environment, the scallop theorem Purcell (1977) describes the importance of breaking reciprocality to achieve net migration via internal cyclic motions.

In nature, there exist a number of microorganisms that crawl on substrates, such as the extracellular matrix and other cells. In contrast to the locomotion of microswimmers, adhesion to the substrate plays an important role in the locomotion of crawling microorganisms. Such crawling motion is often observed in eukaryotic cells, such as Dictyostelium cells and keratocytes.

Many models have been introduced to explain various aspects of the dynamics of crawling cells. Examples include the vertex model Honda (1978); Honda et al. (2004, 2008); Bi et al. (2016), the cellular Potts model Nishimura et al. (2009); Niculescu et al. (2015), the continuous model Marchetti et al. (2013); Nier et al. (2016), the phase field model Ziebert and Aranson (2016); Shao et al. (2010); Taniguchi et al. (2013); Shi et al. (2013); Ziebert et al. (2012); Tjhung et al. (2015), and the particle-based model Newman (2007); Sandersius and Newman (2008); Basan et al. (2011); Zimmermann et al. (2016); Smeets et al. (2016). Here, we construct a model based on the particle-based model.

In our previous study Tarama and Yamamoto (2018), we systematically investigated the cycle of the typical crawling mechanism Ananthakrishnan and Ehrlicher (2007):

1. protrusion at the leading edge,
2. adhesion to the substrate at the leading edge,
3. deadhesion from the substrate at the trailing edge, and
4. contraction at the trailing edge. To clarify the role of this cycle in efficient crawling motion, we introduced a simple mechanical model in which a cell is described by two subcellular elements connected by a viscoelastic bond, which includes an actuator that elongates and shrinks cyclically. The substrate friction of each element switches cyclically between the adhered stick state and the deadhered slip state. By tuning the phase shifts between the actuator elongation and the substrate friction of each element, we demonstrated that the order of the four basic processes of the typical crawling mechanism has a great impact on the crawling distance and efficiency, as well as the crawling direction.

If we consider the extension of the mechanical model to a model cell consisting of far more than two elements or to two- or three-dimensional (2D or 3D) space, adjusting the phase shift of each element “by hand” becomes less feasible. Instead, we need to consider the underlying processes that regulate the actuator elongation and the stick-slip transition of the substrate friction.

The aim of this paper is to develop a basic particle-based model for cell crawling. To this end, we consider a cell crawling on a flat substrate and extend our previous mechanical model to two dimensions. We describe a cell by a set of many subcellular elements connected by viscoelastic bonds Newman (2007). In addition, intracellular chemical reactions are represented by simple 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.

This paper is organized as follows. In the next section, we introduce the model equations that couples cell mechanics and intracellular chemistry. Then, we show that the dependence of the substrate adhesion on the intracellular chemistry determines the direction of the cell crawling in Sect. III. In Sect. IV, we investigate how the cell crawling changes depending on the mechanical and chemical signals that control the substrate adhesion. The impact of the cell shape and size are studied in Sets. V and VI, respectively. In Sect. VII, random crawling motion is realized by random excitation of intracellular chemistry, for which we analyze the traction force multipoles in Sect. VIII. Sect. IX is devoted to the summary and discussion, and this paper concludes with Sect. X.

II Model

First, we introduce our mechanical model of a crawling cell and the RD equations representing intracellular chemical reactions. The choice of RD equations is arbitrary, and we employ a previously introduced model. We then couple the mechanical model and the RD equations, which regulate the intracellular mechanical activities. In particular, we confine ourselves to studying possible couplings between the intracellular chemical and mechanical models.

II.1 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 Fig. 1. Since the typical size of a cell is on the order of ten micrometers, the effect of inertia is negligible. Then, the force balance equation of element $i$ is given by

[TABLE]

where $\bm{v}_{i}$ is the velocity of the element $i$ located at the position $\bm{r}_{i}$. Here, the abbreviations $r=|\bm{r}|$ and $\hat{\bm{r}}=\bm{r}/r$ are used for the relative position $\bm{r}_{ij}=\bm{r}_{j}-\bm{r}_{i}$. The summation is over the set of the elements $\Omega_{i}$ connected to the element $i$ by the viscoelastic bonds.

The first term on the left-hand side of Eq. (1) represents the substrate friction with coefficient $\zeta_{i}(t)$, which changes over time due to intracellular activity. The second term represents intracellular dissipation with the rate $\xi$. The first term on the right-hand side represents intracellular elasticity with the elastic modulus $\kappa$ and free length $\ell_{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^{\rm act}_{ij}(t)$. Here, $\ell^{\rm act}_{ij}(t)$ represents the actuator elongation, from which the force generated by the actuator is calculated as $\bm{f}^{\rm act}_{ij}(t)=-\kappa\ell^{\rm act}_{ij}(t)\hat{\bm{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

[TABLE]

where

[TABLE]

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, which is given by $\bm{f}^{\rm area}_{i}=-\delta U^{\rm area}/\delta\bm{r}_{i}$, where $U^{\rm area}=\sum_{\langle i,j,k\rangle}\sigma/S_{ijk}^{2}$. This potential $U^{\rm area}$ penalizes shrinking of the area of each triangle $\langle i,j,k\rangle$ formed by connected subcellular elements $i$, $j$, and $k$, which is defined by $S_{ijk}=(\bm{r}_{ij}\times\bm{r}_{ik})\cdot\hat{\bm{e}}_{z}/2$ with $\hat{\bm{e}}_{z}$ as the unit vector perpendicular to the 2D substrate.

We scale the system by $L_{0}=10\,\mu\rm{m}$ for length and $T_{0}=1\,\rm{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\,{\rm 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 Appendix B.

II.2 Chemical reaction

In the model Eq. (1), the effects of the intracellular activities are included in the actuator elongation $\ell^{\rm act}_{ij}(t)$ and the change in the substrate friction coefficient $\zeta_{i}(t)$. The former represents the protrusion and contraction processes. The latter corresponds to the adhesion and deadhesion of the cell to the underlying substrate. In actual cells, such cellular activities are caused by various intracellular chemical signals. However, it is not realistic to include all chemical components and their signaling pathways. Therefore, we model the intracellular chemical reactions by simple RD equations.

Here, we employ the RD equations proposed by Taniguchi et al. Taniguchi et al. (2013), which are two-component activator-inhibitor equations:

[TABLE]

where $U_{i}$ and $V_{i}$ represent the inhibitor and activator concentrations for subcellular element $i$, respectively. In the original paper Taniguchi et al. (2013), Eq. (5) were introduced to model the phosphoinositide signaling pathway of Dictyostelium cells, and the activator $U$ and the inhibitor $V$ corresponded to phosphatidylinositol (3,4,5)-trisphosphate (PIP3) and phosphatidylinositol (4,5)-bisphosphate (PIP2) concentrations, respectively. The details of these RD equations are given in Appendix C.

The important property of the RD equations, Eq. (5), is that they are of the 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. Taniguchi et al. claimed that the signaling pathway that they were modeling was excitable, and thus, they considered the parameter region of the excitable case to successfully reproduce the experimental results Taniguchi et al. (2013). Interestingly, similar RD equations were studied by Shao et al. Shao et al. (2010) in the context of cell crawling. However, they assumed a bistable regime to reproduce the steady migration of keratocyte cells.

In this paper, we consider the excitable case with the parameters summarized in Table 3, following the study by Taniguchi et al. Taniguchi et al. (2013). The Laplacian terms in Eq. (5) are calculated by using the moving particle semi-implicit (MPS) method Koshizuka et al. (1998). See also Appendix C for further information.

II.3 Mechanochemical coupling

To combine the cell mechanics, Eq. (1), and the RD equations, Eq. (5), we consider the coupling of the chemical concentrations to the actuator elongation $\ell^{\rm act}_{ij}(t)$ and the substrate friction coefficient $\zeta_{i}(t)$ individually.

Before introducing the coupling, we tested a specified traveling wave that couples to the actuator elongation and the substrate friction change. We introduced a time delay between the elongation and substrate friction change. The results showed that there exist an optimum time delay and optimum wavelength for which the cell exhibits the largest migration distance, as summarized in Appendix D.

II.3.1 Actuator elongation

First, we introduce the coupling between the RD equations and the actuator elongation. In Ref. Taniguchi et al. (2013), actin polymerization was found to be enhanced with increasing PIP3 concentration. Therefore, we presume that the actuator elongation depends on the PIP3 concentration as

[TABLE]

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. (5). Therefore, $\tanh$ is introduced on the right hand side of Eq. (6) to prevent an extremely large elongation. $a$ is a constant denoting sensitivity, and $\ell_{V}$ is the magnitude of the elongation. Here, we set $a=\pi$ and $\ell_{V}=\ell_{ij}$.

II.3.2 Substrate adhesion

Next, we consider the adhesion to the substrate underneath and the deadhesion from it. We model the adhesion/deadhesion processes by the transition of the substrate friction coefficient between the adhered stick state and the deadhered slip state. Here, we consider the dependence of the substrate friction coefficient on both mechanical and chemical signals:

[TABLE]

where $\tau_{\zeta}$ 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. See also Fig. 2 for the plot of these functions.

The function $h_{\zeta}(\zeta)$ is defined by

[TABLE]

where $\zeta_{\rm stick}$ and $\zeta_{\rm slip}$ are the substrate friction coefficients in the adhered stick state and the deadhered slip state, respectively. The small parameter $\epsilon_{\zeta}$ indicates the sharpness of the adhesion-deadhesion transition. Here, we set the transition sharpness to $\epsilon_{\zeta}=\zeta_{\rm slip}/2$. We note that, if there are no change in the substrate friction coefficient, the cell does not exhibit any translational motion.

If we consider an artificial vesicle or droplet sitting on a substrate, its adhesion strength changes depending on the force acting on it Schwarz and Safran (2013). The term $h_{v}(v)$ in Eq. (7) represents this dependence of the cell adhesion to the substrate. Here, instead of the force acting on each subcellular element, we presume that the local velocity changes the adhesion strength through

[TABLE]

where $v^{*}$ is the threshold value. The subcellular element tends to adhere to the substrate ($\zeta=\zeta_{\rm stick}$) if the speed is smaller than the threshold value, i.e., $v<v^{*}$, while the element slips on the substrate ($\zeta=\zeta_{\rm 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. (7) with $A_{v}=1$, was introduced in Ref. Barnhart et al. (2015). As a result of the balance of the two functions, $h_{\zeta}(\zeta)$ 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:

[TABLE]

In Eq. (10), $\sigma_{V}$ stands for the sensitivity, and $V^{*}$ is the threshold concentration. Due to this term, a large value of $V$ prevents strong adhesion if $\sigma_{V}>0$, while large $V$ enhances the adhesion if $\sigma_{V}<0$. However, we are not sure whether PIP3 enhances or diminishes adhesion.

III Crawling by direct and retrograde waves

To determine the sign of $\sigma_{V}$, we numerically integrated Eqs. (1)–(10) for both positive and negative $\sigma_{V}$. See Appendix A for the details of the simulation. Fig. 3 depicts a time series of snapshots of a crawling cell for $\sigma_{V}=2\pi$ in Fig. 3(a) and $\sigma_{V}=-2\pi$ in Fig. 3(b), respectively. Here, we set the threshold in Eq. (10) to $V^{*}=0.5$ throughout this paper.

If $\sigma_{V}$ is positive, the cell moves in the opposite direction to the PIP3 traveling wave, as shown in Fig. 3(a). Interestingly, however, if the sign of $\sigma_{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 Fig. 3(b).

With respect to the migration direction, the traveling wave in the same direction is called the direct wave, while the one in the opposite direction is referred to as the retrograde wave Iwamoto et al. (2014). In this sense, the above crawling motion for positive $\sigma_{V}$ in Fig. 3(a) corresponds to the motion with the retrograde wave, and the one for $\sigma_{V}<0$ in Fig. 3(b) corresponds to the motion with the direct wave. Since the experiments in Ref. Taniguchi et al. (2013) show that cells move in the direction in which PIP3 concentration is increased and thus actin polymerization is enhanced, we set $\sigma_{V}=2\pi$ for the rest of this paper.

IV Mechanical vs. chemical control of adhesion

Now, we study the effect of the mechanical and chemical dependence of substrate friction coefficient on cell crawling. To characterize the cell migration, we measure the migration distance $\Delta R$ of the cell in one cycle, i.e., with one traveling wave, varying the time delay $\tau_{\zeta}$. Varying $A_{v}$ between 0 and 1 tunes the ratio of the mechanical and chemical dependences.

The results are summarized in Fig. 4(a), where different lines correspond to different values of $A_{v}$ as indicated in the legend. Interestingly, the mixing of the mechanical and chemical dependences of the substrate friction may result in larger values of $\Delta R$ than purely mechanical or chemical control. In Fig. 4(a), $\Delta R/L_{\rm cell}$ reaches approximately 0.36 for $A_{v}=0.6$ and $\tau_{\zeta}=0.01$. This value, however, is smaller than the crawling of a Dictyostelium cell, which moves its body length in approximately two cycles Tanimoto and Sano (2014).

V Impact of cellular shape

Thus far, we have assumed a hexagonal cell shape, where the structure of the subcellular elements is given by a perfect hexagonal lattice when the cell is at rest. However, this structure does not describe real cells, which are instead circular or often of more complicated shapes. To elucidate the impact of the cell shape on the crawling motion, we prepare a cell of circular shape as described in Appendix E. We then perform the simulation and measure the migration distance for different values of $A_{v}$. The results are plotted in Fig. 4(b) and are qualitatively the same as those of the hexagonal cell in Fig. 4(a).

VI Impact of cell size

Now, we study the impact of the size of the cell. We prepare a hexagonal cell of size $L_{\rm cell}=1.4$, which is approximately twice the area of the previous cells. We again measure the migration distance $\Delta R$, which is normalized by the cell length $L_{\rm cell}$ to facilitate comparison with the previous results for $L_{\rm cell}=1$. The results are summarized in Fig. 4(c). Qualitatively, the tendency is the same as that of the results in Fig. 4(a) for $L_{\rm 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.

VII Random excitation

In reality, cells change their migration direction over time. In our model, we can reproduce such motion by introducing stochasticity, which may originate from, e.g., the complexity of intracellular processes. Here, we randomly choose one element in every $t=0.15$ and add to that element the stimulus $(\delta U,\delta V)=(-I_{\rm excite},+I_{\rm excite})$ with intensity $I_{\rm excite}=0.75$.

The results are summarized in Fig. 5. Due to the stochasticity, the cell changes its migration direction frequently, as shown in the trajectory of the center-of-mass position in Fig. 5(a). From the snapshots of the cell in Fig. 5(b), we can see that the migration direction depends on the position at which the chemical wave occurs. Because of the excitable nature of the RD equations, a stimulus on the element that has relaxed back to the resting state is more likely to be the origin of the next wave. Therefore, a new wave tends to originate from elements that are near the origin of the previous wave. As a result, the cell tends to maintain the same migration direction for a time, as shown in Fig. 5(a). In Figs. 5(a) and (b), the parameter values are kept the same as in Fig. 3(a).

If the parameter $K_{K}$ is slightly increased from 5 to 5.5, the cell changes its migration direction more frequently, as depicted in Fig. 5(c). Depending on the random stimuli, the cell switches from directional motion to spinning motion as a spiral wave appears, as shown in Fig. 5(d). The spinning motion is rather stable, but the cell can also switch back to directional motion in response to a stimulus. Note that the RD equations maintain their excitable nature at this parameter value.

VIII Traction force multipoles

In many experiments of cell crawling, traction force is measured Style et al. (2014) since it is of fundamental importance for cell motility. Here, we measure the traction force of our model cell. Since traction force is the force that a cell exerts on the substrate underneath, it is described by the interaction between the cell and the substrate: $\bm{f}^{\rm traction}_{i}=\zeta_{i}(t)\bm{v}_{i}$. We calculate the traction force force multipoles, which are defined as follows, to compare with the experimental result Tanimoto and Sano (2014).

First, the traction force monopole is defined by

[TABLE]

where $i$ represents the subcellular element $i$ and the summation runs over the entire cell. The subscript $\alpha$ indicates spatial component $\alpha=1,2$. Here, the traction force multipoles are calculated in the comoving coordinate on the cell, and thus, $\alpha=1$ and 2 represent the components parallel and perpendicular to the centre-of-mass velocity, respectively, to compare with the experimental results in Ref. Tanimoto and Sano (2014). In the numerical simulation of the model crawling cell, the traction force monopole is equal to 0, as shown in Fig. 6, which is consistent with the experimental result Tanimoto and Sano (2014). This is readily understood from the force balance equation, Eq. (1), and the force-free condition, Eq. (2)(a). That is, the traction force monopole vanishes because of the fact that the inertia is negligibly small for crawling cells on top of the force-free condition.

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

[TABLE]

On the one hand, from Fig. 6(a), the diagonal components of the traction force dipole are oscillating around 0. Here, note that the positive and negative force dipoles represent the extensile and contractile force. In the experiment Tanimoto and Sano (2014), only the contractile traction force dipole was observed, which our results fail to reproduce. The reason why our model does not reproduce the contractile force dipole is not clear yet. One possibility is that the friction coefficient of the protrusion process may be different from that of the contraction process, which are set to the same in our current model.

On the other hand, the 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. (12). Such symmetric property of the traction force dipole was also obtained in the experiment Tanimoto and Sano (2014).

Now, we consider the meaning of the symmetric property of the traction force dipole that is obtained in our simulation as well as in the experiment. Actually, this symmetric property of the traction force dipole indicates the torque-free nature of the cell. Here, the torque is defined by

[TABLE]

which, in a two dimensional space, becomes

[TABLE]

by using Eq. (12). Note that, if one describes the force dipole tensor with its invariant, the torque appears as the imaginary part of the eigenvalues.

Finally, in Fig. 6(b), we plot the time evolution of the traction force dipole $M^{(2)}_{11}$ and the traction force quadrupole $M^{(3)}_{111}$, which was also measured in the experiment Tanimoto and Sano (2014). Here the traction force quadrupole is defined by

[TABLE]

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

IX Summary and Discussion

To summarize, we have constructed a mechanochemical model of a cell crawling on a substrate. The mechanical part is described by a subcellular-element model, and the chemical part is described by RD equations. To combine them, we introduce two mechanical activities. One is the actuator of the bond connecting each pair of the subcellular elements, which elongates depending on the intracellular chemical concentration. The other is the substrate friction coefficient of each subcellular element, which shows a sharp transition between the adhered stick state and the deadhered slip state. We consider the dependence of the substrate friction coefficient on both the local velocity and the intracellular chemical concentration. We also introduce a time delay of the substrate friction change.

By using this model, we clarified that the substrate adhesion dynamics affect how the intracellular force is converted cell crawling motion. Depending on the sign of the sensitivity of the substrate friction coefficient to the PIP3 concentration, the model cell exhibited crawling with the retrograde flow or with the direct flow. For the former case, which is consistent with experimental observations, our model showed that there is an optimum time delay and that the combined effect of the mechanical and chemical signals on the substrate friction coefficient can increase the migration distance. We also investigated the impact of the cell shape and the cell size, which led to qualitatively the same results. In addition, we included stochasticity in the RD equations, enabling the cell to change its migration direction and to change its dynamical mode from translational motion to spinning motion. Further, we performed multipole analysis of the substrate traction force, which was qualitatively consistent with the experimental results except the contractile nature of the traction force dipole.

Finally, we discuss some possible extensions of our current model.

Contraction process

In our current model, the protrusion and contraction processes are both modeled by the actuator elongation, $\ell^{\rm act}_{ij}(t)$, of the bond connecting two subcellular elements. The two processes are distinguished by the sign of $\ell^{\rm act}_{ij}(t)$. In this paper, however, we consider only the protrusion process, i.e., $\ell^{\rm act}_{ij}(t)\geq 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 include detailed dynamics of the concentration of adhesion molecules and their diffusion to other subcellular elements. Then, we can discuss detailed structures such as the 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 Fig. 3(a). However, real cells, e.g., Dictyostelium cells, tend to elongate in the direction of motion Maeda et al. (2008); Bosgraaf and Van Haastert (2009). One possible reason that our current model fails to reproduce this elongated shape is that actuator elongation depends only on the absolute value of $V_{i}$. Therefore, this inconsistency may be resolved by, e.g., introducing dependence on the gradient of $V_{i}$.

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.

X Conclusion

In conclusion, the modeling of crawling cells is still a challenging task due to the complexity of intra- and intercellular processes. The force that a cell generates should satisfy the force-free condition, where the total force vanishes. To achieve net migration under the force-free condition, the cell needs to break symmetry.

In our mechanochemical subcellular-element model, the intracellular force acts symmetrically on each pair of subcellular elements; therefore, it naturally satisfies the force-free condition. Symmetry breaking occurs due to the switching of the substrate friction coefficient between the adhered stick state and the deadhered slip state. Therefore, our model clearly distinguishes intracellular force and external force. To control those mechanical activities, we included RD equations representing intracellular chemical reactions.

The RD equations that we employed in this study were introduced to explain the chemical traveling wave observed within Dictyostelium cells Taniguchi et al. (2013). However, a number of chemical reactions occur inside a cell, and which chemical reactions are relevant may depend on the phenomena of interest. Nevertheless, we believe that our model can provide a basic framework for the future construction of mechanochemical models of crawling cells by replacing the RD equations with suitable ones for each specific phenomenon.

Acknowledgements.

This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (17H01083, 19H14673) grant, as well as the JSPS bilateral joint research projects.

Appendix A Simulation method

We integrate Eqs. (1)–(7) numerically by using the Euler method with time increment $\delta t=10^{-5}$. We add the stimulus $(\delta U,\delta V)=(-I_{\rm excite},+I_{\rm excite})$ to Eq. (5) on one of the subcellular elements, which we call the activator element. We set the activator element to the rightmost element denoted by the star in Fig. 1(a), except in the case of random motion, for which an activator element is chosen randomly every time $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_{\rm excite}=0.75$.

Appendix B Parameter values

The typical values of the mechanical parameters that appear in the model Eq. (1) are summarized in Table 1 and Table 2. The length of the cell at rest is set to $L_{\rm cell}=1$, and the total number of subcellular elements is set to $N=91$, which gives the total number of bonds as $N_{\rm bond}=240$, unless otherwise stated. For a larger cell, we set $L_{\rm cell}=1.4$, $N=169$, and $N_{\rm bond}=462$. For a cell consisting of a perfect hexagonal lattice, as in Fig. 1(a), the element distance at rest is set to $\ell_{ij}=L_{\rm cell}/\sqrt{N}$. Then, the mean bond length is $\ell_{0}=L_{\rm bond}/N_{\rm bond}=L_{\rm cell}/\sqrt{N}$, where $L_{\rm bond}$ is the total bond length. The elastic modulus of each bond is set to $\kappa=10$ and the intracellular dissipation rate to $\xi=L_{\rm cell}^{2}/L_{\rm bond}$, which gives the dimensionless parameter $\omega\xi_{\rm cell}/\kappa_{\rm cell}=0.2\pi L_{\rm cell}^{3}$. From the mechanical parameter values of typical living cells in Table 2, the value of this dimensionless parameter is estimated about 0.1–1. The substrate friction coefficient at the deadhered slip state is set to $\zeta_{\rm slip}=L_{\rm cell}^{2}/N$, so that the dissipation by the substrate friction is on the same order as the intracellular dissipation, $\omega\zeta_{\rm cell}/\kappa_{\rm cell}=\omega\xi_{\rm cell}/\kappa_{\rm cell}$. The substrate friction coefficient at the adhered stick state is set to be ten times larger than at the deadhered slip state, i.e., $\zeta_{\rm stick}=10\zeta_{\rm slip}$. The value of $\sigma$ in $U^{\rm area}$ is set to $\sigma=10^{-6}$.

Appendix C Intracellular RD equations

Here, we explain the RD equations proposed by Taniguchi et al. Taniguchi et al. (2013). Taniguchi et al. gathered from their experimental observations of Dictyostelium cells that phosphatidylinositol (3,4,5)-trisphosphate (PIP3) promotes actin polymerization and protrusion of the cellular membrane, and therefore, they considered a signaling pathway around PIP3 including phosphatidylinositol (4,5)-bisphosphate (PIP2), 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):

[TABLE]

where the reaction terms are defined as

[TABLE]

The global couplings are given by

[TABLE]

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.

Some experiments on Dictyostelium cells Taniguchi et al. (2013); Arai et al. (2010); Shibata et al. (2012) revealed that the actin filaments accumulate around the region where the PIP3 concentration increases, and thus, PIP3 promotes actin polymerization, which induces membrane protrusion. Therefore, we assumed that the actuator elongation depends on the PIP3 concentration $V$ as in Eq. (6).

The parameters we used in this paper are summarized in Table 3. Here, the diffusion coefficient $D_{U}=0.48$ corresponds to $0.8\mu\rm{m}^{2}\,(\rm{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 Fig. 7.

To integrate Eq. (17), 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

[TABLE]

where $n_{ij}=(n_{i}+n_{j})/2$, and $n_{i}=\sum_{j\neq i}w(r_{ij})$ is the number of neighboring elements around element $i$. The weight function that we employed is defined by

[TABLE]

where $r_{e}$ is the cutoff length, which is set to $4\ell_{0}$. The normalization factor $\lambda$ 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 $(\delta U,\delta V)=(-I_{\rm excite},+I_{\rm excite})$ to Eq. (17) on the activator element, assuming that the phosphorylation of PIP2 to PIP3 is enhanced for this element.

Appendix D Sinusoidal traveling wave

In this appendix, we show the cell crawling motion induced by a specified intracellular chemical traveling wave.

For simplicity, we consider only a single intracellular chemical component $c(t)$, which we assume travels periodically inside the cell. The chemical component for element $i$ changes as follows:

[TABLE]

where $c_{0}$ is the maximum concentration, which is set to $c_{0}=0.5$ so that $0\leq c_{i}(t)\leq 1$.

For simplicity, we assume a harmonic plane wave that travels in the $x$ direction. Then, the phase of the chemical concentration $\phi_{i}(t)$ changes as follows:

[TABLE]

where $\omega$ 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.

We assume that the cellular activity depends on the value of the intracellular chemical concentration $c_{i}(t)$. The actuator elongation depends on the mean chemical density:

[TABLE]

where $\ell_{c}$ represents the magnitude of the elongation. In addition, the substrate friction coefficient $\zeta_{i}(t)$ of element $i$ changes depending on $c_{i}(t)$. Here, we assume the substrate friction to be described by a two-state function that switches sharply between the adhered stick state and the deadhered slip state as follows:

[TABLE]

where $\psi_{\zeta}$ is a phase shift, $m_{i}$ is an integer that satisfies $2m_{i}\pi<\phi_{i}(t)-\psi_{\zeta}\leq 2(m_{i}+1)\pi$, and $\zeta_{\rm stick}$ and $\zeta_{\rm slip}$ are the values of the friction coefficient in the stick and slip states, respectively. Here, since the internal motion, i.e., the elongation and the stick-slip transition, is perfectly cyclic, we consider the phase shift instead of a time delay. The two are actually equivalent.

We carried out numerical simulations of Eq. (1) together with Eqs. (23)–(26) with the time increment $\delta t=10^{-4}$. In the simulation, we choose the element indicated by the star in Fig. 1(a) as the activator element, and we set $q_{w}>0$ so that the wave travels from right to left. The period of the traveling wave is set to the unit time scale, $T=1$, and thus, $\omega=2\pi$. We let the actuator elongate twice as long as its equivalent length: $\ell_{c}=\ell_{ij}$. For these parameters, we vary the wavenumber $q_{w}$ and the phase shift $\psi_{\zeta}$ and measure the migration distance for one cycle.

The results are summarized in Fig. 8. In Fig. 8(a), the migration distance for one cycle $\Delta R$ is plotted against the phase shift $\psi_{\zeta}$ for different values of $q_{w}$, and time series of a typical crawling cell are depicted in Fig. 8(b) and Fig. 8(c).

Fig. 8(a) shows clearly that the migration direction changes depending on the phase shift $\psi$. The positive $\Delta R$ represents forward motion, i.e., migration towards the activator element. That is, the cell migrates against the intracellular traveling wave, as depicted in Fig. 8(b). In contrast, the negative $\Delta R$ corresponds to backward motion, where the cell migrates in the same direction as the intracellular traveling wave, as shown in in Fig. 8(c). Therefore, there are maximum and minimum values of $\Delta R$, where the cell achieves maximum migration distance in the forward and backward directions, respectively.

The migration distance $\Delta R$ also depends on the wavenumber $q_{w}$. There is an optimal $q_{w}$ around $q_{w}=4\pi/5$ for which the cell can achieve the largest migration distance. That is, in this case, the cell is in a fully stretched state.

Appendix E Circular cell

For cells with a circular shape at rest, we prepared the subcellular elements and the connecting bonds as depicted in Fig. 9(a). The length of the cell was set to $L_{\rm 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_{\rm bond}=267$. We set $\ell_{ij}$ of each bond of the circular cell to the initial element distance, as shown in Fig. 9(a). Then, the mean bond length was calculated to be $\ell_{0}=0.105724$. The other parameters were kept the same as for hexagonal cells.

Fig. 9(b) shows an example of a time series of circular cell crawling obtained numerically from Eqs. (1)–(7). Here, the activator element is set as the one denoted by the star in Fig. 9(a), which is then stimulated by $(\delta U,\delta V)=(-I_{\rm excite},+I_{\rm excite})$ with $I_{\rm excite}=0.75$ every $t=1.5$.