Mathematical modeling of Myosin induced bistability of Lamellipodial fragments

For various cell types and for lamellipodial fragments on flat surfaces, externally induced and spontaneous transitions between symmetric nonmoving states and polarized migration have been observed. This behavior is indicative of bistability of the cytoskeleton dynamics. In this work, the Filament Based Lamellipodium Model (FBLM), a two-dimensional, anisotropic, two-phase continuum model for the dynamics of the actin filament network in lamellipodia, is extended by a new description of actin–myosin interaction. For appropriately chosen parameter values, the resulting model has bistable dynamics with stable states showing the qualitative features observed in experiments. This is demonstrated by numerical simulations and by an analysis of a strongly simplified version of the FBLM with rigid filaments and planar lamellipodia at the cell front and rear. Electronic supplementary material The online version of this article (doi:10.1007/s00285-016-1008-2) contains supplementary material, which is available to authorized users.


Introduction
In a variety of physiological processes such as wound healing, immune response, or embryonic development, crawling cells play a vital role (Ananthakrishnan and Ehrlicher 2007). Cell motility is the result of an interplay between protrusion at the 'front' edge of the cell (w.r.t. the direction of movement), retraction at the rear, as well as translocation of the cell body (Small and Resch 2005). It only occurs when the cell is polarized with a front and a differently shaped rear (Kozlov and Mogilner 2007).
Both protrusion and retraction involve the so-called lamellipodium, a thin, sheet-like structure along the perimeter of a cell, consisting of a meshwork of actin filaments. F-actin is a polar dimer that forms inextensible filaments with a fast-growing plus (barbed) end and a slow-growing minus (pointed) end (Holmes et al. 1990).
The barbed ends abut on the membrane at the leading edge (Mogilner 2009) and have a high probability of polymerization (i.e. elongation of the filament by insertion of new actin monomers), whereas at the pointed ends mostly depolymerization (removal of one monomer) or disassembly of larger parts through severing of the filament occurs. Once a balance between polymerization and depolymerization is reached, each incorporated monomer is being pushed back by newly added monomers. Using the filament itself as a frame of reference, this can be described as movement of monomers from the barbed end towards the pointed end, a process called treadmilling (see Manhart et al. 2015a and the references therein for an overview of the involved processes and proteins). New filaments are nucleated predominantely by branching off existing filaments. The resulting meshwork is an (almost) two-dimensional array of (almost) diagonally arranged actin filaments with decreasing density towards the cell body (Small et al. 1995;Vinzenz 2012).
The lamellipodium is stabilized by the cell membrane (surrounding the entire cell Mitchison and Cramer 1996;Vallotton et al. 2005), adhesions to the substrate (Li et al. 2003;Pierini et al. 2000), cross-linking proteins (Nakamura et al. 2007;Schwaiger et al. 2004) and myosin II filaments , the latter two binding to pairs of filaments. Some of the long filaments from the lamellipodium extend into the region behind, where (through the contractile effect of myosin II) forces are generated which pull the lamellipodium backwards (Small and Resch 2005).
Fish epidermal keratocytes are fast-moving cells with a relatively simple shape (circular, when stationary and crescent-moon-shaped, when moving Lee et al. 1993), which makes them ideal subjects for analysis. Furthermore, they exhibit a lamellipodium with a smooth edge and a fairly uniform distribution of filaments (Lacayo et al. 2007;Small and Resch 2005;Theriot and Mitchison 1991). During the transition from the stationary to the moving state, the lamellipodium in the rear of the cell collapses and the rear bundle is formed, where myosin II generates a contractile force Tojkander et al. 2012;Verkhovsky et al. 1997).
Treatment with staurosporine (a protein kinase inhibitor) results in the formation of completely detached lamellipodial fragments, lacking a cell body, microtubules and most other cell organelles. Remarkably, these fragments can either remain stationary while adopting a circular shape, or can move on their own, adapting their appearance to the same crescent-moon shape as the keratocyte itself (Kozlov and Mogilner 2007;Fig. 1 a A moving keratocyte (right) and a moving cytoplast (left), actin is labelled in green, the nucleus in blue. b, c A moving and a stationary cytoplast (fragment), respectively. The actin network is labelled in red, myosin in green. a, b, c Are reproduced from (Manhart 2011). e Idealization with protruding lamellipodium at the top and lamellipodium collapsed by actin-myosin interaction at the bottom. d Model ingredients of the simplified FBLM (clockwise, starting top left) cross-link stretching, cross-link twisting, filament-substrate adhesion, connection between front and rear by stress fibres, membrane stretching, actin-myosin interaction (color figure online) Verkhovsky et al. 1999) (see Fig. 1a-c). This suggests that the necessary ingredients for movement are all present in the lamellipodium (until it runs out of energy).
Various approaches to continuum mechanical modeling of the lamellipodium exist (see e.g. Kruse et al. 2006;Rubinstein et al. 2005). Bistability results similar to this work have been obtained in Ziebert et al. (2012), where a phase-field approach is used to describe the interplay between the cell shape, the mean orientation of the filaments within the network and the actin-myosin interaction. A strongly simplified model is developed in Kozlov and Mogilner (2007), where bistability has been obtained analytically as the consequence of the properties of a free energy functional containing contributions from the lamellipodium and a possible rear bundle.
This work is based on the FBLM (Manhart et al. 2015a;Ölz and Schmeiser 2010a, b), a two-dimensional, anisotropic, two-phase model derived from a microscopic (i.e. individual filament based) description, accounting for most of the phenomena mentioned above. It describes the actin network in terms of two transversal families of locally parallel filaments, stabilized by transient cross-links and substrate adhesions. In Sect. 2 the FBLM is presented and extended by a model for actin-myosin interaction between the two families. We assume that myosin filaments can connect only when the families are anti-parallel enough and they are described as transient, similar to cross-links. They tend to slide the two families relative to each other, and they are assumed to have a turning effect, making the two families more anti-parallel. The derivation of the additional myosin terms is presented in detail. Readers not consulting (Manhart et al. 2015a;Ölz and Schmeiser 2010a, b) may take this as representative also for the derivation of the adhesion and cross-link models.
The properties of the actin-myosin model are expected to produce the desired bistable behavior. This is demonstrated by numerical simulations in Sect. 8, which indicate the existence of two stable states, a rotationally symmetric nonmoving state and a polarized state, where the cell moves. The moving state is characterized by a more anti-parallel network in the rear of the cell, where actin-myosin interaction is active. Complete collapse of the network and consequential generation of a rear bundle are avoided, since the FBLM is (so far) unable to deal with such topological changes.
The occurrence of bistability is also proven analytically for a strongly simplified model. In Sect. 3 the complexity of the model is reduced in a first step by assuming rigid filaments. Then a planar, translationally invariant lamellipodium is considered in Sect. 4, which reduces the model to a system of three ordinary differential equations. Here we also neglect the effects of branching and capping, assumed to be in equilibrium, as well as filament severing within the modelled part of the lamellipodium, implying a constant actin density there. Bistability is obtained for this model in Sect. 5. Finally, in Sect. 6 a cell (fragment) is replaced by a pair of connected back-toback planar lamellipodia, and the existence of stable stationary (symmetric) as well as moving (polarized) states is proven. The same bistable behavior is observed in the simulations of the full model in Sect. 8. Figure 1 depicts the main components of the simplified version of the FBLM (D and E) together with one keratocyte and three fragments (A-C). The crescent-moon shaped cells and cell fragments are moving, whereas the circularly shaped fragment remains stationary. One can also observe that in moving fragments myosin can predominantely be found at the cell rear. In Fig. 1e, the idealized model obtained in Sects. 3, 4, 5, 6 is illustrated. It can be interpreted as description of lamellipodial sections at the front and at the rear of the cell. The main model ingredients are depicted in Fig. 1d: (e) diagonally arranged filaments (red) together with the membrane (dark green) and arrows indicating inward pulling forces due to stress fibers in the interior of the cell (dashed green line), (d) the cell membrane (green, with arrows indicating the force acting on the barbed ends due to membrane tension), (a) cross-links [blue, producing friction between the filament families and a turning force trying to establish an equilibrium angle-arrows show the forces acting on the cross links due to resistance to stretching and bending above or below a certain angle (visualized by the black dashed line)], (b) myosin filaments (pink, trying to slide the filament families and to make them anti-parallel-straight arrows indicate that myosin moves towards the plus ends along both filaments while curved arrows illustrate the tendency to establish an angle of 180 • between the filaments, analogously to the cross-links), (c) adhesions (yellow, connecting a filament through the membrane with the substrate thus producing friction relative to the substrate).

Adding actin-myosin interaction to the Filament Based Lamellipodium Model (FBLM)
Our starting point is the FBLM as introduced in Ölz and Schmeiser (2010a) (see also Manhart et al. 2015a): where F = F(α, s, t) ∈ R 2 describes the position and deformation of actin filaments in the plane at time t. More precisely, the variable α ∈ A ⊂ R, for some interval A, is a filament label, and s ∈ [−L(α, t), 0] denotes an arclength parameter along filaments, which means that the constraint |∂ s F| = 1 has to be satisfied. Here L(α, t) is the maximal length of filaments in an infinitesimal region dα around α. The filament length density with respect to α and s is given by η(α, s, t), which will be assumed as given (see Manhart et al. 2015a for a dynamic model incorporating polymerization, depolymerization, nucleation, and branching effects). The value s = 0 corresponds to the so called barbed ends of the polar filaments, abutting the leading edge of the lamellipodium. The rear boundary s = −L(α, t) is introduced somewhat artificially since the rear end of the lamellipodium is typically not well defined. By polymerization with speed v(α, t) (also assumed as given in this work), monomers move along filaments in the negative sdirection. Their speed relative to the nonmoving substrate is therefore given by D t F with the material derivative D t = ∂ t − v∂ s . The terms in the first line of (1) correspond to the filaments' resistance against bending with stiffness parameter μ B , to friction relative to the substrate as a consequence of adhesion dynamics with adhesion coefficient μ A , and to the constraint (2) with the Lagrange multiplier λ inext .
The FBLM is actually a two phase model, and F may stand for either of the two families F + or F − . The terms in the second line of (1) describe the interaction between the two families, with the other family indicated by the superscript * . The interaction is the consequence of dynamic cross-linking and leads to a friction term proportional to the relative velocity between the two families and to a turning force trying to push the angle ϕ (cos ϕ = ∂ s F · ∂ s F * ) between crossing filaments to its equilibrium value ϕ 0 , corresponding to the equilibrium conformation of the cross-linker molecule (F ⊥ = (−F y , F x )). The * -quantities corresponding to the other family have to be evaluated at (α * , s * ), determined by the requirement F(α, s, t) = F * (α * , s * , t). It is a basic geometric modeling assumption that the coordinate change (α, s) ↔ (α * , s * ) is one-to-one, wherever the two families overlap. It requires that filaments of the same family do not cross each other and that pairs of filaments of different families cross each other at most once. Finally, the coefficients are given by with constants μ S,T , wherever F crosses another filament, and zero elsewhere. The partial derivative refers to the coordinate transformation introduced above. The FBLM will be extended by the effects of myosin polymers. The basic assumption is that pairs of crossing actin filaments, which lie antiparallel enough, may be connected by a bipolar myosin filament. The modeling is similar to that of cross-links and we will give details of the derivation of the myosin related term (originally presented in Manhart 2011) as an example of the general modeling strategy of the FBLM. We assume the following about myosin molecules: • Myosin filaments connect pairs of crossing actin filaments.
• Due to the motor activity of the myosin heads they "walk" towards the barbed ends of actin filaments with a fixed speed v M . • The rates of creation and breakage of the myosin connections depend on the forces acting on them (see below). • The connections exert twisting forces on the connected actin filaments towards the equilibrium angle π , i.e. the antiparallel state. These forces are caused by the stiffness of myosin filaments. • The connections can be stretched against an elastic restoring force.
The key to the bistability results presented in the next sections is the assumption that the myosin activity (specified below) depends on the angle between crossing actin filaments. This is supported by data reported in Reymann (2012), where it has been shown that myosin can act much more efficiently on antiparallel actin networks as opposed to branched networks. In the original model of the myosin dynamics presented in Manhart (2011), the growth of myosin filaments was also included. Here, we omit myosin size effects for the sake of simplicity. The derivation is done in five steps.
Step 1: Energy contributions Let a filament with label α have a crossing with a filament from the other family with label α * . By s(α, α * , t) we denote the s-position along the α-filament where this crossing occurs. If at time t − a a myosin filament has been attached at this crossing (i.e. the connection has age a at time t), the s-value of the binding site will change due to the aforementioned treadmilling effect (caused by actin polymerization) and the myosin motor activity. Therefore at time t the myosin filament binds to the actin filament at Analogously, the position s * a (α, α * , t) at time t of the binding site on the α * -filament is computed. With the aid of s a and s * a we can now determine the deviations from the unstretched and untwisted equilibrium state of the myosin filament: Note that S M 0 = 0 and T M 0 = ϕ − π , i.e. the myosin filament is unstretchted but possibly bent at the time a connection is created. Now let ρ M (α, α * , a, t) denote the probability density with respect to age a of an active myosin connection between the filaments with label α and α * at time t. Following classical modeling for age-structured populations we assume ρ M to satisfy The rate of breakage ζ may depend on the stretching and twisting forces acting on the myosin filament. We also allow for the possibility that the rate of creation β is affected by how favorable the angle between the actin filaments is. The interpretation of the last factor is that a new connection can only connect a so far unconnected pair of actin filaments. The main Eq. (1) of the FBLM is the result of a variational approach applied to the problem of minimizing a sum of potential energies. Each energy contribution produces forces acting on the actin filament (i.e. bending stiffness, adhesions, etc.). To include new terms, it is therefore necessary to formulate the corresponding energy. For the elastic stretching and twisting of the myosin filaments, we choose the form The integration domain C(t − a) includes all pairs of filaments (α, α * ), which had a crossing at time t − a. The constants κ SM and κ T M are elasticity coefficients.
Step 2: Scaling The key scaling assumption of the original FBLM is that the average lifetime of a cross-link or adhesion is small compared to the average time a monomer spends inside a filament, a ratio denoted by ε. In the limit ε → 0 the non-locality in time is removed from the problem. Furthermore this reduces the effect of adhesions to friction with the substrate and that of the cross-linkers to friction between the filament families (compare Eq. (1)). This raises the question of what to assume about the average lifetime of a myosin filament. In Svitkina et al. (1997) it was observed that in stationary cytoplasts myosin spots stay small and disappear after some time, indicating their transient nature. The situation is quite different for moving cytoplasts, where the myosin spots grow in size over time. Since the aim of this paper is to describe the onset of movement, we use here the same scaling assumption as for cross-links, i.e. that the lifetime of a myosin filament is relatively short. Without going into further detail about the precise reference values for all variables and parameters, we state the resulting energy contributions and equations after scaling: The density equations take the form we can explicitly calculate the solution to the limiting equation. In the following we will replace dependencies on T M 0 by dependencies on ϕ. Ignoring a possible initial time layer, the limiting solution is given by Note that ρ M 0 depends on t and the indices α and α * via the angle between the filaments. The scaled energy contributions take the form The factor 1/ε in front of the stretching contribution is the result of a scaling assumption. It ensures that the effect of myosin filament stretching does not vanish in the limit, although the energy contribution itself does.
Step 3: Variation and macroscopic limit The next step in the derivation is to calculate the variations δU [F, F * ]δ F of the energy contributions, considering ρ M ε given at this stage. After passing to the limit ε → 0, we obtain where D M t F denotes the velocity of a myosin binding site on the actin filament relative to the substrate, with The stiffness parameters μ SM and μ T M result from using the representation of ρ M 0 given in (5) to evaluate the integrals with respect to the myosin age a. They take the ζ(0, ϕ) .
Motivated by the findings in Reymann (2012), we assume that there exists a cutoff angle ϕ < π such that Step 4: Euler-Lagrange Equations The final step in the derivation is to formulate the corresponding Euler-Lange equations. Together with the original terms given in (1), the modified model takes the form with The introduction of μ SM and μ T M is a consequence of the mapping between C and {(α, s) : α ∈ A, s ∈ [−L(α, t), 0]}. Further details of the model derivation can be found in Manhart (2011).
Step 5: Boundary conditions describe the forces acting on the filaments at their barbed ends and at the artificially introduced ends at the boundary of the modeling domain: Thus, there are no torques applied at the ends. The choice of the linear forces f 0 and f L along the leading edge and, respectively, along the artificial boundary will be discussed later.

Rigid actin filaments in the limit of large bending stiffness
We want to derive a simplified model with rigid actin filaments. This is motivated on the one hand by the observation that filaments within the lamellipodium are typically rather straight (Vinzenz 2012). On the other hand stiff filaments can be interpreted as a description of only the outermost part of the lamellipodial region, where filaments are (locally) straight. The resulting model is mathematically much simpler and can be derived by assuming a relatively large bending stiffness μ B . The limit μ B → ∞ will be carried out formally in this section. The solutions of the formal limit 0 = ∂ 2 s η∂ 2 s F of (6), together with the boundary conditions and with the constraint (2), can be written as where s 0 is determined by In other words, F 0 is the center of mass of the filament, and d(ω) its direction. The components of F 0 and the angle ω are still to be determined. The total force balance obtained by integration of (6) with respect to s and using the boundary conditions (8) reads Note that it does not contain μ B and therefore remains valid in the limit. Similarly, the total torque balance is obtained by integration of (6) against (F − F 0 ) ⊥ : This completes the formulation of the rigid filament version of the FBLM. Substitution of (9) into (10) and (11) gives a system of ordinary differential equations for F 0 and ω. Note that coupling with respect to α happens only indirectly through the interaction between the two filament families.

A geometric simplification: the planar lamellipodium
Since in keratocytes the leading edge is rather smooth, we approximate a piece of lamellipodium by an infinite strip, parallel to the x-axis, and invariant to translations and to reflection. For the given data this means that the maximal filament length L and the polymerization speed v are constants. As a further simplification, we assume no filament ends inside the modeled part of the lamellipodium with the consequence η = 1 (and s 0 = −L/2). We assume two families of rigid filaments (9) with giving F ± (α ± , s ± , t) = ±x(t) + α ± ± (s ± + L/2) cos ω(t) The angle between two crossing filaments and the coordinate change between the two families mentioned in Sect. 2 are easily computed: It provides the geometric quantity needed in (3) and (7): This quantity can be interpreted as a measure of the density of crossings, with a maximum at ω = 0 (fully collapsed lamellipodium) and a minimum at ω = π/2 (all filaments are parallel, no crossings). With the planar lamellipodium ansatz, the Eqs. (10) and (11) become independent of α and constitute a system of three ordinary differential equations for the unknowns (x(t), y(t), ω(t)):

Forces at the filament ends-steady protrusion
The membrane stretched around the lamellipodium exerts a force on the polymerizing barbed ends. On the other hand, we assume that the filaments at the rear of the lamellipodium are connected to stress fibres pulling them backwards, another consequence of actin-myosin interaction. Both the membrane force and the stress fibre force will be described as acting in the negative y-direction orthogonal to the leading edge, i.e.
If these forces are modeled as constant, the Eq. (14) for the angle is decoupled from the remaining system. For an analysis of its dynamic behavior, we choose a model for the stiffness coefficients of the actin-myosin connection: with μ SM , μ T M > 0, ϕ 0 < ϕ < π, and with the notation (.) + for the positive part. Bistability can now be obtained with appropriate assumptions on the parameters. The right hand side of (14) can be written as It is a simple exercise to prove: then h(ω) as defined above has three simple zeroes ω 10 , ω 20 , ω 30 , satisfying Theorem 2 Under the assumptions of Lemma 1 and for | f stress − f mem | small enough, the ordinary differential Eq. (14) with the forces given by (15) possesses four stationary solutions ω j , j = 0, . . . , 3 with where ω 0 and ω 2 are unstable, and ω 1 and ω 3 are asymptotically stable.
Again the proof is straightforward. For the stable steady states, the lamellipodium has the constant protrusion speedṡ For the equilibrium angle ω 1 , we typically expect the speed to be positive. It is not affected by actin-myosin interaction. The smaller speed corresponding to ω 3 might actually be negative due to membrane tension and stress fibres, i.e. the second stable state, where the lamellipodium is collapsed by actin-myosin interaction, might be retractive. Finally, the steady states also produce lateral flow with constant speedṡ cos ω 3 , respectively, where in the collapsed state the lateral flow speed produced by polymerization is reduced by actin-myosin interaction.

Coupling of two opposing lamellipodia-bistability
As a caricature of a cell fragment, we consider two back-to-back planar lamellipodia (see Fig. 1e). For notational convenience, the bottom lamellipodium is rotated by 180 • in the mathematical description. Therefore we consider two versions of the system (12)-(14) with unknowns (x, y, ω) and (x,ŷ,ω). The assumption that the total forces exerted on the fragment by membrane tension and by stress fibres vanish, imply that (15) is used in both systems with the same values for f mem and f stress . However, we allow the option that these forces are not constant but regulate the size of the fragment, measured by y +ŷ. We first consider the case of a constant given membrane force and a size dependent force by stress fibres: Case A : f mem = const, f stress = f stress (y +ŷ).
Typically f stress will be an increasing function, but the details are not important for our considerations. Adding the Eq. (13) for y andŷ leads to a closed system of three equations for y +ŷ, ω, andω: with g(ω) = L 2 24 μ A + 4 sin 2 ω cos ω(μ S + μ SM (π − 2ω)) .
We shall prove that with appropriate assumptions on the data, the problem has 4 stable steady states.

Theorem 3 Let the assumptions of Lemma 1 hold, let the function f stress be continuously differentiable with bounded positive derivative, and let f mem , μ A vL, and the Lipschitz constant of f stress be small enough. Then the system (18)-(20) has four stable steady states, satisfying
where the purely cross-link dominated state ω CL and the myosin-influenced state ω MY are given by: Proof From (18) we obtain that steady states have to satisfy This implies, again for stable steady states, The existence of the four steady states is then a consequence of a straightforward perturbation argument. The coefficient matrix in the linearization of (18)-(20) can be written as with positive constants A andÂ, and with 0 < κ = f stress (y +ŷ) 1. A perturbation analysis of the eigenvalue problem for small κ (i.e. formal expansion of eigenvalues in terms of powers of κ and subsequent justification by a contraction argument) gives the eigenvalues which are all negative at the four steady states for small enough κ, because of h (ω 10 ), h (ω 30 ) < 0.
For the steady states the protrusion speed of the fragment is constant and given bẏ For the symmetric steady states (21), (22), the protrusion speeds vanish, hence they describe stationary cells (or fragments). The equilibrium angles in the lamellipodia in this case are either both affected by myosin, (22), or both result only from crosslink activity (21). The asymmetric steady states (23) and (24) describe a protruding, polarized cell. In both cases it consists of a collapsed cell rear, in which myosin is active (ω = ω MY ), and a cell front with a steeper equilibrium angle caused only by cross-link activity (ω = ω CL ). Finally, we also mention the case of a constant stress fibre force and a size dependent membrane force: Case B : f stress = const, f mem = f mem (y +ŷ).
Without going through the details, we note that the qualitative results are the same and a theorem analogous to Theorem 3 can be proven.

Parameter dependencies
The simplifications of the two preceding sections lead to several (testable) statements. This section can be seen as a discussion of the results so far.

Protrusion speed versus adhesion strength
For the coupled lamellipodium of Sect. 6 the protrusion speed is given by (26). For a moving fragment, i.e. situation (23) or (24), we compute Substitution of (25) into the stationary version of (19) (resp. (20)) shows that dω CL /dμ A > 0. On the other hand, the difference of the stationary versions of (19) and (20) gives dω MY /dω CL < 0. Thus, the protrusion speed increases with adhesion strength. Since the results of Theorem 3 are restricted to relatively small adhesiveness, this is in agreement with experimental results (Barnhart et al. 2011) showing two regimes: for smaller adhesion strengths the cell speed is positively correlated with the adhesiveness of the ground. After reaching a maximum speed, the correlation is reversed for larger adhesion strengths. It is not clear if the FBLM in its present form is able to describe the second regime, which might be caused by long-lived focal adhesions, contradicting the assumption of rapid adhesion turnover. It should also be noted that it has been hypothesized (see Barnhart et al. 2011) that the reduction of cell speed on strongly adhesive ground is mainly due to biochemical reasons as opposed to purely mechanical effects. Figure 2 illustrates these considerations.
Onset of myosin activity The inequality (17) gives a criterion for the existence of a myosin influenced stable steady state. Myosin has to be strong enough compared to forces maintaining the branched actin network structure (e.g. Arp2/3, filamin, etc.). Only then is it able-either spontaneously or aided by a pushing force-to initiate the collapse of the network, which precedes bundle formation. Data like those presented in  Reymann (2012) show that myosin activitiy is influenced by actin network structure. The model for myosin activity (16), used for the analyis in this work, can be interpreted as a limiting case, in which the activity is set to zero for unfavorable angles. From the approximation of (17) for ϕ close to π , it is easily seen that if the interval for the myosin angle is small, this needs to be compensated by a large stiffness of the myosin filament in order to initiate symmetry breaking. In our model myosin is necessary for the onset of movement, i.e. polarization of originally symmetric fragments. In Yam (2007) the number of stationary cells spontaneously initiating motility within a given time interval was reduced from 45 to 10 % by inhibiting myosin using blebbistatin, which supports our findings.
Protrusion speed versus myosin strength In the presented model two parameters characterize myosin strength, the elasticity coefficients of myosin filament stretching μ SM and bending μ T M . Whereas μ SM has no influence on the equilibrium angle and therefore also not on the protrusion speed, larger values of μ T M lead to faster cells (easily seen in the expression for h(ω)). This is consistent with experimental results found in Barnhart et al. (2011) where myosin contraction was either inhibited using blebbistatin (leading to slower cells at least for low to medium adhesion strengths) or enhanced using calyculin A (leading to faster cells in all adhesive regimes).

Simulations with the full model
In this section we demonstrate that with the additional term describing myosin within the lamellipodium, the model is able to produce cells/cell fragments that, depending on the initial conditions, will either remain stationary or start moving. In contrast to the simulations presented in Manhart et al. (2015a, b), here the movement is achieved without a continuing external signal and without varying the polymerization speed. In the simulation, we work with the full model (6)-(8) and not with the simplifications introduced in Sects. 3 and 4. However, the qualitative results of Sect. 6 will be reproduced.
Parameter values Parameter values are chosen as in Manhart et al. (2015a) with the following exceptions and additions: we work with a constant filament density η = 1 in parameter space, which means that the filament number remains constant with branching and capping always in equilibrium. No pointed ends appear within the simulation region, which corresponds to a fixed filament length of L = 8 µm. The polymerization speed is fixed at the constant value v = 3 µm min −1 . In Svitkina et al. (1997) it has been observed that myosin speckles that are formed in the lamellipodium drift inwards with time. This indicates that the myosin velocity has to be smaller than the polymerization speed. We therefore chose v M = 1 µm min −1 . Motivated by Reymann (2012) we assume that myosin can only act on actin filaments if the angle between the filaments is larger than ϕ = 100 • . For the stiffness parameters of stretching and twisting the cross-links and myosin good estimates are hard to obtain, since their exact concentration in the lamellipodium is difficult to determine. Motivated by Lemma 1 we chose the myosin twisting force larger than the cross-link twisting force (see Table 1), but of the same order of magnitude. Additionally we increased the bending stiffness by a factor 10 in order to get closer to the analytical case examined in Sect. 3. Here we concentrated on the case of size control via stress fiber forces from the inside, i.e. the membrane force was set to zero. The stress fiber force was chosen to be f stress = μ I P (A − A 0 ) + , where A is the current size of the area surrounded by the lamellipodium, A 0 is the equilibrium inner area and μ I P is the corresponding (inner pulling) force constant. In order not to distort the direction of the filaments, we let the  Manhart et al. (2015a), i.e. 36 discretization nodes in α direction and 9 nodes in s-direction. The timestep used is 2 × 10 −3 min. Figure 3 shows the evolution of the cell in two different numerical experiments over a timespan of 15 min. In both cases identical initial conditions were used, in which the left side of the cell is deformed. This could be caused either by internal fluctuations or by a pushing force from the left, two situations known to lead to symmetry breaking in fragments (Verkhovsky et al. 1999). Both simulations use the same set of parameters, with the exception of the value of the myosin twisting constant, μ T M , which is smaller in the left column than in the right column. In Fig. 3ad it can be observed that even though the filaments on the left initially lie anti-parallel enough for myosin to act on them, myosin does not establish itself there permanently. Lemma 1 gives a possible explanation for this: If myosin is too weak compared to the cross-links, there is no myosin/cross-link-equilibrium, hence the cell reverts back to a purely cross-link dominated state. This situation corresponds to the steady state (21) in Theorem 3. The case of a larger myosin twisting constant is depicted in Fig. 3e-h. In situation case a bundle-precursor is formed on the left, whereas the right stays myosin free. This allows the cell to change to a moving steady state (see the Supplementary Material for a movie) with a collapsed lamellipodium at the back and a non-collapsed lamellipodium at the front, a situation described by the steady states (23) and (24) in Theorem 3. One can also observe the contraction of the rear bundle leading to a more half-moon shaped cell. Clearly if the initial deformations are so small that no myosin can attach, the cell always reverts back to the stationary state. This refers to a situation where the fluctuations or pushing force are too small to cause symmetry breaking.

Discussion and outlook
In this work, we extended the FBLM introduced in Manhart et al. (2015a), Ölz and Schmeiser (2010a) by a description of actin-myosin interaction. The limit of large bending stiffness led to a simplified model for rigid filaments. The additional simplification of a planar lamellipodium reduces the model to a small ODE system for the center of mass of a reference filament and the angle with the membrane. A caricature of a cell fragment has been described by considering two versions of the model coupled by membrane and stress fiber forces. Bistability has been shown in the asymptotic regime of relatively small coupling forces and adhesion strength. Since these results are almost explicit, various parameter dependencies could be obtained. Furthermore, we carried out numerical simulations based on the full FBLM supporting the main analytical result of bistable behaviour.
Our model is able to qualitatively reproduce the observed bistability of cells and cell fragments (Verkhovsky et al. 1999). Since myosin has been shown to be effective only if the angle between the filaments is large enough (Reymann 2012), our specific modeling of the myosin effect to be angle-dependent seems to be a reasonable choice.
In the numerical simulations we have shown that for an initially slightly asymmetric cell, one of two stable steady states is attained, depending on the parameters for myosin Fig. 3 Cell view with clockwise filaments in blue and anti-clockwise filaments in red. Green stars in the lamellipodium mark the area where myosin is active. The black dot marks the origin, the black star the cell's center of mass. The insets show the angles between the filaments of the different families, averaged along each filament, parametrized along the membrane with 0 being at the very right and going counterclockwise. Blue (dashed) lines refer to the clockwise filament family, red (solid) lines to the anti-clockwise filament family. The horizontal green line at 100 • marks the myosin cut-off value ϕ. a-d Time series with μ T M = 1.4 × 10 −2 pN µm. Myosin is initially active in the left part of the cell, but is not strong enough to stay there. Eventually the cell goes back to the stationary cross-link dominated equilibrium. e-h Time series with μ T M = 1.8 × 10 −2 pN µm. Myosin is initially active in the left part of the cell, where a myosin/cross-link equilibrium emerges. Since the right of the cell is unaffected by this, the cell moves to the right. Parameters as in Table 1 (color figure online) and cross-links. If myosin is too weak to exert considerable twisting forces, the effects from cross-links on the local angles between the filaments dominate and the cell reverts back to a symmetric nonmoving shape. On the other hand, if myosin forces are strong enough, a bundle precursor forms at the 'rear' (the location of initial asymmetry), and the cell starts moving. What is apparent in these simulations is that myosin tends to locate at the back of the lamellipodium (away from the membrane). This is also in good agreement with experimental findings. It is remarkable that both for the analytical setting as well as in the numerical simulations, movement is an self-organized behavior, once an initial asymmetry has been established. It is known that keratocytes exhibit little if any chemotaxis, i.e. directed movement to outward cues.
An important issue remains to be addressed: The analytical model, dealing only with pieces of lamellipodium at the front and at the rear of the cell, avoids the transition zone separating the two parts of the lamellipodium with an intact network on the one hand and, on the other hand, compression to a rear bundle. It can be expected that due to lateral flow, actin filaments are drawn into the bundle in these transition regions. However, we expect that this happens with the pointed ends first entering the bundle, as opposed to the simulations presented here. A model supporting the necessary change of orientation in the transition zone is not available so far, and its derivation is the subject of ongoing work.