Two Complementary Signaling Pathways Depict Eukaryotic Chemotaxis: A Mechanochemical Coupling Model

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.


SUPPLEMENTAL FIGURES
. Responses of the modeled cell to a gradient stimulus. (A) Gradient profile as a function of the distance (r) between the cell and the stimulus source. Varying r from approximately 65 to 25 µm results in a 2% to 25% change (p). Here, C and ∆C denote the lowest concentration and increment, respectively. (B) Velocity profile in response to an initial 5% change in the stimulus gradient. Accordingly, the cell takes spontaneous (purple), adaptable (green), and directional motion (pink).

SUPPLEMENTAL VIDEOS
• Video S1 Correlated random migration of a "wildtype" cell.
• Video S2 Chemotaxis behaviors in a "wildtype" cell in response to a lateral gradient stimulus.
• Video S3 Chemotaxis behaviors in a "wildtype" cell in response to an opposite gradient stimulus.
• Video S4 Responses of a cell without FilGAP to a gradient stimulus.
• Video S5 Responses of a FilGAP overexpression cell to a gradient stimulus.
• Video S6 A "wildtype" cell behaves in a local region surrounded by equally spaced obstacles before and after chemoattractant stimuli.
• Video S7 A cell without FilGAP behaves in a local region surrounded by equally spaced obstacles.

Derivation of the chemoattractant field
The production, diffusion, and degradation of a chemoattractant molecule can be described by where C(r, t) is the concentration of the molecule at distance r from the point source at time t, D s is the diffusion coefficient of the chemoattractants, k d is the degradation rate, and k p is the production rate. The solution of this equation has been given by [22]: where λ is the space constant, λ = D s /k d , and τ n is the time constant where τ n = 1/(k d + n 2 πD/L).
Given a baseline concentration (C b ), the percentage change in the chemoattractant concentration (f ) across a small distance 2R (R is the radius of the cell) can be defined as: To achieve a desired f value, the distance between the initial position of the cell and the point source (x 0 ) can be calculated straightforwardly, where γ = C b L/(k d λ) · 1/ sinh(L/λ) and β = e 2R/λ .

Initial signaling processing
Local GPCR occupancy mirrors the local chemoattractant concentration. Hereafter, a balancedinactivation mechanism [11] is introduced to mimic the initial signaling process. This mechanism involves three interacting steps: (i) the local receptor occupancy level ([RL]) drives the production of membranebound species A and cytosolic species I at equal rates, k s ; (ii) the cytosolic species diffuses inside the cell and attaches itself to the membrane at a rate k I and becomes the membrane-anchored species I m ; and (iii) both species A and I m inactivate each other with a rate k i , and A and I m spontaneously degrade at rates δ A and δ I , respectively. The system equations can be written as with the boundary condition

Bidirectional molecular transport
Spatiotemporal regulation of PIP 3 (P 3 ) and PIP 2 (P 2 ) forms the core of the bidirectional molecular transport mechanism, which is described by the following equations: whereR = min{R/R max , 1} andρ = min{ρ/ρ max , 1}. In Eq. S9, the first term on the right-hand side accounts for PIP 3 diffusion, the second accounts for PIP 3 production due to membrane-bound PI3K (Hm) acting on PIP 2 , and the third accounts for PIP 3 diminution due to membrane-bound PTEN acting on PIP 3 . The parameters k PI3K cat (k PTEN cat ) and k PI3k M (k PTEN M ) are based on the steady state levels of PIs.R (ρ) is the normalized factor reflecting the effect of Rac (Rho) activity on PI3K (PTEN) activation, and R max (ρ max ) acts as a constant for Rac (RhoA) activity. IfR (ρ) is greater than R max (ρ max ), the activity of Rac (RhoA) is no longer a limiting factor for PI3K (PTEN) activation, andR (ρ) then equals unity. Eq. S10 describes the PIP 2 dynamics. Similarly, the first term on the right accounts for PIP 2 diffusion, the second accounts for PIP 2 production from PIP 3 via membrane-bound PTEN (T m ), and the third accounts for the reduction of PIP 2 into PIP 3 via membrane-bound PI3K.

Membrane mechanics
The elastic energy stored in the springs due to stretching or compression is given by where l i is the length of the ith spring that is iterated at every step, l 0 is the relaxed (zero-force) length, and K l is the effective stiffness constant of the spring. The elastic energy stored in the spring due to bending is given by where θ i is the angle of the ith spring that is iterated at every step, θ 0 is the relaxed (zero-force) angle, and K b is the spring constant for bending. In addition, an area constraint is also applied to ensure that the cell area is conserved within 1% during the simulation. This is implemented via an energetic penalty as follows: where s and s 0 are the instantaneous and equivalent areas of the cell, respectively, and K s is the penalty coefficient.
The total elastic energy (E m ) of the cell membrane is the sum of all three types of energy: The elastic force acting on the membrane particles is then calculated using the principle of the virtual work as follows: where P i is the position vector of the ith node. The viscous force is given by where γ is the viscosity coefficient of the cellular cytoskeleton, and v ij is the relative displacement rate of the neighboring nodes i and j.
The nodal protrusive force (F pro i ) is related to the local concentration of PIP 3 (P 3 ), and the nodal contractive force (F con i ) is correlated with that of PIP 2 (P 2 ). The general forms of F pro i and F con i are given by where ν pro and ν con are the protrusive and contractive force-concentration transfer factor, respectively, which are derived from the maximum active force generated by a cell, P 3s (P 2s ) is the saturation concentration of PIP 3 (PIP 2 ).

Chemotaxis index
The chemotaxis index (CI) of cells provided a measure of how well cell motion is directed toward an exogenous source. The instantaneous CI for cell i at time t is defined as

MODEL PARAMETERS
Model parameters are within the range of values previously measured and used in former models in the literature (see references in Table S1-S2). Parameters that are not well established were evaluated in the robustness analysis here as well as in our previous work [3].

Parameters for cell mechanics
The spring constants for cell bending (K b ) and stretching (K l ) were adopted from [28] upon the measurements of aspirated cells using a micropipette aspiration technique [29]. The penalty coefficient was set to be at least one order of magnitude larger than K l based on numerical rather than physical considerations. The force generated by a single actin filament can be estimated by the ratchet theory [21], which is a few tenths of the pN/µm. Considering that the density of F-actin at the leading edge is on the order of 10 6 , the protrusive force density at the cell front can be estimated on the order of hundreds of pN. The total amount of contractive force can be experimentally derived from the traction force measurement [16], which is on the order of tens of nN.

Parameters for initial signaling processing
D s and D m are diffusive coefficients for chemoattractants and membrane-bound molecules (i.e., activator and PIs), respectively. It is widely accepted that D s > 100 µm 2 /s, and D m is within the range of 0-1 µm 2 /s [22].
Since the guidance radius of chemoattractant L is on the order of hundreds µm and is generally not seen for distances > 500 µm in vivo [5], L was set to be = 100 µm here. Accordingly, the space constant, λ, is estimated to be 30 µm. Considering that 95% of the molecules are localized within a distance of 3λ from the source [22], the degradation rate of the chemoattractant can thus be estimated as k d = 0.1 s. Given that the concentration limit that does not prevent gradient detection is 0.01-100 nM, C b was set to 10 nM, a reasonable bottom bound of the concentration. Substituting the above parameters (with f = 2%, a minimum percentage change detectable by a chemotactic cell) into Eq. S4, we may estimate k p = 10 nM/s.

Parameters for cytoskeleton remodeling
Estimations of self-delay rates (δ R , δ ρ ), and basal activation rates (I R , I ρ ) for small Rho GTPases were based on the steady-state concentrations of their active forms, first achieved by [9]. and are involved in the bidirectional molecular transport module, which together determine the strength of Rho GTPase-PI feedback loops. A robustness analysis for two parameters was presented in our earlier work [3]. Briefly, only at the intermediate range of parameter combination (i.e., 6 molecules/s < k PI3K cat < 10 molecules/s) does the cell continually develop spontaneous motility.

Parameters in the mechanical-sensing module
Parameterizing the release rate of FilGAP upon stretching was performed by [10] in a minimally reconstituted system [2]. Accordingly, the experimentally derived, time-dependent fluorescence decay rate can be fitted with a two-exponential mixture model [10]. The control of the cytosolic concentration of FilGAP is also affected by the number of lamellipodial nodes (N ). We thus carried out a robustness analysis on N ( Figure S2), ensuring that the temporal evolution of [FilGAP] was not sensitive to N values.

NUMERICAL SIMULATION
The backbone of numerical simulations in our model is an iterative loop. Each time step occurs as follows: 1. At the beginning of the (n + 1)th time step, the system has already reached mechanical equilibrium.
The FDM and MC steps are implanted to derive the concentration fields for both membrane-bound molecules and effector molecules. 2. The protrusive and contractive forces on membrane nodes are updated based on the concentration fields of PIP 3 and PIP 2 , respectively.
3. The cell membrane is driven and deformed by membrane nodal forces.
4. The deformation of the actin network is calculated with the deformed cell membrane, providing the first boundary condition. 5. Once the geometry of the cell is updated, the cytosolic concentration of FilGAP is updated with the bandpass mechanism, and a new time step then begins.

One dimensional finite differences method
As the molecular diffusion on the membrane (circle) of 2D cell is one-dimensional, we use 1D finite difference to solve the diffusion equation on the cell membrane [28] according to Fick's laws. The 1D diffusion flux J of the molecule from node i to i + 1 on the membrane, where the subscript i indicates the ith node on the cell membrane, i is the edge length between node i and i + 1, L i = ( i + i−1 )/2 , D is the diffusivity on the membrane. Therefore, in one simulation time step ∆t, the change of molecular concentration caused by diffusion at the ith node is

MODEL ASSUMPTIONS
Several simplifying assumptions are made in formulating the model, and the major assumptions are specified below.
1. In our model, we introduce a simple balanced-inactivation (BI) mechanism [11] for sensing chemoattractant stimuli. This parsimonious mechanism affords a natural role for achieving a switchlike distribution pattern of G α and G βγ . Other candidate directional sensing mechanisms involve the local modulation of chemoattractants [13] or a phase separation mechanism [4], although they are far more complex in mathematics.
2. The model assumes that the generation of protrusive and contractive forces is in accordance with PIP 3 and PIP 2 , respectively. The relationship between PIP 3 regulation and protrusive force generation is abundant, as PIP 3 may provide membrane binding sites for actin binding proteins, such as Arp2/3 [1]. Moreover, fluorescent imaging data suggest that PTEN contains an N-terminal PIP 2 binding motif, and its deletion completely redistributes the enzyme into the cytosol [7,27]. Since PTEN is an upstream signal transducer for myosin II localization, PIP 2 may regulate contractive force via myosin-II with the help of PTEN [23].
3. We assume that in the bidirectional molecular transport mechanism, the feedback loops between Rho GTPase and PIs are mediated via local recruitment and activation of PI3K and PTEN. However, experimental evidence indicates that PIP 3 also enhances Rac activity by promoting the localized activation of RacGEFs [1]. In earlier modeling works [14], such a feedback effect is included in PIP 3 concentration-dependent rate constants. 4. The translocation behavior of FilGAP between the cytosol and the membrane is not explicitly considered in our model. FilGAP contains an N-terminal PH (pleckstrin homology) domain, a GAP activity region, a spacer, and a C-terminal CC (coiled-coil) domain that mediates protein dimerization [19,18]. It has been well established that the FilGAP C-terminal FLNa-binding site binds to 23 of 24 Ig-like repeats comprising the FLNa dimeric subunit structure (IgFLNa23) [18]. On the other hand,