ORIGINAL RESEARCH article
Sec. Cellular Biochemistry
Two Complementary Signaling Pathways Depict Eukaryotic Chemotaxis: A Mechanochemical Coupling Model
- 1Smart Materials and Advanced Structure Laboratory, School of Mechanical Engineering and Mechanics, Ningbo University, Ningbo Zhejiang, China
- 2Key Laboratory of Microgravity (National Microgravity Laboratory), and Beijing Key Laboratory of Engineered Construction and Mechanobiology, Center for Biomechanics and Bioengineering, Institute of Mechanics, Chinese Academy of Sciences, Beijing, China
- 3State Key Laboratory of Nonlinear Mechanics (LNM) and Beijing Key Laboratory of Engineered Construction and Mechanobiology, Institute of Mechanics, Chinese Academy of Sciences, Beijing, China
- 4School of Engineering Science, University of Chinese Academy of Sciences, Beijing, China
Many eukaryotic cells, including neutrophils and Dictyostelium cells, are able to undergo correlated random migration in the absence of directional cues while reacting to shallow gradients of chemoattractants with exquisite precision. Although progress has been made with regard to molecular identities, it remains elusive how molecular mechanics are integrated with cell mechanics to initiate and manipulate cell motility. Here, we propose a two dimensional (2D) cell migration model wherein a multilayered dynamic seesaw mechanism is accompanied by a mechanical strain-based inhibition mechanism. In biology, these two mechanisms can be mapped onto the biochemical feedback between phosphoinositides (PIs) and Rho GTPase and the mechanical interplay between filamin A (FLNa) and FilGAP. Cell migration and the accompanying morphological changes are demonstrated in numerical simulations using a particle-spring model, and the diffusion in the cell membrane are simulations using a one dimensional (1D) finite differences method (FDM). The fine balance established between endogenous signaling and a mechanically governed inactivation scheme ensures the endogenous cycle of self-organizing pseudopods, accounting for the correlated random migration. Furthermore, this model cell manifests directional and adaptable responses to shallow graded signaling, depending on the overwhelming effect of the graded stimuli guidance on strain-based inhibition. Finally, the model cell becomes trapped within an obstacle-ridden spatial region, manifesting a shuttle run for local explorations and can chemotactically “escape”, illustrating again the balance required in the complementary signaling pathways.
Chemotaxis is characterized by the directed movement of cells along a chemotactic gradient and is crucial for many biological processes, including neuronal patterning (Castellani, 2013), wound healing (Andasari et al., 2018), immune response (Kolaczkowska and Kubes, 2013), and cancer metastasis (Stephens et al., 2008; Polacheck et al., 2013). Even apart from its significant biological function, chemotaxis is an invaluable system for understanding intracellular signaling cascades. For cells to achieve effective movement in response to extracellular guidance, their intracellular signaling network responsible for cytoskeleton remodeling must be differentially and persistently regulated (Petrie et al., 2009; Welf and Haugh, 2011). Studies using modeled organisms, such as neutrophils and Dictyostelium cells, have shown that the molecular mechanisms of chemotaxis are largely shared (Parent, 2004; Van Haastert and Devreotes, 2004), indicating an underlying intrinsic building strategy in physics.
To our knowledge, the spatiotemporal regulation of the biochemical network responsible for eukaryotic chemotaxis can be understood within a dynamic seesaw modeling framework wherein a series of “symmetry breaking” events are involved in different signaling layers (Weiner et al., 2002). At the top layer, cells use uniformly distributed G-protein coupled receptors (GPCRs) to sense the concentration field of a chemoattractant; then, the associated G protein dissociates into Gα and Gβγ subunits, which trigger a fast activation response and a slow inhibition response to achieve directional sensing (Janetopoulos et al., 2004; Levine et al., 2006; Gambardella and Vermeren, 2013). In the middle layer, Rho GTPases (Cdc42, Rac, RhoA) serve as central hubs in transducing signals from the extracellular space to the cytoskeleton (Fukata et al., 2003; Raftopoulou and Hall, 2004). These Rho GTPases are switch-like proteins that are cycled between active membrane-bound (GTP) forms and inactive cytosolic (GDP) forms (Raftopoulou and Hall, 2004), and the existing mutual antagonism effects ensure that the same external signal may determine the features of both the front and the rear of the cell (Onsum and Rao, 2005). At the bottom layer, active Rho GTPases, together with the regulation of phosphatidylinositol-3-kinase (PI3K) and phosphatidylinositol phosphatase (PTEN) (Comer and Parent, 2002; Wang et al., 2002; Billadeau, 2008), contribute to the asymmetric phosphatidylinositol 3, 4, 5-trisphosphate (PIP3) distribution by forming positive/negative feedback loops (Weiner et al., 2002).
Once the spatial effect of such feedback loops has been fully evoked, the cell acts as a seesaw that persists to amplify signaling in a bidirectional manner until an “all-or-none” distribution pattern of molecular components has been achieved (Weiner et al., 2002). For a well-polarized cell, the region of high Cdc42/Rac activity and the activation of Arp2/3 promotes actin branching, thus forming a protrusive pseudopod (Francakoh and Devreotes, 2004; Van Haastert P. J. M., 2010); the region of high RhoA activity and the activation of Rho kinase (ROCK) phosphorylates myosin phosphatase, triggering tail retraction (Wong et al., 2006; Rottner and Stradal, 2011). While the above biochemical components and their interactions are well characterized, the importance of the biomechanical pathways, as represented by FLNa-FilGAP (Nakamura, 2013), is only beginning to be revealed. As the cell extends its pseudopods, FilGAP dissociates from FLNa in actin networks by the myosin-dependent mechanical deformation of the actin network (Nakamura, 2013). FilGAP then relocates to the plasma membrane, where it antagonizes Rac and thus F-actin polarity (Plotnikov et al., 2012). Such a strong dependence between mechanical stretching and FilGAP release has been further proven in reconstitution studies with purified actin cytoskeletal components (Ehrlicher et al., 2011).
A body of mathematical modeling work has been performed for eukaryotic chemotaxis (Devreotes and Janetopoulos, 2003; Danuser et al., 2013). A common feature of these models is the incorporation of diffusion-reaction systems, representing the development of cell polarity, together with a mechanical model of the cell cortex (Holmes and Edelsteinkeshet, 2012). One popular model is based on the level set method (LSM) (Yang et al., 2008). For example, a local excitation, global inhibition biased excitable network (LEGI-BEN) mechanism was incorporated to account for the chemotaxing behavior of cells (Shi et al., 2013). A wave-pinning (WP) mechanism for Rac/Rho signaling was applied to reveal how a cell regulates its shape (Wolgemuth et al., 2010). A layered signaling model was therein proposed to examine how cell shape influences biochemical repolarization in stable keratocyte motion. While these works are helpful as they offer fundamental insights into eukaryotic chemotaxis, the role of biomechanical pathways is largely unknown in those models in the features of chemotaxis.
In this work, we present a mathematical model for eukaryotic chemotaxis with complementary biochemical and biomechanical signaling pathways. We started by using our previous work (Feng et al., 2018), and extended it in two aspects. First, a FLNa-FilGAP pathway was incorporated in a discrete network of actin filaments. Second, the reaction-diffusion in cell membrane is reduced from 2D to 1D for more reasonable and faster computing speed. Our model attempts to incorporate the following experimental findings into an integrative modeling framework: how a cell 1) exhibits correlated random migration in the absence of a chemoattractant stimulus, 2) generates directional and adaptable migration in response to a shallow graded stimulus, and 3) migrates chemotactically within an obstacle-ridden spatial region.
Our 2D model structurally comprises membrane actin and a discrete network of actin filaments (Figure 1A). What we called “the membrane” in the model represents two biological entities: the plasma membrane and the underlying cortex. The membranous viscoelastic element consists of a spring of stiffness k and a dashpot of viscosity η that are connected in parallel. The 2D lamellipod is modeled by N nodes connected by the edges of the Delaunay triangulation. Each lamellipodial element consists of a Hookean spring with stiffness K and a rest length of zero. Our cell migration model strives to capture not only the overall mechanics of cell deformation but also the intracellular signaling events. Therefore, our model can be organized into three modules: the cell mechanics module, biochemistry module, and mechano-sensing module. In the following subsections, we will discuss each of the three modules. As the nucleus and other cytoplasmic organelles play insignificant roles in the phenomenon of interest, they will not be explicitly modeled. Mathematical details and governing assumptions of the model are described in the Supplementary Material.
FIGURE 1. Model components and implications. (A) Two-dimensional cell mesh and the mechanical properties of the elements. The membrane element (blue) consists of a spring of stiffness k and a dashpot of viscosity η that are connected in parallel. The dashpot with a damping coefficient μ connected to the substrate accounts for the viscous dissipation. The lamellipodia are triangulated, and the lamellipodial element consists of an elastic element. FilGAP molecules, marked by colored balls, can either be released from or return to actin-FLNa crosslinking, depending on angular strain. Diffusion-reaction processes are explicitly considered only at the membrane element, implemented by the 1D FDM. (B) Multilayered signaling cascades responsible for cytoskeletal remodeling. 1) Initial signal processing: a cell detects chemoattractants (CAs) by G protein-coupled receptors (GPCRs); the G protein then dissociates into Gα and Gβγ subunits. This layer is closed if no chemotaxis signaling is applied. 2) Small Rho GTPase regulation: spatial regulation of Rac and RhoA by CA is achieved by the localization of GEFs (i.e., PAK1 and Lsc1); releasing FilGAP may downregulate Rac and upregulate RhoA via their antagonistic effects. 3) Bidirectional molecular transport: spatial effects generated from Rho GTPase-PI feedback loops manipulate bidirectional cytoskeletal remodeling. Inserted blue and red boxes show the protrusive and contractive force generation mechanisms, respectively, which are inexplicitly considered. The inserted white plain box shows how such signaling cascades work in physics (a detailed explanation can be found in the discussion).
2.1 Cell Mechanics
The cell membrane is represented by 2D spring-connected particle circle Gallage et al. (2013) (Figure 1A). The motion of each node i on the cell membrane is described by the following force balance equation:
where vi is the velocity vector, μ is the equivalent friction coefficient, and
Cell membrane nodes provide the first boundary condition to lamellipodial nodes, a very good assumption according to the experimental observations that no slippage of the lamellipod is visualized at the leading edge (Paul et al., 2008). The position of the ith lamellipodial node is determined according to the force balance,
where Xj and Xi are the position vector of the edge connecting the ith and jth nodes, and the summation includes all nodes connected to the given node by a Delaunay edge.
2.2 Biochemistry of the Model
The biochemical signaling pathways used in our model are summarized schematically in Figure 1B. Briefly, the model consists of a set of coupled partial differential equations (PDEs) that describe the kinetics, crosstalk, diffusion, and exchange of the signaling molecules and are assigned to three submodules: 1) initial signaling processing, 2) small Rho GTPase regulation, and 3) bidirectional molecular transport. The equations for modules 1) and 3) are in accordance with the balanced inactivation mechanism (Levine et al., 2006) and our previous work (Feng et al., 2018), respectively. The spatiotemporal regulations of active Rac/Rho are described by the following equations:
where G = R and ρ represent the active forms (membrane-bound) of Rac and RhoA, respectively. Gi = Ri, and ρi are the total amounts of the inactive forms of Rac and RhoA, respectively, which can be calculated directly by the conservation law. PG is the activation term and is expressed as
Here, IR and Iρ are the respective baseline activation rates. α and τ are the activation rates of Rac and RhoA by PAK1 (Ep) and Lsc1 (Eρ), respectively, and their concentrations are transmitted from submodule (i). δG is the GAP-mediated baseline inactivation rate. QG represents the FilGAP-mediated antagonism effect between Rac and RhoA and is determined by the first-order ultrasensitivity mechanism as follows (Lipshtat et al., 2010):
Here, AG determined the steepness of this transition, and the parameter b is a measure of the signal impact (Lipshtat et al., 2010).
2.3 Mechano-Sensing Mechanism
The release rate of FLNa-bound FilGAP is mediated by a strain-governed bandpass mechanism (Kang et al., 2015). That is, if an angle change Δθi falls inside range [β1, β2], the FLNa crosslink is in a conformational state of slowly releasing FilGAP; otherwise, the crosslink is in a fast releasing state. The temporal evolution of the number of FilGAPs (m) remaining in a specific network can be described by
where A is a constant and kslow and kfast are slow and fast release constants, respectively. Additionally, sng (Δθi, t) is an indicator function, which is used to switch between slow and fast releasing FilGAP states. The total amount of FilGAP remains conserved throughout the simulations, such that the concentration of cytosolic FilGAP can be derived from the conservation law.
2.4 Computational Setup and Numerical Methods
The cell diameter was set as 10 μm. The mesh parameters used were 421 internal nodes and 45 peripheral nodes for a cell. It is assumed that the cell is chemotactically responsive in a sufficiently large geometry (without boundary effect); otherwise, the first boundary condition is applied to confine cell movement in the presence of obstacles (see Supplementary Material). Since the cell shape and signaling concurrently influence each other, the arising mechano-chemical coupling problem is treated by 1D FDM according to Fick’s laws of diffusion (Merchant et al., 2018a). The diffusion of signaling molecules on the membrane is explicit, while the cytosolic components are assumed to be spatially uniform. A simple Monte-Carlo (MC) method was adopted to evaluate the translocation behaviors of effector molecules (i.e., PI3K/PTEN). Detailed descriptions of the simulation procedures are presented in the Supplementary Material.
3.1 Correlated Random Migration
To illustrate how a cell moves without chemoattractant guidance, we allowed the model cell to move in a simulation run for 900 s with default parameter values (Supplementary Table S1, S2) and repeated the simulations independently for 15 runs. Each simulation was initiated by placing the cell in the same position. The trajectories of cell centroids in each run were recorded and are presented in Figure 2A. The directions of cell migration are completely different among the repeated runs. The dynamics of how the local PIP3 concentration determines the direction of cell movement are presented in Figure 2B (also see Supplementary Video S1). Clearly, once moved, the cell takes a roughly oval-shaped morphology and exhibits a correlated random walk: the direction of future movement is correlated with that of prior movement, and the cell consequently moves with persistence.
FIGURE 2. Simulation of correlated random migration in the absence of an attractant stimulus. (A) Cell centroid trajectories from 15 simulation runs of 900 s each. The cross symbol denotes the starting point. (B) Successive snapshots from a typical simulation. Colored plots are presented of the PIP3 concentrations on the cell membrane and that of FilGAP on the deformed FLNa network. (C,D) Spatiotemporal regulation of PIP3 (C up), PIP2 (C bottom), active Rac (D up), and RhoA (D bottom) during the initial 300 s of the above typical simulation in B. (E,F) Temporal evolution of the cell movement velocity (E) and the cytosolic FilGAP concentration (F).
To interpret the molecular mechanisms of such cellular behaviors, spatiotemporal plots of PIP3/PIP2 and active Rac/RhoA derived from this representative simulation are further demonstrated (Figures 2C,D). Note that the subsequent hotspots appeared at different positions in the PIP3 plot (Figure 2C up), suggesting that the cell generates different pseudopods. More specifically, the distribution of active Rac and RhoA remained uniform during the initial 50 s (Figure 2D). The baseline activities of Rac and RhoA provided a limited “driving force” for evoking cytoskeletal remodeling. In fact, PI3K and PTEN, acting as effectors, were able to translocate from the cytosol to the membrane via PIP3 and PIP2, respectively. After being activated by active Rac and RhoA, the catalytic effect of PI3K leads to the conversion of PIP3 and PIP2 and vice versa for PTEN.
Since membrane diffusion is small, a local increase in PIP3 concentration could be confined locally, resulting in further recruitment of PI3K. A patch of PIP3 is then formed during 30–75 s. As PIP2 has the same diffusivity as PIP3 and their total amounts are constant, the pattern of PIP2 regulation mirrors that of PIP3 (Figure 2C). The regions with high PIP3 and PIP2 are interpreted as the sites where protrusive and contractive force would be enhanced. The accompanying deformation of the actin network also results in the release of FilGAP from the FLNa crosslinks, such that the higher the degree of deformation is, the lower the concentration of FilGAP remains (Figure 2B). After binding with PIP3, FilGAP starts to downregulate the activity of Rac and upregulate that of RhoA, which together results in a local decrease in PIP3 (Figure 2B, 75–175 s). The adjacent regions receive the highest signal amplification and a moderate level of inhibition, defining a narrow range of locations where new pseudopods can be generated.
To further quantify cell motility, temporal evolutions of cell movement velocity and cytosolic FilGAP concentration are presented. The velocity profile displays a transient increase and decrease with a mean value of approximately 4 μm/min (Figure 2E). Moreover, a time lag of 50–100 s exists in the FilGAP response compared with that of velocity (Figure 2E), presumably due to the bandpass mechanism: FilGAP starts to release when the actin network is largely strained and returns to the FLNa network when the angular strain is mostly recovered. Rac and RhoA operate antagonistically through FilGAP, Rac and RhoA signaling and are switched at the leading edge with a time period comparable to that of FilGAP (Figure 2D). Taken together, the basic model reproduces essential aspects of cell motility in the absence of a stimulus with reasonable distributions of the signaling components.
3.2 Model Response to a Gradient Stimulus
To investigate the chemotaxis behaviors of a “wildtype” (WT) cell, a gradient stimulus with default parameter settings was introduced to a spontaneous migrating cell starting at t = 500 s and lasted for an additional 1,600 s. We first tested a case where the initial gradient strength (given in terms of percentage) was only 2% and was applied perpendicular to the polarity the cell had already established (Supplementary Equation S2 describes the percentage as a function of the distance between the cell and the point source). Figure 3A shows the trajectories of cell centroids and a series of successive snapshots of the cell shape colored by PIP3 polarity (also see Supplementary Video S2). Clearly, after receiving this gradient stimulus, the cell gradually changed its migration mode from correlated random migration to directional migration by modulating its PIP3 polarity.
FIGURE 3. Simulation of chemotaxis behaviors in a “wildtype” cell. (A) Trajectory and successive snapshots of cell migration. A gradient stimulus was applied perpendicular to the original migration direction of the cell at t = 500 s when the cell reached the “2” position. The cells were not stimulated on the trajectory located in the white area. (B, C) Spatiotemporal regulation of the molecules responsible for directional sensing, i.e., Gβγ (B, top) and Gα (B, bottom), and for cytoskeletal remodeling, i.e., active Rac (C, top) and PIP3 (C, bottom). (D) Trajectory and successive snapshots of cell migration in response to an opposite gradient stimulus. (E, F) Temporal evolution of the chemotaxis index (CI) in response to different gradient stimuli applied either perpendicular (E) or opposite (F) to the current polarity of the cell. For the definition of CI, please refer to Supplementary Equation S2 in the subsection “Supplementary Material”.
In Figures 3B,C, the spatiotemporal plots of signaling components in a time period before and after receiving directional guidance are further presented. For the balance-inactivation mechanism (Levine et al., 2006), Gβγ and Gα, as activators and inhibitors, respectively, may localize to the membranous regions where the attractant concentration is greatest and weakest, respectively (Figure 3B). The external attractant gradient is thus translated into a spatial bias in Rac activity (Figure 3C, top), which then further regulates PIP3. Specifically, the spatiotemporal plot of PIP3 (Figure 3C, bottom) clearly shows how a new pseudopod evolved from a pre-existing protrusion. Alternatively, we test the case with the same gradient stimulus but in the opposite direction (i.e., at a 180° relative angle). Here, the plots of the derived cell trajectories and successive snapshots show that the cell is not able to reorient itself to the point source (Figure 3D).
By adjusting the relative position of the point source (Supplementary Figure S1A), we further test how sensitive our model is for differentiating gradient strengths. Figure 3E shows the temporal evolution of the chemotaxis index (CI), derived from gradient stimuli, which are applied perpendicularly to the current polarity. Clearly, the cell responds better to the stronger gradients with higher CI values. At the lowest 2% gradient response, the cell movement tends to be aligned but with considerable variability.
We next considered how a spontaneous migrating cell reacts to the gradients from opposite directions. As already illustrated in Figure 3D, the cell is unable to align itself to a 2% gradient, 180° apart (also see Supplementary Video S3). Correspondingly, the derived CI profile starts from 1 and remains at low values (Figure 3F, blue line). Conversely, the profiles increase from the initial value of -1 to maximum value of +1 at higher gradients, indicating that the cell is aligned preferentially along the direction of the gradient. A further comparison is conducted for the cell movement velocity plot (Supplementary Figure S1B). Note that there is a sudden drop in the velocity after receiving a 5% gradient, applied in an opposite manner, indicating that the cell may return to an almost unpolarized state and be repolarized to the external gradient. In other words, the cell performs a reversal rather than a U-turn to react to an opposite gradient stimulus.
3.3 Model Response of Mutant Cells to a Graded Signal
We next considered the chemotaxis behaviors in a number of mutant cells. Applying a 2% graded stimulus perpendicular to a spontaneously migrating FilGAP-abolished cell manifested a complete failure to reorient (Figure 4A and Supplementary Video S4). Spatiotemporal regulation of activator and inhibitor indicated that the cell clearly sensed the stimuli (cf. Figure 4B). However, the amplification of intracellular signaling became too high, causing a freezing state of the initial polarization. In contrast, although a FilGAP-overexpressing cell can respond chemotactically to the same external gradients (Figure 4B and Supplementary Video S5), only a limited number of pseudopods is produced with sluggish motion. The magnitudes of movement velocity and the CI in the late chemotactic stage are presented in the top panel of Figure 4C, wherein the (FilGAP) is varied between these two extremes. Although the cell migrates fast (∼10 μm/min) in the absence of FilGAP, the corresponding CI value is negative, indicating that the cell moves away from the point source. When (FilGAP) increases to 0.12 μM, cellular deformation tends to promote FilGAP release and repress its polarity so that the cell motion is repolarized repeatedly to the external gradient, leading to precise but slow (∼3 μm/min) cell motion. A similar operation is also applied to PI3K mutants by varying (PI3K) from 0 to 0.12 μM while maintaining (FilGAP) at 0.6 μM. As shown in the bottom panel of Figure 4C, the cell is immobilized in the absence of PI3K. In general, the cell movement velocity is increased as (PI3K) increases [the velocity is restricted by the saturation density of F-actin, which is positively correlated with (PIP3)]. Having (PI3K) = 0.12 μM leads to an average velocity of 10.3 μm/min but tends to “freeze” single peaks at the same time, reducing the ability of the cell to respond by turning and dropping off the CI value steeply.
FIGURE 4. Responses of a mutant cell to a gradient stimulus. (A) In the absence of a mechanical inhibition effect (i.e., FilGAP release upon cell deformation), the protrusion front becomes confined, and the cell cannot target itself to the stimulus source (the chemoattractants are released at t = 500 s when the cell reaches the “3” position). (B) With FilGAP overexpression, the cell responds to the graded stimuli precisely with slow velocity. (C) Profiles of average velocity (blue) and the chemotaxis index (red) derived from FilGAP (top) and PI3K mutant cells (bottom).
Taken together, these numerical simulations indicate that the complementary effect of the two signaling pathways is definitely required for a migrating cell that accurately follows the graded stimuli and reaches its target. Consequently, the most effective strategy for chemotaxis is to establish a fine balance between sufficient feedback for autoamplification that enhances the peaks of PI activity and the mechanical inhibition effect that avoids the overly dramatic kurtosis of those peaks.
3.4 Model Responses in the Presence of Obstacles
Under in vivo conditions, chemotactic cells must avoid numerous obstacles, e.g., nonchemotactic cells and foreign debris. While the importance of feedback between shape and biochemistry or the spatial perturbation of chemical field at the obstacle’s vicinity has been well discussed (Grima, 2007; Marée et al., 2012), here we focus on how the effects of complementary signaling pathways enable the cell to cope with such challenge.
In Figure 5A, we illustrate how a “wildtype” cell behaves in a local region surrounded by equally spaced obstacles (also see Supplementary Video S6). In the absence of stimuli (0–20 min), since the cell essentially executes a random walk (Figure 2A), the interactions with obstacles abruptly change its migratory direction and produce a shuttle run. At t = 20 min, the graded stimulus is applied at the outer side of the obstacle cycle. The trajectory plots indicate that the cell smoothly selects a path and navigates from the interval of two obstacles. Figure 5B shows successive snapshots of cell shape changes during spontaneous (points 1, 2 and 3) and chemotactic migration (point 4). When a spontaneously migrating cell encounters obstacles, its leading edge flattens and extends parallel to the wall. Given that the mechanical strain promotes the extra release of FilGAP, this quickly extinguishes the PIP3-rich protrusion fronts and encourages the development of new ones pointing away from the regions of cell-obstacle contact (Figure 5B-2,3). In the presence of chemotactic guidance, however, cell polarity is continually aligned toward the outside of the obstacle region, such that the cell may chemotactically “escape” (Figure 5B-4).
FIGURE 5. Model responses in the presence of obstacles. (A) Cell centroid trajectory during spontaneous (0–500 s) and chemotactic (500–1,500 s) migration within the surrounding steady obstacles (gray circular objects). (B) Snapshots of a cell shape at four time points derived from the above cell trajectory plot. (C) In the absence of a mechanical inhibition effect, the fully polarized cell glides smoothly along the wall of the obstacle that it encounters and escapes from the region between the two adjacent obstacles, leaving a small search area (pink area). (D) In the presence of extremely high inhibition [(FilGAP) ∼1.4 μM], the cell frequently changes its migration direction, achieving a large search area (pink area). (E, F) Temporal evolution of the coverage ratio of FilGAP mutant cells within small (E) and large (F) obstacle surrounding areas. Overexpression of FilGAP can lead to a slow cell moving velocity, resulting in a low searching efficiency, especially regarding a large obstacle surrounding the area.
We next considered how FilGAP mutant cells behave in the presence of obstacles and asked how this translates into biological outcomes. In the absence of FilGAP, as explained in Section “Model response of mutant cells to a graded signal”, the cell becomes fully polarized and migrates in a straight trajectory. When such a cell encounters an obstacle, since no extra FilGAP was released to perturb the polarity already established, the cell smoothly moves along the obstacle wall and squeezes out this local region from the interval of two obstacles (Supplementary Video S7), leaving a small coverage area (Figure 5C, pink area). In the presence of an extremely high mechanical inhibition effect, e.g., (FilGAP) = 1.4 μM, a large coverage area can be achieved because the cell displays the correlated random motion endowed by the obstacle confinement (Figure 5D). To quantify the influence of (FilGAP) on the coverage area, we introduced the coverage ratio, η, defined as the coverage area divided by the total area of the local region. As seen in Figure 5E, the temporal evolution in the absence of FilGAP shows a monotonic increase with time and reaches a maximum at t = 21 s, a time point when the cell escapes from the obstacle area. With an intermediate level of FilGAP, e.g., (FilGAP) = 0.4 μM, the cells move relatively fast, such that the η profile manifests an increase with time, leaving a smaller region to be searched. Note that at later stages, the profile slows down because of the overlap effect. With the overexpression of (FilGAP), e.g., (FilGAP) = 1.2 μM, the cell moves slowly at
In this work, we have proposed a mechanochemical coupling model for eukaryotic chemotaxis. Briefly, the model considers that cell motility is controlled by two signaling pathways. In the first or “biochemical” pathway, the spatial effects of Rho GTPase-PIs mediating positive/negative-feedback loops could be invoked either by chemoattractant stimuli or by endogenous stochastic signals. In the second or “biomechanical” pathway, the interplay of FLNa-FilGAP is incorporated in a discrete network of actin filaments. Elucidating their complementary effects in different conditions obtains a generalized view of eukaryotic chemotaxis. First, the cell uses endogenous signals for the stochastic generation of pseudopods without attractant stimuli. The consequent inhibition effect raised from the strain-mediated FilGAP release ensures the occasional modification of an existing protrusion, allowing the cell to explore the local regions. Second, graded attractant stimuli activate the upstream signaling cascades of Rho GTPase that first induce a positional bias of the pseudopods and then have an overwhelming effect on mechanical inhibition, such that, on average, the extending pseudopods are continually oriented toward the gradient. Third, the cell has become trapped within an obstacle-ridden spatial region, manifesting a shuttle run for local explorations and chemotactically “escaping”. To further explain the above insights, three key issues are discussed below: 1) earlier modeling works and core principles for constructing our model, 2) comparisons with experimental data, and 3) the limitations of our model.
4.1 Earlier Modeling Works and Core Principles for Constructing Our Model
One theoretical difficulty in simulating chemotaxis is how to adequately abstract and simplify the intracellular signaling cascades. Generally, two different approaches exist. One is the “chemical compass” approach, in which the localized intracellular signals at the point of highest receptor occupancy are analogous to a compass needle pointing toward the source of the chemoattractant. Models belonging to this approach are proposed for protein blocking or exogenous activation experiments. For instance, the LEGI mechanism (Janetopoulos et al., 2004) for gradient sensing is based on the fact that the cytoskeleton is frozen by drugs and the wave-pinning mechanism (Mori et al., 2008) on the measurements of exogenous activation of Rac (Lin et al., 2012). Another is the pseudopod-centered approach, which emphasises the endogenous, autocatalytic growth of pseudopods and allows the unprocessed gradient information to bias multiple points in the cycle (Neilson et al., 2011). While persistent effort is being expended in the field of chemotaxis concerning how to integrate the two approaches (Iglesias and Devreotes, 2008; Van Haastert P. J. M., 2010), one similarity is that they both rely on introducing the activity of an inhibitor. Previously, mechanical factors have been shown to act as inhibitors in many theoretical works, yet how they function in conjunction with migratory cell behaviors has not been fully addressed. A pure reaction-diffusion system in a static cell perimeter was proposed where the feedback between Rac and F-actin was combined with the repression of actin polymerization by membrane tension (Wang et al., 2017). A pure mechanical approach was also presented wherein the protrusive rate of the leading edge depends biphasically on the membrane tension upon a force-balance relationship between protrusion and retraction (He and Ji, 2017). Moreover, a mechano-chemical coupling model was proposed for simulating how a neutrophil enters into a narrow channel when the FLNa-FilGAP pathway is implicitly considered (Wu and Feng, 2015). A similar implicit treatment was also adopted by directly introducing hotspots of Rac activity on the cell membrane (Merchant et al., 2018b).
The core principle for constructing our model is to elaborate the above intricate mechanisms into simplified prototypes that can be readily understood. The bottom of Figure 1B shows the bidirectional molecular transport mechanism, which allows a cell to use endogenous signals for the highly ordered extension of pseudopods. Its upstream Rho GTPase regulation and initial signaling processing modules together act as a compass for amplifying extracellular guidance messenger and translate to the bottom module, such that an all-or-none distribution pattern of intracellular signaling components can be achieved. At this stage, the modeled cell can be regarded as a multilayered dynamic seesaw, which only contains biochemical pathways. One recurrent question is that although a cell may easily produce initial polarity and thus directional movement, either by stochastic noise or by guided attractant stimuli, it can only reply to a secondary stimulus with a steeper gradient (Feng et al., 2018). We therefore presume that, in physics, the role of the FLNa-FilGAP pathway is to serve as a brake, as depicted by two springs, into the dynamic seesaw system (Figure 1B, insert). As a result, the seesaw might fluctuate with perturbation, corresponding to the correlated random migration of a cell. However, the dynamic seesaw in the presence of external guidance gradually biases toward the direction of guidance, despite the brake still working, which is essentially analogous to altering the cell movement direction in response to directional cues.
4.2 Comparison With Experimental Data
Our simulations might also provide rational explanations at the molecular level for a set of seemingly controversial measurements at the cellular level. The first set of simulations (Figure 1 and Figure 2) were motivated by the data published by Van Haastert PJ. (2010) and Bosgraaf and Van Haastert (2009), where pseudopodia formation by Dictyostelium cells was analyzed, revealing that pseudopods are frequently formed in two ways. First, pseudopod splitting occurs very frequently, alternating to the right and left at a relatively small angle of ∼55° and leading to a persistent zig-zag trajectory. Second, the cells may also extend pseudopodia from the areas in which the cells are not previously active, described as de novo pseudopodia, inducing a random turn of the cells. In starved Dictyostelium cells, the probability of extending a de novo pseudopod is 10-fold lower than that of pseudopod splitting (Van Haastert PJ., 2010). Upon our simulations in mutant cells, such behaviors could be accounted for as an adaptation strategy for evaluating the efficiency in searching the local environment and could be achieved by adjusting the strength of two signaling pathways. The second set of simulations (Figure 3 and Figure 5) was intended to compare the data published by Andrew and Insall (2007). In their work, the quantitative analysis of pseudopod generation was carried out for chemotaxis in shallow gradients in multiple cell types, including Dictyostelium cells and neutrophils, indicating that 90% of new protrusions were derived from pre-existing ones and rarely emerged from previously pseudopod-free regions of the cell. Upon our simulations, the great difference in the probability of choosing two modes is caused by both the direction and the steepness of the gradient stimulus. Because a cell manifests a proper balance between two signaling pathways, it shows high sensitivity to a shallow gradient stimulus applied from its direction of advancement because such stimuli may continue modifying the existing pseudopod of the cell. However, if the stimulus is applied in an opposite manner, the modification rate may be too slow to guide. In this regard, only a steeper gradient is effective, as its effect is strong enough to produce a new pseudopod directly.
4.3 Limitations of Our Model
First, to avoid introducing too many details at once, our model did not include the signaling cascades for establishing cell-substrate adhesion. While the spatiotemporal regulations of the signaling cascades responsible for cytoskeleton remodeling are widely conserved (Parent, 2004; Van Haastert and Devreotes, 2004), those for cellular adhesion vary significantly from one cell type to another. For instance, rapidly moving cells such as neutrophils (speeds up to 10–2 μm/min) only establish small and transient focal complexes (Holmes and Edelsteinkeshet, 2012). In contrast, slowly moving fibroblasts (speeds of only 1 μm/min) are characterized by strong focal adhesion (Holmes and Edelsteinkeshet, 2012). The establishment of cell-substrate adhesion itself is also a mechano-sensing process (Bershadsky et al., 2006; Giannone and Sheetz, 2006). Further incorporating the models with the assembly and disassembly of cell-substrate adhesion enables our model to provide a possible explanation regarding why there are unique features between different cells.
Second, the current model tends to focus on explaining transient signaling activity on the membrane, for the production and competition of pseudopods, and regarding the role of pseudopod dynamics in chemotaxis. Accordingly, the details of the cytoskeletal remodeling process are omitted, and the signaling activity is implicitly linked to forces causing protrusion and retraction. Cytoskeletal remodeling is likely governed by the polarities of signaling molecule distributions (refer to Supplementary Material), but their intercorrelation is far more sophisticated. A variety of mechanisms, including G-actin treadmilling, dendritic nucleation at the leading edge, and vesicle delivery by microtubules, may function together to regulate the dynamic patterns of actin-myosin flow (Rubinstein et al., 2009; Wolgemuth et al., 2010).
Finally, our model is a simple but efficient numerical tool for mechano-chemical coupling problems. This method could be applicable in the field of collective cell migration (Feng et al., 2018), wherein the position-dependent cell directionality and speed are of greatest concern. Further improvements could be achieved by incorporating the particle-spring model with the immersed boundary (IB) method (Dallon et al., 2009; Wu and Feng, 2015), such that the integrated method could be used to decipher biological processes involving fluid-solid interactions. Nevertheless, the efficient treatment of signaling cascades in our model also requires an accurate account of the cell morphology. In this regard, the CPM method is a more suitable platform (Marée et al., 2012).
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
ML designed the research. LZ and SF carried out all simulations and analyzed the data. ML, SF, LZ, SL, YZ, and LL wrote the article. All authors discussed the results and contributed to the final manuscript.
This work was supported by National Natural Science Foundation of China grants 12002176, 32130061, 11972200 and 11902327, Opening fund of State Key Laboratory of Nonlinear Mechanics (LNM) of Chinese Academy of Sciences, Frontier Science Key Project of Chinese Science Academy grant QYZDJ-SSW-JSC018, and Youth Innovation Promotion Association CAS.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.786254/full#supplementary-material
Andasari, V., Lü, D., Swat, M., Feng, S., Spill, F., Chen, L., et al. (2018). Computational Model of Wound Healing: Egf Secreted by Fibroblasts Promotes Delayed Re-Epithelialization of Epithelial Keratinocytes. Integr. Biol. 10, 605–634. doi:10.1039/c8ib00048d
Andrew, N., and Insall, R. H. (2007). Chemotaxis in Shallow Gradients Is Mediated Independently of Ptdins 3-Kinase by Biased Choices Between Random Protrusions. Nat. Cell Biol. 9, 193–200. doi:10.1038/ncb1536
Bershadsky, A., Kozlov, M., and Geiger, B. (2006). Adhesion-Mediated Mechanosensitivity: a Time to Experiment, and a Time to Theorize. Curr. Opin. Cell Biol. 18, 472–481. doi:10.1016/j.ceb.2006.08.012
Dallon, J. C., Newren, E., and Hansen, M. D. (2009). Using a Mathematical Model of Cadherin-Based Adhesion to Understand the Function of the Actin Cytoskeleton. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 79, 031918. doi:10.1103/PhysRevE.79.031918
Danuser, G., Allard, J., and Mogilner, A. (2013). Mathematical Modeling of Eukaryotic Cell Migration: Insights Beyond Experiments. Annu. Rev. Cell Dev. Biol. 29, 501–528. doi:10.1146/annurev-cellbio-101512-122308
Ehrlicher, A. J., Nakamura, F., Hartwig, J. H., Weitz, D. A., and Stossel, T. P. (2011). Mechanical Strain in Actin Networks Regulates Filgap and Integrin Binding to Filamin a. Nature 478, 260–263. doi:10.1038/nature10430
Feng, S., Zhou, L., Zhang, Y., Lü, S., and Long, M. (2018). Mechanochemical Modeling of Neutrophil Migration Based on Four Signaling Layers, Integrin Dynamics, and Substrate Stiffness. Biomech. Model. Mechanobiol. 17, 1611–1630. doi:10.1007/s10237-018-1047-2
Janetopoulos, C., Ma, L., Devreotes, P. N., and Iglesias, P. A. (2004). Chemoattractant-Induced Phosphatidylinositol 3,4,5-trisphosphate Accumulation Is Spatially Amplified and Adapts, Independent of the Actin Cytoskeleton. Pnas. 101, 8951–8956. doi:10.1073/pnas.0402152101
Kang, J., Puskar, K. M., Ehrlicher, A. J., Leduc, P. R., and Schwartz, R. S. (2015). Structurally Governed Cell Mechanotransduction Through Multiscale Modeling. Sci. Rep. 5, 8622. doi:10.1038/srep08622
Lin, B., Holmes, W. R., Wang, C. J., Ueno, T., Harwell, A., Edelstein-Keshet, L., et al. (2012). Synthetic Spatially Graded Rac Activation Drives Cell Polarization and Movement. Proc. Natl. Acad. Sci. U S A. 109, E3668–E3677. doi:10.1073/pnas.1210295109
Lipshtat, A., Jayaraman, G., He, J. C., and Iyengar, R. (2010). Design of Versatile Biochemical Switches that Respond to Amplitude, Duration, and Spatial Cues. Proc. Natl. Acad. Sci. 107, 1247–1252. doi:10.1073/pnas.0908647107
Marée, A. F. M., Grieneisen, V. A., and Edelstein-Keshet, L. (2012). How Cells Integrate Complex Stimuli: the Effect of Feedback From Phosphoinositides and Cell Shape on Cell Polarization and Motility. Plos Comput. Biol. 8, e1002402. doi:10.1371/journal.pcbi.1002402
Merchant, B., Edelstein-Keshet, L., and Feng, J. J. (2018a). A Rho-Gtpase Based Model Explains Spontaneous Collective Migration of Neural Crest Cell Clusters. Developmental Biol. 444, S262–S273. doi:10.1016/j.ydbio.2018.01.013
Merchant, B., Edelstein-Keshet, L., and Feng, J. J. (2018b). A Rho-Gtpase Based Model Explains Spontaneous Collective Migration of Neural Crest Cell Clusters. Dev. Biol. 444 Suppl 1, S262. doi:10.1016/j.ydbio.2018.01.013
Neilson, M. P., Veltman, D. M., van Haastert, P. J. M., Webb, S. D., Mackenzie, J. A., and Insall, R. H. (2011). Chemotaxis: A Feedback-Based Computational Model Robustly Predicts Multiple Aspects of Real Cell Behaviour. Plos Biol. 9, e1000618. doi:10.1371/journal.pbio.1000618
Paul, R., Heil, P., Spatz, J. P., and Schwarz, U. S. (2008). Propagation of Mechanical Stress Through the Actin Cytoskeleton Toward Focal Adhesions: Model and experiment. Biophysical J. 94, 1470–1482. doi:10.1529/biophysj.107.108688
Plotnikov, S. V., Pasapera, A. M., Sabass, B., and Waterman, C. M. (2012). Force Fluctuations Within Focal Adhesions Mediate Ecm-Rigidity Sensing to Guide Directed Cell Migration. Cell 151, 1513–1527. doi:10.1016/j.cell.2012.11.034
Rubinstein, B., Fournier, M. F., Jacobson, K., Verkhovsky, A. B., and Mogilner, A. (2009). Actin-Myosin Viscoelastic Flow in the Keratocyte Lamellipod. Biophysical J. 97, 1853–1863. doi:10.1016/j.bpj.2009.07.020
Shi, C., Huang, C. H., Devreotes, P. N., and Iglesias, P. A. (2013). Interaction of Motility, Directional Sensing, and Polarity Modules Recreates the Behaviors of Chemotaxing Cells. Plos Comput. Biol. 9, e1003122. doi:10.1371/journal.pcbi.1003122
Wang, F., Herzmark, P., Weiner, O. D., Srinivasan, S., Servant, G., and Bourne, H. R. (2002). Lipid Products of Pi(3)ks Maintain Persistent Cell Polarity and Directed Motility in Neutrophils. Nat. Cell Biol. 4, 513–518. doi:10.1038/ncb810
Wang, W., Tao, K., Wang, J., Yang, G., Ouyang, Q., Wang, Y., et al. (2017). Exploring the Inhibitory Effect of Membrane Tension on Cell Polarization. Plos Comput. Biol. 13, e1005354. doi:10.1371/journal.pcbi.1005354
Weiner, O. D., Neilsen, P. O., Prestwich, G. D., Kirschner, M. W., Cantley, L. C., and Bourne, H. R. (2002). A PtdInsP3- and Rho GTPase-Mediated Positive Feedback Loop Regulates Neutrophil Polarity. Nat. Cell Biol. 4, 509–513. doi:10.1038/ncb811
Wong, K., Pertz, O., Hahn, K., and Bourne, H. (2006). Neutrophil Polarization: Spatiotemporal Dynamics of Rhoa Activity Support a Self-Organizing Mechanism. Proc. Natl. Acad. Sci. 103, 3639–3644. doi:10.1073/pnas.0600092103
Keywords: chemotaxis, cytoskeletal remodeling, mathematical model, biochemical, biomechanical
Citation: Zhou L, Feng S, Li L, Lü S, Zhang Y and Long M (2021) Two Complementary Signaling Pathways Depict Eukaryotic Chemotaxis: A Mechanochemical Coupling Model. Front. Cell Dev. Biol. 9:786254. doi: 10.3389/fcell.2021.786254
Received: 30 September 2021; Accepted: 21 October 2021;
Published: 17 November 2021.
Edited by:Jinglei Hu, Nanjing University, China
Reviewed by:Yuantong Gu, Queensland University of Technology, Australia
Jianshan Wang, Tianjin University, China
Copyright © 2021 Zhou, Feng, Li, Lü, Zhang and Long . This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†These authors have contributed equally to this work