Rotating lamellipodium waves in polarizing cells

Cellular protrusion- and lamellipodium waves are widespread for both non-motile and moving cells and observed for many cell types. They are involved in the cell’s exploration of the substrate, its internal organization, as well as for the establishment of self-polarization prior to the onset of motion. Here we apply the recently developed phase field approach to model shape waves and their competition on the level of a whole cell, including all main physical effects (acto-myosin, cell membrane, adhesion formation and substrate deformation via traction) but ignoring specific biochemistry and regulation. We derive an analytic description of the emergence of a single wave deformation, which is of Burgers/Fisher-Kolmogorov type. Finally, we develop an amplitude equation approach to study multiple competing rotational waves and show how they allow the cell to transition from a non-moving state towards a polarized, steady moving state. Lamellipodial waves are a very general phenomenon observed in many cell types and is a typical phenomenon for animal cells to adhere and move along substrates. The authors present a model showing that the dynamics of these waves can be reproduced with a minimal, well-defined set of parameters.

A nimal cells adhere and move along substrates using a machinery involving the acto-myosin cytoskeleton for force generation [1][2][3] and adhesion protein complexes to link and transfer these forces to the substrate 4 . In many cases, spreading and motion are not monotonous processes in time but involve actin-and shape-wave phenomena [5][6][7] . For instance, periodic lamellipodium contraction waves have been found in adhering fibroblasts, traveling both rearwards and laterally, i.e., along the cell's periphery 8,9 . Lateral waves typically annihilate 10 , hinting at excitable medium type of waves. Similar waves were also found in stationary Xenopus tissue culture cells 11 . A biological reason for the emergence of waves may be that cells try to explore the substrate to find other cells or to select an optimal environment. More recently, protrusion waves have been suggested to be also important for structuring the position of adhesion sites 12,13 . Interestingly, shape waves have also been found in Dictyostelium discoideum, without any contact to the substrate 14 . Finally, non-steady moving cells often display actinrelated shape waves. For instance, wave generation was reported for crawling keratocyte cells 15,16 on highly adhesive substrates, where lamellipodial protrusion waves traveled at the cell's leading edge from side to side with periods of few minutes.
A second important occurrence of cellular shape waves relates to the spontaneous onset of shape polarization and directed cell migration. There, the development of a single protrusion is needed that grows and, when a sufficiently strong front-back asymmetry has been achieved, leads to persistent motion 17 . Instead of this generic scenario, the formation of multiple lamellipodia was reported 18 upon a decrease in membrane tension that apparently impedes the establishment of such a global asymmetry. Similarly, by studying keratocytes from different developmental stages of zebrafish, ref. 19 reported about cells that form several (typically 2-4) protrusions that induce an overall rotational motion of the cell, corresponding to lamellipodia running around the cell's periphery. Here the multiple protrusions are caused by regulatory processes, the rotating cell type having a higher myosin light chain kinase level, a molecule which decreases the lifetime of protrusions by increasing myosin activity. Quite similar multiple protrusion and rotating states have been also found for standard fish keratocytes directly after making first contact to the substrate when arriving from suspension 20 . There, apparently different protrusion waves are coupled and competing, with eventually one winning and taking it all and polarizing the cell.
One can argue that cell protrusion and lamellipodium waves are widespread and observed for many cell types. In contrast to actin polymerization waves that occur also in bulk 7,[21][22][23] , the coupling to both the membrane and the substrate makes the problem much richer. One-dimensional (1D) models for this phenomenon have been proposed in refs. 11,16 , including several regulation steps. In two dimensions, multiple protrusions were described by a simple phenomenological model based on the celledge dynamics controlled by the distance from the cell center 20 . However, to date no model of lamellipodia waves exists on the level of a whole cell. Using the recently developed phase field approach, see refs. [24][25][26][27][28][29][30][31][32][33][34][35][36] for recent comprehensive reviews, we demonstrate the occurrence of protrusion waves and analyze their formation and competition within a minimal physical model. Our study revealed that the onset and competition of rotating lamellipodium waves can be captured in the framework of a model incorporating cell shape dynamics, actin polymerization, substrate deformation, and adhesion. Our numerical analysis of the phase field equations is supplemented by reduced models of cell membrane dynamics where qualitatively similar wave phenomena are observed. We demonstrate that the formation of rotating protrusions can be interpreted as a propagation of shock waves described by a Burgers-like equation for the position of the cell membrane. Furthermore, simple ordinary differential equations for wave competition and the onset of cell polarization can be derived, explaining why standing waves (corresponding to breathing shape modes) are ultimately unstable and how waves trigger the onset of cell polarization and motion. This work hence demonstrates that, while biochemical regulation pathways are obviously relevant to orchestrate the phenomenon and tailor it to the desired function, the emergence of waves only needs the interplay of protrusion, adhesion, and contraction. In addition to migrating cells, our findings may provide additional insights into the onset and complex shape dynamics of active droplets 37,38 .

Results
Rotating lamellipodia in the computational model. The numerical results were obtained using an effectively twodimensional (2D) model framework, whose parts have been developed previously in refs. [27][28][29][30][31] . Namely, the cell shape is described by a phase field ρ(r, t), which evaluates to 1 inside and to 0 outside with a smooth transition region in between these two states defining the membrane position. The dynamics of the actin cytoskeleton is described via a vector field p(r, t), which has a source located at the membrane (modeling actin nucleation) and exerting forces on it (modeling actin polymerization ratcheting, parameter α) provided that adhesive bonds A(r, t) have been formed. Adhesive bonds follow a reaction-diffusion dynamics inside the cell and transfer pushing actin forces to the substrate via traction forces T(r, t). The substrate is deformed (described via the displacement field u(r, t) and Kelvin-Voigt viscoelastic response of shear modulus G), which can lead to adhesive rupture. This "feedback loop" between force exertion, transmission, and bond rupture can induce stick-slip motion, also of parts of the cell, as described previously 28,29 . Details of the model are given in the Methods section.
Although the model was originally designed to describe cell motion, Fig. 1 shows that it also captures a variety of rotational lamellipodium states. While the model is known to describe polarized, moving states-see Supplementary Movie 1-depending on parameters a cell may not be capable of polarizing, or the initial condition may not be sufficiently biased to allow for it. In both cases, transient rotational states emerge, with the cell finally settling to a non-motile circularly symmetric state in the former case or polarizing with a single lamellipodium and adopting a persistent motion in the latter. The simplest case is one rotating shape wave as demonstrated in  Figure 2 displays the traction distribution of a polarized cell (Fig. 2a) compared to a cell displaying two counter-rotating lamellipodia (Fig. 2b) (similar to Fig. 1e-h). It is evident that both, the traction distribution and the polarization (as already shown in Fig. 1) have local maxima close to the shape deformations. Hence, they can be interpreted as local lamellipodia, similar to those occurring in the experiments [18][19][20] .
Since the cell shape, the actin distribution, as well as the adhesion build-up and substrate deformation are all coupled, lamellipodium waves are expected to be in nonlinear competition. In fact, typically all but one wave is damped in the long term, with a dominant wave gaining in amplitude and polarizing the whole cell. For cells displaying a single shape wave, the wave travels along the cell's periphery for several turns until it eventually polarizes the whole cell. One can conclude that the wave phenomena are transient and that they finally induce directed motion. This suggests a coupling of deformation modes and the translational mode, see the discussion of wave competition below. Figure 2c shows a phase diagram of the cellular dynamics obtained when varying the substrate stiffness G vs. the actinassociated pushing rate, α. For soft substrates, the cell polarizes immediately without noticeable initial oscillation. In turn, if the substrate is too stiff, all oscillatory modes are strongly dampened and the cell remains stationary. For slow actin dynamics, cells display a very slow, ameboid-like gliding behavior, as its effect is too weak to polarize the cell, see Supplementary Movie 7. For too high values of α, the cell is unable to maintain its integrity, with (lamellipodium) fragments beginning to split off from the main cell body. In between these boundaries in parameter space, rotational states with a diversity of modes exist: single rotating lamellipodia, double co-rotating and counter-rotating lamellipodia, and even higher modes. This corroborates that the wave phenomena emerge in the transition region between conditions where cells are easily polarizable and conditions where cells can not polarize, as an additional route to break the symmetry and initiate cell motion. Overall, the number of protrusions tends to increase with the cell's size, although the detailed dependence is not simple.
Fourier analysis of rotating states. To obtain additional insight, we extracted the Fourier modes of the cellular shape as a function of time. The position of the membrane was identified with the location of the ρ = 1/2-isoline. The membrane was then discretized using quadratic interpolation to produce discrete Cartesian coordinates, which were transformed into radial coordinates (ϕ, r(ϕ)) relative to the center of mass (c.o.m.) determined by r c:o:m: = R ρr dxdy= R ρ dxdy. The Fourier modes as defined bỹ were then determined numerically performing a Fast-Fourier Transform (FFT) of r(ϕ). The leading order Fourier modes are shown in Fig. 3a for the case of a single rotating wave and in Fig. 3b for the case of two counter-rotating waves. For the single wave case, the leading Fourier modes monotonically increase until the rotating wave is strong enough to polarize the cell. For the two wave case, the plot shows a periodic pattern-as the two waves separate on one side of the cell and collide on the opposite side. Eventually, one wave wins the competition, the weaker becomes damped and from this point on the Fourier modes resemble the single wave case. In the process of polarization, the cell typically forms a transient elongated elliptical shape, as demonstrated in Fig. 1c. Finally, it eventually stabilizes into the steady gliding shape, as seen in Fig. 1d. The build-up of this transient elongated shape is why the second mode dominates the Fourier spectrum prior to cell polarization.
In the following, we will investigate in more detail the formation of the waves, how these waves interact and compete with each other, and how the rotational states eventually polarize the cell. Formation of rotational waves. To understand the formation process of localized lamellipodium waves, we performed an asymptotic reduction of the full computational model. We approximate the phase field for a quasi-circular cell close to the stationary state as where r 0 (ϕ, t) is the location of the cell membrane, see also ref. 39 . As shown in Methods, we can approximate the stationary values of the adhesion field and the actin polarization close to the cell's boundary by A s ≈ a nl ρ/s, and Using these approximations, the dynamic phase field equation can then be reduced by evaluating the corresponding solvability condition to a single 1D partial differential equation for the membrane deformation, δr 0 , from a reference value R 0 , i.e., δr 0 (ϕ, t) = r 0 (ϕ, t) − R 0 : where δV = 1 2 R 2π 0 δr 0 þ R 0 ð Þ 2 dϕ À πR 2 0 is related to the volume conservation of the cell (corresponding in the effective 2D model to the conservation of the cell's contact area).α is an effective parameter containing the actin polymerization rate, the pushing force it exerts, and the myosin-mediated contractile forces in a combined fashion (see Methods, Eq. (27), and subsequent discussion). The equation has a diffusive term and a contribution from volume conservation. The third term leads to linear growth of perturbations, while the nonlinear term~(∂ ϕ δr 0 ) 2 saturates the growth. Finally, the pushing due to actin polymerization leads to a source term attempting to globally increase δr 0 .
This Burgers-like equation suggests that the lamellipodium waves are a kind of shock wave traveling around the cell's periphery as the result of some small local perturbation in δr 0 . In fact, a numerical integration of Eq. (3), using as initial condition a small narrow peak, is shown in Fig. 4a and Supplementary Movie 8: two rotating waves develop and travel in opposite directions; the front of the wave is where j∂ ϕ δr 0 j is the largest. The diffusion term diminishes the back end of the wave, while δV ensures the conservation of the overall volume of the cell. The qualitative behavior is the same as for the full phase field model shown in Fig. 4b, where the cell radius r 0 was interpolated from a simulation, see the Fourier analysis of rotating states discussed above. The deformation dynamics of lamellipodia is caused by the simultaneous action of molecular motors and the ratcheting of actin, pushing the membrane out radially in a, in general, nonuniform fashion. This is countered by a restoring force resulting from the conservation of area. Since actin is oriented roughly normal to the membrane, once a lamellipodium protrusion forms the actin along the sides of the lamellipodium has a significant component in the azimuthal direction. This effect is present in Eq. (3) as the last term as j∂ ϕ δr 0 ð Þj is greatest along the lamellipodium sides. The azimuthal component then forces the lamellipodium to begin rotating with the lamellipodium sides becoming the wave fronts: the actin along one side forces it clockwise and the other counterclockwise, effectively splitting the protrusion into two separate counter-rotating waves. The most generic scenario observed is hence a pair of counter-rotating waves originating from a single lamellipodium.
By applying a Cole-Hopf transformation, δr 0 = ðD ρ =αÞlogðWÞ, the Burgers-like Eq. (3) can be reduced to a Fisher-Kolmogorov-Petrovsky-Piskunov-type equation 40 μδV ÀαÞ. From its structure, this equation is known to allow for moving fronts connecting stable and unstable equilibria in certain parameter regimes. Further analysis of this equation is detailed in Methods. Especially, the existence of transient traveling wave solutions around r 0 ≃ R 0 is found by phase plane analysis, supporting the numerical findings.
Competition of rotational waves. While the previous section explains the formation of the rotating waves, the question remains how different rotational modes are coupled together, as well as how the waves affect the cell's velocity and vice versa. We here propose a phenomenological dimensionless reduced system of ordinary differential equations to explain the dynamics of wave competition, as well as how waves trigger cell polarization.
For simplicity, consider two rotating waves with complex amplitudes R 1 and R 2 , each with corresponding frequency ω 1 and ω 2 . The waves are assumed to be small perturbations of the circular cell, i.e., r 0 (ϕ, t) Then, as detailed in Methods, a system of amplitude equations can be derived by Here ϵ is the supercriticality parameter, ζ describes amplitude saturation, and γ is the wave coupling parameter. For γ > 1, only one wave survives competition, the surviving wave being determined by the initial conditions. As we obtain γ = 2 (see Methods), this is indeed the case for the system at hand, explaining the behavior observed in Fig. 3b). As shown in the Methods section, the supercriticality parameter ϵ is proportional to the (square of the) actin polymerization parameter β of the phase field model: in fact, actin pushes the membrane outwards radially and causes the waves to grow in size. The amplitude saturation parameter ζ in turn is proportional to the parameter d 0 in the underlying model, which describes the rupturing of adhesions due to high cell traction: the higher the rupturing rate, the lower the force that can be transferred and consequently the lower the protrusion force and the size of the protrusion. Note that this parameter was not significant for the formation of the waves, but it is for the competition of multiple waves. Now let us consider the time evolution for the cell's center of mass velocity. It is convenient to formally rewrite the velocity vector, v ¼ v xx þ v yŷ , as a "complex velocity,"V ¼ v x þ iv y . First, we take the velocity to be uncoupled to the rotating waves. Motivated by the observation that in the full model cells exhibit a subcritical transition to the moving state 27,28 , we suggest the following generic equation for the onset In fact, this equation has a subcritical transition toward motion, if a, b, c > 0 and (b/c) 2 − 4a/c > 0. In this case, there are three equilibria: V 0 = 0 and Both V 0 (stationary cell) and V + (polarized cell) are stable, while V − is an unstable solution. Now, the symmetry-i.e., the fact that the translation mode has the angular dependence exp(iϕ)-suggests the following leading order coupling of the center of mass velocity and the rotational waves For our sakes, we are only concerned with the real components of the rotating waves. Results from Eqs. (9)-(11) can also be compared to simulations using the full computational model in a similar fashion as described above for the Fourier mode analysis: we interpolated r 0 and determined the location of the peak of the jth wave in polar coordinates (r j , ϕ j ). Its amplitude can then be evaluated as R j = exp(iϕ j )(r j − R 0 )/R 0 . Finally, the cell's center of mass velocity can be evaluated by taking the time derivative of the center of mass position, v ¼ _ r c:o:m: . Let us first consider the case of one wave rotating with frequency ω 1 = ω. As shown in Fig. 5a as the solid curves, the rotating wave perturbs the cell, inducing oscillations of the cell's center of mass velocity, as seen in the terms ∝q 1 , q 2 in Eq. (9). Eventually, the velocity stabilizes to V + since, in turn, it dampens the rotating wave as expressed by the terms ∝w in Eqs. (10) and (11). In the full computational model, the dashed curves, the cell has a transient elongated elliptical shape when it polarizes, accompanied by an increasing velocity, and the cell eventually settles in the steady moving state. We can also consider the case of two counter-rotating waves such that ω 1 = −ω 2 = ω, as shown in Fig. 5b. Qualitatively, the counter-rotating case is similar to the single wave case: as before, the rotating waves disturb the cell causing oscillations in the velocity. The two waves compete until only one survives and the other is damped, with the initial conditions determining which one of the waves survives. Eventually, the cell polarizes and its velocity dampens the last remaining rotating wave (blue curve). Figure 5a, b directly compare the dynamics of the reduced amplitude equations to the . Although the oscillations differ in amplitude and are more complex in the full model-due to the neglect of higher modes in the reduced model-the general behavior is well captured. Finally, we can consider the case of a single rotating wave with ω 1 = ω and using a smaller ϵ value. This case provides insight into the ameboid-like motility as found in Fig. 2b. Note that this motility state is common for low values of α (one has ϵ / α 2 as we used β ∝ α, see Methods and Supplementary Table 1). Often the oscillations for this state are not distinct rotating waves about a roughly circular cell but instead cause shape deformations of the whole cell. It is then instructive to measure the strength of the oscillation by tracing the second moment of ρ with respect to the center of mass: R r À r c:o:m: ð Þ 2 ρðx; yÞdxdy. In this case, the rotational oscillation is not capable to drive the cell's velocity beyond the unstable node (V − ). Instead, eventually both the center of mass velocity and the shape wave coexist and oscillate out of phase with respect to each other, as seen in Fig. 5c. Therefore, a possible explanation for the ameboid-like motion is that weak rotational oscillations continue to exist, forcing the cell to periodically change speed and direction, similar to the meandering core of a spiral wave 41 .

Discussion
We used a whole-cell model that incorporates all relevant physical processes-the dynamics of the actin-myosin network, the formation of adhesive bonds to the substrate, and the substrate deformation via cellular traction forces-to study shape waves of spreading and polarizing cells. Although we neglected the specific biochemistry and regulation pathways present within the cell, rotating lamellipodium waves similar to those seen experimentally have been found and we could investigate in depth their occurrence, their competition, and their triggering of directed motion.
Our study suggests that the rotating states occur only within a limited range of the parameters-as exemplified by the actin propulsion rate and the substrate stiffness-and that the emergence of these waves enlarges the parameter range wherein cells are able to polarize. It also shows that the emergence of these waves only requires protrusion, adhesion, and feedback via traction-mediated breaking of adhesive bonds. One can anticipate a roughly linear increase in the number of protrusions with an increase in cell size. However, owing to nonlinear interaction between the protrusion waves, the cell size dependence is more complex and requires an in-depth computational study, which is beyond the scope of this work.
Our analysis is presently limited to a 2D description. This approximation is justified by the fact that the lamellipodium thickness is much smaller than the spatial scale of the protrusions and that the cell body is typically not participating in the front dynamics. One may hence expect that a three-dimensional consideration of the phenomenon will not change the qualitative conclusions on the physical mechanisms associated with rotating lamellipodia waves. A delicate point, however, is to what extent the spreaded area is conserved or changed by material flow from the third dimension. Hence, a 3D analysis may bring additional insights into the spatial and temporal organization of lamellipodium waves but at this point is still computationally prohibitive and left for future study.
Another important question is to what extent external stimuli and/or internal regulation processes lead to similar shape waves and rotating states or whether they only orchestrate the ones induced by the described physical mechanism. To answer this question, as demonstrated in ref. 42 , the approach can be generalized to account for the currently discussed regulation pathways, like Rho/Rac guanosine triphosphate (GTP)-ases 43 or Ena/VASP proteins 16 .
The developed analytical reduced description for how perturbations in the cell's periphery lead to the formation of rotating waves will stimulate further analysis. Not only is its agreement with the full numerics very reasonable but, importantly, such equations of Burgers-or Fisher-Kolmogorov-Petrovsky-Piskunov-type, respectively, put the cellular shape waves in the mathematically wellstudied framework 44 of waves in reaction-diffusion systems or excitable media and allow future analytical in-depth studies.
Finally, the proposed amplitude equations allowed us to rationalize competition of cellular shape waves and on a phenomenological level also how these waves can trigger-via generic symmetry-imposed couplings-the onset of the cell's center of mass velocity. This framework is again capable of reproducing the dynamics seen in the full cell model semi-quantitatively. It may also provide insight into the onset of rotating states in related systems like active droplets 45 and the so-called "circus movements" 46,47 found in blebbing blastomere cells from, e.g., fish embryos, where single blebs travel around the cell's periphery for long times.
While the present work focused on the dendritic actin near the membrane, future work could also incorporate actin waves experimentally seen in the bulk of the cell 21,23 and study how they induce and/or affect rotational states. Furthermore, similar concepts and ideas can be used for the description of active poroelastic droplets modeling Physarum polycephalum dynamics 48,49 . Another intriguing aspect is the possibility of synchronization of shape waves of different, close-by cells due to mechanical and chemical interactions, calling for dedicated experiments.

Methods
Model. As developed in refs. 27,29,31 , we use the following fields to describe effectively the 2D cell: first, an Allen-Cahn phase field, with ρ(r, t) = 1 inside the cell and ρ = 0 outside, augmented by a (2D) volume conservation to model the cell shape; second, a vector field p(r, t) describing the actin orientation and degree of ordering. And third the density of engaged adhesive bonds A(r, t) that is coupled to the deformable substrate.
The phase field equation reads with the phase field parameter given by Here δ ¼ 1 2 is the stationary point and the second term introduces 2D contact area conservation with target area V 0 ¼ πR 2 0 for some fixed radius R 0 , with the stiffness of the constraint given by μ. As explained in more detail before 27 , the term ∝σ models actin network contraction by myosin molecular motors with constant rate σ and the last term models the motion ("advection") of the membrane due to the ratcheting of actin filament polymerization, with velocity α. Only adhering filaments can transfer the force to the substrate, hence the proportionality to A.
The equation modeling actin dynamics reads The first term accounts for diffusion of actin filaments/elasticity in the ordered state. The second term is a source and models that actin is generated at the membrane with (polymerization) rate β. It includes the effect of membrane tension feedback on polymerization as developed in ref. 31 : The constant T 1 sets the strength of the membrane tension, c ¼ À∇ Á ∇ρ ∇ρ j j is the local curvature of the cell membrane and δS ¼ SðρÞÀS 0 S 0 is (proportional to) the excess surface area, where SðρÞ = R ∇ρ j jdxdy and S 0 = 2πR 0 þ 4π ffiffiffiffiffiffiffiffi 2D ρ p log 2 are, respectively, the instantaneous and the reference surface area, the latter being defined as S(ρ s ) with ρ s from Eq. (19) below. The next two terms in Eq. (14) are sinks and model the degradation/depolymerization of actin filaments inside the cell with the typical timescale τ 1 and suppress actin outside the cell, respectively. The final term coarsely models the myosin motors' ability to form antiparallel actin bundles reducing the actin polarization and breaks the symmetry between the front and rear of the cell 27 .
The engaged adhesion bonds are modeled by an equation of reaction-diffusion type as The first term is diffusion and the second term models attachment, a linear contribution in the presence of actin, and a nonlinear one modeling collective effects. The third term restricts the amount of adhesions by excluded volume, ∝s, and substrate-dependent detachment. The detachment rate accounts for the fact that the substrate is soft and deformed by the cell via the traction forces it exerts on it. In Eq. (16), we made the assumption that adhesions will break where the local substrate deformation directly underneath the cell, |u z=H | = |u(x, y, z = H)|, exceeds a critical value u c . The substrate is modeled as a 3D isotropic homogeneous visco-elastic solid of Kelvin-Voigt type, introducing the shear modulus G, a viscosity η (modeling viscous losses due to adhesion rupture), and the effective height of the substrate, H. This treatment has been thoroughly discussed in refs. 29,36 , yielding the displacement at the top of the substrate, u z=H , due to the traction force distribution exerted by the cell, defined as Here the first two terms in the bracket are the counter forces related to actin polymerization and membrane tension preserving the membrane surface area. The final term is of frictional origin and ensures that the net force on the substrate is zero, as it should, the cell being an isolated self-propelled object.
Numerical method. The system of equations Eqs. (12), (14), and (15), including the displacement field u z=H , were solved numerically using a highly parallelized algorithm implemented on a graphics processing unit, using Compute Unified Device Architecture (CUDA). As the phase field takes care of the deformable and/ or moving boundaries, the problem can be solved on a double periodic square domain using classical operator splitting 50 , pseudo-spectral method, and FFT. For a more comprehensive review of the numerical scheme, see ref. 36 . The phase field curvature was computed (using finite differences) only in a tube around the cell membrane 31 . Typical parameter values are listed in Supplementary Table 1. Field distributions for a stationary cell are shown Supplementary Fig. 1.
The qualitative motility behavior was tested to be consistent when refining the spatial mesh or the time step. Simulations were started with a small amount of random noise added within the circular region modeling the cell, where the initial value of the phase field is ρ(x, y, t = 0) = 1. The polarization of the cell and/or the formation of rotating waves was consistent irregardless of the initial noise; the noise realization, however, affects the direction of polarization or, in the case of multiple waves, which wave dominates and survives.
Derivation of reduced descriptions of wave formation. Consider first the homogenous 1D case of the phase field equation. Rewriting it as ∂ t ρ ¼ L ρ ½ρ using the self-adjoint differential operator L ρ it is easy to show that is the solution of L ρ ρ s Â Ã ¼ 0, i.e., ρ s solves the stationary 1D phase field equation [for ease of notation, one definesr = r À R 0 ð Þ =ð2 ffiffiffiffiffiffiffiffi 2D ρ p Þ]. If we assume a large cell (R 0 = ffiffiffiffiffi ffi D ρ p ) 1), it can be also shown that Eq. (19) is the leading order approximation for the stationary solution of the 2D homogeneous phase field equation. Now consider the full system of equations. For simplicity, we will ignore the diffusion terms of p and A, take δS = 0, and define δV ¼ R ρdxdy À πR 2 0 : We assume that the terms involving d(u) and a 0 in the adhesion equation are negligible for the stationary case; we will include and discuss them later. Therefore, we can approximate the stationary solution as A s % a nl s ρ. For the region inside the cell (ρ = 1), the stationary values then are p s = 0 and A s ≈ a nl /s, while outside the cell (ρ = 0) one trivially has p s = 0 and A s = 0.
The actual region of interest is the one near the cell membrane (r j j ( 1). Expanding Eq. (19), one gets the stationary values in this region as ρ s ≈ 1/2 and hence A s % a nl 2s . Furthermore, we can approximate p s by setting ∂ t p s = 0 and taking the scalar product of Eq. (21) with ∇ρ to obtain which can be solved to get approximately definingβ as the effective polymerization vs. degradation rate of the actin filaments. These approximations for the stationary values of ρ, A, and p are compared to the numerically obtained values in the full model in Supplementary Fig. 2. We now generalize Eq. (19) to non-circular cells with radius r 0 (ϕ, t), as already given in Eq. (2). Assuming small perturbations of the stationary state, we can use A ≈ A s and p ≈ p s and write Eq. (20) as L ρ ½ρ ¼ 0 with Now, in order for a nontrivial solution of Eq. (25) to exist, its solvability condition h∂ r ρ; L ρ ½ρi ¼ 0 must be satisfied, where the inner product is defined by A Brdr dϕ and ∂ r ρ is the solution to the linearized Eq. (20). In the large cell limit and using Eq. (2), a straightforward calculation shows that the solvability condition is satisfied given the following equation is fulfilled with the effective actin propulsion parameter Note that this parameter contains all main non-equilibrium effects of the cytoskeleton: polymerization vs. degradation (β), actin ratcheting, and force generation against the membrane (α) provided actin adheres (a nl ), as well as actomyosin contraction (σ). Finally, introducing δr 0 (ϕ) = r 0 (ϕ) − R 0 as the radial perturbation from R 0 , Eq. (26) can be rewritten in the form given by Eq. (3). This Burgers-like equation can be transformed into the Kolmogorov-type Eq. (4) as indicated in the main text.
Derivation of equations describing wave competition. Having demonstrated that rotating waves can form from a perturbed circular cell, we now assume their existence and deduce equations for how the amplitudes of several waves as well as the cell's center of mass velocity couple. In order to do so, we include the terms of the adhesion Eq. (15) that were previously neglected (but still neglect diffusion). Similar to the previous section, we define and accordingly write Eq. (22) as Taking the previously obtained result, A s ¼ a nl s ρ, we know that it fulfills L A ½A s ¼ 0. Now we consider a cell with two harmonic rotating waves, with small nondimensional amplitude R i and angular frequency ω i . We can hence write the cell's radius r 0 as r 0 ðϕ; tÞ ¼ R 0 1 þ R 1 ðϕ; tÞe Àiω 1 t þ R 2 ðϕ; tÞe Àiω 2 t Â Ã þ c:c:; with R 1 j j; R 2 j j ( 1. Furthermore, we will assume Aðr; ϕ; tÞ % a nl s ρðr; ϕ; tÞ and pðr; ϕ; tÞ % Àβ∇ρðr; ϕ; tÞ, with ρ given by Eq. (2).
The step-like detachment of adhesive bonds, d(u), as defined in Eq. (16) and entering Eq. (30) is triggered by an increase of substrate deformation, |u z=H | 2 , caused in turn by the traction force exerted by the cell. As evident from Fig. 1a, the traction distribution is significant primarily at the front of the traveling waves. This can also be seen by explicitly evaluating |T| 2 ≈ |−ρAp| 2 using the approximations listed above, showing that the traction is high, due the perturbations of the cell's radius r 0 , only near the front where actin polarization is high. We can therefore approximate Eq. (16) by d(u) ≈ d 0 |∂ ϕ r 0 | 2 |∇ρ|, which, as we assume harmonic traveling waves and extract only slowly-varying contribution, further simplifies to dðuÞ % d 0 R 1 e Àiω 1 t þ R 2 e Àiω 2 t j j 2 ∇ρ j j where R Ã i denotes the complex conjugate of R i . Using all the approximations motivated above, Eq. (30) reads and integration yields the condition It is then straightforward to show, using Eqs. (2), (31), and (33), that Eq. (35) is satisfied if the amplitudes R 1 and R 2 follow the coupled amplitude equations given by Eqs. (5) and (6)  Code availability. Code used to generate the data shown here are available upon reasonable request from the corresponding author.

Data availability
Data that supports the findings of this study are available upon reasonable request from the corresponding author.