A mathematical model coupling polarity signaling to cell adhesion explains diverse cell migration patterns

Protrusion and retraction of lamellipodia are common features of eukaryotic cell motility. As a cell migrates through its extracellular matrix (ECM), lamellipod growth increases cell-ECM contact area and enhances engagement of integrin receptors, locally amplifying ECM input to internal signaling cascades. In contrast, contraction of lamellipodia results in reduced integrin engagement that dampens the level of ECM-induced signaling. These changes in cell shape are both influenced by, and feed back onto ECM signaling. Motivated by experimental observations on melanoma cells lines (1205Lu and SBcl2) migrating on fibronectin (FN) coated topographic substrates (anisotropic post-density arrays), we probe this interplay between intracellular and ECM signaling. Experimentally, cells exhibited one of three lamellipodial dynamics: persistently polarized, random, or oscillatory, with competing lamellipodia oscillating out of phase (Park et al., 2017). Pharmacological treatments, changes in FN density, and substrate topography all affected the fraction of cells exhibiting these behaviours. We use these observations as constraints to test a sequence of hypotheses for how intracellular (GTPase) and ECM signaling jointly regulate lamellipodial dynamics. The models encoding these hypotheses are predicated on mutually antagonistic Rac-Rho signaling, Rac-mediated protrusion (via activation of Arp2/3 actin nucleation) and Rho-mediated contraction (via ROCK phosphorylation of myosin light chain), which are coupled to ECM signaling that is modulated by protrusion/contraction. By testing each model against experimental observations, we identify how the signaling layers interact to generate the diverse range of cell behaviors, and how various molecular perturbations and changes in ECM signaling modulate the fraction of cells exhibiting each. We identify several factors that play distinct but critical roles in generating the observed dynamic: (1) competition between lamellipodia for shared pools of Rac and Rho, (2) activation of RhoA by ECM signaling, and (3) feedback from lamellipodial growth or contraction to cell-ECM contact area and therefore to the ECM signaling level.

Then s 2,3 are associated with two polarized states while s 4 is an apolar co-existence state. In region I, the only stable solution is s 4 . In regions II and III, solutions s 2 and s 3 respectively are stable. In Regime IV, only s 1 is stable.
2.2. Model 1a. Using R I = ρ I = 1 as before, and taking n = 1, B E = 0 and a linear assumption for b ρ , we arrive at the model equations Technically, the transition from polarized to oscillatory behavior shown in Figure 3 results from a bifurcation to heteroclinic cycles rather than a canonical Hopf bifurcation. (But the conclusions are similar.) In the oscillatory regime of this model, there are two unstable polarized states. The oscillation is a heteroclinic cycle connecting these two states. As the bifurcation is approached from the oscillatory side, the period of oscillation becomes infinitely long, and at the bifurcation the system spends infinite time at one or the other polarized states, representing the transition to stability.

2.4.
Model 1c. This model augments Model 1b to consider the influence of Rac and Rho effector molecules in the protrusion and contraction process. Here w k represents effectors of Rac and c k effectors of Rho in lamellipod k. We assume that these directly influence protrusion and contraction terms (B E and L E ) and that their dynamics are described by linear ODE's. We retain the same Rac and Rho equations as in Model 1b and include the following to describe dynamics of E k , w k , c k .
Here, 1/ 2 is a timescale variable associated with the speed of these reactions. asymptotically reduces to Model 1b provided the GTPase effector dynamics (e.g. ROCK / WASP) are fast ( 2 ). In this case a quasi steady state approximation reduces Model 1c to Model 1b and bifurcation analysis of the full Model 1c equations (not shown) closely matches that of Figure 3. In that case, Model 1c has the same problem that GTPase dynamics are required to be slow to generate oscillations. When 2 , the effectors determine the timescale of feedbacks, and consequently, GTPases need not be slower than ECM dynamics (e.g. it can be the case that < 1).

Model 2
We set R I = ρ I = 1, and Hill coefficients n = 3, and use the assumed b ρ (E k ) expression. Setting A E = 0, and using assumed forms of B E and L E , we get (for k = 1, 2 and j = k):

Model 3
We use the conservation of GTPases combined with the bistable GTPase model (2) to obtain the equations: 4.1. Model 3b. In Figure 5d, we demonstrate the influence of ECM feedback on Rac activation. To do so, we consider a modification of Model 3 given by: Here, the parameter γ ER represents the magnitude of ECM feedback on the rate of Rac activation. Setting γ ER = 0, reduces this variant to the original Model 3.
4.1.1. Model 3b parameter set. All parameters are the same as in Model 3 with γ ER = 0.5. This parameter set was used to generate the dashed borders in Figure 5d. Note that the strength of ECM feedback to Rho activation is γ E = 2. This parameter set represents a scenario where ECM influences both Rac and Rho, but where the influence on Rho is dominant.

Parameter Selection
In our models, a slow negative feedback interacts with a bistable system to generate a relaxation oscillator. This suggests a strategy to find appropriate parameter regimes: first pick parameters for the bistable submodel that provide a reasonable range of behavior (spanning the low monostable, bistable, and high monostable states); then parametrize the negative feedback subsystem so that the full system swings across those three regimes.
For example, in Model 3, bistability stems from the conserved GTPase dynamics. Our experience with such models stems from both dimension-carrying versions [3,4], and from a theoretical analysis of the Rac-Rho system [5]. To first decouple the ECM dynamics, we considered b ρ as a fixed parameter. We reduce the dimensionality of parameter space by scaling the variables. Here we scaled time by the typical GTPase inactivation time 1/δ (equivalent to setting δ = 1 in the full model) and GTPase concentration by the "IC50" level (at which GTPase rate of activation is inhibited down to 50% of its maximal level. (This is equivalent to setting the Hill function IC50 parameters to 1.) Past experience suggests that the system exhibits bistability when the total amount of GTPases is roughly double the value of the IC50 parameter. Thus we chose the dimensionless total concentration of GTPase to be R T = ρ T = 2. Finally, we selected b R and the range of b ρ to span the bistable domain (as noted above) using bifurcation analysis of the Rac / Rho system (Figure 5a). This fully parameterized the GTPase submodel.
We next parameterized the negative feedback subsystem so that the resulting range of values of E i (max and min) suffice to force the parameter b ρ (E i ), now considered a function of E i , to traverse the bistable region in Figure 5a. To do so, k E and γ E were chosen to determine the maximum and minimum values of b ρ . We scaled the ECM variables by the Hill function parameter E 0 (level of ECM that leads to half-maximal ECM-feedback-induced Rho activation rate). Then the ECM feedback in Eqn. (4.1c) switches on and off in the proximity of E j ≈ 1, determining the required range of interest for the variables E j to ensure b ρ (E i ) spans the range [k E , k E + γ E ]. The maximum and minimum values of E i are related to the quantities (k R + γ R )/k ρ and k R /(k ρ + γ ρ ). These are set to ensure that E j 's operate in a reasonable range about the "switching value" of 1. The Hill parameters R 0 and ρ 0 (half-max levels of Rac and Rho for their respective effects on ECM in Eqn. (4.1b)) were adjusted with the same goal in mind, to ensure that (E j ) is able to nearly achieve their required maximum and minimum values. Here we also used the behavior of the decoupled bistable subsystem, namely the maximum and minimum values for Rac and Rho in the two steady states. The midpoint, (max-min)/2, was a reasonable guess for appropriate values of R 0 and ρ 0 . Finally, some small adjustments were made to this basic parameter set to obtain the final results. With this basic simple strategy it was easy to find large parameter swaths that produced the three regimes of behavior and the models were robust to parameter variations.