A Stochastic Modelling Framework for Single Cell Migration: Coupling Contractility and Focal Adhesions

Abstract

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, I 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.

1. Introduction

Cell migration is essential for embryogenesis, wound healing, immune surveillance and progression of diseases, such as cancer metastasis. During embryogenesis, coordinated collective migration to particular sites is essential to the development of an organism. Tissue repair and immune response rely on directed movement of cells following external cues, which are produced after the tissue layer is damaged and infected with pathogens. Likewise, cancer cells migrate away from the tumour into surrounding tissue and distant organs to form metastases, which is the leading cause of death among cancer patients.

Cell motility is a cyclical process, involving morphological changes of the cell body and adhesive contacts with the underlying substrate [1]. The cycle can be divided into three 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 [2]. 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 signalling molecules, the most prominent of which are Rho GTPases [3,4]. 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 modelling 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 [5,6] for extensive reviews on whole cell motility models and [7] for a review on modelling of its critical components). For example, free boundary and phase-field models of steadily migrating cells in [8,9] (and its extension in [10]), [11,12,13] are able to reproduce cell morphology as a result of mechanical and biochemical interactions. Models of cell migration in [14,15] 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. [16,17] 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 [18,19,20,21,22]. There stochasticity is driven by a Gaussian process, although there is evidence that the paths of the migrating cells follow a non-Gaussian process [23,24,25]. Another way to include randomness is based on velocity jump process, as proposed by [26] in the context of cell motility. Extending the model and including cell–substrate interactions, population migration models have been proposed in [27,28]. Here, however, our focus is solely on single cell migration.

In this paper and the follow-up works, I strike a middle ground. By proposing a minimal cell representation including a few cellular structures essential for cell motility, I 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, I do not make any prior assumptions about the distributions that the events and their arrival times follow. Rather, I 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 Markov process (PDMP) [29,30]. Our simulations reproduce experimental observations, such as superdiffusive scaling of the squared displacement [23,24,31], biased migration in the presence of an external cue, contact guidance [32] and directed movement due to asymmetric contractility (and in the absence of guidance cues) [33,34]. Thus, our approach illustrates how detailed treatment of adhesion cluster dynamics can translate into stochastic cell motility description in a mathematically consistent manner.

Due to the minimal character of the model and our bottom-up approach, this work serves as a basis for modelling cell–cell collisions in [35] and population migration (Chapter 7 in [36]). In particular, the latter, mesoscopic scale model can be derived from our microscopic model presented here (Chapter 7 in [36]). The derived model describes population migration via a (probability) density function, and captures the aforementioned intracellular activity. Therefore, the work presented here falls, in part, into a broader framework called Kinetic Theory of Active Particles (KTAP) [37]. Several multiscale approaches utilized this framework to model chemo-, haptotaxis and movement along fibre tracts [27,28,38]. Although in these works subcell dynamics has been considered via a description of receptor binding, interaction of subcellular structures, essential for generation of cell migration, has been neglected. Here, I describe such interactions in detail and show that it can be wrapped into PDMP theory, which can then be used to derive kinetic description more rigorously. As the derivation is out of scope, I leave its demonstration for the forthcoming work.

This article is organized as follows. In Section 2, I briefly describe the cell motility cycle and the relevant agents. I then introduce the minimal cell representation and describe the deterministic cell motion between the cycle steps. In Section 3, I construct a stochastic model of adhesion events, which signify transitions of cycle stages. In Section 4 I combine the deterministic and stochastic components of the migration cycle to obtain a well-defined piecewise deterministic Markov process.At the end of this section, I briefly remark how on relation of the process to the KTAP. In Section 5, I 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.

2. 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 signalling pathways responsible for cell polarization [4]. 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 [3]. Furthermore, FA maturation depends on the applied tension and occurs concomitantly with actomyosin bundle formation [39]. These bundles, called stress fibres (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, I make the following assumptions. First, FA unbinding leads to reconfiguration of SFs and the cell movement, whereas assembly of new FAs leads to only to the former. 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, I consider only interactions of SFs and FAs. Moreover, I 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. As a result that the repolarization of migrating cells occurs frequently as an outcome of intricate biochemical activity, then, in order to keep the model tractable, I do not explicitly model cell polarity. Instead, (de)adhesion frequency is indicative of (rear)front. By neglecting the polymerization process, dynamic alterations to actin cytoskeleton (i.e., cytoskeletal remodelling) is simplified to SF reconfiguration.

2.1. 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 centre at time t be $\mathbf{x}\left(t\right)\in {\mathbb{R}}^2$. Suppose there are M equally spaced adhesion sites ${\mathbf{x}}_i\left(t\right)\in R_{cell}{\mathbb{S}}^1$, $i=1,\dots M$ on a cell circumference with constant relative distance (see Remark 1 in Section 2.3). Let $\mathbf{Y}\left(t\right)\in {\left\{0,1\right\}}^M$ be a vector of focal adhesion states at time t, i.e., $Y_i\left(t\right)=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 ${\mathbf{x}}_n\left(t\right)\in {\mathsf{\Omega }}_{cell}:=\left\{(x,y)\in {\mathbb{R}}^2\phantom{.}|\phantom{.}{x}^2+y^2\le R_{cell}^2\right\}$ (in a cell’s reference frame with origin at $\mathbf{x}$), i.e., all SFs are connected at ${\mathbf{x}}_n$. Since stress fibres behave like Hookean springs on extension, but readily buckle under compression (see Figure 2b for illustration) [40], then, inspired by Guthardt Torres et al. [41], the force ${\mathbf{F}}_i$ at focal adhesion i is given by:

$$\begin{array}{c}{\mathbf{F}}_i=\left\{\begin{array}{c}\left(T_i+EA\frac{L_i-L_0}{L_0}\right){\mathbf{e}}_i,\phantom{aabc}{L}_0<L_i\hfill \\ T_i{\mathbf{e}}_i,\phantom{asdasdasdaabcsaaa}{L}_c\le L_i\le L_0\hfill \\ \frac{L_i-L_c+\delta }{\delta }{T}_i{\mathbf{e}}_i,\phantom{Asdaasssda}{L}_c-\delta \le L_i<L_c\hfill \\ 0\phantom{sdaasdAADSADAdA}{L}_i<L_c-\delta ,\hfill \end{array}\right.\hfill \end{array}$$

(1)

where $T_i$ is the magnitude of the contractile force due to myosin motors, $EA$ is the one-dimensional Young’s modulus, $L_0$ and $L_c$ are, respectively, rest and critical lengths, $L_i=\parallel {\mathbf{x}}_n-{\mathbf{x}}_i\parallel$, ${\mathbf{e}}_i=\frac{{\mathbf{x}}_n-{\mathbf{x}}_i}{L_i}$ is the unit vector along the $i^{\mathrm{th}}$ SF and $\delta$ is a small positive constant. The first case in Equation (1) is due to the Hookean behaviour of SFs upon extension and myosin tension generation. Furthermore, stress fibre laser ablation experiments [42,43,44] revealed that the initial instantaneous response (elastic behaviour 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 fibre. Combined with stress fibre buckling, the second case in Equation (1) is obtained. Deguchi et al. [45] also found that SF contraction ceased after reaching a certain critical length. This implies that $F_i:={\mathbf{F}}_i·{\mathbf{e}}_i$, the magnitude of ${\mathbf{F}}_i$ along ${\mathbf{e}}_i$, is zero when $L_i<L_c-\delta$. For technical reasons, I assume $F_i$ is piecewise continuous (see Figure 2a for profile of $F_i$)—hence the last cases in Equation (1). For simplicity, it is also assumed that myosin generated force $T_i$ may vary between SFs, but is otherwise constant.

Note that since ${\mathbf{x}}_i\left(t\right)\in R_{cell}{\mathbb{S}}^1$ and FA sites are equally spaced, then in polar coordinates we have:

$$\begin{array}{ccc}{\mathbf{x}}_i\left(t\right)=R_{cell}{(cos\left({\theta }_i\right),sin\left({\theta }_i\right))}^T,\hfill & \hfill & {\theta }_i\left(t\right):={\theta }_1\left(t\right)+(i-1)\frac{2\pi }{M},\hfill \end{array}$$

and so ${\mathbf{F}}_i={\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1)$.

Since the force at ${\mathbf{x}}_n$ due to $i^{\mathrm{th}}$ SF is $-{\mathbf{F}}_i$, then the net force at ${\mathbf{x}}_n$ is

$$\begin{array}{c}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y}):=-\sum _{i=1}^M{Y}_i{\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1),\hfill \end{array}$$

(2)

then, assuming negligible inertial effects (due to the viscous nature of cytoplasm) and constant $\mathbf{Y}$:

$$\begin{array}{c}{\beta }_{cell}{\dot{\mathbf{x}}}_n=\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y}),\hfill \end{array}$$

(3)

where ${\beta }_{cell}$ is the drag coefficient in the cytoplasm (see Appendix B for details).

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 centre [46,47]. 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.
• Paul et al. [48] demonstrated that application of force, originating from nuclear region, on FAs by star-like SF arrangement, results in cells acquiring morphologies typical for motile cells. Since the forces applied on FAs by SFs determine which one ruptures, then it also influences the motion of a cell (due to retraction). Since I am primarily interested in cell migration, it is justified to assume that this architecture represents a realistic situation. Furthermore, Oakes et al. [49] found that modelling SFs embedded in contractile networks, where only SFs actively contract, yields a behaviour mimicking their experimental results—the cytoskeletal flow was directed along the stress fibres. 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).
• 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 normalized to a circle. Möhl et al. [46] 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.

Without loss of generality, I assume that our system starts to evolve at $t=0$ and ${\mathbf{x}}_n\left(0\right)\in {\mathsf{\Omega }}_{cell}$. Note that if ${\mathbf{x}}_n\left(0\right)\in {\mathsf{\Omega }}_{cell}$, then ${\mathbf{x}}_n\left(t\right)\in {\mathsf{\Omega }}_{cell}$ for $t>0$, proved below. This result shows that our representation of the cell does not yield unphysical outcome, where stress fibres move outside of the cell. That is, provided initial conditions at $t=0$ are physically well-defined (stress fibres are contained in the cell, i.e., ${\mathbf{x}}_n\left(0\right)\in {\mathsf{\Omega }}_{cell}$), then they will remain so.

Proposition 1.

Let ${\mathbf{x}}_n\in \partial {\mathsf{\Omega }}_{cell}$, ${\theta }_1\in [0,2\pi )$, and $\mathbf{Y}\in {\left\{0,1\right\}}^M$ be arbitrary. Let $\mathbf{n}$ be outward unit normal at ${\mathbf{x}}_n$. Then $\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\mathbf{n}\le 0$ with equality sign if and only if $\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})=0$.

Proof. 

The proposition above is obviously true, and it follows from the fact that

$$\begin{array}{c}{Y}_i(-{\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1))·\mathbf{n}=Y_i\parallel {\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1))\parallel (-{\mathbf{e}}_i·\mathbf{n})=\frac{Y_i\parallel {\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1))\parallel }{L_i}\left({\mathbf{x}}_i-{\mathbf{x}}_n\right)·\mathbf{n}\le 0,\hfill \end{array}$$

since ${\mathbf{x}}_i·\mathbf{n}<0$ and ${\mathbf{x}}_n·\mathbf{n}>0$. □

Corollary 1.

Let ${\mathbf{x}}_n\left(0\right)\in {\mathsf{\Omega }}_{cell}$ be arbitrary and let ${\theta }_1\left(t\right)$ be given. Suppose ${\mathbf{x}}_n\in C^1\left([0,\infty )\right)$ is a solution of Equation (3). Then, ${\mathbf{x}}_n\left(t\right)\in {\mathsf{\Omega }}_{cell}$$\forall t>0$ and $\forall \mathbf{Y}\left(t\right)\in {\left\{0,1\right\}}^M$.

Proof. 

Due to Equation (3), it suffices to show that $\forall {\mathbf{x}}_n\in \partial {\mathsf{\Omega }}_{cell}$ we have $\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\mathbf{n}\le 0$, which follows from Proposition 1. □

2.2. 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, I 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 $\tau >0$ be the time of the next adhesion event, be it binding or unbinding. Let $\mathbf{Y}\left(0\right)\in {\left\{0,1\right\}}^M$, $\mathbf{x}\left(0\right)\in {\mathbb{R}}^2$, and ${\mathbf{x}}_n\left(0\right)\in {\mathsf{\Omega }}_{cell}$ be arbitrary. Then, $\mathbf{Y}\left(t\right)=const.$ for $t\in [0,\tau )$.

I assume ${\theta }_1(t=0)=0$. Since the FA sites are equally spaced, it is sufficient to consider ${\theta }_1\left(t\right)$ only.

2.2.1. Focal Adhesion Binding

Following an FA binding, I 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\in [0,\tau )$:

$$\begin{array}{cc}\dot{\mathbf{x}}\hfill & =0\hfill \\ {\dot{\mathbf{x}}}_n\hfill & ={\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})\hfill \\ {\dot{\theta }}_1\hfill & =0.\hfill \end{array}$$

(4)

2.2.2. 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 [46]. Oakes et al. [49] found that the cytoskeleton behaves like an elastic solid on timescales up to one hour. Provided the time $\tau$ between adhesion events is small enough, the following is justified.

The force $\mathbf{F}$ along the radial vector $\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)$ is acting on the cell centre, thereby inducing translational motion (see Figure 3). On the other hand, the rotational motion is produced due to $\mathbf{F}$ acting along the tangential vector $\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)$. The radial and tangential components of the force $\mathbf{F}$ are given by:

$$\begin{array}{cc}{F}_r\hfill & :=\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\hfill \\ F_{\mathit{\phi }}\hfill & :=\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right),\hfill \end{array}$$

where ${\mathbf{x}}_n=(x_{n,1},x_{n,2})$ and

$$\begin{array}{ccc}\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)=\frac{{\mathbf{x}}_n}{\parallel {\mathbf{x}}_n\parallel },\hfill & \hfill & \widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)={\left(-\frac{x_{n,2}}{\parallel {\mathbf{x}}_n\parallel },\frac{x_{n,1}}{\parallel {\mathbf{x}}_n\parallel }\right)}^T.\hfill \end{array}$$

Note that the characteristic Reynolds number $Re$ is given by

$$\begin{array}{c}Re=\frac{\rho ·s·L}{\nu }\sim {10}^{-6}-{10}^{-4},\hfill \end{array}$$

Where I assumed the surrounding fluid is water (with corresponding values for density $\rho$ and viscosity $\nu$, and that characteristic cell speed s and size L are 0.1–1 $\mathsf{\mu }$m/s, L = 10–50 $\mathsf{\mu }$m, respectively. Thus, neglecting inertia, we have:

$$\begin{array}{cc}\dot{\mathbf{x}}\hfill & ={\beta }_{ECM}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\hfill \\ {\dot{\mathbf{x}}}_n\hfill & ={\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})\hfill \\ {\dot{\theta }}_1\hfill & ={\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right),\hfill \end{array}$$

(5)

where ${\beta }_{ECM}$ and ${\beta }_{rot}$ are, respectively, translational and rotational drag coefficients in the ECM (see Appendix A for derivations). Note that drag depends on the body shape, which changes frequently in migrating cells. However, for simplicity, the coefficients are assumed to be constant.

2.3. Specification of Kinematics

It is convenient to transform the system above into nondimensional form. In order to do so, I define the following scales, characterized by their respective constants:

• The spatial and cell length scales are defined by cell radii $R_{cell}$.
• The time scale is set by FA disassociation rate $k_{off}^0$, since FA unbinding of leads to locomotion.
• The force scale is defined by the characteristic force $F_b$ usually observed at an FA.

The values of these constants are given in Appendix B. Whence the new variables are defined:

$$\begin{array}{cccccc}\tilde{\mathbf{x}}\hfill & :=\frac{\mathbf{x}}{R_{cell}},\hfill & {\tilde{\mathbf{x}}}_n\hfill & :=\frac{{\mathbf{x}}_n}{R_{cell}},\hfill & \tilde{t}\hfill & :=k_{off}^{0}t\hfill \end{array}$$

and transform ${\mathbf{F}}_i$ from Equation (1):

$$\begin{array}{c}{\tilde{\mathbf{F}}}_i:=\frac{{\mathbf{F}}_i}{F_b}=\left\{\begin{array}{c}\left({\tilde{T}}_i+\widetilde{EA}\frac{{\tilde{L}}_i-{\tilde{L}}_0}{{\tilde{L}}_0}\right){\tilde{\mathbf{e}}}_i,\phantom{abc}{\tilde{L}}_0<{\tilde{L}}_i\hfill \\ {\tilde{T}}_i{\tilde{\mathbf{e}}}_i,\phantom{asd1sdasdaabcsaa}{\tilde{L}}_c\le {\tilde{L}}_i\le {\tilde{L}}_0\hfill \\ \frac{{\tilde{L}}_i-{\tilde{L}}_c+\tilde{\delta }}{\tilde{\delta }}{\tilde{T}}_i{\tilde{\mathbf{e}}}_i,\phantom{sdaasdssda}{\tilde{L}}_c-\tilde{\delta }\le {\tilde{L}}_i<{\tilde{L}}_c\hfill \\ 0\phantom{sdaasdAADSADAdA}{\tilde{L}}_i<{\tilde{L}}_c-\tilde{\delta },\hfill \end{array}\right.\hfill \end{array}$$

(6)

where

$$\begin{array}{cccccccc}{\tilde{L}}_i\hfill & =\frac{L_i}{R_{cell}},\hfill & {\tilde{L}}_0\hfill & =\frac{L_0}{R_{cell}},\hfill & {\tilde{L}}_c\hfill & =\frac{L_0}{R_{cell}},\hfill & \tilde{\delta }\hfill & =\frac{\delta }{R_{cell}},\hfill \\ {\tilde{T}}_i\hfill & =\frac{T_i}{F_b},\hfill & \widetilde{EA}\hfill & =\frac{EA}{F_b}\hfill & {\tilde{\mathbf{e}}}_i\hfill & =\frac{{\tilde{\mathbf{x}}}_n-{\tilde{\mathbf{x}}}_i}{{\tilde{L}}_i},\hfill & {\tilde{\mathbf{x}}}_i\hfill & =\frac{{\mathbf{x}}_i}{R_{cell}}.\hfill \end{array}$$

Note that we have ${\tilde{\mathbf{x}}}_n\in {\tilde{\mathsf{\Omega }}}_{cell}:=\left\{(x,y)\in {\mathbb{R}}^2\phantom{,}|\phantom{,}{x}^2+y^2\le 1\right\}$ and ${\tilde{\mathbf{x}}}_i\in {\mathbb{S}}^1$.

Let

$$\begin{array}{cccccccc}\tilde{\mathbf{F}}\hfill & :=\mathbf{F}/F_b,\hfill & {\tilde{\beta }}_{cell}\hfill & :=\frac{k_{off}^0{R}_{cell}}{F_b}{\beta }_{cell},\hfill & {\tilde{\beta }}_{ECM}\hfill & :=\frac{k_{off}^0{R}_{cell}}{F_b}{\beta }_{ECM},\hfill & {\tilde{\beta }}_{rot}\hfill & :=\frac{k_{off}^0}{R_{cell}{F}_b}{\beta }_{rot}.\hfill \end{array}$$

To complete the specification of cell kinematics between adhesion events, another discrete variable is introduced $\mu \left(t\right)\in \left\{0,1\right\}$:

$$\mu =\{\begin{array}{c}1,\mathrm{if}\mathrm{the}\mathrm{last}\mathrm{event}\mathrm{was}\mathrm{unbinding}\hfill \\ 0,\mathrm{if}\mathrm{the}\mathrm{last}\mathrm{event}\mathrm{was}\mathrm{binding}.\hfill \end{array}$$

Then, plugging in Equations (4) and (5) the rescaled quantities and dropping tildes, it follows that the system evolves according to the following ODE system between the FA events:

$$\begin{array}{cc}\dot{\mathbf{x}}\hfill & =\mu {\beta }_{ECM}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\hfill \\ {\dot{\mathbf{x}}}_n\hfill & ={\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})\hfill \\ {\dot{\theta }}_1\hfill & =\mu {\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)\hfill \end{array}$$

(7)

Remark 1.

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 neighbourhood around each adhesion site, in which no other site is present, and the size of this neighbourhood 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\pi /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.

3. FA Event Model

In the previous section a model of cell motion between FA events has been constructed. Following [50], here I 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 [50], as I 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.

3.1. 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. The following convention for enumerating reactions is adopted: 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 $\mathbf{Y}\left(t\right)$ be defined as before.

Let $a_j(\mathbf{y},t)dt$ be the probability, given $\mathbf{Y}\left(t\right)=\mathbf{y}\in {\left\{0,1\right\}}^M$, $\mu \left(t\right)$, $(\mathbf{x}\left(t\right),{\mathbf{x}}_n\left(t\right),{\theta }_1\left(t\right))$ and time t, that reaction j will occur in the time interval $[t,t+dt)$. For clarity, the dependence of the rate $a_j(\mathbf{y},·)$ on $(\mathbf{x}(·),{\mathbf{x}}_n(·),{\theta }_1(·))$ and $\mu (·)$ is suppressed. I assume that the rate $a_j$ satisfies the following:

$$a_j(\mathbf{y},t)=\{\begin{array}{c}0,\mathrm{if}j\mathrm{is}\mathrm{odd}\mathrm{and}{y}_{(j+1)/2}=1\hfill \\ 0,\mathrm{if}j\mathrm{is}\mathrm{even}\mathrm{and}{y}_{j/2}=0\hfill \\ \ne 0,\mathrm{else}.\hfill \end{array}$$

(8)

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(\mathbf{y},t)\ne 0$ for at least one $j\in \left\{1,\dots ,2M\right\}$ (for each FA site $i\in \left\{1,\dots ,M\right\}$, either $a_{2i-1}$ or $a_{2i}$ is nonzero). For consistency of mathematical theory, we must consider the entirety of state space, including the case of FA (un)binding when it is already (un)bound. However, due to Equation (8), we effectively rule out such scenario.

Lemma 1.

Let $\mathbf{Y}\left(t\right)=\mathbf{y}$. Then the probability that no FA event occurs in the time interval $[t,t+dt)$ is $1-{\sum }_{j=1}^{2M}{a}_j(\mathbf{y},t)dt+o\left(dt\right)$.

Proof. 

Using the definition of $a_j$, the probability that reaction j does not happen is $1-a_j(\mathbf{y},t)dt$. Then,

$$\begin{array}{c}\prod _{j=1}^{2M}\left(1-a_j(\mathbf{y},t)dt\right)=1-\sum _{j=1}^{2M}{a}_j(\mathbf{y},t)dt+o\left(dt\right),\hfill \end{array}$$

is the probability that no FA reaction occurs is. □

Let $K(\tau ,j|t,\mathbf{y})d\tau$ be the probability, given $\mathbf{Y}\left(t\right)=\mathbf{y}$ and $(\mathbf{x}\left(t\right),{\mathbf{x}}_n\left(t\right),{\theta }_1\left(t\right))$ at time t, that the next reaction will occur in the time interval $[t+\tau ,t+\tau +d\tau )$ and will be reaction j. Here, again, the dependence on $\mathbf{x}$, ${\mathbf{x}}_n$, ${\theta }_1$ is suppressed for clarity.

Proposition 2.

Let $\tau >0$ and $\mathbf{Y}\left(t\right)=\mathbf{y}$. Then,

$$\begin{array}{c}K(\tau ,j|t,\mathbf{y})=a_j(\mathbf{y},t+\tau )exp\left(-{\int }_t^{t+\tau }\sum _{j^{\prime }=1}^{2M}{a}_{j^{\prime }}(\mathbf{y},{\tau }^{\prime })d{\tau }^{\prime }\right).\hfill \end{array}$$

Proof. 

Let $P\left(\tau \right|t,\mathbf{y})$ denote the probability that no reaction occurs in the time interval $[t,t+\tau )$, given $\mathbf{y}$ (and $\mathbf{x}$, ${\mathbf{x}}_n$, ${\theta }_1$) at time t. Then, by Lemma 1:

$$\begin{array}{cc}\hfill P(\tau +d\tau |t,\mathbf{y})& =P\left(\tau \right|t,\mathbf{y})\left(1-\sum _{j=1}^{2M}{a}_j(\mathbf{y},t+\tau )d\tau +o\left(d\tau \right)\right)\Rightarrow \hfill \\ \hfill \frac{P(\tau +d\tau |t,\mathbf{y})-P(\tau |t,\mathbf{y})}{d\tau }& =-P\left(\tau \right|t,\mathbf{y})\sum _{j=1}^{2M}{a}_j(\mathbf{y},t+\tau )+P\left(\tau \right|t,\mathbf{y})\frac{o\left(d\tau \right)}{d\tau }.\hfill \end{array}$$

Letting $d\tau \to 0$ we obtain the following ODE:

$$\begin{array}{cc}\frac{d}{d\tau }P\left(\tau \right|t,\mathbf{y})\hfill & =-P\left(\tau \right|t,\mathbf{y})\sum _{j=1}^{2M}{a}_j(\mathbf{y},t+\tau ).\hfill \end{array}$$

Since $P\left(0\right|t,\mathbf{y})=1$, the solution $P\left(\tau \right|t,\mathbf{y})$ is given by:

$$\begin{array}{c}P\left(\tau \right|t,\mathbf{y})=exp\left(-{\int }_t^{t+\tau }\sum _{j=1}^{2M}{a}_j(\mathbf{y},{\tau }^{\prime })d{\tau }^{\prime }\right).\hfill \end{array}$$

We have then:

$$\begin{array}{c}K(\tau ,j|t,\mathbf{y})=P\left(\tau \right|t,\mathbf{y})a_j(\mathbf{y},t+\tau )=a_j(\mathbf{y},t+\tau )exp\left(-{\int }_t^{t+\tau }\sum _{j^{\prime }=1}^{2M}{a}_{j^{\prime }}(\mathbf{y},{\tau }^{\prime })d{\tau }^{\prime }\right).\hfill \end{array}$$

(9)

 □

Let $K_{time}\left(\tau \right|t,\mathbf{y})d\tau$ be the probability that the next reaction will occur in the time interval $[t+\tau ,t+\tau +d\tau )$, given $\mathbf{Y}\left(t\right)=\mathbf{y}$ and $(\mathbf{x}\left(t\right),{\mathbf{x}}_n\left(t\right),{\theta }_1\left(t\right))$ at time t.

Let $K_{index}\left(j\right|\tau ,t,\mathbf{y})$ be the probability that the index of the next reaction is j given $\mathbf{Y}\left(t\right)=\mathbf{y}$, $(\mathbf{x}\left(t\right),{\mathbf{x}}_n\left(t\right),{\theta }_1\left(t\right))$ at time t and given that the reaction will occur at time $t+\tau$.

By elementary probability theory (using the definition of conditional probability), we know that

$$\begin{array}{c}K(\tau ,j|t,\mathbf{y})d\tau =K_{index}\left(j\right|\tau ,t,\mathbf{y})K_{time}\left(\tau \right|t,\mathbf{y})d\tau .\hfill \end{array}$$

Due to Equation (9), we see that:

$$\begin{array}{cc}\hfill K_{index}\left(j\right|\tau ,t,\mathbf{y})& =\frac{a_j(\mathbf{y},t+\tau )}{a_0(\mathbf{y},t+\tau )}\hfill \\ \hfill K_{time}\left(\tau \right|t,\mathbf{y})& =a_0(\mathbf{y},t+\tau )exp\left(-{\int }_t^{t+\tau }{a}_0(\mathbf{y},{\tau }^{\prime })d{\tau }^{\prime }\right),\hfill \end{array}$$

(10)

where

$$\begin{array}{c}{a}_0(\mathbf{y},t)=\sum _{j=1}^{2M}{a}_j(\mathbf{y},t),\hfill \end{array}$$

and $a_0\ne 0$ due to Equation (8).

Obviously,

$$\begin{array}{cc}\sum _{j=1}^{2M}{K}_{index}\left(j\right|\tau ,t,\mathbf{y})\hfill & =1\hfill \\ {\int }_0^{\infty }{K}_{time}\left(\tau \right|t,\mathbf{y})d\tau \hfill & =1.\hfill \end{array}$$

Thus, if T is (random) time until the next reaction, then its probability density function given by $K_{time}$, its survival function $S\left(s\right)$ is given by (without loss of generality, suppose that $t=0$):

$$\begin{array}{c}\mathbb{P}(T>s)=S\left(s\right)=exp\left(-{\int }_0^s{a}_0(\mathbf{y},{\tau }^{\prime })d{\tau }^{\prime }\right),\hfill \end{array}$$

(11)

and its (cumulative) distribution function is given by $1-S\left(s\right)$ (one can check this by differentiating the distribution function, given by $1-S\left(s\right)$, with respect to s). Note that the distribution of a random variable is uniquely determined by its distribution function.

Using the proof of Proposition 2 one has the following:

Proposition 3.

Let $\tau >0$ and let $\widehat{K}\left(\tau \right|t,\mathbf{y})$ be the probability of more than one FA event occurring in the time interval $[t+\tau ,t+\tau +d\tau )$, given the state of the system at time t. Then $\widehat{K}\left(\tau \right|t,\mathbf{y})=o\left(d\tau \right)$ as $d\tau \to 0$.

Proof. 

By the proof of Proposition 2:

$$\begin{array}{c}\widehat{K}\left(\tau \right|t,\mathbf{y})=P\left(\tau \right|t,\mathbf{y})o\left(d\tau \right),\hfill \end{array}$$

since, following the definition of $a_j$, the probability of more than one reaction occurring in time interval $[t,t+d\tau )$ is $o\left(d\tau \right)$. □

Proposition 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 corresponds 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 [51,52], simultaneous events could lead to ambiguity in determining the motile state of the cell.

3.2. Combining the Cell Motility and the FA Event Model

With the results of the previous section I can now formally state the cyclical mesenchymal cell motility model as a stochastic process (see Figure 4).

Let $t=0$, $\mathbf{x}\left(0\right)$, ${\mathbf{x}}_n\left(0\right)$, ${\theta }_1\left(0\right)$, $\mu \left(0\right)$ be given and $\mathbf{Y}\left(0\right)={\mathbf{y}}^0$.

• The time $T_1$ of the FA event is chosen such that $\mathbb{P}(T_1>s)=S\left(s\right)$.
• The system evolves according to Equation (7) for $t\in [0,T_1)$.
• At time $t=\tau$, the index j of the FA event is chosen with probability $K_{index}\left(j\right|T_1,0,{\mathbf{y}}^0)$ $\mathbf{Y}$ and $\mu$ jump to new values:

  $$\begin{array}{c}\mathbf{Y}(t=\tau )=\{\begin{array}{c}{\mathbf{y}}^0+{\widehat{\mathbf{e}}}_i,\phantom{a}i=(j+1)/2,\mathrm{if}j\mathrm{is}\mathrm{odd}\hfill \\ {\mathbf{y}}^0-{\widehat{\mathbf{e}}}_i,\phantom{a}i=j/2,\mathrm{else}\hfill \end{array},\hfill \\ \mu (t=\tau )=\{\begin{array}{l}0,\mathrm{if}j\mathrm{is}\mathrm{odd}\\ 0,\mathrm{else}\end{array},\hfill \end{array}$$

  where ${\widehat{\mathbf{e}}}_i\in {\mathbb{R}}^M$ is the standard basis vector. Note that due to Equation (8), we always have $\mathbf{Y}(t=\tau )\in {\left\{0,1\right\}}^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$ the time $T_2$ of the FA event is chosen such that

  $$\begin{array}{c}\mathbb{P}(T_2>s|T_1)=exp\left(-{\int }_{T_1}^{T_1+s}{a}_0(\mathbf{y},{\tau }^{\prime })d{\tau }^{\prime }\right).\hfill \end{array}$$
• The system evolves according to Equation (7) for $t\in [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 it will be shown that this process is well-defined.

4. Piecewise Deterministic Process

In this section, I briefly overview a class of piecewise deterministic processes, first introduced by Davis [29]. I 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).

4.1. PDMP Overview

Let A be countable and let $\mathsf{\Gamma }\subset {\mathbb{R}}^d$ be open. Let ${\mathbf{X}}_t\in \mathsf{\Gamma }$ and let ${\mathbf{H}}_{\nu }:\mathsf{\Gamma }\to {\mathbb{R}}^d$ for $\nu \in A$.

Let $(\mathsf{\Omega },\mathcal{F},{\left({\mathcal{F}}_t\right)}_{t\ge 0},\mathbb{P})$ be a filtered probability space, where $\mathsf{\Omega }$ is a sample space, $\mathcal{F}$ is a $\sigma$-algebra on $\mathsf{\Omega }$, ${\left({\mathcal{F}}_t\right)}_{t\ge 0}$ is a (natural) filtration and $\mathbb{P}$ is a probability measure. Let $E:=\left\{(\nu ,\mathbf{\xi }):\nu \in A,\mathbf{\xi }\in \mathsf{\Gamma }\right\}$ and let $(E,\mathcal{E})$ be a Borel space. For details, see Chapter 2 in [30].

I can define the piecewise deterministic process on the state space $(E,\mathcal{E})$ (for a more detailed general treatment see Davis [30]) by the following objects (here, I first provide a constructive definition. The verification of the conditions and their explicit representation corresponding to our case of cell motility is postponed for the sake of clearer exposition):

I

Vector fields $({\mathbf{H}}_{\nu },\nu \in A)$ such that for all $\nu \in A$ there exists a unique global solution ${\mathbf{X}}_t\in \mathsf{\Gamma }$ to the following equation:

$$\begin{array}{cc}\frac{d}{dt}{\mathbf{X}}_t\hfill & ={\mathbf{H}}_{\nu }\left({\mathbf{X}}_t\right)\hfill \\ {\mathbf{X}}_0\hfill & \in \mathsf{\Gamma }.\hfill \end{array}$$

(12)

Let ${\varphi }_{\nu }:[0,\infty )\times \mathsf{\Gamma }\to \mathsf{\Gamma }$ denote the flow corresponding to Equation (12), i.e.,

$$\begin{array}{c}{\varphi }_{\nu }(t,{\mathbf{X}}_0)={\mathbf{X}}_t.\hfill \end{array}$$

II

A measurable function $a_0:E\to {\mathbb{R}}_{+}$ such that the function $s\mapsto a_0(\nu ,{\varphi }_{\nu }(s,{\mathbf{X}}_0))$ is integrable.

III

A transition measure $Q:\mathcal{E}\times E\to [0,1]$, such that for fixed $C\in \mathcal{E}$, $(\nu ,\mathbf{\xi })\mapsto Q(C;(\nu ,\mathbf{\xi }\left)\right)$ is measurable for $(\nu ,\mathbf{\xi })\in E$ and $Q(·;(\nu ,\mathbf{\xi }\left)\right)$ is a probability measure for all $(\nu ,\mathbf{\xi })$ on $(E,\mathcal{E})$.

Let $({\nu }^0,{\mathbf{X}}^0)\in E$ at time $t=0$ be given. Let a survival function S be defined similarly as in Equation (11):

$$\begin{array}{c}S(t,(\nu ,\mathbf{X})):=exp\left(-{\int }_0^t{a}_0(\nu ,{\varphi }_{\nu }(s,\mathbf{X}))ds\right).\hfill \end{array}$$

(13)

Let $T_1$ be the first jump time such that

$$\begin{array}{c}\mathbb{P}(T_1>t\phantom{,}|\phantom{,}({\nu }^0,{\mathbf{X}}^0))=S(t,({\nu }^0,{\mathbf{X}}^0)),\hfill \end{array}$$

and let $({\nu }^1,{\mathbf{X}}^1)$ be distributed according to the probability law $Q(·,{\varphi }_{{\nu }^0}(T_1,{\mathbf{X}}^0))$. Then, the motion of $({\nu }_t,{\mathbf{X}}_t)$ for $t\le T_1$ is given by:

$$({\nu }_t,{\mathbf{X}}_t)=\{\begin{array}{c}({\nu }^0,{\varphi }_{{\nu }_0}(t,{\mathbf{X}}^{\mathbf{0}})),\phantom{asd}t<T_1,\hfill \\ ({\nu }^1,{\mathbf{X}}^1),\phantom{aasdasdds}t=T_1.\hfill \end{array}$$

At time $t=T_1$ the next jump time $T_2$ is distributed such that

$$\begin{array}{c}\mathbb{P}(T_2-T_1>s\phantom{,}|\phantom{,}({\nu }_{T_1},{\mathbf{X}}_{T_1}))=S(s,({\nu }_{T_1},{\mathbf{X}}_{T_1})).\hfill \end{array}$$

The value of the process at the jump time $T_2$ is determined by the measure $Q(·,{\varphi }_{{\nu }_{T_1}}(T_2,{\mathbf{X}}_{T_1}))$ and the process continues in a similar way. Thus, we have a well-defined piecewise deterministic process [29].

Theorem 1.

[29] The process ${({\nu }_t,{\mathbf{X}}_t)}_{t\ge 0}$ is a homogeneous Markov process.

4.2. Cell Motility and PDMP

In this section I show that the cyclical cell motility model described in Section 3.2 is a well-defined PDMP.

One can show that ${\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1)$ satisfies the Lipschitz condition for $({\mathbf{x}}_n,{\theta }_1)\in {\mathsf{\Omega }}_{cell}\cup [0,2\pi ):=D_{cell}$ (the restriction to the interval $[0,2\pi ]$ is due to the periodic dependence on ${\theta }_1$ in the definition of ${\mathbf{x}}_i$. See Section 2). Furthermore, one can show that ${\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})$ and $\mu {\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)$, given by

$$\begin{array}{cc}\hfill {\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})& =-{\beta }_{cell}^{-1}\sum _{i=1}^M{Y}_i{\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1)\hfill \\ \hfill \mu {\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)& =-\mu {\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \sum _{i=1}^M{Y}_i{\mathbf{F}}_i({\mathbf{x}}_n,{\theta }_1)·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)\hfill \end{array}$$

also satisfy the Lipschitz condition for $({\mathbf{x}}_n,{\theta }_1)\in D_{cell}$ and arbitrary $\mu \in \left\{0,1\right\}$, $\mathbf{Y}\in {\left\{0,1\right\}}^M$.

Proposition 4.

Let $\mathbf{x}\left(0\right)={\mathbf{x}}_0$, $({\mathbf{x}}_n\left(0\right),{\theta }_1\left(0\right))\in D_{cell}$. Let $\mu \in \left\{0,1\right\}$, $\mathbf{Y}\in {\left\{0,1\right\}}^M$. Then there exists a unique solution of the system

$$\begin{array}{cc}\dot{\mathbf{x}}\hfill & =\mu {\beta }_{ECM}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\hfill \\ {\dot{\mathbf{x}}}_n\hfill & ={\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})\hfill \\ {\dot{\theta }}_1\hfill & =\mu {\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right),\hfill \end{array}$$

(14)

for $t>0$.

Proof. 

Note that since the evolution of $\mathbf{x}$ is decoupled from the other two equations, it is sufficient to prove the claim for the following subsystem:

$$\begin{array}{cc}{\dot{\mathbf{x}}}_n\hfill & ={\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})\hfill \\ {\dot{\theta }}_1\hfill & =\mu {\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \mathbf{F}({\mathbf{x}}_n,{\theta }_1,\mathbf{Y})·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)\hfill \end{array}$$

(15)

Since the right hand side of this system is Lipschitz on $D_{cell}$ and $({\mathbf{x}}_n\left(0\right),{\theta }_1\left(0\right))\in D_{cell}$, then there exists a unique solution of the subsystem Equation (15) for time $t\le t_{D_{cell}}$, where $t_{D_{cell}}=inf\left\{t^{\ast }>0\phantom{,}|\phantom{,}{\mathbf{x}}_n\left(t^{\ast }\right)\notin {\mathsf{\Omega }}_{cell}\right\}$ is the exit time from $D_{cell}$. By Corollary 1, we see that $t_{D_{cell}}=\infty$. □

Let $A:=\left\{1,2,\dots ,2^{M+1}\right\}$ and let $\mathbf{\alpha }:A\to \left\{0,1\right\}\times {\left\{0,1\right\}}^M$ be a bijection. This is simply a mapping such that $\mathbf{\alpha }\left(\nu \right)=(\mu ,\mathbf{Y})\in \left\{0,1\right\}\times {\left\{0,1\right\}}^M$ corresponds to a particular cell motion and FA states (recall that the former can either be moving or stationary).

Let $(\mathbf{x}\left(0\right),{\mathbf{x}}_n\left(0\right),{\theta }_1\left(0\right))\in \mathsf{\Gamma }:={\mathbb{R}}^2\times {\mathsf{\Omega }}_{cell}\times [0,2\pi )$ and denote ${\mathbf{X}}_t=(\mathbf{x}\left(t\right),{\mathbf{x}}_n\left(t\right),{\theta }_1\left(t\right))$. Moreover, let ${\mathbf{H}}_{\nu }:\mathsf{\Gamma }\to {\mathbb{R}}^5$ be such that

$${\mathbf{H}}_{\nu }\left(\mathbf{X}\right):=\left(\begin{array}{c}\hfill {\alpha }_1\left(\nu \right){\beta }_{ECM}^{-1}\mathbf{F}({\mathbf{x}}_n,{\mathbf{\alpha }}_2\left(\nu \right),{\theta }_1)·\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\widehat{\mathbf{r}}\left({\mathbf{x}}_n\right)\hfill \\ \hfill {\beta }_{cell}^{-1}\mathbf{F}({\mathbf{x}}_n,{\mathbf{\alpha }}_2\left(\nu \right),{\theta }_1)\hfill \\ \hfill {\alpha }_1\left(\nu \right){\beta }_{rot}^{-1}\parallel {\mathbf{x}}_n\parallel \mathbf{F}({\mathbf{x}}_n,{\mathbf{\alpha }}_2\left(\nu \right),{\theta }_1)·\widehat{\mathit{\phi }}\left({\mathbf{x}}_n\right)\hfill \end{array}\right),$$

(16)

where ${\alpha }_1$ and ${\mathbf{\alpha }}_2$ are, respectively, the first and the second components of $\mathbf{\alpha }$. In this work, the subscript index in ${\mathbf{\alpha }}_2$ always denotes the second component of $\mathbf{\alpha }$).

Let the probability $(\mathsf{\Omega },\mathcal{F},{\left({\mathcal{F}}_t\right)}_{t\ge 0},\mathbb{P})$ and state space $(E,\mathcal{E})$ be defined as in the previous section.

The objects (I,II,III) described in Section 4.1 can now be specified.

I

By Proposition 4 we see that for all $\nu \in A$, there exists a unique global solution to (12).

II

Note that in our case the rate function $a_0$ is given by (recalling Section 3.1):

$$\begin{array}{c}{a}_0(\nu ,{\mathbf{X}}_t)=a_0({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_t)=\sum _{j=1}^{2M}{a}_j({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_t).\hfill \end{array}$$

Here, I abuse the notation introduced in Section 3.1: $a_j({\mathbf{\alpha }}_2\left(\nu \right),t)=a_j({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_t)=a_j(\nu ,{\mathbf{X}}_t)$ for $j=0,\dots ,2M$. Thus, for the integrability condition to be satisfied, I assume that each probability rate function $a_j$ is integrable along the solution of Equation (12). 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 Equation (8).

III

In our case, the measure $Q(·;(\nu ,\mathbf{\xi }\left)\right)$ is given by (recalling Section 3):

$$\begin{array}{cc}Q(\left\{\eta \right\}\times d{\mathbf{\xi }}^{\prime };(\nu ,\mathbf{\xi }))={\delta }_{\mathbf{\xi }}\left(d{\mathbf{\xi }}^{\prime }\right)\sum _{j=1}^M\hfill & {\delta }_{{\alpha }_1\left(\eta \right),0}\frac{a_j^{+}({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{a_0({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_j,1}\prod _{i\ne j}^M{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_i,{\mathbf{\alpha }}_2{\left(\nu \right)}_i}\hfill \\ & +{\delta }_{{\alpha }_1\left(\eta \right),1}\frac{a_j^{-}({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{a_0({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_j,0}\prod _{i\ne j}^M{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_i,{\mathbf{\alpha }}_2{\left(\nu \right)}_i},\hfill \end{array}$$

(17)

where $\delta$ is the Kronecker delta function, ${\delta }_{\mathbf{\xi }}(·)$ is the Dirac measure at $\mathbf{\xi }$, $a_j^{+}=a_{2j-1}$ and $a_j^{-}=a_{2j}$ correspond to, respectively, the binding and unbinding probability rates at FA site j, and ${\mathbf{\alpha }}_2{(·)}_i$ is the $i^{\mathrm{th}}$ component of the vector ${\mathbf{\alpha }}_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 (10): the probability density function of the jump time $T_{k+1}$, given that $T_k\le t<T_{k+1}$, is given by:

$$\begin{array}{cc}{K}_{time}(·|t,{\mathbf{\alpha }}_2\left({\nu }_t\right))\hfill & =a_0({\mathbf{\alpha }}_2\left({\nu }_t\right),{\mathbf{X}}_t)exp\left(-{\int }_0^{·}{a}_0({\nu }_t,{\varphi }_{{\nu }_t}(s,{\mathbf{X}}_t))ds\right)\hfill \\ \hfill & =a_0({\mathbf{\alpha }}_2\left({\nu }_{T_k}\right),{\mathbf{X}}_t)exp\left(-{\int }_0^{·}{a}_0({\nu }_{T_k},{\varphi }_{{\nu }_{T_k}}(s,{\mathbf{X}}_t))ds\right),\hfill \end{array}$$

which corresponds to the survival function given by Equation (13).

I now turn our attention to the measure Q in Equation (17). The components of ${\mathbf{X}}_t$ do not jump, and vary continuously in time, i.e., if $T_k$ is the jump time, then ${\mathbf{X}}_{T_k^{-}}={\mathbf{X}}_{T_k}$ (see Section 3.2), hence the Dirac measure ${\delta }_{\mathbf{\xi }}(·)$ at $\mathbf{\xi }$ in Equation (17). Clearly, such transition of the continuous component ${\mathbf{X}}_t$ of the PDMP at a jump time is probabilistically independent of the transition of the discrete component $\nu$. 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, there is only one FA event at a jump time. Hence, for the probability of transition $\nu \to \eta$ to be nonzero, the vectors of FA states ${\mathbf{\alpha }}_2\left(\nu \right)$ and ${\mathbf{\alpha }}_2\left(\eta \right)$ must differ only by one component. Consider the following example to illustrate the jump mechanism.

Example. Let $M=4$, $T_k\le t<T_{k+1}$ and suppose $T_{k+1}$, ${\mathbf{X}}_t$ are given. Let $\nu ,\eta \in A$ be such that $\mathbf{\alpha }\left(\nu \right)=({\mu }_{\nu },{\mathbf{Y}}_{\nu })$ and $\mathbf{\alpha }\left(\eta \right)=({\mu }_{\eta },{\mathbf{Y}}_{\eta })$, where

$${\mathbf{Y}}_{\nu }=\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 1\hfill \end{array}\right],{\mathbf{Y}}_{\eta }=\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 1\hfill \\ \hfill 0\hfill \\ \hfill 1\hfill \end{array}\right].$$

Then, by Equation (10), the probability of the transition $\nu \to \eta$ is given by:

$$\begin{array}{cc}{\delta }_{{\mu }_{\eta },0}{K}_{index}\left(3\right|T_{k+1}-t,t,{\mathbf{Y}}_{\nu })\hfill & ={\delta }_{{\mu }_{\eta },0}\frac{a_3({\mathbf{\alpha }}_2\left(\nu \right),{\varphi }_{\nu }(T_{k+1}-t,{\mathbf{X}}_t))}{a_0({\mathbf{\alpha }}_2\left(\nu \right),{\varphi }_{\nu }(T_{k+1}-t,{\mathbf{X}}_t))}\hfill \\ \hfill & ={\delta }_{{\alpha }_1\left(\eta \right),0}\frac{a_3({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_{T_{k+1}})}{a_0({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_{T_{k+1}})}\hfill \\ \hfill & ={\delta }_{{\alpha }_1\left(\eta \right),0}\frac{a_2^{+}({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_{T_{k+1}})}{a_0({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_{T_{k+1}})}.\hfill \end{array}$$

Clearly, the transition ${\mathbf{Y}}_{\nu }\to {\mathbf{Y}}_{\eta }$ corresponds to the binding event at FA site 2, explaining the Kronecker delta term (see Section 2.3 and Section 3.2). Now, consider the sum in Equation (17) for this example. We see that

$$\begin{array}{c}\prod _{i\ne j}^M{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_i,{\mathbf{\alpha }}_2{\left(\nu \right)}_i}\ne 0,\mathrm{for}j=2\mathrm{only}.\hfill \end{array}$$

We therefore obtain

$$\begin{array}{c}Q(\left\{\eta \right\}\times d{\mathbf{\xi }}^{\prime };(\nu ,{\mathbf{X}}_{T_{k+1}}))={\delta }_{{\mathbf{X}}_{T_{k+1}}}\left(d{\mathbf{\xi }}^{\prime }\right){\delta }_{{\alpha }_1\left(\eta \right),0}\frac{a_2^{+}({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_{T_{k+1}})}{a_0({\mathbf{\alpha }}_2\left(\nu \right),{\mathbf{X}}_{T_{k+1}})}{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_2,1}.\hfill \end{array}$$

Note that if at time t the vector of FA states is given by ${\mathbf{Y}}_{\nu }$, then there are M possible FA state vectors into which a transition can occur with nonzero probability:

$$\begin{array}{c}\left\{\left[\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ 0\hfill \end{array}\right],\left[\begin{array}{c}1\hfill \\ 0\hfill \\ 0\hfill \\ 1\hfill \end{array}\right],\left[\begin{array}{c}0\hfill \\ 1\hfill \\ 0\hfill \\ 1\hfill \end{array}\right],\left[\begin{array}{c}0\hfill \\ 0\hfill \\ 1\hfill \\ 1\hfill \end{array}\right]\right\}.\hfill \end{array}$$

Similarly as with the rate function $a_0$, we can derive Equation (17) from the principles we established before.

Proposition 5.

The transition probability measure $Q(·,(\nu ,\mathbf{\xi }\left)\right)$ is given by Equation (17) for each $(\nu ,\mathbf{\xi })\in E$.

Proof. 

Let $(\nu ,\mathbf{\xi })\in E$, $\left\{\eta \right\}\times d{\mathbf{\xi }}^{\prime }\in \mathcal{E}$. Let $(N,\mathsf{\Xi })$ and $(N_{-},{\mathsf{\Xi }}_{-})$ be E-valued random variables before and after the jumps. Then,

$$\begin{array}{cc}Q(\left\{\eta \right\}\times d{\mathbf{\xi }}^{\prime };(\nu ,\mathbf{\xi }))\hfill & =\mathbb{P}\left((N,\mathsf{\Xi })\in \left\{\eta \right\}\times d{\mathbf{\xi }}^{\prime }|(N_{-},{\mathsf{\Xi }}_{-})=(\nu ,\mathbf{\xi })\right)\hfill \\ \hfill & =\mathbb{P}\left(\left\{\eta \right\}\times d{\mathbf{\xi }}^{\prime }|(\nu ,\mathbf{\xi })\right),\hfill \end{array}$$

where the random variables were omitted for notational convenience. Then we have:

$$\begin{array}{c}\mathbb{P}\left(\left\{\eta \right\}\times d{\mathbf{\xi }}^{\prime }|(\nu ,\mathbf{\xi })\right)=\mathbb{P}\left(d{\mathbf{\xi }}^{\prime }\right|\left\{\eta \right\},(\nu ,\mathbf{\xi }))\mathbb{P}\left(\left\{\eta \right\}\right|(\nu ,\mathbf{\xi }))={\delta }_{\mathbf{\xi }}\left(d{\mathbf{\xi }}^{\prime }\right)\mathbb{P}\left(\eta \right|(\nu ,\mathbf{\xi })),\hfill \end{array}$$

since $\mathsf{\Xi }={\mathsf{\Xi }}_{-}$ a.s., by construction of the process. Note that if $(H,\mathbf{Z})$ is a random variable, then, due to $\mathbf{\alpha }$ being a bijection:

$$\begin{array}{c}\left\{H=\eta \right\}=\left\{{\mathbf{\alpha }}^{-1}\left(H\right)={\mathbf{\alpha }}^{-1}\left(\eta \right)\right\}=\left\{\mathbf{\alpha }\left(H\right)=\mathbf{\alpha }\left(\eta \right)\right\}.\hfill \end{array}$$

Since $\mathbf{\alpha }$ is a bijection, we have

$$\begin{array}{cc}\mathbb{P}\left(\eta \right|(\nu ,\mathbf{\xi }))\hfill & =\mathbb{P}\left(({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right))|({\alpha }_1\left(\nu \right),{\mathbf{\alpha }}_2\left(\nu \right)),\mathbf{\xi }\right)\hfill \\ \hfill & =\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right),\hfill \end{array}$$

since, by construction of the cell motility process, the new FA state ${\mathbf{\alpha }}_2\left(\eta \right)$ is determined independently of whether a cell was moving or not (represented by ${\alpha }_1\left(\nu \right)\in \left\{0,1\right\}$) and the new motility state ${\alpha }_1\left(\eta \right)$ 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, only one of $j=1,\dots ,2M$ possible (binding and unbinding) FA events occurs. Thus, for $j,j^{\prime }\in \left\{1,\dots ,2M\right\}$ and $j\ne j^{\prime }$ the events “reaction j occurs” and “reaction $j^{\prime }$ occurs” are mutually exclusive. We then have, by the definition of conditional probability:

$$\begin{array}{c}\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)=\sum _{j=1}^{2M}\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|j,{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\mathbb{P}\left(j|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right),\hfill \end{array}$$

where

• $\mathbb{P}\left(j|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)$ is the probability that the FA event j occurs, given ${\mathbf{\alpha }}_2\left(\nu \right)$ and $\mathbf{\xi }$.
• $\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|j,{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)$ is the probability of a jump into cell state $\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)\right)$, given ${\mathbf{\alpha }}_2\left(\nu \right)$ and $\mathbf{\xi }$, and that the FA event j occurred.

Let $j\in \left\{1,\dots ,M\right\}$ and $j^{+}=2j-1$ and $j^{-}=2j$. Due to Equation (10) we have:

$$\begin{array}{cc}\mathbb{P}\left(j^{\pm }|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\hfill & =\frac{a_j^{\pm }({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{a_0({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}.\hfill \end{array}$$

(18)

Furthermore,

$$\begin{array}{cc}\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|j^{+},{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\hfill & ={\delta }_{{\alpha }_1\left(\eta \right),0}{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_j,1}\prod _{i\ne j}^M{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_i,{\mathbf{\alpha }}_2{\left(\nu \right)}_i}\hfill \\ \mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|j^{-},{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\hfill & ={\delta }_{{\alpha }_1\left(\eta \right),1}{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_j,0}\prod _{i\ne j}^M{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_i,{\mathbf{\alpha }}_2{\left(\nu \right)}_i}.\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\sum _{j=1}^{2M}\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|j,{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\hfill & \mathbb{P}\left(j|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\hfill \\ \hfill & =\sum _{j=1}^M\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|j^{+},{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\mathbb{P}\left(j^{+}|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\hfill \\ \hfill & \phantom{\sum _{j=1}^M}+\mathbb{P}\left({\alpha }_1\left(\eta \right),{\mathbf{\alpha }}_2\left(\eta \right)|j^{-},{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\mathbb{P}\left(j^{-}|{\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi }\right)\hfill \\ \hfill & =\sum _{j=1}^M{\delta }_{{\alpha }_1\left(\eta \right),0}\frac{a_j^{+}({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{a_0({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_j,1}\prod _{i\ne j}^M{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_i,{\mathbf{\alpha }}_2{\left(\nu \right)}_i}\hfill \\ \hfill & \phantom{\sum _{j=1}^M}+{\delta }_{{\alpha }_1\left(\eta \right),1}\frac{a_j^{-}({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{a_0({\mathbf{\alpha }}_2\left(\nu \right),\mathbf{\xi })}{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_j,0}\prod _{i\ne j}^M{\delta }_{{\mathbf{\alpha }}_2{\left(\eta \right)}_i,{\mathbf{\alpha }}_2{\left(\nu \right)}_i}.\hfill \end{array}$$

 □

Remark 2.

As mentioned in the introduction, one can derive the description of cell migration on the mesoscopic scale. If we let $f(t,\nu ,\mathbf{X})$ to be the density of cells with FA state ${\alpha }_2\left(\nu \right)$ and positions $\mathbf{X}$ at time t, then:

$$\begin{array}{cc}\frac{\partial }{\partial t}f(t,\nu ,\mathbf{X})\hfill & +{\nabla }_{\mathbf{X}}·\left({\mathbf{H}}_{\nu }\left(\mathbf{X}\right)f(t,\nu ,\mathbf{X})\right)\hfill \\ & =-a_0(\nu ,\mathbf{X})f(t,\nu ,\mathbf{X})+\sum _{\eta }{a}_0(\eta ,\mathbf{X})q(\nu ;\eta ,\mathbf{X})f(t,\eta ,\mathbf{X}),\hfill \end{array}$$

where $q(\nu ;\eta ,\mathbf{X})$ is the probability that FA state will jump from ${\alpha }_2\left(\eta \right)$ to ${\alpha }_2\left(\nu \right)$, given that the continuous component has value $\mathbf{X}$. Although it may seem daunting due to the dimensionality of variables, it is still amenable for simplifications and macroscopic scaling. For example, by assuming that subcellular (FA binding/unbinding) events occur at the fastest time scale, changes to cell position at an intermediate scale and appreciable changes to the density function at the slowest scale, one can obtain parabolic-like scaling. See Chapter 7 in [36] for details.

5. Adhesion Kinetics

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

5.1. Unbinding Rate

Consider the unbinding rate $a_j^{-}$ of FA adhesion site $j\in \left\{1,\dots ,M\right\}$ and let $\mathbf{y}\in {\left\{0,1\right\}}^M$, $\mathbf{\xi }=(\mathbf{x},{\mathbf{x}}_n,{\theta }_1)\in \mathsf{\Gamma }$. Following Bell [53], the bond disassociation rate under applied force is given by:

$$\begin{array}{c}{a}_j^{-}(\mathbf{y},\mathbf{\xi })=k_{off}^0{e}^{\parallel {\mathbf{F}}_j({\mathbf{x}}_n,{\theta }_1)\parallel /F_b}{y}_j,\hfill \end{array}$$

(19)

where $k_{off}^0$ is the FA disassociation rate under no load, ${\mathbf{F}}_i$ is the force applied at the FA, given by Equation (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 (8)). Clearly, the function in Equation (19) is integrable. Here, I 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 3.

In the context of cell migration and within the framework of our model, I 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 [54]. Hence the force dependence of the unbinding rate $a_j^{-}$. Recalling the definition of ${\mathbf{F}}_i$ in Equation (1), such indirect biochemical mediation can be included by considering mediators of the force $T_i$. In order to keep the model tractable, I omit the interaction between RhoA and myosin motors.

5.2. Binding Dynamics

Consider the binding probability rate $a_j^{+}$ of the FA adhesion site $j\in \left\{1,\dots ,M\right\}$ and let $\mathbf{y}\in {\left\{0,1\right\}}^M$, $\mathbf{\xi }=(\mathbf{x},{\mathbf{x}}_n,{\theta }_1)\in \mathsf{\Gamma }$. The probability rate $a_j^{+}$ is given by:

$$\begin{array}{c}{a}_j^{+}(\mathbf{y},\mathbf{\xi })=k_{on,j}\left(\mathbf{\xi }\right)(1-y_j),\hfill \end{array}$$

where $k_{on,j}:\mathsf{\Gamma }\to {\mathbb{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, I focus on two major mediators of FA formation: Rac activity and contractile forces.

5.2.1. Rac Dependence

For simplicity, I 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 [54], then local Rac activity correlates with the concentration of the chemical cues. Conversely, local Rac activity negatively correlates with chemorepellent.

Let $Q_{cue}:{\mathbb{R}}^2\to {\mathbb{R}}^{+}$ denote the concentration of a cue in the spatial domain and let $q:{\mathbb{R}}^{+}\to {\mathbb{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,

$$\begin{array}{c}{k}_{on,j}\left(\mathbf{\xi }\right)\propto q\left(Q_{cue}(\mathbf{x}+{\mathbf{x}}_j)\right),\hfill \end{array}$$

where I recall that ${\mathbf{x}}_j$ is the position of the $j^{\mathrm{th}}$ FA site.

For example, one can take $Q_{cue}\left(\mathbf{x}\right)$ to be the density of the ECM (or chemoattractant) at $\mathbf{x}\in {\mathbb{R}}^2$ and take $q\left(x\right)=x$. Then, the probability of binding is simply proportional to the ECM (or chemoattractant) density. This corresponds to a bi-molecular reaction rate, which depends on the number of one of the reactants (FA) and on the concentration (ECM or a chemoattractant) of the other. Although a more complex function q can be considered, such as those in [55,56,57], in order to keep the minimal character of the model, I opted for simple linear relation. Moreover, following Model 4 in [56] and assuming no feedback with phosphoinositides, then in a steady state we have:

$$\begin{array}{cc}R\hfill & =\overline{R}\frac{\alpha }{{\delta }_R}C+\overline{R}({\widehat{I}}_R+Q_{cue})\hfill \\ \rho \hfill & =\frac{\overline{\rho }}{{\delta }_{\rho }}\frac{{\widehat{I}}_{\rho }}{1+{\left(\frac{R}{a_2}\right)}^n}\hfill \\ C\hfill & =\frac{\overline{C}}{{\delta }_C}\frac{{\widehat{I}}_C}{1+{\left(\frac{R}{a_1}\right)}^n},\hfill \end{array}$$

where $R,\rho ,C$ denote the concentrations of Rac, RhoA, Cdc42 and the rest are constant parameters (see [56] for details). That is, there is a linear dependence of Rac on an external cue.

5.2.2. Force Dependence

The enlargement of nascent (immature) adhesions is concurrent with their maturation into focal adhesions [39,58]. Thus, enlargement and maturation are synonymous. While the initial stage of adhesion growth is force-independent [58], maturation occurs in a force-dependent manner [39,59,60].

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 effect of applied force diminishes [59].

That is, the study by Stricker et al. [59] 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 $\mathsf{\mu }$m² (see e.g., [61,62]).

Choi et al. [58] showed that nascent adhesions assemble at a rate of ∼1.3 ${\min }^{-1}$ = 0.021 s⁻¹ reaching a size of ∼0.2 $\mathsf{\mu }$m². Furthermore, it was shown that the formation of these adhesions is independent of the surrounding environment’s mechanical properties (such as fibronectin density and stiffness) [58] and force generation within the cell (i.e., myosin II activity) [58,63,64].

Let $k_{on}^0$ be the force-independent FA maturation (binding) rate. Since FAs are larger in size than immature adhesions, which assemble at a rate of 0.021 s⁻¹, then their assembly is slower and hence I take $k_{on}^0=0.01$ s⁻¹. It is now needed to find a function that could represent force dependence of FA maturation rate. Denote this function $k_{force}:{\mathbb{R}}^{+}\to {\mathbb{R}}^{+}$.

It satisfies the following:

• $k_{force}\left(0\right)=k_{on}^0$, i.e., if there is no force applied, the rate is $k_{on}^0$.
• If the applied force is below the characteristic force $F_b$, then $k_{force}$ 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_{force}$ should equal the unbinding rate, i.e., I 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 diminishes [59]. Thus, $k_{force}$ should plateau around $F_b$. I also assume that for large applied forces $k_{force}$ plateaus at $k_{on}^0$, since exceeding loads rupture integrin bonds frequently and impede stable maturation.

The following form of $k_{force}$ satisfies the conditions above:

$$k_{force}\left(F\right)=\{\begin{array}{c}\frac{k_{off}^{0}e+k_{on}^0}{1+exp\left(-{\gamma }_1(F-F_1^{\ast })/F_b\right)}+k_{on}^0-ϵ,\phantom{a}F\le F_b\hfill \\ \frac{k_{off}^{0}e+k_{on}^0}{1+exp\left({\gamma }_2(F-F_2^{\ast })/F_b\right)}+k_{on}^0,\phantom{a}\mathrm{else}\hfill \end{array},$$

(20)

where $F_1^{\ast }=F_b/4$ and $F_2^{\ast }=5F_b/4$ are the midpoints of the sigmoid functions (see Figure 5). To find the values of ${\gamma }_1,{\gamma }_2$ and $ϵ$, see Appendix B.

Remark 4.

Chan and Odde [65] 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 Equation (19)). Consequently, the force dependent binding rate $k_{force}$ 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\in \left\{1,\dots ,M\right\}$ is given by:

$$\begin{array}{c}{a}_j^{+}(\mathbf{y},\mathbf{\xi })=q\left(Q_{cue}(\mathbf{x}+{\mathbf{x}}_j)\right)k_{force}\left(\parallel {\mathbf{F}}_j({\mathbf{x}}_n,{\theta }_1)\parallel \right)(1-y_j).\hfill \end{array}$$

(21)

6. Numerical Simulations

Here, I show the results of simulating cell trajectories under different scenarios, which represent various experimental settings, namely:

• Uniform environment with no cues.
• Non-uniform environment with external cue gradient and uneven myosin motor activity within a cell.
• Striped ECM architecture.

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 [31,66]. Note that if $a_j^{\pm }\sim O\left(1\right)$, then the rate function is $a_0\sim O\left(M\right)$. 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 the same initial conditions for all cells are employed. Namely, at $t=0$ the conditions are:

• $\mathbf{x}$ is at the origin, ${\mathbf{x}}_n$ is uniformly distributed on a circle with radius $R_{cell}$, and $\theta =0$.
• Each FA is in (un)bound state and each cell is in moving state with probability $1/2$.

I simulate trajectories of $n_{cell}:=56$ cells for time $t_{end}:=600$ min. The time interval is divided into $n_{time}:=5000$ intervals, at the end of which the cell centroid positions $\mathbf{x}$ are recorded. For details on the parameters and numerical methods used for simulations, see Appendices Appendix B and Appendix D, respectively. For details on the data analysis, see Appendix C.

6.1. 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}$, $q=1$ in Equation (21).

The cell trajectories with varying number of adhesion sites, depicted in Figure 6a–c, show no clear trend and resemble those of a Brownian motion. Indeed, fitting the mean-squared displacement $msd\left(t\right)$ to the curve $\widehat{msd}\left(t\right)={\beta }_0{t}^{\overline{\beta }}$ (see Appendix C for details), we see that the exponent is $\overline{\beta }\approx 1$ for the three cases (see Figure 6d–f). This suggests that the cell motion has diffusive characteristics in this scenario. In Figure 6g–i, we see the simulated distribution of speeds. The average speeds $s_{av}$ (and the interquartile range, $s_{IQR}$, see Appendix B for details) and the parameters ${\beta }_0$, $\overline{\beta }$ are shown in

Table 1. Parameters computed from the simulations.

M	8	16	32
$s_{av}$, $\mathsf{\mu }$m/min	1.7595	2.4845	3.6047
$s_{IQR}$, $\mathsf{\mu }$m/min	1.4557	2.217	3.1787
$\overline{\beta }$, 1	1.0683	1.0035	1.0552
${\beta }_0$, $\mathsf{\mu }$m2/${\min }^{\overline{\beta }}$	2.1493	3.3971	4.6308
$\overline{r}$, 1	0.0483	0.0408	0.0497

. We see that as M increases, the cell motion becomes progressively faster and more diffusive (since $\overline{\beta }\sim 1$, the slope ${\beta }_0$ is a measure of diffusivity. See Appendix C), which is expected for a dominantly amoeboid type of motility. As a result of $\overline{\beta }\sim 1$, I can estimate the diffusion coefficient $D={\beta }_0/4$. The obtained values are lower, but within an order of magnitude estimated by Liu et al. [31], who found that $D\approx 2.7\mathsf{\mu }$m²/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 computed values of the average speed $s_{av}$ fall in the range reported by Liu et al. [31], who found the speeds to be in the range from $0.2\mathsf{\mu }$m/min to $7\mathsf{\mu }$m/min. Although there are very high speeds observed in Figure 6i, which seem to be outliers, speeds as high as $25\mathsf{\mu }$m/min have been observed [67]. As expected, rose plots of normalized velocities show no bias in any particular direction in Figure 6j–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 [68]. Expectedly, Figure 7 illustrates that the cells strongly deviate from straight-path migration (Figure 7 (left); see also time average of the directionality ratio $\overline{r}$ 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 6d–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. [23], Liang et al. [24] and Liu et al. [31] showed that the displacement scales superdiffusively. In these studies it was experimentally found that the time scaling went as $\sim t^{\overline{\beta }}$, where $\overline{\beta }\approx 1.2-1.5$. The primary reason why, in our case, we have diffusive behaviour is that, given a certain state of adhesion sites, there is a complete circular symmetry of the probability rates $a_j^{\pm }$ with respect to ${\mathbf{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 ${\mathbf{x}}_n$ by $\pi$ radians. Hence, somewhat reminiscent of a random walk, we have 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\sim U(-\delta ,\delta )$ is uniformly distributed on the interval $(-\delta ,\delta )$ with $\delta \in (0,1)$. Then, on average, the rates are unmodified (the multiplication factor $1+u$ for each $j=1,\dots ,M$ of every cell is computed at the beginning of simulations and is held fixed thereafter). This way, I not only simulate a non-homogeneous binding rate (and hence, for example, non-homogeneous Rac activity) within a cell, but also simulate otherwise completely identical copies of cells. Such a modification, where one does not explicitly apply a directed, predefined bias can be referred to as chemokinesis [69].

The effect of modifying the rates $a_j^{+}$ with $\delta =0.05,0.1,0.15$ can be seen in Figure 8. The cell trajectories, depicted in Figure 8a–c, show that the motion consists of periods with relatively regular path intermingled with highly irregular and random movement. In Figure 8d–f, we see that the rate modification leads to a superdiffusive time scaling of the mean-squared displacement, as the exponent $\overline{\beta }$ becomes larger than one and falls within the experimentally observed range of values [23,24,31]. Moreover, we see that as $\delta$ increases, so does $\overline{\beta }$, and the increase of the latter is more pronounced for a larger number of adhesion sites M (see also

Table 2. Parameters computed from the simulations with varying $\delta$.

M	8			16			32
$\delta$	0.05	0.1	0.15	0.05	0.1	0.15	0.05	0.1	0.15
$\overline{\beta }$, 1	0.9859	1.3184	1.4084	1.0086	1.3505	1.5581	1.1014	1.4299	1.5639
$s_{av}$, $\mathsf{\mu }$m/min	1.7656	1.7768	1.7557	2.4918	2.5009	2.5021	3.5818	3.5735	3.6104
$s_{IQR}$, $\mathsf{\mu }$m/min	1.3964	1.4319	1.472	2.24	2.2322	2.2355	3.2136	3.1108	3.0269
${\beta }_0$, $\mathsf{\mu }$m2/${\min }^{\overline{\beta }}$	3.0846	0.6934	0.5534	3.4517	0.7103	0.3543	4.1257	0.9536	0.5716
$\overline{r}$, 1	0.0452	0.0519	0.0597	0.0440	0.0513	0.0587	0.0522	0.05	0.0623

). 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 $\overline{\beta }$. However, the distribution of speeds for the corresponding values of M is virtually identical to the case with the unmodified probability rates (Figure 8g–i and Table 2). The uniform distribution of normalized velocities also remained unchanged (Figure 8j–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 [66].

Note that the increased values of $\overline{\beta }$ indicate that the cells explore a larger surface area [68]. However, other indicators of motion persistence are not affected significantly (Figure 9), although migration paths become slightly straighter, as indicated by increased values of $\overline{r}$ (Table 2). These results suggest that symmetry breaking of adhesion binding may allow cells to explore larger area without introducing velocity correlations (Figure 9d–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 [23,24,31]. Nevertheless, this highlights a potential mechanism of anomalous diffusion. In the following, I examine whether our model gives biologically consistent results in the case of externally induced polarization.

6.2. External Cue Gradient

I 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, I let the functions $Q_{cue}$ and q to have the following form (recall Equation (21)):

$$\begin{array}{cc}\hfill Q_{cue}\left(\mathbf{x}\right)& =\{\begin{array}{c}1+{\delta }_E{x}_2,\mathrm{if}{x}_2\ge 0\hfill \\ 1,\phantom{+{\delta }_E{x}_{2}a}\mathrm{else}\hfill \end{array}\hfill \\ \hfill q\left(Q_{cue}\left(\mathbf{x}\right)\right)& =Q_{cue}\left(\mathbf{x}\right),\hfill \end{array}$$

where ${\delta }_E$ represents the gradient magnitude and $x_2$ is the second component of $\mathbf{x}$. Here, for simplicity I 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, I simulate, respectively, hapto- or chemotactic migration.

In the presence of a cue gradient, we see that the cell trajectories, shown in Figure 10a–c, 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 10d–f, with the exponent $\overline{\beta }>1$ for all cases. Notice that as the number of adhesion sites M increases, so does $\overline{\beta }$ for the same ${\delta }_E$ (see

Table 3. Parameters computed from the simulations with varying ${\delta }_E$.

M	8			16			32
${\delta }_E$	0.05	0.1	0.15	0.05	0.1	0.15	0.05	0.1	0.15
$\overline{\beta }$, 1	1.2551	1.5051	1.52	1.3405	1.6963	1.7545	1.5427	1.7722	1.8569
$s_{av}$, $\mathsf{\mu }$m/min	1.8136	1.9133	2.0365	2.5235	2.5972	2.6089	3.5819	3.4218	3.3074
$s_{IQR}$, $\mathsf{\mu }$m/min	1.7384	1.5461	1.6434	2.0398	2.1999	2.1180	2.8058	2.9067	3.0175
${\beta }_0$, $\mathsf{\mu }$m2/${\min }^{\overline{\beta }}$	1.0697	0.4625	0.6845	1.0312	0.3496	0.3654	0.8319	0.3412	0.2242
$\overline{r}$, 1	0.0523	0.0607	0.0693	0.053	0.08	0.097	0.0726	0.1	0.1223

). 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 Table 1 and Table 2 vs. Table 3 and Figure 6 and Figure 8a–c vs. Figure 10a–c). In Figure 10g–i, we see that the evolution of time-averaged exponents ${\beta }_{av}\left(t\right)$ (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 ${\beta }_{av}\left(t\right)$ at $\overline{\beta }$. Curiously, a similar behaviour has also been observed by Dieterich et al. [23].

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 11d–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 12d–f). We also observe that an external signal has a stronger impact on motion persistence for higher number of adhesions due to relative increases of $\overline{r}$ 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 11g–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. In addition, for a smaller number of adhesion sites, the effects of increasing the cue gradient have a more noticeable effect on the ratios of events (see Figure 11g–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 Equation (21)) 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 ${\mathbf{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. In addition, 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 Equation (1)) and the force dependence of the FA binding rate (see Equation (20)) are necessary for directed migration resulting from biased adhesion formation in the presence of an external signal.

6.3. Fibrillar Architecture of ECM

The ECM topography is another important determinant of directed cell migration. In particular, the spatial distribution of the ECM fibres guides the motility by inducing cell shape alignment along the adhesive cues, resulting in characteristic directed movement along the fibre tracts [69]. Such guided migration is called contact guidance [32,69]. Ramirez-San Juan et al. [32] showed that contact guidance can be modulated by micrometer scale variations of interfibre spacing. Inspired by this study, I simulate how subcellular scale fibre 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:

$$\begin{array}{cc}\hfill Q_{cue}\left(\mathbf{x}\right)& =\{\begin{array}{c}1,\hspace{1em}\mathrm{if}\mathbf{x}\in {\mathsf{\Omega }}_{{\delta }_G}\hfill \\ 0.01,\phantom{,}\mathrm{else}\hfill \end{array}\hfill \\ \hfill q\left(Q_{cue}\left(\mathbf{x}\right)\right)& =Q_{cue}\left(\mathbf{x}\right),\hfill \end{array}$$

where ${\mathsf{\Omega }}_{{\delta }_G}$ represents the stripe pattern, ${\delta }_G=0.15,0.25,0.35$ represents the spacing between stripes such that the distances between them is ${\delta }_G{R}_{cell}$ (Figure 15). The stripe width is taken to be $0.25R_{cell}$. Similarly as in [32], 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 16a–c), and the displacements are aligned with the fibre pattern (Figure 16d–f). However, increasing the spacing does not simply lead to a greater adhesion alignment along the horizontal direction, as can be observed in Figure 16g–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 16h).

Such binding (and unbinding) pattern leads to fluctuating movement along northwest-southeast and northeast-southwest axis, with the resulting net migration pattern shown in Figure 16b. Similarly, the binding pattern shown in Figure 16i with more frequent events along the equator corresponds to a mixture of diagonal and horizontal movements (Figure 16c), as larger interfibre 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 [70], where it was found that adhesion alignment determines contact guidance. I also found that the average speeds were lower than in previous scenarios (Table 2 and Table 3): $1.52\mathsf{\mu }$m/min, $0.94\mathsf{\mu }$m/min and $0.87\mathsf{\mu }$m/min corresponding to, respectively, ${\delta }_G=0.15,0.25,0.35$. Interestingly, the average speeds reported in [32] were $\sim 0.6\mathsf{\mu }$m/min, although in that study the speeds were nearly constant for varying fibre pattern.

In Figure 17, I 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, I obtain the expected guidance of cell movement (Figure 16a–f) and the geometric constraint of adhesion sites (Figure 16g–i) by the fibrillar ECM pattern, in agreement with [32,69]. 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 [32,69], increasing the interfibre distance does not necessarily lead to greater alignment of cell migration along the ECM fibres in our simulations: the values of the guidance parameter G, defined as in [32] (See Appendix C for details), were found to be $0.64,0.70,0.64$, corresponding to, respectively, ${\delta }_G=0.15,0.25,0.35$. 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.

6.4. Asymmetric Contractility

I 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 [52]. 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 [34]. Here, I 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 (1) denotes the force generated by myosin motors at an adhesion site i. Instead of taking it constant, it varies with the radial position of an FA. Namely, let $T_i:[0,2\pi )\to {\mathbb{R}}_{+}$ be defined as:

$$T_i\left(\theta \right)=\{\begin{array}{c}(1+{\delta }_{myo})T_{i_0},\mathrm{if}\pi <\theta +(i-1)\frac{2\pi }{M}<2\pi \hfill \\ T_{i_0},\phantom{aasdaasad}\mathrm{else}\hfill \end{array},$$

where ${\delta }_{myo}=0.35,0.40,0.45$, $T_{i_0}$ is the constant value used in previous simulations (see Appendix B), and $\theta +(i-1)\frac{2\pi }{M}$ is the radial position of the $i^{\mathrm{th}}$ FA. Thus, contractile forces south of cell equator are larger by $35\%,40\%,45\%$ for corresponding values of ${\delta }_{myo}$. I 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 19a–c shows, as expected, the trajectories of cells maintaining north-south 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 [34]. Consequently, higher values of ${\beta }_{av}\left(t\right)$ are calculated, as shown in Figure 19d–f. In particular, for ${\delta }_{myo}=0.45$, we see that the time scaling of the mean-squared displacement is close to ballistic (see

Table 4. Parameters computed from the simulations with varying ${\delta }_{myo}$.

M	8			16			32
${\delta }_{myo}$	0.35	0.40	0.45	0.35	0.40	0.45	0.35	0.4	0.45
$\overline{\beta }$, 1	1.3072	1.6325	1.7759	1.1311	1.7905	1.8524	1.4650	1.8892	1.9353
$s_{av}$, $\mathsf{\mu }$m/min	1.6230	1.5639	1.5103	2.3742	2.3798	2.3516	3.5861	3.7414	3.7701
$s_{IQR}$, $\mathsf{\mu }$m/min	1.3827	1.4534	1.4607	2.0273	2.0186	2.2693	3.1370	3.2738	3.1152
${\beta }_0$, $\mathsf{\mu }$m2/${\min }^{\overline{\beta }}$	0.7767	0.3659	0.2493	2.3088	0.2893	0.3037	1.0713	0.9894	1.2787

for values of $\overline{\beta }$).

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 ${\delta }_{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 19g–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 ${\delta }_{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 ${\delta }_{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 ${\delta }_{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 ${\delta }_{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 ${\delta }_{myo}=0.35$. Specifically, I varied the buckling length $L_0$ and stiffness $EA$ such that $x=(1+{\delta }_x)x^0$ corresponds to the modified value of the parameter $x\in \{L_0,EA\}$, where $x^0$ corresponds to the default value given in Appendix B. In Figure 21a, I 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 ${\delta }_{L_0}=-0.18,-0.27$, respectively, leads to the opposite results (Figure 21b,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 5.

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. [34] 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. [71] showed that directed collective cell migration of neural crest cells requires greater contractility at the rear of the clump.

7. Discussion and Outlook

In this paper I constructed a stochastic model of cell migration using a minimal representation of cellular structures, essential for crawling, such as stress fibres 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. Nonetheless, several biological assumptions were imposed. The most principle of them are:

• The migration cycle is reduced to two steps: FA unbinding leads to movement and cytoskeletal reconfiguration, while binding to the latter only (Section 2).
• Cell shape is assumed to be spherical and constant (Section 2.1).
• The number of adhesion sites (occupied and unoccupied) is constant (Section 2.3).
• Adhesion binding and unbinding depend on the force applied by stress fibres and position of the cell, i.e., we omit intracellular biochemical interactions (Section 5).

We refer the reader to the referenced sections for detailed justifications.

Having specified the coupling between SFs and FAs, as well as between the cellular environment and FAs, I performed numerical simulations. I showed that our model is able to reproduce experimental observations, such as: superdiffusive scaling of the mean-squared displacement [23,24,31] (Figure 8); biased motility in the presence of external cue (Figure 10); contact guidance [32] (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 fibre tracts (Figure 11 and Figure 16d–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, I obtained directed movement in the opposite direction (Figure 20). Here, asymmetric FA dynamics (and so front-rear polarity) were also obtained, but as a consequence of locally induced contractile activity, consistent with [34].

Our simulation results in various settings suggest that the cell speeds follow a gamma distribution (Figure 6g–i, Figure 8g–i and Figure 11a–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 (Table 1, Table 2, Table 3 and Table 4). Note that faster and diffusive amoeboid movement is characterized by an increased number of weaker adhesions with high turnover and contractility [66]. Thus, the aforementioned correlation is also consistent with experimental observations. I 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 fibre tracts, where cell reshaping takes place [32]. 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 remodelling is simple, the resulting numerical simulations reproduce some important aspects of migrating cells observed in experiments. 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 fibres, 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. I 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. I believe this can be done with the framework provided by vertex-based models [72]: a more complex contractility apparatus can be described via active cable network model [41] 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 signalling pathways [5,56,73], 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 cell–cell 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 [74]. 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 [35]. 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.