One dimensional cell motility patterns

During migration cells exhibit a rich variety of seemingly random migration patterns, which makes unraveling the underlying mechanisms that control cell migration a daunting challenge. For efficient migration cells require a mechanism for polarization, so that traction forces are produced in the direction of motion, while adhesion is released to allow forward migration. To simplify the study of this process cells have been studied when placed along one-dimensional tracks, where single cells exhibit both smooth and stick-slip migration modes. The stick-slip motility mode is characterized by protrusive motion at the cell front, coupled with slow cell elongation, which is followed by rapid retractions of the cell back. In this study, we explore a minimal physical model that couples the force applied on the adhesion bonds to the length variations of the cell and the traction forces applied by the polarized actin retrograde flow. We show that the rich spectrum of cell migration patterns emerges from this model as different \emph{deterministic} dynamical phases. This result suggests a source for the large cell-to-cell variability (CCV) in cell migration patterns observed in single cells over time and within cell populations. The large heterogeneity can arise from small fluctuations in the cellular components that are greatly amplified due to moving the cells' internal state across the dynamical phase transition lines. Temporal noise is shown to drive random changes in the cellular polarization direction, which is enhanced during the stick-slip migration mode. These results offer a new framework to explain experimental observations of migrating cells, resulting from noisy switching between underlying deterministic migration modes.


I. INTRODUCTION
Eukaryote cell migration, whereby cells crawl actively over an external substrate, is a subject of great interest for biological processes such as development and cancer progression.Adhesion-based cell motility involves the orchestration of a large number of cytoskeletal proteins: Typically, the cell needs to break its symmetry (polarize), and produce traction forces in the direction of polarization (balanced by drag or friction forces).The traction forces are mediated by adhesion between the cell and the external substrate, but these adhesions have to detach at the trailing end of the cell in order to allow the cell to migrate forward.When observing freely migrating cells (i.e.not guided by an external gradient of any kind), one is often struck by the large variation in their migration patterns [1][2][3][4][5], both for a single cell over time, and for a (seemingly identical) cell population.The origin of the large Cell-to-cell variability (CCV), or phenotypic, population heterogeneity, observed during cell migration is not understood [6], and is usually ascribed to the inherent noise of cellular systems [7,8].
To simplify the study of the complex process of cell motility, cells can be confined to move along onedimensional tracks, either on flat adhesive stripes [9,10], linear grooves [11] and channels [12,13], or thin fiber [14,15].In addition to being a simple geometry for the study and analysis of the cell motion in general, such confined motion appears also in-vivo [16], for example when cells move along axonal fibers [17,18], or cancer invades in confined spaces between tissues [19].
The experiments listed above have shown that iso-lated cells on one-dimensional tracks exhibit the following stereotypical behaviors [9,10]: (i) Non-migrating and unpolarized, by remaining quiescent or elongating symmetrically, (ii) Undergoing spontaneous symmetry breaking, polarization and migrating smoothly, (iii) As in (ii) but exhibiting stick-slip migration.The stick-slip motility mode is characterized by protrusive motion at the cell front, coupled with an overall elongation of the cell and followed by rapid retraction of the cell back.In both (ii,iii) the cell motility can be highly persistent, or undergo sporadic direction changes.The appearance of this variety of migration behaviors within a uniform cell population, as well as the switching of cells between these modes, remains an open puzzle and is the aim of this work.Clearly the diversity of the migration patterns is beyond modeling the cell motility as a simple random walk process [4,[20][21][22].
We present a theoretical model that describes the cell motility in a highly simplified, and coarse-grained manner.Nevertheless, the model contains two key, and strongly coupled, components: a slip-bond adhesion module at the cell back, and a cellular polarization module.We find that these components are sufficient to drive the entire spectrum of observed motility patterns, and explain the transitions between them.Our work therefore demonstrates how a minimal model gives rise to a rich variety of deterministic migration patterns.In experiments these deterministic patterns may underlie the migration of cells, but further confounded by noise.
The stick-slip motion is shown here to have an underlying deterministic oscillatory behavior, separated from smooth migration by a bifurcation line.Unlike previ-ous treatments of stick-slip dynamics of cells, which focused on the protrusion-retraction cycles at the cell edge [23,24], we couple the stick-slip adhesion at the cell back to the overall cell polarization, through the dependence of the polarization on the cell length.This dependence, which is an inherent property of the UCSP model [25], was not previously explored.Deterministic oscillations in the speed of migrating dendritic cells, for example, were related to competition for finite resources that directly affected the acto-myosin polarization mechanism, but did not involve length oscillations [26].Furthermore, the dendritic-cell oscillations depend on a specific macropinocytosis process, while here we obtain deterministic oscillatory (stick-slip) migration patterns that are driven by adhesion dynamics that are much more general across cell types.

II. MODEL
The model is introduced by three parts, with increasing levels of complexity and realism.In this manner we expose the motility patterns and the key components that drive them.
The first part describes a cell that is constantly polarized, with a constant protrusive activity at the leading edge, and slip bond adhesions at the rear [24].This part allows us to expose the oscillatory stick slip behavior through the dynamics of the cell length and adhesion concentration at the rear.
The second part adds a self-polarization model to the polarized cell [25].This part couples the dynamics in cell length to the protrusive activity, and introduced a critical polarization length scale.However the model simplifies the cell to have a single leading edge, where all the protrusive activity is concentrated.
In the third part the model is extended to be symmetric, such that the protrusion and adhesion dynamics acts on both edges of the cell.This part outlines the conditions for symmetry breaking and the role of noise in choosing a migration direction.

Model description
Consider a cell of length l that migrates along a linear track.The two ends of the cell, the front and back, denoted by x f and x b , are connected by a spring (Fig. 1A).The stiffness of the spring k, represents the effective elasticity of the cell cytoplasm and membrane.Such a mechanical coupling between the front and rear was recently demonstrated experimentally [27].
In this part of the model, the motion of the cell is considered to be already polarized, such that the actin treadmilling from the front to the back occurs at a constant velocity v, and produces a constant protrusive force at the front where α describes the strength of the coupling of the actin retrograde flow and the effective friction force generated by adhesions which grip the sliding filaments.This is a simplified representation of the "clutch" mechanism [23], that converts the sliding of the actin to a protrusive force that pushes on the membrane.This coupling is dependent on the adhesion strength, and we therefore expect a term of the form: r/(r+r 0 ) to multiply the r.h.s. of Eq.1, where r is the ratio between the binding/unbinding rates of the cell-substrate adhesion molecules r = k on /k 0 of f , and r 0 quantifies the the cell-substrate adhesion saturation.Such a term accounts for the loss of traction force when the adhesion diminishes r → 0. We do not explicitly describe here the catch-bond property of these adhesions, as we are not interested in the stick-slip dynamics of the leading edge, but rather wish to focus on the stick-slip events on the whole cell scale.For the rest of this paper we effectively work in the limit of r 0 r (the effects of r 0 are shown in Fig. 23).
The pushing force is balanced by a drag force which is proportional to the speed of the moving front, and a restoring force due to the global cell elasticity where γ represents the effective resistance to the motion of the cell front due to the friction generated by the contact of adhesion molecules with the substrate.As in Eq.( 1), we expect a term of the form r/(r + r 0 ) to multiply the first term on the r.h.s. of Eq.2, but we neglect it by choosing to work in the limit of r 0 r.The parameter k represents the stiffness of the spring which describes the effective elasticity of the cell, and the parameter x 0 represents the rest length of the cell.
Equating ( 1) and ( 2) due to the force balance at x f yields which allow to obtain the equation of motion for the moving cell front At the back, x b , the pulling force due to the cell elasticity is balanced by a friction force which results from the stretching of bound linkers, which model slip bond adhesions where n is the mean number of bound linkers, κ is the spring constant of the linkers, and ∆x is the average displacement of the linkers (Fig. 1C).The number of bound linkers evolves dynamically and obeys the following kinetics where N is the total number of linkers, k on is the basal attachment rate, and k of f is a detachment rate which depends exponentially on the force applied on each bound linker [23,24] where k 0 of f is the basal detachment rate and f s represents the susceptibility of the linkers to the applied force [24].The stretching force per bound linker is The average stretch of each linker ∆x depends on the dissociation rate of the linkers (8,7) and can be written as ∆x = ẋb k of f (10) By combining (8,9,10) into the force balance between ( 5) and ( 6) we obtain where κ is the effective spring constant of the linkers, and we substituted for the average stretching of each linker: ∆x = ẋb /k of f .By reorganizing (11) we obtain the equation of motion of the moving back Combining ( 4) and (12), and changing coordinates to Next, the system (13,14) is rescaled by the time and length scales of k −1 of f ,x 0 , and by the total number of adhesion sites N , as well as rescaled by the parameters Finally, we set α = 1 and γ = 1 and remove the tilde signs in Eqs.(15,16), to obtain Note that the dynamics is now captured by two ODEs, where the spatial component appears only as the total cell length variable.These reduced equations have four structural parameters (r, k, κ, f s ) and the parameter describing the strength of the actin treadmilling flow (v).

Results
It is convenient to plot the resulting dynamics on the k −r phase diagram shown in Fig. 2A.We find that above a critical value of the cell stiffness k (to the right of the blue solid line), the spring is too stiff to allow for large length changes that enables the spring to store and release large forces.In this regime we get a single, stable fixed point that corresponds to smooth cell motion, and no stick-slip behavior (for the phase-diagram as function of the adhesion saturation parameter r 0 see Appendix A).
Below the critical stiffness, we find that for small values of r the system corresponds to smooth motion (stable fixed point), of a short and fast-moving cell (Fig. 2H-J).
Above a transition line (solid blue line in Fig. 2A), the fixed point undergoes a Hopf bifurcation, which marks the transition from smooth motion to stick-slip motion (limit cycle dynamics)(Fig.2E-G).In this regime there are large length oscillations, as well as large oscillations in the number of bound adhesion linkers at the cell back (Fig. 3A,B).These oscillations in n occur due to the large values of the force per linker that is reached (Fig. 3C), which leads to catastrophic, avalanche-like detachment events at the back.Within the stick-slip regime, we find that the duration of the limit-cycle increases with increasing r, i.e. the dynamics slow down with increasing substrate adhesiveness (Fig. 3D).
As r increases there is a second Hopf bifurcation, whereby the fixed point is stable again.However, there is a narrow region of bistability, where the phase-space is separated by a separatrix (solid black line in Fig. 3B), such that the stable fixed point coexists with the limit cycle.Within this regime, noise can induce transitions from stick-slip (limit cycle) to smooth motion (stable fixed point), as demonstrated in Fig. 4.
For increasing r the separatrix grows until it meets the limit cycle, which marks the transition to flows that all lead to the single stable fixed point.The dynamics in this regime correspond to a smooth motion of a slow-moving

Comparison to Experiments
We now compare in detail the model's stick-slip limitcycle dynamics (Fig. 2E-G) to experiments (see Appendix D for experimental methods).In Fig. 5 we analyze the stick-slip migration mode of a C6 glioma cell (seeded on laminin-coated lines of 5µm width), by following the cell length and the intensity of membrane-bound vinculin, using total internal reflection fluorescence (TIRF) microscopy.As a component of the adhesion complex, we denote the total intensity of vinculin in the rear of the nucleus as a measure of the adhesion strength at the cell back in the model (n).We find that the experimentally observed dynamics during a stick-slip cycle give rise to a limit-cycle in n, l-phase space that has the same qualitative features as the model predicts (Fig. 2F).
The robust features of the migration modes predicted by the model, can be observed in experiments (Fig. 6).In Fig. 6a we show a kymograph of a motile cell moving along a one-dimensional stripe, exhibiting stick-slip behavior.At a certain point along the trajectory the cell increased both its rate of stick-slip events, and its overall migration velocity.In addition, the amplitude of the stick-slip events decreased in size.All of these features are captured by the model, if the average adhesion strength between the cell and the substrate (r) decreases abruptly along the trajectory, and the cell moves from high to low r within the stick-slip regime (Fig. 2A).The cell-substrate adhesion strength is sensitive to the surface properties along the stripe, which may vary along its length.The backwards growth of the lamellipodia at the cell back, during each stick-slip cycle, are described by the more elaborate model below.
In Fig. 6b we show a cell migrating smoothly, and then abruptly decreasing in length and increasing in speed.This behavior is captured by the model assuming that the adhesion strength r decreases such that the cell jumps from high to low r as shown in (Fig. 2A,D,J).
Note that we can not exclude that the abrupt change in behavior observed for the cells in Fig. 6 originates from a change in some other internal parameter of the cell.In addition, we observe highly dynamic stick-slip behavior at the leading edge of the cell, which we do not describe by our current model [23].

Model description
We now complement the length-adhesion model described above, by incorporating a model for the spontaneous self-polarization of the cell.We use a modified version of the scheme developed in [25] (Fig. 7A), whereby in addition to the advection of proteins that enhance the actin treadmilling, such as myosin-II, we also consider an inhibitor of actin polymerization.The inhibitor protein is free to diffuse in the cytoplasm, and is also advected by the actin treadmilling flow.An example for an inhibitor of local actin polymerization, that is advected by the actin flow, is Arpin [25,28].We furthermore assume that the time-scale of the redistribution of this inhibitor across the cell is much shorter than the timescale of changes in the actin treadmilling speed.This assumption is corrob-orated by recent measurements of myosin-II redistribution time when advected by actin (∼ 10sec) [27], while the stick-slip dynamics of the cell are over timescales of tens of minutes [9].We can therefore use the steady-state distribution of the concentration of this inhibitor c(x), which is given by an exponential function (Fig. 7B) where v is the effective instantaneous actin treadmilling velocity, D the effective diffusion coefficient of the inhibitor and c tot is the total amount of inhibitor molecules in the cytoplasm (see Appendix B).
To complete the model [25], the net treadmilling flow is given by the difference between the actin polymerization flows created at the two ends of the cell, which are inhibited by the local concentration of the inhibitor at the front and rear where β gives the scale of the actin flow in the cell, and the polymerization activity at each end is diminished by the inhibitor concentration at each end, given by a simple Hill function: c(x b,f ) = cs cs+c(x b,f ) (Fig. 7B), where c s is the concentration at which the effect of the inhibitor saturates.
After rescaling (see Appendix C), using Eqs.19,20, we obtain the following implicit equation for the net actin treadmilling flow In the persistent regime (sufficiently large c s ), we find that as l increases above a critical value Eq.( 21) undergoes a pitchfork bifurcation, where spontaneous actin treadmilling appears, corresponding to a polarized and motile cell (Fig. 7C).The critical length is given by which provides a lower bound for the coupling strength β > D/c, below which there is no spontaneous motility.A lower critical length for polarization was also obtained in [29].
The equations of motion are now extended beyond the basic model (Eqs.17,18) to include the dynamics of the actin treadmilling flow where v ss (l) is the steady-state solution to Eq.21, which is the length-dependent actin flow speed, and δ is the rate at which the actin flow adjusts to length changes, i.e. the rate at which v approaches v ss (l).For a cell that has a rest-length that is longer than the critical length for polarization, l c < 1, the resting state of the cell is polarized.Since the velocity saturates to a constant value for l > l c (Fig. 7C), the results are qualitatively similar to those we obtained for a constant, length-independent, v (Fig. 3).
We will therefore focus on the more interesting case where l c > 1, and the rest-length of the cell is smaller than the critical length for polarization.In this case we find another fixed point that corresponds to v = 0 where the cell is at its rest length, and stationary.

Results
In Fig. 8A we plot the k − r phase diagram for the simpler case of δ → ∞, such that the actin flow speed is given directly by the solution v ss (l) (Eq.20,Fig. 7C).In comparison with the case of fixed actin flow (Fig. 2A), there are several new phases, despite the overall similar shape of the phase diagram.
For very low values of r, we find that there is a transition line below which the cell is non-motile (Brown region in Fig. 8A), which corresponds to all the flows in the n − l phase-space leading to the single fixed point at v = 0 (Eq.26).Above this transition line, with increasing r, there are phases of coexistence between smooth motion or stick-slip, and non-motile behavior.We plot in Fig. 9 the dynamics for different points of increasing value of r, for fixed k = 0.8 (vertical dashed line in Fig. 8A).Two flows are demonstrated for each case (green and purple paths and corresponding kymographs in Fig. 9), exposing the coexistence of either smooth motion (E-H, Q-T) or stick-slip (I-L) with the non-motile phase.The bifurcations of the stable and unstable solutions are denoted by their cell length (and the limit cycle amplitude), as well as the limit cycle periods, in Fig. 8B,C (along the same k = 0.8 line).
There is even a thin region of tri-stability (Fig. 8A), which we demonstrate in Fig. 8D-F, by inducing a transition between the three phases by adding noise at specific times.For a fuller exploration of the k − r phase diagram see Appendix E.
Next, we explore the effect of a finite value of δ, where the treadmilling velocity of the actin does not instantaneously adjust to its length-dependent value v s s(l) (Eq.25).
In Fig. 10A we plot the shift in the stick-slip transition line to lower values of k, for decreasing values of δ.The region of no-motility is now pushed to lower values of k.The reason for this is shown in Fig. 10B-D: when the length drops below l c after the slip event, the velocity of the actin does not drop to zero instantaneously (as in the case when δ → ∞).There is therefore a region of the phase-diagram where the length recovers and increases above l c , allowing the stick-slip limit-cycle to survive.The "stall duration", during which l < l c and the cell is almost stalled, increases with increasing value of δ (Fig. 10E).
Part 3: Self Polarized symmetric cell with a dynamic protrusion

Model description
In the model we developed so far, we did not allow for the actin polymerization to produce local traction and protrusive forces at both ends, but rather assumed that the competition between the two ends gives rise to a single leading edge (or to an unpolarized cell) (Fig. 7A).
FIG. 8. A) k − r phase diagram.Brown region correspond to a regime where the cell is stationary.Striped texture region corresponds to a regime where there is no motility due to the unstable limit cycle created by the discontinuous vector field.Blue region correspond to a bi-stable region of polarized stick-slip cells and stationary cells.Orange/Green regions correspond to bi-stable regions of polarized cells with a constant length and stationary cells.Red region corresponds to a region of tri-stability between polarized cells with a constant length, stick slip motility and no motility.Blue/Red/Brown curves are the hopf/saddle node/hopf bifurcation transition lines.Black dashed line is the r cross section for k = 0.8.Black dots correspond tp the points (k, r) = (0.8, 0.4), (0.8, 0.7), (0.8, 1.65), (0.8, 5), (0.8, 8).Red star correspond to the point (k, r) = (0.8, 7.7).B) Maximal and minimal amplitude of the l state variable as a function of r for the limit cycles in the vector field along the cross section of k = 0.8.Dashed blue curve correspond to the unstable discontinuous limit cycle.Solid Black curve corresponds to the stable limit cycle and the cell length.Black dashed curve is the unstable limit cycle.C) The time period of the limit cycles.Blue dashed line indicate the discontinuous unstable limit cycle.Black solid line indicate the stable limit cycle (period of stick slip).Black dashed line indicates the unstable limit cycle in the tri-stable region.D-F) The tri-stable region.Solution of Eqs.23 We now extend the model to describe the local protrusive forces that are produced by actin polymerization at the two ends of the cell (Fig. 11A).
The dynamics of the moving front and back depend on their respective direction of motion: when they protrude there is a drag force that is given by a constant drag coefficient, while when they retract the friction depend on the slip-bond activity of the adhesion (Fig. 11B).We therefore have the following equations of motion for the two ends of the cell where the friction coefficients are given by where Θ is the theta Heaviside function.This way the friction coefficient depends on the local motion of each end.This is illustrated in Fig. 11B, where the friction at the back is due to the slip-bonds when the back end of the cell is retracting ẋb > 0 and Θ ( ẋb ) = 1.Note that the force that pulls on the adhesions during retraction (in the exponential function in Eq.28) is not only given by the elastic spring force (as in the previous models presented above), but is countered by the local force of actin polymerization that is pushing the membrane Similarly, the dynamics of the adhesions are given by actin is now calculated separately on both sides of the cell, given by The distribution of the inhibitor along the cell is affected by the competition between the treadmilling actin from both ends, so its given by Eq.19 with:

Results
The solution of Eq.31 is shown in Fig. 11C.The solution shows that the cell remains unpolarized for l < l c , and the actin tredmilling flows are equal at the two ends of the cell, given by v f /b β l c+l .A polarized cell for l > l c has a higher treadmilling flow from the front, and in the limit of large l the treadmilling velocity approaches these limiting values: v * f = β and v * b = d/c.We find that there a new length scale in the system, which determines the ability of the cell to elongate and reach polarization.For the the unpolarized cell we equate the protrusive forces at both ends to the restoring force of the spring which yields the polarization length By equating the critical length l c (22) to the polarization length l p (33) we obtain a critical value of β, above which the cell polarizes (Fig. 12A,B) (34) As for the previous model, the interesting behavior is for a cell that has a rest-length that is smaller than l c > 1.For such a system we demonstrate the dynamics for different values of β (above and below β c ), in Fig. 12D-L.For β < β c the cell elongates symmetrically (Fig. 12D-F), but does not polarize (no spontaneous symmetry breaking).For β > β c we find regimes of smooth motion (β = 8, Fig. 12G-I) and stick-slip (β = 11, Fig. 12J-L) for parameters that correspond to these phases in the phase diagram of the stationary treadmilling system (Fig. 2).
When the cell is within the stick-slip regime, it may end up after each slip event with a cell length that is shorter than l c .Since in the symmetric model the protrusive force acts locally at both ends, the cell elongates symmetrically until l > l c and it repolarizes.This is demon-strated in Fig. 11B.In Fig. II we show the effect of the parameter c on the stick-slip migration, where increasing c decreases the retrograde flow and increases the time it takes the cell to repolarize following each slip event.
Furthermore, since the cell is effectively stalled when its length is shorter than l c , the time it takes for the cell to expand back to l c after each stick-slip event determines the overall migration velocity.This is demonstrated in Figs.21,22.During each stick-slip event, since the cell loses polarity, it may repolarize in a new direction.This occurs when we add random noise to the equations of motion: the system is solved as a stochastic differential equation where additive noise is added to the actin flow v f and v b .For different amplitudes of the noise term we plot in Fig. 13A the probability p that the cell changes its direction of motion per stick-slip event.We find that for each level of noise there is a critical value of the actin response rate δ (Fig. 30) that sets a threshold between p = 0 and p = 0.5.During each stick-slip event, when l < l c , the actin flow velocities on both sides approach each other exponentially with a time-scale of δ −1 .As v f /b approach within the noise amplitude, they cross each other (Fig. 13B,C), and the side that has the higher flow when l = l c determines the new direction of cell motion.The dynamics around the critical δ where the turning probability increases strongly, are demonstrated in Fig. 14.

COMPARISON TO EXPERIMENTS
When comparing the results of our fully symmetric model to the experimental observations of motile cells along one-dimensional tracks, we need to note that the cellular system is very noisy.Nevertheless, our model exposes distinct motility patterns that should underlie the noisy motion of the cells.
First, we wish to show that the main classes of cell behaviors, which is observed within a uniform population of identical cells (Fig. 15) can be qualitatively spanned by the model if, for example, the actin treadmilling activity (model parameter β) is variable within the population.
Next, we focus on the stick-slip migration pattern, and compare several examples of such observed cell trajectories to the dynamics predicted by the model (Fig. 16).We find that there is strong similarity between the dynamics exhibited by the cell and in the model, regarding the oscillations in the cell length and velocity.The limit-cycle trajectory predicted by the model can be observed in the experimental data.For slightly different set of parameters, we find that the model gives rise to other stick-slip regimes observed in the experiment: a cell with lowfrequency of small-amplitude stick-sip events (Fig. 17), and a cell with high-frequency of large-amplitude stickslip events (Fig. 18).
In Fig. 19 we give a wide range of different cell migrations observed in experiments, demonstrating that the model can produce a similar range of cell migrations.In Fig. 19A we show a cell that seems to be oscillating between being polarized to the right and to the left, such that it hardly migrates over time.In our model a similar behavior may occur if the cell is in the stick-slip regime, and following each slip event it is shorter than the critical length, allowing for noise to induce direction changes (as explained in Figs.13,14).
In Fig. 19B a kymograph is shown of a cell that is initially at rest, and elongates symmetrically.At some later time it breaks the right-left symmetry and moves in a persistent manner.This behavior is well described by our model, where the rest-length of the cell is below the polarization length, but the conditions are such that β > β c and it elongates symmetrically until it reaches a length that is longer than the critical length l c .At a later time point we decreased the value of the surface adhesion strength r, leading to a transition from long-slow to short-fast migration mode, which seems to fit this cellular behavior (similar to Fig. 6).
In Fig. 19C we plot a kymograph of a cell that is persistently moving on a 1D track in the stick-slip regime.During each stick-slip event the progression of the cell's leading edge is slightly modulated, while at the back there are periodic extensions of lamellipodia that undergo large retractions.This behavior corresponds very well to the stick-slip behavior in our model.
Finally, in Fig. 19D we present a kymograph of a cell that is first observed to be migrating by slow stick-slip events.At some point the cell stops and rounds up, as it enters mitosis, from which two daughter cells emerge, both migrating with similar speeds, but one seems to be moving smoothly while the other exhibits some stickslip cycles.We demonstrate a similar chain of migration changes using our model, by modulating the parameters that determine the actin treadmilling activity (β) and adhesion strength (r).These comparisons serve to show that the model can describe the complex migration patterns of cells.The parameters used in the model may correspond to a unique migration mode, or correspond to regimes where different migration modes coexist.The transitions in the cell behavior from one migration mode to another may therefore correspond to the dynamics of the internal parameters of the cell.In addition, noise in the dynamical variables (such as in the actin treadmilling speed), can drive the enhanced polarization changes observed during stick-slip motion.

CONCLUSION
Cell motility involves a large number of cellular components, connected by a complex network of interactions.In addition, the system is noisy, which makes the analysis of cell migration a daunting task.In the present work we developed a theoretical model aimed at exposing the motility patterns of cells moving along a one-dimensional track.One dimensional motion is both simpler to analyze experimentally, and describe theoretically, as well as highly relevant to cells migrating within tissues during development and cancer progression.
We find that the coupling between the cell length and the slip-bonds of adhesion molecules, through the springlike elasticity of the cell, provides the basic mechanism of stick-slip during cell migration.Furthermore, we show that when cells migrate and their polarization is maintained by the overall actin treadmilling flow, the cell length couples strongly to the polarization state of the cell: cells below a critical length do not polarize.While persistent smooth and stick-slip migration modes arise naturally in our model, we predict that stick-slip events may allow for noise-induced direction changes, as the cell recoils to less than the critical length at each stick-slip cycle.
The minimal model that we present recovers the rich variety of cellular migration patterns, which we then compare to experiments on migrating cancer cells (glioma).These comparisons validate the model and demonstrate how it can provide the framework for understanding the complex migration patterns exhibited by migrating cells.
Our model makes detailed predictions regarding the dependence of the cellular migration modes on the coarse-grained parameters of the model, which describe the cell's mechanics (k), actin-polymerization activity (β) and surface adhesion (r).We predict a reentrant smooth migration regime as function of the surface adhesion, with stick-slip occurring only at the intermediate regime (Figs.2,8,10).
Despite the fact that cell motility is very noisy, our models help expose the underlying deterministic patterns of motility [26] that drive the cellular motion, which may be partially masked by noise.We show that the noise can also drive dramatic transitions between the different motility patterns, since the cell may reside close to the transition lines, or in a regime of coexistence, between such modes.These results ofer a new explanation for the large "phenotypic heterogeneity" (CCV) that is observed in cell migration experiments, even under well controlled conditions and monoclonal cell population.These results offer a new framework to explain experimental observations of migrating cells, resulting from noisy switching between underlying deterministic migration modes.Small fluctuations in the cellular components that are greatly amplified due to shifting the cells' internal state between coexisting migration modes, and across phase transition lines.The richness in migration modes that our model predicts for cells with identical or slightly different internal states offers a new paradigm to explain "phenotypic heterogeneity" in cell migration.
Our results explain naturally many observations of cells that grow symmetrically, and migrate along 1D tracks.In addition, the model makes predictions that can motivate future experimental exploration.Finally, having a model for single-cell motility can serve as the basis for the description of collective cell migration [30], which goes beyond treatments of the cells as simplified self-propelled particles.In Eqs.1,2 the parameters α and γ should be multiplied by a factor of r/(r + r 0 ), but throughout the main text (Fig. 2) we treated the limit of r 0 → 0.Here we show the effects of a finite value of r 0 .Fig. 23 demonstrates that as the value of r 0 increases the stick-slip region along the k −r phase diagram decreases, and the cells migrate with a slower speed.Consider a generic polarity cue which follows an advection diffusion transport along the length of the cell l.The polarity cue acts on both ends of the cell, and inhibits actin polymerization.
The derivation starts by dividing the cue into two populations: 1) bound cue proteins c b which are advected by the retrograde flow v, and 2) unbound cue proteins c f which diffuse with a diffusion coefficient D.
The bound/unbound proteins can attach/detach to/from the advected actin, with rates kon / koff respectively.The total concentration of the polarity cue in the cell is considered as conserved, i.e, c tot = c b (x, t)+c f (x, t) ∀ {t, x}.
The advection-diffusion transport equations of c b and c f are given by At the limit of fast exchange, in which the kinetic rates kon , koff By considering mass conservation we obtain B7) and therefore the concentration profile is given by Eq.19.
Appendix C: Rescaling the velocity equation in the asymmetric self-polarized model Next, we rescale Eq.( 20) by the time and length scale k −1 of f and x 0 , such that → v, and Appendix D: Experimental methods 1. Cell culture C6 rat glioma cells were obtained from ATCC and grown in High Glucose DMEM supplemented with 10 % heat inactivated FBS (HI-FBS) and glutamine (Invitrogen).Graded brain tumor specimens were obtained with informed consent, as part of a study protocol approved by the SingHealth Centralised Institutional Review Board A (Singapore).From one of those specimen, derived human glioma propagating cells NNI21 (generous gift from Carol Tang lab, National Neuroscience Institutte Singapore were isolated and cultured as described previously [31] as tumor spheres in high glucose DMEM/F12 (1:1) supplemented with sodium pyruvate, non-essential amino acid, penicillin/streptomycin, glutamine, B27 supplement (Invitrogen), bFGF (20ng/ml), EGF (20ng/ml) (PeproTech), and heparin (5g/ml) (Sigma).For transfection and migration assays, NNI21 were cultured as monolayers on laminin (10 g/ml) coated petri dishes for 3-5 days before transfection.HGPCs transfections were performed with a Neon electroporator (Invitrogen) as per manufacturers recommendations using vinculin-mCherry (gift from P. Kanchanawong, Mechanobiology Institute, National University of Singapore, Singapore) as an adhesion rapporteur and GFP-Plasma Membrane (GFP-PM, Clonetech) for plasma membrane staining.

Micropatterning
5µm lines were printed on glass coverslips using deep UV photopatterning technique as directed in [32].
Briefly, 25-mm glass coverslips were plasma-treated for 5 min and incubated for 1h at room temperature (RT) with poly-l-lysinegraftedpolyethylene glycol (0.1 mg/ml, pLL-PEG, SuSoS) diluted in Hepes [10 mM (pH 7.4)].After washing in water, the pLL-PEGcovered coverslip was placed with the polymer brush facing downward onto the chrome side of a quartz photomask for photolithography treatment (7-min ultraviolet-light exposure).Subsequently, the coverslip was removed from the mask and coated with laminin (10 µg/ml) (Invitrogen) diluted in dPBS for 1h at 37C.Cells were seeded on the patterns and incubated at 37 C. Cell imaging typically started within the following hour.

Microscopy
Phase contrast of live specimens were performed on a Leica AM TIRF MC system equipped with temperature, humidity, and CO2 control.Long term imaging was done using a 10X objective (Leica HCX PL FLU-OTAR 10x/0.30NAPH1 Objective).Acquisitions were typically obtained over a period varying from 2 to 12 h (1 image/30 s).Black dashed line indicate the unstable limit cycles.Black solid lines indicate the stable limit cycles.E-G) Maximal amplitude of the n state variable as a function of r for the limit cycles in the vector field along the cross section ofk = 1.1/0.8/0.5 respectively.Black dashed curves correspond to the unstable limit cycles.Solid Black curves correspond to the stable limit cycles.Dashed blue curve corresponds to the n value of the fixed point in the region where l < x0 (l < 1).Black dot corresponds to the transition between a finite and infinite period of the discontinuous unstable limit cycle (the n cross section in which the unstable limit cycle crosses the section where n = r

FIG. 1 .
FIG. 1.The simplified model.A) Illustration of a the physical model of a cell migrating along a linear track.x b and x f represent the back and front part of the cell of the cell which are connected by a spring with a stiffness k. x0 is the rest length of the spring and l is the length of the cell.v is the velocity of the actin retrograde flow which is assumed to be constant.At the front acts a protrusion force (red arrow) and a drag force (teal arrow).At the back acts a friction force due to the slip bonds (blue arrow).B) The physical model of the stick slip adhesion at the back x b .Stochastic linkers with stiffness κ attach with an average rate of kon and detach with a length dependent rate of kof f (l) (Eq.8) at the back.The linkers stretch on average with a displacement of ∆x (Eq.10).C) The physical model of the protrusion force acting at the front x f : F 1 p and F 2 p (Eqs.1,2) are forces that act on/from the sliding actin filaments, which are balanced by catch adhesions (green short lines).

FIG. 2 .
FIG. 2. A) k − r phase diagram.Blue region correspond to the polarized stick-slip cells.Orange/Green regions corresponds to the polarized cells with a constant length.Red region corresponds a region of bi-stability between polarized cells with a constant length and stick slip.Blue/Red curves are the hopf/saddle node bifurcation transition lines.Black dashed line is the r cross section for k = 0.8.Black dots correspond to the points (k, r) = 0.8, 1, (k, r) = 0.8, 5 and (k, r) = 0.8, 8.4.Red star correspond to the point (k, r) = (0.8, 7.7).B-D) The dynamics at (k, r) = 0.8, 1. B-D) The dynamics at (k, r) = 0.8, 5. H-I) The dynamics at (k, r) = 0.8, 1. Blue/Orange curves at panels B,E,H correspond to the time series of the cell length and adhesion concentration at the rear.Red curve on panels C,F,I represent the trajectory in l, n phase space.Black curves on panels D,G,I are the kymographs where the black curves are cell edges x b and x f .Parameters: fs = 5, κ = 20, v f = 10.

FIG. 3 .
FIG. 3. Stability analysis along the line k = 0.8 (vertical black dashed line in Fig.2A).A) The maximal and minimal length as a function of r.B) The maximal and minimal force applied on a linker at the back part of the cell as a function of r.C) The maximal and minimal adhesion concentration at the rear as a function of r.D) The time period of the limit cycle as a function of r.Solid Black curves indicate the stable limit cycle.Dashed black curves indicate the unstable limit cycle.

FIG. 5 .
FIG. 5. Comparison of the basic stick-slip model (Eqs.17,18) to the stick-slip dynamics observed in experiments (Patientderived Human glioma propagating cell NNI2 transfected with fluorescently taged Vinculin and seeded on laminincoated lines of 5µm width and imaged every 110 sec).A) Snap-shots from the movie, where the cells have fluorescently labeled vinculin (represents the extent of adhesion complexes).The region that defines the cell back is to the left of the vertical dashed line, i.e. behind the nucleus.B) A kymograph of the cell migration, with the stick-slip events marked.C) The dynamics of the cell length and the total amount of vinculin signal at the cell back region, as function of time.The black line on both graphs indicates the stick-slip cycle used to plot the phase-space limit-cycle shown in (D).

FIG. 7 .
FIG. 7. The stick-slip UCSP model.A) Illustration of a the physical model.Unlike the simple model of Fig.1, the treadmilling velocity is now dependent on the cell length v(l), due to the advection of an inhibitory cue.The concentration of the inhibitory cue is represented by the colorbar on the right.The steady-state concentration of inhibitory cue, c(x), is shown in (B, upper panel), which is an exponential function due to a balance between diffusion and advection (Eq.19).The effect of c(x) on the local actin polymerization rate is given by a Hill function (Eq.20) ((B) lower panel): where the inhibitory cue concentration is low the local actin polymerization rate is large.C) The steady state retrograde flow speed as a function of cell length (solution of Eq. 21).Red solid line represent the stable solution.Red dashed line represent the unstable solution.Black dashed line notes the critical length lc of polarization.
FIG.8.A) k − r phase diagram.Brown region correspond to a regime where the cell is stationary.Striped texture region corresponds to a regime where there is no motility due to the unstable limit cycle created by the discontinuous vector field.Blue region correspond to a bi-stable region of polarized stick-slip cells and stationary cells.Orange/Green regions correspond to bi-stable regions of polarized cells with a constant length and stationary cells.Red region corresponds to a region of tri-stability between polarized cells with a constant length, stick slip motility and no motility.Blue/Red/Brown curves are the hopf/saddle node/hopf bifurcation transition lines.Black dashed line is the r cross section for k = 0.8.Black dots correspond tp the points (k, r) = (0.8, 0.4), (0.8, 0.7), (0.8, 1.65), (0.8, 5), (0.8, 8).Red star correspond to the point (k, r) = (0.8, 7.7).B) Maximal and minimal amplitude of the l state variable as a function of r for the limit cycles in the vector field along the cross section of k = 0.8.Dashed blue curve correspond to the unstable discontinuous limit cycle.Solid Black curve corresponds to the stable limit cycle and the cell length.Black dashed curve is the unstable limit cycle.C) The time period of the limit cycles.Blue dashed line indicate the discontinuous unstable limit cycle.Black solid line indicate the stable limit cycle (period of stick slip).Black dashed line indicates the unstable limit cycle in the tri-stable region.D-F) The tri-stable region.Solution of Eqs.23,24 with noise of amplitude (δl, δv) = (0.5, 0.11) injected at t = 3.55 and (δl, δn) = (0.5, 0.05) injected at t = 16.2, to demonstrate the transition between the phases (across the separatrix lines).Blue/Orange/Red curves in D correspond to the dynamics of length/adhesion concentration and actin retrograde flow.Gray dashed lines indicate the time point in which the noise was injected.Green curve in E is the trajectory in the l − n − v phase space.Black solid lines are the separatrices.F displays the kymograph which corresponds to D-E.Parameters: fs = 5,κ = 20,β = 11, c = 3.85, d = 3.85.
FIG. 9.The dynamics along the line of constant k = 0.8 (vertical dashed line in Fig.8A).A-D) r = 0.4.E-H) r = 0.7.I-L) r = 1.65.M-P) r = 5.Q-T) r = 8.Blue/Orange/Red curves in panels A,E,I,M,Q correspond to the time series of the cell length, adhesion concentration and actin retrograde flow respectively (bold/thin lines corresponds to green/purple trajectory).Green and purple lines on panels B,F,J,N,R demonstrate the trajectories in the bi-stable l−n−v phase space regime.Black solid curves are the separatrices.Panels D,G,I display the corresponding kymographs.Parameters: fs = 5, κ = 20, β = 11,c = 3.85,d = 3.85.

FIG. 11 .
FIG. 11.The symmetric model.A) Illustration of the physical model.x b and x f are the back and front part of the cell of the cell which are connected by a spring with a stiffness k. x0 is the rest length of the spring and l is the length of the cell.v f and v b are the steady state velocities of actin retrograde flow from both ends of the cell, which are coupled to the gradient of an inhibitory cue.The concentration of the inhibitory cue is represented by the colorbar on the left.At both ends of the cell acts a protrusion force (red arrow) which is proportional to the retrograde flow speed (red arrow), a drag force (teal arrow) when the edge is moving forward and a friction force due to the slip bonds (blue arrow) when the edge is moving backwards.B) Demonstration of the friction model in Eq.28: Upper panel displays in blue the derivative of the moving rear ẋb .Below in orange is the Heaviside theta function of the time series of ẋb .Right panel displays the corresponding kymograph.C) The steady state retrograde flow speed as a function of cell length (Eq.31).Red solid line represent the stable solution.Red dashed line represent the unstable solution.Black dashed line notes the critical length of polarization lc.Parameters: β = 19,c = 9,d = 40.

FIG. 12 .
FIG. 12.The polarization length.A) The critical length lc (blue) and polarization length lp (orange) as a function of the coupling coefficient β.Gray dashed lines indicate the sample points of β = 4, 8, 11.B) The steady state velocity profile as a function of length for β = 4, 8, 11 (purple, teal, yellow curves).C) The critical coupling coefficient as a function of cell elasticity k.The purple, teal, yellow dashed curves denote the critical values of k for the values of β used in (B,D-L).The black dashed line denotes the critical value of β for k = 0.8 (which is the value of k used in D-L).D-F) Time series, kymograph, and phase space for β = 4. G-I) Time series, kymograph, and phase space for β = 8.J-L) Time series, kymograph, and phase space for β = 11.In panels D,G,I upper blue curve is the time series of the cell length.Dashed black/gray lines are the critical/polarization lengths.Lower red/blue curves represent the time series of the actin retrograde flow speed at the front/back respectively.In panels E,H,K the purple/teal/yellow colors represent β = 4, 8, 11.In panels F,I,L red/blue curves represent the trajectories at the l − n f − v f /l − n b − v b phase spaces (upper/lower part of the surface).Parameters: κ = 20, fs = 5, r = 5, k = 0.8, c = 3.85, d = 3.85.

FIG. 16 .
FIG. 16.Stick-slip migration analysis in experiments and the model.A-C) Experiment: C6 glioma cells seeded on laminincoated lines of 5 µm width (imaged every 30 sec).Kymograph correspond to total time of 3 hours.A) Kymograph of a cell with stick-slip migration.B) Upper panel is the timedependent length of the cell.Lower panel is the cell speed (from the trajectory of the geometric center).C) Phase-space trajectory of the cell.D-F) Model results for a cell in the stick-slip regime: A) kymograph.B) length and velocity time series.C) phase space trajectory.parameters: β = 14, c = 6, d = 20, r = 4, k = 0.8 and δ = 120.

ACKNOWLEDGMENTS
NSG acknowledges that this work is made possible through the historic generosity of the Perlman family.NSG is the incumbent of the Lee and William Abramowitz Professorial Chair of Biophysics and this research was supported by the Israel Science Foundation (Grant No. 1459/17).NCG and PM acknowledge support by IFOM starting package, and the Italian Association for Cancer Research (AIRC), Investigator Grant (IG) 20716 to NCG.

FIG. 20 .
FIG. 20.Time series of the cell length (dark upper panels) and the retrograde flow speed (lower panels) of the symmetric system for different values of the parameter c: A) for c = 3.85.B) for c = 7.7 and C) for C = 11.55.Black dashed line in the upper panels is the critical length of polarization lc.Red/Blue curves in the lower panels represent the actin flow at the front/back.
Appendix B: Concentration profile of the polymerization inhibitor c(x)

D δx 2
, v δx , along any segment δx ∈ [0, L] , (Eqs.B1,B2) reduce to ∂c(where c(x, t) ≡ c b (x, t).At steady state, (Eq.B3) takes the form of 0 The solution of (Eq.B4) readsc(x) = c 0 e − vx D + c 1 (B5)and the coefficients are obtained by applying a no-flux boundary condition at x f and x b and then rescale x 0 x f → x f , x 0 x b → x b and ctot csx0 → c which gives v ss (x f , x b ) = =β change coordinates to l = x f − x b , and obtain Eq.21.

FIG. 24
FIG.24.A) k −r phase diagram as presented in Fig.9.Black dashed lines correspond to the r cross section for k = 0.5, 0.8, 1.1 and k cross sections of r = 0.5, 1.5, 3.5, 5, 8. Red dots correspond to intersections (dynamics at the intersections are shown in Figs.25,26,27).B-C) The time period of the limit cycles in the model along the cross sections of k = 1.1/0.8/0.5 respectively.Black dashed line indicate the unstable limit cycles.Black solid lines indicate the stable limit cycles.E-G) Maximal amplitude of the n state variable as a function of r for the limit cycles in the vector field along the cross section ofk = 1.1/0.8/0.5 respectively.Black dashed curves correspond to the unstable limit cycles.Solid Black curves correspond to the stable limit cycles.Dashed blue curve corresponds to the n value of the fixed point in the region where l < x0 (l < 1).Black dot corresponds to the transition between a finite and infinite period of the discontinuous unstable limit cycle (the n cross section in which the unstable limit cycle crosses the section where n = r 1+r .Parameters: fs = 5,κ = 20,β = 11, c = 3.85, d = 3.85.
FIG.24.A) k −r phase diagram as presented in Fig.9.Black dashed lines correspond to the r cross section for k = 0.5, 0.8, 1.1 and k cross sections of r = 0.5, 1.5, 3.5, 5, 8. Red dots correspond to intersections (dynamics at the intersections are shown in Figs.25,26,27).B-C) The time period of the limit cycles in the model along the cross sections of k = 1.1/0.8/0.5 respectively.Black dashed line indicate the unstable limit cycles.Black solid lines indicate the stable limit cycles.E-G) Maximal amplitude of the n state variable as a function of r for the limit cycles in the vector field along the cross section ofk = 1.1/0.8/0.5 respectively.Black dashed curves correspond to the unstable limit cycles.Solid Black curves correspond to the stable limit cycles.Dashed blue curve corresponds to the n value of the fixed point in the region where l < x0 (l < 1).Black dot corresponds to the transition between a finite and infinite period of the discontinuous unstable limit cycle (the n cross section in which the unstable limit cycle crosses the section where n = r 1+r .Parameters: fs = 5,κ = 20,β = 11, c = 3.85, d = 3.85.