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: μ v i = F i elas + F i vis + F i pro + F i cont + F i drag (1) where v _i is the velocity vector, μ is the equivalent friction coefficient, and F i elas , F i vis , F i pro , F i cont and F i drag are the elastic energy force, viscous energy force, protrusive force, contractive force, and drag force, respectively. Expressions for each term are provided in the Supplementary Material.

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, ∑ j X j − X i = 0 (2) where X _j and X _i 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.

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: ∂ G ∂ t = D m ∇ 2 G + P G ( R , ρ ) G i − δ G G + Q G G , (3) where G = R and ρ represent the active forms (membrane-bound) of Rac and RhoA, respectively. G _i = R _i , and ρ _i are the total amounts of the inactive forms of Rac and RhoA, respectively, which can be calculated directly by the conservation law. P _G is the activation term and is expressed as P R = I R + α E P and P ρ = I ρ + τ E ρ . (4) Here, I _R and I _ρ are the respective baseline activation rates. α and τ are the activation rates of Rac and RhoA by PAK1 (E _p ) and Lsc1 (E _ρ ), respectively, and their concentrations are transmitted from submodule (i). δ _G is the GAP-mediated baseline inactivation rate. Q _G 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): Q R = − Q ρ = − 1 1 + A G ⁡ exp − b [ FilGAP ] ⋅ [ PIP 3 ] (5) Here, A _G determined the steepness of this transition, and the parameter b is a measure of the signal impact (Lipshtat et al., 2010).

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 [β ₁, β ₂], 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 d m i d t = A ⋅ m i 0 − m i ∑ i m i 0 − ∑ i m i − k slow + ( k fast − k slow ) ⋅ s n g Δ θ i , t (6) where A is a constant and k _slow and k _fast 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.

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.

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 PIP₃ 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.

To interpret the molecular mechanisms of such cellular behaviors, spatiotemporal plots of PIP₃/PIP₂ 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 PIP₃ 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 PIP₃ and PIP₂, respectively. After being activated by active Rac and RhoA, the catalytic effect of PI3K leads to the conversion of PIP₃ and PIP₂ and vice versa for PTEN.

Since membrane diffusion is small, a local increase in PIP₃ concentration could be confined locally, resulting in further recruitment of PI3K. A patch of PIP₃ is then formed during 30–75 s. As PIP₂ has the same diffusivity as PIP₃ and their total amounts are constant, the pattern of PIP₂ regulation mirrors that of PIP₃ (Figure 2C). The regions with high PIP₃ and PIP₂ 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 PIP₃, FilGAP starts to downregulate the activity of Rac and upregulate that of RhoA, which together results in a local decrease in PIP₃ (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.

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 PIP₃ 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 PIP₃ polarity.

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 PIP₃. Specifically, the spatiotemporal plot of PIP₃ (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.

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 (PIP₃)]. 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.

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.

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 PIP₃-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).

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 ∼ 4 μ m / min , such that the η profile manifests a slower increase. We then carried out another set of simulations within a larger obstacle region with a radius of r = 40 μm. Clearly, the cell with an intermediate level of FilGAP, e.g., (FilGAP) = 0.4 μM, exhibits the highest searching efficiency, demonstrating that the magnitude of the movement velocity becomes a limiting factor. Taken together, the results in this section imply that the complementary effects of the FLNa-FilGAP and PI-Rho GTPase pathways allow a cell to efficiently search its surrounding area when chasing the stimulus out of the obstacle region.

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.

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.

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