Cellular Tango: how extracellular matrix adhesion choreographs Rac-Rho signaling and cell movement

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–47]. Experimental evidence suggests that ECM signaling strongly enhances the activation of RhoA [48, 49], although it likely also has positive effects on Rac activity [45, 50].

Integrins cluster into focal adhesions (FAs), 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–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 and 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.

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, more biophysical detail, and enhanced computational tools.

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

$$\begin{align}& \frac{\partial R}{\partial t}={A}_{R}(\rho ){R}_{I}-{I}_{R}R+D{\triangle}R,\\ & \frac{\partial {R}_{I}}{\partial t}=-{A}_{R}(\rho ){R}_{I}+{I}_{R}R+{D}_{I}{\triangle}{R}_{I},\end{align} \tag{ 2.2a }$$

$$\begin{align}\frac{\partial \rho }{\partial t}& ={A}_{\rho }(R,E){\rho }_{I}-{I}_{\rho }\rho +D{\triangle}\rho ,\\ \frac{\partial {\rho }_{I}}{\partial t}& =-{A}_{\rho }(R,E){\rho }_{I}+{I}_{\rho }\rho +{D}_{I}{\triangle}{\rho }_{I},\end{align} \tag{ 2.2b }$$

where

$$\begin{equation}{A}_{R}(\rho )=\frac{{b}_{R}}{1+{\rho }^{n}},\quad {A}_{\rho }(R,E)=\frac{{b}_{\rho }(E)}{1+{R}^{n}},\quad {I}_{R}={I}_{\rho }=\delta .\end{equation} \tag{ 2.2c }$$

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

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

$$\begin{equation}\frac{\partial E}{\partial t}={\epsilon}(a(R,\rho ,E)-d(R,\rho ,E)E),\end{equation} \tag{ 2.3a }$$

$$\begin{equation}\text{and}\;\quad {b}_{\rho }(E)={k}_{E}+{\gamma }_{E}\frac{{E}^{n}}{({E}_{0}^{n}+{E}^{n})}.\end{equation} \tag{ 2.3b }$$

The Rho activation rate b_ρ (E) is affected by ECM signaling [48, 49]. In (2.3b), k_E , is a basal Rho activation rate, and γ_E is a parameter governing the strength of feedback from ECM adhesion to Rho activation. Figure 1(a) describes the above Rac-Rho-ECM feedback.

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

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. We keep these assumptions, but simplify to an elementary equation for E that has a rate of increase dependent on R and a decay rate dependent on ρ

$$\begin{equation*}\frac{\partial E}{\partial t}={\epsilon}(a(R)-d(\rho )E).\end{equation*}$$

(We have dropped the two-lamellipod 'species-competition' form of the original model in [43] where a term quadratic in E had appeared. This simplifies our starting-point for the ECM equation.) Results for model I are provided in the SI, and will not be the focus of the main text.

Our main modification of model I was to depict more accurate biophysical detail for the ECM signaling via integrin bonds. In the integrin bond model variants, we take into account a force F(R, ρ) exerted on adhesion bonds with basic form

$$\begin{equation}F(R,\rho )={\beta }_{\rho }\frac{\rho }{1+\rho }-{\beta }_{R}\frac{R}{1+R},\end{equation} \tag{ 3.1a }$$

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

$$\begin{equation}\frac{\partial E}{\partial t}={\epsilon}\left(K({E}_{t}-E)-d(F,E)E\right).\end{equation} \tag{ 3.1b }$$

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.

We can interpret E_t as the maximal local density of bound integrin bonds. E_t is related to the available contact area with the ECM, due to, for instance, greater distance between pillars in [42], or any other experimental manipulation that increases ECM binding opportunity. So, higher values of E_t are associated with more available cell-ECM contact.

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]. We first considered the case of slip-bonds (model II). We later determined that experimental data supports the presence of catch-bond [59], so we considered those (model III). We compare the two below, and include model II results in the SI.

3.2.1. Model II: slip bond

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

$$\begin{equation}d(F,E)={k}_{0}\enspace \mathrm{exp}\left(\frac{F}{p(E+{E}_{\mathrm{s}})}\right).\end{equation} \tag{ 3.2 }$$

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 prevents blowup as E → 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₀, whereas when F/E ≈ p, the breakage rate is d ≈ k₀ e¹ ≈ 2.7k₀. 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

$$\begin{equation}d={k}_{0}\enspace \mathrm{exp}\left(\frac{F}{p(E+{E}_{\mathrm{s}})}\right)+{k}_{0c}\enspace \mathrm{exp}\left(-\frac{F}{p(E+{E}_{\mathrm{s}})}\right).\end{equation} \tag{ 3.3 }$$

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.0004, k_0c ≈ 55, appropriate to α₅ β₁ 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 equation (3.1b), together with either of the two force-dependent bond dissociation rates (3.2) or (3.3) (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₀ (level of adhesion at which ECM-signaling to Rho turns on in equation (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 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.

Zoom In Zoom Out Reset image size

Briefly, all three models use the basic Rac-Rho equation (2.2a). Model I used generic ECM dynamics (2.3). 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 in the SI.

We can pinpoint the cause of oscillations that show up later in the full model by first considering the interactions of Rac and Rho on their own, and how feedback from the ECM affects that behavior. Consider the well-mixed version of (2.2a) where spatial gradients are ignored. The inactive Rac and Rho can be eliminated using conservation,

$$\begin{equation*}{R}_{T}=R+{R}_{I},\qquad {\rho }_{T}=\rho +{\rho }_{I},\end{equation*}$$

and the remaining pair of ordinary differential equations (ODEs) for R, ρ, (2.2a) and (2.2b) is a mutual-antagonistic system. This system is bistable for some range of parameter values [58], as shown by the S-shaped curve in 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.

Zoom In Zoom Out Reset image size

Model I is the simplest of the three variants that we studied, facilitating some straightforward experimentation. It is also closest to the original version described in the ODE model of [42, 43]. In the SI, we demonstrate the time-dependent cycles and the bifurcation behavior of the well-mixed variant of model I. We also display preliminary 1D spatial results using model I in SI figure 1.

The three distinct geometries considered are shown in figure 1. (a) A 1D domain, as in figure 1(b). 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–63]. (b) A static 2D circular domain, as illustrated in figure 1(c). 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 1(d). In all these cases, we assume Neumann boundary conditions, since no proteins leak out of the cell edges.

We represent results by kymographs, a popular tool for demonstrating spatial behavior over time. Such plots are commonly used to quantify waves of actin, GTPases, or other cell properties in slices or around the perimeter of a cell [64–67]. Examples are shown in figures 1(e) and 5, 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 −π < x < π around the cell perimeter.

Zoom In Zoom Out Reset image size

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, figure 5(A). Waves that reflect from one edge to the opposite edge show up as a set of peaks and troughs, figure 5(B). 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 Δx/Δt is the wave speed along the boundary, figure 5(C). We also show an example of a random pattern in figure 5(D). 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 11 for an example where all variables are displayed.

In correspondence with Park et al [42], we explore how model behavior changes when γ_E (rate of ECM-feedback to Rho activation) and b_R (basal rate of Rac activation) are varied. We also vary the constant E_t (the maximal adhesion size). Our results, shown as arrays of kymographs for each model, can be compared to the two-parameter bifurcation diagrams in figure 5(d) in [43] (for 'hybrid model 3' in that paper). We concentrate on results for model III below. Similar results for models I and II are given in the SI figures 1 and 2.

Conclusions from the spatial 1D results, including those of model I and II in the supplementary information, are as follows. (1) Such models have overall concurrence with the simplified two-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 behavior in an otherwise polarizable, multistable system [58]. (3) The details of the ECM mechanism, and the ECM equation affect the breadths of parameter regimes, but otherwise do not significantly change the overall conclusions. The three models tested all share regimes of uniform, polar, and oscillatory behaviors 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 7, 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 7(A) 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₀, it signals to activate Rho. If at the same time, the adhesion at the back is lower than E₀, Rho is not sufficiently active to inhibit Rac, allowing Rac activity to invade the rear of the cell, resulting in a front-back polarity flip. If however adhesion at the back exceeds E₀, Rho is upregulated, which inhibits Rac from moving to the rear, maintaining a persistent cell polarity.

Zoom In Zoom Out Reset image size

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

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

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

We carried out parameter sweeps for three values of E_t . We investigate whether the circular geometry in 2D introduces new phenotypes. Results are shown in figure 8 for the catch-slip bond model III, and a comparison to the slip-bond model II is provided in SI figure 4.

Zoom In Zoom Out Reset image size

The catch-slip bond model (III) still promotes substantial regimes of polarized persistent behavior. But it also produces an exotic variety of patterns for E_t = 1000, including standing oscillations and zig-zag waves, along with random-looking transition patterns. See SI figure 3 for a few snapshots of the evolving patterns in figures 8(A)–(D). 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. In most cases, back-front oscillatory patterns are transient and evolve into a single spiral around the cell edge. For higher values of γ_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 (γ_E ). Increasing b_R and decreasing γ_E both lead to higher Rac and thus wider Rac fronts.

Overall, the behavior of model III in 2D shares some features of model II (comparison in SI figure 4). Certain features of the analogous 1D results persist: regimes of low Rac, high Rac, oscillatory, or persistent polar patterns are evident as before. However, the slip-bond (II) version, has a smaller regime of polarization, and more spiral waves.

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 the 1D results. For example, in model III (figure 8) 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.

We next solve the model PDEs over a similar range of parameters, but in a 2D deforming domain, representing the top-down view of the cell shape. In 1D, Buttenschön et al [68] adapted a GTPase model for domain deformation that affects the cell 'volume', leading to dilution that influences the outcomes. Here our 2D simulations CPM is constructed to approximately preserve the total amount of each GTPase in the cell, as explained in the supplementary information. The total volume (area in 2D) of a cell is relatively constant.

Results for model III are shown below, and a comparison with model II behavior is given in SI figure 5. Briefly, domain deformation was tracked using a custom-built cellular Potts model (CPM) calculation with an in-built PDE solver, as described in the SI. (As of this writing, open-source packages such as Morpheus [69] and CompuCell3D [70] 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. [71]), whereas high levels of Rho lead to inwards contraction as in [72], bypassing the explicit representations of actin, myosin, and other cellular components that were included in previous work [26]. 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 supplementary information.

Using the same parameter sweeps as before we show model III results in figure 9, with selected sample cell-shape time sequences in figure 10. (Compare models II and III in SI figure 5.) As before, kymographs track the Rac activity level along the cell edge (dotted white curve in figure 1(d)). Because the domain deforms, and its perimeter length fluctuates slightly with time, the right edge of each kymograph in panels of figure 9 appears slightly ragged.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By comparing patterns in the 2D deforming domain (figure 9) to the identical conditions in the 2D static domain (figure 8), we can identify the effect that domain dynamics has on internal patterns. For the catch-slip bond model (III), 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 γ_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.

Regimes of spiral waves in a static domain tend to become front-to-back oscillations in the deforming domain, especially in the slip-bond model (II), SI figure 5, for E_t = 2000, 3000 where the slanted yellow-purple bands that represent spiral waves in the static case are replaced by the 'honeycomb' kymograph patterns that depict front-back oscillations. For the slip-bond with E_t = 3000, a larger regime of polarization emerges in a band at low γ_E that was previously only spiral waves. The cell deformation appears to dampen spirals in favor of either front-to-back cycles or stable polarity.

The catch-bond model produces much less regular oscillations than the slip-bond model (SI figure 5). 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 (γ_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 γ_E makes the cell more sensitive to weak signals from the ECM, allowing for quick oscillations.

We selected three examples of interesting dynamics (indicated by circles in figure 9) 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 figure 10, 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 figure 9 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 tends to break symmetry, damping the spirals that arrive at the domain boundary. An elongated domain then favors 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.

Among the phenotypes, we see varying degrees of periodicity and spatial regularity in oscillations. In figure 11, we illustrate some time-steps in between front-back switches for the cell labeled B in figure 9, and visualize Rac, Rho and the ECM fields. (See also movies linked to caption in the same figure.) As explained in Figure 7, 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 11 and its accompanying movie. As the cell protrudes to the right, it creates adhesions at the front (175 MCS). At the back of the cell, forces are high due to the presence of Rho. With the catch-bond dynamics, these forces stabilize the adhesion at the back (175 MCS). However, at some point, the forces breach the threshold that breaks adhesions at the back (200 MCS). Then, the ECM-induced Rho activation at the front results in repolarization (225 MCS). Interestingly, some mature adhesions are long-lived, even at the back of the cell (t ⩾ 250 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 phenotypes we observe.

Zoom In Zoom Out Reset image size

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].

While we employed nondimensionalized parameters, SI table 3, we can define a spatial and temporal scales of our ultimate simulations based on properties of melanoma cells in [3, 42]. Details are provided in the SI. Briefly, we use known mean melanoma cell diameter to assign a width (2/3 μm) to a CPM pixel. We then use the approximate speed of persistent melanoma cells, 0.1 μm min⁻¹ [3], to translate Monte Carlo steps to time steps (t₀ = 12s per MCS). We then determine the real time associated with oscillations (roughly one front-back flip in 12 min). In [42], we estimated an oscillating melanoma cell to flip its front every 14 min. So, our time scale of oscillations fits well to this data. We then use the time and spatial scales so determined to quantify the rate of diffusion for inactive GTPases (≈ 6.5 μm² s⁻¹), which is a reasonable ball-park estimate. Hence, while our main focus has been on patterns arising from the PDE models, our results are at appropriate biological scales.

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.