A stochastic modeling framework for single cell migration: coupling contractility and focal adhesions

The interaction of the actin cytoskeleton with cell-substrate adhesions is necessary for cell migration. While the trajectories of motile cells have a stochastic character, investigations of cell motility mechanisms rarely elaborate on the origins of the observed randomness. Here, guided by a few fundamental attributes of cell motility, we construct a minimal stochastic cell migration model from ground-up. The resulting model couples a deterministic actomyosin contractility mechanism with stochastic cell-substrate adhesion kinetics, and yields a well-defined piecewise deterministic process. Numerical simulations reproduce several experimentally observed results, including anomalous diffusion, tactic migration, and contact guidance. This work provides a basis for the development of cell-cell collision and population migration models.

steps: protrusion of the leading edge, assembly of adhesions to the substrate at the front and disassembly in the rear, and contraction of the cell body, thereby producing locomotion [1]. This type of crawling movement requires transmission of contractile forces, generated within the cell cytoskeleton, through substrate adhesions. Such mechanical interactions between various cellular structures is attained by integrating numerous signaling molecules, the most prominent of which are Rho GTPases [66], [69]. While there are other modes of motility relying on, for example, flagellar activity (e.g. spermatozoa or E.coli) or rolling in the bloodstream (e.g. leukocytes), here the focus is on the crawling type of movement.
Due to cell motility being a highly complex and not fully elucidated process, mathematical modeling of the corresponding phenomena is a challenging task. Numerous approaches have been developed, each being able to capture a certain aspect of the process (see [34], [90] for extensive reviews on whole cell motility models and [53] for a review on modeling of its critical components). For example, free boundary and phase-field models of steadily migrating cells in [10], [60] (and its extension in [49]), [65], [71], [75] are able to reproduce cell morphology as a result of mechanical and biochemical interactions. Models of cell migration in [14], [89] explored emergence of various motility modes due to mechanical coupling of intracellular components and the substrate. In their hybrid motility model, Marée et al. [50], [51] explored the mechanochemical interaction in detail by considering the Rho GTPase signaling circuit. A common feature of these continuum models is that they are not able to reproduce experimentally observed random paths of migrating cells. Stochastic motility models of eukaryotic cells have been proposed, for example, in [5], [19], [32], [73], [80]. There stochasticity is driven by a Gaussian process, although there is evidence that the paths of the migrating cells follow a non-Gaussian process [20], [47], [74]. Another way to include randomness is based on velocity jump process, as proposed by [61] in the context of cell motility. Extending the model and including cell-substrate interactions, population migration models have been proposed in [21], [40]. Here, however, our focus is solely on single cell migration.
In this paper and the follow-up works, we strike a middle ground. By proposing a minimal cell representation including a few cellular structures essential for cell motility, we aim at reproducing stochastic migratory paths in various experimental settings. In our model, stochasticity arises as a result of mechanochemical coupling between the cell cytoskeleton and the substrate through adhesive contacts. Specifically, the random events under consideration are the (de)adhesion events of the cell migration cycle, whereas deterministic locomotion and contraction occur between arrival of the events. Here, we do not make any prior assumptions about the distributions that the events and their arrival times follow. Rather, we consider, in detail, the major determinants of adhesion dynamics and derive the complete stochastic description. This is in contrast to most mathematically well-posed models on stochastic cell migration, where it is assumed that the motion follows a Gaussian process. Thereby, the construction of our cell motility model will result in a piecewise deterministic process [16], [17]. Our simulations reproduce experimental observations, such as superdiffusive scaling of the squared displacement [20], [47], [48], biased migration in the presence of an external cue, contact guidance [67], and directed movement due to asymmetric contractility (and in the absence of guidance cues) [84], [87]. Thus, our approach illustrates how detailed treatment of adhesion cluster dynamics can translate into stochastic cell motility description in a mathematically consistent manner. Moreover, due to the minimal character of the model and our bottom-up approach, this work serves as a basis for modeling cell-cell collisions in [81] and population migration, which will be presented in follow-up papers.
This article is organized as follows. In Section 2 we briefly describe the cell motility cycle and the relevant agents. We then introduce the minimal cell representation and describe the deterministic cell motion between the cycle steps. In Section 3 we construct a stochastic model of adhesion events, which signify transitions of cycle stages. In Section 4 we combine the deterministic and stochastic components of the migration cycle to obtain a well-defined piecewise deterministic Markov process. In Section 5 we specify the kinetics of adhesion events. Numerical simulations are performed in Section 6. A discussion of the results and future outlook are presented in Section 7.

The cell motility model
The cell migration cycle begins with protrusion of the leading edge as a result of actin polymerization (Figure 1a). The polymerization process in lamellipodia is mediated by the Arp2/3 complex, which acts downstream of signaling pathways responsible for cell polarization [69]. Next, the protrusions are stabilized due to formation of focal adhesions (FAs) in the lamellae (regions behind the lamellipodia), which link the actin cytoskeleton to the extracellular matrix (ECM). An FA is a multiprotein integrin-based adhesion cluster, which matures in a Rho GTPase dependent manner [66]. Furthermore, FA maturation depends on the applied tension and occurs concomitantly with actomyosin bundle formation [26]. These bundles, called stress fibers (SFs), generate contractile forces due to non-muscle myosin II motors. Due to increased tension at cell rear, FAs rupture. Finally, deadhesion leads to cell body translocation due to cytoskeletal contraction.
In order to construct the mathematical model, we make the following assumptions. First, FA unbinding leads to remodeling of the SF configuration (and of the entire cytoskeleton) and the cell movement, whereas assembly of new FAs leads to restructuring only. Second, FA events need not occur in the order described above. Several adhesions might be assembled (disassembled) before deadhesion (adhesion) occurs. Note also that while the contractile machinery is important, the dynamic instability of adhesions is what drives the migratory process, for stable FAs prevent retraction. Thus, we consider only interactions of SFs and FAs. Moreover, we do not consider the actin polymerization process and simplify the migration cycle to two steps: after FA assembly occurs, a cell does not move, but reconfigures SFs; after disassembly, a cell does both. Neglecting the polymerization process and the reduction to binding/unbinding events can be justified by the fact that one of the major consequences of the leading edge protrusions is promotion of FA assembly. Because the repolarization of migrating cells occurs frequently as an outcome of intricate biochemical activity, then, in order to keep the model tractable, we do not explicitly model cell polarity. Instead, (de)adhesion frequency is indicative of (rear)front.

The cell representation
Consider the situation in Figure 1b. The disk represents a cell. Let the radius be R cell and let the position of the center at time t be x(t) ∈ R 2 . Suppose there are M equally spaced adhesion sites x i (t) ∈ R cell S 1 , i = 1, . . . M on a cell circumference with constant relative distance (see Remark in Section 2.3). Let Y(t) ∈ {0, 1} M be a vector of focal adhesion states at time t, i.e. Y i (t) = 0, 1 correspond to unbound and bound FA at node i, respectively. Since the traction stresses are oriented inward, transmitted to ECM by FAs, and generated by contractile SFs, then the FAs on the circumference must be one of the ends of SFs. Suppose the other end of all SFs at time t is at the position x n (t) ∈ Ω cell := {(x, y) ∈ R 2 | x 2 + y 2 ≤ R 2 cell } (in a cell's reference frame with origin at x), i.e. all SFs are connected at x n . Since stress fibers behave like Hookean springs on extension, but readily buckle under compression [57], then, inspired by Guthardt Torres et al. [33], the force F i 4 at focal adhesion i is given by: where T i is the magnitude of the contractile force due to myosin motors, EA is the onedimensional Young's modulus, L 0 and L c are, respectively, rest and critical lengths, is the unit vector along the i th SF, and δ is a small positive constant. The first case in (2.1) is due to the Hookean behavior of SFs upon extension and myosin tension generation. Furthermore, stress fiber laser ablation experiments [38], [42], [72] revealed that the initial instantaneous response (elastic behavior due to the SF length dependence in the first case) is followed by slower contraction due to myosin activity (force dependence on T i ) in the remaining portion of the fiber. Combined with stress fiber buckling, we obtain the second case in (2.1). Deguchi et al. [18] also found that SF contraction ceased after reaching a certain critical length. This implies that F i = 0 when L i < L c − δ. For technical reasons we assume F i is piecewise continuous -hence the last cases in (2.1). We also assumed for simplicity that myosin generated force T i may vary between SFs, but is otherwise constant.
Note that since x i (t) ∈ R cell S 1 and FA sites are equally spaced, then in polar coordinates we have: and so F i = F i (x n , θ 1 ).
Since the force at x n due to i th SF is −F i , then the net force at x n is then, assuming negligible inertial effects (due to the viscous nature of cytoplasm) and constant Y: where β cell is the drag coefficient in the cytoplasm.
The representation of a cell in such a way is justified for the following reasons: • The traction stresses are largely applied on the cell periphery and their magnitude decays rapidly towards the center [54], [82]. Thus, the cell body SF ends are at or near mechanical equilibrium. Since contractile forces are generated by SFs, then a cell body SF end must be balanced by all other SFs (due to the equilibrium). Hence, it is reasonable to have a single connecting node of radial SFs which is either at mechanical equilibrium (for stationary cells) or tends to it. Figure 2: (a) Magnitude of F i in red. The blue dashed line corresponds to the profile of F i if we were to treat the fiber as a Hookean spring with constant EA/L 0 . (b) Schematic representation of stress fiber contraction. As the fiber contracts below the rest length L 0 , buckling occurs. As myosin mediated contraction causes the fiber to contract below the rest length L 0 , buckling occurs due to lack of resistance to compression. Below the critical length L c , the fiber ceases to contract due to vanishing interfilament distance. Modified from [57].
• Paul et al. [63] demonstrated in their active cable network model, combined with force application originating from nuclear region on FAs by star-like SF arrangement, results in cells acquiring morphologies typical for motile cells. Since the distribution of forces applied on FAs determines which one ruptures, then it also influences the motion of a cell (due to retraction). Since we are primarily interested in cell migration it is justified to assume that this architecture represents a realistic situation. Furthermore, Oakes et al. [59] found that modeling SFs embedded in contractile networks, where only SFs actively contract, yields a behavior mimicking the experimental results -the cytoskeletal flow was directed along the stress fibers. In the same study, the authors concluded that it is appropriate to treat an SF as a 1D viscoelastic contractile element, which also justifies neglecting inertia in Equation (2.2).
• Since motile cells assume a wide variety of cell shapes and continuously remodel their actin cytoskeleton, one can view this representation as a cell shape normalization (it is implicitly assumed that a cell volume remains constant). That is, Figure 1b depicts a cell and forces applied on FAs normalized to a circle. Möhl et al. [54] applied the shape normalization technique to a timelapse series data of migrating keratinocytes and demonstrated that this allows consistent analysis of FA dynamics, actin flow and traction forces. In view of their results, a particular cell traction force map and FA configuration normalized to a circle can be effectively captured by our representation.
Note that x n (t) ∈ Ω cell for t > 0, proved below.
Proof. The proposition above is obviously true, and it follows from the fact that since x i · n < 0 and x n · n > 0.

The cell migration cycle
Recall that during the migration cycle, deadhesion leads to cell body translocation, while adhesion binding does not. In both cases actomyosin contractility leads to reconfiguration of the cytoskeleton. Here we show how our cell representation can describe the reconfiguration and cell body motion following binding and unbinding events. Without loss of generality assume that an event occurred at t = 0. Let τ > 0 be the time of the next adhesion event, be it binding or unbinding. Let Y(0) ∈ {0, 1} M , x(0) ∈ R 2 , and x n (0) ∈ Ω cell be arbitrary. Then, Y(t) = const. for t ∈ [0, τ ).
We assume θ 1 (t = 0) = 0. Since the FA sites are equally spaced, it is sufficient to consider θ 1 (t) only.

Focal adhesion binding
Following an FA binding, we suppose that a cell becomes stationary (i.e. the cell centroid remains constant). However, a newly formed FA and the associated SF lead to cytoskeletal reshaping. Thus, we have the following system of ODEs for t ∈ [0, τ ):

Focal adhesion unbinding
Following an unbinding event, cytoskeletal contraction leads to cell body movement. Due to the circular geometry, the contractile forces induce both rotational and translational motion.
Note that the bound focal adhesions are able to slide for short distances [54]. Oakes et al. [59] found that the cytoskeleton behaves like an elastic solid on timescales up to one hour. Provided the time τ between adhesion events is small enough, the following is justified.θ Figure 3: Force diagram showing transmission of internally generated contractile forces into translational and rotational motion.r andφ are radial and angular unit vectors, respectively.θ 1 is the angular velocity, F is a net contractile force, F r and F ϕ are radial and tangential components of F, x and R cell are cell center and radius, respectively.
The force F along the radial vectorr(x n ) is acting on the cell center, thereby inducing translational motion (see Figure 3). On the other hand, the rotational motion is produced due to F acting along the tangential vectorφ(x n ). The radial and tangential components of the force F are given by: where x n = (x n,1 , x n,2 ) and Note that the characteristic Reynolds number Re is given by where we assumed the surrounding fluid is water (with corresponding values for density ρ and viscosity ν, and that characteristic cell speed s and size L are 0.1 − 1µm/s, L = 10 − 50µm, respectively. Thus, neglecting inertia, we have:ẋ where β ECM and β rot are, respectively, translational and rotational drag coefficients in the ECM (see Appendix A for derivations).

Specification of kinematics
It is convenient to transform the system above into nondimensional form. In order to do so, we define the following scales: • The spatial and cell length scales are defined by cell radii R cell .

8
• The time scale is set by FA disassociation rate k 0 of f , since FA unbinding of leads to locomotion.
• The force scale is defined by the characteristic force F b .
The constants are to be specified later. Whence we define the new variables: and transform F i from (2.1): Note that we have x n ∈ Ω cell := {(x, y) ∈ R 2 | x 2 + y 2 ≤ 1} and x i ∈ S 1 . Let To complete the specification of cell kinematics between adhesion events, we introduce another discrete variable µ(t) ∈ {0, 1}: µ = 1, if the last event was unbinding 0, if the last event was binding.
Then, plugging in (2.4), (2.5) the rescaled quantities and dropping tildes, it follows that the system evolves according to the following ODE system between the FA events: Suppose that just before an event occurs at time t = 0, the cell is in state I. If at time t = 0 (de)adhesion occurs, the cell jumps into state the (II )II and the system evolves according to (2.7) until the next event occurs at time t = τ , after which the cycle begins anew. The scenarios can be characterized as "run" and "tumble" phases in the bottom and top panels, respectively. (b) Schematic representation of the FA positions projected on cell's circumference at t = 0 − and t = τ − in the top and bottom panels, respectively.
Remark. Our assumption on constant relative distance between FA sites stems from two slightly weaker assumptions: 1) total number of adhesion sites (occupied and unoccupied) is constant; 2) there is a neighborhood around each adhesion site, in which no other site is present, and the size of this neighborhood is the same (and constant) for each site. Figure 4 how it reflects on their peripheral motion. This assumption implies that in each line segment of size 2π/M (with M = 8) there is only one FA site present, which may correspond to bound (in red) or unbound FA (in gray). See also Appendix B on our discussion on the number of FA sites.

FA event model
In the previous section we constructed a model of cell motion between FA events. Following [29], here we construct a stochastic model describing the random adhesion/deadhesion events and their arrival times. The discussion here differs from the standard approach of the Gillespie algorithm in [29], as we do not assume that the propensity functions vary inappreciably between the reactions. Moreover, it provides a connection to the theory of PDMPs, as the forms of the objects, necessary to define a piecewise deterministic process (see the next section), follow from the derivations here.

Focal adhesion events
Since there are M FAs and since each FA can participate in only two reactions (binding and unbinding), then there are 2M total possible reactions. We adopt the following convention for enumerating reactions: reaction j corresponds to a binding reaction of the FA site i = (j + 1)/2 if j is odd; otherwise reaction j corresponds to an unbinding reaction of the FA site i = j/2. Let Y(t) be defined as before.
That is, if the FA is (un)bound, the probability of the (un)binding reaction is zero; if the FA is (un)bound, the probability of (binding) unbinding is nonzero. This implies that a j (y, t) = 0 for at least one j ∈ {1, . . . , 2M }. 1 Lemma 3.1. Let Y(t) = y. Then the probability that no FA event occurs in the time Proof. Using the definition of a j , the probability that reaction j does not happen is 1 − a j (y, t)dt. Then, the probability that no FA reaction occurs is: Let K(τ, j|t, y)dτ be the probability, given Y(t) = y and (x(t), x n (t), θ 1 (t)) at time t, that the next reaction will occur in the time interval [t + τ, t + τ + dτ ) and will be reaction j. Here, again, we suppress for clarity the dependence on x, x n , θ 1 .
Proof. Let P (τ |t, y) denote the probability that no reaction occurs in the time interval [t, t + τ ), given y (and x, x n , θ 1 ) at time t. Then, by Lemma 3.1: Letting dτ → 0 we obtain the following ODE: a j (y, t + τ ).
Since P (0|t, y) = 1, the solution P (τ |t, y) is given by: We have then: Let K time (τ |t, y)dτ be the probability that the next reaction will occur in the time interval [t + τ, t + τ + dτ ), given Y(t) = y and (x(t), x n (t), θ 1 (t)) at time t.
Let K index (j|τ, t, y) be the probability that the index of the next reaction is j given Y(t) = y, (x(t), x n (t), θ 1 (t)) at time t and given that the reaction will occur at time t + τ .
Due to equation (3.2), we see that: and a 0 = 0 due to equation (3.1). Obviously, Thus, if T is (random) time until the next reaction, then its probability density function given by K time , its survival function S(s) is given by (without loss of generality, suppose that t = 0): and its (cumulative) distribution function is given by 1 − S(s). 1 Note that the distribution of a random variable is uniquely determined by its distribution function.
Using the proof of Proposition 3.2 one has the following: Let τ > 0 and letK(τ |t, y) be the probability of more than one FA event occurring in the time interval [t + τ, t + τ + dτ ), given the state of the system at time t. ThenK(τ |t, y) = o(dτ ) as dτ → 0.
1 One can check this by differentiating the distribution function, given by 1 − S(s), with respect to s.
Proof. By the proof of Proposition 3.2: since, following the definition of a j , the probability of more than one reaction occurring in Proposition 3.3 implies that we can neglect the case when more than one FA event occurs at the event time. Thus, an FA event (binding or unbinding) unambiguously correspond to a switch in motility state. If this were not the case and the probability of two FA events at the same time were not negligible, then binding and unbinding of distinct FAs could occur simultaneously. Since the cell becomes motile after unbinding only, simultaneous events could lead to ambiguity in determining the motile state of the cell.

Combining the cell motility and the FA event model
With the results of the previous section we can now formally state the cyclical mesenchymal cell motility model as a stochastic process (see Figure 4 and Supplementary Video 1).
• The time T 1 of the FA event is chosen such that P(T 1 > s) = S(s).
• At time t = τ , the index j of the FA event is chosen with probability K index (j|T 1 , 0, y 0 ) Y and µ jump to new values: where e i ∈ R M is the standard basis vector. Note that due to (3.1), we always have Y(t = τ ) ∈ {0, 1} M , since the probability of (un)binding of (un)bound FA is zero.
• The cycle starts anew with initial time t = T 1 and initial values of other variables at this time: starting at t = T 1 we choose the time T 2 of the FA event such that a 0 (y, τ )dτ .
• The system evolves according to (2.7) for t ∈ [T 1 , T 1 + T 2 ) and so on.
One sees that the cyclical process described above is a Markov process, since the evolution of the system depends only on the current state. This completes the formal specification of the model. In the following we will show that this process is well-defined.

Piecewise deterministic process
In this section we briefly overview a class of piecewise deterministic processes, first introduced by Davis [16]. We then show how the deterministic equations, describing motion between stochastic focal adhesion events, can be combined to yield a well-defined piecewise deterministic Markov process (PDMP).

PDMP overview
Let (Ω, F, (F t ) t≥0 , P) be a filtered probability space, where Ω is a sample space, F is a σ-algebra on Ω, (F t ) t≥0 is a (natural) filtration, and P is a probability measure. Let We can define the piecewise deterministic process on the state space (E, E) (for a more detailed general treatment see Davis [17]) by the following objects 1 : I Vector fields (H ν , ν ∈ A) such that for all ν ∈ A there exists a unique global solution X t ∈ Γ to the following equation: Let (ν 0 , X 0 ) ∈ E at time t = 0 be given. Let a survival function S be defined similarly as in equation (3.4): Let T 1 be the first jump time such that and let (ν 1 , X 1 ) be distributed according to the probability law Q(·, φ ν 0 (T 1 , X 0 )). Then, the motion of (ν t , X t ) for t ≤ T 1 is given by: The value of the process at the jump time T 2 is determined by the measure Q(·, φ ν T 1 (T 2 , X T 1 )) and the process continues in a similar way. Thus, we have a well-defined piecewise deterministic process [16]. 16]). The process (ν t , X t ) t≥0 is a homogeneous Markov process.

Cell motility and PDMP
In this section we show that the cyclical cell motility model described in Section 3.2 is a well-defined PDMP.
Then there exists a unique solution of the systeṁ for t > 0.
Proof. Note that since the evolution of x is decoupled from the other two equations, it is sufficient to prove the claim for the following subsystem: Since the right hand side of this system is Lipschitz on D cell and (x n (0), θ 1 (0)) ∈ D cell , then there exists a unique solution of the subsystem (4.4) for time Let A := 1, 2, . . . , 2 M +1 and let α : A → {0, 1} × {0, 1} M be a bijection. This is simply a mapping such that α(ν) = (µ, Y) ∈ {0, 1} × {0, 1} M corresponds to a particular cell motion and FA states (recall that the former can either be moving or stationary). Let where α 1 and α 2 are, respectively, the first and the second components of α 1 . Let the probability (Ω, F, (F t ) t≥0 , P) and state space (E, E) be defined as in the previous section.
We now specify the objects (I,II,III) described in Section 4.1.
I By Proposition 4.2 we see that for all ν ∈ A, there exists a unique global solution to (4.1).
II Note that in our case the rate function a 0 is given by (recalling Section 3.1) 2 : Thus, for the integrability condition to be satisfied, we assume that each probability rate function a j is integrable along the solution of equation (4.1). An exact form of the rates a j satisfying this condition will be given in the subsequent section. Note that a 0 is nonzero, which follows from (3.1).
III In our case, the measure Q(·; (ν, ξ)) is given by (recalling Section 3): where δ is the Kronecker delta function, δ ξ (·) is the Dirac measure at ξ, a + j = a 2j−1 and a − j = a 2j correspond to, respectively, the binding and unbinding probability rates at FA site j, and α 2 (·) i is the i th component of the vector α 2 (·).
The justification for choosing the functions above stems from our deductions in Section 3.1. In particular, the rate function a 0 in (II) is due to (3.3): the probability density function of the jump time T k+1 , given that T k ≤ t < T k+1 , is given by: which corresponds to the survival function given by (4.2).
We now turn our attention to the measure Q in (4.6). The components of X t do not jump, and vary continuously in time, i.e. if T k is the jump time, then X T − k = X T k (see Section 3.2), hence the Dirac measure δ ξ (·) at ξ in (4.6). Clearly, such transition of the continuous component X t of the PDMP at a jump time is probabilistically independent of the transition of the discrete component ν. Hence we have the product of the Dirac measure with the sum, which is a discrete measure for the transition of the discrete component.
By Proposition 3.3, there is only one FA event at a jump time. Hence, for the probability of transition ν → η to be nonzero, the vectors of FA states α 2 (ν) and α 2 (η) must differ only by one component. Consider the following example to illustrate the jump mechanism. Example.
Then, by equation (3.3), the probability of the transition ν → η is given by: Clearly, the transition Y ν → Y η corresponds to the binding event at FA site 2, explaining the Kronecker delta term (see Sections 2.3, 3.2). Now, consider the sum in (4.6) for this example. We see that We therefore obtain Note that if at time t the vector of FA states is given by Y ν , then there are M possible FA state vectors into which a transition can occur with nonzero probability: Similarly as with the rate function a 0 , we can derive equation (4.6) from the principles we established before.
variables before and after the jumps. Then, where we omitted the random variables for notational convenience. Then we have: s., by construction of the process. Since α is a bijection, we have 1 since, by construction of the cell motility process, the new FA state α 2 (η) is determined independently of whether a cell was moving or not (represented by α 1 (ν) ∈ {0, 1}) and the new motility state α 1 (η) is determined only by which FA event took place (binding or unbinding), regardless of whether a cell was previously moving or not. Note that when a jump occurs, then, by Proposition 3.3, only one of j = 1, . . . , 2M possible (binding and unbinding) FA events occurs. Thus, for j, j ∈ {1, . . . , 2M } and j = j the events "reaction j occurs" and "reaction j occurs" are mutually exclusive. We then have, by the definition of conditional probability: is the probability that the FA event j occurs, given α 2 (ν) and ξ.
given α 2 (ν) and ξ, and that the FA event j occurred.

Adhesion kinetics
While we introduced the probability rates of binding and unbinding events, we have not yet fully specified them. Here we elaborate on how such quantities unambiguously correspond to the relevant subcell dynamics by providing the precise functional forms of the propensity functions.

Unbinding rate
Consider the unbinding rate a − j of FA adhesion site j ∈ {1, . . . , M } and let y ∈ {0, 1} M , ξ = (x, x n , θ 1 ) ∈ Γ. Following Bell [7], the bond disassociation rate under applied force is given by: where k 0 of f is the FA disassociation rate under no load, F i is the force applied at the FA, given by equation (2.1), and F b is a characteristic force scale. The last factor y j simply indicates that only bound FAs can unbind (thus satisfying equation (3.1)). Clearly, the function in (5.1) is integrable. Here we neglect the fact that the force may be applied at the FA (and consequently at the transmembrane receptors) at an angle and assume for tractability of the model that it is applied normally to the FA (hence consider only magnitude).
Remark. In the context of cell migration and within the framework of our model, we only consider FA disassembly on the cell periphery (including the lamellae). The primary cause of such FA unbinding has mechanical, rather than biochemical nature due to the cell contractile mechanism applying load to the adhesions. Although it is known that the Rho family of GTPases (in particular its member RhoA) mediates disassembly of FAs, their effect is indirect: the activity of myosin motors, which generate contractile forces in SFs, is regulated by RhoA [68]. Hence the force dependence of the unbinding rate a − j . Recalling the definition of F i in equation (2.1), we can include such indirect biochemical mediation by considering mediators of the force T i . In order to keep the model tractable, we omit the interaction between RhoA and myosin motors.

Binding dynamics
Consider the binding probability rate a + j of the FA adhesion site j ∈ {1, . . . , M } and let y ∈ {0, 1} M , ξ = (x, x n , θ 1 ) ∈ Γ. The probability rate a + j is given by: where k on,j : Γ → R + is the effective binding rate at FA site j. The last term (1−y j ) simply indicates that only unbound FAs can bind. Whereas unbinding can be viewed effectively as a bond rupturing under applied tension, a binding reaction, or focal adhesion assembly and maturation, is a highly regulated process. Due to the complexity of the FA assembly process, we focus on two major mediators of FA formation: Rac activity and contractile forces.

Rac dependence
For simplicity, we assume that the probability of FA formation is directly proportional to local Rac concentration. Consider the case of chemotactic cell migration. Leading edge protrusions preferentially form in the direction of a chemoattractant. Since Rac is required for formation of lamellipodium and FA formation [68], then local Rac activity correlates with the concentration of the chemical cues. Conversely, local Rac activity negatively correlates with chemorepellent. Let Q cue : R 2 → R + denote the concentration of a cue in the spatial domain and let q : R + → R + denote the Q cue dependent concentration of Rac. Clearly, q is an increasing function for the case of an attractant and a decreasing function for a repellent. Then, where we recall that x j is the position of the j th FA site.
For example, we can take Q cue (x) to be the density of the ECM (or chemoattractant) at x ∈ R 2 and take q(x) = x. Then, the probability of binding is simply proportional to the ECM (or chemoattractant) density. 1 Although a more complex function q can be considered, such as those in [30], [35], [58], in order to keep the minimal character of the model, we opted for simple linear relation. Moreover, following Model 4 in [35] and assuming no feedback with phosphoinositides, then in a steady state we have: where R, ρ, C denote the concentrations of Rac, RhoA, Cdc42, and the rest are constant parameters (see [35] for details). That is, there is a linear dependence of Rac on an external cue. Note that the enlargement of nascent adhesions is concurrent with their maturation into focal adhesions. Thus, enlargement and maturation are synonymous. While the initial stage of adhesion growth is force-independent [13], maturation occurs in a force-dependent manner [26], [78], [88]. However, during such a force-dependent maturation, the positive correlation between the adhesion size and the applied tension exists only in the initial stages of maturation. As FAs increase in size, the force dependence plateaus [78].

Force dependence
That is, the study by Stricker et. al. [78] showed that for mature FAs there is no correlation between applied force and FA size. One can consider an adhesion site as mature when its size reaches ∼ 1µm 2 (see e.g. [6], [77]).
Let k 0 on be the force-independent FA maturation rate. Since FAs are larger in size than nascent adhesions, which assemble at a rate of 0.021s −1 , then their assembly is slower and hence we take k 0 on = 0.01s −1 . We now want to find a function that could represent force dependence of FA maturation rate. Denote this function k f orce : R + → R + .
It satisfies the following: • k f orce (0) = k 0 on , i.e. if there is no force applied, the rate is k 0 on .
• If the applied force is below the characteristic force F b , then k f orce is greater than the unbinding rate, i.e. it is more likely that an FA increases in size than that it ruptures.
• If the characteristic force F b is applied, the rate k f orce should equal the unbinding rate, i.e. we assume that there is a dynamic equilibrium in some sense.
• If the applied force is larger than F b , then the unbinding rate dominates the binding rate. Note that as FA increases in size, the force dependence plateaus [78].
Thus, k f orce should plateau around F b . We also assume that for large applied forces k f orce plateaus at k 0 on , since exceeding loads rupture integrin bonds frequently and impede stable maturation.
The following form of k f orce satisfies the conditions above: where F * 1 = F b /4 and F * 2 = 5F b /4 are the midpoints of the sigmoid functions (see Figure  5). To find the values of γ 1 , γ 2 and , see Appendix B.
Remark. Chan and Odde [11] showed that adhesion site dynamics depends on substrate stiffness. In particular, they showed that for a stiff substrate the transmembrane bonds rupture more frequently, compared to the case with softer substrate under the same load, since the critical load is reached faster. This mechanism provides means for a cell to assess the surrounding rheology. Within the context of our model, this means that the force F b is smaller for the stiffer substrate, thus increasing the disassociation rate for the same load (see (5.1)). Consequently, the force dependent binding rate k f orce also changes for the stiffer ECM. In this case, the curves in Figure 5 will shift to the left. Therefore, our model provides an opportunity to include mechanosensitivity of migrating cells by considering the relevant dynamics for individual FAs.
Therefore, the binding propensity rate a + j of an adhesion j ∈ {1, . . . , M } is given by: 2. Non-uniform environment with external cue gradient and uneven myosin motor activity within a cell.
Note that the total number of adhesion sites M is a free parameter, which differs from cell to cell. Nevertheless, it is a crucial parameter, determining whether the motility type is amoeboid or mesenchymal. Amoeboid motility is characterized by a large number of weak adhesions, high turnover, and more contractile cell body. On the other hand, mesenchymal migration relies on fewer, but stronger focal adhesions with slower turnover and lower overall contractility. The cells with the former motility type are faster and have more diffusive motion [48], [62]. Note that if a ± j ∼ O(1), then the rate function is a 0 ∼ O(M ). Therefore, adhesion events occur more frequently for increasing M , implying higher adhesion turnover. Thus, by varying M , our model is capable of explaining this particular aspect of the difference between the two migration types.
For all scenarios we employ the same initial conditions for all cells. Namely, at t = 0 the conditions are: • x is at the origin, x n is uniformly distributed on a circle with radius R cell , and θ = 0.
• Each FA is in (un)bound state and each cell is in moving state with probability 1/2.
We simulate trajectories of n cell := 56 cells for time t end := 600min. We divided the time interval into n time := 5000 intervals, at the end of which we recorded the cell centroid positions x. For details on the parameters and numerical methods used for simulations, see Appendix B and D, respectively. For details on the data analysis, see Appendix C.

Uniform environment
The results of the simulation with uniform spatial cue Q cue are presented in Figure 6. Due to the absence of spatial variation of Q cue , we take q = 1 in equation (5.3).
The cell trajectories with varying number of adhesion sites, depicted in Figure 6 (a-c), show no clear trend and resemble those of a Brownian motion. Indeed, fitting the meansquared displacement msd(t) to the curve msd(t) = β 0 tβ (see Appendix C for details), we see that the exponent isβ ≈ 1 for the three cases (see Figure 6 (d-f)). This suggests that the cell motion has diffusive characteristics in this scenario. In Figure 6 (g-i) we see the simulated distribution of speeds. The average speeds s av and the parameters β 0 ,β are shown in Table 1. We see that as M increases, the cell motion becomes progressively faster and more diffusive 1 , which is expected for a dominantly amoeboid type of motility. Becauseβ ∼ 1, we can estimate the diffusion coefficient D = β 0 /4. The obtained values are lower, but within an order of magnitude estimated by Liu et. al. [48], who found that D ≈ 2.7µm 2 /min. Interestingly, the gamma distribution gives a very good fit to the simulated data for various values of M , suggesting that cell speeds are gamma distributed. Moreover, the obtained values of the average speed s av fall in the range reported by Liu et. al. [48], who found the speeds to be in the range from 0.2µm/min to 7µm/min. Although there are very high speeds observed in Figure 6 (i), which seem to be outliers, speeds as high as 25µm/min have been observed [25]. As expected, rose plots of normalized velocities show no bias in any particular direction in Figure 6 (j-l). Along with time scaling of the squared displacement, persistence of motion can be measured by directionality ratio (distance between cell centroids divided by path length) and velocity autocorrelation [31]. Expectedly, Figure 7 illustrates that the cells strongly deviate from straight-path migration (Figure 7 (left); see also time average of the directionality ratior in Table 1) and the deviation directions are uncorrelated (Figure 7 (right)). The rapid decay in Figure 7 (right) also suggests that correlations in time become stationary very fast. The increased oscillations in Figure 7 (right) towards the end of the observation window are due to decreased number of observations (see Appendix C).  Although our results show in Figure 6 (d-f) the mean-squared displacement scales diffusively (i.e. linearly) with time, this is not consistent with the reported results. For example, Dietrich et. al. [20], Liang et. al. [47], and Liu et. al. [48] showed that the displacement scales superdiffusively. In these studies it was experimentally found that the time scaling went as ∼ tβ, whereβ ≈ 1.2 − 1.5. The primary reason why, in our case, we have diffusive behavior is that, given a certain state of adhesion sites, there is a complete circular symmetry of the probability rates a ± j with respect to x n variable. Due to this symmetry, then, the probability of moving in one direction is exactly the same as the probability of moving in the opposite direction if we rotate x n by π radians. Hence, somewhat reminiscent of a random walk, we obtain a diffusive time scaling of the squared displacement. Moreover, this symmetry of the probability rates is somewhat idealistic, since it implies that the signaling activity relevant for adhesion dynamics is homogeneous within a cell. One of the ways to break this symmetry, is to multiply each binding probability rate a + j by 1 + u, where u ∼ U (−δ, δ) is uniformly distributed on the interval (−δ, δ) with δ ∈ (0, 1). Then, on average, the rates are unmodified 1 . This way, we not only simulate a non-homogeneous binding rate (and hence, for example, nonhomogeneous Rac activity) within a cell, but also simulate otherwise completely identical copies of cells. Such a modification, where we do not explicitly apply a directed, predefined bias can be referred to as chemokinesis [64].   The effect of modifying the rates a + j with δ = 0.05, 0.1, 0.15 can be seen in Figure 8. The cell trajectories, depicted in Figure 8 (a-c), show that the motion consists of periods with relatively regular path intermingled with highly irregular and random movement. In Figure 8 (d-f) we see that the rate modification leads to a superdiffusive time scaling of the mean-squared displacement, as the exponentβ becomes larger than one and falls within the experimentally observed range of values [20], [47], [48]. Moreover, we see that as δ increases, so doesβ, and the increase of the latter is more pronounced for a larger number of adhesion sites M (see also Table 2). This is due to the fact that as each adhesion site is modified independently, the variance of the modified rates of a cell grows with the number of FAs, which corresponds to increased cell polarization, and hence more prominent persistent motion resulting in higher values ofβ. However, the distribution of speeds for the corresponding values of M is virtually identical to the case with the unmodified probability rates (Figure 8 (g-i) and Table 2). The uniform distribution of normalized velocities also remained unchanged (Figure 8 (j-l)). These results suggest that in the absence of spatial cues, the distribution of speeds for a given adhesiveness (represented by the total number of adhesions M ) remains invariant under symmetry breaking of adhesion binding, while the diffusion type (normal vs. anomalous) does not. Thus, the adhesion number and its turnover is a major determinant of the cell speed, which is consistent with [62]. Note that the increased values ofβ indicate that the cells explore a larger surface area [31]. However, other indicators of motion persistence are not affected significantly ( Figure  9), although migration paths become slightly straighter, as indicated by increased values ofr ( Table 2). These results suggest that symmetry breaking of adhesion binding may allow cells to explore larger area without introducing velocity correlations (Figure 9(d-f)).
As cell polarization is required for migration even in the absence of external signals, it is not surprising that our results show that an imbalance of adhesion formation within a cell leads to experimentally observed superdiffusive scaling of the squared displacement [20], [47], [48]. Nevertheless, this highlights a potential mechanism of anomalous diffusion. In the following, we examine whether our model gives biologically consistent results in the 30 case of externally induced polarization.

External cue gradient
We first investigate how cell trajectories are varied in the presence of an external cue gradient. If a cue, for example, is a chemoattractant, then it is well known that adhesion formation in a cell is biased in the direction of the attractant. Thus, to simulate such biased migration, we let the functions Q cue and q to have the following form (recall equation (5.3)): where δ E represents the gradient magnitude and x 2 is the second component of x.
Here, for simplicity we took the identity function for q and a linear cue gradient in the y coordinate. This cue can represent, for example, density of ECM or concentration of a chemoattractant. Thus, we simulate, respectively, hapto-or chemotactic migration.  In the presence of a cue gradient, we see that the cell trajectories, shown in Figure 10 (ac), exhibit a clear trend in the direction of an increasing concentration. The corresponding plots of the mean-squared displacements show the superdiffusive time scaling in Figure 10 (d-f), with the exponentβ > 1 for all cases. Notice that as the number of adhesion sites M increases, so doesβ for the same δ E (see Table 3). Together with the trajectory plots in Figure 10, our results suggest that in the presence of an external gradient, the taxis becomes more prominent and a cell more sensitive to a cue for increasing number of FAs. Moreover, comparing with the case of a uniform environment, we see that although the amoeboid motility is more diffusive in the absence of external cues, it is also more regular and directed when a cue gradient is present (see Tables 1, 2 vs. Table 3 and Figures 6, 8(a-c) (a) vs. 10(a-c)). In Figure 10 (g-i) we see that the evolution of time-averaged exponents β av (t) (see Appendix C) have three phases. Following the rapid increase in the first phase, there is a gradual decrease in the rate of change in the second phase, followed by stabilization of β av (t) atβ. Curiously, a similar behavior has also been observed by Dieterich et al. [20]. The distribution of speeds again remained invariant and the average speeds are very close to the cases with no external cues (see Table 3). However, the frequency of normalized velocities (see Figure 11 (d-f)) show, as expected, that the cell velocities are aligned with the cue gradient. Accordingly, we see that persistent motion emerges: directionality ratio increases compared to unbiased migration (Table 3) and the velocities become correlated  Figure 12(d-f)). We also observe that an external signal has a stronger impact on motion persistence for higher number of adhesions due to relative increases ofr and the degree of velocity autocorrelation. Recall that in the presence of, for example, a chemotactic cue, a cell polarizes so that its adhesion dynamics is aligned with the gradient. In particular, adhesions are preferentially formed at the front (where the chemoattractant concentration is larger), and preferentially ruptured at the back. We can see in Figure 11 (g-i), that our simulation results reproduce such polarized dynamics: the ratio of binding to unbinding events is larger (smaller) than unity in the northern (southern) part of the cells, where the cue is stronger (weaker) relative to the cell centroid. Also, for a smaller number of adhesion sites, the effects of increasing the cue gradient have more noticeable effect on the ratios of events (see 11 (g-i)). This is simply due to the reduced density of adhesion sites, which leads to larger relative difference in the concentration of the cue between them. From Figure 13 we can asses the effect of an external cue Q cue on the the binding rate a + i (omitting the force dependence for clarity), since the rate is proportional to Q cue . Together with Figure 13, the simulations illustrate that directed tactic migration, resulting from biased adhesion formation, follows from the local information about the external cue. That is, the spatial dependence of the FA binding rate is solely due to the local concentration of an external cue (see (5.3)) and no central mechanism for gradient determination was utilized to bias adhesion formation. Consequently, migration along the gradient of an external cue is achieved without its explicit "computation" by the cell.
Along with external cue, force dependence of the binding rate is also important for directed migration and without it, the cells do not exhibit biased migration (data not shown). Figure 14 illustrates how the dependence fits into the migration cycle (recall Figure 1a). For the directed migration to occur, at the time of FA disassociation x n must be preferentially in the rear (Figure 14 step 2). After FA unbinding the increased force at the rear FAs due to extended SFs promotes binding there (Figure 14 step 3). Note that since cell body translocation occurs only after an unbinding event, formation of new FA in the prospective rear of the cell does not lead to backwards movement. Also, due to the external signal more FAs tend to be at the front than at the rear. Thus, the pulling force exerted by the front on the rear tends to be larger than the opposite and hence the cell moves preferentially in the direction of the gradient. Without the signal, of course, movement becomes unbiased, as shown in the previous section. This suggests that the SF length dependence of the forces (see (2.1)) and the force dependence of the FA binding rate (see (5.2)) are necessary for directed migration resulting from biased adhesion formation in the presence of an external signal. Figure 14: The force dependence of the binding rate and the biased adhesion formation during the migration cycle. Side view schematic of the cell is illustrated, where (un)bound FAs are shown as (white)black circles. 1) Initial configuration. 2) Unbinding leads to cell translocation and motion of x n within the cell. 3) Increased force on the cell rear (due to its dependence on SF extension) promotes FA association due to force dependence k f orce of the binding rate (see Section 5.2.2), after which the cycle begins anew.

Fibrillar architecture of ECM
The ECM topography is another important determinant of directed cell migration. In particular, the spatial distribution of the ECM fibers guides the motility by inducing cell shape alignment along the adhesive cues, resulting in characteristic directed movement along the fiber tracts [64]. Such guided migration is called contact guidance [64], [67]. Ramirez-San Juan et. al. [67] showed that contact guidance can be modulated by micrometer scale variations of interfiber spacing. Inspired by this study, we simulate how subcellular scale fiber spacing influences cell motility, and whether such ECM architecture yields migration patterns characteristic of contact guidance.
Similar to the case with an external cue gradient, the functions Q cue and q have the following form: where Ω δ G represents the stripe pattern, δ G = 0.15, 0.25, 0.35 represents the spacing between stripes such that the distances between them is δ G R cell (Figure 15). The stripe width is taken to be 0.25R cell . Similarly as in [67], these dimensions are chosen so that a cell is spread on multiple stripes. The simulation results, shown in Figure 16, indicate that the cell motility has characteristics of contact guidance. Namely, the trajectories show preferential horizontal cell movement (Figure 16 (a-c)), and the displacements are aligned with the fiber pattern (Figure 16 (d-f)). However, increasing the spacing does not simply lead to a greater adhesion alignment along the horizontal direction, as can be observed in Figure 16 (g-i). Rather, it is the combination of the ECM pattern and the radial position of FAs that gives rise to, for example, definite x-shaped adhesion binding patterns (Figure 16 (h)). Such binding (and unbinding) pattern leads to fluctuating movement along northwestsoutheast and northeast-southwest axis, with the resulting net migration pattern shown in Figure 16 (b). Similarly, the binding pattern shown in Figure 16 (i) with more frequent events along the equator corresponds to a mixture of diagonal and horizontal movements (Figure 16 (c)), as larger interfiber spacing precludes FA binding at the poles and facilitates adhesion along, as well as across the stripes in x-shaped pattern (see also Figure 15 (right) for illustration of a characteristic FA configuration). On the other hand, smaller spacing also leads to horizontal movement, but with more frequent vertical displacement across the stripes (Figure 16). These results are also in line with conclusions made by in [41], where it was found that adhesion alignment determines contact guidance We also found that the average speeds were lower than in previous scenarios (Tables 2, 3): 1.52µm/min, 0.94µm/min, and 0.87µm/min corresponding to, respectively, δ G = 0.15, 0.25, 0.35. Interestingly, the average speeds reported in [67] were ∼ 0.6µm/min, although in that study the speeds were nearly constant for varying fiber pattern. In Figure 17 we illustrate the characteristic adhesions pattern and the profiles of the FA binding rate corresponding to ECM architecture in Figure 15 (right). Assuming that there is a mechanical equilibrium for simplicity, we see that the adhesion pattern on the cell's periphery reflect the structure of the environment, since low values of Q cue translate into low probability of focal adhesion binding. Alternatively, if the cell is positioned as in Figure  17 (bottom, left), then the adhesion pattern is modified accordingly. Thus, we see that our assumption about constant relative distance of FAs does not preclude the characteristic cell adhesion patterns to reflect environmental inhomogeneities (see also Figure 13).
Altogether, our simulations of contact guidance are, for the most part, consistent with the observations reported in the literature. In particular, we obtain the expected guidance of cell movement (Figure 16 (a-f)) and the geometric constraint of adhesion sites ( Figure  16 (g-i)) by the fibrillar ECM pattern, in agreement with [64], [67]. Nevertheless, since our model does not explicitly take into account morphological changes in cell shape (recall that in our model cell shape is normalized to a circle; see Section 2) and since cell shape control is essential to contact guidance [64], [67], increasing the interfiber distance does not necessarily lead to greater alignment of cell migration along the ECM fibers in our simulations 1 . Moreover, in the case when the total number of adhesion sites is very low, the stripes are too narrow, and the separation between them is large, then it might occur that all adhesion sites "miss" the stripes, although the cell is spread over multiple stripes. In this case, the probability that any adhesion binds to the substrate is low, which is not biologically consistent. To remedy these shortcomings, the model needs to be extended in order to accommodate strong changes to cell morphology.

Asymmetric contractility
We now investigate how cell motility is influenced by asymmetrical contractile forces in a cell. Along with preferential adhesion formation, due to, for example, a chemotactic gradient, formation of cell rear by increased actomyosin contractile activity serves as an alternative mechanism by which a directed migration can be induced in the absence of such gradient [15]. In particular, local stimulation of contractility leads to directed motility in the direction opposite to the stimulated area, even in the absence of response to chemotactic stimuli [87]. Here we show that our model is also capable of capturing such directed movement, triggered by breaking myosin mediated contractile symmetry.
Recall that T i in equation (2.1) denotes the force generated by myosin motors at an adhesion site i. Instead of taking it constant, we let it vary with the radial position of an FA. Namely, let T i : [0, 2π) → R + be defined as: where δ myo = 0.35, 0.40, 0.45, T i 0 is the constant value used in previous simulations (see Appendix B), and θ + (i − 1) 2π M is the radial position of the i th FA. Thus, contractile forces south of cell equator are larger by 35%, 40%, 45% for corresponding values of δ myo . We should, therefore, expect in our simulations that the northern part of a cell becomes the front due to the imposed contractile symmetry breaking, and that cells will move accordingly (see Figure 18 for illustration).
Indeed, Figure 19 (a-c) shows, as expected, the trajectories of cells maintaining northsouth polarity corresponding to, respectively, front and rear. Since the asymmetry of myosin forces remained during the simulations, the cell's north-south polarity also persisted, resulting in the cell movement that was highly directed along this axis, consistent with [87]. Consequently, we obtain higher values of β av (t), as shown in Figure 19 (d-f). In particular, for δ myo = 0.45, we see that the time scaling of the mean squared displacement is close to ballistic (see Table 4 for values ofβ).   Moreover, increasing the number of FAs leads to more polarized, directed migration. As in the previous cases, neither speed averages (Table 4) nor their distribution (data not shown) changed significantly for a given number of focal adhesions. Interestingly, for δ myo = 0.35, the binding is relatively more frequent in the rear (i.e. south of equator) than in the front, and unbinding is relatively more frequent in the front (i.e. north of equator) than in the back (Figure 19 (g-i)). This suggests, then, that cells were preferentially moving in the southern direction. However, as can be seen in Figure 20, this is not the case. Although movements southwards are more frequent in this situation (due to the above-mentioned event frequencies), the speeds are lower than northward movements: the ratios of the average speeds directed north to the average speeds directed south were found to be 1.0165, 1.0181, and 1.0858 corresponding to, respectively, M = 8, 16, 32. The net effect is northward movement. For higher values of δ myo , we see that the unbinding is, expectedly, more frequent in the rear, while binding is preferentially in the front. These adhesion frequency patterns also illustrate the significance of the force dependence of the FA binding rate. Recalling Figure 5, we see that, for δ myo = 0.4, 0.45 (corresponding to T i = 1.018F b , 1.054F b ), the binding rate dominates unbinding north of equator due to greater SF extension (see Figure 18 for an illustration) leading to increased contractile force. Since the expected adhesion pattern is reversed for δ myo = 0.35 (corresponding to T i = 0.981F b ) and yet the cells migrate northwards, it may suggest that there is a threshold value of δ myo , above which cells can migrate in a certain direction solely by asymmetric contractility, and/or below which cells must additionally bias adhesion formation to do so. This prompted us to investigate whether varying mechanical properties of SFs can yield the expected adhesion pattern for lower degree of asymmetry, corresponding to δ myo = 0.35. Specifically, we varied the buckling length L 0 and stiffness EA such that x = (1 + δ x )x 0 corresponds to the modified value of the parameter x ∈ {L 0 , EA}, where x 0 corresponds to the default value given in Appendix B. In Figure 21(a) we see that reducing the buckling length L 0 by 27% leads to the expected adhesion pattern, while reducing it by 18% leaves it largely unchanged. However, decreasing and increasing stiffness when δ L 0 = −0.18, −0.27, respectively, leads to the opposite results (Figure 21(b,c)). This suggests that if SFs are less prone to buckling and less stiff, lower degree of myosin induced contractile asymmetry may be required to drive directed migration.
Remark. Another way to induce contractile asymmetry is, for example, to decrease the myosin force T i north of the cell's equator. Then, again, the south of the cell equator is more contractile. However, the simulated trajectories show southward directed movement (data not shown), contrary to what we should expect. Therefore, merely inducing contractile asymmetry is not sufficient. For the expected directed migration to occur, there must be a local increase of contractile forces above some critical level in the prospective cell rear, rather than a local decrease of contractility in the prospective front. Interestingly, Yam et al. [87] were able to initiate directed movement by increasing local actomyosin contraction, while locally decreasing the contractile activity did not lead to migration initiation. More recently, Shellard et al. [76] showed that directed collective cell migration of neural crest cells requires greater contractility at the rear of the clump.

Discussion and outlook
In this paper we constructed a stochastic model of cell migration using a minimal representation of cellular structures, essential for crawling, such as stress fibers and focal adhesions. Using this representation, and observing that FA assembly and disassembly events of the migration cycle lead to different migratory outcomes, we obtained the equations describing deterministic cell motion between the random occurrence of FA events. After introducing the rates of FA binding and unbinding, we obtained the remaining necessary objects to define a piecewise deterministic Markov process: the distribution of interarrival times and of the next FA event. Note that the forms of these distributions have been derived, rather than simply postulated. Having specified the coupling between SFs and FAs, as well as between the cellular environment and FAs, we performed numerical simulations. We showed that our model is able to reproduce experimental observations, such as: superdiffusive scaling of the mean-squared displacement [20], [47], [48] (Figure 8); biased motility in the presence of external cue ( Figure 10); contact guidance [67] (Figure 16). In these cases, the obtained results followed solely due to asymmetric, dynamic instability of FAs in direct response to environmental stimuli. Specifically, it is only the biased FA assembly rate that drives biased cell motility along the cue gradient or the fiber tracts ( Figures 11  and 16 (d-i)). That is, preferred velocities were not imposed or chosen in any way, but simply followed as a consequence of front-rear polarity, as the cell front is characterized by preferential FA binding and the rear by unbinding.
Another characteristic of directed migration is the asymmetric contraction of actomyosin bundles. By increasing the force generation of myosin motors in the prospective rear, we obtained directed movement in the opposite direction ( Figure 20). Here asymmetric FA dynamics (and so front-rear polarity) was also obtained, but as a consequence of locally induced contractile activity, consistent with [87].
Our simulation results in various settings suggest that the cell speeds follow a gamma distribution (Figures 6 (g-i), 8 (g-i), 11 (a-c)). Furthermore, the number of adhesion sites seems to be a determinant of the gamma distribution, as its parameters are similar under different settings and given number of FAs. These results suggest that cell speeds are independent of biased FA formation, i.e. the bias only alters the directionality and not the speed. It is also interesting to see a correlation between the number of adhesion sites and diffusivity (Table 1), as well as average speed (Tables 1-4). Note that faster and diffusive amoeboid movement is characterized by an increased number of weaker adhesions with high turnover and contractility [62]. Thus, the aforementioned correlation is also consistent with experimental observations. We note that our model is not fit to take into account motility strongly relying on cell shape control, which is required, for example, in highly mobile cells. However, the simulations reproduce migration along fiber tracts, where cell reshaping takes place [67]. Our results suggest, then, that adhesion along the tracts is sufficient to produce such migration patterns.
Although the model of the internal contractile machinery driving cell locomotion and cytoskeletal remodeling is simple, the resulting numerical simulations explain several experimental observations. Moreover, the cyclical nature of cell motility is captured with our piecewise deterministic model. While migration of a crawling cell is accompanied by changes in its shape, dynamic coupling of cell-substrate adhesions and contractile machinery, i.e. focal adhesions and stress fibers, represent another side of the coin. Numerous sophisticated phase-field or free boundary models that produce realistic morphology of motile cells, often do not emphasize this coupling (and stochasticity) during the migration cycle. We attempted to remedy this issue in our model, and showed with our simulations that numerous aspects of cell migration can be explained without detailed account of cell shape changes. Nevertheless, shape control is essential for a more complete description of the phenomena. We believe this can be done with the framework provided by vertex-based models [24]: a more complex contractility apparatus can be described via active cable network model [33] and a more detailed account of mechanical forces (e.g. protrusions due to actin polymerization) can be done as in [14]. Together with models of RhoGTPases signaling pathways [34], [35], [56], the most significant drawbacks of our approach (including rigid rotation of the SF structure) can be overcome. The presented framework of piecewise deterministic motility process can also be extended to three-dimensional setting, as neither the event interarrival time distribution, nor the transition measure rely on the particular features of migration in a plane.
However, given the relative simplicity of the stochastic model and its ability to explain a handful of the experimental observations, it is possible to extend the model to include cellcell collisions in the context of contact inhibition of locomotion (CIL). Here, the collisions lead to cessation of locomotion and to repolarization, such that the formation of new adhesions at the site of contact is inhibited, while contractility is stimulated [70]. Within the framework of our model this can be implemented in a straightforward manner: collisions cause a switch to a non-moving state and the FA probability rates are modified according to contact location, as was done in our work in [81] Yet another extension is obtained by treating cells as particles and using kinetic theory, yielding equations governing population migration. Thus, we can achieve a multiscale description of cell motility. Due to limitations in size and scope, cell-cell collisions and population migration will be treated in forthcoming works.

Acknowledgement
The author acknowledges support of the German Academic Exchange Service (DAAD).

Appendix A Equations of cell motion
In our model, using common, lab's reference frame will yield the same governing relations, because the involved forces are determined by relative position of cellular structures. Below, we show why this is the case and provide a more detailed explanation regarding the equation of motions for x, x n , θ 1 presented in Sections 2.1-2.2.
Let x n = x + x n and x i = x + x i , where primes indicate the corresponding variables in the lab's reference frame (recall that x i is the position of the i th FA in cell's reference frame). Then, in this frame, the length of the i th SF L i and the unit vector e i along it are given by respectively. Thus, F i (x n , θ 1 ) = F i (x n , θ 1 ), where F i is the force applied by the i th SF at the i th FA. Note that the force at x n (or x n ) due to the i th SF is −F i (x n , θ 1 ) (or −F i (x n , θ 1 )) by action-reaction principle. Therefore, net force F at x n is F (x n , θ 1 , . Neglecting inertia, we have Now, let us examine the equations of motion after FA unbinding, stated in Section 2.2.2, but in lab's reference frame. In this frame, the radial unit vectorr (x n ) from the cell center x is given by (see Figure 3 in the manuscript for illustration) Analogously, the tangential unit vectorφ (x n ) is given bŷ Note that the tangential component F ϕ of the force F at x n induces rotational motion, while the radial component F r of the force F at x n induces translational motion. These components are given by: Neglecting rotational inertia, we then have where the right hand side in the first (second) line is the torque due to tangential component of the force F (F) at x n (x n ). Because of low Reynolds number, we also have due to translational motion induced by the radial component of the force F at x n . In the common reference frame, the following system of ODEs holds after unbinding (using the definition of x n ): which is equivalent to (2.5). Using the common reference frame becomes even less convenient when we formulate and analyze our stochastic process of cell motility. Moreover, our approach in the main text does not contradict the formulation with the single reference frame, and is equivalent to it.

Appendix B Parameter assessment
Note that the length of myosin mini-filaments is ∼ 0.3µm [79] and the interfilament distance is ∼ 1µm in an uncontracted fiber [3]. Assuming vanishing interfilament distance (see Figure 2 for illustration and [57] for a review on actomyosin contraction mechanism), then the proportion of the minifilaments to the initial, uncontracted SF length is 0.3 1+0.3 = 0.23. Kassianidou et. al. [38] showed that the retraction length scales linearly with the initial length. If the interfilament distance vanishes, the front myosin motor cannot perform a power stroke and step forward, which renders the single motor and the entire minifilament unable to apply contractile stress. Therefore, taking the initial length to be ∼ R cell , we estimate the critical length L c = 0.2R cell . Interestingly, Deguchi et. al. [18] found that stress fibers contract, on average, to 20% of their original length. As stress fibers generally span more than half of a cell, and since it was found that there is a preexisting strain [18] in them, we estimate L 0 = 1.1R cell .
If we take R cell = 25µm, and assuming a SF has at rest the length of R cell , we can estimate the number of myosin minifilaments in a fiber to be 25µm/1.3µm ≈ 20. It was estimated in [83] that there are 10-30 myosin motors in each minifilament. As each motor produces a force of 2 − 10pN [23], [55], [57], we then estimate T i = 4nN . Balaban et. al. [6] found that focal adhesions apply a constant stress of ∼ 5.5nN/µm 2 over an area of ∼ 1µm 2 on an elastic substrate. Thus, we take F b = 5.5nN . Assuming a preexisting strain was 0.1 when these measurements were taken, and since T i = 4nN , we then estimate the one dimensional Young's modulus EA = 15nN .
Using Stokes' Law for drag in the low Reynolds number regime, the drag coefficient β ECM can be estimated as: where η ECM is the dynamic viscosity of the environment. Assuming that the viscosity η ECM is higher than that of water, and taking into account that the contact between cell surface and the substrate increase the effective friction, we estimate β ECM ≈ 10 − 10 3 N ·s m 2 × R cell . Similarly, due to the low Reynolds number, the rotational drag coefficient β rot is given by: In order to obtain estimates for the drag coefficient β cell one needs to have an estimate of the cytoplasm viscosity. Assessing the effective cytoplasmic viscosity of migrating cells is a challenging task, since the viscoelastic properties of the cytoskeleton (which, among other things, consists of polymer networks) are highly dynamic due to constant remodeling and spatiotemporal mediation of the rheology by various signaling pathways. Particularly, the actin network bundle size and cross-linkers influence the viscoelastic properties [27]. Furthermore, the effective viscosity experienced by an experimental probe (or a protein) in polymer solutions depends not only on the type of material properties of the fluid, but also on the size of a probe 1 [37]. Inferring that the body being dragged in the cell is the nucleus with radius R nucleus , we estimate: where η cell is the cytoplasm viscosity. Note that a focal adhesion is a cluster of transmembrane receptors (integrins) linking the substrate with the cytoskeleton, which is always under tension. The cluster also includes adapter proteins, which interlinks these receptors with the cytoskeleton (see e.g. reviews [9], [28]). Thus, if there is no load on a focal adhesion, then, since the cytoskeleton is always under tension, such a focal adhesion can be treated as a complex of independent  bonds to the substrate without a link to SF. In the absence of a load, the average cluster lifetime T lif e is given by [22]: where N is the number of bonds in a cluster, H N is the N th harmonic number, k 1 0 and k 1 on are, respectively, unbinding and binding rates of integrins under no load. Note that in the absence of a load, (re-)binding of individual integrins is an independent event, which bears no relation to the FA, since tension is required for an FA to form and sustain itself. We then estimate: k 0 of f = 1/T lif e | k 1 on =0 = k 1 0 /H N .
Li et al. [46] found that k 1 0 = 0.012s −1 for α 5 β 1 -integrin binding to fibronectin. For N = 10 3 − 10 4 we estimate that k 0 of f = 0.05s −1 − 0.1s −1 . The rationale for simulating M = 8, 16, 32 FAs is the following. Note that the number of cell-substrate adhesions is higher than the number of focal adhesions we chose for our simulations. However, not all adhesions are directly involved in translocating the cell body, during which large traction forces are applied to the substrate through focal adhesions (which are fewer in number than immature, less stable focal complexes/points and nascent adhesions). Moreover, detachment of focal adhesions that leads to translocation, is primarily the result of contractile tension applied by ventral stress fibers, as opposed to transverse arcs and dorsal stress fibers [39]. The latter two have primarily structural role, while the former is fundamental to rear retraction [12], [39]. Thus, the number of focal adhesions that are directly involved in cell body translocation is controlled by the number of the associated ventral stress fibers, which are also the most significant source of traction force applied to the substrate due to large tension within them [39], [77]. Although reports of ventral stress fiber numbers are elusive, visual inspection of the fluorescence images in, for example, [36], [43], [77] (or any other appropriate study) suggests that simulations with the chosen number of (ventral) fibers (and focal adhesions) is realistic. Moreover, diameter of focal adhesions d ∼ 1 − 5µm [26]. Assuming that the separation between focal adhesions is comparable to their size, and taking the cell radius to be R cell = 25µm (as in our simulations), we see that the upper range of possible number of adhesions on the cell circumference is 2πR cell /2d ≈ 16 − 80. We reiterate that this number is an estimate of focal adhesions attached to ventral stress fibers, and it underestimates the total number of focal adhesions that a cell employs, since significant number of them are attached to other types of stress fibers and may also be present within the cell body and not on its periphery. We performed simulations with M = 64 focal adhesions and did not find any added insight.
The values of γ 1 , γ 2 and , mentioned in 5.2.2, can be found as follows: Suppose F ≤ F b . Then, Since k f orce (0) = k 0 on and k f orce (F b ) = k 0 of f e, then: It follows that is given as the solution of the following equation: Similarly, since k f orce (F b ) = k 0 of f e, we find: The values of γ 1 , γ 2 , and are fixed for a value of F b .

6:
Set initial condition X s ← List X [last].

10:
Compute X s i := G(X s ; h) and a 0 (ν t , X s i ).