Membrane waves driven by forces from actin filaments

Membrane waves propagating along the cell circumference in a top down view have been observed with several eukaryotic cells (Döbereiner et al 2006 Phys. Rev. Lett. 97 10; Machacek and Danuser 2006 Biophys. J. 90 1439–52). We present a mathematical model reproducing these traveling membrane undulations during lamellipodial motility of cells on flat substrates. The model describes the interplay of pushing forces exerted by actin polymerization on the membrane, pulling forces of attached actin filaments on the cell edge, contractile forces powered by molecular motors across the actin gel and resisting membrane tension. The actin filament network in the bulk of lamellipodia obeys gel flow equations. We investigated in particular the dependence of wave properties on gel parameters and found that inhibition of myosin motors abolishes waves in some cells but not in others in agreement with experimental observations. The model provides a unifying mechanism explaining the dynamics of actin-based motility in a variety of systems.


Introduction
Experiments have shown traveling membrane protrusions or retractions in many cell types such as mouse embryonic fibroblasts, fly wing disc cells, mouse T cells and Dictyostelium discoideum [1,3,[6][7][8]. These waves can be identified as the dynamics of the cell contour in a top down view of the cell on the substrate. The contour dynamics consists sometimes of homogeneous oscillations rather than propagating structures [8]. Membrane waves are believed to play a central role in cell motility, endocytosis and probing the extracellular matrix [9]. In the damped liquid environment of the cell, these propagating waves are maintained by active forces from the cytoskeleton. Actin polymerization produces most of the driving force for membrane protrusion [10]. Polymerization also creates a network of actin filaments providing support for the forces acting on the membrane. The network of actin filaments is cross-linked by various proteins. Among them, myosin II molecular motors not only cross-link but can also contract the network. Some of the membrane waves disappear upon inhibition of contraction of the actin gel [1,3,4,7,11], and others exist without myosin activity [1][2][3][4][5][6]8]. This suggests that the membrane waves are driven by the protrusive force of actin polymerization and the contractile stress produced by myosin II molecular motors contributes in some cases only.
Several theoretical approaches have been introduced to explain the mechanism responsible for membrane wave propagation. In the approach of [12][13][14], membrane waves are driven by curvature-mediated activation of actin polymerization. A non-decaying membrane wave was produced by the coupling between membrane shape and protrusive-contractive forces of actin-myosin and the addition of only convex actin activators in this model [13]. The intrinsic curvature values that can be induced or are preferred by membrane proteins suggest that the mechanism applies to ruffle-like waves [14]. The observation of a transition between homogeneous oscillations with vanishing curvature and waves of the cell contour induced by signaling mechanisms suggests that different dynamic regimes of one system defined by different parameter values produce both the phenomena. The homogeneous oscillations are unlikely to be produced by a mechanism relying on curvature.
Here we propose a mechanism for membrane wave generation that does not require any spontaneous curvature of actin polymerization activators. The dynamics of the cell boundary, and therefore cell shape, are determined by the interaction of actin filaments with cell membrane, cross-linking of actin filaments, polymerization-and myosin-driven retrograde flow of the actin gel and load due to membrane tension. In our model, while transiently membranebound filaments are pulling back the membrane, the polymerizing actin network pushes on the cell membrane from within (see figure 1). Membrane tension equilibrates spatially very fast, The actin network consists of two parts: a highly cross-linked actin meshwork (actin gel) and free actin filaments (semiflexible region (SR)). The membrane and the gel boundary are characterized in Cartesian coordinates by (x(s), y(s)) and (x g (s), y g (s)), respectively, with s being the contour length of the membrane or actin gel boundary. In the SR, filaments undergo cycles of attachment (k a ) and detachment (k d ): while detached filaments with contour length l d polymerize and push against the membrane ( f d ), attached filaments with contour length l a can either pull or push the membrane ( f a ). Both pushing and pulling forces depend on the contour length of the filaments as well as the orientation of the filaments (θ ) and the normal distance to the membrane (z). R ≈ l a [1 − l a /4 p ] cos θ is the projected equilibrium end-to-end distance of the polymer onto the membrane normal. but temporarily varies with the cell's projected area. This introduces a global spatiotemporal coupling in the dynamics of distant regions on the cell boundary.

The model
The model is schematically represented in figure 1. The cell is characterized in Cartesian coordinates by x(s) and y(s) with s measuring the cell's contour length. The membrane at 4 each small interval (s, s + ds) is pushed or pulled by actin filaments which are anchored in a highly cross-linked actin gel with a linear density n l . The gel forms due to cross-linking of actin filaments. The newly polymerized filaments are not yet cross-linked and therefore behave like an ensemble of individual filaments. We call the filamentous range a semiflexible region (SR) since its properties are dominated by the semiflexible filament behavior. It takes some time till that critical degree of cross-linking is reached, which turns the network into a gel. Newly polymerized filament parts move backwards during that time owing to retrograde flow. That turns the temporal transient into a spatial gradient of the degree of cross-linking. We define the location of the critical concentration of bound cross-linkers as the gel boundary. The gel boundary advances with a rate set by cross-linking and effectively shortens the free length of filaments in the SR. Since polymerization occurs at the leading edge membrane only, we obtain the lamellipodium structure shown in figure 1 with an SR at the leading edge and a gel further in the bulk. The complete protrusion is often described as consisting of the posterior lamellum with highly cross-linked and bundled actin filaments and the anterior lamellipodium, which is much less bundled. We do not identify the gel with the lamellum and the SR with the lamellipodium. However, we assume the gel boundary to be in the lamellipodium.
Cross-linkers dissociate from filaments deep in the bulk of the gel and diffuse back to the lamellipodium leading edge, where they can bind again. Solving the equations for this reaction-diffusion process explains the dependence of the gel boundary velocity on system parameters and variables [24]. The velocity increases with the free length of actin filaments with a characteristic length scalel and saturates at the maximum value v max g , which depends on the concentration of available cross-linkers [24]: The SR at each interval (s, s + ds) is formed by populations of attached and detached filaments with linear densities n a (s) and n d (s) and average contour lengths l a and l d , respectively. Filaments bind transiently via linker proteins to the membrane. The filament's attachment rate to the linker protein on the membrane is k a . The filament-linker complex exerts force f a normal to the membrane during attachment. The linker proteins are modeled as springs with spring constant k l and zero equilibrium length. The force f a depends on k l , the filament's contour length l a as well as on the distance between gel and membrane z. The filament-linker complex has a nonlinear force-extension relation which we approximate by a piece-wise linear function [20]. Let R ≈ l a [1 − l a /4 p ] cos θ be the projected equilibrium end-to-end distance of the polymer onto the membrane normal. The elastic response of filaments experiencing small compressional forces (z R ) is approximated by a spring constant k = 12k B T 2 p /l 4 a [21]. For small pulling forces (z R ), the linker-filament complex acts as a spring with an effective constant k eff = k l k /(k l + k ). In the strong force regime, the force-extension relation of the filament is highly nonlinear and diverges close to full stretching [22]. Therefore, only the linker will stretch out. The complete force-extension relation is captured by −k l (z − l a cos θ ) − k eff (l a cos θ − R ), z l a cos θ. ( Attached filaments are usually under tension and detach with a force-dependent rate k d [23] Here, k 0 d is the spontaneous detachment rate and δ ∼ 2.7 nm is one-half of the actin monomer size.
The pushing force due to actin polymerization near the edge is the main type of active force at the membrane. Detached filaments polymerize at subsecond timescales and exert normal protrusive force f d on the membrane. The membrane hinders the polymerization. The higher the filament pushing against the membrane is, the lower the polymerization velocity is. According to [43], v p decreases exponentially with f d : with v max p = k on δG the saturation polymerization velocity in the absence of any obstacle. k on is the monomer assembly rate and G is the actin monomer concentration. The pushing force f d of semiflexible polymers in the presence of an obstacle has been extensively studied in [19]. For a stiff polymer such as actin with persistence length p = 15 µm and contour length l d p , the pushing force can be approximated by f d = f cfd (ζ ), in which f c = k B T p /l 2 d is the Euler buckling force and the scaling parameter ζ is given by ζ = p (l d − z)/l 2 d . For a small compression (ζ 0.2), the dimensionless forcef reads We now write the set of equations for the dynamics of the coupled system consisting of attached and detached filaments in the SR, the boundary of the cross-linked actin gel and the membrane [16]. The average lengths of attached and detached filaments in the interval (s, s + ds) are denoted by l a and l d , respectively. The filament lengths shrink with velocity v g and grow only in the detached state by polymerization with velocity v p . The time evolution of the average number of attached filaments n a is described by a constant attachment rate k a and a stress-dependent detachment rate k d . The gel boundary characterized by (x g , y g ) advances with velocity v g and moves backward due to retrograde flow v r . Furthermore, the membrane resists motion with a drag force (coefficient η), and resists bending deformations defined by the local membrane curvature κ(s) due to membrane tension S. All these dynamic processes are captured by the following set of equations: Here primed quantities indicate derivative with respect to contour length s and the prefactors in equations (6d)-(6g) are geometrical factors explained in the supplementary data (available from stacks.iop.org/NJP/14/115002/mmedia). We assume that capping of existing filaments and nucleation of new filaments balance, and the total number of attached and detached filaments n l (s) at the interval (s, s + ds) is conserved. Consequently, the average number of detached filaments is given by defines the local velocity of gelation and points in the direction normal to the gel boundary. Membrane tension imposes an opposing force on growing actin filaments at the cell's leading edge. Tension gradients in the membrane at any point equilibrate within milliseconds in a fluid membrane such as lipid membranes. Hence, on the scale of seconds and minutes relevant to cell motility, tension is constant along the cell boundary but can change with time [27,28]. We hypothesize that membrane tension changes linearly with the cell's projected area: in which S 0 and S 1 are constants and A 0 is the cell's preferred area and deviations from this area produce an effective extra membrane tension. This implies a spatial long-range coupling of distant points on the membrane and may lead to different morphodynamic patterns on the membrane. We choose A 0 = 60 µm 2 , which gives a cell radius of the order of 4.4 µm. This is the typical size of a D. discoideum cell. In addition, κ(s) is the signed local membrane curvature defined as where and are the first and the second derivatives with respect to s. In most motile cells, the actin network is simultaneously transported away from the leading edge in a process known as retrograde flow [34]. The flow is driven by the net force exerted by the SR on the leading edge membrane and contractile forces of myosin motors (µ). Friction of the flowing actin gel against the intracellular interface of cell adhesion complexes (ξ ) and gel viscosity η g damp retrograde flow. A semi-analytical solution of the corresponding gel equations in the radial direction provides the expression for retrograde flow v r as [24] v r = a 0 7 with a 0 ret and b 0 ret written as Here, L is the gel depth and h 0 is the lamellipodium height. The term a 0 ret is the contribution of gel contraction to retrograde flow. It increases with motor activity. Friction between the actin gel and the adhesion complexes or stress fibers and gel viscosity slow it down. Myosin II is found at lower levels throughout the lamellipodium than in the rear of polarized motile cells. The high concentrations in the cell rear, equation (10a), cause higher flow at the cell near, which is consistent with experiments.
We use this expression of retrograde flow calculated in a one-dimensional geometry as an approximation for the radial flow in our model cell in order to simplify the model. Exact calculation of retrograde flow would, of course, require the simulation of the two-dimensional gel equations [29,31].

Lateral traveling waves
We integrate numerically the set of differential equations (6) and combine time behavior with a stability analysis. The cell has a preferred radius R close to but larger than √ A 0 /π . The cell shape is a circle initially. The distribution of actin filaments (n 0 l ) and saturation polymerization speed (v max p ) over the cell boundary is assumed to be uniform. Linear stability analysis illustrates that within certain ranges of polymerization, attachment, detachment and cross-linking, the starting fixed point is unstable. We choose that unstable state as the initial condition for the simulation. The cell contour reaches its asymptotic behavior after some transient (see, e.g., figure 3). The parameters µ and ξ are among important parameters governing the dynamic modes of the model described by equation (6) as can be seen in figure 2 (see [15,16,24] for bifurcation diagrams in dependence on SR parameters). The parameter µ corresponds to myosin activity and the friction coefficient ξ is proportional to the cell-substrate adhesion site density.
In the regime with only unstable fixed points and a stable limit cycle, the membrane position oscillates, spatially desynchronizes and the cell exhibits lateral wave patterns propagating on the circumference. Since distribution of actin filaments is homogeneous on the cell edge, there is no net cell locomotion. The cell is spread on the substrate, has a roughly circular shape and the membrane periodically undergoes global shape changes.
Examples of the membrane normal velocity map and curvature map are shown in figures 3(a) and (b). The wave pattern is also visible in the other dynamic variables, such as polymerization velocity v p , radius of cell and gel boundaries, net normal velocity of the gel boundary, average length of attached and detached filaments (l a , l d ) and total pulling/pushing forces of bound/unbound filaments (n a f a , n d f d ) (see figures 3(c)-(j)). The cell and gel boundaries are measured with respect to the fixed origin at (0, 0) (the laboratory frame of reference). Interestingly, the wave propagation is also observed in the net velocity of the gel boundary v g − v r as well as in the radius of the gel boundary (figures 3(c) and (j)). In this The leading edge membrane velocity and gel boundary velocity are in antiphase during the oscillations shown in figure 4. When the membrane is pushed forward, the opposite force acts on the gel boundary, drives retrograde flow and tips the balance between flow and crosslinking to net inward. When the membrane is pulled inward, the opposite force slows down retrograde flow and the gel boundary moves outward. Oscillations and the alternating forces arise from an interplay of free filament length dynamics and attachment and detachment of filaments. The free filament length shortens during membrane retraction, since the gel boundary advances and the membrane retracts. Filaments become stiffer due to this shortening, which increases the force exerted by polymerizing filaments on the leading edge membrane. Attached filaments can stand that force for a while, but when it reaches a critical strength, they are ripped off the membrane. With little left holding it back, the membrane jerks forward and the phase of membrane protrusion starts. The total force on the membrane is pushing, now. The  [44] filament lengths grow and forces decrease, since longer filaments are less stiff. Therefore also the detachment rate decreases and filaments re-attach. That increases the pulling forces, which after some time cause the transition back to the retraction phase. Not all of the inward forces acting on the leading edge membrane are exerted by attached filaments. Membrane tension contributes a substantial part. However, the switches in the balance between pushing and pulling forces arise from the interaction between filament attachment and length dynamics.
Our model provides a unifying oscillation mechanism across different systems. The mechanism has been most clearly shown experimentally with oil drops by Trichet et al [50]. Trichet et al describe filament detachment and re-attachment as the processes in the core of the oscillation mechanism as in our model. A decrease of filament attachment by vasodilatorstimulated phosphoprotein (VASP) can lead to an onset of oscillatory motion both in experiments and in the model [50]. The amount of pushing and pulling forces acting on the oil drop is indicated by its degree of deformation from the spherical force free shape toward a kiwi shape [50]. The total force is proportional to the velocity. That allows us to directly read the phase relation between forces and velocity from movies of oil shape and motion such as supplementary movie 2 of [50]. Oil drops show increasing deformation in the slow velocity phase [50], indicating the rise of both n d f d and n a f a . The release of these forces indicated by relaxation to a spherical shape coincides with an increase in velocity. Later, onset of deformation coincides with a decrease in velocity. This phase relation between forces and velocity agrees with our model results. The results of modeling the velocity time course of drop motion are presented in [17]. We show the phase relation between forces and velocities for drops in figure 5. As in the experiments, both n d f d and n a f a increase during the slow velocity phase, a fast increase of the velocity coincides with a drop in the amount of forces and filament detachment. The high-velocity phase is terminated by re-attachment of filaments and an increase in both n d f d and n a f a . The same phase relation applies to the oscillation mechanism of the cell simulations reported here. Merely, all velocity values are shifted toward negative velocities due to gel contraction and retrograde flow, and the values of n a f a are smaller than n d f d due to membrane tension.
An example of wave propagation with faster protrusion-retraction events as well as larger wave amplitude is shown in figure 6. In this example, we determined the wave period The membrane's retraction and protrusion velocities are almost of the same order of magnitude as the gel edge velocity, although out of phase: as the membrane retracts (v < 0), the gel boundary advances forward (v g − v r > 0) and vice versa. The transition from retractive to a protrusive state is the consequence of a significant increase in propulsive force (n d f d ) and is associated with a sudden drop in polymerization velocity v p . All the parameters are the same as in figure 3. T = 21.79 ± 4.09 at a given point i on the cell membrane as the maximum of the Fourier spectrum and subsequent averaging over all points. Similarly, we determined the spatial wavelength as time average λ = 10.42 ± 6.62 µm. In contrast to previous patterns, the wavelength is here much shorter than the cell circumference (about 1/3).
We see in figure 6 more clearly than before that protrusions and retraction events are organized in lateral waves along the cell membrane. Notably, the bands of positive (protrusion) or negative velocity (retraction) make an angle with respect to the temporal axis. This indicates that protrusion or retraction are not happening simultaneously, but rather local protrusion or retraction activity propagates with a definite speed in either direction across the cell boundary. Lateral waves propagating in opposite directions collide and annihilate each other. Interestingly, membrane waves exist even in the absence of molecular motor contractile forces (µ = 0), as can be seen in figure 2. We also looked at the dependence of wave period on myosin activity µ and adhesion strength ξ . Although the dependence of wavelength on ξ as well as µ was not significant (data not shown), reducing myosin contractile stress µ within the oscillatory regime increases the wave period significantly (figures 7(a) and (b)).
The shape of protrusion or retraction events is sensitive to the attachment/detachment rates of the filaments to the membrane. Increasing k a generates more regular wave patterns on the membrane and increases the wavelength as well (see figure 8). The retraction velocity of the membrane is faster at high values of k a due to the higher number of attached filaments. On the other hand, fast membrane retraction buckles and shortens the unbound actin filaments generating a strong propulsive force which quickly dominates the resisting membrane tension and force of attached filaments and ultimately the membrane expands again (see figure 8(b)). The type and amplitude of wave pattern depend also on the cell size and the saturation length of the cross-linking velocity. In biological terms, that saturation length can be increased by decreasing the cross-linker binding rate [24]. Increasing the cell size weakens global coupling since relative area changes by protrusions are smaller. Simulations in figure 9 were carried out with a cell with a ten times larger preferred area A 0 = 600 µm 2 (the typical size of a fish keratocyte),l increased to 400 nm and all the other parameters are the same as in figure 6. The resulting wave pattern is completely different and the wave amplitude is larger (∼1 µm).
Because membrane tension is constant along the cell boundary, it effectively couples protrusion and retraction events that take place in spatially distinct regions of the cell. Thus, membrane tension provides for a global coupling in membrane dynamics. Figure 10 shows a pattern typical of systems with global coupling [49]. It exhibits two different regions oscillating in synchrony (phase clusters) but with a phase difference between the regions. Due to an analogy   to phase clusters in systems of discrete coupled oscillators, these patterns are called cluster patterns [49].
Experiments on fibroblast and keratocyte cells show that cells spreading on a highly adhesive substrate exhibit traveling waves around the edge [26,31], but a lack of cell-substrate adhesion can lead to a 'frustrated' instability in which wave propagation is blocked along the membrane but the cell edge oscillates [26]. Within our model parameters, in the limit of weak cell-substrate adhesion strength (ξ → 0), the retrograde flow coefficient b 0 ret increases significantly although a 0 ret approaches the constant value µL/8η (see equation (10)). To  determine the effect of adhesion strength, we fixed the parameters as in figure 6, and selected b 0 ret to be a large value such as 800 nm pN −1 s −1 , which corresponds to weak adhesion. In addition, we also decreased S 1 to the value S 1 = 10 −5 pN nm −2 , which allows larger variations in the cell's area. We observed that, in agreement with experiment, the formation of protrusion and retraction events is stopped and all the points along the membrane show synchronized oscillations (see figure 11).

Discussion
The protrusions of cells spreading and resting or motile with well-developed lamellipodia consist of the posterior lamellum with highly cross-linked and bundled actin filaments and the anterior lamellipodium with a network of individual filaments polymerizing against the leading edge membrane. The circumference of the protrusions exhibits a variety of wave patterns. Alternating retractions and protrusions comprise in some cases the whole depth of the lamellipodium (see [1,3] and possibly also [7]) and have in other cases an amplitude smaller than this depth [6][7][8]. The large-amplitude oscillations may involve cyclic (partial) loss of the lamellipodium and regeneration by nucleation of filaments and are therefore not covered by this study. They will be discussed in a future publication [52]. Our simulations here apply to depth oscillations of existing lamellipodia. We have previously shown [16] that the oscillation mechanism described above and the excitability linked to it describe the experimental results reported by Machacek et al quantitatively [8], including the bifurcations between different wave regimes. There is also qualitative agreement between simulations and waves reported by Döbereiner et al [6] and Enculescu and Falcke [18]. Here, we implemented the model into a two-dimensional cell geometry, coupled it to actin gel dynamics with retrograde flow and added dynamic membrane tension as new features. We find that wave patterns also exist under these conditions. They exist for a wide range of gel parameter values (figure 2). The bifurcation diagram also suggests an explanation for the observation that in some cell types waves disappear upon inhibition of myosin and in others they are not affected. If myosin contraction is large and waves are observed (large ξ range in figure 2(a)), inhibition of myosin will terminate waves. However, if myosin contraction is small and waves are observed (intermediate ξ range in figure 2(a)), they exist for µ = 0 also.
There have been alternative explanations for cell contour wave patterns, e.g. by autocatalytic nucleation [2] or polymerization waves traveling through the bulk of the cell [51]. Membrane waves and velocity oscillations have been observed in many other cells and systems to which these two mechanisms were not applied. The similarity in molecular constituents across these different systems and cell types (at the least in function) strongly suggests that it should be possible to define a unifying concept. The mathematical model we suggest is our proposition for such a unifying concept, since it explains a large variety of experimental observations: the hopping motion of Listeria bacteria [20], the oscillatory motion of oil drops ( figure 5 and [15]), the oscillatory motion of protein-coated beads [18], waves in the cell contour [16,18] and the force-velocity relation of fish keratocytes [25] in a quantitative way. The crucial ingredients are the attachment/detachment dynamics of filaments and the dynamics of the free filament length in combination with the semiflexibility of F-actin. Attachment/detachment explains why we can observe resistance from the filaments to substantial pulling loads in some reconstituted systems such as beads and oil drops and at the same time actin polymerization-driven propulsion and velocity oscillations by the interplay of forces and attachment/detachment of filaments to/from the pushed surface [15,18]. These processes are also the core of the oscillation mechanism reported by Dickinson and Purich [53]. A simple change in the binding parameters explains the completely different force balance at the leading edge of fish keratocytes [25]. The inclusion of filaments bound to the propelled surface explains why often the net force of propulsion is much smaller than the pushing force to be expected from the observed number of filaments: some are pushing, some are holding back. The filament free length dynamics arise from the simple fact that it takes some time for newly polymerized filaments to become cross-linked into a gel. The sensitive dependence of mechanical properties of filaments on their free length [19,21] is a type of feedback to forces and detachment rates able to produce robustly all the nonlinear dynamic regimes observed in actin-based propulsion and morphodynamics [15]. Its relevance is strongly confirmed by the force velocity measurements and the elastic properties of the SR [25], as well as recent direct measurements of the free length close to the leading edge of lamellipodia by electron tomography [54][55][56][57].
In addition, the model explains qualitative experimental observations, which we have not linked in a quantitative way in explicit studies to specific experiments. In some cells exhibiting alternations between protrusion and retraction, the retrograde flow oscillates also [1]; in others it does not [1,2]. Similarly, the model may exhibit oscillations with constant retrograde flow [15] or an oscillating one (figure 3). In this study we explain in a qualitative way the different responses of wave patterns to myosin inhibition.
We think that the very mathematical realization of the model is not its most important aspect. It is rather the spectrum of biological processes and interactions defining it that has the explanatory and predictive power discussed here.