Spiral actin-polymerization waves can generate amoeboidal cell crawling

The phase field is an auxiliary real-valued positive field defining the interior and exterior of the cell. The interior is given by $\psi \to 1$, whereas the exterior corresponds to $\psi \to 0$. The phase field continuously varies between these two values. The dynamics of the phase field is given by [27, 28, 31, 32]

$$\begin{eqnarray}&&{{\partial }_{t}}\psi =f\left( \psi ,\delta \right)+{{D}_{\psi }}\Delta \psi -\beta \mathbf{p}\cdot \nabla \psi .\end{eqnarray} \tag{ 1 }$$

The zeros of the function f determine the values of the pure phases, $\psi =0$ and $\psi =1$. Explicitly, we use the cubic form

$$\begin{eqnarray}&&f\left( \psi ,\delta \right)=\kappa \psi \left( 1-\psi \right)\left( \psi -\delta \right),\end{eqnarray} \tag{ 2 }$$

such that the phase field relaxes into the state $\psi =0$ for an initial value ${{\psi }_{0}}<\delta$ and into the state $\psi =1$ in the opposite case. The parameter κ determines the time scale on which the phase field reaches these values. Finally, the value of δ depends on the phase field itself through the constraint

$$\begin{eqnarray}&&\delta =\frac{1}{2}+\epsilon \left( \mathop{\int }^{}{{\text{d}}^{2}}\mathbf{r}\;\psi \left( \mathbf{r} \right)-{{V}_{0}} \right)\end{eqnarray} \tag{ 3 }$$

and serves to maintain a constant cell area ${{V}_{0}}$, see [28]. Here, the parameter $\epsilon >0$ gives the stiffness of the constraint. In principle, its value can depend on the membrane surface tension and elasticity. However, different values of essentially amount to rescaling ${{V}_{0}}$ and ${{D}_{\psi }}$, such that the exact value of is not important. Since cell motion takes place in an overdamped regime, we neglect inertial terms in the dynamic equations.

The effective diffusion term in equation (1) determines the width of the interface between phase $\psi =0$ and $\psi =1$, that is, the membrane and the value of the constant ${{D}_{\psi }}$ is related to the interfacial energy. For simplicity, we neglect a possible contribution of the membrane's bending stiffness to the phase-field dynamics. The third term in the right-hand side of equation (1) couples the phase field to the polymerization dynamics of actin filaments. Here, we assume that filaments pointing with their barbed end towards the cell boundary push the boundary further outwards, whereas filaments pointing with their pointed end towards the cell boundary pull it inwards. The action of molecular motors has indeed been suggested to induce pulling on the membrane [35]. For simplicity we consider a situation in which molecular motors are absent. We will show below that in our solutions the polarization field at the membrane always points outwards, and, hence, is responsible for pushing the membrane forward. The value of β determines the strength of the coupling between the filaments and membrane, and depends on the substrate's properties and the actin polymerization rate, see [30, 31] for details.

Actin assembly is an intricate process depending on the state of the nucleotide bound to an actin monomer and the presence of accessory proteins that can speed up or slow down the addition and removal of monomers at the filament ends [4]. We drastically simplify this dynamics and assume that the filaments grow on average at the barbed end and shrink at the pointed end. In the absence of accessory proteins, though, the average growth rate at the barbed end of an actin filament is lower than the average disassembly rate at the pointed end, such that the filaments shrink on average [36–38]. As a consequence, for example, capping proteins that prevent disassembly at the pointed end or proteins that act as an elongation factor at the barbed end are necessary to maintain net filament growth. To keep our system reasonably simple, we do not consider such proteins explicitly, but take the growth speed ${{v}_{\text{a}}}$ at the barbed end to be larger than the disassembly speed ${{v}_{\text{d}}}$ at the pointed end, ${{v}_{\text{a}}}>{{v}_{\text{d}}}$. To prevent filaments from growing indefinitely, we also introduce a rate ${{k}_{\text{d}}}$ at which filaments instantaneously degrade completely. Under physiological conditions, the spontaneous nucleation of actin filaments is a rare event [4] and we neglect it in the following.

These processes can be accounted for in a kinetic theory that describes the system in terms of a filament density depending on the filament position, orientation, and length [17]. As we show in appendix A, the explicit length dependence can be integrated out. Furthermore, by making a moment expansion with respect to the orientational degree of freedom, the system can be expressed in terms of the total actin density T and the average orientation p at a position r [39]. Here, the length of p scales with the filament density. The dynamic equations now read

$$\begin{eqnarray}&&{{\partial }_{t}}T=D\Delta T-{{k}_{\text{d}}}\psi T-{{v}_{\text{a}}}\psi \nabla \cdot \mathbf{p}+\alpha \psi {{n}_{\text{a}}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\partial }_{t}}\mathbf{p}=D\Delta \mathbf{p}-{{k}_{\text{d}}}\psi \mathbf{p}-{{v}_{\text{a}}}\psi \nabla T.\end{eqnarray} \tag{ 5 }$$

The respective first terms in the right-hand sides of equations (4) and (5) are effective diffusion terms capturing the effects of fluctuations in the system. The following terms result, respectively, from the filament degradation and polymerization. The last term on the right-hand side of equation (4) stems from the nucleation of new filaments by active nucleators. Here, α is a measure of the nucleation rate. The degradation and polymerization terms are confined to the cell interior by multiplication with the phase field ψ.

We do not describe the distribution of actin monomers explicitly and assume that they form a reservoir. In fact, the situation in living cells is somewhat similar as there are proteins that sequester actin monomers, such that the concentration of actin monomers that can be incorporated into filaments is roughly constant [4].

What remains is to specify the dynamics of the NPFs. We take their number to be conserved, and there is a spontaneous rate ${{\omega }_{\text{a}}}$ of the NPF activation. In addition, there is a rate ${{\omega }_{1}}{{n}_{\text{a}}}$ at which active NPFs activate inactive NPFs. We neglect spontaneous NPF inactivation, but consider inactivation mediated by actin filaments. There is experimental evidence for this process in human neutrophils [14] and D. discoideum [40]. The corresponding rate is given by ${{\omega }_{\text{d}}}T$, and the activation and inactivation processes are restricted to the cell interior. Finally, NPFs can diffuse with respective diffusion constants ${{D}_{\text{a}}}$ and ${{D}_{\text{i}}}$. As active NPFs are often linked to larger protein complexes or to the membrane, we take ${{D}_{\text{a}}}$ to be much smaller than ${{D}_{\text{i}}}$. Explicitly, dynamic equations for the NPF densities read:

$$\begin{eqnarray}&&{{\partial }_{t}}{{n}_{\text{a}}}={{D}_{\text{a}}}\Delta {{n}_{\text{a}}}+{{\omega }_{\text{a}}}\psi \left( 1+{{\omega }_{1}}n_{\text{a}}^{2} \right){{n}_{\text{i}}}-{{\omega }_{\text{d}}}\psi T{{n}_{\text{a}}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{n}_{\text{i}}}={{D}_{\text{i}}}\Delta {{n}_{\text{i}}}-{{\omega }_{\text{a}}}\psi \left( 1+{{\omega }_{1}}n_{\text{a}}^{2} \right){{n}_{\text{i}}}+{{\omega }_{\text{d}}}\psi T{{n}_{\text{a}}}.\end{eqnarray} \tag{ 7 }$$

This completes the specification of the model, which is also illustrated in figure 1. For our numerical solutions, we use a dimensionless form of dynamic equations (1)–(7), which is given in appendix B. We scale space by $\lambda =\sqrt{{{D}_{\text{i}}}\omega _{\text{a}}^{-1}}$ and time by $\tau =\omega _{\text{a}}^{-1}$. Dimensionless parameters are indicated by tildes.

Zoom In Zoom Out Reset image size

Since in equations (6) and (7) the diffusion terms are not directly coupled to the field ψ, NPFs can be transported outside the cell. For this reason, in our numerical integration scheme, NPFs having leaked out of the cell domain are removed in every time step by setting ${{n}_{\text{a}}}{{\left| {} \right.}_{\psi <\bar{\psi }}}=0$ and ${{n}_{\text{i}}}{{\left| {} \right.}_{\psi <\bar{\psi }}}=0$, where we choose $\bar{\psi }={{10}^{-2}}$. Since the total number of NPFs is a conserved quantity, we also add a term to the NPF densities compensating for their loss outside the cell, ${{n}_{\text{a}}}={{n}_{\text{a}}}+\psi {{N}_{\text{a}}}/\mathop{\int }^{}{{\text{d}}^{2}}\mathbf{r}\;\psi \left( \mathbf{r} \right)$ and ${{n}_{\text{i}}}={{n}_{\text{i}}}+\psi {{N}_{\text{i}}}/\mathop{\int }^{}{{\text{d}}^{2}}\mathbf{r}\;\psi \left( \mathbf{r} \right)$. Here, ${{N}_{\text{a}}}=\mathop{\int }^{}{{\text{d}}^{2}}\mathbf{r}\;{{n}_{\text{a}}}\left| _{\psi <\bar{\psi }} \right.$ and ${{N}_{\text{i}}}=\mathop{\int }^{}{{\text{d}}^{2}}\mathbf{r}\;{{n}_{\text{i}}}\left| _{\psi <\bar{\psi }} \right.$ are the respective amounts of degraded active and inactive NPFs outside the cell domain. We checked different ways for compensating for 'escaped' NPFs. The systems generated the states we report below practically independently of the method used. Similarly, the actin density T and orientation p can leak out of the cell domain. However, we will choose D = 0 in our simulations. As a consequence, this effect does not alter the solutions qualitatively and so we do not remove these fields explicitly.

The cytoskeleton can only generate force dipoles such that the total momentum flux through the plasma membrane of a cell must vanish, if external forces are absent. In the overdamped regime considered here, this amounts to the sum of all internal forces acting on the cytoskeleton being zero. In the case considered here, the cytoskeleton can exchange momentum with the environment through interactions with the membrane and with the substrate it moves on. In addition, there might be a friction force originating from interactions between the surrounding fluid and the cell, but this force is negligibly small. While in our description forces do not appear explicitly, we can still interpret our dynamic equations in this respect. To this end we assume that the cytoskeleton adheres tightly to the substrate and neglect any actin flux relative to the substrate.

Let us consider force balance in a situation close to a sharp interface limit, that is, ${{D}_{\psi }}$ is very small. At the membrane, three forces need to be considered. First of all, surface tension $\gamma \propto {{D}_{\psi }}$ leads to a force density ${{\mathbf{f}}_{\text{ten}}}$ along the membrane. Here, ${{\mathbf{f}}_{\text{ten}}}=\gamma \chi \mathbf{n}$, with χ being the membrane curvature and $\mathbf{n}$ denoting the outward normal $\mathbf{n}=-\nabla \phi /|\nabla \phi |$, see [27]. Dissipative (friction) forces are captured by $\xi \mathbf{v}$, where $\mathbf{v}$ is the (local) membrane velocity and ξ is the effective friction coefficient. Finally, there is the force ${{\mathbf{f}}_{\text{pol}}}$ exerted by the cytoskeleton on the membrane. This is given by ${{\mathbf{f}}_{\text{pol}}}=\xi \beta \mathbf{nn}\cdot \mathbf{p}$. We will see below in the numerical solutions that actin is only pushing on the membrane, that is, p is pointing outwards. This is in contrast to systems containing molecular motors generating contractile stresses in the cytoskeleton [28, 41]. Locally, we arrive at the following kinematic condition for the membrane

$$\begin{eqnarray}&&\xi \mathbf{v}={{\mathbf{f}}_{\text{ten}}}+{{\mathbf{f}}_{\text{pol}}}.\end{eqnarray} \tag{ 8 }$$

We now return to our discussion of the global force balance. As stated above, the sum over all external forces acting on the cytoskeleton must vanish. Under the assumptions made above, we only consider the traction and membrane forces on actin. Explicitly, we have

$$\begin{eqnarray}&&{{\mathbf{F}}_{\text{tract}}}+\mathop{\oint }^{}\text{d}\ell \ {{\mathbf{f}}_{\text{mem}}}=0,\end{eqnarray} \tag{ 9 }$$

where the integral is along the membrane and ${{\mathbf{f}}_{\text{mem}}}=-{{\mathbf{f}}_{\text{pol}}}={{\mathbf{f}}_{\text{ten}}}-\xi \mathbf{v}$. Let us emphasize that the traction forces are confined to the position of the membrane. This is because we do not consider any active stresses generated by the interaction of molecular motors with the actin network [31]. Indeed, the slime mold D. discoideum does not form focal adhesions and traction forces are essentially restricted to the outer boundaries [42].

For large enough values of ${{v}_{\text{a}}}$ and ${{\omega }_{\text{d}}}$, the system can spontaneously break the circular symmetry of the stationary state and evolve into a state of persistent motion, see figure 2 and supplementary movie S01, available at stacks.iop.org/njp/16/055007/mmedia. In such a state, the active NPFs localize at the leading edge and, hence, also the filament density T. The polarization field has topological defects at the maximum actin densities, see figure 3 and appendix B, such that the polarization field points into the direction of motion at the leading edge, whereas it points into the opposite direction at the trailing edge. However, in contrast to the states of persistent motion reported in [43], the fields do not form a solitary wave, but change incessantly: the NPF-rich region periodically dissolves and reforms closer to the cell center, from where it then moves towards the leading edge. Note that the density of active NPFs can have a varying number of maxima such that the number of point defects of the polarization field is not conserved.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These periodic changes in the distributions lead to oscillatory variations in the propulsion force, see figure 3(f), inset. As a result, the speed of the cell exhibits profound periodic oscillations similar to that observed in D. discoideum [44]. Periodic changes of the migration speed have also been observed in neutrophils [45], but there might be a consequence of $\text{C}{{\text{a}}^{2+}}$ oscillations. The average propulsion velocity is largely independent of the polymerization velocity ${{v}_{\text{a}}}$. Instead, it depends sensitively on the rates of actin nucleation and of nucleator binding and unbinding. This behavior can be rationalized by considering a dimensionless form of equations (1)–(7). As we show in appendix B, for the chosen parameters, the actin density T and polarization p are enslaved by the density ${{n}_{\text{a}}}$ of active nucleators. This property also suggests that details of the actin polymerization process are not important for the overall system behavior. Indeed, other polymerization dynamics in combination with the NPF dynamics similar to the one used here have been found to spontaneously produce waves [14, 16, 18, 19]. Even though the polymerization velocity does not have a major impact on the wave's propagation velocity or the cell's migration speed, which are set by the rates of nucleator activation and inactivation as well as the actin nucleation rate, changes in ${{v}_{\text{a}}}$ can lead to the qualitative changes in the type of motion, as we will discuss now.

Decreasing the growth velocity from ${{\tilde{v}\,}_{\text{a}}}=0.39$ to ${{\tilde{v}\,}_{\text{a}}}=0.27$, the state changes qualitatively: the fragment is still moving, however, its trajectory is no longer straight; see figure 4(a). Instead it consists of curved segments that are interspersed by points at which the cell discontinuously changes direction. We would like to emphasize that this apparently random motion results from deterministic dynamic equations (1)–(7) without addition of noise.

Zoom In Zoom Out Reset image size

We can characterize the trajectory by calculating its root mean squared displacement (RMSD) $\sqrt{\left\langle {{r}^{2}}\left( t \right)\right\rangle }={{\left( \underset{T\to \infty }{\mathop{\lim }}\,\;\frac{1}{T}\cdot \mathop{\int }^{}_{0}^{T}{{\left| \mathbf{r}\left( t\prime +t \right)-\mathbf{r}\left( t\prime \right) \right|}^{2}}\;\text{d}t\prime \right)}^{1/2}}$; see figure 4(b). On a short time scale $t\lesssim 1.7\omega _{\text{a}}^{-1}$, the RMSD grows essentially linearly with time. Then there is a cross-over to a region where the RMSD grows with an exponent close to but a little less than 1/2. The deviations from 1/2 might be due to finite size effects. We conclude that, within the accuracy of our calculations, the migration is diffusive in this state and can be described as a random walk. This conclusion is further supported by the velocity autocorrelation function, which decays rapidly. Furthermore, the curved segments have an equal probability of being curved in one or the other direction, making the random walk unbiased.

To obtain a better understanding of the processes underlying the formation of an irregularly moving state, let us consider the actin distribution T during an event of a discontinuous change in the direction of motion, see figure 5 and supplementary movie S02. Again the actin field is concentrated in a fraction of the cell domain. In contrast to actin 'blobs' in the gliding states, it now lacks a lateral symmetry axis, though. Instead it spirals along the cell membrane leading to curved segments of the trajectory (figure 5). Similar to the gliding states, in the irregularly moving state, too, the blobs have only a finite life-time. New blobs are then generated in the cell interior, which again leads to spiral motion. We could not detect any correlation between the orientations of two subsequent spirals. Consequently, the trajectory can change discontinuously. An additional feature of the dynamics is that a blob can break up into two or more parts, see figures 5(b), (c). This typically leads to changes in the curvature radius of the trajectories' curved segments.

Zoom In Zoom Out Reset image size

The actin blobs develop into spirals if the cell size ${{V}_{0}}$ is chosen to be sufficiently large, see below. We conclude that the irregular state essentially results from the intrinsic non-trivial spiral dynamics. Such dynamics is known from excitable systems with delayed inhibitor production [46] or more than two components [47], as well as from oscillatory media [33, 48].

To obtain a more systematic understanding of the emerging states, we construct the phase diagram as a function of two most crucial parameters: the nucleator detachment rate ${{\omega }_{\text{d}}}$ and the actin polymerization velocity ${{v}_{\text{a}}}$, see figure 6. Before describing the phase diagram let us comment on the dependence of the system on the nucleation rate α and the nucleator activation rate ${{\omega }_{\text{a}}}$. As the actin density follows the dynamics of the active nucleators closely, changing the value of α has a similar effect as changing the value of ${{\omega }_{\text{d}}}$. Indeed, in the dynamic equations (1)–(7) only the product ${{\omega }_{\text{d}}}T$ appears and $T\sim \alpha$. In contrast, increasing the value of ${{\omega }_{\text{a}}}$ has a similar effect as decreasing the value of ${{\omega }_{\text{d}}}$.

Zoom In Zoom Out Reset image size

We classify the states in terms of the actin dynamics inside the cell. As noted in the previous sections, the actin distribution can be stationary and axially symmetric, see figure 7(a), it can organize into spirals, or form blobs of high concentration having a finite lifetime. The spiral states can be further divided into stable and unstable spirals, meaning that they incessantly break up, disappear, and nucleate again. Two examples of the former are shown in figures 7(b), (c); see also supplementary movies S04 and S05. The blobs can either be axially symmetric as discussed above for the gliding state or not. In the non-symmetric case, they move along the cell boundary, see figures 7(d)–(f) and supplementary movie S03. The appearance and disappearance of blobs is periodic in time not only for the axially symmetric blobs as we had seen above, but apparently also in the general case. From our numerical solutions it is not clear if the transition from spirals to blobs moving along the membrane is a bifurcation or a continuous transition. Let us also note that we did not detect regions of co-existence between the different states.

Zoom In Zoom Out Reset image size

Different states of the internal actin dynamics correspond to different mobility patterns of the cell. In a stationary state, the cell's center of mass does obviously not move. In the case of stable spirals, the center of mass does not move either. The spiral arm moves along the membrane and induces a protrusion travelling along the membrane similar to observations made on spreading fibroblast [49, 50] and in drosophila cells with silenced Pak3, a homolog of p21-activated kinase [51]. Interestingly, Pak3 depleted cells show an increased actin-polymerization activity. In our case, however, protrusions are small and more work is needed to establish the relation between the polymerization waves in adhering cells and the ones studied here. Unstable spirals produce the irregular orbits presented above. Moving away from the spirals' stability threshold, the corresponding effective diffusion constant increases. For unstable blobs, the center of mass moves along trajectories reminiscent of Lissajous figures, see figure 7(f), inset. As the period necessary for circling along the cell perimeter approaches the lifetime of a blob, the trajectory approaches a circle. Finally, the axially symmetric blobs are associated with directional motion as discussed above.

The cell dynamics also depends on the typical cell volume ${{V}_{0}}$ if all other parameters are fixed, such that the dynamics of the nucleator and actin fields spontaneously generate waves. Below a critical value, the cell remains immobile in a stationary state. Beyond this value, the cell moves directionally as in the phase indicated by the blue squares in figure 6. As the cell size is increased further, the system settles first into the randomly moving dynamic state. Alternatively, the system can assume the state with pinned center of mass and fluctuation of the cell boundary. For large cell volumes, the actin density formed a well-developed spiral, see supplementary movie S06.

In addition to the physical parameters discussed so far, the phase field parameters also have an impact on the states generated by the system. An obvious effect of increasing the diffusion constant ${{D}_{\psi }}$ is a broadening of the interface, whereas the limit ${{D}_{\psi }}\to 0$ yields the sharp interface limit. In addition, larger values of ${{D}_{\psi }}$ also lead to a decreasing area of the cell: for the parameter values given in table B1, an increase of ${{\tilde{D}\,}_{\psi }}$ from 0.002 to 0.01, the area changed from 0.037 ${{\lambda }^{2}}$ to 0.03 ${{\lambda }^{2}}$. Importantly, though, we could not detect a change in the state due to changes in ${{D}_{\psi }}$.

This is different for changes in the value of β, which describes the coupling of the actin field to the interface and also depends on the actin polymerization rate. With increasing β we observe the same sequence of states as for increasing ${{v}_{\text{a}}}$ and ${{\omega }_{\text{d}}}$: unstable spirals turn into periodically appearing and vanishing blobs, first without and then with axial symmetry. Finally, for decreasing values of , fluctuations in the cell area around the value ${{V}_{0}}$ increase. Below a critical value of , the circular symmetric shape of the cell is stable. In this case, the cell volume varies periodically corresponding to a breather state of the actin network.