Minimal Model of Directed Cell Motility on Patterned Substrates

TL;DR

This study presents a minimal biophysical model demonstrating how substrate patterning and cell adhesion parameters influence the direction of cell motility, with potential applications in cell sorting and guidance.

Contribution

It introduces a minimal model that links substrate adhesiveness and cell parameters to directional cell movement, revealing new insights into cell motility control.

Findings

Cells with low adhesion rates move perpendicular to stripes.

High adhesion rate cells move parallel to stripes.

Varying actin polymerization affects movement directionality.

Abstract

Crawling cell motility is vital to many biological processes such as wound healing and the immune response. Using a minimal model we investigate the effects of patterned substrate adhesiveness and biophysical cell parameters on the direction of cell motion. We show that cells with low adhesion site formation rates may move perpendicular to adhesive stripes while a those with high adhesion site formation rates results in motility only parallel to the substrate stripes. We explore the effects of varying the substrate pattern and the strength of actin polymerization on the directionality of the crawling cell; these results have applications in motile cell sorting and guiding on engineered substrates.

Figures

Tables

Table 1. Table 1. Physical parameter values used in numerical simulations.

parameter value description $k_0$ 2.5 wave number of substrate pattern ${\overline{a}}_0$ $.0025-.25$ linear attachment rate of adhesion sites ${\overline{a}}_{nl}$ 1.5 effective collective (nonlinear) attachment rate of adhesion sites $\overline{s}$ 1 local saturation of adhesion sites $\eta$ 10 dissipation in the adhesive layer $G$ .15 substrate stiffness $\overline{\alpha }$ $.5$ effective propulsion strength $b$ 5 sharpness of breaking function $U_c$ $\sqrt{5}$ critical extension to break adhesive contacts

Table 2. Table 2. Physical parameter values.

parameter value description $R_0$ 3 radius of cell $\beta$ 3 creation of $p$ at interface ${\tau }_1$ 10 degradation of $p$ inside cell $k_0$ 2.5 wave number of substrate pattern $\alpha$ 2 advection of $\rho$ by $p$ $\eta$ 10 dissipation in the adhesive layer $G$ .15 substrate stiffness ${\overline{a}}_0$ $.0025-.25$ linear attachment rate of adhesion sites ${\overline{a}}_{nl}$ 1.5 effective collective (nonlinear) attachment rate of adhesion sites $\overline{s}$ 1 local saturation of adhesion sites $d$ 1 detachment rate of adhesion sites $b$ 5 sharpness of breaking function $U_c$ $\sqrt{5}$ critical extension to break adhesive contacts

Equations

$$\frac{d}{dt}x$$

$$\frac{d}{dt}y$$

$$A$$

$$\frac{d}{dt}{A}_0$$

$$-\eta \frac{d}{dt}U$$

$$\epsilon \ddot{x}+\dot{x}=f(\dot{x},\dot{y},A,y)$$

$$\epsilon \ddot{y}+\dot{y}=f(\dot{x},\dot{y},A,y),$$

$${\partial }_t\rho =D_{\rho }\Delta \rho -(1-\rho )(\delta -\rho )\rho -\alpha Ap\cdot (\nabla \rho ),$$

$${\partial }_{t}p=-{\tau }_1^{-1}p-{\tau }_2^{-1}(1-{\rho }^2)p-\beta \nabla \rho -\gamma [(\nabla \rho )\cdot p]p.$$

$${\partial }_{t}A=D_A\Delta A+\rho (a_0{p}^2+a_{nl}{A}^2)-s{A}^3-d(U)A.$$

$$d(U)=\frac{1}{2}(1+\tanh [b(U^2-U_c^2)]),$$

$$\frac{d}{dt}\mathbf{U}=-\frac{1}{\eta }(G\mathbf{U}+\mathbf{V}),$$

$$\rho (x,y,t)=\rho (x-x_0(t),y-y_0(t)),$$

$$\rho (x,y)=\exp (-(x^2+y^2)/R_0^2).$$

$$V_x\int ({\partial }_x\rho {)}^2=\alpha \int Ap\cdot (\nabla \rho ){\partial }_x\rho .$$

$$V_y\int ({\partial }_y\rho {)}^2=\alpha \int Ap\cdot (\nabla \rho ){\partial }_y\rho .$$

$$p(x,y,t)=p(x-x_0(t),y-y_0(t))$$

$$-V\cdot \nabla p=-{\tau }_1^{-1}p-\beta \nabla \rho .$$

$$p(x-{\tau }_1{V}_x,y-{\tau }_1{V}_y)$$

$$=-{\tau }_1\beta \nabla \rho (x,y).$$

$$p(x,y)=-{\tau }_1\beta \nabla \rho (x+{\tau }_1{V}_x,y+{\tau }_1{V}_y).$$

$$\frac{d}{dt}{A}_0={\bar{a}}_0-d(U)A_0+{\bar{a}}_{nl}{A}_0^2-\bar{s}{A}_0^3,$$

$${\bar{a}}_0=a_0\frac{⟨\rho p^2⟩}{⟨\rho ⟩},{\bar{a}}_{nl}=a_{nl}\frac{⟨{\rho }^3⟩}{⟨\rho ⟩},\bar{s}=s\frac{⟨{\rho }^3⟩}{⟨\rho ⟩}.$$

$$-\eta {\partial }_{t}U=GU+\bar{\alpha }{A}_0.$$

$$A=A_0(t)\left(\frac{1}{2}(1+sign(\sin (k_{0}y)))\right).$$

$$f$$

$$-6e^{-(1/12)k_0^2{R}_0^2}(2k_0{R}_0^2{\tau }_1{V}_y\cos (k_0(({\tau }_1{V}_y)/3+y_0))$$

$$+(-6R_0^2+k_0^2{R}_0^4+8{\tau }_1^2(V_x^2+V_y^2))\sin (k_0(({\tau }_1{V}_y)/3+y0)))$$

$$-2e^{-(3/4)k_0^2{R}_0^2}(6k_0{R}_0^2{\tau }_1{V}_y\cos (k_0({\tau }_1{V}_y+3y_0))+(-6R_0^2+9k_0^2{R}_0^4$$

$$+8{\tau }_1^2(V_x^2+V_y^2))\sin (k_0({\tau }_1{V}_y+3y_0))))$$

$$g$$

$$+k_0^2{R}_0^4+4{\tau }_1^2(2V_x^2+3V_y^2))\cos (k_0(({\tau }_1{V}_y)/3+y_0))$$

$$+3k_0{R}_0^2(-24R_0^2+9k_0^2{R}_0^4+4{\tau }_1^2(2V_x^2+3V_y^2))\cos (k_0({\tau }_1{V}_y+3y_0))$$

$$+{\tau }_1{V}_y(-3R_0^2+4{\tau }_1^2(V_x^2+V_y^2))(-3e^{(3k_0^2{R}_0^2)/4}\pi$$

$$+12e^{(2k_0^2{R}_0^2)/3}\sin (k_0(({\tau }_1{V}_y)/3+y_0))+4\sin (k_0({\tau }_1{V}_y+3y_0)))).$$

$$\zeta =f(\zeta ,\dot{y},y)$$

$$\dot{y}=g(\zeta ,\dot{y},y).$$

Full text

MSM: Department of Mathematics, The Pennsylvania State University, University Park, PA 16802, USA. The work of MSM was supported by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program and by NSF grant DMS-1405769.

LB: Department of Mathematics, The Pennsylvania State University, University Park, PA 16802, USA. The work of LB was supported by NSF grant DMS-1405769.

ISA: Departments of Biomedical Engineering, Chemistry and Mathematics, The Pennsylvania State University, University Park, PA 16802, USA. The research of I.S.A. was supported by the US Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), Materials Science and Engineering Division.

1. Introduction

Eukaryotic cell motility is crucial to many biological processes ranging from wound healing [1] and the immune response [2] to cancer metastasis [3]. The underlying biophysical mechanisms leading to persistent cell motion are generally understood: actin treadmilling (polymerization) drives protrusions at the cell front [4, 5] while adhesion complexes transfer traction forces to the substrate [6] and myosin motors produce contractions at the cell rear [7, 8]. However the interactions of these mechanisms with external stimuli, e.g., varying substrate properties, remains relatively unexplored.

Prototypical cells for experiments, and subsequently, mathematical models, are keratocytes (e.g., harvested from fish scales). In homogeneous environments, once individual keratocytes initiate motion they exhibit characteristic crescent profiles and maintain essentially constant shape, speed, and direction [9, 10]. Moreover as the characteristic cell length/width is two orders of magnitude larger than the height while motile, keratocytes are amenable to 2D models and thus may be considered the simplest cells for development of mathematical models (for a review on advances in 3D modeling techniques see, e.g., [11]).

For a general overview of both biological aspects of actin driven cell motility and of several modeling approaches we recommend the survey [12]. In particular, in recent years both free boundary and phase-field models have been extremely successful in replicating, explaining, and predicting experimental results (see, e.g., [13, 14, 15, 16, 17, 18, 19]), see for review [20]. We highlight that recent mathematical analyses have elucidated biological mechanisms [21, 22, 23, 24], and numerical simulations have described a wide range of behaviors ranging from motility initiation via stochastic fluctuations [25] to capturing various modes of motility such as stick-slip and bipedal motions [26, 27, 28].

The focus of our study is the effect of variable substrate properties on the direction of the cell’s motion. In experimental settings the substrate may be coated non-homogeneously with fibronectins which allow for the creation of adhesion complexes, resulting in regions of variable adhesiveness. In [29] microcontact printing of regions enriched with (or depleted of) fibronectins resulted in alternating substrate stripes of high and low adhesiveness with period smaller than the cell size. In this setting the authors of [29] found that keratocytes exhibited directed motion parallel to the adhesive stripes.

We briefly mention other works which study the effect of non-homogeneous substrates. If the cell’s size is smaller than the width of the adhesive stripe then both experiment and modeling [30, 28] observe that the cell may be contained within adhesive regions: when the cell encounters a region of low adhesiveness it will change its direction in order to remain on regions of high adhesiveness. Moreover, it has been observed that cell morphologies can be controlled when placed on specific adhesive geometries [31] and directionality can be controlled on substrates with varying stiffness [32, 33]. As such, controlling and predicting the motility of keratocytes on engineered patterned substrates has direct applications to cell screening and sorting both in biomimetics [34] and experimental settings.

In this work we study the directionality of a motile cell on a patterned substrate with alternating stripes of high and low adhesiveness; here the cell size is assumed to be larger than the stripe width so that the cell spans several regions of both high and low adhesiveness corresponding to the experimental setup of [29]. We stress that the dynamics of adhesion site formation are, in general, complex [35, 36, 37, 38]. In [28] a phase-field model which included both dynamics of the adhesion site formation as well as substrate deformation (e.g., due to traction forces generated by the cell) was introduced. The authors introduced patterned adhesiveness and reproduced experimental results observed in [29]. Surprisingly, the authors of [28] also obtained for certain parameter regimes that the cells move perpendicular (rather than parallel) to substrate adhesion stripes. More specifically, when the attachment rate of adhesion complexes is sufficiently high (e.g., when adhesive stripes had high adhesiveness), numerical simulations reproduce experimental results. When the attachment rate of adhesion complexes is low (e.g., when adhesive parts of stripes have low adhesiveness), numerical simulations show that cells move perpendicular to substrate stripes; this striking behavior has thus far remained unobserved in experiments and merits additional numerical investigation.

Due to the complexity of the phase-field model in [28] there is no clear, simple mechanism to explain and differentiate parallel/perpendicular motions. To elucidate the effects of substrate adhesiveness patterns and biophysical parameters on the resultant direction of motion, we derive from the full system a minimal model of cell motility on patterned substrates which is capable of reproducing both parallel and perpendicular (to stripes) directed motions. This minimal model reveals the underlying mechanical processes which may give rise, in particular, to perpendicular motions. Moreover the reduction allows for efficient numerical experiments over a wide range of parameters revealing the dependence of direction of motion on biophysical (e.g., actin polymerization strength) and substrate properties (e.g., stripe sizes, adhesiveness).

2. Results

2.1. Description of the reduced model

In [28] a 2D phase-field model of cell motility was introduced to describe the onset and persistence of cell motility as well as a broad range of cell morphologies. The model contains four coupled differential equations: a phase-field equation describing the location of the cell membrane, an evolving vector field representing the effect of the actin filament network, a scalar equation for density of adhesion sites, and a Kelvin-Voigt visco-elastic equation for the deformation of the substrate.

As exploration of the full phase-field model is complicated we propose the following simplified system of differential and algebraic equations which track the location of the center of the cell $(x,y)$ as well as the effective adhesion of the cell to the substrate $A$ and effective substrate deformation $U$ (see Figure 1 for an illustration of the model):

[TABLE]

Here, $V_{x},V_{y}$ are the center of mass velocities, $d(U):=\frac{1}{2}(1+\tanh[b(U^{2}-U_{c}^{2})])$ is an effective cut-off function which destroys adhesions if the substrate deformation $U$ exceeds some critical value $U_{c}$, $k_{0}$ is the wavenumber of microprinted stripe pattern. Equation (3) encodes the striped substrate by creating inhomogeneity in the $y$-direction.

System (1)-(5) is derived from the phase-field model in [28] after several simplifying assumptions; a detailed derivation can be found in the Supplementary Materials. In this reduced model we aim to replicate both the numerically observed emergence of perpendicular motion to the stripes as well as the experimentally observed parallel motions whereby we may study the robustness of both as well as understand the biophysical mechanisms which give rise to each.

We note that the system (4)-(5) is decoupled from equations (1)-(2) so that $A=A(t)$ enters (1)-(2) as a time varying coefficient depending on physical parameters. Since we will focus on the effect of varying initial conditions and physical parameters on the motion of the cell, the decoupling yields efficient numerical simulation. Indeed, we can compute $A_{0}$ and $U$ independently form equations for $V_{x},V_{y}$, and thus to reduce overall computation time and effectively decrease the search over initial conditions by two dimensions. See Section 4 for a detailed description of the numerical methods used.

Although this decoupling is a vast simplification to the physical model, subsequent numerical study shows that we retain sufficient structure for meaningful results; future work may investigate the fully coupled system by replacing (5) by the vector equation $-\eta\frac{d}{dt}{\bf U}=G{\bf U}+{\bf V}$, where ${\bf V}=(V_{x},V_{y})$.

2.2. Smaller driving force required for motion on striped substrate

We first consider the effect of a striped substrate on the minimal driving force required for persistent cell motion. As observed in [28] if we take $\bar{a}_{0}$ sufficiently small, e.g., $\bar{a}_{0}=.0025$, then the system (4)-(5) tends to a limit cycle. On the other hand taking $\bar{a}_{0}$ large, say $\bar{a}_{0}=.25$, (4)-(5) tends to an equilibrium, suggesting the existence of a supercritical Hopf bifurcation in the parameter $\bar{a}_{0}$; the point of bifurcation is $\bar{a}_{0}\approx.1$. For simplicity we take $\bar{a}_{0}=.25$ so that the pair $(A_{0},U)$ tends to an equilibrium $(A_{\infty},U_{\infty})$. Following [28], we define $\kappa:=8\alpha\beta\tau_{1}^{2}A_{\infty}/(81R_{0}^{2})$, the normalized driving force in the cell. In [28] a subcritical onset of motion was observed on homogeneous substrates: for $\kappa<\kappa_{c}\approx.746$ there is no persistent motion of the cell while for $\kappa>\kappa_{c}$ there is some finite minimal velocity of the cell.

Upon inclusion of substrate stripes the minimal driving force required for the onset of persistent motion decreases: $\kappa_{c}\approx.395$. As expected, this motion is parallel to the substrate stripes. This predicts that microcontact printing may allow for directed cell motion even for cells which cannot sustain persistent motion on homogeneous substrates. These results are analogous to observations in [39] where it was shown that the distance traveled by a motile cell is larger when restrictions to the cell geometry are imposed (e.g., diameter of the background matrix is made smaller than the cell diameter). Moreover results of [40] show that micro-contact printing may dominate chemical cues for contraction driven motile cells (e.g., fibroblasts).

2.3. Robustness of vertical trajectories decreases as adhesion formation rate increases

Since $f(0,V_{y},A,y)=0$ then $V_{x}=0$ is always a trivial solution of $V_{x}=f(V_{x},V_{y},A,y)$. We thus first investigate the existence of persistent motion perpendicular to the substrate pattern in the reduced system (2)- (5) with $V_{x}=0$. Numerical analysis shows that persistent motion perpendicular to stripes is possible over a large range of physical parameters. For subsequent results we assume to use parameter values as in Table 1 unless otherwise mentioned.

Surprisingly our numerical analysis also shows persistent vertical trajectories exist for all values of $\bar{a}_{0}$ (in contrast to results obtained in the full PDE model). First taking the effective linear attachment rate $\bar{a}_{0}=.0025$, the $A_{0}$-$U$ system exhibits a limit cycle (i.e., stick-slip behavior), which in turn results in large oscillations in the $y$ velocity, see Figure 1. Taking $\bar{a}_{0}=.25$, then the $A_{0}$-$U$ system tends to equilibrium but a purely perpendicular motion still exists. Oscillations in the $y$ velocity result only from the non-homogeneity of the substrate and thus are smaller than in the case of oscillating $A_{0}$ and $U_{0}$. Since persistent motion perpendicular to stripes is observed for all values of $\bar{a}_{0}$, in order to corroborate our results with the full PDE simulations, we aim to understand the robustness of perpendicular trajectories via a stability analysis of the vertical motion as a function of $\bar{a}_{0}$. That is, we consider the long time direction of motion if the initial $x$-velocity $V_{x}$ is chosen to be non-zero. Due to the algebraic dependence of $V_{x}$ in (1) classical linear stability analysis techniques are not available. Thus, to study stability of vertical trajectories, we exhaustively search initial conditions numerically in order to see long time asymptotics, for details see Section 4.

We define $\theta$ to be the angle (in degrees) of the initial velocity of the cell from the positive vertical axis, see Figure 1. For each $\bar{a}_{0}$ we exhaustively search initial conditions $y_{0}$ in order to determine the maximal value of $\theta$ which gives rise to persistent vertical motion. We define $\max\theta$ to be the maximal $\theta$ over all possible initial conditions. The dependence of $\max\theta$ on $\bar{a}_{0}$ is shown in Figure 2.

We observe that $\max\theta$ is largest where $a_{0}$ is small, indicating that vertical motion is more robust in this regime. This corresponds with the results in [28]; in particular $\max\theta$ is largest in the regime that the cell undergoes stick-slip motion (i.e., limit cycles in the $A_{0}-U$ system).

When $\bar{a}_{0}$ is sufficiently large, the pair $(A_{0},U)$ tends to an equilibrium $(A_{\infty},U_{\infty})$. The equilibria $A_{\infty}$ and $U_{\infty}$ both increase with $\bar{a}_{0}$. We observe that $(A_{\infty},U_{\infty})$ remains relatively constant for all $\bar{a}_{0}<.3$ (e.g., $.65\leq A_{\infty}\leq.7$). However, $(A_{\infty},U_{\infty})$ begins to very rapidly increase for $\bar{a}_{0}>.3$ (e.g., $A_{\infty}\approx.88$ for $\bar{a}_{0}=.4$). Interestingly $U_{\infty}$ is approximately $U_{\infty}\approx R_{0}$ when $\max\theta=0$. Although substrate variations of the order of magnitude of the cell size are unphysical (an artifact of our approximations), we conclude that large substrate variations may lead to stabilization of parallel motion whereas small substrate variations allow for perpendicular motion to stripes. This hypothesis agrees with previous numerical and experimental evidence which shows that cells may overcome variations in substrate stiffness provided the substrate is sufficiently stiff [28, 30].

2.4. Aspect ratio of substrate stripes sorts cells depending on actin polymerization strength

We investigate the motion of cells on substrates with striped patterns of adhesiveness where the ratio of adhesive stripe width to non-adhesive stripe width is not necessarily equal. Let $L_{1}$ be the width of the adhesive stripe and $L_{2}$ to be the width of the non-adhesive stripe so that $L:=L_{1}+L_{2}$ is the period of the substrate pattern. We investigate the effect of varying the ratio $L_{1}/L$ on $\max\theta$.

Here, we keep $\bar{a}_{0}$ constant ($\bar{a}_{0}=.0025$) and vary $L_{1}$. Since $\alpha$ is a key physical parameter measuring the strength of actin polymerization, we additionally investigate how changing $\alpha$ changes these data. The results are summarized in Figure 3.

It is clear that if $L_{1}=0$ or $L_{1}=L$ then the cell cannot create a biased directionality and $\max\theta=0$. This is expected since if the adhesive stripe is too small then the cell cannot develop a sufficient number of adhesive bonds to the substrate to initiate persistent motion. On the other hand if $L\approx 1$ the substrate is entirely adhesive and the substrate is effectively homogeneous.

Interestingly we observe a monotonic dependence between the actin polymerization strength $\alpha$ and the percentage of adhesiveness of the substrate which results in perpendicular motion to stripes. If the substrate is predominantly non-adhesive then the cell requires high actin polymerization strength to generate perpendicular motion and vice-versa. This suggests that different effectiveness of internal biophysical parameters may lead to different behaviors of the cell depending on the shape of the patterned substrate, providing evidence for cell sorting and directed cell motility.

2.5. Memory increases the basin of attraction of persistent perpendicular motion

Thus far, the basin of attraction for vertical motion has been bounded by $\max\theta<40$. However in full PDE simulations even with motion initialized parallel to the stripes, the cell switches directions and begin moving perpendicular to the stripes (provided $a_{0}$ is sufficiently small) so that $\max\theta\approx 90$. We predict that this discrepancy arises from the lack of memory in the reduction (1) - (5). That is, we neglect both shape deformation and directional memory (or inertia) which are present in the full PDE model. In particular, there is finite time relaxation of the actin vector field which is neglected in the reduced equations.

To include these memory effects we modify (1)-(2):

[TABLE]

[TABLE]

where $\varepsilon$ is an effective memory (or )inertial) coefficient characterizing the time of relaxation of the cell’s directional inertia. This coefficient is a function of both shape deformations and relaxation of the actin polymerization field. Mathematically it is also important to note that $\varepsilon>0$ provides regularization so that (6)-(7) absent in (1)-(2). We compare numerical results with those results obtained for $\varepsilon=0$, see Figure 2. Again, our interest is to study the dependence on $\max\theta$ as a function of $\bar{a}_{0}$.

As (6)-(7) is a singular perturbation of (1)-(2), it is natural that the behavior drastically changes qualitatively. However, it is surprising that the value of $\max\theta$ becomes approximately piecewise constant in numerical studies. Moreover the maximum value of $\bar{a}_{0}$ for which $\max\theta>0$ decreases as $\varepsilon$ increases: a cell with longer memory is less likely to overcome the nonadhesive stripes and exhibit perpendicular motion. However we observe that even for $\varepsilon=.005$ the positive values of $\max\theta$ are very large: $\max\theta\approx 77$, agreeing with simulations in [28]. This suggests that memory of shape deformations as well as persistence of the actin network may be correlated with persistent motion perpendicular to the stripes.

3. Conclusions

We have presented a reduced system derived from the full phase-field model presented in [28]. The reduced system is much less computationally expensive to simulate while still retaining the crucial solution behaviors. We reproduce both perpendicular and parallel persistent cell motions on periodically striped substrates. Our numerical study indicates that the robustness of vertical and horizontal motions can be quantified via consideration of long time asymptotics over all initial conditions (measured by $\max\theta$). These simulations indicate that low adhesion site formation rates are necessary for persistence of perpendicular motions. Physically we expect that if the number of adhesion sites is too high then the forces “pulling” the cell back towards the adhesive part of the substrate will prevent the cell from overcoming the arrest caused by the nonadhesive part of the substrate.

Moreover, we exhibited a inverse, monotonic dependence of actin polymerization strength on the amount of adhesive substrate stripes allowing for perpendicular motion. We predict that carefully engineered substrates could potentially use this correlation in applications of cell sorting and directed cell motility.

Finally, by including memory (e.g. due to finite cell deformation relaxation time) we obtain a singularly perturbed system which shows better agreement with full phase-field simulations: for small values of $\bar{a}_{0}$ the initial velocity may be almost parallel to stripes and still we observe persistent perpendicular motion over long time. These results suggest that shape deformation and directional cell memory correlate with the ability of the cell to move perpendicular to the substrate stripes. We believe that this methodology of reduction of a full PDE model to finite dimensional models (e.g., systems of ordinary differential equations) may be applicable to a wide range of physical models and may provide new insights via both analytical and numerical study of the resultant, simplified systems.

4. Methods

We solve system (1) - (5) using a forward finite difference scheme. In particular $V_{x}$ is defined via an algebraic equation, as opposed to a differential equation. As such to solve the implicitly defined velocities $V_{x}$ and $V_{y}$ we require a nonlinear function solver. We use a predefined function in Scilab (fsolve), which is based on an iterative Powell hybrid method, initializing the solver at the solution of the previous time step.

In order to conduct an exhaustive search through all initial conditions we first simulate the longtime dynamics of (4)-(5) to determine either the value of the equilibrium point or the stable limit cycle. In the former case, we assign the limiting equilibrium value as the initial conditions for $A_{0}$ and $U$. In the latter case, we may fix any point on the limit cycle (for consistency, we choose the point on the limit cycle where $A_{0}$ is maximal).

When viewed as an algebraic system, (1)-(2) has non-unique solutions $(V_{x},V_{y})$ for fixed $y$. Thus for a fixed value of $y_{0}\in[0,k_{0}/2\pi]$ we first compute all admissible initial velocities $V_{x}(0),V_{y}(0)$ and simulate the long time dynamics of all to deduce if either (i) the motion becomes eventually horizontal, (ii) the motion becomes eventually vertical, or (iii) there is no continuous in time velocities which solve (1)-(2) for all time. Since the equations (1)-(2) are nonlinearly coupled, case (iii) may arise, for example, as the result of bifurcations in the $V_{x}$-$V_{y}$ solution plane as $A$ or $y$ vary. We recall that there are generically non-unique solutions to (1)-(2), so if such a discontinuity occurs, it is deemed unphysical and so we omit this scenario from our analysis.

We additionally note that solving (6)-(7) requires small time steps $\Delta t\ll\varepsilon$ to ensure convergence as it is a singularly perturbed system. In general we may take any initial velocities $V_{x}(0),V_{y}(0)$, however we restrict ourselves to initial conditions which are compatible with initial velocities computed in the case that $\varepsilon=0$ so that results may be compared: given an initial $y_{0}$ we initialize $V_{x}(0),V_{y}(0)$ to be admissible solutions to the system (1)-(2).

5. Appendix

5.1. Phase-field model

In [41] a 2D phase-field model of cell motility was introduced to describe the onset and persistence of cell motility as well as a broad range of cell morphologies. This model contains two PDEs. First, a phase-field parameter $\rho$ describes the location of the cell membrane (i.e., $\rho\approx 1$ on the interior of the cell and $\rho\approx 0$ outside of the cell). Second, the motion of the cell membrane is actively driven by a vector field $p$ which models the averaged orientation field of the actin filament network.

The coupling of these equations reflects two main experimental observations: (i) the nucleation of branches of actin filaments near the cell membrane via Wiskott-Aldrich syndrome proteins (WASP) activation of the ARP2/3 complex and (ii) presence of actin filaments allows for the creation of adhesive contacts and subsequent transfer of momentum of the polymerizing actin network to the substrate driving the cell membrane. Additionally the model enforces approximate volume preservation. The full PDE model introduced in [41] is:

[TABLE]

[TABLE]

For our subsequent analysis we assume that $\tau_{2}^{-1}=\gamma=0$. In particular, the term containing $\gamma$ accounts for symmetry breaking due to myosin driven contraction in the rear of the cell; since our subsequent analysis initializes data with non-zero velocities, this symmetry breaking is not required for capturing persistent motion. Moreover it is mentioned in [28] that self-sustained cell motion is possible in the system (8)-(9) even without this term. We refer the interested reader to [28, 41] for a full description of the model.

In [41] it is assumed for simplicity that the friction generated from adhesion complexes is homogeneous so that $A$ is constant. To account for the complex interaction between cell and substrate it is necessary to include the dynamics of adhesion site formation, see [28]:

[TABLE]

These adhesion contacts describe integrin complexes which (through a series of intermediate proteins such as zyxin, talin, and vinculin) engage both the substrate and the cytoskeleton.

In the final term of (10), the function $d(U)$ describes the coupling of $A$ to the average substrate deformation ${\bf U}(t)$, with $U=|{\bf U}|$:

[TABLE]

where ${\bf U}$ satisfies

[TABLE]

where $\eta$ is the effective viscous friction of the substrate, $G$ is an effective spring constant and ${\bf V}$ is the velocity of the cell’s center of mass. That is, the substrate is viewed as a Kelvin-Voigt visco-elastic material and adhesion sites are broken if deformations exceed the threshold $U_{c}$.

As all attachment/detachment rates are effective parameters, they incorporate both characteristics of the adhesion complex as well as the substrate preparation. Thus, spatial inhomogeneity in the substrate may be introduced through coefficients, e.g., spatial dependence of $a_{0}=a_{0}(y)$.

5.2. Derivation of the reduced system

We reduce the system (8) to a two-dimensional system for the location of the center of the cell $(x,y)$. All assumptions are analogous to those made in [28]. The biggest difference is that we do not assume that $V_{y}\ll 1$ and as such we have more complex coupling between all equations. We first assume that the cell has fixed circular shape for all time:

[TABLE]

where

[TABLE]

In particular, note that under the assumption (13) we have $\delta\equiv\frac{1}{2}$. We denote the velocity of the center of the cell $V=(V_{x},V_{y})=(\dot{x}_{0},\dot{y}_{0})$. Multiplying (8) by $\partial_{x}\rho$, and integrating over the domain yields

[TABLE]

Likewise multiplying (8) by $\partial_{y}\rho$ and integrating we have

[TABLE]

We likewise assume that

[TABLE]

to rewrite equation (9):

[TABLE]

Under the assumption that $|V|$ is small:

[TABLE]

Thus to first order

[TABLE]

To simplify the equation for $A$ we first consider the problem on a homogeneous substrate (i.e., $a_{0}$ is constant) and consider the following ansatz for the adhesion sites density: $A=A_{0}(t)\rho(x,y,t)+\delta A_{1}(x,y,t)$ where $\delta$ is much smaller than the size of the cell. That is, we assume that $A$ is essentially spatially constant on the interior of the cell. Plugging this expansion into (10) and integrating over the domain we have to leading order:

[TABLE]

where

[TABLE]

Similarly we can simplify the equation for ${\bf U}$ to a scalar equation under the assumption that $U$ and $V$ are essentially co-linear. As in [28] we make the approximation $V(t)\approx\langle\alpha A_{0}p\rangle\approx\bar{\alpha}A_{0}$, where $\bar{\alpha}$ is a numerical constant:

[TABLE]

Finally, in order to simulate the patterned substrate we approximate $A$ to have the form

[TABLE]

For subsequent numerical simulations, we use the first two terms of the Fourier expansion of (25) for simplicity and regularity.

By plugging in (21), (25) into (15)-(16), we derive the equations $V_{x}=f(V_{x},V_{y},A,y)$, $V_{y}=g(V_{x},V_{y},A,y)$, where $f$ and $g$ are defined

[TABLE]

For all simulations we use parameter values from Table 2 unless otherwise noted.

If we assume that $V_{y}=0$ we of course recover equation (14) in [28]. We note that due to homogeneity of the substrate in the $x$ direction, equation $V_{x}=f(V_{x},V_{y},A,y)$ is algebraic; letting $\zeta=V_{x}$:

[TABLE]

[TABLE]

Importantly (26)-(27) may not be in general solvable (or have a continuous in time solution) for all values of $\dot{y},y$. This leads to difficulty of both numerical simulations as well as prevents the use of classical linear stability analysis.

We highlight here that the reduction described above is similar to the one conducted in [28] however we additionally incorporate expansions for the $V_{y}$ component. Moreover the equation for the effective adhesion $A_{0}$ is derived using the asymptotic expansion $A=A_{0}(t)\rho(x,y,t)+\delta A_{1}(x,y,t)$ which gives rise to the effective coefficients (23).