Plasticity of cell migration resulting from mechanochemical coupling
Abstract
Eukaryotic cells can migrate using different modes, ranging from amoeboidlike, during which actin filled protrusions come and go, to keratocytelike, characterized by a stable morphology and persistent motion. How cells can switch between these modes is not well understood but waves of signaling events are thought to play an important role in these transitions. Here we present a simple twocomponent biochemical reactiondiffusion model based on relaxation oscillators and couple this to a model for the mechanics of cell deformations. Different migration modes, including amoeboidlike and keratocytelike, naturally emerge through transitions determined by interactions between biochemical traveling waves, cell mechanics and morphology. The model predictions are explicitly verified by systematically reducing the protrusive force of the actin network in experiments using Dictyostelium discoideum cells. Our results indicate the importance of coupling signaling events to cell mechanics and morphology and may be applicable in a wide variety of cell motility systems.
https://doi.org/10.7554/eLife.48478.001Introduction
Eukaryotic cell migration is a fundamental biological process that is essential in development and wound healing and plays a critical role in pathological diseases, including inflammation and cancer metastasis (Ridley et al., 2003; Roussos et al., 2011; Montell, 2003). Cells can migrate using a variety of modes with a range of corresponding morphologies. The repeated extensions and retractions of pseudopods in amoeboidlike cells, for example, result in a constantly changing morphology and random migration while keratocytelike cells have a stable and broad actinrich front, a nearconstant shape, and move in a persistent fashion (Webb and Horwitz, 2003; Keren et al., 2008). Furthermore, many cells do not have a unique migration mode and can switch between them, either as a function of the extracellular environment or upon the introduction of a stimulus (Paul et al., 2017; Bergert et al., 2012; Charras and Sahai, 2014; Liu et al., 2015; Petrie and Yamada, 2016; Miao et al., 2017). This plasticity is currently poorly understood and is thought to play a role in pathological and physiological processes that involve cell migration, including cancer metastasis (Friedl and Alexander, 2011).
A key step in cell migration is the establishment of an asymmetric and polarized intracellular organization where distinct subsets of signaling molecules, including PAR proteins, Rho family GTPases and phosphoinositides, become localized at the front or back of the cell (Jilkine and EdelsteinKeshet, 2011; Rameh and Cantley, 1999; Goldstein and Macara, 2007; Raftopoulou and Hall, 2004; Rappel and EdelsteinKeshet, 2017). In the absence of directional cues, this symmetry breaking can be a spontaneous and dynamic process with waves of cytoskeletal and signaling components present on the cell cortex (Vicker, 2002; Weiner et al., 2007; Whitelam et al., 2009; Case and Waterman, 2011; Gerisch et al., 2012; Allard and Mogilner, 2013; Gerhardt et al., 2014; Barnhart et al., 2017). Addressing this spontaneous symmetry breaking and the role of waves have generated numerous theoretical studies (Jilkine and EdelsteinKeshet, 2011; Meinhardt, 1999; Otsuji et al., 2007; CsikászNagy et al., 2008; Beta et al., 2008; Mori et al., 2008; Xiong et al., 2010; Knoch et al., 2014; Miao et al., 2017). Most models, however, study cell polarity in the context of biochemical signaling and do not consider cell movement or deformations originated from cell mechanics. This may be relevant for nonmotile cells including yeast (Park and Bi, 2007; Slaughter et al., 2009) but might not be appropriate for motile cells where the coupling between intracellular pathways and cell shape can be crucial in determining the mode of migration (Camley et al., 2013; Camley et al., 2017). Furthermore, most of these models only focus on one specific migration mode and do not address transitions between them. Therefore, it remains an open question how cell mechanics, coupled to a biochemical signaling module, can affect spontaneous cell polarity and can determine transitions between cell migration modes.
Here we propose a novel model that couples an oscillatory biochemical module to cell mechanics. Our choice of the biochemical model was motivated by recent findings that the selforganized phosphatidylinositol (PtdIns) phosphate waves on the membrane of Dictyostelium cells exhibit characteristics of a relaxation oscillator (Arai et al., 2010). Our model is able to generate amoeboidlike, keratocytelike, and oscillatory motion by varying a single mechanical parameter, the protrusive strength, without altering the biochemical signaling pathway. We determine how the transitions depend on these parameters and we show that keratocytelike motion is driven by an emergent traveling wave whose stability is determined by the mechanical properties of the cell. Finally, we experimentally obtain all three migration modes in wildtype Dictyostelium discoideum and explicitly verify model predictions by reducing the actin protrusive force using the drug latrunculin B. Our model provides a unified framework to understand the relationship between cell polarity, motility and morphology determined by cellular signaling and mechanics.
Models and results
Model
Our twodimensional model is composed of two modules: a biochemical module describing the dynamics of an activatorinhibitor system which works in the relaxation oscillation regime, and a mechanical module that describes the forces responsible for cell motion and shape changes (Figure 1a). Our biochemical module consists of a reactiondiffusion system with an activator $A$ (which can be thought of as PtdIns phosphates and thus upstream from newlypolymerized actin; Gerhardt et al., 2014; Miao et al., 2019) and an inhibitor $R$ (which can be thought of as the phosphatase PTEN). This activator and inhibitor diffuse in the cell and obey equations that reproduce the characteristic relaxation oscillation dynamics in the PtdIns lipid system (Arai et al., 2010; Matsuoka and Ueda, 2018; Fukushima et al., 2018):
where ${D}_{A}$ and ${D}_{R}$ are the diffusion coefficients for $A$ and $R$, respectively. In these expressions, $F(A)$ is the selfactivation of the activator with a functional form that is similar to previous studies: $F(A)=[{k}_{a}{A}^{2}/({K}_{a}^{2}+{A}^{2})+b]({A}_{t}A)$ (Miao et al., 2017). The activator is inhibited by $R$ through the negative feedback $G(R)={d}_{1}+{d}_{2}R$ while $R$ is linearly activated by $A$ (see Materials and methods). The timescale of the inhibitor $\tau $ is taken to be much larger than the timescale of the activator, set by $k}_{a$. Finally, to ensure robustness to stochasticity, we add uniformly distributed spatial white noise terms $\u27e8{\zeta}_{i}(\mathbf{\mathbf{r}},t){\zeta}_{j}({\mathbf{\mathbf{r}}}^{\prime},{t}^{\prime})\u27e9=\sigma {\delta}_{ij}\delta (\mathbf{\mathbf{r}}{\mathbf{\mathbf{r}}}^{\prime})\delta (t{t}^{\prime})$.
Nullclines for this system are shown in Figure 1b, where we have chosen parameters such that the fixed point is unstable and the system operates in the oscillatory regime. As a result of the separation of timescales for $A$ and $R$, the dynamics of $A$ and $R$ are characteristic of a relaxation oscillator (inset of Figure 1b): $A$ reaches its maximum quickly, followed by a slower relaxation phase during which the system completes the entire oscillation period.
To generate cell motion, we couple the output of the biochemical model to a mechanical module which incorporates membrane tension and protrusive forces that are proportional to the levels of activator $A$ and normal to the membrane, similar to previous studies (Shao et al., 2010; Shao et al., 2012) (see Materials and methods and Figure 1a). To accurately capture the deformation of the cell in simulations, we use the phase field method (Shao et al., 2010; Ziebert et al., 2012; Shao et al., 2012; Najem and Grant, 2013; Marth and Voigt, 2014; Camley et al., 2017; Cao et al., 2019). Here, an auxiliary field $\varphi $ is introduced to distinguish between the cell interior ($\varphi =1$) and exterior ($\varphi =0$), and the membrane can be efficiently tracked by the contour $\varphi =1/2$. Coupling this field to the reactiondiffusion equations can guarantee that noflux boundary conditions at the membrane are automatically satisfied (Kockelkoren et al., 2003). The evolution of the phasefield is then determined by the force balance equation:
where $\xi $ is a friction coefficient, $\u03f5$ is the boundary width of the phase field, and $H(\varphi )$ is a Hamiltonian energy including the membrane tension, parameterized by $\gamma $ and area conservation (see Materials and methods). The first term on the right hand side describes the actin protrusive force, parameterized by $\eta $, and acts on the cell boundary since $\nabla \varphi $ is nonzero only in a region with width $\u03f5$ formulates the dependence of the protrusive force on the activator levels and is taken to be sigmoidal: $M(A)={A}^{n}/({A}^{n}+{A}_{0}^{n})$, where $n$ is a Hill coefficient. As initial conditions, we use a disk with radius $r$ with area $S=\pi {r}^{2}$ and set $A=R=0$. Default parameter values for our model are estimated from experimental data and given in Table 1.
Cell motion is quantified by computing the velocity of the center of mass of the cell. Furthermore, we record the trajectory of the center of mass and compute its average curvature $\u27e8\kappa \u27e9=\int \kappa (l)\mathit{d}l/L$, where $k(l)$ is the local curvature, and $L$ is the total length of the trajectory. These quantities can be used to distinguish between different migration modes (see Results and Materials and methods).
Computational results
We first examine the possible migration modes as a function of the protrusive strength $\eta $ for fixed area $S$, parametrized by the radius $r$ of the disk used as initial condition, and default parameters. As shown in Figure 2, there are three distinct cell migration modes. When $\eta $ is small, activator waves initiate in the interior and propagate to the cell boundary. However, the protrusive force is too small to cause significant membrane displacement, as also can be seen from the trajectory in Figure 2b. Consequently, the cell is almost nonmotile and the activator and inhibitor field show oscillatory behavior (Figure 2a I and b and Video 1).
As $\eta $ increases, an activator wave that reaches the boundary can create membrane deformations, leading to the breaking of spatial homogeneity. This wave, however, is competing with other traveling waves that emerge from random positions. Consequently, the cell exhibits transient polarity, moves in constantly changing directions, and displays amoeboidlike migration (Figure 2a II and b and Video 2).
When $\eta $ is increased further, protrusions generated by activator waves reaching the cell boundary become even larger. As a result of the coupling between the waves and membrane mechanics, a single traveling wave will emerge within the cell, characterized by a broad and stationary band of high levels of activator. This wave pushes the membrane forward in a persistent direction with constant speed and the cell will adopt a steady keratocytelike morphology, even in the presence of noise (Figure 2a III and b and Video 3).
The transition from oscillatory dynamics to amoeboidlike unstable cell motion can be understood by considering the coupling between the traveling waves and membrane motion. In Materials and methods we show that these traveling waves, that emerge naturally in systems of relaxation oscillators (Kopell and Howard, 1973; Keener, 1980; Murray, 2002) are stable as long as the activator front can ‘outrun’ the inhibitor’s spreading speed. This condition results in a minimal wave speed ${c}_{\text{min}}$ that depends on ${D}_{R}$ and $\tau $ (see Materials and methods). The activator wave pushes the membrane outward and will keep propagating as long as the boundary can keep up with the wave speed. In our model, the membrane is pushed outward by a protrusive force resulting in a speed approximately given by ${v}_{b}\sim \alpha \eta /\xi $, where $\alpha $ is the boundaryaveraged value of $M(A)$, that depends on mechanical parameters and is independent of biochemical parameters (see Materials and methods and Figure 2—figure supplement 1). This implies that if the speed of the membrane is less that the minimum speed of the activator wave (${v}_{b}<{c}_{\text{min}}$) there will be no significant membrane motion. On the other hand, when ${v}_{b}>{c}_{\text{min}}$, a traveling wave can be selected by matching the wave speed and the boundary speed:
The above equation indicates that when $\eta >\xi {c}_{\text{min}}/\alpha $, the cell will break its symmetric shape through traveling waves that deform the membrane.
In simulations, the critical value of $\eta $, ${\eta}_{c,1}$, for which the oscillating cell becomes amoeboidlike can be determined by slowly increasing $\eta $ and computing the centerofmass speed of the cell, ${v}_{CM}$. Oscillatory cells are then defined as cells with a vanishing center of mass speed (see also Materials and methods). The transition between oscillatory and amoeboidlike cells is shown in Figure 2c where we plot ${v}_{CM}$ as a function of protrusion strength $\eta $ for cells with as initial condition either homogeneous solutions in $A$ and $R$ that are perturbed with noise (blue curve) or asymmetric distributions of $A$ and $R$ (red curve). For both initial conditions, the speed shows a subcritical bifurcation at a critical value of ${\eta}_{c,1}$, above which a nonzero cell speed emerges. Furthermore, these simulations reveal that, as argued above, ${\eta}_{c,1}$ depends on both ${D}_{R}$ and $\tau $ (Figure 2—figure supplement 2).
Once a traveling wave is able to generate membrane deformations, why does it not always result in stable, keratocytelike motion with a single traveling wave in the cell’s interior? Notice that in our model, if the spatial extent between the cell front and the back, $d$, is larger than the wavelength of the activator wave $\lambda $ (schematically shown in Figure 1a), a new wave will be generated behind the original one. This wavelength can be approximated by $\lambda \approx \sqrt{2D\tau}$ (where we have taken ${D}_{A}={D}_{R}=D$ for simplicity) such that stable keratocytelike cells driven by a single wave are only possible when $d<\lambda $. This frontback distance, however, depends on the balance between protrusive force and membrane tension at the traveling wave’s lateral ends. Increasing values of $\eta $ result in a broader front and therefore smaller values of $d$ (Figure 2d). As a consequence, as $\eta $ is increased, waves propagating within the cell eventually become stable, resulting in a single, propagating wave and a cell with a keratocytelike morphology. For a discussion on the role of tension we refer to Materials and methods and Figure 2—figure supplement 3.
To distinguish between keratocytelike and amoeboidlike cells we compute the average curvature of the center of mass trajectory $\u27e8\kappa \u27e9$ (see Materials and methods). For keratocytelike cells, which move in a more persistent way, $\u27e8\kappa \u27e9$ will take on small values while for amoeboidlike cells $\u27e8\kappa \u27e9$ will become large. This is shown in Figure 2e where we plot $\u27e8\kappa \u27e9$ as a function of protrusion strength (red curve). The transition from unstable to stable polarity, and therefore keratocytelike motion, happens smoothly in a narrow region of $\eta $ and the critical value ${\eta}_{c,2}$ can be defined as the point for which $\u27e8\kappa \u27e9=0.02\mu {m}^{1}$. Alternatively, and as in the experiments, we can compute the standard deviation of angles of trajectory points, taken at fixed intervals (see Materials and methods). As illustrated by the blue curve in Figure 2e, this standard deviation also decreases rapidly as $\eta $ is increased and thus provides an alternative method to distinguish the two cell migration modes.
Our analysis also implies that for equal protrusive strengths larger cells will be less stable. These cells have a larger area which allows for the nucleation of a new wave front which destabilizes the cell. For equal size cells, those with a smaller protrusion strength would have a larger fontback distance (Figure 2d) and therefore more space to generate a new wave, potentially destabilizing them. Consequently, decreasing the protrusion would destabilize larger keratocytelike cells and only smaller keratocytelike cells will remain.
From the above analysis, it becomes clear that the protrusive strength and size of initial disk, and thus cell area, are critical parameters in determining the stability of the polarity established by interactions between traveling waves and moving boundaries. In simulations, we therefore determine the phase diagram in the $(\eta ,S)$ space by systematically varying the initial radius (with step size 0.5 µm) and protrusive force (with step size 0.25 pN) while keeping all other parameters fixed (Figure 2f). We constrain our cell area to be within the physiologically relevant range with a initial radius between r = 4µm and r = 8µm, corresponding to an area between S ~ 50µm_{2} and S ~ 200µm_{2} (see Materials and methods and Figure 2—figure supplement 4 for an extension of the phase space to larger values of $r$). As stated above, there are three distinct phases, corresponding to the three different cell migration modes of Figure 2a. The transition from oscillatory to amoeboidlike motility occurs at small $\eta $ and is independent of cell size, as predicted in Equation 4. The transition from unstable to stable polarity is quantified by ${\eta}_{c,2}$, which increases for increasing values of $r$ and thus $S$.
The latter transition depends on parameters that affect either $d$ or $\lambda $ (Figure 2g). For example, we can reduce the timescale of the inhibitor $\tau $ to half the value reported in Table 1 (${\tau}_{0}$) which leads to a decrease of $\lambda $ and should lead to larger values of ${\eta}_{c,2}$. Secondly, for increasing values of the membrane tension the transition occurs for larger values of $\eta $. This can be understood by realizing that an increase in the membrane tension will reduce the cell’s deformability. Therefore, the curvature of the cell’s front will decrease, which will increase the frontback distance and thus the critical protrusion strength. Finally, increasing value of the friction coefficient $\xi $ should, according to Equation 4, lead to a decrease in the membrane speed. Since the biochemical wave speed $c$ is unchanged, the transition from amoeboidlike to keratocytelike motion should occur for larger values of $\eta $. All those predictions are confirmed in our simulations, as shown in Figure 2g. Importantly, we find the speed of the keratocytelike cell in all the situations is linearly dependent on $\eta /\xi $, and independent of other parameters, as predicted by Equation 4 (Figure 2h).
In summary, our model predicts that a sufficient decrease of $\eta $ can destabilize keratocytelike cells, resulting in cells that employ amoeboidlike migration, and can transform keratocytelike and unstable cells into oscillatory cells. Furthermore, for decreasing protrusive force, the keratocytelike cells should have a reduced speed and cell size.
Experimental results
To test our model predictions, we carry out experiments using wildtype Dictyostelium cells (see Materials and methods). In vegetative foodrich conditions, during which food is plentiful, most cells migrate randomly using amoeboidlike motion. Starvation triggers cellcell signaling after which cells become elongated and perform chemotaxis. However, we found that starving cells for 6 hr at sufficiently low density is enough to prevent cellcell signaling. Under these conditions, the majority of cells still moves as amoeboidlike cells but a significant fraction, approximately 20–50%, migrates in a keratocytelike fashion (Video 4). These cells adopt a fanshaped morphology and move unidirectionally, as was also observed in certain Dictyostelium mutants (Asano et al., 2004). Employing these low density conditions, we can alter the cell’s protrusive force by interrupting actin polymerization using the drug latrunculin B, an inhibitor of actin activity. Our model predicts that as the concentration of latrunculin increases, keratocytelike cells are more likely to switch to unstable or oscillatory cells. Furthermore, the size and speed of the remaining keratocytelike cells should decrease.
A snapshot of starved cells is shown in Figure 3a, before (left panel) and after exposure to latrunculin (right panel). Higher magnification plots of the amoeboidlike cells (top two panels) and of a keratocytelike cell are shown to the right. In these panels, the actin distribution is visualized with the fluorescent marker limEGFP. As can be seen by comparing the snapshots in Figure 3a, the number of keratocytelike cells decreases after the exposure of latrunculin. This decrease is due to keratocytelike cells becoming unstable and switching to the amoeboidlike mode of migration. We quantify the percentage of keratocytelike cells, as well as the speed and shape, as a function of time for different concentrations of latrunculin for at least 100 cells. As shown in Figure 3b, the percentage of keratocytelike cells decreases upon the introduction of latrunculin. Furthermore, this decrease becomes more pronounced as the concentration of latrunculin is increased (inset Figure 3b), consistent with our model prediction.
To further verify the model predictions, we quantify the cell area $S$ and cell speed for the keratocytelike cells. Both the area (Figure 3c) and speed (Figure 3d) decrease after the introduction of latrunculin. This decrease becomes more significant for larger concentrations of latrunculin, again consistent with our predictions. The effect of latrunculin is further shown in Figure 3—figure supplement 1 where we plot the speed of the amoeboid and keratocytelike cells for different areas as a function of latrunculin concentration. Speed decreases for increasing latrunculin concentration, indicating that the protrusive force is reduced in the presence of latrunculin.
Finally, our model predicts that, with a relaxed area conservation, a sufficient reduction of protrusive strength results in the appearance of oscillatory cells with oscillating basal area size (Figure 2—figure supplement 5). Indeed, after the exposure to latrunculin, a small fraction of cells are observed to display oscillatory behavior characterized by repeated cycle of spreading and contraction, resulting in a basal surface area that oscillates, similar to the engineered oscillatory cells of Miao et al. (2017) (Figure 3e). Cell tracking reveals that these cells originate through a transition from the unstable, amoeboidlike state to the oscillatory state. Interestingly, the observed oscillation in surface area is often very regular (Figure 3e). We find that the average period of different cells is largely independent of the latrunculin concentration, ranging from 6.2 ± 0.9 min for 1 µM to 6.6 ± 0.9 min for 4 µM, while the coefficient of variation, defined as the ratio of the standard deviation and the mean, for single cell periods varies between 0.31 ± 0.12 (1 µM) and 0.20 ± 0.10 (4 µM; Figure 3—figure supplement 2).
Our experimental results can be summarized by identifying the migration mode of cells for different values of the cell area and the latrunculin concentration, thus constructing a phase diagram that can be compared to the computational one (Figure 2f). Due to celltocell variability, we find a distribution of migration modes for each point in this phase diagram. This is shown in Figure 3f where we plot, using a separate color scale for each migration mode, the percentage of keratocytelike, amoeboidlike, and oscillatory mode of migration for different values of the cell area and latrunculin concentration. The resulting experimental phase diagrams agree well with the computational phase diagram presented in Figure 2f. Without latrunculin, many cells migrate using the keratocytelike mode. Exposing cells to latrunculin, and thus reducing the protrusive force, results in a shift from keratocytelike towards amoeboidlike cells. Furthermore, for the maximum value of the latrunculin concentration, almost all keratocytelike cells are destabilized and a small proportion of cells (≈ 9%) are in the oscillatory mode. Thus, our experimental results are in good agreement with the model predictions.
Discussion
In this paper, we propose a simple but unified paradigm to understand cell migration and cell morphology. As in previous modeling studies (Nishimura et al., 2009; Nishimura et al., 2012; Miao et al., 2017; Miao et al., 2019), our model displays different migration modes. These modes can be induced by varying the protrusive force which is attractive since the switching of these modes can occur on a timescale that is shorter than gene expression timescales (Figure 3b), suggesting that a single model with conserved components should be able to capture all three modes. Importantly, our model predictions are verified in experiments using wildtype Dictyostelium cells which, under our conditions, exhibit all three migration modes. These experiments show that upon the introduction of latrunculin the speed of keratocytelike cells decrease. This is perhaps not surprising since latrunculin inhibits actin polymerization which can be expected to result in smaller cell speeds. Our experiments also show, however, that the area of the moving keratocytelike cells decreases in the presence of latrunculin. In addition, and more importantly, our experiments demonstrate that transitions between the migration modes can be brought about by reducing the protrusive strength of the actin network. These nontrivial effects of latrunculin are consistent with the predictions of our model and are also captured in our experimental phase diagrams (Figure 3f).
Key in our model is the coupling of traveling waves generated through biochemical signaling and cell mechanics. Our main finding is that cell migration is driven by the traveling waves and that persistent propagation of these waves result in keratocytelike cells with a broad and stable front. This stable front is only present if the frontback distance is smaller than the biochemical wavelength. Reducing the protrusive strength results in a larger frontback distance, resulting in unstable, amoeboidlike migration. For even smaller values of the protrusive strength, cells display oscillatory behavior. For these values, the membrane speed is smaller than the minimum biochemical wave speed.
Several recent computational studies have addressed switching between different migration modes. For example, Nishimura et al. have presented a model that includes actin and cortical factors, control factors of actin polymerization, and have shown that feedback between cell shape deformations and the spatially distributed control factors can result in amoeboidlike motion (Nishimura et al., 2009; Nishimura et al., 2012) . Furthermore, changing the rate of polymerization as well as the threshold of polymerization in the model can result in transitions between amoeboidlike and keratocytelike cells (Nishimura et al., 2009; Nishimura et al., 2012). In addition, Miao et al. have proposed a model that can generate all three migration modes observed in experiments of engineered Dictyostelium cells (Miao et al., 2017; Miao et al., 2019). This model contains an excitable network and the different migration modes can be generated by altering the threshold of this network. Our current model is distinct from these studies in several ways. First, the migration mode transitions in our model are induced by the mechanical module with the same biochemical components, while in other models the transitions are generated by changing the dynamics of the biochemical signaling pathways from, for example excitable to oscillatory (Miao et al., 2017), or the threshold of actin polymerization (Nishimura et al., 2009; Nishimura et al., 2012). In addition, the biochemical module in our model is much simpler and only contains an activator and inhibitor while the model of Miao et al. (2017) requires additional feedback from a postulated polarity module. Of course, our work does not exclude the existence of this polarity module but it shows that we can explain the observed cell morphologies and movement within a minimal framework of coupling two biochemical components and cell mechanics. Second, based on the earlier measurements of Arai et al. (2010), our model assumes that the biochemical module operates as a relaxation oscillator. As a result, and in contrast to Nishimura et al. (2009) and Nishimura et al. (2012), our model is able to generate oscillatory cells. This is also in contrast to the model of Miao et al. (2017), which uses nested excitable networks. Note however, that in our model we can tune the negative feedback to make the biochemical module operate in the excitable regime. We have explicitly verified that qualitatively similar migration modes and transitions are observed if our model is excitable (Figure 2—figure supplement 6). In the excitable version of our model, however, the oscillations in the nonmotile mode are, in general, less regular and periodic than the ones obtained in the relaxation oscillator version. From the statistical features of the periods obtained from oscillatory cells in experiments, it is likely that the cellular signaling dynamics can be most accurately described by relaxation oscillation models. As a final distinction, we point out that the biochemical and mechanical module in the model of Miao et al. (2017) are solved separately on a 1D ring while in Nishimura et al. (2009) and Nishimura et al. (2012), the biochemical reactions are coupled to cell deformation in a lattice model. As a result, the keratocytelike cells in Miao et al. (2017) display constant excitations at the front that travel along the membrane in the lateral direction rather than stationary activator bands, as observed in the experiments and in our model (Figure 3a(v)). Furthermore, the keratotytelike cells in Nishimura et al. (2009) and Nishimura et al. (2012) have similar shape to the ones generated in our model but have less persistent directionality.
Several future extensions of our study are possible. First, our current study is restricted to two dimensional geometries while actual cell motion is of course three dimensional. Extending our model to 3D would allow us to relax the area conservation constraint and should result in cells for which the basal surface area show clear oscillations that are coupled to extensions away from the surface (Figure 2—figure supplement 5). Second, we have ignored fluid flow within the cytosol, which may play a role in signaling and polarity formation. Including fluid flow, which adds considerable computational complexity to the model (Shao et al., 2012), will be part of future extensions. Third, it should be possible to couple the biochemical model to upstream chemotaxis pathways, allowing it to address directed motion or more complex pathways. In addition, it should be possible to consider multiple parallel and excitable pathways which may regulate cell motility in chemotaxis (Tanabe et al., 2018) and models with more molecular details (Matsuoka and Ueda, 2018; Fukushima et al., 2018). Fourth, it would be interesting to compare wave dynamics obtained in our model with waves observed in giant Dictyostelium cells (Gerhardt et al., 2014). In addition, alternative biochemical models in which parameters determine the qualitatively different dynamics can be studied (Miao et al., 2017). Furthermore, our study predictions may also be verified in other cell types. For example, we predict that overexpression of actin in fast moving cells should result in cells migrating with keratocytelike morphologies while disturbing actin polymerization in keratocytes could lead to unstable migration. Finally, it would be interesting to determine how the feedback between mechanical and biochemical modules can potentially help understand other cell migration processes.
Materials and methods
Full model
Request a detailed protocolOur model for the cell boundary and cell motion is detailed in earlier studies (Shao et al., 2010; Shao et al., 2012; Camley et al., 2014). Briefly, we model the cell boundary as an interface with tension, driven in this study by activator $A$ at the front. Cell motion obeys the overdamped force balance equation ${\mathit{\bm{F}}}_{act}+{\mathit{\bm{F}}}_{mem}+{\mathit{\bm{F}}}_{area}+{\mathit{\bm{F}}}_{fric}=0$ where ${\mathit{\bm{F}}}_{act}$ is the active force proportional to the activator concentration; ${\mathit{\bm{F}}}_{mem}$ describes the membrane tension (line tension since we are modeling a 2D cell); ${\mathit{\bm{F}}}_{area}$ represents area conservation to prevent cells from expanding or shrinking indefinitely and ${\mathit{\bm{F}}}_{fric}$ is a friction force. The active force from the activator is governed by ${\mathit{\bm{F}}}_{act}=\eta M(A)\widehat{\mathit{\bm{n}}}$, where $\widehat{\mathit{\bm{n}}}=\nabla \varphi /\nabla \varphi $ is the outwardpointing normal direction of the membrane, and $M(A)={A}^{3}/({A}^{3}+{A}_{0}^{3})$ where ${A}_{0}$ represents a threshold value for activation of protrusive force. The membrane tension force is computed using the functional derivative (Camley et al., 2014) ${\mathit{\bm{F}}}_{mem}=\frac{\delta {H}_{tension}(\varphi )}{\delta \varphi}\nabla \varphi /{\delta}_{\u03f5}$, with ${\delta}_{\u03f5}=\u03f5{\nabla \varphi }^{2}$ and
Here, $G(\varphi )$ is a double well potential with minima at $\varphi =1$ and $\varphi =0$. As in our earlier work (Camley et al., 2017) we neglect membrane bending and we have verified that it does not qualitatively change the results. We implement area conservation as ${\mathit{\bm{F}}}_{area}={B}_{S}(\int \varphi {d}^{2}\mathit{\bm{r}}{S}_{0})\widehat{\mathit{\bm{n}}}$ where ${B}_{S}$ represents the strength of the area conservation and ${S}_{0}$ is the prescribed area size determined by initial cell radius $r$ through ${S}_{0}=\pi {r}^{2}$. The friction is ${\mathit{\bm{F}}}_{fric}=\xi \mathit{\bm{v}}$ so that $\mathit{\bm{v}}$ is obtained from the force balance equation: $\mathit{\bm{v}}=({\mathit{\bm{F}}}_{act}+{\mathit{\bm{F}}}_{mem}+{\mathit{\bm{F}}}_{area})/\xi $. Note that the friction coefficient takes into account the interaction with the substrate, and that fluid drag can be ignored (Del Alamo et al., 2007). The motion of the phase field $\varphi $ is then determined by the advective equation $\partial \varphi /\partial t=\mathit{\bm{v}}\cdot \nabla \varphi $. Finally, coupling the phase field equations to the reactiondiffusion equations presented in the main text, we arrive at the full equations:
Through the coupling of $\varphi $ to the reactiondiffusion equations, all reaction and diffusion processes are constrained to be inside the auxiliary field (Kockelkoren et al., 2003). The membrane tension parameter is similar to the one used in earlier studies (Shao et al., 2010; Shao et al., 2012; Camley et al., 2014) and taken from Simson et al. (1998). The parameters of the biochemical module are estimated from experiments (Arai et al., 2010; Gerhardt et al., 2014), such that the minimum wave speed in simulations is approximately 0.12µm/s and the wavelength is about 15µm. Note, however, that we can rescale time constants to make simulations more efficient. Therefore, these values are obtained after increasing $\tau $ and time by a factor of 2 and 5, respectively.
Numerical details
Request a detailed protocolThe parameters used for numerical simulations are listed in Table 1. Equations are evolved in a region with size of ${L}_{x}\times {L}_{y}$ = 50 × 50 µm with discrete grids of n×m = 256×256 and periodic boundary conditions are used. Equation 8 is discretized using the forward Euler method with ${\partial}_{t}\varphi =({\varphi}^{(n+1)}{\varphi}^{(n)})/\mathrm{\Delta}t$. Derivatives are calculated using finite difference formulas: ${\partial}_{x}\varphi =({\varphi}_{i+1,j}{\varphi}_{i1,j})/(2\mathrm{\Delta}x)$ and ${\partial}_{x}^{2}\varphi =({\varphi}_{i+1,j}+{\varphi}_{i1,j}2{\varphi}_{i,j})/\mathrm{\Delta}{x}^{2}$, with similar equations for the derivatives in the $y$direction. Equations 6 and 7 are discretized using the forward Euler scheme with ${\partial}_{t}(\varphi A)={\varphi}^{(n)}({A}^{(n+1)}{A}^{(n)})/\mathrm{\Delta}t+{A}^{(n)}({\varphi}^{(n+1)}{\varphi}^{(n)})/\mathrm{\Delta}t$. The diffusion terms $\nabla \cdot (\varphi \nabla A)$ are also approximated using finite difference. The xterm, for example, reads $[({\varphi}_{i+1,j}+{\varphi}_{i,j})({A}_{i+1,j}{A}_{i,j})/(2\mathrm{\Delta}x)({\varphi}_{i,j}+{\varphi}_{i1,j})({A}_{i,j}{A}_{i1,j})/(2\mathrm{\Delta}x)]/\mathrm{\Delta}x$. The white noise terms are simulated as Wiener processes with $\zeta (t)\mathrm{\Delta}t=\sqrt{\sigma \mathrm{\Delta}t}N(0,1)$. As initial condition for $\varphi $, we use a disk $\varphi =\frac{1}{2}[1+\mathrm{tanh}(3(r{r}_{0})/\u03f5)]$, where ${r}_{0}$ is the prescribed radius and $A=R=0$. The activator and inhibitor concentration outside the boundary is 0. To implement this boundary condition, we solve Equations 6 and 7 only in region ${\u03f5}_{0}$ away from $\varphi =1/2$ which is $\varphi >\chi =1/2+1/2\mathrm{tanh}(3{\u03f5}_{0}/\u03f5)\approx 0.0025$, and leave $A=R=0$ outside this region. Here, we have taken ${\u03f5}_{0}=2\mu m$.
Equations are parallelized with CUDA and simulated using GPUs. Typical simulations speeds on a highend graphics board are less than one minute for 100 s of model time.
Identification of migration modes
Request a detailed protocolIn our simulations, we track the motion of the center of mass of the cell which results in a cell trajectory. Oscillatory cells are defined as cells with a vanishing center of mass speed (see also Figure 2c). To distinguish amoeboidlike from keratocytelike cells we can use one of two strategies. The first one is identical to the one used in the experiments and uses the position of the center of mass at discrete time intervals. We then compute the angle between the vector connecting consecutive points and an arbitrary fixed axis and compute the standard deviation of this distribution. Keratocytelike cells, which move more persistently and thus have straighter trajectories, will have a smaller standard deviation than amoeboidlike cells (Figure 2b). The result is shown in Figure 2e where we plot the standard deviation for 100 consecutive intervals, separated by 3 s as a function of protrusion strength (blue curve). Alternatively, taking advantage of the high temporal resolution of the computational tracks, we can use the curvature of the cell trajectory to identify the transition between amoeboidlike and keratocytelike cells. This is also shown in Figure 2e where we plot the average curvature $\u27e8\kappa \u27e9=\int \kappa (l)\mathit{d}l/L$ (red curve) as a function of protrusion strength. Here, $\kappa (l)$ is the local curvature, and $L$ is the total length of the trajectory. We can then define keratocytelike cells as those cells with a standard deviation smaller than 2° or with an average curvature smaller than 0.02 µm^{1}.
To identify the migration mode in the experiments we also track the centroid of the cell. However, since the temporal resolution in the experiments are much lower than that of our simulations, we cannot employ the average curvature identification described above. Instead, we use the angle method and compute angles between two successive positions, separated by 30 s. The standard deviation of this angle is then computed for five consecutive pairs and keratocytelike cells are defined by having a standard deviation less than 25°. To determine whether cells can be considered oscillatory we compute the cell area and evaluate three criteria. First, we compute the coefficient of variation COV (ratio of the standard deviation and the mean) of the area oscillations (computed using five consecutive frames separated by 30 s). Second, we determine the maximum (peaks) and minimum values (valleys) in a time trace of the cell area (see Figure 3e) and quantify the ratio of the difference between the average peak and valley and the average peak: P=(<peak><valley>)/<peak>. Third, we evaluate the total amount of time $T}_{tot$ an oscillation is present. Oscillatory cells are then defined as cells with COV>5%, P>19%, and $T}_{tot$ > 40 min.
Speed of keratocytelike cells
Request a detailed protocolThe local cell boundary velocity can be approximated as $\mathit{\bm{v}}=\widehat{\mathit{\bm{n}}}\eta M(A)/\xi $. Since, for keratocytelike cells, the front is almost flat (see main text Figure 2), the cell speed can be approximated as $v=\alpha \eta /\xi $, where $\alpha $ is the average of $M(A)$ across the boundary: $\alpha ={\int}_{\u03f5}^{\u03f5}M(A)\varphi \mathit{d}x/(2\u03f5)\approx 0.55$. Thus, the speed of keratocytelike cells only depends on the protrusive force and the friction but not on the tension or on the inhibitor timescale. We can verify this in simulations by changing the tension $\gamma $ and the inhibitor’s timescale $\tau $ while keeping other parameters fixed. The results are shown in Figure 2—figure supplement 1 which demonstrates that the cell speed changes little if these parameters are varied.
Traveling waves
Request a detailed protocolOur relaxation oscillator system exhibits traveling waves with a minimum speed that depends on the timescale of the inhibitor and its diffusion constant. Our model can be written as ${\partial}_{t}A=D{\nabla}^{2}A+f(A,R),{\partial}_{t}R=D{\nabla}^{2}R+g(A,R)$, where $f(A,R)$ and $g(A,R)$ can be found from Equations 6 and 7. First we consider the case of $\tau \to \mathrm{\infty}$ and ${D}_{R}=0$ so that the inhibitor is uniformly distributed and constant: $R={R}_{0}$. The relevant equation now is ${\partial}_{t}A=D{\partial}_{x}^{2}A+f(A;{R}_{0})$. When ${R}_{0}$ is in a proper range, there are three steady states $A={A}_{1,2,3}$, with ${A}_{1}<{A}_{2}<{A}_{3}$ and $A={A}_{1,3}$ stable and $A={A}_{2}$ unstable. For a given ${R}_{0}$, we seek for solutions of wave form $A(z)=A(xct)$. The excitable version of this system has been extensively studied and for the cubic reaction term $k(u{u}_{1})({u}_{2}u)(u{u}_{3})$ there is a stable traveling wave solution that connects $u$ with ${u}_{1}$ and ${u}_{3}$ and has a wave speed $c=\sqrt{kD/2}({u}_{1}2{u}_{2}+{u}_{3})$ (Murray, 2002). Likewise, our model ${\partial}_{t}A=D{\partial}_{x}^{2}A+f(A;{R}_{0})$, which has the same structure as the excitable system with a cubic reaction equation, has a stable traveling wave with speed $c\approx w({A}_{1}2{A}_{2}+{A}_{3})$ which depends on ${R}_{0}$ through ${A}_{1,2,3}$. Here, $w$ is a constant that only depends on the diffusion coefficient ${D}_{A}$ and reaction rates (cf. $\sqrt{kD/2}$ for the cubic reaction equation).
For nonzero values of $\tau $ and ${D}_{R}$, the relaxation phase, characterized by the accumulation of inhibitor, sets in behind the activator’s wave front. From the above analysis, it can be deduced that a necessary condition for a stable wave front is a constant profile of the inhibitor $R={R}_{0}$ in the width of the front. For ${D}_{R}=0$, this can only be the case when $\tau $ is large so that the reaction rate of the inhibitor is slow compared to the rate of the activator. Thus, there needs to be a separation of timescales and the system has to obey relaxation dynamics. Furthermore, for ${D}_{R}>0$, the activator wave front is only stable if it can outrun the inhibitor’s spreading speed which is proportional to $\sqrt{{D}_{R}/\tau}$. Therefore, the traveling wave is only stable if both the reaction and diffusion of the inhibitor are slow enough. In other words, stable waves will have a speed that exceeds a minimum value which depends both on ${D}_{R}$ and $\tau $. In simulations, we can change ${c}_{\text{min}}$ by changing the inhibitor’s diffusion coefficient ${D}_{R}$ and the timescale $\tau $. As expected, larger ${D}_{R}$ and smaller $\tau $ leads to larger ${c}_{\text{min}}$, and consequently larger ${\eta}_{c,1}$ (Figure 2—figure supplement 2a and b).
The role of tension in morphology
Request a detailed protocolTension is important to maintain the unidirectional movement of keratocytelike cells. In our model, at the two lateral ends of the traveling wave the morphology is determined by a balance between the protrusive force and the tension which is determined by the local curvature and the parameter $\gamma $. When the protrusive force increases, a larger membrane curvature is necessary to balance the protrusive force, resulting in a flatter front and a decreased frontback distance. If $\gamma $ is not large enough, the traveling wavefront can have a turning instability: the cell will no longer migrate along a straight path and will make a turn (Camley et al., 2017). An example of this instability and the resulting motion is shown in Figure 2—figure supplement 3a where we decrease the surface tension by a factor of 2 (from $\gamma $=2pN/µm to $\gamma $=1pN/µm). Similar to Camley et al. (2017), the cell can also be destabilized by increasing the diffusion of the activator and inhibitor. An example of a simulation showing this can be found in Figure 2—figure supplement 3b.
Parameter variations
Request a detailed protocolWe have examined how the model results change when certain parameters are varied. For example, we can extend the $(\eta ,r)$ phase space to larger values of $r$ as shown in Figure 2—figure supplement 4a. For cell sizes that are beyond ones observed in experiments, we find that the critical protrusive force ${\eta}_{c,2}$ saturates. For these large values of $r$, a dominant wave forms at the front of the cell and new waves that are generated at the back of the cell are not strong enough to break this dominant wave’s persistency. The cell will move persistently in the direction of the dominant wave, with smaller waves repeatedly appearing at the back, as shown in the snapshot presented in Figure 2—figure supplement 4a. For keratocytelike cells with a single wave, $d$ will saturate to the wavelength $\lambda \approx $ 13 µm, as shown in Figure 2—figure supplement 4b. Finally, we have verified that changing the Hill coefficient in $M(A)$ does not change the phase diagram and thus the transitions in a qualitative fashion. This is shown in Figure 2—figure supplement 4c where we plot the boundaries between the different migration modes in the phase diagram for three different values of the Hill coefficient.
Varying area conservation
Request a detailed protocolReal cells are threedimensional objects in which the changing of the basal surface area will be compensated by the morphological changes away from the substrate. Our model represents the cell as a twodimensional object and therefore includes an area conservation term. Making the strength of this area conversation large in simulations allows us to define and sample the $(\eta ,r)$ phase space (Figure 2b). To determine how this area conservation term affects the observed dynamics, we reduce the area conservation parameter from ${S}_{B}=10$ to ${S}_{B}=0.1$. As shown in Figure 2—figure supplement 5, all migration modes are still present. For small values of $\eta $, cells are oscillatory and, compared to large values of ${S}_{B}$, exhibit measurable oscillations in cell size (Figure 2—figure supplement 5a–c). Furthermore, increasing the value of the protrusive strengths results in amoeboidlike motion while even larger values of $\eta $ lead to keratocytelike cells (Figure 2—figure supplement 5d).
Excitable model
Request a detailed protocolOur biochemical model is based on relaxation oscillation dynamics. However, it is straightforward to consider an excitable version of the model. For this, we take ${c}_{2}=30,\sigma =0.1\mu {M}^{2}\mu {m}^{2}/s$ and keep all the other parameters same as listed in Table 1. For the excitable version of our biochemical module, we find similar migration mode transitions as found in the main text as well as a qualitatively similar phase diagram (Figure 2—figure supplement 6a and b). Specifically, we also find a nonmotile, amoeboidlike, and a keratocytelike mode (Figure 2—figure supplement 6a). In the excitable case, however, perturbations are required to initiate waves and movement. As a consequence, the patterns of activator are more noisy than for the oscillation model and the nonmotile cells do not exhibit oscillatory dynamics. Finally, the transition between nonmotile cells and amoeboidlike cells is also subcritical (Figure 2—figure supplement 6c).
Experiments
Request a detailed protocolWildtype AX2 cells were transformed with the plasmid expressing limEdeltacoilGFP. Cells were kept in exponential growth phase in a shaker at 22 °C in HL5 media with hygromycin (50 µg/mL). On the day before the experiment, cells were diluted to low concentration (1–2×10^{5} cells/mL) to stop the exponential growth. After 15 h18h, cell concentration reached 2–5×10^{5} cells/mL and 10^{5} cells were plated in a 50 mm round chamber with glass bottom (WillCo). After 15 min, cells attached to the substrate and HL5 was replaced with 7 mL DB (5 mM Na_{2}HPO_{4}, 5 mM KH_{2}PO_{4},200µM CaCl_{2}, 2 mM MgCl_{2}, pH6.5). Specified amount was diluted in 500µL of DB and added to the sample using a pipette at 6 hr of starvation.
Differential interference contrast (DIC) images are taken every 30 s in six fields of view across the sample using a 10x objective from 5h45 to 6h30 after beginning of starvation. Cell centroids, area, minor and major axis are tracked using Slidebook 6 (Intelligent Imaging Innovations). Statistical analysis of trajectories is performed in MATLAB (2018a; The Mathworks). Experimental data presented before latrunculin B are for cells between 5h45 to 6 hr of starvation, whereas effect of the drug is quantified on cells from 6h15 to 6h30. Speeds are measured using a time interval of 1 min. For each concentration, data are collected on three different days, resulting in N = 333 , N = 315, N = 385, N = 567, and N = 399 for 0, 1, 2, 3, and 4 nM latruculin, respectively. P values are computed with the Wilcoxon rank sum test. Fluorescent images (488 nm excitation) are captured with a 63x oil objective using a spinningdisk confocal Zeiss Axio Observer inverted microscope equipped with a Roper Quantum 512SC cameras.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files.
References

Traveling waves in actin dynamics and cell motilityCurrent Opinion in Cell Biology 25:107–115.https://doi.org/10.1016/j.ceb.2012.08.012

Keratocytelike locomotion in amiBnull Dictyostelium cellsCell Motility and the Cytoskeleton 59:17–27.https://doi.org/10.1002/cm.20015

AdhesionDependent wave generation in crawling cellsCurrent Biology 27:27–38.https://doi.org/10.1016/j.cub.2016.11.011

A bistable mechanism for directional sensingNew Journal of Physics 10:083015.https://doi.org/10.1088/13672630/10/8/083015

Periodic migration in a physical model of cells on micropatternsPhysical Review Letters 111:158102.https://doi.org/10.1103/PhysRevLett.111.158102

Cell motility dependence on adhesive wettingSoft Matter 15:2043–2050.https://doi.org/10.1039/C8SM01832D

Physical influences of the extracellular environment on cell migrationNature Reviews Molecular Cell Biology 15:813–824.https://doi.org/10.1038/nrm3897

Actin and PIP3 waves in giant cells reveal the inherent length scale of an excited stateJournal of Cell Science 127:4507–4517.https://doi.org/10.1242/jcs.156000

PIP3 waves and PTEN dynamics in the emergence of cell polarityBiophysical Journal 103:1170–1178.https://doi.org/10.1016/j.bpj.2012.08.004

The PAR proteins: fundamental players in animal cell polarizationDevelopmental Cell 13:609–622.https://doi.org/10.1016/j.devcel.2007.10.007

Waves in excitable mediaSIAM Journal on Applied Mathematics 39:528–548.https://doi.org/10.1137/0139043

Plane wave solutions to ReactionDiffusion equationsStudies in Applied Mathematics 52:291–328.https://doi.org/10.1002/sapm1973524291

Signaling networks and cell motility: a computational approach using a phase field descriptionJournal of Mathematical Biology 69:91–112.https://doi.org/10.1007/s0028501307044

Orientation of chemotactic cells and growth cones: models and mechanismsJournal of Cell Science 112:2867–2874.

Wave patterns organize cellular protrusions and control cortical dynamicsMolecular Systems Biology 15:e8585.https://doi.org/10.15252/msb.20188585

Bordercell migration: the race is onNature Reviews Molecular Cell Biology 4:13–24.https://doi.org/10.1038/nrm1006

Wavepinning and cell polarity from a bistable reactiondiffusion systemBiophysical Journal 94:3684–3697.https://doi.org/10.1529/biophysj.107.120824

Cortical factor feedback model for cellular locomotion and cytofissionPLOS Computational Biology 5:e1000310.https://doi.org/10.1371/journal.pcbi.1000310

NonBrownian dynamics and strategy of amoeboid cell locomotionPhysical Review E 85:041909.https://doi.org/10.1103/PhysRevE.85.041909

A mass conserved reactiondiffusion system captures properties of cell polarityPLOS Computational Biology 3:e108.https://doi.org/10.1371/journal.pcbi.0030108

Central roles of small GTPases in the development of cell polarity in yeast and beyondMicrobiology and Molecular Biology Reviews 71:48–96.https://doi.org/10.1128/MMBR.0002806

Cancer cell motility: lessons from migration in confined spacesNature Reviews Cancer 17:131–140.https://doi.org/10.1038/nrc.2016.123

Multiple mechanisms of 3D migration: the origins of plasticityCurrent Opinion in Cell Biology 42:7–12.https://doi.org/10.1016/j.ceb.2016.03.025

Cell migration: rho GTPases lead the wayDevelopmental Biology 265:23–32.https://doi.org/10.1016/j.ydbio.2003.06.003

The role of phosphoinositide 3kinase lipid products in cell functionJournal of Biological Chemistry 274:8347–8350.https://doi.org/10.1074/jbc.274.13.8347

Mechanisms of cell polarizationCurrent Opinion in Systems Biology 3:43–53.https://doi.org/10.1016/j.coisb.2017.03.005

Computational model for cell morphodynamicsPhysical Review Letters 105:108104.https://doi.org/10.1103/PhysRevLett.105.108104

Structural mosaicism on the submicron scale in the plasma membraneBiophysical Journal 74:297–308.https://doi.org/10.1016/S00063495(98)777872

Symmetry breaking in the life cycle of the budding yeastCold Spring Harbor Perspectives in Biology 1:a003384.https://doi.org/10.1101/cshperspect.a003384

Transformation from spots to waves in a model of actin pattern formationPhysical Review Letters 102:198103.https://doi.org/10.1103/PhysRevLett.102.198103

Model for selfpolarization and motility of keratocyte fragmentsJournal of the Royal Society Interface 9:1084–1092.https://doi.org/10.1098/rsif.2011.0433
Decision letter

Aleksandra M WalczakReviewing Editor; École Normale Supérieure, France

Detlef WeigelSenior Editor; Max Planck Institute for Developmental Biology, Germany

Jayson PauloseReviewer
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Plasticity of cell migration resulting from mechanochemical coupling" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Jayson Paulose (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
The manuscript addresses the question of how cells switch between different migration modes by coupling mechanics and signalling. All three reviewers are generally positive about the paper. Yet they have certain concerns that would improve both the readability and strengthen the conclusions. We think it is instructive to present the reviewers comments to you, since they are for the most part convergent and they highlight what the reviewers appreciated about this work. Please note that we are aware that the answers to some of the comments are in the supplementary information, yet since the reviewers had a hard time finding them, it is likely that general readers would find this difficult as well.
During the postreview discussion among the reviewers, additional concrete suggestions emerged:
1) It would be helpful if the authors provided a precise definition of the amoeboid > keratocyte transition in the model, and of the amoeboid phase in the experiments, and to move the relevant information (including Figure 2—figure supplement 2A and a similar panel for the amoeboid > keratocyte transition in the model) to the main text and figures.
2) Producing a "migrationmode diagram" from the existing experimental data to provide a more concrete comparison of model prediction to data would be really helpful. The reviewers are not asking for time consuming new experiments.
Reviewer #1:
In this work, the authors address the question of how cells can switch between different migration modes, namely amoeboidlike, characterised by rapid motility induced by unstable protrusions, and keratocytelike, characterised by stable morphology and persistent motion. To this aim, they propose a theoretical model which takes into account (i) a signalling module, made of an activator and an inhibitor, modelled as a reactiondiffusion system; (ii) a mechanical module, modelled by a phasefield method, describing the cell shape deformations as induced by tension, bending, coupling to the signalling module and external friction with the surrounding medium. Their analysis shows that: (i) increasing the coupling between the two modules allows the cells to switch between different modes of migration; (ii) as far as the activator front can outrun the inhibitor's spreading, cells are able to migrate; (iii) while the transition from a nonmigrating oscillatory behaviour to the amoeboidlike mode does not depend on the cell size, that between the amoeboidlike and keratocytelike does. To validate their results, they perform experiments on cells treated with different amounts of an actinpolymerisation inhibitor, which allows to control the feedback between signalling and mechanics. Their findings show that qualitatively the theoretical predictions are recovered. I think that this paper is beautiful and it represents a nice step forward in the general understanding of cell migration and mechanochemical coupling. The manuscript is in general clear and well written, the literature well represented and the authors already addressed in the paper most of the major concerns I would have in principle. Yet, I list below the points I think would definitely improve the clarity and quality of the manuscript.
1) The authors define a reactiondiffusion system made of an activator and inhibitor, respectively called A and R, which works in a relaxation oscillation regime. I do understand that this choice is made in order to observe alternated phases for the activator and the authors discuss their model reduction in the Discussion. Yet, I wonder whether the inhibitor presence is really necessary and a more coarsegrained model could actually reproduce the same results.
2) I have a general feeling that the model results could be better discussed in the main text. Most of the points about the variation of parameters are described in the supplementary materials and the reasoning behind the main text is hard to follow. To make the manuscript clearer, I would also move some of the supplementary figures in the (already existing) main figures by condensing/rearranging the panels and giving them proper titles (e.g., "Protrusive force effect", etc.), once the theory discussion in the main text has been expanded.
3) The authors model cell mechanics with a phase field method. The cell shape is described by an auxiliary field that defines the boundary of the cell changing as a function of mechanochemical forces. I have two main questions about this point. First, I wonder why the authors did not use a more realistic continuous activeviscous shell description. In the context of the phase field model, any viscous dissipation of the membrane is completely neglected and I wonder whether this would only shift the migrationmode switch appearance or would strongly modify the results. Also, in general, the parameters in the phasefield model are quite coarse grained compared to those actually measurable experimentally (e.g., surface tensions). Can the authors clarify this point? Second, the coupling between the mechanical and signalling modules is taken to be strongly nonlinear via a sigmoidal function of the activator. The authors analyse in depth the dependency of the cell migration modes on the intensity of this coupling. I do think it would be worth to investigate also whether and how changing the "Hill coefficient" could lead to motility switches in the model, as this would give more information about the mathematical description of the properties of the mechanochemical coupling in this context and possible cooperative effects at the biochemical level.
4) To my knowledge, there are not universally accepted values for the bending modulus but just very wide ranges of parameters (if I am wrong, I would be very happy to have the reference precisely quantifying this). Could the authors describe whether they investigated the role of this parameter (I have looked into the supplementary materials but did not find any mention about it)? In principle, since bending prevents changes in the mean curvature, I would expect a smoothening or even a disappearance of the switch between the different migration modes.
5) The authors perform experiments on cells where their protrusive force is modulated by interfering with actin polymerisation. Their results all agree qualitatively with their model, i.e., decreasing the protrusive force destabilises keratocytelike migration and allows the switch to the amoeboidlike migration. This has also effects on other quantities as frontback distance, speed and size as qualitatively found by the model. Yet, I do not see anywhere any quantitative comparison between model and experiments. Can the authors clarify this and show an experimental "migrationmode diagram"? A quantitative comparison would be very interesting as it would require in principle a more general definition of protrusive force in the model, as interfering with actin polymerisation alters all mechanical parameters such as tension and bending.
Reviewer #2:
The manuscript by Cao et al. presented a model which incorporates a biochemical module to cell mechanics. The mathematical model produced different migration modes, such as amoeboidlike and keratocytelike motility, depending on strengths of mechanical forces. The authors also evaluated the model by some experiments using Dictyostelium cells. It is a prevailing idea that different migration modes are exchangeable and arise from a common mechanism. This study revealed that a single mechanical parameter, the protrusive strength, enables the transition of different migration modes. However, the author's group has previously constructed similar models for cell boundary and motility. The biochemical model is based on a relaxation oscillator by Arai et al. (2010). Thus, this study likely combined the two existing models. Also, experimental evaluation is not conclusive. My major concerns are described below.
An important finding in this study is that the transition of migration modes depends on the protrusive strength (η) and cell size (r). But it is not clear what is the definition of each migration mode produced by the model. Previous work characterized the amoeboid as well as keratocytelike motility of Dictyostelium cells (Takagi et al., PLoS One, 2008: Asano et al., 2004). Does the model here recapitulate these features of real cell behaviors? Also, I am not sure how the activator is maintained at the front of a keratocytelike cell even after the inhibitor is produced.
The model was evaluated by adding latrunculin to reduce the protrusive strength in Figure 3 but there are some issues to be addressed. First, the criteria for each migration mode is missing in the manuscript.
Second, the keratocytelike cells may differ between the model and the experiment. The model predicted that Factin asymmetrically localize at the cellular front but the keratocytelike cell shown in Figure 3A has Factin at the entire cortex. Also, the model should be evaluated by drawing the phasediagram in the (η,r) space of Dictyostelium cells shown in Figure 3A. The protrusive strength could be quantified by the fluorescence intensity of limEGFP.
Third, the authors claim that the model predicts decreases in cell area size and frontbackdistance of keratocytelike cells after latrunculin treatment in the third paragraph of the subsection “Experimental Results” but I do not understand why. I thought that a decrease in protrusion strength increases frontbackdistance shown in Figure 2—figure supplement 3A.
Fourth, in addition to the study by Miao et al. which is deeply discussed in this manuscript, Nishimura et al. previously published mathematical analysis for the transition of different migration modes (Nishimura, Ueda and Sasai, 2009; 2012). The authors should discuss these works as compared with the current model.
Fifth, the current model still can operate in the excitable regime of the PtdIns phosphate module in Figure 2—figure supplement 6. The recent study showed that multiple pathways including PtdIns phosphate and sGC regulate random cell motility through their distinct excitable properties (Tanabe et al., 2018). Moreover, the model presented by Arai et al. was revised by Fukushima et al. and a new feature of PtdIns phosphate dynamics via the mutual inhibition between PIP3 and PTEN was reported to ensure polarity formation. These points should be discussed so that the presented model could be more realistically considered.
Reviewer #3:
The authors report combined numerical and experimental work on cell motility driven by the coupling between biochemical waves arising from a reactiondiffusion equation and the mechanics of cell shape. The main advance lies in coupling chemical oscillations of an activatorinhibitor system to a mechanical model of the cell: the activator pushes on the cell boundary to change its shape and position, which in turn modifies the spatial domain of the chemical reaction and influences the observed patterns. From simulations and analytical considerations of the waves supported by the reactiondiffusion system, three distinct spatiotemporal patterns are identified which are associated with specific motility and morphology signatures observed in migrating cells. Transitions between these phases are driven by various factors, of which the cell area and the strength of pushing are the most significant. The predictions of the model are tested in an experiment using Dictyostelium cells, in which a reduction of the pushing strength is observed to change the dominant migration mode and the cell morphologies in a manner consistent with the model predictions.
The work appears correct, complete, and wellpresented overall except for the major issues below. As far as I can judge, the mechanism does stand apart from previous works in reducing the complexity and number of components needed to recover the range of cell behaviors observed, and in highlighting the interplay between chemical dynamics and isotropic mechanics (without additional polarization fields). This interplay could be relevant to a variety of other systems which rely on the feedback between reactiondiffusion chemical dynamics and cell mechanics, and the relative simplicity/generality of the model presented here makes it potentially applicable to these other situations. Therefore, I find that the results are significant and of broad interest. I have two major concerns about the methodology which I would like to see addressed before publication.
1) The nature of the transition from amoeboidlike to keratocytelike behavior in the model is not adequately quantified. The transition is described as being from "unstable to stable polarity" and the procedure for classifying experimental trajectories is provided (subsection “Experiments”), but I did not see a rigorous definition of the phase in terms of quantities measured in simulations. How is η_{c,2} determined, and is there a sharp transition in some quantity akin to that seen in the cell speed for the stationary > amoeboid transition (Figure 2—figure supplement 2A)? This information is needed to judge whether the amoeboid > keratocyte transition is indeed a sharp one in the rη plane (as suggested in the Abstract and in Figure 2B), or is rather a smooth crossover.
2) I would like to see more information about how parameter values were chosen in the model. To what extent are the parameter values trying to capture the specifics of the experiments, as opposed to being representative values which qualitatively reproduce the behavior observed in experiments? Specifically, the mechanical parameters are taken from Camley et al. (2014) in subsection “Full model”, but that reference appears to be a computational study rather than an experimental one – were the model parameters used in the earlier study connected to measurements of cell mechanics? If not, what is their source, and are they reasonable? The biochemical parameters are stated as being informed by experiments, but if so, how do the authors explain the large discrepancy between the oscillation period observed in stationary cells (~few minutes, Figure 3—figure supplement 2) and the oscillation period in the simulations (~20 s, Figure 2—figure supplement 5C)?
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Plasticity of cell migration resulting from mechanochemical coupling" for further consideration at eLife. Your revised article has been favorably evaluated by Detlef Weigel as the Senior Editor, a Reviewing Editor, and three reviewers.
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below. I am copying the remaining comments of the reviewers below; these should not be too hard to address.
Reviewer #1:
In the revision, the authors have greatly addressed most of the concerns I raised in the previous review. In my opinion, the manuscript has really improved and it is now much clearer. Yet, I still have two major concerns about the comparison between the experimental and the theoretical migrationmode diagrams and about a check on the forcebalance equation. To be clear, I do really appreciate the effort the authors made in building the experimental migrationmode but I think it would be best if they improved it as far as possible to show a straightforward comparison between theory and experiments.
– I understand that inferring the relationship between the "molecular effect" of latrunculin and the protrusive force by fluorescence analysis is currently difficult or not possible. Yet, I am missing why more quantitative comparisons cannot be tried whatsoever. For example, for a given value of the cell area (averaged over a bin for the experimental data), one could analyse the average velocity of the centre of mass varying latrunculin in the experiments and varying the protrusive force in the theory. I think that this would give a first hint on the relationship between the latrunculin concentration and the protrusive force and how robust this relationship is across different values of the cell area. The same could be carried out for the average curvature. Of course, this is only a suggestion: it could be highly possible that the outcome of this analysis is quite noisy and no relationship can be inferred. In this case, then, I would only focus on merging the three experimental migrationmode diagrams into one and using a thresholding/averaging procedure to make it more uniform and directly comparable to the theoretical one.
– As check, I went through the table of parameters and got a bit confused by the units of measurement of the parameters in the force balance equation.
For example, if I use the units of the friction coefficient reported in the table and plug them into the force balance equation (Equation 8 – Materials and methods), I obtain pN/μm^{3} on the LHS, which does not match the protrusive force term and the area conservation term on the RHS (note I assumed time in s, forces in pN, velocity in μm/s and dimensionless phase field). If I am not missing anything, I think that some of the forces should be rescaled by the disk surface to obtain the "pointbypoint force balance" and Bs, the parameter constraining the cell area, should have dimensions pN/μm^{2}. Could the authors clarify this point or tell me what I am missing? For the concentration module, should not the noise intensity have units μM^{2}/s^{2} (instead of μM^{2}/s)? The other units should be fine instead (but please have a final check).
https://doi.org/10.7554/eLife.48478.020Author response
The manuscript addresses the question of how cells switch between different migration modes by coupling mechanics and signalling. All three reviewers are generally positive about the paper. Yet they have certain concerns that would improve both the readability and strengthen the conclusions. We think it is instructive to present the reviewers comments to you, since they are for the most part convergent and they highlight what the reviewers appreciated about this work. Please note that we are aware that the answers to some of the comments are in the supplementary information, yet since the reviewers had a hard time finding them, it is likely that general readers would find this difficult as well.
During the postreview discussion among the reviewers, additional concrete suggestions emerged:
1) It would be helpful if the authors provided a precise definition of the amoeboid > keratocyte transition in the model, and of the amoeboid phase in the experiments, and to move the relevant information (including Figure 2—figure supplement 2A and a similar panel for the amoeboid > keratocyte transition in the model) to the main text and figures.
We have now provided a precise definition of the transitions between the different migration modes. This is discussed in the main text and in more detail in the new subsection of Materials and methods, “Identification of migration modes”. Also, we have added panels to Figure 2, which now show the transition from oscillatory to amoeboidlike cells (panel C), and from amoeboidlike to keratocytelike cells (panel E).
2) Producing a "migrationmode diagram" from the existing experimental data to provide a more concrete comparison of model prediction to data would be really helpful. The reviewers are not asking for time consuming new experiments.
We thank the editors and reviewers for this suggestion. We have now added such a migrationmode diagram (new Figure 3F) which shows the percentage of cells using the different migration modes in the latrunculin concentration vs. cell area space. Importantly, the qualitative features of these phase diagrams agree very well with the computational phase diagram: large protrusive strength (corresponding to experiments carried out in the absence of latrunculin) results in large numbers of cells migrating in the keratocytelike mode. Introducing latrunculin leads to a decrease in protrusive strength and the destabilization of some of the keratocytelike cells into amoeboidlike cells. For even larger values of latrunculin, we find a larger percentage of cells that show oscillatory behavior, fully consistent with our numerical predictions.
Reviewer #1:
[…] The manuscript is in general clear and well written, the literature well represented and the authors already addressed in the paper most of the major concerns I would have in principle. Yet, I list below the points I think would definitely improve the clarity and quality of the manuscript.
1) The authors define a reactiondiffusion system made of an activator and inhibitor, respectively called A and R, which works in a relaxation oscillation regime. I do understand that this choice is made in order to observe alternated phases for the activator and the authors discuss their model reduction in the Discussion. Yet, I wonder whether the inhibitor presence is really necessary and a more coarsegrained model could actually reproduce the same results.
The coupling between traveling waves and cell mechanics is a crucial element in our model: waves propagate towards the membrane and continue to move as long as the membrane can “keep up” with the wave. Therefore, the ability of the biochemical model to generate propagating waves is necessary. To the best of our knowledge, a minimal model (either excitable or oscillatory) for robust traveling waves with a finite wavelength contains at least two components and these components can be thought of as the activator and inhibitor.
2) I have a general feeling that the model results could be better discussed in the main text. Most of the points about the variation of parameters are described in the supplementary materials and the reasoning behind the main text is hard to follow. To make the manuscript clearer, I would also move some of the supplementary figures in the (already existing) main figures by condensing/rearranging the panels and giving them proper titles (e.g., "Protrusive force effect", etc.), once the theory discussion in the main text has been expanded;
We thank the reviewer for this suggestion. We have substantially rearranged the text and have modified the figures. Specifically, we have added text to the Results section explaining some of the theoretical reasoning and have added 5 new panels to Figure 2. These panels display typical trajectories of our computational cells (panel B) and show the transitions from oscillatory to amoeboidlike cells and from amoeboidlike to keratocytelike cells (panels C and E). Also, we have moved the supplementary figure which shows the effect of protrusion force on keratocytelike cell shape to Figure 2D.
Furthermore, we have moved the supplementary panel showing the effects of parameter variations on the boundaries in the phase diagram to Figure 2G. With these modifications, we hope that the model results are easier to follow.
3) The authors model cell mechanics with a phase field method. The cell shape is described by an auxiliary field that defines the boundary of the cell changing as a function of mechanochemical forces. I have two main questions about this point. First, I wonder why the authors did not use a more realistic continuous activeviscous shell description. In the context of the phase field model, any viscous dissipation of the membrane is completely neglected and I wonder whether this would only shift the migrationmode switch appearance or would strongly modify the results. Also, in general, the parameters in the phasefield model are quite coarse grained compared to those actually measurable experimentally (e.g., surface tensions). Can the authors clarify this point? Second, the coupling between the mechanical and signalling modules is taken to be strongly nonlinear via a sigmoidal function of the activator. The authors analyse in depth the dependency of the cell migration modes on the intensity of this coupling. I do think it would be worth to investigate also whether and how changing the "Hill coefficient" could lead to motility switches in the model, as this would give more information about the mathematical description of the properties of the mechanochemical coupling in this context and possible cooperative effects at the biochemical level.
We appreciate the comment of the reviewer and would like to stress that our goal was to create a model that was as simple as possible and able to capture the critical role of interplay between signaling waves, cell shape and mechanics. Adding fluid flow adds considerable complexity to our computational model and will be our next step. This addition can certainly affect the cell mechanics and might lead to cell moving speed different from Equation 4. However, we believe that it should not strongly modify the main results because migration mode transitions are mostly determined by matching waves to cell shape.
The biochemical parameters (including diffusion constants and reaction rates) are estimated to recapitulate the wave propagation measured in Gerhardt et al. (2014) and Arai et al. (2010). The membrane tension is the same as we have used in our earlier studies (Shao et al., 2010, Shao et al., 2012, Camley et al., 2014) and is estimated from an earlier experimental study (Simson et al., 1998). We now cite this study.
At the suggestion of the reviewer, we have investigated the effect of changing the Hill coefficient, n. These new simulation results are reported as a figure supplement (Figure 2—figure supplement 4C) where we plot the phase diagram of migration modes for different values of n. The simulations reveal that there is no qualitative difference for different values of n. In particular, for fixed cell size there is still a transition between oscillatory and amoeboidlike cells for a critical protrusion strength η_{c,1} and a transition from amoeboidlike to keratocytelike cells for a larger protrusion strength (η_{c,2}). Furthermore, the transition from oscillatory to amoeboidlike cells is independent of cell size for all values of n while η_{c,2} shows the same qualitative behavior as a function of cell size. The shifting of the η_{c,1} and η_{c,2} curves can be simply attributed to the fact that the integrated value of the sigmoidal function M(A) within the cells increases for increasing values of n.
4) To my knowledge, there are not universally accepted values for the bending modulus but just very wide ranges of parameters (if I am wrong, I would be very happy to have the reference precisely quantifying this). Could the authors describe whether they investigated the role of this parameter (I have looked into the supplementary materials but did not find any mention about it)? In principle, since bending prevents changes in the mean curvature, I would expect a smoothening or even a disappearance of the switch between the different migration modes.
We agree with the reviewer that there is experimental uncertainty about the bending modulus. We have followed up on the suggestion of the reviewer and have tested the effects of bending by increasing the modulus 2.5 fold (to 5 pN µm) and by setting it equal to zero. We found that the migration mode transitions are virtually unchanged. This is in line with our previous results (e.g., Camley et al., 2017) and, for simplicity, we have now removed the bending energy from the model.
5) The authors perform experiments on cells where their protrusive force is modulated by interfering with actin polymerisation. Their results all agree qualitatively with their model, i.e., decreasing the protrusive force destabilises keratocytelike migration and allows the switch to the amoeboidlike migration. This has also effects on other quantities as frontback distance, speed and size as qualitatively found by the model. Yet, I do not see anywhere any quantitative comparison between model and experiments. Can the authors clarify this and show an experimental "migrationmode diagram"? A quantitative comparison would be very interesting as it would require in principle a more general definition of protrusive force in the model, as interfering with actin polymerisation alters all mechanical parameters such as tension and bending.
As we already mentioned above, we now include such a migrationmode diagram (Figure 3F) and discuss its results in the main text. This diagram agrees well with the computational predictions.
Reviewer #2:
[…] The author's group has previously constructed similar models for cell boundary and motility. The biochemical model is based on a relaxation oscillator by Arai et al. (2010). Thus, this study likely combined the two existing models. Also, experimental evaluation is not conclusive. My major concerns are described below.
We thank the reviewer for their comments and suggestions. We would like to stress that our current model is more than a simple combination of two models. The reviewer correctly points out that we have presented studies of deforming and moving cells before (e.g., Camley et al., 2013, 2014, Shao et al., 2010 and 2012) and that the biochemical model is based on the important work by Arai et al. However, the novel and most important result of our current study is that coupling between a biochemical module and a mechanical module can produce different cell migration modes.
An important finding in this study is that the transition of migration modes depends on the protrusive strength (η) and cell size (r). But it is not clear what is the definition of each migration mode produced by the model. Previous work characterized the amoeboid as well as keratocytelike motility of Dictyostelium cells (Takagi et al., PLoS One, 2008: Asano et al., 2004). Does the model here recapitulate these features of real cell behaviors?
We thank the reviewer for pointing this out. We now explicitly include a definition of each migration mode in the main text as well as in a new subsection in the Materials and methods section. In simulations, we can define oscillatory cells as cells with a vanishing center of mass velocity while amoeboidlike and keratocytelike cells can be distinguished using either the angles of the cell trajectory or the timeaveraged curvature of this trajectory. Examples of these cell trajectories are now shown in Figure 2B. Since an amoeboidlike cell has an unstable polarity, its trajectory is more curved, resulting in a high average curvature. The trajectory of a keratocytelike cell, on the other hand, is characterized by long and straight segments and thus much lower average curvature. As we now show in new panels Figure 2C and E, these definitions clearly capture the transitions between the different migration modes.
Also, I am not sure how the activator is maintained at the front of a keratocytelike cell even after the inhibitor is produced.
The activator can be maintained because our activatorinhibitor system produces traveling wave. Viewed in a frame moving with the wave speed the front is stable and the inhibitor is always distributed at the back of the activator wave front.
The model was evaluated by adding latrunculin to reduce the protrusive strength in Figure 3 but there are some issues to be addressed. First, the criteria for each migration mode is missing in the manuscript.
We agree with the reviewer and have now added explicit definitions for the criteria for each migration mode.
Second, the keratocytelike cells may differ between the model and the experiment. The model predicted that Factin asymmetrically localize at the cellular front but the keratocytelike cell shown in Figure 3A has Factin at the entire cortex. Also, the model should be evaluated by drawing the phasediagram in the (η,r) space of Dictyostelium cells shown in Figure 3A. The protrusive strength could be quantified by the fluorescence intensity of limEGFP.
Please keep in mind that the variable A in our model does not represent Factin. Rather, it can be thought of as an upstream component such as a PtdIns phosphate. Therefore, the distribution of A is not supposed to exactly correspond to the experimental distribution of Factin. In addition, as can be seen by overlaying the fluorescence image of Figure 3A(v) with the corresponding DIC image, limEGFP is not located at the back membrane of the keratocytelike cell. This is further clarified in Author response image 1 in which we plot the fluorescence image of Figure 3A(v) and the cell outline in red (determined form the DIC image). This image shows that the maximum fluorescence intensity occurs away from the membrane. As shown in experiments (see, e.g., Asano et al., 2008) actin leads PIP3 at the front of the cell while it is trailing PIP3 at the back of the cell. This would imply that the distribution of PtdIns phosphates, including PIP3, will not extend to the membrane but is bounded by the actin distribution. We are currently pursuing further investigations to quantify the distributions of PIP3 and actin in migrating cells.
Reviewer 3, along with the other reviewers and the editors, suggests to construct an experimental phase diagram similar to Figure 2F. Our new panel (Figure 3F) provides such a phase diagram.
We thank the reviewer for the suggestion to quantify the protrusive strength using the fluorescence intensity. Unfortunately, our cells display a relatively large heterogeneity in limEGFP expression levels which makes it not feasible to quantify the protrusive strength.
Third, the authors claim that the model predicts decreases in cell area size and frontbackdistance of keratocytelike cells after latrunculin treatment in the third paragraph of the subsection “Experimental Results” but I do not understand why. I thought that a decrease in protrusion strength increases frontbackdistance shown in Figure 2—figure supplement 3A.
We thank the reviewer for pointing this out and we apologize for the confusion. Our simulations are carried out for constant cell size S_{0}. This assumption, necessary to prevent cells from shrinking or expanding indefinitely and an immediate consequence of the fact that we are modeling a 2D cell, prevents us from making a direct comparison of the frontback distance in simulations and experiments. We have now removed the statement the reviewer was referring to. However, the fact that the cell size in the experiments decreases upon the introduction of latrunculin is consistent with our model. Larger cells have more “room” to initiate a new traveling wave and should therefore be less stable than smaller cells. Thus we can expect that smaller keratocytelike cells will remain fan shaped after a latrunculin treatment while larger keratocytelike cells will become amoeboidlike cells. In the main text, we now only have statement that the cell size will decrease after latrunculin treatment.
Fourth, in addition to the study by Miao et al. which is deeply discussed in this manuscript, Nishimura et al. previously published mathematical analysis for the transition of different migration modes (Nishimura, Ueda and Sasai, 2009; 2012). The authors should discuss these works as compared with the current model.
We thank the reviewer for pointing out these interesting references. We now compare these studies to our model in the Discussion section.
Fifth, the current model still can operate in the excitable regime of the PtdIns phosphate module in Figure 2—figure supplement 6. The recent study showed that multiple pathways including PtdIns phosphate and sGC regulate random cell motility through their distinct excitable properties (Tanabe et al., 2018). Moreover, the model presented by Arai et al. was revised by Fukushima et al. and a new feature of PtdIns phosphate dynamics via the mutual inhibition between PIP3 and PTEN was reported to ensure polarity formation. These points should be discussed so that the presented model could be more realistically considered.
We again thank the reviewer for pointing out these references. They are now included in the revised manuscript, where they are also discussed.
Reviewer #3:
[…] I have two major concerns about the methodology which I would like to see addressed before publication.
1) The nature of the transition from amoeboidlike to keratocytelike behavior in the model is not adequately quantified. The transition is described as being from "unstable to stable polarity" and the procedure for classifying experimental trajectories is provided (subsection “Experiments”), but I did not see a rigorous definition of the phase in terms of quantities measured in simulations. How is η_{c,2} determined, and is there a sharp transition in some quantity akin to that seen in the cell speed for the stationary > amoeboid transition (Figure 2—figure supplement 2A)? This information is needed to judge whether the amoeboid > keratocyte transition is indeed a sharp one in the rη plane (as suggested in the Abstract and in Figure 2B), or is rather a smooth crossover.
We thank Dr. Paulose for this comment, which was also brought up by the other reviewers. As we mentioned in our response to the editors and reviewer 1 and 2, we now include a description of our definitions of the three different migration modes. These definitions are based on the cell trajectories: in simulations, oscillatory cells are defined as cells with a vanishing center of mass velocity while amoeboidlike and keratocytelike cells are defined based on either the average local curvature of the trajectory or based on the standard deviation of the angles between points along this trajectory. In experiments, oscillatory cells are determined by evaluating three criteria, as detailed in Materials and methods. These definitions are also used in the new panels Figure 2C and E which show the transition between oscillatory amoeboidlike cells and between amoeboidlike and keratocytelike cells in simulations, respectively.
2) I would like to see more information about how parameter values were chosen in the model. To what extent are the parameter values trying to capture the specifics of the experiments, as opposed to being representative values which qualitatively reproduce the behavior observed in experiments? Specifically, the mechanical parameters are taken from Camley et al. (2014) in subsection “Full model”, but that reference appears to be a computational study rather than an experimental one – were the model parameters used in the earlier study connected to measurements of cell mechanics? If not, what is their source, and are they reasonable? The biochemical parameters are stated as being informed by experiments, but if so, how do the authors explain the large discrepancy between the oscillation period observed in stationary cells (~few minutes, Figure 3—figure supplement 2) and the oscillation period in the simulations (~20 s, Figure 2—figure supplement 5C)?
Dr. Paulose correctly points out that Camley et al. (2014) is a theoretical study. This study cites other studies which are used to estimate the tension. We apologize for this incomplete referencing and now include the original experimental reference (Simson et al., 1998).
Dr. Paulose is completely correct in noting that the oscillation period in the simulations (~20s) is much shorter (approximately by a factor of 10) than the one observed in experiments (~200s). The oscillation period is determined by the timescale of the inhibitor (τ) as well as the overall time scale of the model. To increase computational efficiency we have chosen τ to be small (10 s). However, rescaling this value by a factor of 2 while rescaling all timedependent equations be a factor of 5 will result in an oscillation period that is 20x10=200s, consistent with the experimental value, and migration speeds and wavelengths that are comparable to experimental values.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below. I am copying the remaining comments of the reviewers below; these should not be too hard to address.
Reviewer #1:
[…]
– I understand that inferring the relationship between the "molecular effect" of latrunculin and the protrusive force by fluorescence analysis is currently difficult or not possible. Yet, I am missing why more quantitative comparisons cannot be tried whatsoever. For example, for a given value of the cell area (averaged over a bin for the experimental data), one could analyse the average velocity of the centre of mass varying latrunculin in the experiments and varying the protrusive force in the theory. I think that this would give a first hint on the relationship between the latrunculin concentration and the protrusive force and how robust this relationship is across different values of the cell area. The same could be carried out for the average curvature. Of course, this is only a suggestion: it could be highly possible that the outcome of this analysis is quite noisy and no relationship can be inferred. In this case, then, I would only focus on merging the three experimental migrationmode diagrams into one and using a thresholding/averaging procedure to make it more uniform and directly comparable to the theoretical one.
We thank the reviewer for this suggestion. We have quantified the average velocity of cells as a function of the latrunculin concentration for various cell areas (new Figure 3—figure supplement 1). This analysis shows that the average speed clearly decreases as the latrunculin concentration increases, indicating that adding latrunculin decreases the protrusive force.
Concerning the suggestion of the reviewer to merge the three experimental migrationmode diagrams into one we respectfully disagree. In our simulations, the migration mode is uniquely determined by initial conditions, resulting in a single phase diagram. In the experiments, on the other hand, we encounter a distribution of different migration modes for a given area and given latrunculin concentration. This makes it difficult to merge the diagrams and therefore we prefer to plot the percentage of cells for the different modes in three separate diagrams.
– As check, I went through the table of parameters and got a bit confused by the units of measurement of the parameters in the force balance equation.
For example, if I use the units of the friction coefficient reported in the table and plug them into the force balance equation (Equation 8 – Materials and methods), I obtain pN/μm^{3} on the LHS, which does not match the protrusive force term and the area conservation term on the RHS (note I assumed time in s, forces in pN, velocity in μm/s and dimensionless phase field). If I am not missing anything, I think that some of the forces should be rescaled by the disk surface to obtain the "pointbypoint force balance" and Bs, the parameter constraining the cell area, should have dimensions pN/μm^{2}. Could the authors clarify this point or tell me what I am missing? For the concentration module, should not the noise intensity have units μM^{2}/s^{2} (instead of μM^{2}/s)? The other units should be fine instead (but please have a final check).
We thank the reviewer for pointing out the unit problem. We have corrected the units in the table as follows: tension γ, pNµm; ξ, pNs/µm; B_{S}: pN/µm^{2}. The noise intensity sigma has units of µM^{2}µm^{2}/s, because the delta function δ(tt’) has units of 1/s while δ(rr’) has units of 1/m^{2}.
https://doi.org/10.7554/eLife.48478.021Article and author information
Author details
Funding
National Science Foundation (PHY1707637)
 WouterJan Rappel
Human Frontier Science Program (LT000371/2017C)
 Elisabeth Ghabache
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Brian A Camley for many useful discussion. This work was supported by the National Science Foundation under grant PHY1707637 and HFSP number LT000371/2017 C.
Senior Editor
 Detlef Weigel, Max Planck Institute for Developmental Biology, Germany
Reviewing Editor
 Aleksandra M Walczak, École Normale Supérieure, France
Reviewer
 Jayson Paulose
Publication history
 Received: May 15, 2019
 Accepted: October 2, 2019
 Version of Record published: October 18, 2019 (version 1)
Copyright
© 2019, Cao et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,968
 Page views

 286
 Downloads

 19
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Physics of Living Systems
Schooling in fish is linked to a number of factors such as increased foraging success, predator avoidance, and social interactions. In addition, a prevailing hypothesis is that swimming in groups provides energetic benefits through hydrodynamic interactions. Thrust wakes are frequently occurring flow structures in fish schools as they are shed behind swimming fish. Despite increased flow speeds in these wakes, recent modeling work has suggested that swimming directly inline behind an individual may lead to increased efficiency. However, only limited data are available on live fish interacting with thrust wakes. Here we designed a controlled experiment in which brook trout, Salvelinus fontinalis, interact with thrust wakes generated by a robotic mechanism that produces a fishlike wake. We show that trout swim in thrust wakes, reduce their tailbeat frequencies, and synchronize with the robotic flapping mechanism. Our flow and pressure field analysis revealed that the trout are interacting with oncoming vortices and that they exhibit reduced pressure drag at the head compared to swimming in isolation. Together, these experiments suggest that trout swim energetically more efficiently in thrust wakes and support the hypothesis that swimming in the wake of one another is an advantageous strategy to save energy in a school.

 Physics of Living Systems
When a fish beats its tail, it produces vortices in the water that other fish could take advantage of to save energy while swimming.