Excitable actin dynamics and amoeboid cell migration

Amoeboid cell migration is characterized by frequent changes of the direction of motion and resembles a persistent random walk on long time scales. Although it is well known that cell migration is typically driven by the actin cytoskeleton, the cause of this migratory behavior remains poorly understood. We analyze the spontaneous dynamics of actin assembly due to nucleation promoting factors, where actin filaments lead to an inactivation of these factors. We show that this system exhibits excitable dynamics and can spontaneously generate waves, which we analyze in detail. By using a phase-field approach, we show that these waves can generate cellular random walks. We explore how the characteristics of these persistent random walks depend on the parameters governing the actin-nucleator dynamics. In particular, we find that the effective diffusion constant and the persistence time depend strongly on the speed of filament assembly and the rate of nucleator inactivation. Our findings point to a deterministic origin of the random walk behavior and suggest that cells could adapt their migration pattern by modifying the pool of available actin.


Introduction
The ability of cells to migrate is one of their most fascinating characteristics. During mesenchymal migration, cells persistently polarize and adhere to the substrate, which leads to persistent directional motion [1,2]. In contrast, during amoeboid migration, cells frequently change their polarization and hence their direction of motion. They also adhere less strongly to the substrate than cells during mesenchymal migration. Amoeboid migration can be observed for the soil amoeba Dictyostelium discoideum and for immune cells, for example, dendritic cells. The random walk performed during amoeboid migration is an important aspect of immune cells' task to scan the organism for pathogens. The origin of the random polarization changes during amoeboid migration is largely unknown [3] and it is not clear to what extent cells can control the characteristics of their random walk. Molecular noise is an obvious candidate for generating random migration [4,5]. The processes involved in generating migration are indeed subject to noise due to the stochastic nature of molecular reactions. However, these stochastic events take place on length and time scales that are small compared to those characteristic of cellular random walks. It is not obvious how cells could influence the strength of this noise and hence their migration behavior. Fluctuating external cues could also generate random walks. Indeed, cells respond to a multitude of external signals, notably, chemical or mechanical gradients, and adapt their migration accordingly.
Here, the cells have a certain degree of control as they can tune the strength of their responses. However, cellular random walks have been observed in the absence of external cues [6][7][8][9]. Finally, there is the possibility that cells generate internal polarization cues, which would give them the maximal possible control over their behavior. In this context, spontaneous actin polymerization waves have been proposed to provide such internal cues [10]. Actin is an important constituent of the cytoskeleton, which drives cell migration. It assembles into linear filaments-called F-actin-with two structurally different ends. This structural polarity of actin filaments is exploited by molecular motors that transform the chemical energy released during hydrolysis of adenosine-triphosphate (ATP) into mechanical work. The assembly and disassembly of F-actin is regulated by various cofactors. For example, formins and the Arp2/3 complex nucleate new filaments. Actin depolymerizing factor (ADF)/cofilin, on the other hand, can promote their disassembly. Interestingly, there is evidence for feedback between the actin cytoskeleton and the activity of these regulatory cofactors. For example, nucleation promoting proteins have been reported to be less active in regions of high F-actin density [11,12]. Such a feedback can lead to spontaneous actin polymerization waves [13][14][15][16][17]. Such waves are present during migration [13,15,18,19], and theoretical analysis has shown that they can be sufficient to cause cell motility [10,18,20,21].
From a physical point of view, spontaneous actin polymerization waves are akin to waves in excitable media. Early indications of this connection were given in [13,14,17]. Further support came from the observation that actin polymerization waves exhibit a refractory period [15,22]. More recently, the actin cytoskeleton of D. discoideum was shown to be poised close to an oscillatory instability [19]. The dynamics of excitable systems is exemplified by the Fitz-Hugh-Nagumo system, which is a very much simplified version of the Hodgkin-Huxley equations describing action potentials traveling along the axons of nerve cells.
In this work, we analyze the description of actin polymerization waves proposed in Ref [16]. We clarify its connection to the FitzHugh-Nagumo system and characterize the waves it generates. Furthermore, we use a phase-field approach [23,24] to study the impact of actin polymerization waves on cell migration. Here, the phase field is an auxiliary field that distinguishes between the inside and outside of a cell. We analyze in detail a recently introduced current for confining proteins to the cell interior [10]. Finally, we explore the relation between the system parameters and the characteristics of the random walks generated by chaotic polymerization waves.

Actin dynamics
In this section, we present the description of the actin cytoskeleton developed in Refs [10,16,21]. In Ref [16], the basic mechanism for the actin-nucleator dynamics, see Eqs (1)-(4) below, was presented and studied in fixed geometries. The coupling to a dynamic phase field, which represents the cell interior, was introduced in Ref [21]. There, the nucleator current in presence of a phase field had a form that led to strong leaking of nucleators from the cell interior. This was remedied in Ref [10], which focused on experiments and lacked a detailed study of the dynamic equations, which is the purpose of the present article. After establishing the dynamic equations, we discuss their relation to the FitzHugh-Nagumo model (FHN) and show that oscillations and waves emerge spontaneously in our system. Finally, we characterize the shape, length and propagation velocity of these waves.

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][26][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 r 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 r c indicates that changes in the polarization are linked to the actin polymerization current v a c: when actin is polymerizing from an actindense 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.
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 . Blue circles represent inactive nucleators. They are spontaneously activated at rate ω 0 , 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 . https://doi.org/10.1371/journal.pone.0246311.g001 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, R 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 ω 0 . 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 o À 1 0 and space by ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi D i =o 0 p . We use the same notation for the rescaled parameters as in Eqs (1)-(4), such that the non-dimensionalization corresponds to setting ω 0 = 1 and D i = 1. Unless noted otherwise, we use in the following the parameter values given in Table 1.

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 ! 1. The remaining dynamic equations become where we have used n i = n tot − n a .
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 ðvÞ ¼ v À 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 condi- 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. 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 3 . 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.  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. https://doi.org/10.1371/journal.pone.0246311.g002 Let us take a closer look at the nullclines. Analogously to _ w ¼ 0 for the FHN system, _ c ¼ 0 yields a linear relation between c and n a and the n a -nullcline exhibits the characteristic S-shape of _ 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 0 , 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 _ 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 which is always positive as the number of active nucleators is bounded from above by the total number of nucleators, n tot � n a . This proves the above statement.
If the fixpoint is unstable against small perturbations, the system exhibits oscillations as mentioned above, see Fig 2D and 2E. In case, (c 0 , 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, lim n a !1 cðn a Þ ¼ À 1 and lim n a !0 þ cðn a Þ ¼ þ1. To determine whether the n a -nullcline is monotonously decreasing, we consider the positive roots of the derivative c 0 = @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 on 2 tot > 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 which is always positive. In conclusion, the fixpoint (c 0 , n a,0 ) is unstable and the system oscillates for on 2 tot > 27 and if n 0 a < n a;0 < n þ a . We now turn to a linear stability analysis of the fixpoint. For the dominating growth exponent s of the perturbation, we find ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ða À k d Þ 2 À 4ao d n a;0 q 2 ; where a ¼ À 1 À 3on 2 a;0 þ 2on tot n a;0 À o d c 0 only depends on k d /α. By increasing the nucleation rate α while keeping k d /α = const, the nullcline remains unaffected. For k d > a the real part of the eigenvalue becomes negative, leading to a stationary state. Thus, k d < a is the last condition for the presence of oscillations in our system. The oscillation frequency ω F close to the instability can be estimated from the imaginary part of the growth exponent s, which gives o F ¼ IðsÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ao d n a;0 p .

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 2 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  Table 1.
https://doi.org/10.1371/journal.pone.0246311.g003 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. 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 0 = αn a /k d , p(x) = p 0 = 0, and n i,0 = n tot − n a,0 with ð1 þ on 2 a;0 Þn i;0 À o d c 0 n a;0 ¼ 0: 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 0 + δ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, dc ¼ P 1 n;m¼À 1ĉ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 d dtĉ nm ¼ À iv a ðq x;npx;nm þ q y;mpy;nm Þ À k dĉnm þ an a;nm The solutions to these equations are of the formĉ / 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 I(s nm )6 ¼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 λ 0 of the unstable mode determines the wave length of the emerging pattern. This remains true in a large region beyond the instability, see  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. 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] n a ða 1 ; a 2 ; a 3 ; a 4 ; xÞ ¼ a 1 2 e À a 2 x ð1 þ tanh ½a 3 x�Þð1 À 2xÞ where a 1 to a 4 are variational parameters. We constrain a 2 and a 3 to vary in the intervals [5,15] and [30,50], respectively, whereas a 4 can take on the values 2, 3, 4; we do not impose any  Table 1. https://doi.org/10.1371/journal.pone.0246311.g004 constraints on a 1 . Note that the test function (23) does not fulfill the periodic boundary condition. However, since a 2 , a 3 � 1, n a (a 1 , a 2 , a 3 , a 4 , ±1/2) � 0. In our ansatz, the active-nucleator density n a increases according to the exponential polynomial x a 4 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 3 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 3 , 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 1 to a 4 , 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 ≲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.
Stationary patterns. In addition to planar traveling waves, the dynamic Eqs (1)-(4) can also produce stationary patterns, see Fig 3C and 3D. These Turing patterns appear if v a ≳ 1 and consist either of 'blobs' of high or low active nucleator densities or of labyrinthine stripes of high active nucleator density. These structures can coexist in the same system. Since our focus in this work is on actin waves, we refrain from discussing these states further.

Cell motility from actin polymerization waves
Having analyzed the intrinsic actin dynamics, we now turn to a characterization of cell migration patterns emerging form spontaneous actin waves. We start by introducing a phase-field approach for describing the cellular domain. It contains a novel current for confining the nucleators to the cell interior. We then describe migration patterns and study the dependence of their characteristics on the system parameters.

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 C = δ. Conservation of the cell area can be achieved  Table 1. by adjusting the value of δ as described in Eq (26): The actual cell area is given by R A CdA 0 , it's target area by A 0 . If the cell is bigger than A 0 , 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 C(r) = 0.5. The term proportional to D C 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 C. Conservation of the nucleators is an important aspect of these dynamic equations. Simply multiplying the corresponding transport term by C 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 @ t n ¼ DðCΔn À nΔCÞ ð27Þ

¼ rðDCrnÞ À rðDnrCÞ: ð28Þ
This term evidently conserves the total particle number. It can be interpreted as a combination of scaling the diffusion constant with C 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 / C.
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 4n j � (n j−1 − 2n j + n j+1 )h −2 , 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 C j > 0, see Fig 6A. 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  Table 1. https://doi.org/10.1371/journal.pone.0246311.g006

PLOS ONE
corresponding sources to the cell interior through multiplication with C. 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 rC(r)d 2 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 actinnucleator 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 hr 2 (t)i is given by hr 2 (t)i = 4Dt+ 2(vτ) 2 (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  Table 1.
https://doi.org/10.1371/journal.pone.0246311.g007 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. 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.

Discussion
In this work, we have shown that a deterministic, self-organized system describing the actin assembly dynamics in living cells is capable of generating cellular random walks akin to amoeboid migration [10,21]. We elucidated its relation to excitable systems by a comparison with the FitzHugh-Nagumo system and characterized in detail spontaneously emerging traveling waves. We recall that the wave propagation speed is independent of the actin polymerization velocity v a , such that the waves are driven by the nucleator dynamics and not the actin dynamics.
By coupling the actin dynamics to a phase field, we studied the impact of the spontaneous actin dynamics on cell migration. In this context, we introduced a new expression for the nucleator current in presence of a phase field, such that nucleators are confined to the cell interior. In other phase-field studies of cell migration, conservation of particle numbers is typically not an issue and all material leaving the cell interior is simply quickly degraded [24]. If nucleators are not conserved, for example, by replacing the concentration of inactive nucleators n i by a constant, then the density of active nucleators diverges and waves are absent from the system. In Ref [21], nucleators that had leaked out of the system were reintroduced into the cell by homogenously distributing them in the cell interior. In contrast, the current −D(C r n − n r C) used in this work acts locally. All phases reported in Ref [21] are recovered and also the topologies of the phase spaces are the same in both systems with one notable exception: whereas in the present work erratic migration occurred for larger values of v a and ω d than for persistent migration, it was the opposite in Ref [21].
By analyzing the mean-squared displacement of the simulated cells, we characterized their persistent random walks in terms of a diffusion constant, a persistence time, and the cell speed. We linked these effective parameters to the actin-polymerization speed v a and the strength ω d of the negative feedback of actin on nucleator activity. It showed that these parameters had a strong effect on the effective diffusion constant and the persistence time, whereas the cell speed varied only by a factor of two. This suggests that by changing the pool of available actin monomers, cells can control important aspects of their random walks. This might allow notably cells of the immune system patrolling an organism for pathogens to adapt their behavior to the tissue they reside in.
In this work, we have neglected the effects of molecular motors, which can generate contractile stresses in actin networks. Their effects can be included into the description [28], which is in particular important when stress fibers are present [36]. The description we considered rather applies to cells that do not adhere to a substrate, like the immature dendritic cells studied in Refs [9,10]. Still, the migration of immature dendritic cells depends on the presence of molecular motors [9]. Further work is necessary to address the influence of molecular motors on actin polymerization waves. Similarly, the effects of hydrodynamic flows on these waves [37] remain to be studied.
Furthermore, it will be interesting to study in future work collective cell migration driven by spontaneous actin-polymerization waves. Previous phase-field studies revealed how steric interactions between cells can lead to collective migration [38,39] and how topographic surface structures influence this behavior [40]. In the context of our work, one might expect interesting synchronization phenomena between actin waves in different cells.

Appendix: Numerical implementation of the dynamic equations
A self-written CUDA program was used to solve the nondimensional dynamic equations efficiently on graphics processing units (GPUs). The system was discretized into 256 points on each axis in both 1D and 2D. The protein densities and the phase field were updated using the explicit midpoint rule with adaptive step size control. Fourier transformations were used to compute spatial derivatives, a method that has a higher accuracy than a finite difference scheme. The code and a more detailed description of the algorithm are available at [41].

Appendix: Wave profile
In this appendix, we determine the actin and nucleator densities for a wave traveling at velocity v.
The homogenous solution c h (x) to this equation can be written as where l ¼ k d Lv=ðv 2 À v 2 a Þ and analogously l a ¼ k d Lv a =ðv 2 À v 2 a Þ. In the above equation, the amplitude of the homogenous solution is fixed by the conditions c h (0) = c 0 and c 0 h ð0Þ ¼ v 0 . The solution to the in-homogenous Eq (37) with the source term SðxÞ is obtained by the method of variation of constants. We write c 0 = AS(x) and v 0 = AS 0 (x), where A is the Wronskian of our system and arrive at the full solution cðxÞ ¼ aLl where l ¼ v=ðv 2 À v 2 a Þ and analogously for λ a . The solution corresponds to a fraction of vÀ v a 2v of the scaled nucleator density decaying on a lengthscale of L À ¼ 1 lÀ l a and a fraction of vþv a 2v decaying with L þ ¼ 1 lþl a . The decaying part of the actin wave can be fitted perfectly with the single parameter a vÀ v a 2v e À l À x þ vþv a 2v e À l þ 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 @ t n a ¼ D a @ 2 x n a þ ð1 þ on 2 a ÞðN À n a Þ À o d cn a ð40Þ 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, with n tot being the average total nucleator density. 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).
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: where n i ðxÞ ¼ n tot À D a n a ðxÞ À ð1 À D a ÞvL e vL À 1 is the distribution of inactive nucleators. This non-linear integro-differential equation can be solved using the variational ansatz of Sect Wave solutions.
Supporting information S1 Video. Example of a traveling wave solution to Eqs (1)-(4) in two dimensions with periodic boundary conditions for v a = 0.44, ω d = 0.32. Other parameters as in Table 1