# Cellular tango: How extracellular matrix adhesion choreographs Rac-Rho
signaling and cell movement

Elisabeth G. Rens  Current: Delft Institute of Applied Mathematics, Delft
University of Technology, Delft, the Netherlands, Previous: Department of
Mathematics, University of British Columbia, Vancouver, Canada  and  Leah
Edelstein-Keshet  Department of Mathematics, University of British Columbia,
Vancouver, Canada

###### Abstract.

The small GTPases Rac and Rho are known to regulate eukaryotic cell shape,
promoting front protrusion (Rac) or rear retraction (Rho) of the cell edge.
Such cell deformation changes the contact and adhesion of cell to the
extracellular matrix (ECM), while ECM signaling through integrin receptors
also affects GTPase activity. We develop and investigate a model for this
three-way feedback loop in 1D and 2D spatial domains, as well as in a fully
deforming 2D cell shape. The model consists of reaction-diffusion equations
solved numerically with open-source software, Morpheus, and with custom-built
cellular Potts model simulations. We find a variety of patterns and cell
behaviors, including persistent polarity, flipped front-back cell polarity
oscillations, and random protrusion-retraction. We show that the observed
spatial patterns depend on the cell shape, and vice versa. The cell stiffness
and biophysical properties also affect patterning, overall cell migration
phenotypes.

##  1\. Introduction

Cell migration is a vital phenomenon that occurs in health and disease,
including wound healing and cancer metastasis. Signals that promote cell
migration include chemical gradients, topographic or mechanical cues or
adhesion gradients [1, 2, 3]. Properties of a cell’s internal regulatory
system, and the various aspects of its environment, together determine the
spatiotemporal distribution of intracellular signaling, and the resultant cell
response [4]. The onset, direction, persistence and type of cell migration has
important biological significance. During wound healing, cells need to start
moving in the right direction to properly close a wound. In cancer,
spontaneous persistent migration is a major determinant of metastasis. Cell
mutations and cell response to changes in the environmental cues affect
whether cell migration is normal, as in wound healing and development, or
defective, as in cancer metastasis.

Cell motility is regulated by signaling networks in the cell, where proteins
of the Rho-family GTPases, called Rac and Rho are central hubs [5, 6, 7, 8].
These proteins can exist in an activated, membrane bound form, or inactive
form that freely moves in the cytosol [9]. Their (in)activation is regulated
by GEF/GAP proteins. Rac is generally localized in protrusions, where it
promotes actin assembly to further push the cell membrane out [10]. Rho is
generally activated at the trailing ends of the cell, where it promotes cell
contraction via ROCK, a kinase that recruits and activates myosin on the actin
cytoskeleton [11, 12]. Together, these two opposing actions of Rac and Rho on
the cytoskeleton regulates the location of protrusive or contractile
activities [13, 14, 15].

One typical, and relatively well-studied, example of a Rac/Rho pattern is two
opposing gradients, consistent with persistently polarized cell motion [4].
However, the spatial distribution of Rac/Rho is not always this simple.
Complex spatial and temporal patterning have been observed. For instance, Rho
exhibits transient bursts of activity at the trailing edge [16]. Also, waves
of Rho from the front to the back of the cell have been observed [17]. The
activity at multiple sites of protrusion is often associated with high Rac,
but bands of intense Rho activity have been also been observed in active
protrusions [16]. It is suggested that Rac and Rho actually cyclically
interchange their activity in protrusions [18]. All this shows that the
dynamics and distribution of Rac-Rho are non-trivial and depend on the exact
internal signaling network, cell geometry and environmental inputs.

Mutations can affect protein conformations, and their binding and activation
rates, and this could bias the competition between Rac and Rho and thus the
expected cell behavior. Signals from the environment also play a significant
role. Cells are surrounded by the extracellular matrix (ECM), a network of
proteins and fibers [19]. The mechanical architecture of the ECM, a topic of
increasing interest in cancer and wound healing, is an important determining
factor for cell phenotype [2]. Signals from outside of the cell, such as
external forces [20, 21], cell-ECM bonds [22, 23] or growth factors [24], are
channeled to GTPases like Rac, Rho and Cdc42, and promote or inhibit their
activation or inactivation. Cell protrusions interact with the ECM, causing
external cues to change local GTPase activity, which by diffusion also affects
the signaling profile of the whole cell.

The interplay between ECM and cell signaling has received much attention in
recent years. For instance, one study [22] looked at cell migration on
patterned adhesive islands. Their experiments show that the precise local
pattern determines where Cdc42 was activated. The authors went on to show that
with the right adhesive pattern, cell migration could be reversed. Another
study used different sizes of and spacing between ECM islands to show that the
onset of Rac waves correlated with local adhesion, affecting the cells’
orientation of motion [25].

Various mathematical models have considered GTPase dynamics within a cell but
most of these studies only focused on cell polarization, such as [26, 27, 28,
29, 30, 31, 32, 33, 34, 35] and others. The role of cell-ECM adhesion has been
mathematically modelled and studied in detail in the context of cell spreading
and directional migration (see e.g. [36, 37, 38, 39, 40] and many others), but
not so much towards GTPase patterning [41]. So, what still remains elusive is
how various, more diverse phenotypes in Rac-Rho patterning can emerge. Here we
focus specifically on the patterning in the context of ECM signaling.

Even though there are many complex interactions between a cell and the ECM, it
was proposed that experimental observations of melanoma cells [42] can be
explained by a minimal model of Rac-Rho-ECM [42, 43] where adhesion promotes
Rho activation, Rac and Rho are mutually antagonistic, Rac upregulates
adhesion and Rho down-regulates adhesion. Both these papers motivate our own
work. There, the dynamics were studied in a 2-compartment model, representing
the front and the back of a cell.

In this paper, we will explore the effect of this ECM feedback mechanism on
the full spatiotemporal dynamics of Rac and Rho. We focus our plan on
addressing the following questions (1) Can we find the same dynamic cell
polarity patterns as described in [42] in the full spatial models? What are
minimal requirements to get the oscillatory patterns of behaviour observed in
melanoma cells in [42]? (2) Do realistic force-dependent adhesion dynamics
give rise to the same patterning as in [42]? (3) What are the possible GTPase
patterning dynamics in the spatial model? How does the geometry affect those
dynamics? (4) How does the feedback between GTPase activity and cell
protrusion/retraction affect the resulting pattern dynamics and the cell
behaviour? (5) How do biophysical properties such as cell stiffness affect the
patterns and resulting cell phenotype?

Figure 1. A schematic diagram of the models considered in this paper. (a) The
well-mixed mutual inhibition model, for Rac -Rho without and later with
feedback from the extracellular matrix (ECM). (b)-(d) spatially distributed
model variants where PDEs describe the basic interactions distributed across
the cell. Rac and Rho then have both active and inactive forms that diffuse at
distinct rates. ECM signaling is distributed spatially as well. The domain
geometry considered is (b) a 1D strip, (c) a static circular domain and (d) a
full 2D version of the model with cell shape captured by a cellular Potts
computation. In all three cases, we represent the results in a kymograph, with
the spatial variable horizontal. For 2D geometries, the spatial variable is
the perimeter of the full 2D domain and $-\pi\leq x\leq\pi$.

###  1.1. The extracellular matrix (ECM)

The ECM is a meshwork of fibrous proteins such as collagen surrounding cells.
It presents a complex topography and adhesive environment on which cells
crawl, pull, deform, and remodel [19, 44]. As cells pull and exert forces on
the ECM, the resulting mechanical tension creates stimuli that affect the
cells.

Cells attach to ECM via integrin bonds, that are distributed over the cell
surface. A migrating cell creates new integrin attachments to the ECM as the
cell front expands to new regions. The rear of the cell detaches, breaking
some of the existing bonds. See Figure 2 for a cartoon that illustrates this
concept. Engagement of integrin receptors then signals to GTPases [45, 46,
47]. Experimental evidence suggests that ECM signaling strongly enhances the
activation of RhoA [48, 49], although it likely also have positive effects on
Rac activity [50, 45].

Figure 2. Side view of cell with new cell-ECM adhesion bonds forming in front
and old ones rupturing in back

###  1.2. Integrin bonds

Integrins cluster into focal adhesions (FA’s), creating local hot-spots of
adhesion. Structural proteins within FAs enable binding of the integrins to
the cytoskeleton, which stabilize the adhesion and allows cells to transmit
cytoskeletal forces to the integrins. The biophysical properties of integrins
have been elucidated in recent years [51]. The binding-unbinding kinetics of
integrin bonds are influenced by mechanical force. For a long time, it was
thought that the detachment rate of integrins increases with force (slip-
bond), but recent experiments revealed that certain integrins [51] actually
behave like catch-bonds. Applying force to a catch-bond causes tightening and
reduces the rate of unbinding up to some force threshold, beyond which the
bonds start to break. Hence, the lifetime of a catch-bond is non-monotonic,
and maximal for some intermediate force magnitude.

Experimental [52] and modeling ([53, 54, 55] and many others) papers have
investigated the relationship between applied force and adhesion bond growth
and breakage. In this work, we followed the ideas of Novikova & Storm [56],
who approximated a focal adhesion as a cluster of catch-slip bonds of which
the cluster size is optimal under finite force (catch) and ruptures above a
certain force threshold (slip). In this paper, we compare this model with a
pure slip-bond model.

In our ultimate deforming cell model, we keep track of parts of the cell that
advances, forming fresh integrin bonds and establishing ECM signaling.
Retraction of parts of the cell result in forces pulling on these adhesion
bonds, and tearing them off, which is also included in our ultimate cell
deformation simulations.

###  1.3. Experimentally observed phenotypes

Much of the work reported here was motivated by the experiments of [42] on
aggressive melanoma cells (1205Lu) adhering to topographical surfaces
mimicking ECM. These substrates were arrays of nano-posts of varying densities
and anisotropy, coated with the adhesion protein fibronectin (FN). Our work
here follows on the foundations built by the theoretical work on [43]. We
briefly summarize the key experimental observations.

The cells in [42] expressed fluorescent tagged pleckstrin homology (PH) domain
of Akt, an indicator of the activity of PI3K, a kinase associated with the
activity of Rac at the “cell front”. This region corresponded with protruding
lamellipodia. The location of the cell front was tracked over time under
various conditions. The dynamics were classified into three main types,
persistent, random, or oscillatory. Persistent cells had a stable “front”, and
tended to be polarized. In oscillatory cells, the PI3K “hot spot” switched
back and forth across the cell, so that front and back polarity exchanged in
some rhythmic process. In the “random cells”, the region of “frontness” jumped
from one site to another along the cell edge, with no clear periodic pattern.

In [42], every experiment produced some fraction of each of the three types,
with relative proportions dependent on the type of manipulation. The adhesion
of the cell and its access to the ECM could be changed by modifying post
density. On sparser posts arrays, cells penetrate between the posts and attach
to the underlying ECM. Higher post arrays (hence lower cell-ECM contact) led
to a greater proportion of random cells and less persistently polarized cells
[42].

The model in [43] could account for major experimental observations in [42]
with a minimal model based on Rac-Rho mutual inhibition, opposed effects of
Rac and Rho on cell expansion (Rac) and contraction (Rho), together with
feedback from the ECM to the GTPases. Our paper proceeds to explore this idea
with greater spatial resolution and with enhanced computational tools.

##  2\. Modeling background

We briefly summarize background information for the GTPase models,
assumptions, and common formulation on which this paper rests. Details for
these models can be found in [57, 58].

Each GTPase acts like a “molecular switch”, where only the active “ON” state,
bound to the plasma membrane, has downstream influence. However, limited
availability of inactive GTPase can affect the rate of activation. On the
timescale of interest (seconds, minutes), the total amount of a given GTPase
is roughly constant in the cell. Thus, for GTPase $j$, a pair of partial
differential equations (PDEs) is used to keep track of the distribution of
active ($G^{j}$) and inactive ($G_{I}^{j}$) forms:

(2.1a) $\displaystyle\frac{dG^{j}}{dt}$
$\displaystyle=A_{j}G_{I}^{j}-I_{j}G^{j}+D^{j}\triangle G^{j},$ (2.1b)
$\displaystyle\frac{d{G_{I}}^{j}}{dt}$
$\displaystyle=-A_{j}G_{I}^{j}+I_{j}G^{j}+D_{I}^{j}\triangle G_{I}^{j}.$ Here
$A_{j},I_{j}$ are rates of GTPase activation and inactivation, that depend on
interactions with other GTPases and the ECM. Inactive GTPase resides in the
cytosol, and diffuses faster than active GTPase, so that $D^{j}\ll D_{I}^{j}$.

We also track the adhesion of the cell to the ECM. A typical equation for this
ECM adhesion variable is,

(2.1c) $\frac{dE}{dt}=\epsilon(a-dE).$

The ECM rate of increase, $a$, and decay $d$ will be important GTPase-
dependent terms describe below, and $\epsilon$ is a small parameter signifying
slower dynamics. The ECM does not diffuse. The set of equations (2.1) for the
case of Rac ($R$) and Rho ($\rho$) lead to a total of five PDES. The specific
assumptions are described next.

###  2.1. Rac-Rho mutual antagonism

For the mutually antagonistic Rac and Rho, we follow the basic assumptions in
[58, 42, 43] to arrive at the GTPase equations

(2.2a) $\displaystyle\frac{dR}{dt}$
$\displaystyle=A_{R}(\rho)R_{I}-I_{R}R+D\triangle
R,\qquad\frac{dR_{I}}{dt}=-A_{R}(\rho)R_{I}+I_{R}R+D_{I}\triangle R_{I},$
(2.2b) $\displaystyle\frac{d\rho}{dt}$
$\displaystyle=A_{\rho}(R,E)\rho_{I}-I_{\rho}\rho+D\triangle\rho,\qquad\frac{d\rho_{I}}{dt}=-A_{\rho}(R,E)\rho_{I}+I_{\rho}\rho+D_{I}\triangle\rho_{I},$
where (2.2c) $A_{R}(\rho)=\frac{b_{R}}{1+\rho^{m}},\quad
A_{\rho}(R,E)=\frac{b_{\rho}(E)}{1+R^{m}},\quad I_{R}=I_{\rho}=\delta.$

The crosstalk of Rac and Rho is represented in the activation rates, whereas
inactivation rates have been taken as constant.

###  2.2. Rac-Rho ECM feedback

We add two-way feedback from ECM to Rho and from Rac and Rho to the ECM by
assuming that

(2.3a) $\displaystyle\frac{dE}{dt}$
$\displaystyle=\epsilon(a(R,\rho,E)-d(R,\rho,E)E),$ (2.3b)
$\displaystyle\mbox{and }\quad b_{\rho}(E)$
$\displaystyle=k_{E}+\gamma_{E}\frac{E^{m}}{(E_{0}^{m}+E^{m})}.$

The Rho activation rate $b_{\rho}(E)$ is affected by ECM signaling [48, 49].
In (2.3)b, $k_{E}$, is a basal Rho activation rate, and $\gamma_{E}$ is a
parameter governing the strength of feedback from ECM signals to Rho
activation.

The basic assumptions made in Eqs. (2.2c) and (2.3)b will be a common basis
for all models discussed in this paper. Details of the assumptions for the ECM
equation (2.3)a will be discussed in what follows.

##  3\. Model variants

The three model variants we consider differ only in the way that the ECM
reacts to the downstream effects of Rac and Rho. In all three cases, we will
see that Rac feedback is positive, and Rho feedback is negative, but the
latter two models are formulated with a greater detail about integrin bond
formation and breakage.

###  3.1. Model I: the original version

In the original model [42, 43], it was reasoned that cell-ECM contact
increases when Rac causes the cell to spread, and decreases when Rho causes
cell contraction. This idea was modeled with a assumptions for (2.1c) that:

(3.1) $a=a(R)=K+\gamma_{R}\frac{R^{n}}{(R_{0}^{n}+R^{n})},\quad
d=d(\rho)=k_{P}+\gamma_{\rho}\frac{\rho^{n}}{(\rho_{0}^{n}+\rho^{n})}.$

We refer to this version as Model I.

###  3.2. Integrin bonds and ECM dynamics

In the integrin bond model variants, we take into account a force $F(R,\rho)$
exerted on adhesion bonds with basic form

(3.2a) $F(R,\rho)=\beta_{\rho}\frac{\rho}{1+\rho}-\beta_{R}\frac{R}{1+R},$ and
a lower bound of 0. Rho-driven contraction pulls on adhesion bonds. When Rac
dominates and leads to local protrusion, there is no force on the bonds, given
that the cell “rolls over” those adherent sites. The ECM equation is then
taken to be (3.2b) $\frac{dE}{dt}=\epsilon\left(K(E_{t}-E)-d(F,E)E\right).$

So, $a=a(E)=K(E_{t}-E)$ and $d(F,E)$ is the force-dependent integrin bond
breakage rate, where the amount of force depends on Rac and Rho. In the
absence of force, the ECM variable settles into its steady state level,
$E_{t}$, on a timescale of $1/(\epsilon K)$. We can interpret $E_{t}$ as the
maximal local density of bound integrin bonds. Rac and Rho activity creates
force that affects the bond breakage, and impacts ECM signaling. To arrive at
reasonable assumptions about $d(F,E)$ we used the integrin biophysics model of
[56].

Figure 3. Bifurcation diagrams for (3.2b) showing the steady state level of
adherent integrin bonds, $E$, as a function of applied force, $F$. Solid lines
indicate the stable branches and dotted lines are the unstable branches.
$E_{0}$ is the level of adhesion at which Rho activation turns on, see (2.3).
(A) Steady-state $E$ for slip-bond case, with $d(F,E)$ as in (3.3). Higher
force leads to lower adhesion. (B) as in (A) but for the catch-slip-bond,
$d(F,E)$ given by (3.4). Higher force leads to higher adhesion, until the
force reaches a maximum and adhesions break down completely. In each case, two
values of the maximal adhesion $E_{t}$ are shown.

####  3.2.1. Model II: Slip bond

The rate of bond dissociation of slip-bonds is well-approximated by

(3.3) $d(F,E)=k_{0}\exp\left(\frac{F}{p(E+E_{s})}\right).$

The total force is distributed over all local integrin bonds at a given
adhesion site, and $F/E$ is a force per bond that leads to rupture. The small
correction $E_{s}$ in the denominator prevent blowup as $E\to 0$, ensuring
that small bonds under low force can grow.The parameter $p$ is a reference
value of force per bond, making the term in large braces dimensionless. When
$F=0$, the bond breakage rate is $d=k_{0}$, whereas when $F/E\approx p$, the
breakage rate is $d\approx k_{0}e^{1}\approx 2.7k_{0}$. Hence $p$ sets the
reference force per bond at which adhesions detach at 2.7 times their basal
rate of detachment. Larger $p$ means that greater force per bond is needed to
cause the same bond breakage rate.

####  3.2.2. Model III: Catch-slip bond

In this case, we take

(3.4)
$d=k_{0}\exp\left(\frac{F}{p(E+E_{s})}\right)+k_{0c}\exp\left(-\frac{F}{p(E+E_{s})}\right).$

Catch-bonds tend to grow stronger under the influence of force up to some
limit, so that their lifetime is maximal under some optimal force. Beyond that
point, for larger applied force, the bonds break. We adopt the values
$k_{0}\approx 0.0004,k_{0c}\approx 55$, appropriate to $\alpha 5-\beta 1$
integrin [56].

####  3.2.3. ECM-force bifurcation plot

To gain some intuition, we compared the force-dependent adhesion for catch-
slip versus slip-bonds. To do so, we isolated Eqn. (3.2b), together with
either of the two force-dependent bond dissociation rates (3.3) or (3.4)
(leaving out the Rac-Rho dependence of force). We plotted the steady state
adhesion $E$ as a function of force $F$ in Figure 3. We also compare the value
$E_{0}$ (level of adhesion at which ECM-signaling to Rho turns on in Eqn.
(2.3)) on the same plots. Figure 3 shows that in the case of a slip-bond, the
adhesion cluster decreases as force increases and quickly ruptures above some
force threshold, implying that (in the full model) signaling to Rho is only
sustained at low force. In the case of catch-slip bonds, the adhesion clusters
increase as force increase up to some threshold for rupture. In that case,
signaling to Rho is most strongest at intermediate force magnitude. When there
is more ECM contact (compare $E_{t}=2000$ with $E_{t}=1000$), the clusters are
larger and more long-lived, as more force is required for rupturing the bonds.

###  3.3. Model summary

Briefly, all three models use the basic Rac-Rho equations (2.2) and (2.3),
with a generic ECM dynamics. Model II assumes slip bond dynamics, and Model
III takes slip-catch bond dynamics for the ECM equation. The three variants
are summarized in Table 1.

Model Type | Description | Equations | Ref
---|---|---|---
Original (I) | Generic ECM eqn. | (2.2), (2.3), (3.1) | [42, 43]
Integrin bonds (II) | Slip-bond | (2.2),(2.3)b,(3.2), (3.3) | [56]
Integrin bonds (III) | Slip-catch bond | (2.2),(2.3)b,(3.2), (3.4) | [56]
Table 1. Summary of the three models considered in this paper. All models use
the basic set of Rac-Rho equations, (2.2), but the specific assumptions about
the ECM dynamics differ, as indicated in this table.

##  4\. Preliminary results: spatially uniform case

There are several submodels and simplified cases that are of interest. We
briefly survey some previous results and link to the fuller spatial models.

###  4.1. Well-mixed Rac-Rho model

A simple reduced variant of the model for only Rac and Rho (with no ECM,
$b_{\rho}$= constant, and no spatial terms) is the well-mixed version of (2.2)
where spatial gradients are ignored. In that case, it is common to assume that
the total GTPase of each type, $G_{T}=R_{T},\rho_{T}$ is constant in the cell
on the timescale of interest, so that

$R_{T}=R+R_{I},\quad\rho_{T}=\rho+\rho_{I}.$

Hence, the inactive forms can be eliminated from the system, resulting in a
pair of ordinary differential equations (ODEs), (11.1). This mutual-antagonism
system is is bistable for some range of parameter values [58], as shown in the
left panel of Figure 4. Bistability implies hysteresis, and slow negative
feedback from some influence can lead to an excursion around the hysteresis
loop that results in oscillations. The ECM in the Rac-Rho-ECM model plays the
role of this additional feedback.

Figure 4. Left: The well-mixed Rac-Rho system (11.1) on its own is bistable,
leading to an S-shaped bifurcation diagram with respect to a parameter such as
$b_{R}$ or $b_{\rho}$ (not shown). Slow negative feedback that shifts a
parameter (for example, the effect of the ECM on $b_{\rho}$) can, in
principle, lead to cycling around the hysteresis loop, and emergence of limit
cycle oscillations. Right: The well-mixed Rac-Rho-ECM (Model I) bifurcation
diagram. Some of the Rac-Rho bifurcation structure is preserved, but we see a
Hopf bifurcation at some intermediate value of the parameter $b_{R}$. There is
a very narrow range of $b_{R}$ where two stable limit cycles coexist, as shown
in the two zoom insets.

###  4.2. Well-mixed Rac-Rho-ECM model (I)

Once ECM is coupled to the Rac-Rho system in for Model I, interesting dynamics
emerge. For example, the well-mixed variant of the Rac-Rho-ECM Model I, whose
equations are given by (11.2) has a regime of oscillations sandwiched between
states of high and low Rac. A typical bifurcation diagram for (11.2) is shown
on the right panel of Figure 4.

Figure 5. The target volume of a CPM cell is linked to its internal Rac-Rho
GTPase levels for a well-mixed Model I Rac-Rho-ECM signaling. The cell
contracts and expands accordingly in this oscillatory regime. Top left:
schematic diagram. Top right: the time dependent variables in the cell. Bottom
row, left to right, the level of Rac (on the color scale at right) in cell
shape snapshots at $t=20,24,26,28,32,34,46$. Produced with the open source
software Morpheus [59] using the file OneLamelShape.xml.

The well-mixed model is a coarse approximation of the cell signaling, but we
can visualize its basic implications by linking the GTPase state to a cell’s
“contact area”. To do so, we simulated the well mixed Model I, Eqs. (11.2), in
a cellular Potts Morpheus simulation, and assumed that the cell target area is
high when Rac dominates and low when Rho dominates (Schematic, top left). A
representative example of results is shown in panels of Figure 5. We see the
relaxation oscillations typical of such systems (top right panel), and the
fluctuations in cell contact area as well as Rac activity level over once
cycle in the time sequence of shapes. (Rho activity peaks whenever Rac
activity drops.)

In the paper [42, 43], the authors considered a well-mixed model with two
lamellipods, rather than a single compartment as we show in Figure 5. Each
lamellipod had an internal Rac-Rho-ECM model (similar, though not identical to
our Model I). The two were coupled by a constant total pool of Rac and of Rho,
and additional terms for competition of the lamellipod size were added. This
model could then show regimes of persistent polarity, oscillations, as well as
oscillations. A similar example is included in the Appendix Figure 16.
Oscillations and Hopf bifurcations are also possible in the well-mixed
variants of Models II and III.

##  5\. Setup for the spatially distributed models

We next consider the full spatially distributed (PDE) models. Now there are
five PDEs in each of the model variants, to include the PDEs for inactive Rac
and Rho. These cannot be eliminated as they diffuse at rates different from
the active forms. The full equations are presented in the Appendix for each
case.

###  5.1. Spatial geometries

We consider three distinct geometries, as shown in the schematic, Figure 1.
(a) A 1D domain, as in Figure 1b. This is interpreted as a transect across the
diameter of the cell, with ends at opposite cell edges. This 1D model can also
be interpreted as the geometry of a cell that is confined to a narrow channel
as in [60]or moving along thin “ECM” fibers or nano-wires, as in [61, 62, 63].
(b) A static 2D circular domain, as illustrated in Figure 1c. This represents
a cell that is “frozen” by an inhibitor of the cytoskeleton, or stuck to an
adhesive island. In this case, the signaling can take place, but there is no
net motion or change of shape. (c) A fully deforming 2D cell, as shown in
Figure 1d. In all these cases, we assume Neumann boundary conditions, since no
proteins leak out of the cell edges.

###  5.2. Kymographs

In each case, and for each of the model variants I, II, III, we represent
results by a kymograph, as shown in Figure 1e, with time increasing up the
vertical axis, and position along the horizontal axis. In 1D, the position is
distance along the cell diameter, $0<x<L$. For the 2D cells, the model is
solved on a 2D domain but results are summarized by similar kymographs of
activity along the cell edge. The“position” is taken along the (normalized)
circumference of the shape, so that $-\pi<x<\pi$ around the cell perimeter.

Excitable systems that are spatially distributed are known to sustain standing
waves, target waves and spiral waves (2D) as well as waves that reflect from
domain boundaries. A single standing wave that has high Rac at one domain
boundary and high Rho at the opposite end corresponds to a polarized cell.
Spiral waves can be recognized by the rotation of the high Rac activity along
the domain edge, which in our kymographs, appears as a series of bands where
the slope of the band $\Delta x/\Delta t$ is the wave speed along the
boundary. Waves that reflect from one edge to the opposite edge show up as a
set of peaks and troughs. In all cases, we display the Rac activity on a color
scale from low (black or purple) to high (bright orange and yellow). Rho
activity (not shown) is the reciprocal; high Rac implies low Rho and vice
versa. See Figure 15 for an example where all variables are displayed.

###  5.3. Parameter sweeps

In correspondence with Park et al [42], we explore how model behaviour changes
when $\gamma_{E}$ (rate of ECM-feedback to Rho activation) and $b_{R}$ (basal
rate of Rac activation) are varied. We also vary the constant $K$ in Model I
(basal rate of ECM signaling increase), and the corresponding constants
$E_{t}$ in Models II and III (the maximal adhesion size).

Our results, shown as arrays of kymographs for each model, can be compared to
the 2-parameter bifurcation diagrams in Figure 5d in [43] (for “Hybrid Model
3” in that paper).

##  6\. Results: 1D spatial domain

###  6.1. Model I in 1D

The spatially distributed Model I is given by the system (11.3). Results shown
in Figure 6, consist of multiple kymographs of $R(x,t)$ (with the same color
and axes convention.) These results agree, broadly speaking, with previous
work (“Hybrid” model in [42] Fig 2c in and [43] Fig 5d, Model 3). For a fixed
value of $K$, we find regimes of low Rac (along the low $b_{R}$, vertical
region of the parameter plane), of high Rac (along the low $\gamma_{E}$,
horizontal region of the parameter plane). Wedged between these are regimes of
front-back oscillation and/or persistent polarization. Increasing $\gamma_{E}$
at intermediate values of $b_{R}$ results in higher frequencies of
oscillation. As $K$ increases, greater swaths to the parameter plane
correspond to polarized states, and the wedge of the oscillatory regime
shrinks.

Figure 6. Model I regimes: Behaviour of Rac for the full 1D spatially
distributed version of Model I, Eqs. (11.3), showing parameter sweeps for the
basal Rac activation rate $b_{R}$ and the ECM-induced Rho activation rate
$\gamma_{E}$. Kymographs oriented as in Fig 1e, with $0\leq x\leq 3$
(horizontally) and $0<t<1000$ (vertically upwards). The parameters were varied
as follows: $1\leq b_{R}\leq 7$ (in steps of 1.0), $\gamma_{E}=1,2,4,6,8$,
$K=0.1,0.5,1.5,1.75,2,2.5$. Other parameter values:
$R_{T}=3.0,\rho_{T}=6.0,\delta=1.0,k_{E}=2,\epsilon=0.001,E_{0}=1.5,n=3,\gamma_{\rho}=10,\gamma_{R}=5,k_{\rho}=0.45,\rho_{0}=2.4,R_{0}=1.0$.
Heatmap is purple for low and yellow for high Rac activity levels. Approximate
regime borders (green) added manually. Produced with Morpheus file
Model1RacRhoECMPDEsIn1D.xml

###  6.2. Model II: Slipbond model in 1D

We carried out a similar parameter sweep for the spatially distributed Model
II, Eqs. (11.3) but with ECM as in (11.4) and $d$ given by (3.3) for slip-bond
adhesion sites. Results are shown in Figure 7, but with variation in the
maximal adhesion size, $E_{t}$ between panels. If the maximal size of
adhesions is too small, $E_{t}\approx 100$, no interesting dynamics is
observed, and the domain remains spatially uniform. For $E_{t}\geq 500$, we
see the onset of oscillations. Higher values of $b_{R}$ delays the onset of
the cycles, and increasing $\gamma_{E}$ leads to period doublings. We see no
polar patterns, as yet, but these emerge once maximal adhesion size is tuned
to $E_{t}\geq 2000$, and are even more evident by $E_{t}=3000$. Finally, by
$E_{t}=5000$, the oscillations are no longer present in the analogous regimes,
and the cell is either homogeneous or polarized. Interestingly, a few
kymographs demonstrate a polarity reversal in one or another run.

Overall, the results of Model II share similarities with results of Model I,
with a few subtle differences. There are similar “wedge-shaped” regimes. The
behaviour near the low $b_{R}$ or low $\gamma_{E}$ regimes are still uniform.
There is a similar tendency to settle into a polar state when the adhesions
become more dominant (increasing $K$ in Model I versus increasing $E_{t}$ for
Model II). As before, the frequency of cycles increases with $\gamma_{E}$, and
decreases with $b_{R}$.

Figure 7. Model II regimes: As in Fig 6 but for the slip-bond integrin model
(Model II). The parameters were varied as follows: $1\leq b_{R}\leq 7$ (in
steps of 1.0), $\gamma_{E}=1,2,4,6,8$, max adhesion size, $E_{t}$: varied from
$100,500,1000,2000,3000,5000$. Other parameter values:
$R_{T}=3.0,\rho_{T}=6.0,\delta=1.0,k_{E}=2,\epsilon=0.001,E_{0}=300,n=3,\rho_{0}=2.4,R_{0}=1.0,\gamma_{E}=4,K=10,E_{T}=1000,k_{0}=5,E_{s}=100,p_{s}=0.35,\beta_{R}=1200,\beta_{\rho}=1600$.
Produced by Morpheus xml file Model2SlipBondRacRhoECMPDEsIn1D.xml

###  6.3. Model III: Catch-Slip bond model in 1D

Results for Model III are shown in Figure 8. Overall, we observed a
substantial decrease in the size of the oscillatory parameter regime. The runs
tended to produce polarized patterns of Rac activity, occasionally after a few
transient reversals. We found that the oscillatory phenotype only emerges when
the back of the cell is in the slip-regime of force, meaning that all adhesion
bonds are broken. Only then it is possible for Rac to be sufficiently
activated at the back for the front to switch.

###  6.4. Significance

Conclusions from the spatial 1D results are as follows. (1) Such models have
overall concurrence with the simplified 2-compartment models in [42, 43], but
allow for finer detail on both the dynamics and the spatial distributions. (2)
In each of the cases tested, the ECM variable is responsible for emergence of
oscillatory behaviour in an otherwise polarizable, multistable system [58].
(3) The details of the ECM mechanism, and the ECM equation affect the breadths
of some parameter regimes, but otherwise do not significantly change the
overall conclusions. The three models tested all share regimes of uniform,
polar, and oscillatory behaviours in similar swaths of parameter space. What
is more important is the feedback from Rac to amplifying ECM, from Rho to
depressing ECM, and from ECM to activating Rho.

In Figure 9 we summarize the mechanisms by which adhesion can cause
oscillations or persistence in model II and III. The main idea is that the
right level of adhesion is necessary for oscillations to take place. Figure 9A
shows profiles for Rac, Rho and force $F$ in a 1D polarized cell. Whether the
cell polarization persists or exhibits a front-back flip depends on the level
of adhesion, which is itself determined by the magnitude of the force and type
of adhesion bonds. In general, if the adhesion at the front is higher than
$E_{0}$, it signals to activate Rho. If at the same time, the adhesion at the
back is lower than $E_{0}$, Rho is not sufficiently active to inhibit Rac,
allowing Rac activity to invade the rear of the cell, allowing a front-back
polarity flip. If however adhesion at the back exceeds $E_{0}$, Rho is
upregulated, which inhibits Rac from moving to the rear, maintaining a
persistent cell polarity.

The slip-bond and catch-slip bond models lead to different levels of adhesion
(Figure 9), but both have regimes consistent with oscillations. Referring to
Figure 3 helps to understand the basis for this behaviour. In case of slip-
bonds, adhesion is highest at the front and lowest at the back, (Figure 9B,D),
as adhesion decrease with force (Figure 3A). If the force at the back of the
cell is above the rupture threshold (Figure 3A), oscillations can occur
(Figure 9B). If the force at the back of the cell is too low, polarity is
persistent (Figure 9D).

For catch-slip bonds, adhesion is optimal under higher forces (Figure 3B), and
hence, maximal at some distance from the cell front, either near the middle
(Figure 9C), or at the back (Figure 9E). The conditions for
oscillations/persistence for the catch-slip bond are similar to the slip-bond,
but with different force thresholds (Figure 3B).

Figure 8. Model III: the catch-slip bond model. As in Figs. 6 and 7 but for
the case where $E$ has the catch-slip bond representation. Similar trends are
observed, but for a smaller oscillatory regime. Parameters as in Fig 7 with
the following exceptions:
$p_{s}=0.08,\beta_{R}=1000,k_{0}=\exp(-7.78),k_{0c}=\exp(4.02)$. Produced with
Morpheus xml file Model3CatchSlipBondRacRhoECMPDEsIn1D.xml Figure 9. Schematic
explanation of how the slip vs catch-slip integrin dynamics lead to
oscillations or persistence. (A) Typical Rac/Rho and force profiles in a 1D
static cell. Highest forces are at the back while low forces are at the front.
Increasing or decreasing the total force magnitude leads to high ($>E_{0}$) vs
low ($<E_{0}$) adhesion, as shown in Figure 3. (B-C) Adhesion profiles (light
blue curve) consistent with oscillations, or (D-E) persistence. (B) Under high
force, slip-bonds break at the back of the cell and grow at the front. This
upregulates Rho at the front, leading to Rac inhibition. Hence, Rac moves to
the back. The pattern repeats, resulting in front-back oscillation. (C) As in
(B), but for the catch-bond model: here adhesion is elevated at some distance
away from the front. (D) If the force magnitude is too low, slip-bonds grow at
the back of the cell. This reinforces existing Rho activity at the back,
causing the polar Rac/Rho profiles to persist. (E) In case of catch-bonds,
intermediate force at the back causes maximal adhesion at the back, which as
in (D) upregulates Rho and reinforces persistence.

##  7\. Results: 2D domains

From here on, we drop the Original Model I, as it is phenomenological, and
concentrate on a comparison of Models II and III. We also restrict attention
to the more interesting parts of the parameter regimes.

###  7.1. A circular 2D static domain

Figure 10. Models II (left) and III (right) in a 2D static circular domain.
Kymographs show Rac activity along the rim of the circle. Parameter values and
colour scheme as in previous figures.

We simulate the model equations in a fixed circular domain with no flux
(Neumann) boundary conditions and initial conditions as described in the
Appendix. The kymographs show dynamics along the circular perimeter. The full
2D patterns consisted of spiral, oscillating or static standing waves.

We carried out parameter sweeps for three values of $E_{t}$ and we added a few
larger values for $\gamma_{E}$. We investigate whether the circular geometry
in 2D introduces new phenotypes. Results are shown in Figure 10 for Model II
(left column) and Model III (right column) for three values of the adhesion
size $E_{t}$.

Overall, the behaviour of both Model II and III in 2D shares certain features
of the analogous 1D results: regimes of low Rac, high Rac, oscillatory, or
persistent polar patterns are evident as before. The catch-slip model (III)
still promotes more polarized persistent behaviour than the slip-bond (II)
version, as before. However, there are several notable differences. (1) for
$E_{t}=1000$, both models, but especially Model III produce a more exotic
variety of patterns, including standing oscillations and zig-zag waves, along
with random-looking transition patterns. (2) In place of back and forth 1D
cycles, we see traveling waves of Rac activity that get wider as $E_{t}$
increases. These represent a rotation of Rac along the rim of the circular
domain, as would occur when a spiral-wave pattern forms inside. (3) In most
cases, back-front oscillatory patterns are transient and evolve into a single
spiral around the cell edge. For higher values of $\gamma_{E}$, the spirals
sometimes change direction from counter-clockwise to clockwise, or vice versa.
Spirals are more likely for higher values of $b_{R}$, the basal Rac activation
and lower values of ECM-Rho activation ($\gamma_{E}$). Increasing $b_{R}$ and
decreasing $\gamma_{E}$ both lead to higher Rac and thus wider Rac fronts.

We find that parameters that resulted in a persistent phenotype in 1D can
become spirals in 2D, most notable in the slip-bond model. We also find
intermediate phenotypes, where spirals are followed by a few transient front-
back oscillations, that then evolve into a spiral again. Without carefully
observing the pattern, this behavior looked seemingly random at first.

Finally, we note that the boundary between regimes is less sharp in 2D than in
1D results. For example, in Model III for $E_{t}=2000$, the oscillatory and
persistent regimes appear to bleed into one another. This may stem from
coexistence of multiple steady states (migratory phenotypes) in a given
transition zone in parameter space. This kind of coexistence has already been
discussed in a simpler model of Rac and Rho alone in [58] where patterns were
far simpler (uniform low, high, or polar). The 2D geometry appears to
accentuate this possibility.

In conclusion, the 2D geometry gives rise to more complex patterns, such as
spirals. In a motile cell, spirals could be indicative of circular motion. Or,
cell shape changes may transform spiraling into front-back oscillations or
stabilize the spiral into a persistent front. Consequently, we asked how the
patterning changes as we allow the 2D domain to deform.

###  7.2. A 2D deforming domain

Models II and III, are now solved over a similar range of parameters, but in a
2D deforming domain, representing the top-down view of the cell shape.
Briefly, domain deformation was tracked using a custom-built cellular Potts
model (CPM) calculation with an in-built PDE solver, as described in the
Appendix. (As of this writing, open-source packagew such as Morpheus [59] and
CompuCell3D [64] do not yet have a PDE solver for a deforming CPM cell.)

To link GTPase activity to the evolving cell shape, we assumed that high Rac
activity at the cell edge promotes local outwards protrusion of the edge (see,
e.g. [65]), whereas high levels of Rho lead to inwards contraction as in [66],
bypassing the explicit representations of actin, myosin, and other cellular
components that were included in previous work [26]. In distinction with the
coarse-grained results of Figure 5 where the entire cell spreads or contracts,
here the spatial distributions of the GTPases is resolved in detail. Hence
protrusion and contraction of the cell edge are localized to spots on the
perimeter where there is high Rac or high Rho, and we can track both cell
polarization, change of overall shape, and motility. Details of the custom
built CPM computation are provided in the Appendix.

Figure 11. As in Figure 10 but for a 2D deforming cell. Parameters as in Table
3. The same two models and parameter sweeps are shown. The shapes of three
sample cells (circled above) are shown in Fig. 12. Movies can be found here
https://imgur.com/a/Btw4z9H Figure 12. Sample cell shapes and Rac activity in
a brief time sequence. The three parameter settings correspond to circles in
Fig 11. (A) Persistent cell. (B) Oscillatory cell. (C) Random cell
protrusions. The kymographs (right) show the long-time behaviour for each
case, and are magnified views of the circled images in Fig 11. Note that the
expanding perimeter of the oscillatory cell is seen in the kymograph for (B).
Movies can be found here https://imgur.com/a/Btw4z9H

Using the same parameter sweeps as before for Models II and III, we show
results in Fig 11, with selected sample cell-shape time sequences in Fig. 12.
As before, kymographs track the Rac activity level along the cell edge (dotted
white curve in Figure 1d). Because the domain deforms, and its perimeter
length fluctuates slightly with time, the right edge of each kymograph in
panels of Fig 11 appear slightly ragged.

By comparing patterns in the 2D deforming domain (Figure 11) to the identical
conditions in the 2D static domain (Fig. 10), we can identify the effect that
domain dynamics has on internal patterns. Regimes of spiral waves in a static
domain tend to become front-to-back oscillations in the deforming domain. This
is evident in the Slip-bond model (II) for $E_{t}=2000,3000$ where the slanted
yellow-purple bands that represent spiral waves in the static case (Fig. 10)
are replaced by the “honeycomb” kymograph patterns that depict front-back
oscillations (Fig. 11). For the slip-bond with $E_{t}=3000$, a larger regime
of polarization emerges in band at low $\gamma_{E}$ that was previously only
spiral waves. The cell deformation appears to dampen spirals in favour of
either front-to-back cycles or stable polarity.

The slip-bond model produces much more regular oscillations than the catch-
bond model. It is easy to note that the frequency of oscillation is reduced at
higher ECM $E_{t}$ (compare $E_{t}=2000$ with $E_{t}=3000$). On the other
hand, the ECM-Rho feedback strength ($\gamma_{E}$) increases the frequency.
So, the feedback strength and ECM binding rate, both of which regulate the
amount of signaling, have opposite effects on the oscillation frequency. The
explanation is that for higher values of $E_{t}$, the cells need to apply more
force to break the bonds at the back of the cell to allow a front-back switch
to occur. Higher values of $\gamma_{E}$ makes the cell more sensitive to weak
signals from the ECM, allowing for quick oscillations.

For the catch-slip bond model (III), at $E_{t}=1000$, results in static and
deforming cell domains are somewhat similar, except for the corner at low
values of both $b_{R}$ and $\gamma_{E}$, where spirals (in the static domain)
turn into cycles (deforming domain). Polar patterns still tend to dominate at
$E_{t}=2000,3000$, but we find an increased swath of irregular and random
oscillations for $E_{t}=2000$ when the cell can deform.

We selected three examples of interesting dynamics (indicated by circles in
Fig. 11) to track in the full 2D simulations. These include (A) a persistent
polar case, (B) front-back oscillations, and (C) random activity pattern. The
cell shape dynamics are shown for a short time sequence in Fig. 12, together
with a more detailed view of the corresponding (long-term) kymographs, to
demonstrate the overall outcomes. In (A) we see a cell that polarizes, and
starts to migrate directionally. The polarity persists over time. For (B),
high Rac activity oscillates between two cell ends. This causes the cell to
elongate. Hence, the perimeter of the cell also increases with time, as seen
in the gradually expanding envelope of the kymograph. This cell fails to have
significant net migration. In (C), the zone of active Rac continues to move
around the cell irregularly, so protrusion/retraction is undirected, and the
cell fails to migrate. The kymograph demonstrates many interpenetrating waves
along the cell edge, without clear periodicity.

We can understand some of the results of Fig. 11 from the actual cell shapes.
First, for a cell with initial spiral internal dynamics, the Rac-induced edge
protrusion can result in slight elongation of the cell. This tend to break
symmetry, damping the spirals that arrive at the domain boundary. An elongated
domain then favours the oscillations that dominated the previous 1D results,
and the oscillations, in turn make the cell longer and contribute to self-
amplification. Cell edge deformation can also push spirals into a single
front, entering the persistent regime. This is visible in the slip-bond model
with $E_{t}=3000$, where many patterns were spirals in the static 2D domain,
but have become persistent in the deforming domain. Most notably the catch-
bond model in combination with shape deformations causes irregular
oscillations. In some cases, these oscillations ultimately evolve into a
single persistent front.

In conclusion, the 2D deformable simulations show that cell motility can
change the phenotypical outcome expected from static 2D simulations. This can
explain the huge variability in the three phenotypes of the same cell types on
the same kind of ECM, something observed also in [42].

##  8\. Results: effect of cell stiffness on patterning dynamics

Finally, we asked how the biophysics of the cell and the properties of the
domain affect the internal patterning of the signaling activity. To do so, we
reran parameter sweeps for variants of the CPM simulations with distinct CPM
parameters. (1) We decreased the feedback from Rac/Rho to cell edge
protrusion/retraction. This manipulation would be related to the
responsiveness of actin polymerization at the front (activation of WASp and
availability of Arp2/3 downstream of Rac) or the downstream recruitment of
ROCK and myosin activation downstream of Rho). (2) We also changed the
stiffness of the cell. For this purpose, we modified the relative values of
the CPM area and perimeter constraint and cell-medium interfacial energy [67].
(3) Some cell types have far more rapid signaling responses than others. For
example, GTPase responses in neutrophils occur on a timescale of seconds,
whereas in fibroblasts, HeLa cells, and melanoma, the timescale is far longer
(many minutes, hours). To investigate this, we reduced the number of PDE
iterations per MCS by a factor of 2. This is equivalent to slowing down the
signaling kinetics relative to cell motion.

Results of these tests, shown in Figure 13 demonstrate that biophysical
parameters such as cell stiffness greatly affect internal reaction-diffusion
patterning. We highlight a few of these simulations by small panels labelled
(A) to (J) in Figure 13. We demonstrate the full cell-shape dynamics of
several typical cases in Figure 14 (See also video links in the caption of
Figure 13.)

As before, we observe persistence (A) and very regular oscillations (F). But,
we now also observe patterning in what were formerly uniform Rac parameter
regimes. In addition, we find various interesting intermediate phenotypes. For
instance, we observe a cell that has irregular oscillations (B,G,H,I), where
the period of the front-back oscillations vary. Some cells exhibit a kind of
stick-slip motion (J), where Rac transiently appears at the back of the cell
and briefly halts the forward motion, but the cell remains overall persistent.
We also observe phenotypes that oscillate and spiral simultaneously or
consecutively (D,E). In these, transient spirals and oscillations briefly push
the cell edge outwards.

Among the phenotypes, we see varying degrees of periodicity and spatial
regularity. For example, in Figure 15, we show a time series of a persistently
moving cell (labeled (I) in Figure 13) that suddenly flips its front. We show
some time-steps right before the front-back switch and visualize Rac, Rho and
the ECM field. As explained in Figure 9, a front-back switch can occur
provided adhesion is sufficiently high at the front and, simultaneously,
sufficiently low at the back. This is seen in the 2D moving cell in Figure 13
and its accompanying movie. As the cell moves, it creates adhesions at the
front (150 MCS). At the back of the cell, forces are high due to the presence
of Rho. Due to the catch-bond dynamics, these forces stabilize the adhesion at
the back (250 MCS). However, at some point, the forces breach the threshold
that breaks adhesions at the back (350 MCS). Then, the ECM-induced Rho
activation at the front results in repolarization (450 MCS). Interestingly,
some mature adhesions are very long-lived, even at the back of the cell
($t\geq$350 MCS). Such repolarization takes place periodically, every time
there is sufficient adhesion breakage at the rear of the cell. The precise
combination of GTPase dynamics and biophysical cell parameters, determines how
often and where adhesions tend to break, and thus when and where the Rac front
will appear. This results in the various of phenotypes we observe.

Figure 13. Rac signaling in 2D deforming cell with distinct biophysical
properties. Parameters as in Table 3, except with
$n_{i}=400,N=2000,\lambda_{A}=1.5,\lambda_{\epsilon}=10,P_{0}=400,\lambda_{P}=3,J=20000,\beta_{G}=0.25$.
Movies can be found here https://imgur.com/a/Btw4z9H Figure 14. Three
phenotypes from the circled examples in Fig 13. Movies can be found here
https://imgur.com/a/Btw4z9H Figure 15. Rac, Rho signaling and ECM adhesion
during cell motility of example (I) from Figure 13. Movie can be found here
https://imgur.com/a/p8zxH64

##  9\. Discussion

In this paper, we formulated a spatially-distributed set of models for
signaling and ECM feedback, motivated by the experiment-theory work of [42,
43]. We added more detail in several ways: (1) We considered a fully spatio-
temporal (PDE) model in place of the original 2-compartment models in [42,
43]. (2) We revised a generic ECM-feedback model to one based on the
biophysics of integrin adhesion bonds, allowing for more direct experimental
validation. (3) We investigated how geometry in 1D, 2D, and a deforming 2D
cell affect the GTPase patterning dynamics. For the deforming cell
simulations, we used the Cellular Potts formalism. Previous work of this type
includes [26] (a Cdc42-Rac-Rho model with F-actin barbed ends) and [66] (a
single-GTPase toy model with F-actin feedback).

By explicitly formulating the roles of local cell forces, integrin dynamics,
and cell shape changes, we could predict observed phenotypes, new phenotypes,
and phenotypical switches. While many details of our models differ from those
of [42, 43], we find similar basic regimes of persistent, oscillatory, and
random behaviour, together with new intermediate behaviours that were not
previously apparent. We examined the models in a hierarchy of well-mixed, 1D,
2D, and deforming 2D. This helps to build up some basic understanding before
exploring the full complexity of the spatial patterning.

As previously observed, Rac-Rho mutual antagonism results in bistability [57],
and feedback from other influences (here the ECM) allows for oscillations and
other dynamics in single cells [68, 43, 69] or in a cell sheet [70]. The
precise nature of the feedback has a limited influence, as we observed in our
three model variants. The key idea is that the relative Rac/Rho activity
levels determine the magnitude of contact with the ECM, whereas the slow
feedback from the ECM signaling leads the system “around a hysteresis loop” by
gradually tuning the Rac-Rho competition. As described in [58, 71], there are
regimes in which several types of steady state patterns coexist (e.g. polar
and oscillatory, in our case). This shows that cells can have multiple
repertoires of behaviour, even when their parameters are the same. A single
“genotype”, that corresponds to a given set of kinetic rates, can give rise to
multiple phenotypes.

We found that geometry has significant impact on the dynamics of patterns. For
example, in 2D domains, we see regular and irregular internal waves such as
fronts and spiral waves. Spirals have been intensively studied in pattern-
forming chemical and biological systems since the 1970s [72, 73]. Analysis of
such patterns by various mathematical, geometric, and physics-based methods
[74, 75, 76] is typically restricted to static domains. It is well-known that
bistable reaction-diffusion PDEs with slow negative feedback can give rise to
traveling waves [77], and that addition of a conservation condition can cause
those waves to stall [78]. Presence of spiral-waves inside cells have been
connected to feedback between actin and other components [79, 80, 81, 82, 68].
A survey of traveling waves in actin dynamics more generally can be found in
[83].

Our most interesting observation is the interplay between a deforming domain
and the internal reaction-diffusion dynamics. The regimes of behaviour in
static 2D domains shift when the patterns interact with an evolving domain
boundary. To some extent, this is expected from previous models of static and
deforming 2D cells. For example, static and deforming domains are compared in
[84] for several basic models of a single GTPase coupled to negative feedback
from F-actin. Dynamic patterns of a Rho network on a static cell shape were
modeled and analyzed mathematically in [17]. In [85], a circuit of Cdc42, Rac,
Rho, and phosphoinositides was simulated on static cell shapes; it was shown
that certain shapes can reverse polarization gradients. In one of the original
fully deforming domain simulations [86], it was shown that domain deformation
can accelerate internal RD dynamics of these intermediates. Neumann boundary
conditions used to depict sealed cell edges set constraints on the level
curves of internal concentrations: those level curves must be orthogonal to
the boundary, so small boundary changes dramatically affect the internal
distributions. A full mathematical exploration of PDEs in deforming domains is
still in its infancy, and our observations provide motivation for further
analysis of such patterning dynamics.

We note that to obtain our results, we set the diffusion length of the active
forms of GTPase to be 1/10 of those of the inactive form. This assured that
the length scale of the Rac-front would be of the right width relative to cell
diameter, leading to front-back cycles instead of waves. Decreasing the
diffusion length of the active form (or increasing cell size) generates
multiple interacting waves. The kymographs of the Rac behavior at the cell
edge in the 2D static domains were similar, but the wave patterning in the
cell interior changed. In the deforming cell, this leads to even more
irregular phenotypes than shown here, as multiple co-existing waves hit the
cell edge in different locations.

We found that cell deformation can reinforce patterns that were only weakly
stable on a static domain. For example, cells that elongate tend to become
more stably polarized; the Rac front tends to concentrate at one end of the
oval, rather than slowly spiraling around. There is then a positive feedback
between polarization and further polarization, favouring persistence. In a
different parameter regime, cell elongation can also evolve spirals into
front-back oscillations.

In correspondence with the two-compartment model, we find the persistent and
oscillatory phenotype, and cells with either uniform high Rac or uniform low
Rac. We find that by increasing the feedback between ECM and Rho, persistent
cells become oscillatory. Further increasing the feedback parameter increases
the frequency of the oscillations. By extending our model to static 2D
domains, we found a new phenotype within the oscillatory regime: spiraling of
Rac around the cell edge, that may turn around to spiral in the opposite
direction. By simulating cell edge deformations, we also find intermediate
phenotypes: cells that oscillate with irregular frequencies, and persistent
cells that occasionally attempt but fail to re-polarize, and remain
persistent. We also find cells that create and retract protrusions at random
new locations (the random phenotype). Such behaviours cannot be detected in
the 1D (or even static 2D) simulations. Finally, we observed that spontaneous
pattern formation in motile cells occurred in parameter regimes where static
simulation predicted uniform solutions.

We found that changes in biophysical properties such as cell stiffness can
affect the internal patterning dynamics. So, the same signaling network and
rates, but slightly different biophysical parameters result in very different
overall cell behaviour. This suggests that cell biophysics is an important
factor in cell migratory behaviour that should be more closely examined in the
context of intracellular signaling. See [34] for an example of this type.

We developed our Models II and III based on integrin bond dynamics modeled by
Novikova and Storm [56]. That is, adhesions are depicted as clusters of catch-
slip bonds that stabilize with intermediate force, and break beyond some force
threshold. Many other models for catch-bond dynamics exist, with slightly
different assumptions and greater complexity [52]. The advantage of the
version we used is in its simplicity and the fact that it was fit to
experimental data for force-lifetime curves of $\alpha 5-\beta 1$ integrin
bonds [56]. This type of integrin is expressed by many cells and binds to
fibronectin. Future data for force-lifetime curves of other cell types, other
ECMs, annd integrin types could be used to refine parameters and equations.
Other interesting extensions include modeling adhesions as mixtures of
different integrin bonds (e.g. combination of slip and catch-slip bonds [87]).
Such model extensions could help to further explain the response of
intracellular signaling to structural and spatial variations in the ECM and
integrin expression of cells [48, 88].

Real cells exhibit a wide range of dynamic behaviour. Examples include actin
waves, lamellipodial ruffles, random filopodial protrusions, circular cell
motility, oscillations in situ, and many more. The underlying mechanisms are
often unclear, and the distribution of signaling proteins in such dynamics is
only rarely characterized. Here, we have shown that various exotic dynamics
are inherent in even simple subsets of signaling networks. Such internal
dynamics can affect a cell’s ability to migrate as well as its invasive
potential.

Our results also contribute to further research on ECM-cell signaling. We have
expanded on the work of [42] with more highly resolved models, finding overall
agreement. We are able to actually simulate the class of random polarity and
motility that previous models did not resolve. Finally, the specific details
of the model could be altered to generalized to other sources of feedback or
to other hypotheses about the influence of ECM signaling.

Our models have several limitations that result from simplifications or
assumptions we made. First, we have adopted the hypothesis from [42] that the
Rac-Rho module accounts for major aspects of melanoma cell migratory phenotype
control. This patently ignores thousands of components that provide additional
feedback, input, or fine-tuning. We also directly linked Rac and Rho to cell
edge dynamics, ignoring the delays in recruitment and assembly of actin, as
well as the activation of myosin motors. We have also ignored the
inhomogeneities of the cell and its thickness in the third dimension. All
these are common simplifications.

Another simplification is our expression for the traction forces of the cell.
We employed a simple phenomenological relation between traction force and
levels of Rac/Rho in the cell. More detailed force models, as in [89, 36]
could be used to increase accuracy and relate other biophysical parameters
(stiffness, membrane tension) to the forces that are applied on the integrins.

Results in the paper provide examples of possible outcomes. The kymographs for
moving cells are for one realization only, and tend to vary somewhat from one
simulation to another. Both the initial conditions and also the stochastic
nature of the CPM affect the given outcome. As noted, there are regions of
parameter space in which multiple behaviours coexist, and while we display
examples of this type, we made no attempt to fully characterize these.

With the basic behaviour mapped out into regimes in parameter space, our next
investigation (ongoing) is to understand implications to real cells and
compare to new experiments. The full model can now be used to study the effect
of local variations in the ECM (or post density arrays) that lead to directed
cell motility. The ability of cells to migrate towards topographical cues
(topotaxis) as well as the effect of cell and substrate stiffness provide a
rich new set of phenomena to explore computationally. Cell migration also
affects the deposition and degradation of ECM [90], providing yet another set
of phenomena for investigation in the future.

##  10\. Acknowledgements

While at UBC, EGR was funded by a Natural Sciences and Engineering Research
Council of Canada (NSERC) Discovery Grant 41870 to LEK. We thank Andreas
Buttenschön for his FastVector library and general help with optimizing the
CPM code written by EGR. We are grateful to Lutz Brusch, Jörn Starruß and
Andreas Deutsch for the open software Morpheus, and to Bard Ermentrout for
XPPauto. We are grateful to Zachary Pellegrin for his development of a CPM
reaction-diffusion solver. We thank Andre Levchenko and JinSeok Park for early
discussions about this project. We thank Andreas Buttenschön for his
FastVector library and general help with building the CPM code.

##  11\. Appendix

###  11.1. Full Model Equations

For convenience, we gather full model equations in this section.

####  11.1.1. The well-mixed Rac-Rho system

(11.1a) $\displaystyle\frac{dR}{dt}$ $\displaystyle=$
$\displaystyle\frac{b_{R}}{1+\rho_{a}^{3}}(R_{T}-R)-\delta R,$ (11.1b)
$\displaystyle\frac{d\rho}{dt}$ $\displaystyle=$
$\displaystyle\frac{b_{\rho}}{1+R_{a}^{3}}(\rho_{T}-\rho)-\delta\rho.$

####  11.1.2. The well-mixed Rac-Rho-ECM Model I

(11.2a) $\displaystyle\frac{dR}{dt}$
$\displaystyle=\frac{b_{R}}{(1+\rho^{3})}(R_{T}-R)-\delta R,$ (11.2b)
$\displaystyle\frac{d\rho}{dt}$
$\displaystyle=\left(k_{E}+\gamma_{E}\frac{E^{m}}{E_{0}^{m}+E^{m}}\right)\left(\frac{1}{1+R^{3}}\right)(\rho_{T}-\rho)-\delta\rho,$
(11.2c) $\displaystyle\frac{dE}{dt}$
$\displaystyle=\epsilon\left[\left(k_{R}+\gamma_{R}\frac{R^{n}}{R_{0}^{n}+R^{n}}\right)-E\left(k_{\rho}+\gamma_{\rho}\frac{\rho^{n}}{\rho_{0}^{n}+\rho^{n}}\right)\right].$

####  11.1.3. Spatially distributed Rac-Rho-ECM Model I

In each of the spatial model variants, we can no longer eliminate the total
GTPase, since now conservation leads to an integral constraint,

$R_{T}=\int_{\Omega}(R+R_{I})d\Omega=\mbox{constant},$

and similarly for Rho (where $\Omega$ is the domain representing the cell).
Hence we keep the full PDEs for active and inactive Rac, and active and
inactive Rho.

(11.3a) $\displaystyle\frac{dR}{dt}$
$\displaystyle=\frac{b_{R}}{(1+\rho^{m})}R_{I}-\delta R+D\Delta
R,\quad\frac{dR_{I}}{dt}=-\frac{b_{R}}{(1+\rho^{m})}R_{I}+\delta R+D_{I}\Delta
R_{I},$ (11.3b) $\displaystyle\frac{d\rho}{dt}$
$\displaystyle=\frac{b_{\rho}(E)}{1+R^{m}}\rho_{I}-\delta\rho\,+D\Delta\rho,\quad\frac{d\rho_{I}}{dt}=-\frac{b_{\rho}(E)}{1+R^{m}}\rho_{I}+\delta\rho\,\,+D_{I}\Delta\rho_{I},$
(11.3c) $\displaystyle\frac{dE}{dt}$ $\displaystyle=\epsilon(a(R)-d(\rho)E),$
where the functions appearing above are (11.3d) $\displaystyle b_{\rho}(E)$
$\displaystyle=k_{E}+\gamma_{E}\frac{E^{m}}{(E_{0}^{m}+E^{m})},$ (11.3e)
$\displaystyle a(R)$
$\displaystyle=K+\gamma_{R}\frac{R^{n}}{(R_{0}^{n}+R^{n})},$ (11.3f)
$\displaystyle d(\rho)$
$\displaystyle=k_{P}+\gamma_{\rho}\frac{\rho^{n}}{(\rho_{0}^{n}+\rho^{n})}.$

We assume that $D_{R}=D_{\rho}=D$, since GTPases have similar sizes. The
active form diffuses much slower than the inactive form, so $D<<D_{I}$.

####  11.1.4. Spatially distributed Models II and III

We use the same Rac -Rho equations as in (11.3)a,b, d, with a new PDE for the
ECM variable:

(11.4) $\frac{dE}{dt}=\epsilon\left(K(E_{t}-E)-d(F,E)E\right).$

where the force on an integrin bond is

(11.5)
$F(\rho,R)=\mbox{max}\left(0,\frac{\beta_{\rho}\rho}{(1+\rho)}-\frac{\beta_{R}R}{(1+R)}\right).$

For the slip-bond (Model II) we use

(11.6) $d(F,E)=k_{0}\exp\left(\frac{F}{p(E+E_{s})}\right).$

whereas for the slip-catch bond (Model III) we use

(11.7)
$d(F,E)=k_{0}\exp\left(\frac{F}{p(E+E_{s})}\right)+k_{0c}\exp\left(-\frac{F}{p(E+E_{s})}\right).$

Altogether, the model structure for the active forms in Models II and III is

(11.8a) $\displaystyle\frac{dR}{dt}$
$\displaystyle=\frac{b_{R}}{(1+\rho^{m})}R_{I}-\delta R+D\Delta R,$ (11.8b)
$\displaystyle\frac{d\rho}{dt}$
$\displaystyle=\left(k_{E}+\gamma_{E}\frac{E^{m}}{(E_{0}^{m}+E^{m})}\right)\cdot\left(\frac{1}{1+R^{m}}\right)\rho_{I}-\delta\rho+D\Delta\rho,$
(11.8c) $\displaystyle\frac{dE}{dt}$
$\displaystyle=\epsilon\left(K(E_{t}-E)-d(F(R,\rho),E)E\right),$

together with two additional PDEs for the inactive GTPases, $R_{I},\rho_{I}$,
as before in (11.3).

###  11.2. Parameter values

We employed nondimensionalized parameters (see Table 3). We can define a
spatial and temporal scale, by comparing with melanoma cells in [3, 42]. These
cells have a diameter of roughly 30 $\mu$m [42] and our CPM cells have a
diameter of 45 pixels. So 1 CPM pixel can be assigned a width of 2/3 $\mu$m.
Given that one CPM pixel is 0.075 in non-dimensional units, our spatial scale
is then given by $L_{0}\approx 9\mu$m. We estimate the speed of persistent
melanoma cells in [3] as 0.1 $\mu$m/min, similar to melanoma migration speeds
reported elsewhere, of 5-40 $\mu$m/hr [91, 92]. The simulated persistent cell
in Fig. 12A has a speed of roughly 0.02 $\mu$m/MCS. During one MCS, we perform
1 timestep of the reaction-diffusion equations. So, one MCS is associated with
$t_{0}$=12 s. These temporal and spatial scale result in a diffusion rate of
$\frac{L_{0}^{2}}{t_{0}}\approx 6.5\mu$m2/s, for the inactive form of the
GTPases. This gives realistic values for our parameters, and they can be
adapted to link with different types of cells and conditions.

###  11.3. Kymographs in 2D

In both static and deforming 2D, we produce kymographs by tracking the Rac
level on the cell edge. We identified the cell edge using the edge-detection
function bwboundaries in Matlab. To align the Rac front to the left of the
kymograph, we first identify the location of highest Rac relative to the cell
centre (See methods in [42]) and rotate the image by that angle.

###  11.4. Implementation details, Morpheus

####  11.4.1. Morpheus software and files

Results | Figure | Morpheus xml file
---|---|---
Model I in 1D | Fig. 6 | Model1RacRhoECMPDEsIn1D
Model II in 1D | Fig. 7 | Model2SlipBondRacRhoECMPDEsIn1D
Model III in 1D | Fig. 8 | Model3CatchSlipBondRacRhoECMPDEsIn1D
One cell, WM Model I | Fig. 5 | OneLamelShape
2Lam Model I | Fig. 16 | Melanoma2LamelCPMRac-Rho-ECMWellMixed
Table 2. A list of the Morpheus files used to prepare a number of figures. WM=
well-mixed, 2Lam = two -lamellipod.

Morpheus is an open-source software package for multiscale simulations of
single and interacting cells [59]. Morpheus solves differential equations, and
simulates the cellular-Potts model to depict cells. Morpheus produces small
.xml files that allow users to exactly reproduce simulation results.

####  11.4.2. Cell shapes with well-mixed internal Rac-Rho-ECM

For the cell shape simulation with well-mixed internal variables, we linked
the internal chemistry to the target area $A_{0}$ in the Hamiltonian,

(11.9) $H=\lambda_{A}(A-A_{0})^{2}+\dots$

One cell compartment: For the single-cell compartment shown in Figure 5, we
used the well-mixed Rac-Rho-ECM Model I and assumed that

$A_{0}=400+1000\cdot R-100\cdot\rho.$

Two cell compartments: We also ran Morpheus simulations in which a single cell
was depicted by by a pair of geometric objects as in [93]. This is directly
motivated by [43] where two “lamellipods” compete for “frontnes”. As example
of this type is shown in Figure 16. The levels of active Rac and Rho,
$R_{1},R_{2},\rho_{1},\rho_{2}$ in each lamellipod, and the ECM contact area
of each, $E_{1},E_{2}$ are given by the system (11.2) but the compartments are
coupled by assuming that,

(11.10) $R_{T}=R_{I}+R_{1}+R_{2},\quad\rho_{T}=\rho_{I}+\rho_{1}+\rho_{2}.$

Since inactive GTPases $R_{I},\rho_{I}$ diffuse rapidly, it is reasonable to
assume that their levels are the same in both compartments. Hence, we can
eliminate $R_{I},\rho_{I}$ from the well-mixed Model I ODE system. We assumed
that the target area of each compartment was given by

$A_{0}=400+1500R-300\rho.$

The adhesive energy between the two compartments is kept purposely low to
ensure that they adhere to form a single “cell”. As shown in Figure 16, the
two lamellipods oscillate in antiphase manner, with high Rac in one
corresponding to low Rac in the other.

Figure 16. The two-lamellipod model. Top left: a schematic diagram indicating
how the various interactions are modeled based on Eqs. (11.2) of Model I. Top
right: A plot of the variables over time in each of the two lamellipods,
created with Morpheus file Melanoma2LamelCPMRac-Rho-ECMWellMixed.xml. The red
arrow spans the time-frame displayed across the bottom row. Bottom: Morpheus
simulation showing the “shape” of a cell consisting of the two compartments.
The color scheme gives the Rac level in each compartment. The areas of the
“front” and “back” lamellipods evolve as the levels of GTPases oscillates.
First Rac is high on the left, and then it shifts to being high on the right.
Produced with Morpheus file Melanoma2LamelCPMRac-Rho-ECMWellMixed.xml

###  11.5. Implementation details, custom-built CPM simulations

The CPM is an energy-based cell shape computation, where the energy is
represented by a Hamiltonian [94, 67]. In a single-cell model, the Hamiltonian
has typical terms that include the following:

(11.11) $H=\lambda_{A}(A-A_{0})^{2}+JP+\lambda_{P}(P-P_{0})^{2}.$

These terms include energetic penalties for cell area that deviates from a
“rest area” $A_{0}$, for perimeter deviation from the value $P_{0}$. There is
also an energy associated with cell-medium contact, depicted by the term $JP$.

A cell is represented by some set of contiguous pixels (in 2D for our paper).
During each Monte-Carlo step (MCS), the cell attempts a number of
protrusions/retractions at pixels along its boundary, and accepts all attempts
that decrease the Hamiltonian. Otherwise, the probability of accepting a
change that increases $H$ is given by a Boltzman probability distribution that
depends on the ratio of energy gain to some energy fluctuation (“thermal
energy” term). More information about the CPM formalism and a comparison with
other methods of simulating cells is provided in [67, 89].

####  11.5.1. Imposed cell elongation

We added cell elongation to avoid circular motion and to keep cell trajectory
straight with a persistent front. (See [33] for circular motion in a simple
GTPase cell-shape model.) To impose cell elongation we add an additional term
to the Hamiltonian

(11.12) $H_{e}=\lambda_{e}\cdot(100\cdot\xi-\xi_{0})^{2},$

where $\xi$ is the eccentricity of the cell. To approximate the eccentricity
[95], we used the inertia tensors $I(S)$ of the cells, defined as

(11.13) $I=\left(\begin{matrix}\sum(y-c_{y})^{2}&-\sum(x-c_{x})(y-c_{y})\\\
-\sum(x-c_{x})(y-c_{y})&\sum(x-c_{x})^{2}\end{matrix}\right),$

where the sum goes over all the coordinates within the cell, and
$(c_{x},c_{y})$ is the center of mass of the cell. The two eigenvalues $a,b$,
with $a>b$, of this matrix approximate the relative sizes of the long and
short axes of the cell. The eccentricity of a cell is defined as
$\xi=\sqrt{1-\left({\frac{b}{a}}\right)^{2}}$. A circular cell has $\xi\approx
0$, while an elongated cell has $\xi\approx 1$. This elongation constraint
deviates from the typical cell length constraint (see e.g. [95]), which only
affects the long and not the short axis. Using the old constraint, a cell
tended to round up. Our new elongation constraint prevented this and assured
that the cell long and short axes were significantly different.

####  11.5.2. Forces due to Rac-Rho activities

Technically, we implement the Rac-Rho dependent edge forces through the CPM
Hamiltonian, as done, e.g., in [26, 89]. That is, we modify the local change
in the Hamiltonian (equivalent to a local force) as follows

(11.14) $\displaystyle\Delta
H=\beta_{G}\cdot\begin{cases}\beta_{R}\frac{R_{a}}{1+R_{a}}-\beta_{P}\frac{\rho_{a}}{1+\rho_{a}},&\text{if
retraction},\\\
-\beta_{R}\frac{R_{a}}{1+R_{a}}+\beta_{P}\frac{\rho_{a}}{1+\rho_{a}},&\text{if
extension}.\end{cases}$

Here, $\beta_{G}$ sets the relative change in energy ($\approx$ force)
contribution of the GTPases compared to terms in the general Hamiltonian.

####  11.5.3. Forces needed to overcome adhesion bonds:

To make the cells less motile when their ECM attachment is higher, we
implement an additional ECM-dependent yield energy, $H_{0}$, (a kind of
dissipative energy) that the cell needs to overcome to make a movement. To do
so, we used the probability function

(11.15) $P(\Delta H)=\begin{cases}1&\textrm{if}\;\Delta H+H_{0}<0\\\
e^{-(\Delta H+H_{0})/T}&\textrm{if}\;\Delta H+H_{0}\geq 0,\end{cases}$

where $H_{0}$ depends on the cell-ECM adhesivity as follows:

(11.16) $H_{0}=\begin{cases}0&\text{if extension}\\\ \lambda_{N}E&\text{if
retraction}.\end{cases}$

####  11.5.4. Initial conditions

For the 2D simulations, we use a circular domain of diameter $L$=3. We imposed
uniform noise in active Rac throughout the cell with a maximal level of 0.8.

In 1D and static 2D simulations, we used grid size $dx=0.1$ to solve the PDEs,
but for the deforming cells, we increased the resolution by setting
$dx=0.075$.

We initiate the simulation with a cell with uniform adhesion:
$N(\vec{x})=N_{u}$, where $N_{u}$ is some small value, chosen here as 50.
Whenever the cell protrudes, we assume it binds to the ECM and set
$N(\vec{x})=N_{0}=50$. After a cell retracts, the ECM attachment is set to
zero.

####  11.5.5. Redistribution of fields after cell moves

After a protrusion, the newly created cell site obtains the GTPase levels of
the pixel that was “copied” to the target site. After a retraction, the GTPase
level formerly at a cell site is set to zero. After each MCS, we normalize the
GTPase levels (over the whole cell) to ensure conservation of the total GTPase
of each type.

parameter | description | value
---|---|---
space/time |  |
$dx$ | lattice site (pixel) width | 0.075
$dt$ | PDE integration time-step | 0.00125
$n_{i}$ | number of PDE integration time steps | 800
$N$ | number of Monte Carlo Steps in CPM | 1000
CPM |  |
$A_{0}$ | target area | 1600
$\lambda_{A}$ | strength of area constraint | 0.5
$P_{0}$ | target perimeter | 200
$\lambda_{P}$ | strength of perimeter constraint | 12
$e_{0}$ | target eccentricity | 90
$\lambda_{e}$ | strength of eccentricity constraint | 20
$J$ | adhesive energy | 40000
$n_{r}$ | neighbourhood radius for $J\cdot P$ | 3
$\xi(r)$ | scaling factor for $J\cdot P$ | 18
$\beta_{G}$ | GTPase force strength | 0.5
$\lambda_{N}$ | adhesion strength | 0.3
$T$ | cellular temperature | 200
Rac/Rho |  |
$D_{I}$ | diffusion of inactive Rac and Rho | 1
$D$ | diffusion of active Rac and Rho | 0.1
$\delta$ | Rac, Rho inactivation rate | 1
$b_{R}$ | Rac activation rate | 2
ECM-Rho coupling |  |
$k_{E}$ | basal Rho activation rate | 2
$\beta_{R}$ | protrusive force of Rac | 1200 or 1000
$\beta_{P}$ | contractile force of Rho | 1600
$E_{0}$ | half-max adhesion for Rho activation | 300
$\gamma_{E}$ | adhesion dependent Rho activation rate | 4
Adhesion model |  |
$r$ | integrin binding rate | 10
$E_{T}$ | maximum adhesion size | 1000
$k_{0}$ | integrin unbinding rate due to slip | 5 or 0.0004
$k_{0c}$ | integrin unbinding rate due to catch | 0 or 55
$\epsilon$ | time-scale | 0.001
$p_{s}$ | scaling factor force | 0.35 or 0.08
$E_{s}$ | small adhesion for scaling | 100
$E_{n}$ | nascent adhesion at new protrusion | 50
Table 3. Dimensionless parameters for the CPM simulations. The value options
indicated by ”value1 or value2” signify the different values chosen for Model
II (slip) or Model III (catch-slip).

## References

  * [1] Houk AR, Jilkine A, Mejean CO, Boltyanskiy R, Dufresne ER, Angenent SB, et al.  Membrane tension maintains cell polarity by confining signals to the leading edge during neutrophil migration.  Cell. 2012;148(1-2):175–188.
  * [2] Charras G, Sahai E.  Physical influences of the extracellular environment on cell migration.  Nature Reviews Molecular Cell Biology. 2014;15(12):813–824.
  * [3] Park J, Kim DH, Kim HN, Wang CJ, Kwak MK, Hur E, et al.  Directed migration of cancer cells guided by the graded texture of the underlying matrix.  Nature materials. 2016;.
  * [4] Ridley AJ, Schwartz MA, Burridge K, Firtel RA, Ginsberg MH, Borisy G, et al.  Cell migration: integrating signals from front to back.  Science. 2003;302(5651):1704–1709.
  * [5] Ridley AJ.  Rho family proteins: coordinating cell responses.  Trends in cell biology. 2001;11(12):471–477.
  * [6] Etienne-Manneville S, Hall A.  Rho GTPases in cell biology.  Nature. 2002;420(6916):629.
  * [7] Burridge K, Wennerberg K.  Rho and Rac take center stage.  Cell. 2004;116(2):167–179.
  * [8] Etienne-Manneville S.  Cdc42-the centre of polarity.  Journal of cell science. 2004;117(8):1291–1300.
  * [9] Bourne HR, Sanders DA, McCormick F.  The GTPase superfamily: a conserved switch for diverse cell functions.  Nature. 1990;348(6297):125–132.
  * [10] Burridge K, Doughman R.  Front and back by Rho and Rac.  Nature cell biology. 2006;8(8):781.
  * [11] Maekawa M, Ishizaki T, Boku S, Watanabe N, Fujita A, Iwamatsu A, et al.  Signaling from Rho to the actin cytoskeleton through protein kinases ROCK and LIM-kinase.  Science. 1999;285(5429):895–898.
  * [12] Machacek M, Hodgson L, Welch C, Elliott H, Pertz O, Nalbant P, et al.  Coordination of Rho GTPase activities during cell protrusion.  Nature. 2009;461(7260):99–103.
  * [13] Guilluy C, Garcia-Mata R, Burridge K.  Rho protein crosstalk: another social network?  Trends in cell biology. 2011;21(12):718–726.
  * [14] Parri M, Chiarugi P.  Rac and Rho GTPases in cancer cell motility control.  Cell Commun Signal. 2010;8(23):10–1186.
  * [15] Sailem H, Bousgouni V, Cooper S, Bakal C.  Cross-talk between Rho and Rac GTPases drives deterministic exploration of cellular shape space and morphological heterogeneity.  Open biology. 2014;4(1):130132.
  * [16] Pertz O, Hodgson L, Klemke RL, Hahn KM.  Spatiotemporal dynamics of RhoA activity in migrating cells.  Nature. 2006;440(7087):1069–1072.
  * [17] Bolado-Carrancio A, Rukhlenko OS, Nikonova E, Tsyganov MA, Wheeler A, Garcia-Munoz A, et al.  Periodic propagating waves coordinate RhoGTPase network dynamics at the leading and trailing edges during cell migration. 2020 jul;9:e58165.  Available from: https://doi.org/10.7554/eLife.58165.
  * [18] Tomar A, Schlaepfer DD.  Focal adhesion kinase: switching between GAPs and GEFs in the regulation of cell motility.  Current opinion in cell biology. 2009;21(5):676–683.
  * [19] Sheetz MP, Felsenfeld DP, Galbraith CG.  Cell migration: regulation of force on extracellular-matrix-integrin complexes.  Trends in cell biology. 1998;8(2):51–54.
  * [20] Lessey EC, Guilluy C, Burridge K.  From mechanical force to RhoA activation.  Biochemistry. 2012;51(38):7420–7432.
  * [21] Poh YC, Na S, Chowdhury F, Ouyang M, Wang Y, Wang N.  Rapid activation of Rac GTPase in living cells by force is independent of Src.  PloS one. 2009;4(11):e7886.
  * [22] Lee G, Han SB, Kim DH.  Cell-ECM contact-guided intracellular polarization is mediated via Lamin A/C dependent nucleus-cytoskeletal connection.  Biomaterials. 2020;p. 120548.
  * [23] Plotnikov SV, Pasapera AM, Sabass B, Waterman CM.  Force fluctuations within focal adhesions mediate ECM-rigidity sensing to guide directed cell migration.  Cell. 2012;151(7):1513–1527.
  * [24] Aoki K, Nakamura T, Fujikawa K, Matsuda M.  Local Phosphatidylinositol 3,4,5-Trisphosphate Accumulation Recruits Vav2 and Vav3 to Activate Rac1/Cdc42 and Initiate Neurite Outgrowth in Nerve Growth Factor-stimulated PC12 Cells.  Mol Biol Cell. 2005;16(5):2207–2217.
  * [25] Xia N, Thodeti CK, Hunt TP, Xu Q, Ho M, Whitesides GM, et al.  Directional control of cell motility through focal adhesion positioning and spatial control of Rac activation.  The FASEB journal. 2008;22(6):1649–1659.
  * [26] Marée AFM, Jilkine A, Dawes A, Grieneisen VA, Edelstein-Keshet L.  Polarization and Movement of Keratocytes: A Multiscale Modelling Approach.  Bull Math Biol. 2006;68:1169–1211.
  * [27] Jilkine A, Edelstein-Keshet L.  A Comparison of Mathematical Models for Polarization of Single Eukaryotic Cells in Response to Guided Cues.  PLoS Comput Biol. 2011;7(4):e1001121.
  * [28] Holmes WR, Lin B, Levchenko A, Edelstein-Keshet L.  Modelling Cell Polarization Driven by Synthetic Spatially Graded Rac Activation.  PLoS Comput Biol. 2012 06;8(6):e1002366.
  * [29] Kopfer KH, Jäger W, Matthäus F.  A mechanochemical model for rho GTPase mediated cell polarization.  Journal of Theoretical Biology. 2020;504:110386.
  * [30] Wang W, Tao K, Wang J, Yang G, Ouyang Q, Wang Y, et al.  Exploring the inhibitory effect of membrane tension on cell polarization.  PLoS computational biology. 2017;13(1):e1005354.
  * [31] Tao K, Wang J, Kuang X, Wang W, Liu F, Zhang L.  Tuning cell motility via cell tension with a mechanochemical cell migration model.  Biophysical Journal. 2020;.
  * [32] Cao Y, Ghabache E, Rappel WJ.  Plasticity of cell migration resulting from mechanochemical coupling.  Elife. 2019;8:e48478.
  * [33] Camley BA, Zhao Y, Li B, Levine H, Rappel WJ.  Crawling and turning in a minimal reaction-diffusion cell motility model: Coupling cell shape and biochemistry.  Physical Review E. 2017;95(1):012401.
  * [34] Zmurchok C, Collette J, Rajagopal V, Holmes WR.  Membrane Tension Can Enhance Adaptation to Maintain Polarity of Migrating Cells.  Biophysical Journal. 2020 Oct;119(8):1617–1629.  Available from: https://doi.org/10.1016/j.bpj.2020.08.035.
  * [35] Copos C, Mogilner A.  A hybrid stochastic–deterministic mechanochemical model of cell polarization.  Molecular biology of the cell. 2020;31(15):1637–1649.
  * [36] Rens EG, Merks RM.  Cell shape and durotaxis explained from cell-extracellular matrix forces and focal adhesion dynamics.  Iscience. 2020;23(9):101488.
  * [37] Deshpande VS, Mrksich M, McMeeking RM, Evans AG.  A Bio-mechanical Model For Coupling Cell Contractility With Focal Adhesion Formation.  Journal of the Mechanics and Physics of Solids. 2008;56(4):1484 – 1510\.  Available from: http://www.sciencedirect.com/science/article/pii/S002250960700172X.
  * [38] Ronan W, Deshpande VS, Mcmeeking RM, Mcgarry JP.  Cellular Contractility And Substrate Elasticity: A Numerical Investigation Of The Actin Cytoskeleton And Cell Adhesion.  Biomechanics and Modeling in Mechanobiology. 2014;13(2):417–435.
  * [39] Vernerey FJ, Farsad M.  A Mathematical Model Of The Coupled Mechanisms Of Cell Adhesion , Contraction And Spreading.  Journal of Mathematical Biology. 2014;68(4):989–1022.
  * [40] Bangasser BL, Shamsan GA, Chan CE, Opoku KN, Tüzel E, Schlichtmann BW, et al.  Shifting the optimal stiffness for cell migration.  Nature communications. 2017;8(1):1–10.
  * [41] Hastings JF, Skhinas JN, Fey D, Croucher DR, Cox TR.  The extracellular matrix as a key regulator of intracellular signalling networks.  British journal of pharmacology. 2019;176(1):82–92.
  * [42] Park J, Holmes WR, Lee SH, Kim HN, Kim DH, Kwak MK, et al.  Mechanochemical feedback underlies coexistence of qualitatively distinct cell polarity patterns within diverse cell populations.  Proceedings of the National Academy of Sciences. 2017;114(28):E5750–E5759.
  * [43] Holmes WR, Park J, Levchenko A, Edelstein-Keshet L.  A mathematical model coupling polarity signaling to cell adhesion explains diverse cell migration patterns.  PLoS computational biology. 2017;13(5):e1005524.
  * [44] Larsen M, Artym VV, Green JA, Yamada KM.  The matrix reorganized: extracellular matrix remodeling and integrin signaling.  Current opinion in cell biology. 2006;18(5):463–471.
  * [45] Costa P, Scales TM, Ivaska J, Parsons M.  Integrin-specific control of focal adhesion kinase and RhoA regulates membrane protrusion and invasion.  PLoS One. 2013;8(9):e74659.
  * [46] Huveneers S, Danen EH.  Adhesion signaling–crosstalk between integrins, Src and Rho.  Journal of cell science. 2009;122(8):1059–1069.
  * [47] Lawson CD, Burridge K.  The on-off relationship of Rho and Rac during integrin-mediated adhesion and cell migration.  Small GTPases. 2014;5(1):e27958.
  * [48] Danen EH, Sonneveld P, Brakebusch C, Fässler R, Sonnenberg A.  The fibronectin-binding integrins $\alpha$5$\beta$1 and $\alpha$v$\beta$3 differentially modulate RhoA–GTP loading, organization of cell matrix adhesions, and fibronectin fibrillogenesis.  The Journal of cell biology. 2002;159(6):1071–1086.
  * [49] Park JH, Ryu JM, Han HJ.  Involvement of caveolin-1 in fibronectin-induced mouse embryonic stem cell proliferation: Role of FAK, RhoA, PI3K/Akt, and ERK 1/2 pathways.  Journal of cellular physiology. 2011;226(1):267–275.
  * [50] Arthur WT, Petch LA, Burridge K.  Integrin engagement suppresses RhoA activity via a c-Src-dependent mechanism.  Current Biology. 2000;10(12):719–722.
  * [51] Kong F, García AJ, Mould AP, Humphries MJ, Zhu C.  Demonstration of catch bonds between an integrin and its ligand.  Journal of Cell Biology. 2009;185(7):1275–1284.
  * [52] MacKay L, Khadra A.  The bioenergetics of integrin-based adhesion, from single molecule dynamics to stability of macromolecular complexes.  Computational and structural biotechnology journal. 2020;18:393–416.
  * [53] Nicolas A, Geiger B, Safran SA.  Cell Mechanosensitivity Controls The Anisotropy Of Focal Adhesions.  Proceedings of the National academy of Sciences of the United States of America. 2005;101(34):12520–12525.
  * [54] Besser A, Safran SA.  Force-induced Adsorption And Anisotropic Growth Of Focal Adhesions.  Biophysical Journal. 2006;90(10):3469–3484.
  * [55] MacKay L, Khadra A.  Dynamics of mechanosensitive nascent adhesion formation.  Biophysical journal. 2019;117(6):1057–1073.
  * [56] Novikova EA, Storm C.  Contractile fibers and catch-bond clusters: a biological force sensor?  Biophysical journal. 2013;105(6):1336–1345.
  * [57] Holmes WR, Mata MA, Edelstein-Keshet L.  Local Perturbation Analysis: A Computational Tool for Biophysical Reaction-Diffusion Models.  Biophysical journal. 2015;108(2):230–236.
  * [58] Holmes WR, Edelstein-Keshet L.  Analysis of a minimal Rho-GTPase circuit regulating cell shape.  Physical Biology. 2016;13:046001.
  * [59] Starruß J, de Back W, Brusch L, Deutsch A.  Morpheus: a user-friendly modeling environment for multiscale and multicellular systems biology.  Bioinformatics. 2014;30(9):1331–1332.
  * [60] Lin B, Holmes WR, Wang J, Ueno T, Harwell A, Edelstein-Keshet L, et al.  Synthetic spatially graded Rac activation drives directed cell polarization and locomotion.  PNAS. 2012;109(52):E3668–E3677.
  * [61] Sharma P, Ng C, Jana A, Padhi A, Szymanski P, Lee JS, et al.  Aligned fibers direct collective cell migration to engineer closing and nonclosing wound gaps.  Molecular biology of the cell. 2017;28(19):2579–2588.
  * [62] Padhi A, Alabduljabbar DS, Kapania RK, Nain AS.  Contractile Cytoplasmic Fragments Exhibit ROCK-dependent Oscillatory Motion.  American Society for Cell Biology. 2020;.
  * [63] Padhi A, Danielsson BE, Alabduljabbar DS, Wang J, Conway DE, Kapania RK, et al.  Cell Fragment Formation, Migration, and Force Exertion on Extracellular Mimicking Fiber Nanonets.  Advanced Biology. 2021;2000592.
  * [64] Cickovski TM, Huang C, Chaturvedi R, Glimm T, Hentschel HGE, Alber MS, et al.  A framework for three-dimensional simulation of morphogenesis.  IEEE/ACM Transactions on Computational Biology and Bioinformatics. 2005;2(4):273–288.
  * [65] Zhan H, Bhattacharya S, Cai H, Iglesias PA, Huang CH, Devreotes PN.  An excitable Ras/PI3K/ERK signaling network controls migration and oncogenic transformation in epithelial cells.  Developmental Cell. 2020;54(5):608–623.
  * [66] Liu Y, Rens EG, Edelstein-Keshet L.  Spots, strips, and spiral waves in models for static and motile cells.  arXiv preprint arXiv:190910504. 2019;.
  * [67] Magno R, Grieneisen VA, Marée AF.  The biophysical nature of cells: potential cell behaviours revealed by analytical and computational studies of cell surface mechanics.  BMC biophysics. 2015;8(1):1–37.
  * [68] Holmes W, Carlsson A, Edelstein-Keshet L.  Regimes of wave type patterning driven by refractory actin feedback: transition from static polarization to dynamic wave behaviour.  Physical biology. 2012;9(4):046005.
  * [69] Buttenschön A, Liu Y, Edelstein-Keshet L.  Cell size, mechanical tension, and GTPase signaling in the single cell.  Bulletin of mathematical biology. 2020;82(2):1–33.
  * [70] Zmurchok C, Bhaskar D, Edelstein-Keshet L.  Coupling mechanical tension and GTPase signaling to generate cell and tissue dynamics.  Physical biology. 2018;15(4):046004.
  * [71] Zmurchok C, Holmes WR.  Simple Rho GTPase Dynamics Generate a Complex Regulatory Landscape Associated with Cell Shape.  Biophysical Journal. 2020;.
  * [72] Winfree AT.  Spiral waves of chemical activity.  Science. 1972;175(4022):634–636.
  * [73] Tyson JJ, Murray J.  Cyclic AMP waves during aggregation of Dictyostelium amoebae.  Development. 1989;106(3):421–426.
  * [74] Winfree AT.  Varieties of spiral wave behavior: An experimentalist’s approach to the theory of excitable media.  Chaos: An Interdisciplinary Journal of Nonlinear Science. 1991;1(3):303–334.
  * [75] Cross MC, Hohenberg PC.  Pattern formation outside of equilibrium.  Reviews of modern physics. 1993;65(3):851.
  * [76] Keener JP.  A geometrical theory for spiral waves in excitable media.  SIAM Journal on Applied Mathematics. 1986;46(6):1039–1056.
  * [77] Rinzel J, Terman D.  Propagation phenomena in a bistable reaction-diffusion system.  SIAM Journal on Applied Mathematics. 1982;42(5):1111–1137.
  * [78] Mori Y, Jilkine A, Edelstein-Keshet L.  Wave-pinning and cell polarity from a bistable reaction-diffusion system.  Biophys J. 2008;94(9):3684–3697.
  * [79] Vicker MG.  Reaction–diffusion waves of actin filament polymerization/depolymerization in Dictyostelium pseudopodium extension and cell locomotion.  Biophysical chemistry. 2000;84(2):87–98.
  * [80] Whitelam S, Bretschneider T, Burroughs NJ.  Transformation from spots to waves in a model of actin pattern formation.  Physical review letters. 2009;102(19):198103.
  * [81] Bretschneider T, Anderson K, Ecke M, Müller-Taubenberger A, Schroth-Diez B, Ishikawa-Ankerhold HC, et al.  The three-dimensional dynamics of actin waves, a model of cytoskeletal self-organization.  Biophysical journal. 2009;96(7):2888–2900.
  * [82] Bernitt E, Döbereiner HG, Gov NS, Yochelis A.  Fronts and waves of actin polymerization in a bistability-based mechanism of circular dorsal ruffles.  Nature communications. 2017;8(1):1–8.
  * [83] Allard J, Mogilner A.  Traveling waves in actin dynamics and cell motility.  Current opinion in cell biology. 2013;25(1):107–115.
  * [84] Liu Y, Rens EG, Edelstein-Keshet L.  Spots, strips, and spiral waves in models for static and motile cells.  arXiv preprint arXiv:190910504. 2019;.
  * [85] Eroumé K, Vasilevich A, Vermeulen S, de Boer J, Carlier A.  On the influence of cell shape on dynamic reaction-diffusion polarization patterns.  Plos one. 2021;16(3):e0248293.
  * [86] Marée AFM, Grieneisen VA, Edelstein-Keshet L.  How Cells Integrate Complex Stimuli: The Effect of Feedback from Phosphoinositides and Cell Shape on Cell Polarization and Motility.  PLoS Comput Biol. 2012 03;8:e1002402.
  * [87] Novikova EA, Storm C.  Evolving roles and dynamics for catch and slip bonds during adhesion cluster maturation.  Physical Review E. 2021;103(3):032402.
  * [88] Danen EH, van Rheenen J, Franken W, Huveneers S, Sonneveld P, Jalink K, et al.  Integrins control motile strategy through a Rho–cofilin pathway.  The Journal of cell biology. 2005;169(3):515–526.
  * [89] Rens EG, Edelstein-Keshet L.  From energy to cellular forces in the Cellular Potts Model: An algorithmic approach.  PLoS Computational Biology. 2019;15(12).
  * [90] Park D, Wershof E, Boeing S, Labernadie A, Jenkins RP, George S, et al.  Extracellular matrix anisotropy is determined by TFAP2C-dependent regulation of cell collisions.  Nature materials. 2020;19(2):227–238.
  * [91] Li R, Hebert JD, Lee TA, Xing H, Boussommier-Calleja A, Hynes RO, et al.  Macrophage-secreted TNF$\alpha$ and TGF$\beta$1 influence migration speed and persistence of cancer cells in 3D tissue culture via independent pathways.  Cancer research. 2017;77(2):279–290.
  * [92] Howard JD, Moriarty WF, Park J, Riedy K, Panova IP, Chung CH, et al.  Notch signaling mediates melanoma–endothelial cell communication and melanoma cell migration.  Pigment cell & melanoma research. 2013;26(5):697–707.
  * [93] Schnyder SK, Molina JJ, Tanaka Y, Yamamoto R.  Collective motion of cells crawling on a substrate: roles of cell shape and contact inhibition.  Scientific reports. 2017;7(1):1–14.
  * [94] Graner F, Glazier JA.  Simulation of biological cell sorting using a two-dimensional extended Potts model.  Physical review letters. 1992;69(13):2013.
  * [95] Zajac M, Jones GL, Glazier JA.  Simulating convergent extension by way of anisotropic differential adhesion.  Journal of Theoretical Biology. 2003;222(2):247–259.