The dynamic equations

Amoeboid cell migration is driven by the actin cytoskeleton, which is mostly concentrated in the actin cortex, a layer beneath the plasma membrane. The cortex thickness is a few hundred nanometers [25–27] and thus much smaller than the lateral extension of a cell (>10 μm). In this work, we aim at describing the actin cytoskeleton adjacent to the substrate and thus use a two-dimensional geometry.

We use the continuum description of Refs [10, 16, 21] for the actin dynamics, where the actin density is captured by the field c. The alignement of actin filaments can lead to (local) orientational order in the system. This effect is captured by the orientational order parameter p, which is similar to the nematic order parameter of liquid crystals. In the dynamic equations, all terms allowed by symmetry up to linear order and up to first order in the derivatives are considered, such that

Here, v_a is the average polymerization speed and k_d an effective degradation rate, see Fig 1. Instead of the phenomenological approach used here, Eqs (1) and (2) can also be obtained by coarse-graining a kinetic description [21]. In this way, one sees that the term v_a ∇ p in Eq (1) describes changes of the actin density resulting form the addition or removal of actin monomers at the ends of actin filaments. The term v_a ∇ c indicates that changes in the polarization are linked to the actin polymerization current v_a c: when actin is polymerizing from an actin-dense region into an actin-sparse region, the polarization of the actin network grows in the direction opposite to the actin-density gradient. Note, that this description neglects flows of the actin network [28] that could, for example, be generated by molecular motors. We also neglect a possible diffusion term that would account for fluctuations in the actin dynamics. We have checked that our results are not affected qualitatively for sufficiently small diffusion constants.

Fig 1

Schematic representation of the actin dynamics captured by Eqs (1)–(4).

Blue circles represent inactive nucleators. They are spontaneously activated at rate ω₀, a process that is often associated with membrane binding. The activation rate is enhanced by already active nucleators, represented by green circles, which is captured by the parameter ω. Active nucleators generate new actin filaments (red) at rate α. The latter grow at velocity v_a and spontaneously disassemble at rate k_d. Furthermore, actin filaments attract factors that inactivate nucleators. This complex process, which can involve several different proteins in a cell, is captured by the rate ω_d.

The last term of Eq (1) is a source term that describes nucleation of new actin filaments. For the conditions present in cells, new actin filaments hardly form spontaneously. Instead, specialized proteins assist in this process. Examples are members of the formin family or the Arp2/3 complex. These proteins can be in an active or an inactive state and their spatial distribution in a cell can change with time. In this way they can contribute essentially to orchestrating the organization of the actin cytoskeleton. We introduce the densities n_i and n_a to describe these actin nucleation promoting factors—’nucleators’ for short -, where the indices refer to the inactive and active forms, respectively. Active nucleators generate new actin at a rate α, hence the form of the last term in Eq (1).

The dynamic equations for the fields n_a and n_i capture their transport by diffusion and their activation and inactivation dynamics. On the time scales that are relevant for the dynamics we study in the remainder of this work, nucleator synthesis and degradation can be neglected. Consequently, the dynamic equations should conserve the number of nucleating proteins, ∫_A(n_a + n_i)dA = An_tot = const, where A is the cell area adjacent to the substrate. We write

The diffusion constants for active and inactive nucleators are D_a and D_i, respectively. Spontaneous activation of nucleators occurs at rate ω₀. There is some experimental evidence for a positive feedback of nucleator activation [29], such that active nucleators promote the activation of further nucleators. Recently, it was reported that the molecular network underlying this positive feedback involves the Rho activating Guanine nucleotide exchange factor (GEF) GEF-H1 [30]. The small guanosine triphosphatase (GTPase) Rho in turn activates actin nucleating factors of the formin family. We capture this effect by the parameter ω. Nucleator deactivation can occur spontaneously. Furthermore, it has been proposed that nucleator deactivation can be induced by factors that are recruited by actin filaments [11, 12, 29, 31, 32]. We assume that the latter dominates [29] and neglect spontaneous deactivation. Actin induced deactivation is controlled by the parameter ω_d.

To fully determine the dynamics of the fields c, p, n_a, and n_i, Eqs (1)–(4) have to be complemented by boundary conditions. In this section, we use periodic boundary conditions to study the intrinsic actin dynamics. Later we will add the presence of the cell membrane through a phase field, see Sect Cell motility from actin polymerization waves.

In the following we use a non-dimensionalized version of the dynamic equations. We scale time by $$ {\omega }_0^{-1}$$ and space by $$ \sqrt{D_i/{\omega }_0}$$. We use the same notation for the rescaled parameters as in Eqs (1)–(4), such that the non-dimensionalization corresponds to setting ω₀ = 1 and D_i = 1. Unless noted otherwise, we use in the following the parameter values given in Table 1.

Table 1

Nondimensional parameter values used in this work unless indicated otherwise.

Parameter	Meaning	Value
Da	Diffusion constant of active nucleators	4 ⋅ 10−2
va	Effective actin polymerization speed	0.1–0.6
kd	Effective filament degradation rate	176
ω	Cooperative binding strength of nucleators	6 ⋅ 10−3
ωd	Detachment rate of active nucleators	0.1 − 0.6
α	Actin polymerization rate	588
ntot	Average total nucleator density	700
L	System length	1.3
Ng	Number of grid points per dimension	256
$$	Time scale	91.6 s
$$	Length scale	63.5 μm
DΨ	Phasefield relaxation / surface tension coefficient	5 ⋅ 10−3
κ	Phasefield timescale modifier	118
ϵ	Area conservation strength	8
β	Actin-membrane interaction coefficient	5.75 ⋅ 10−3
A0	Mean cell area	0.083

The length and time scales are chosen such that the ensuing dynamics is comparable to that of immature dendritic cells [10].

Spatially homogenous solutions

Consider the case of homogenous protein distributions. The constraint on the nucleator density thus is n_a + n_i = n_tot = const, where n_tot is the average total nucleator density. According to Eq (2), the polarization field is decoupled from the other fields and will tend to zero, p → 0, for t → ∞. The remaining dynamic equations become

Eqs (5) and (6) are reminiscent of the FitzHugh-Nagumo (FHN) system [33, 34]. In its general form, the latter is given by [35]:

Eq (7) describes generation of the ‘carrier’ w by the ‘driver’ v and degradation of w with rate a. Here, ϵ ≪ 1 is a small parameter, such that the dynamics of w occurs on longer time scales than the one of v. The second equation captures inhibition of v by w and I is an external stimulus. Finally, f(v) describes a feedback of v on its own production: in general, it promotes generation of v for small values of v, whereas it inhibits its production for larger values of v.

A typical specific choice of f is $$ f\left(v\right)=v-\frac{v^3}{3}$$. In that case, the system essentially depends only on the parameter a and the external stimulus I, because variations in ϵ do not affect the dynamics qualitatively as long as ϵ ≪ 1. Although the stimulus can depend on time, for the time being, we consider the case of constant I. Information about the asymptotic behavior can be obtained by analyzing the nullclines in phase space, that is, the curves defined by the respective conditions $$ \dot{v}=0$$ and $$ \dot{w}=0$$ in the (v, w)-plane. Intersections of the two nullclines correspond to fixpoints of which there are either one or three. In the latter case, the system is bistable as two fixpoints are stable against small perturbations, whereas the third is unstable, see Fig 2A.

Fig 2

Phase space diagrams for spatially homogenous dynamics.

A-C) Phase space for the FitzHugh-Nagumo Eqs (7) and (8) with a = 2, I = 0 (A), a = 0.04, I = 2 (B), and a = 0.4, I = 2 (C). D-F) Phase space for the dynamic Eqs (5) and (6) with k_d = 5, α = 50 (D), k_d = 50, α = 400 (E), and k_d = 80, α = 400 (F). Other parameters as in Table 1. In each case, the nullclines are shown in red, the vector fields as blue arrowheads and an example trajectory in black. For the FHN equations, the diagrams show a bistable case (A), a limit cycle (B) and an excitable case (C). For Eqs (5) and (6) we present limit cycles (D, E) and an excitable case (F). For these equations, there is no bistable case.

In the case that there is one fixpoint, it can be stable or unstable against small perturbations. If it is unstable, the system exhibits a limit cycle and asymptotically oscillates, see Fig 2B. In the opposite case, the FHN system can present excitable dynamics, that is, even though the fixpoint is stable against small perturbations, sufficiently large perturbations induce an ‘excursion’ in phase space, before returning to the fixpoint, see Fig 2C. This behavior can be observed, when the intersection of the two nullclines is left to the minimum or right to the maximum of the v-nullcline. If the intersection is between the two extrema, the system spontaneously oscillates, see Fig 2B.

The similarity between the actin-nucleator dynamics, Eqs (5) and (6), and the FHN system becomes evident when choosing c = w, n_a = v, ϵ = α, a = k_d/α, I = n_tot, and f(v) = −v + ωIv² − ωv³. The two dynamical systems differ in that the term −w of Eq (8) corresponds to −ω_d vw in Eq (6). Lastly, in contrast to v and w in the FHN system, which can take any real value, we now have w ≥ 0 and 0 ≤ v ≤ n_tot. Note that, in the FHN system, I is an external signal and can depend on time, while the corresponding term n_tot in the actin-nucleator system is a constant.

From the comparison between the actin-nucleator dynamics and the FHN system, we see that the actin-nucleator dynamics is driven by the nucleators, whereas actin is the carrier providing negative feedback. This is in agreement with experimental observations [17, 22]. The similarity between the two systems suggests that the actin-nucleator dynamics can also show oscillations as well as excitable behavior. This is indeed the case as we discuss now. We consider the case, where α is not a small parameter.

Let us take a closer look at the nullclines. Analogously to $$ \dot{w}=0$$ for the FHN system, $$ \dot{c}=0$$ yields a linear relation between c and n_a and the n_a-nullcline exhibits the characteristic S-shape of $$ \dot{v}=0$$. The nullclines of our system intersect exactly once in the region c ≥ 0 and n_a ≥ 0, such that there is only one fixpoint (c₀, n_a,0), independently of the parameter values. To see this, note first that the c-nullcline is a straight line through the origin. Now consider the function c(n_a) defined by the nullcline $$ {\dot{n}}_a$$. If there were parameter values for which three intersection points existed, then there would be some tangent to c(n_a) with a negative y-intercept c_y. However, for any value n_a ≥ 0 the value c_y is given by

If the fixpoint is unstable against small perturbations, the system exhibits oscillations as mentioned above, see Fig 2D and 2E. In case, (c₀, n_a,0) is stable, the system can amplify a finite perturbation, but will eventually return to the fixpoint, see Fig 2F. Before performing a linear stability analysis of the fixpoint, we first develop a physical picture of the necessary conditions for an instability.

The fixpoint can only be unstable, when the n_a-nullcline c(n_a) exhibits two extrema for n_a > 0. Explicitly, the nullcline is given by

Consequently, $$ \underset{n_a\to \infty }{lim}c\left(n_a\right)=-\infty $$ and $$ \underset{n_a\to 0^{+}}{lim}c\left(n_a\right)=+\infty $$. To determine whether the n_a-nullcline is monotonously decreasing, we consider the positive roots of the derivative c′ = ∂c/∂n_a. They are determined by

This equation always has one negative real solution. The two other roots are real only if the discriminant of the polynomial is negative. This leads to $$ \omega n_{\text{tot}}^2>27$$. In that case, two of the three real roots take the form

The value of $$ n_a^{+}$$ is always positive and $$ n_a^{-}$$ is always negative, because the argument of the sine function takes values between π/6 and π/2. The third root is

We now turn to a linear stability analysis of the fixpoint. For the dominating growth exponent s of the perturbation, we find

Wave solutions

After having analyzed the dynamic Eqs (1)–(4) for spatially homogenous fields, we now turn to the general case and study the system in a domain of size L² with periodic boundary conditions in the x- and y-direction. Then, the system can generate a variety of spatially heterogeneous solutions, including planar traveling waves and stationary patterns, see Fig 3 and S1 and S2 Videos. For information about the numerical approach we used for solving the dynamic equations, see Appendix: Numerical implementation of the dynamic equations. In the following we will determine the parameter region in which these patterns exist and characterize the shape of planar waves.

Fig 3

Snapshots of solutions for the actin concentration c to Eqs (1)–(4) in two dimensions with periodic boundary conditions.

A, B) Travelling planar waves for D_a = 0.04, ω_d = 0.28, v_a = 0.2 (A) and D_a = 0.04, ω_d = 0.32, v_a = 0.44 (B). Green arrows indicate the direction of motion. The disclinations in (B) might heal after very long times. C, D) Stationary Turing patterns for D_a = 0.04, ω_d = 0.45, v_a = 6.0 (C) and D_a = 0.21, ω_d = 0.42, v_a = 9.5 (D). For different initial conditions a pure hexagonal pattern of blobs can appear. All other parameters as in Table 1.

Linear stability analysis

We start our analysis by investigating the stability of the homogenous steady state against small spatially heterogeneous perturbations. The homogenous state is characterized by c(x) = c₀ = αn_a/k_d, p(x) = p₀ = 0, and n_i,0 = n_tot − n_a,0 with

As shown above there is only one positive solution n_a,0 ≤ n_tot to this equation, such that there is a unique homogenous stationary state.

Consider c(x, y, t) = c₀+ δc(x, y, t) and similarly for the fields p, n_a, and n_i. Linearizing the dynamic equations with respect to the steady state and expressing the perturbations in terms of a Fourier series, $$ \delta c={\sum }_{n,m=-\infty }^{\infty }{\widehat{c}}_{nm}{e}^{-i(q_{x,n}x+q_{y,m}y)}$$ and similarly for δ p, δn_a, and δn_i with q_x,n = 2πn/L and q_y,m = 2πm/L, leads to

The solutions to these equations are of the form $$ \widehat{c}\propto {\mathrm{e}}^{s_{nm}t}$$ etc, where s_nm are the growth exponents of the modes (n, m). If s_nm > 0, then a heterogeneous steady state emerges. If instead, ℜ(s_nm)>0 and ℑ(s_nm)≠0, then an oscillatory state, that is, either a standing or a traveling wave, can be expected.

Our numerical solutions indicate that all instabilities in our system are super-critical such that there is no coexistence of different states that are not linked by a symmetry transformation. Close to the instability, the wavelength λ₀ of the unstable mode determines the wave length of the emerging pattern. This remains true in a large region beyond the instability, see Fig 4. The wave length depends only weakly on the actin assembly velocity v_a, Fig 4A and 4B, and not on the nucleator inactivation rate ω_d, Fig 4D and 4E. It increases with the diffusion constant D_a, Fig 4C, and decreases with the cooperativity parameter ω, Fig 4F.

Fig 4

Wavelength as a function of system parameters.

Orange dots represent values obtained from numerical solutions in two spatial dimensions with periodic boundary conditions (L = 1.0), blue lines are the results of a linear stability analysis, see Sect Linear stability analysis. Parameter values are ω_d = 0.44 (A), ω_d = 0.48 (B), v_a = 0.46 and ω_d = 0.43 (C), v_a = 0.32 (D), v_a = 0.48 (E), and v_a = 0.46 and ω_d = 0.43 (F). All other parameters as in Table 1.

In contrast to the wave length, we only get a poor estimate of the wave’s propagation velocity from the linear stability analysis. In the following we use a variational ansatz to determine the wave form and propagation velocity of plane waves.

Wave form

We start by rewriting the dynamic Eqs (1)–(4). First of all, we combine the equations for the actin density c and the polarization p to obtain one equation for the density. Furthermore, we exchange n_i for N = n_a+ n_i. Finally, we consider solutions in a reference frame moving with the wave velocity v. We use periodic boundary conditions with period Λ and thus arrive at

where we have scaled space by Λ, such that the period is equal to 1, see see Appendix: Wave profile.

Eqs (20) and (21) are linear and can be solved as soon as n_a is known, see Dependence of migration characteristics on parameter values. To solve the nonlinear Eq (22) we make the following ansatz for a right-moving wave in the interval [−1/2, 1/2]

where a₁ to a₄ are variational parameters. We constrain a₂ and a₃ to vary in the intervals [5, 15] and [30, 50], respectively, whereas a₄ can take on the values 2, 3, 4; we do not impose any constraints on a₁. Note that the test function (23) does not fulfill the periodic boundary condition. However, since a₂, a₃ ≫ 1, n_a(a₁, a₂, a₃, a₄, ±1/2) ≈ 0.

In our ansatz, the active-nucleator density n_a increases according to the exponential polynomial $$ x^{a_4}{\mathrm{e}}^{a_{2}x}$$ at the front of the wave. In this region actin is nucleated and increases correspondingly. The trailing region of the wave is defined by a decrease of the active nucleator density according to 1+ tanh(a₃ x). This decrease results from a threshold actin concentration beyond which nucleator inactivation occurs at a higher rate than nucleator activation. Due to the large value of a₃, the nucleator density drops sharply to zero and also the actin density decays exponentially in the trailing region. The corresponding decay length is v/k_d, see see Appendix: Wave profile.

After solving the linear Eqs (20) and (21), we calculate an error by integrating the difference between the left and the right hand sides of (22) over the whole period:

Minimizing the error yields values for the variational parameters a₁ to a₄, v, and Λ.

In Fig 5A, we compare a solution obtained by the variational ansatz and by numerically solving the dynamic Eqs (1)–(4). The agreement is very good with the largest deviations being present at the front of the wave. Similarly, the parameter dependence of the wave speed is reproduced well by our variational ansatz, Fig 5B and 5C. The wave speed is essentially independent of the actin polymerization speed v_a as long as $$ v_a\lesssim 1$$, which is consistent with our earlier remark that the wave dynamics is driven by the nucleator activity rather than actin assembly. Furthermore, the wave speed increases with the parameter ω_d describing nucleator inactivation by actin. Indeed, as ω_d increases, nucleators are more rapidly inactivated, such that they become available for activation at the wave front.

Fig 5

Shape and velocity of traveling waves in one dimension.

A) Actin and active nucleator concentrations c (green) and n_a (orange) for v_a = 0.8 and ω_d = 0.35. Dots are from a numerical solution, solid lines are obtained from the variational ansatz Eq (23) and the solution Eq (39). B, C) Wave speed as a function of the actin polymerization velocity v_a (B) and the nucleator inactivation parameter ω_d (C). All other parameters as in Table 1.

Phase-field dynamics

Similar to previous work on cell motility, we use a phase-field approach to define the dynamic cell shape [23, 24]. A phase field is an auxiliary scalar field with values ranging between 0 and 1, which are called the pure phases of the system. We treat values of 0 as being outside of the cell and values of 1 as being inside. The phase-field dynamics is given by [23, 24]

The term proportional to κ derives from a free energy with minima at the pure phases. They are separated by an energy barrier at Ψ = δ. Conservation of the cell area can be achieved by adjusting the value of δ as described in Eq (26): The actual cell area is given by ∫_A ΨdA′, it’s target area by A₀. If the cell is bigger than A₀, then δ > 0.5 such that the overall cell area shrinks and vice versa. For sufficiently large values of κ, the transition between the two pure phases is sharp.

The transition region between the two pure phases determines the position of the cell membrane. Specifically, we implicitly define the location of the cell membrane by all positions r with Ψ(r) = 0.5. The term proportional to D_Ψ accounts for interfacial tension between the two pure phases and thus the surface tension of the membrane. For cells, surface tension of the membrane dominates its bending energy [24], which we neglect. Finally, the term proportional to β describes the interaction strength between the phase field and the actin network. The interaction is always directed along the polarization vector, such that the membrane can be pushed outwards or pulled inwards [24]. In our solutions we do not observe pulling to the inside.

The dynamics of the actin network and the nucleators is confined to the cell interior by multiplying the dynamic Eqs (1)–(4) by Ψ. Conservation of the nucleators is an important aspect of these dynamic equations. Simply multiplying the corresponding transport term by Ψ violates conservation of the total nucleator amount and also leads to nucleators leaking out of the cell interior [21]. Here, we choose a different option and instead modify the nucleator current at the position of the membrane. For a particle density n, we write

This term evidently conserves the total particle number. It can be interpreted as a combination of scaling the diffusion constant with Ψ and introducing an inwards flux proportional to D at the membrane. This suggests that the expression is efficient for keeping the nucleators inside the cell. This is indeed the case as can be seen by solving for the stationary state of Eq (28), which is given by n ∝ Ψ.

In this context, it is also instructive to look at the discretized version of the right hand side of Eq (28). Using the discretized Laplacian △n_j ≡ (n_j−1 − 2n_j + n_j+1)h⁻², where h is the discretization length, we get in one dimension:

From this expression it is evident that nucleators can hop only to a site j inside the cell, i.e., with Ψ_j > 0, see Fig 6A.

Fig 6

Polymerization waves in presence of a phase field.

A) Schematic comparison of the discretized diffusion in absence (left) and presence (fight) of a phase-field, see Eq (29). B) Phase diagram of migration patterns as a function of the actin growth velocity v_a and nucleator inactivation parameter ω_d. C-E) Example trajectories with cell outlines drawn at 8 equidistant points in time for v_a = 0.34 and ω_d = 0.45 (diffusive migration, C), v_a = 0.22 and ω_d = 0.38 (random walk with straight segments, D), and v_a = 0.46 and ω_d = 0.43 (random walk with curved segments, E). Scale bars correspond to a length of 0.3. Other parameters as in Table 1.

In presence of the phase field, the dynamic equations are

For actin, the diffusion current can be neglected as argued above, such that its dynamics is unaffected by the modified diffusion introduced in Eq (28). The coupling of the actin density c and the polarization field p to the phase field is thus obtained by simply restricting the corresponding sources to the cell interior through multiplication with Ψ. However, since the actin concentration is not a conserved quantity and rapidly degraded in the absence of nucleators, we chose the degradation term to act also outside the cell interior to get rid of any actin that might have left the cell.

Actin-wave induced cell trajectories

In Fig 6B we show the phase diagram of the different dynamic patterns of the phase field’s center r_c = ∫rΨ(r)d² r as a function of the parameters v_a and ω_d. Five different dynamic states can be distinguished. Below a critical value of ω_d, waves do not emerge in the system and the center settles into a stationary state. The critical value of ω_d depends only weakly on v_a. There is a second critical value, such that the center r_c is again stationary if ω_d is larger than this critical value.

Close to the critical values of ω_d, the actin-nucleator system forms a spiral wave, see S3 Video. These spirals are symmetric and do not deform the phase field. They spin around a fixed point, which coincides with the center r_c. Since the dynamic equations are isotropic, solutions with clockwise or counter-clockwise rotations coexist. As the value of ω_d is, respectively, further increased or decreased, the spiral loses its symmetry. In this case, the motion of the center r_c becomes erratic and can be described as a random walk.

Three different types of random walks can be identified. First, the center r_c can exhibit diffusive dynamics, see Fig 6C and S4 Video. Second, it can perform a random walk, where straight segments along which the cell moves with constant velocity alternate with segments of diffusive motion, see Fig 6D and S5 Video. Also in the third type of random walk the cell center changes between two states, namely, diffusive or curved motion, see Fig 6E and S6 Video. Along the curved segments, the radius of curvature typically varies, but there are special cases, for which the radius of curvature along the curved segments is constant and the same for all segments. Note that for all kinds of random walk trajectories, the direction of motion after a diffusive segment is uniformly distributed. Similarly, the handedness of a curved segment is uncorrelated with that of the preceding segment.

For the erratic trajectories, the actin-nucleator dynamics is chaotic. For the persistent random walk, axisymmetric waves emanate from a center with a fixed position within the cellular domain. During the diffusive states, we observe spiral wave chaos instead. In the states corresponding to curved segments, the waves are not axisymmetric, which leads to ‘protrusions’ of the membrane and a turning of the cell axis. In case of the diffusive trajectories, the actin-nucleator dynamics exhibits spiral chaos. The deterministic dynamic equations are thus able to replicate salient migration features of searching cells [10].

Dependence of migration characteristics on parameter values

The random walks discussed above fall into the class of persistent random walks. For a persistent random walk, the velocity of the walker has a finite time autocorrelation, that is, its magnitude and direction persist for a characteristic time τ. Note that there are several realizations of a persistent random walk. In a run-and-tumble process, the walker exhibits periods during which it moves along straight lines with constant speed. These periods are interrupted by events during which the walker essentially does not move but changes its direction. Another possibility is that the direction of motion and the speed varies constantly in a smooth way. Inbetween these extremes, the segments of a run-and-tumble motion shows continuous changes of the velocity. In all cases, the mean square displacement 〈r²(t)〉 is given by 〈r²(t)〉 = 4Dt+ 2(vτ)²(e^−t/τ − 1). Here, v is the mean velocity of the persistent period and D is the diffusion constant describing the effective diffusive behavior on very long time scales. In the following we study, how the effective parameters τ, v, and D depend on our system parameters.

As shown in Fig 7A–7C, the persistence time τ, the speed v and the diffusion constant D initially increase with v_a and then decrease for larger values of v_a. The non-monotonous behavior of these quantities is a consequence of two competing effects. To see this, let us first recall that the wave speed does not increase with increasing v_a, Fig 5B. However, the polarization of the actin network does increase in this case as can be read of directly from Eq (2). Consequently, the interaction between the actin field and the membrane gets stronger and the membrane deformations are more pronounced. At the same time, the pronounced membrane deformations feed back on the actin waves, which are getting less regular. Thus, the cell polarization is less efficient, such that the periods of persistent migration are effectively shortened. At the same time, the migration speed decreases during these periods. This is confirmed by the mean instantaneous speed of the cell centers, which are very similar to the effective speed v, see Fig 7B.

Fig 7

Effective parameters of random walk trajectories.

A-C) Diffusion constant D (A), speed v (B), and persistence time τ (C) as a function of the actin polymerization speed v_a. D-F) As (A-C), but as a function of the nucleator inactivation parameter ω_d. Values were measured by fitting a persistent random walk model to the mean square displacement (MSD) of the respective trajectories. In (B) and (E), also the mean speed measured directly on the trajectories is shown (orange squares). Other parameters as in Table 1.

As a function of the parameter ω_d, we observe a transition from a persistent to a diffusive random walk. Below the transition, the parameters v and D increase with ω_d. In contrast, the value of τ depends non-monotonically on ω_d; it first increases and then decreases. Above a critical value of ω_d, we find τ = 0. For these values, the diffusion constant varies only slightly with ω_d and is two orders of magnitude smaller than for the persistent random walks. The dependence of v on ω_d is linear for the persistent random walks. Note that the values of v obtained from fitting the mean square displacement for τ ≈ 0 are not meaningful. The mean instantaneous velocity is again very similar to v for the persistent random walks. In the diffusive regime, it still grows linearly with ω_d. This is in line with the wave velocity, which increases with ω_d, see Fig 5C.

Actin density

The actin density c and the polarisation field p are given by Eqs (1)–(4), which in one spatial dimension and after non-dimensionalization read

Deriving Eq (34) with respect to time, we can eliminate the field p and obtain a linear equation for c with an inhomogeneity proportional to n_a:

This is the equation for a wave with speed v_a, internal friction with 2k_d and a driving proportional to $$ k_d^2$$. The source of the wave depends on n_a and its time derivative. We will assume that the active nucleators move as a solitary wave with velocity v, that is, n_a(x, t) = n(x − vt).

In the reference frame moving with the nucleation wave speed v and normalized by the wavelength L, Eq (36) becomes

The homogenous solution c_h(x) to this equation can be written as

The solution to the in-homogenous Eq (37) with the source term $$ S\left(x\right)=\frac{\alpha k_d{L}^2}{v^2-v_a^2}(1-\frac{v}{k_{d}L}{\partial }_x)n_a$$ is obtained by the method of variation of constants. We write c₀ = AS(x) and v₀ = AS′(x), where A is the Wronskian of our system and arrive at the full solution

The solution corresponds to a fraction of $$ \frac{v-v_a}{2v}$$ of the scaled nucleator density decaying on a lengthscale of $$ L_{-}=\frac{1}{\lambda -{\lambda }_a}$$ and a fraction of $$ \frac{v+v_a}{2v}$$ decaying with $$ L_{+}=\frac{1}{\lambda +{\lambda }_a}$$. The decaying part of the actin wave can be fitted perfectly with the single parameter $$ a(\frac{v-v_a}{2v}{e}^{-{\lambda }_{-}x}+\frac{v+v_a}{2v}{e}^{-{\lambda }_{+}x})$$. Note that the nucleation rate α has no effect on the shape of the wave, but only affects its amplitude.

The solution for the polarization field p is obtained by solving Eq (35) for p(x, t)≡p(x − vt).

Total nucleator density

We now rewrite the dynamic Eqs (1)–(4) for the active and inactive nucleator concentrations n_a and n_i in terms of the total nucleator concentration N = n_a+ n_i and n_a. In one spatial dimension and after non-dimensionalization, we have

In the reference frame of the traveling wave, (41) becomes

Integrating once and determining the integration constant by integrating once more over the entire system, we arrive at a first order equation for the total amount of nucleators,

Eq (43) implies that with a homogenous total nucleator concentration N = const = n_tot, gradients in n_a also vanish. Thus, a heterogeneity in the total nucleator concentrations is necessary to observe waves and wave propagation requires nucleator transport.

Furthermore, D_a is a measure for how far active nucleators can diffuse around the bulk of the wave while bound before detaching, on a time scale proportional to the wave period τ, thus affecting the wave length. D_a needs to be sufficiently smaller than D_i to create a length scale difference large enough to enable the formation of the bulk of the wave and maintain the imbalance in total nucleator concentration, otherwise the constant distribution of proteins is the only solution (as the wave length grows too large, or the imbalance shrinks too much to be supported).

The solution to Eq (43) is given by

From this equation we see that there are no waves, when D_a = 1 (= D_i).

Active nucleator density

Using the solutions for c, Eq (39), and N, Eq (44), we arrive at a single equation for the distribution of the active nucleators in the reference frame moving at the wave speed v: