Altering the threshold of an excitable signal transduction network changes cell migratory modes.

The diverse migratory modes displayed by different cell types are generally believed to be idiosyncratic. Here we show that the migratory behaviour of Dictyostelium was switched from amoeboid to keratocyte-like and oscillatory modes by synthetically decreasing phosphatidylinositol-4,5-bisphosphate levels or increasing Ras/Rap-related activities. The perturbations at these key nodes of an excitable signal transduction network initiated a causal chain of events: the threshold for network activation was lowered, the speed and range of propagating waves of signal transduction activity increased, actin-driven cellular protrusions expanded and, consequently, the cell migratory mode transitions ensued. Conversely, innately keratocyte-like and oscillatory cells were promptly converted to amoeboid by inhibition of Ras effectors with restoration of directed migration. We use computational analysis to explain how thresholds control cell migration and discuss the architecture of the signal transduction network that gives rise to excitability.

INTRODUCTION

Diverse cells display different migratory modes. For example, amoeba move by repeatedly extending and retracting pseudopods, keratocytes glide with a single broad anterior protrusion, fibroblasts slowly project filopodia and lamellipodia with strong attachment, and some cells display oscillatory behavior[1-6]. Transitions between these migratory behaviors can be important, for example as cancer cells become metastatic[7,8]. While these migratory modes likely depend on a constellation of expressed genes, the connections among different genetic profiles, various types of protrusions, and cell migratory modes are not understood.

In amoeboid cells, random migration requires activation of the same signal transduction networks that are triggered by extracellular cues[9-11]. These networks contain multiple compensatory pathways involving Ras super family GTPases, phosphoinositides, and kinases and phosphatases[12-15]. Coordinated increases or decreases in these activities start spontaneously and propagate in waves along the cell cortex that are spatiotemporally correlated with actin waves, which are thought to drive pseudopodia formation at the edge of the cell[11,16-26]. Thus an understanding of the relationship among waves of signal transduction activities, cytoskeletal events, and cellular protrusions may reveal the essential differences between distinct migratory modes.

We have proposed that the signal transduction components involved in cell migration are linked in a signal transduction excitable network, or STEN. In addition to wave propagation, evidence for excitability includes all-or-none responses to supra-threshold stimuli and the existence of a refractory period[11,17]. In the classical example of the action potential excitability derives from self-amplifying positive and delayed negative feedback loops comprised of voltage-gated ion channels, whereas in cell migration these loops most likely involve biochemical interactions. Computational models based on excitability have led to realistic simulations of random and directed cell migration[11,16,27-32]. These models suggest that small changes in thresholds can alter wave properties and lead to distinct patterns of behavior, but this has not been experimentally tested.

Here, we hypothesize that the spectrum of cell migratory modes observed in various cells arises from the different set points, or thresholds, of STEN. Since components in the STEN undergo highly coordinated transient changes during network activation, “clamping” one component near the level it achieves during activation (or inactivation) might alter the excitability of the entire network, offering an excellent opportunity to test the idea.

RESULTS

Acute perturbations of STEN components control migratory profiles

We used a chemically inducible dimerization system in Dictyostelium to clamp PIP2 at low levels or the activities of Ras/Rap GTPases at high levels, as would be expected to transiently occur during STEN activation[33-38]. To decrease PIP2, we co-expressed myristoylated FKBP and the yeast PIP2-specific phosphatase Inp54p[39,40] fused to FRB. Upon addition of rapamycin, mCherry-FRB-Inp54p was rapidly and irreversibly recruited to the plasma membrane (Fig. 1a). Several PIP2 biosensors with different apparent affinities, including PHPLCδ and derivatives of nodulin[41], suggested that PIP2 levels decreased rapidly, reaching a new lower steady-state level within 15 minutes (Fig. 1b,Supplementary Fig. 1a). After rapamycin addition, the speed of the cell population increased to 7.3±2.7 from 4.2±2.0 μm/min, but not when inactive Inp54pD281A was recruited (Supplementary Fig. 1b,c, Video 1). Similar increases in motility occurred when we recruited RasCQ62L(∆CAAX), Rap1G12V(∆CAAX), or the GEF domain of Rap1 GEF, GbpD[42] (Supplementary Fig. 1d–f).

Figure 1

Acute clamping of PIP2 at lowered levels triggers cell migratory mode transitions

(a) Left, confocal images of myr-GFP-FKBP-FKBP (green) and mCherry-FRB-Inp54p (red) before and after rapamycin treatment. Right, temporal profile of normalized intensity of cytosolic mCherry (mean±s.d., n=25 cells). (b) Left, confocal images of PIP2 biosensors PHPLCδ-YFP and GFP-NodulinQ534P before and after rapamycin treatment. Right, normalized ratios of membrane to cytosol intensity of PHPLCδ-YFP (black), GFP-NodulinQ534P (red), and GFP-Nodulin (blue) following rapamycin treatment at time 0 (mean±s.e.m, n=21 cells for each biosensor). Initial data are shown in Supplementary Fig. 1a. (“Nodulin” refers to the Nlj6-like nodulin domain of the Arabidopsis Sec14-nodulin protein AtSfh1[41]. NodulinQ534P is a mutant with defects in self-oligomerization and thus diminished affinity for PIP2 compared with Nodulin.) (c) Time lapse phase-contrast images showing the transition of an amoeboid cell to a fan-shaped (top) or an oscillatory cell (bottom). Scale bars in all images represent 10 μm. (d) Centroid tracks showing random movement of amoeboid (left), fan-shaped (middle), and oscillatory (right) cells. Each track lasts 10 min and was reset to the same origin. Insets show color-coded (1 min intervals) outlines of a cell. Velocity (μm/min) and directedness are 4.2±1.6 and 0.39 for amoeboid, 11.7±1.8 and 0.91 for fan-shaped, and 7.5±2.0 and 0.40 for oscillatory cells (mean +/− s.d, n=50 cells for each). (e) Temporal profiles of normalized areas of 10 control (left) and experimental cells (right). Dashed lines indicate the addition of rapamycin at time 0. Based on the corresponding movies, amoeboid, fan-shaped, and oscillatory modes are colored as blue, green, and red respectively. (f) Fractions of each migratory mode in 10-min time windows in a cell population of ~100 cells before and after Inp54p recruitment. 0.5 μM (top) or 5 μM (bottom) rapamycin were added at time 0.

Careful examination revealed that initially amoeboid cells switched to fan-like or oscillatory behavior (Fig. 1c, Supplementary Fig. 2, 3, Video 2). “Fans” resemble gliding keratocytes, with broad protrusions along a wide anterior zone and a flat contractile back. “Oscillators” spread globally, contracted strongly, and migrated briefly before repeating the cycle. The appearance of these distinct populations could also be detected with computer algorithms which defined fans, or oscillators, as cells that either migrated perpendicular to their long axis or displayed a large variance in area, respectively (Supplementary Fig. 4). The speed and directedness were much higher in the fan-shaped than in the amoeboid cells (Fig. 1d). The oscillatory cells were also somewhat faster than the amoeboid since they often glided into the spreading phase. The heterogeneity within the population is illustrated in Figure 1e. Some cells transitioned persistently to fans or oscillators while others switched between modes. Retrospective examination of wild-type cells revealed that about 1% of cells displayed fan-like or oscillatory behaviors, but this increased to ~10% when expressing cytosolic Inp54p. Following recruitment, the probability of the fan-shaped and oscillatory modes grew to about 20% and 50% respectively during 40 min (Fig. 1f). At a lower dose of rapamycin, fans were more prominent than oscillators at early time points but the same final distribution was reached (Fig. 1f), suggesting that the fan-shaped cells represent an intermediate phenotype.

We used LimE∆coil-YFP (LimE) to monitor actin polymerization. Three-dimensional kymographs, or T-stacks, showed that newly polymerized F-actin was initially confined to patches but was localized in broad zones after the addition of rapamycin (Fig. 2a, Supplementary Video 3). Time-lapse videos revealed the spatiotemporal F-actin dynamics in different modes of migration (Fig. 2b,d,f, Supplementary Video 4). In kymographs of the amoeboid cells narrow regions of LimE intensity appeared stochastically with low coordination (Fig. 2c). In fans, broad, sometimes V-shaped, regions centered in the direction of migration appeared periodically, rarely extending to the back of the cell (Fig. 2e). In oscillators, LimE intensity was highly coordinated with cell area, increasing and decreasing during the spreading and contracting phases respectively (Fig. 2g,h). LimE usually occupies 10–40% of the cell periphery in amoeboid cells, 40–70% in fans, and a broad range including 0 and 100% in oscillators (Figure 2i). None of these intensity patterns was observed when we monitored cell periphery using recruited mCherry-FRB-Inp54p (Supplementary Fig. 5a–d). Similar LimE patterns were observed in cells after recruiting RasCQ62L(∆CAAX), Rap1G12V(∆CAAX), and GbpD-GEF (Supplementary Video 5). These observations suggest that the lateral dimension of actin-filled projections plays a key role in cell migratory behavior.

Figure 2

Actin profiles show distinct spatiotemporal patterns in amoeboid, fan-shaped, and oscillatory cells

(a) T-stacks of LimE∆coil-YFP (LimE) in initially amoeboid cells before and after recruiting inactive (left) or active (right) Inp54p. Red dashed line indicates the addition of rapamycin. (b, d, f) Time-lapse confocal images and color-coded overlays of the YFP fluorescence intensities at 2 sec intervals of an amoeboid (b), fan-shaped (d), and oscillatory cell (f) expressing LimE. Scale bars represent 10 μm. (c, e, g) Kymographs of cortical LimE intensity in the same amoeboid (c), fan-shaped (e), and oscillatory cells (g). Durations of the kymographs are 5 min (c, e) and 10 min (g). (h) Temporal profiles of cell area (black) and normalized mean cortical LimE intensity (green) of the oscillatory cell above. Normalization for changes in cell shape is based on the kymograph of stably recruited mCherry-FRB-Inp54p intensity in the same cell as shown in Supplementary Fig. 5d. (i) Relative frequency of different LimE sizes (percentage of cell perimeter) appearing in each mode of migration (>2000 time points combined from 10 kymographs of each mode automatically thresholded and quantified in MATLAB).

STEN activities are required for the transition in migratory modes

We examined whether activation of the STEN coincided with the large protrusions in fans and oscillators. As representatives of STEN activity, we used RBD-YFP, PHcrac-YFP, and PTEN-GFP. RBD, which monitors multiple Ras activities, and PHcrac, which reports PIP3 levels, are found preferentially at the front of migrating cells, at the tips of newly projected protrusions, in the active zones of propagating waves, and associate transiently with the plasma membrane during global stimulation. PTEN displays complementary localization and kinetics. In oscillatory cells, the “front” biosensors were recruited to the membrane globally during the spreading phase and returned to the cytosol when the cell contracted whereas the “back” marker, PTEN, moved oppositely (Fig. 3a–c, Supplementary Fig. 5e–g, Video 6). In migrating fans RBD and PHcrac were recruited to the front while PTEN located at the rear (Fig. 3a–c).

Figure 3

A signal transduction network mediates the transitions of cell migratory modes

(a, b, c) Dynamics of signal transduction biosensors RBD-YFP (a), PHcrac-YFP (b), and PTEN-GFP (c). In each panel, top shows time-lapse confocal images of cells transitioning from amoeboid to oscillatory mode after rapamycin addition at time 0; bottom left shows kymograph of cortical intensity of each biosensor in an oscillatory cell; bottom middle shows temporal profile of cell area (black) and normalized mean cortical signaling intensity (green) in the corresponding oscillatory cell (normalization for changes in cell shape based on the kymograph of stably recruited mCherry-FRB-Inp54p intensity in the same cell as shown in Supplementary Fig. 5e–g.); bottom right shows confocal image of a fan-shaped cell. Scale bars in all images represent 10 μm. (d) Percentage (mean±s.d., n=3experiments, >500 cells in each experiment, ***P<0.0001, one-way ANOVA with post-hoc Tukey HSD test) of fan-shaped plus oscillatory cells 10 min before (blue bars) and 10–20 min (red bars) after the addition of rapamycin at time 0. (e) Distributions of different sizes of RBD-GFP (top) and PHcrac-YFP (bottom) activity in cells of wt, pkbA-/pkbR1-, and pkbA-/pkbR1- overexpressing PKBA. Cells were treated with 5 μM LatrunculinA for 30 min. (f) Distributions of different sizes of PHcrac-YFP activity in wt cells, as well as wt cells expressing myr-FKBP-FKBP plus PKBA-mCherry-FRB, treated with DMSO or rapamycin in addition to 5 μM LatrunculinA for 30 min. More than 3000 cells from 3 independent experiments were analyzed in each condition in e and f.

Next we sought to determine the extent to which the transition among migratory modes depended on changes in gene expression or STEN activities. Addition of 1 mg/ml cycloheximide did not block the induced migratory mode transitions (Supplementary Video 7). Previous studies have shown that heterotrimeric G-proteins activate Ras and Rap small GTPases, leading to the activation of PI3K and TORC2, which converge on PKBs activities to organize actin polymerization. Here loss of PKBA and PKBR1 (pkbA-/pkbR1-) blocked the induced mode transitions (Figure 3d). Inhibition of neither PI3K nor TORC2 activities was sufficient to block the transitions, but combined inhibition with LY294002 plus PP242, or by LY294002 treatment of cells lacking TORC2 subunit PiaA,nearly prevented the formation of fans and oscillators. Interference with small GTPases by expressing dominant negative Rap1 (DNRap1), or double deletion of RasC and RasG (rasC-/rasG-) diminished the transitions, and combined perturbations blocked the responses. G-protein function was not required since Gγ deletion did not prevent the induced transitions. These suggest that the same STEN activities involved in chemoattractant stimulation are required for the enlarged protrusions underlying fans and oscillators.

Given that PKBs link the activity of the STEN to the cytoskeletal network, it is not surprising that they are required for the mode transitions, but these kinases have also been reported to negatively regulate RasC[36]. To access the effects of PKBs without interference of the cytoskeleton, we treated cells with LatrunculinA to inhibit actin polymerization and monitored the activity of Ras and PI3K directly, using RBD and PHcrac. Cells lacking PKBs (pkbA-/pkbR1-) showed enhanced activities of both biosensors, which were suppressed by expressing PKBA (Fig. 3e). Similarly, elevated RBD activities were observed when LY294002 and PP242 were used to inhibit upstream regulators of PKBs (Supplementary Fig. 6a–f). In addition, recruitment of PKBA in wild type cells with the dimerization system suppressed Ras-mediated PI3K activities (Fig. 3f, Supplementary Fig. 6g,h). Thus the negative regulation by PKBs extends beyond RasC to other Ras small GTPases and PI3K. Interestingly, during chemotactic stimulation, activation of PKBs is delayed relative to the increase in Ras and PI3K activities[36,37,43], suggesting they could act as a delayed negative feedback regulator of the excitable network.

Migratory profiles are dictated by properties of propagating STEN waves

To examine the changes in wave behavior that accompanied the transitions to different migratory modes we monitored STEN and cytoskeletal activities, using biosensors PHcrac and LimE, respectively, on the basal surface of cells by TIRF microscopy (Fig. 4a,b and Supplementary Video 8, 9). As previously described in amoeboid cells, LimE and PHcrac wave activities started at a point and propagated a few microns before extinguishing. In oscillators, these activities began in a larger region, expanded in a wave that momentarily covered the basal surface of the cell, then extinguished from the center outward just before the cell collapsed. In fans, the expansion was similar except that it was not symmetric and did not extend to the rear of the cell and therefore pushed out a wide band on the front. This was often followed immediately by another asymmetric wave which kept the cell gliding in the same direction. Note that LimE and PHcrac maintain the same relative distribution pattern (Fig. 4c, Supplementary Video 10) as previously published[44]. Consistent with the increased width and range, the waves in fans and oscillators move faster (Fig. 4d).

Figure 4

Properties of the STEN waves are altered independent of cytoskeleton

(a, b) Merged time-lapse TIRFM images of each mode of migration. PHcrac-YFP (green) and mCherry-FRB-Inp54p (red) are shown in a, and LimE-YFP (green) and mCherry-FRB-Inp54p (red) in b. The mCherry-FRB-Inp54p signals serve as membrane markers following recruitment. The white arrows point to new waves initiating at the front of fan-shaped cells. (c) Time-lapse TIRFM images of PHcrac-YFP (green) and LimE-RFP (red) in the same oscillatory cell. (d) Scatter plot of wave speed in each mode of migration (n=31 cells for each mode and the horizontal lines indicate median, P=0.00019 comparing amoeboid and fan-shaped, P<0.0001 comparing amoeboid and oscillator, Wilcoxon-Mann-Whitney rank sum test). (e) PHcrac-YFP activities in LatrunculinA-treated cells. Left, confocal images of PHcrac-YFP in the same cells before and after Inp54p recruitment. Right, Distributions of different sizes of PHcrac-YFP activity in cells expressing mCherry-FRB-Inp54p(D281A) or mCherry-FRB-Inp54p (in addition to myr-FKBP-FKBP) treated with DMSO or rapamycin for 30 min. More than 3000 cells from 3 independent experiments were analyzed in each condition. (f, g) Kymographs of cortical PHcrac-YFP (f) or PTEN-GFP (g) intensity in cells treated with 5 μM LatrunculinA. The red dashed lines indicate time 0 at which rapamycin was added. Solid vertical lines on the left of each kymograph represent 20 μm. Scale bars in all images represent 10 μm.

Next we assessed whether the increased STEN wave activities require an intact cytoskeleton. In cells immobilized with LatrunculinA, PHcrac still exhibited occasional flashes at the cell membrane. After Inp54p recruitment these patches became more frequent and larger (Fig. 4e, Supplementary Video 11). Kymographs of individual cells show that PHcrac as well as PTEN began to exhibit periodic waves along the periphery (Fig. 4f,g). Thus the switch in the dynamic behavior of the STEN is independent of the actin cytoskeleton or the normally accompanying morphological changes.

The altered wave properties can be traced to changes in the threshold for STEN activation

The expanded STEN waves could indicate that the threshold to activate the excitable network is lowered. To assess the sensitivity of the STEN, we monitored the translocation of PHcrac in response to chemoattractant in the same cells before and after Inp54p recruitment. Figure 5a shows that the same set of cells, which did not respond to low amounts of cAMP, did respond to the same dose after Inp54p recruitment, while the responses to saturating levels of stimulus were nearly identical (Fig. 5b). Further, the concentration of cAMP eliciting a half-maximal response shifted from about 1 nM to 200 pM after lowering PIP2 levels (Fig. 5c), which was not observed when recruiting inactive Inp54pD281A (Fig. 5d). Recruiting RasCQ62L (∆CAAX) also slightly left-shifted the dose-response curve to cAMP (Supplementary Fig. 7). These results suggest that the perturbations that bring about expanded waves and migratory mode transitions also increase the excitability of the STEN.

Figure 5

The threshold for STEN activation is lowered

(a, b) Responses of PHcrac-YFP to global cAMP stimulations in the same cells before (blue box) and after (red box) Inp54p recruitment. Time-lapse confocal images of PHcrac (green) and mCherry-FRB-Inp54p (red) are shown on the top. cAMP was added at time 0 (316 pM cAMP in a and 100 nM cAMP in b). Scale bars represent 10 μm. Temporal profiles of normalized mean cytosolic PHcrac intensity are shown at the bottom, with blue line indicating before and red after rapamycin treatment (mean±s.e.m., n=26 cells for a and n=18 cells for b). (c, d) Normalized PHcrac responses (drop of cytosolic intensity) to different doses of cAMP before (blue) and after (red) rapamycin. c, recruiting Inp54p. d, recruiting inactive Inp54p. Response in each individual condition was normalized to that upon 100 nM cAMP stimulation (mean±s.e.m., from the lowest to the highest cAMP concentrations, n=24, 21, 26, 21, 25, 18 cells in c, and n=21, 21, 23, 23, 21, 23 cells in d, *P=0.046, **P=0.002, ***P<0.0001, two-tailed t-test).

Simulations altering the threshold of excitable network capture the entire spectrum of cell behaviors

To see the effect of perturbations on an excitable system, we simulated an activator-inhibitor system coupled to a polarization module using stochastic, partial-differential equations on a 2-D surface (Methods and Supplementary Fig. 8a,b). The model was parameterized to recreate the pattern of activity seen in amoeboid cells. Simulations recreated waves of activity seen by TIRF which propagate outward and extinguish relatively quickly (Fig. 6a, Supplementary Video 12). When we only lowered the threshold for triggering excitable behavior, we observed more robust waves that traveled faster and farther. The frequency of firings also increased leading to more waves being present at any time (Fig. 6b). After additional reduction of the threshold, nearly the whole field underwent synchronized periods of activity and inactivity (Fig. 6c).

Figure 6

Simulations of the excitable network with varying thresholds successfully capture various cell migratory modes

(a, b, c) Two-dimensional simulations showing waves of activity propagating outward and extinguishing, with gradually decreasing threshold from a to c. Green and red represent the activator (F) and inhibitor (R) activities, respectively. Yellow denotes the region where the two overlap. (d, e, f) Kymographs (top) and simulated cell morphologies (bottom) give rise to amoeboid, fan-shaped and oscillatory behaviors, respectively. The cell shapes are color-coded from blue to red, showing the same cell at subsequent times; the colored crosses denote a fixed origin. (g) Kymograph (top) of the response to a gradual threshold decrease over a larger time scale. The threshold was held constant until t=100 (dotted line), after which it was lowered exponentially. The lower plots show the maximum (blue) and mean (red) levels of activity. (h) Fraction of cells in one of the three migratory modes as a function of time, with slowly (left) and rapidly (right) decreasing threshold obtained from simulations of 45–55 cells. (i) Normalized response curve to different levels of external stimuli (cAMP) for higher (blue) and lower (red) thresholds. Normalization is with respect to the 100 nM cAMP response.

We used a 1-D model of the excitable system coupled to a viscoelastic cell model in the level set framework to determine the effect of these changes on cell migratory modes[27]. In these simulations, initially randomly extended projections gradually became localized toward the front as the cells polarized (Fig. 6d, Supplementary Video 13). As we lowered the threshold, the waves that guided projections moved faster and spread further to cover roughly half of the cell and were triggered periodically while the back remained quiescent (Fig. 6e, Supplementary Video 13), resulting in a fan-shaped cell. When the threshold was further reduced, multiple initiations appeared and the waves rapidly encompassed the whole cell causing it to spread. This symmetric response was closely followed by a retraction as the waves extinguished, before starting the next cycle (Fig. 6f, Supplementary Video 13).

The migratory mode transitions were observed in a single simulation by gradually lowering the threshold (Fig. 6g). Plotting the mean and maximum intensity of either the activator or the inhibitor of the excitable network over the perimeter as a function of time revealed behavior that coincided with different modes: the mean activity of amoeboid cells is low, but the maximum is high and fluctuates; in fans, the maximum is persistently high, but the mean increases as the size of the active zone expands; and in oscillators the maximum and mean levels of activity always coincide (Fig. 6g, Supplementary Fig. 8c). During the transition from amoeboid to fan-shaped, the front of the cell expanded outward while the back of the cell retracted (Supplementary Video 14). As the cell transitioned from fan-shaped to oscillatory, the activity spread instantly around the entire perimeter as the cell spread and then disappeared causing the contraction (Supplementary Video 14).

To recreate cell-to-cell heterogeneity, we continuously lowered the threshold as above but randomized its initial and final values in different cells. At each time point, different migratory modes were observed among the cell population. As was observed experimentally, the predominant mode shifted in time from amoeboid to fan-shaped to oscillatory (Fig. 6h). Curiously, we occasionally saw “reverse” shifts such as ones from oscillatory to fan-shaped cells in the simulation, similar to those observed in cells (compare Fig. 1e and Supplementary Fig. 8c). We also altered the rate at which we lowered the threshold and observed that the transition away from amoeboid cells was faster for higher rates of reduction (Fig. 6h) as seen experimentally (Fig. 1f). Finally, we simulated the dose-response behavior of a heterogeneous population of cells before and after lowering the threshold. The curves showed a left shift (Fig. 6i) matching the experimental observations (Fig. 5c).

Conversion of innately fan-shaped or oscillatory cells to amoeboid mode recovers directed migration

We examined pikI- cells lacking a major PIP-5-kinase[45] (PI5K), which were reported to have decreased PIP2 levels and to fail to perform chemotaxis. We observed many fan-shaped or oscillatory cells in the population (Fig. 7a, Supplementary Video 15). Further, cytoskeletal and signal transduction biosensors exhibited identical patterns in these oscillatory cells as in those induced by acute perturbations (Fig. 7b). As this suggested the presence of excess STEN activity, we reasoned that inhibition of the coupling of these events to the cytoskeleton could convert the cells to the amoeboid mode. Indeed, addition of PI3K inhibitor caused the cells to round up initially. Remarkably, within about 30 min more than 90% cells were moving in the amoeboid mode (Fig. 7c, Supplementary Video 15).

Figure 7

Innately fan-shaped or oscillatory cells can be converted to amoeboid mode

(a, c, d, e) Phase-contrast images of pikI- (a, c) and amiB- cells (d, e) treated with (c, e) or without (a, d) LY294002 (LY). Color-coded outlines (1 min apart) of several cells were imposed on top of the phase images, with magenta outlining the current cells. Scale bars represent 25μm. (b) T-stacks of merged confocal images of oscillatory pikI- cells expressing RBD-YFP and LimE-RFP (left) or PTEN-GFP and LimE-RFP (right). Blue arrows indicate the spreading phases during the oscillation. (f) Percentage of fan-shaped cells in amiB- before and after 30 μM LY treatment, and after 20 μM PP242 treatment (mean±s.d., n=3 experiments, >100 cells in each experiment, *P=0.0043, **P=0.0002, one-way ANOVA with post-hoc Tukey HSD test). (g, h) Chemotaxis assays before (red) and after (blue) LY treatment. Number of cells migrating through pores into folic acid filled chambers and controls are shown (mean±s.e.m., n=3 experiments, >7500 cells in each experiment, *P=0.0036, 0.0033, 0.0013 in g, *P=0.0019 and ***P<0.0001 in h, two-tailed t-test). g, pikI-; h, amiB-.

A previously reported mutant, amiB-[46], consists of almost entirely fan-shaped cells when starved (Fig. 7d). The migration profile of these cells is similar with those we produced with our acute perturbations. As above we treated amiB- fan-shaped cells with the PI3K inhibitor. The amiB- cells initially rounded up, but within 30 min, greater than 70% of amiB- cells recovered and continued to move in the amoeboid mode for a couple of hours (Fig. 7e,f, Supplementary Video 16). TORC2 inhibitor PP242 had a similar effect (Fig. 7f).

Consistent with previous results we observed that both pikI- and amiB- cells failed to carry out chemotaxis over the range of folic acid gradients. We speculated that the excess STEN activities in these cells might impair directed migration and that conversion of the cells to amoeboid modes of migration might restore their ability to move up chemotactic gradients. Indeed, after treatment with the PI3K inhibitor, both pikI- and amiB- cells displayed greatly improved chemotaxis to every folic acid concentration tested (Fig. 7g,h).

DISCUSSION

Our results show that various cell migratory modes arise from different thresholds of an excitable signal transduction network. In simulations small decreases in the threshold of an excitable network were sufficient to expand the range of waves to cover increasing fractions of the cell perimeter, and this led to dramatic changes in protrusions and cell migratory modes. Our experimental results demonstrate a causal chain of events where changing the threshold alters the range of propagating STEN waves and in turn cytoskeletal activities and the nature of cellular protrusions, which promotes transitions among migratory modes.

We propose that local regions of the cell cortex transition between inactive and active states which we designate as B and F, respectively, since “back” and “front” markers such as PTEN and PI3K, tend to associate with these regions (Fig. 8a). We propose a refractory state, R, follows the F-state. The B- and F-states are mutually inhibitory, creating a positive feedback loop. The F- and R-states are related through a delayed negative feedback loop. In resting cells, most of the cortex is in the B state. Once initiated, waves propagate unidirectionally because diffusion of F-state components triggers activation in adjoining B but not R regions. Lowering the threshold for activation promotes the F-state, thereby increasing the speed and range of the propagating waves, leading to fans and oscillators (Fig. 8b,c). Bear in mind that the cartoon in Figure 8b is a one-dimensional representation of the more-or-less circular waves propagating on the cell cortex as seen in the TIRF videos and in the two-dimensional simulations. The coupling of the F-state waves to cytoskeleton we have proposed explains the expanding cup-like structures observed by lattice light sheet microscopy of randomly migrating Dictyostelium cells[47,48].

Figure 8

Propagating waves of STEN activity control cell migratory modes

(a) Snapshot of the cortical membrane showing regions where the STEN is in B- (grey), F- (green), and R- (red) states. Positive feedback, depicted as mutual inhibition between F- and B states, and negative feedback from R- to F-states, together with diffusion of network components, advances the wave unidirectionally (black arrow). (b) During amoeboid migration (left), network activities initiate locally and spread in waves that dissipate over a short distance. Since the F-state organizes cytoskeletal activities, this pattern controls the formation of pseudopodia. When the threshold of the network is lowered (middle), the waves initiate at larger zones and spread faster and further around the cell before dissipating. When the threshold is further lowered, initiated waves propagate even further to swiftly cover the nearly entire membrane. Color scheme is same as in a. (c) Schematic representation of STEN and cytoskeletal activities in amoeboid, fan-shaped, and oscillatory cells. (d) Proposed molecular architecture of STEN. PIP2 levels antagonize while GbpD(GEF), Rap1, and RasC promote the positive feedback loop (green) of the excitable network. PKBs serve in a delayed negative feedback loop (red) through a set of their substrates including Sca1 and PI5K, and promote F-actin through other sets of substrates.

Our results provide important guidelines for the molecular architecture that confers excitability to the network (Fig. 8d). Our studies implicate RasC, Rap1, GbpD, PIP2, and PKBR1/PKBA in the feedback loops that regulate the transitions between states of the cell cortex. They suggest that activation of Ras/Rap-related activities, or lowering PIP2, promotes positive feedback. The PKBs appear to have dual roles: On one hand, they are strong candidates for the delayed negative feedback that resolves the activation since STEN activities were elevated in cells lacking both PKBs and were suppressed upon recruitment of PKBA (Fig. 3e,f and Supplementary Fig. 6). These are consistent with previous reports suggesting that PKBs indirectly inhibit Ras GEF Aimless[36] and directly activate PI5K[43] which would increase PIP2 levels. On the other hand, likely through a different set of substrates, the PKBs are required to link the elevated STEN activities to the cytoskeleton since their inhibition blocked expanded protrusions (Fig. 3d).

We found that decreasing PIP2 levels effectively lowered the threshold of the STEN even in the absence of F-actin (Fig. 4e–g). PIP2 has been reported to play diverse roles in cytoskeletal organization and regulation[49-54]. We speculate that lowering PIP2 levels might function by promoting the activation of multiple Ras and Rap proteins. Activated Ras signals were detected at the entire periphery of induced spreading cells (Fig. 3a), and acutely increasing Ras or Rap activities had a similar effect as lowering PIP2 levels (Supplementary Video 5). Furthermore, combined inhibition of several Ras and Rap GTPases blocked the migratory mode transitions induced by decreasing PIP2 (Fig. 3d).

We propose that directional sensing could be achieved if an extracellular cue lowered the threshold at the front and raised it at the back. Thus, there would be more and larger protrusions at the front compared with the back, which could bring about guidance with high sensitivity. Of course, directed cell migration would be impaired at either end of the dynamic range. In fact, manipulations of amiB- and pikI- cells showed that amoeboid outperformed oscillatory and keratocyte-like fan-shaped cells in chemotaxis (Fig. 7g,h). Inhibition of PI3K is effective in restoring the chemotactic ability of pikI- cells, implying that the system can still function with lowered PIP2 and PIP3 levels. Thus what is important for proper chemotaxis is an optimal level of STEN activity and cytoskeletal coupling, rather than simply the activity or amount of a certain molecule.

The STEN wave scheme of cell migration has important implications. First, it can explain why diverse cell types display predominantly different migratory modes, since the resting set point of a cell is expected to be determined intrinsically by the expression profiles of network components. Second, it would allow cells to rapidly adjust behavior to various environments, using nearly the same components, which might underlie the plasticity of migrating cells as during cancer metastasis. Third, it enables a cell to integrate multiple chemical and mechanical cues, if they all fed into controlling the threshold of the same STEN. Finally, the natural heterogeneity in the set points of the STEN among individual cells within a population can lead to different behaviors and a range of responses to the same cue, facilitating sorting and patterning.

METHODS

Cells and plasmids

Wild-type Dictyostelium discoideum cells of the AX2 strain, obtained from the R. Kay laboratory (MRC Laboratory of Molecular Biology, UK), were used in this study. The piaA- and pkbA-/pkbR1- cells were described previously[37, 43]. rasC-/rasG- and amiB- were obtained from DictyBase stock center. Gγ- cells were obtained from the M. Ueda laboratory (Osaka University, Japan), and pikI- from the R. Kay laboratory. pikI- cells were grown on Klebsiella aerogenes lawn on SM plate, while wild-type and other gene deletion cell lines were cultured axenically in HL5 at 22°. In addition, pkbA-/pkbR1- cells were supplemented with heat-killed Klebsiella aerogenes during culture. Growth-stage cells were used for most experiments (except for cAMP stimulation experiments). For TIRF microscopy, starvation in suspension for 1-2 hours was carried out in DB at 2 × 107 cells/ml, in order to reduce photosensitivity, which allowed for prolonged and frequent imaging. Cells were used within 2 months of thawing from frozen stocks.

To make membrane-anchored FKBP constructs, sequences encoding the first 150 amino acids at the N-terminus of PKBR1 (denoted as “myr”), as well as tandem FKBPs from ARGENT™ Regulated Heterodimerization Kit (Version 2.0), were linked with a linker sequence of AGTGCTGGTGGT and cloned into pCV5 and also pDM358. Sequences encoding mCherry, FRB, and 1-331 amino acids of Inp54p (or inactivate form Inp54pD281A) were linked and cloned into pCV5 as well as pB18. The linker sequence between mCherry and FRB is GGAGCAAGTGGA, and two repeats of AGTGCTGGTGGT were used between FRB and Inp54p. Similarly, mCherry-FRB-RasCQ62L(4-555bp) was constructed with 5 repeats of AGTGCTGGTGGT between FRB and RasC, as with mCherry-FRB-Rap1G12V(1-546bp) and mCherry-FRB-GbpD(GEF) (1-1761bp of GbpD). PKBA-mCherry-FRB was constructed using two repeats of AGTGCTGGTGGT as linker between PKBA and mCherry. Biosensors PHplcδ-YFP, PTEN-GFP, PHcrac-YFP, and j-YFP were cloned into pCV5. Conventional RBD-GFP was used in Fig. 3e and Supplementary Fig. 6, while the longer version RBD (51–220 amino acids of Raf)-YFP[55] was used in Fig. 3a to enhance expression and reduce toxicity. PHcrac-YFP/LimE-RFP in the same pDM358 vector with due expression cassette was built using the shuttle vector pDM344. myr-FKBP-FKBP/myc-Rap1S17N (pDM358) was built similarly for experiments in Fig. 3d where dominant negative Rap1 was overexpressed.

For experiments where cells overexpressed two proteins, as presented in Fig. 1c–f, 3d, Supplementary Fig. 1b–f, 2, 3, and Supplementary Video 1, 2, the combination of pB18 and pDM358 were used. The cells were first transformed with mCherry-FRB-X in pB18 by electroporation and selected with G418 (20 μg/ml). Then the selected cells were transformed again with myr-FKBP-FKBP in pDM358 and cultured with G418 (20 μg/ml) and hygromycin (50 μg/ml).

For experiments where cells overexpressed three proteins, as presented in Fig. 1b, 2, 3a–c,f, 4, 5, Supplementary Fig. 1a, 5, 6g,h, 7 and Supplementary Video 3–9, 11, two alternative transformation systems were used. One is the combination of biosensor (pCV5), mCherry-FRB-X (pCV5), and myr-FKBP-FKBP (pCV5), and the cells were cultured with G418 (20 μg/ml). The other is the combination of biosensor (pCV5), mCherry-FRB-X (pCV5), and myr-FKBP-FKBP (pDM358), and the cells were cultured with both G418 (20 μg/ml) and hygromycin (50 μg/ml).

For experiments where cells overexpressed four proteins, as presented in Fig. 4c and Supplementary Video 10, FRB-Inp54p (pCV5), myr-FKBP-FKBP (pCV5), and PHcrac-YFP/LimE-RFP (pDM358) were used, and the cells were cultured with both G418 (20 μg/ml) and hygromycin (50 μg/ml).

For expereiments with inducible RasCQ62L in Supplementary Fig. 2d–e, we used the pDM359 vector, obtained from the DictyBase stock center. Cells were cultured in hygromycin (50 μg/ml) and RasCQ62L expression was induced with Doxycycline (10 μg/ml) for 12 hours before imaging.

Microscopy

Growth-phase cells were placed in an eight-well cover slip chamber and allowed to adhere for 10-15 min. Then the media was replaced with 450 μl DB. In order to add 5μM rapamycin in most experiments, firstly 1 μl stock (10 mM in DMSO) was diluted with 200 μl DB, then 50 μl of the solution was added to the 450 μl in the chamber drop wise. Imaging usually started no later than 30 min after DB replacement. For experiments presented in Fig. 3d, cells were incubated with 40 μM LY294002 or 20 μM PP242 or both for 80–100 min in DB before imaging started. And for experiments presented in Fig. 3e,f, Fig. 4e–g and Supplementary Video 11, cells were treated with 5 μM LatrunculinA in DB for about 30 min.

Zeiss LSM780 single-point laser-scanning microscope (Zeiss AxioObserver with 780-Quasar confocal module; 34-channel spectral, high-sensitivity gallium-arsenide phosphide detectors, GaAsP), was used for confocal image acquisition.

Total Internal Reflection Fluorescence (TIRF) microscopy was carried out with a Nikon Eclipse TiE microscope illuminated by an Ar laser (YFP) and diode laser (mCherry). Images were acquired by a Photometrics Evolve EMCCD camera controlled by Nikon NIS-Elements.

Phase illumination on a Zeiss Observer.Z1 inverted microscope equipped with a 20×/0.3 air objective or a 40×oil objective was used for phase image acquisition.

Image analysis

Images were analyzed with ImageJ (NIH) and MATLAB (MathWorks, Natick, MA).

was measured by manually selecting the cytosolic portion of the cells (Fig. 1a, 5a,b). The responses to cAMP stimulation (Fig. 5c,d and Supplementary Fig. 7) were calculated from the decrease in the mean cytosolic fluorescence. Each value point came from 1 minus the mean peak value (within 30 sec after cAMP) of PHcrac cytosolic intensity in individual condition, and then was normalized by the response to 100 nM cAMP.

To get the in Fig. 1b and supplementary Fig. 1a, the cell image was first converted to binary after the threshold was adjusted to cover all the pixels of the cell, yielding image A. Then “Erode” was used to generate image B and “Erode” was applied a second time to generate image C. Image B was subtracted from image A by “Image Calculator” to get Image D. Image D was multiplied with the original cell image to get the pixels of membrane portion. Image C was multiplied with the original cell image to get the cytosol portion. The mean intensity of the membrane portion and cytosol portion were measured and the ratio was calculated.

in Fig. 1e were measured by using the “Analyze Particle” function on threshold adjusted cell images, and were normalized according to each mean cell area before rapamycin addition.

in Fig. 1d were generated using the “Find Edges” and “Time-Lapse Color Coder” in ImageJ. The overlays in Fig. 2b,d,f, and Supplementary Fig. 5a,c were generated by “Time-Lapse Color Coder” in Image J. Cell outlines in Fig. 7a,c,d,e, and Supplementary Video 8, 9 were obtained manually and imposed on the images.

The ImageJ 3D Viewer plugin was used to stack the frames from videos to create the in Fig. 2a and Fig. 7b, as previously described[11]. The resampling factor is set to 1 to avoid blurring of activities between frames.

The presented in Fig. 2–4 and Supplementary Fig. 5, 6h were created using the image processing toolbox of MATLAB (MathWorks, Natick, MA) as previously described[11]. Colors were assigned linearly in the same fashion across all kymographs, with blue indicating the lowest intensity and red the highest.

presented in Fig. 1d, and Supplementary Fig. 1 were created using the Gradientech Tracking Tool, by clicking the centroid of each cell in each frame (Videos were acquired every 12 sec).

presented in Fig. 4d was measured in ImageJ. The TIRF images were adjusted to stand out edges of propagating waves, then the distance between two fixed locations along the propagating direction was measured, as well as the time that the wave edge crossed these two locations. Such measurement was carried out three times for each cell, and each data point in the figure represents the averaged value. The lagging edges of the waves were measured for oscillators.

Distribution of biosensor activity sizes

The sizes of the activity patches, as presented in Fig. 3e,f, 4e and Supplementary Fig. 6, were computed using a custom MATLAB script. The fluorescent images were first segmented using adaptive threshold, followed by morphological opening (MATLAB commands: adaptthresh, imbinarize, imerode, bwareaopen, imdilate, imfill and imclearborder, were applied in that order). The round, latrunculin-treated cells were then identified using the circular Hough transform (MATLAB command: imfindcircles). The intensity around the perimeter was then computed, smoothed with a moving average filter of length three, and compared to the mean and standard deviation of the intensity of the interior region. A patch was defined as a contiguous region around the perimeter that is at least two standard deviations higher than the interior’s mean. The “activity size” was computed by calculating the fraction of the perimeter occupied by total patches.

Assignment of migratory modes

Each fraction in Fig. 1f was obtained by manually following every cell appearing in the same field (354.25μm×354.25μm) in a 10-min time window. Fan-shaped cells resemble fish keratocytes, with protrusions persistently taking over one side of the cell and migrate almost in a straight line. Oscillatory cells exhibit repeatedly spreading and crunching morphologies. The rest were assigned as amoeboid mode. Based on corresponding videos, migratory modes of every cell, including the ones that migrated in or out of the field, are identified and prorated according to the time spent in the window. A rougher quantification was used for Fig. 3d, 7f, Supplementary Fig. 2, 3 in which a cell was identified as a responding cell (fan-shaped or oscillatory) as long as it showed any spreading morphology in selected time window, and in which fans and oscillators were not separated.

Automatic assignments were carried out in MATLAB (MathWorks, Natick, MA) as detailed in Supplementary Fig. 4. Fans were defined as cells that migrated perpendicular to their long axis, oscillators as cells displayed a large variance in area, and amoeboid cells as cells with small variance in area.

Excitability of PHcrac to cAMP stimulation

Cells associated with Fig. 5 and Supplementary Fig. 7 were developed in suspension, by pulsing with 50–100 nM cAMP every 6 min at 2 × 107 cells/ml for about 4 hrs. Then cells in 450 μl DB was plated in one 8-well chamber, and incubated with 5μM LatA for about 30 min. Confocal imaging was performed during which 50 μl cAMP solution (10-fold dilution) was added. Images were acquired every 2.5 sec. Then in order to get rid of cAMP all media in this well was carefully removed, followed by immediately replenishing fresh DB media with 5 μM LatrunculinA and 5 μM rapamycin. The chamber and the cells could not be perturbed in order to keep observing the same cells in the same field. After incubation for 10-15 min, confocal imaging was performed again with stimulation by the same concentration of cAMP. The same process was repeated using other wells of cells to get data of other concentrations of cAMP stimulation.

Cell migration assay

amiB- and pikI- chemotaxis assays were carried out using clear-view migration plates Incucyte, Ann Arbor, MI. Briefly, cells were ceded at a density of 4×103 cells in the upper chamber of the cell migration plate and allowed to chemotaxis to folic acid at the indicated concentrations in the lower chamber for 3h. Chemotaxis was scored using Incucyte Zoom software, Essen Biosciences, Ann Arbor, MI.

Simulation of excitable signal transduction system

The excitable network is described by three interacting states (Fig. 8a). State F is mutually inhibitory to state B and initiates a negative feedback loop through state R. The system can be described by the following three partial differential equations: Note that if the coefficients for the equation for B are large, then we can approximate this state by its quasi-steady-state value Replacing this into the equation for F leads to the second order system where and , resulting in a two-species activator-inhibitor system as previously considered[16]. Both the components in this system can diffuse spatially, with diffusion coefficients D and D respectively. The signals U, U and P are inputs to the excitable system. They refer to the basal level of activation, the stochastic component and the polarization component respectively (explained in more detail below). In this case there is no external stimulus[27], as we are interested in the system response to intrinsic noise. However, for the dose-response – where we simulated the response to different levels of cAMP, the U term was replaced with the cAMP concentration. The stochastic component is modeled as a zero mean, white noise process with a controllable variance: to modulate the amount of noise in the system.

Additionally, we incorporate a polarity mechanism given by P = Z – W, where the individual components are implemented as local excitation, global inhibition mechanisms[27]: For simplicity, we let DW be sufficiently high that W is spatially independent. The polarization module is activated by signal σpro, which represents the stress from actin polymerization and is proportional to R. The complete system is shown in Supplementary Fig. 8a.

Varying the slope of the R-nullcline gives rise to different cell behaviors

If we fix P, set each of the partial equations to zero and solve for R as a function of F we obtain the two nullclines of the activator-inhibitor system. The parameter g determines the slope of the R-nullcline (Supplementary Fig. 8b). The difference between the equilibrium state and the minimum of the F-nullcline roughly dictates the threshold of the system, where decreasing the value of g causes a decrease in threshold and vice versa. In amoeboid cells (slope: game=5.5), if noise is sufficient to perturb the steady state across the threshold, a large excursion in the phase plane takes place before the system settles back to equilibrium (Supplementary Video 13). The contribution of the time- and spatially-dependent P causes the R-nullcline to oscillate around the steady state as the excursions take place. It subsequently causes the threshold to lower at the front of cells and increase at the back thus making it unlikely for the back to get excited, resulting in the cell becoming polarized in one unique direction. Fan-shaped cells manifest at a lower threshold (slope: gfan=4.6), with strong polarity at the front, causing the nullcline to shift across the minimum resulting in oscillations at the front (Supplementary Video 13). Reducing this slope further (slope: gosc=2.2) leads to oscillators as the R-nullcline shifts completely to the right of the minimum of F and creates an unstable equilibrium state for the whole cell, in which case, the slightest perturbation drives the system into periodic large excursions, without the system ever reaching steady state again, resulting in whole cell expansions and contractions (Supplementary Video 13).

To see the whole spectrum of cell behaviors on one kymograph (Fig. 6g), we held the slope constant at 6 until t=100, after which it was lowered exponentially so that, at t=600, the slope was 2, thereby encompassing the required thresholds for amoeboid cells, fans and oscillators, and enabling us to visualize their transition. To determine the probabilities of each migratory mode at different time points (Fig. 6h), we used the random number generator in MATLAB (Mathworks, Natick, MA) to obtain random slopes between 6 and 2, and then categorized the responses at each time point. The dose response curve (Fig. 6i) was similarly simulated where the higher thresholds were obtained by generating random slope values between 6.5 and 5.5 (blue), and lower threshold values between 5 and 2 (red), and then normalizing the peak response to each threshold by the maximum simulated response at 100 nM cAMP.

Simulating cell morphology

To determine the effect of these changes on morphology we simulated cell behavior using level set methods (LSM), as previously described[56]. Briefly, in the LSM the cellular boundary is defined as the zero-level set of a potential function φ(x,t), x∈R[2] whose evolution of the cell is described by the equation: where v(x,t) describes the local velocity of the potential function. To obtain this velocity we apply different stresses on the cell and use a viscoelastic mechanical model of the cell to determine the local velocity. In our case we use: where σtot is the total stress applied on the cell, xmem and xcor are the local displacements of the membrane and cortex, respectively, and K, D and B are viscoelastic components of the cell describing the elasticity (K) and viscosity (D) of the membrane, and the viscosity of the (B) of the cytoplasm. The velocity in equation (1) is given by v=∂xmem/∂t.

The total stress (σtot) acts normal to the cell boundary and incorporates surface tension (σten= γκ(x), where γ is the local cortical tension and κ is the local curvature), volume conservation (σvol = karea(A(t)−A), where A is the surface area enclosed by the cell boundary either at time t or initially) and protrusive forces (σpro= σ0R(θ)).

Parameters used in these simulations are provided in Supplementary Table 1. The model and all simulations are implemented using MATLAB. The PDEs for the signaling system were solved by representing the cell boundary as a one-dimensional system, discretized in space using 314 points. The two dimensional simulations were carried out discretizing the system similarly on a 300×300 point lattice, using the same parameters as the one dimensional system. Spatial diffusion terms, which contain the second derivatives, were approximated by central differences in space, subsequently converting the partial differential equations to ordinary differential equations. The time step for simulation was set to 0.01 seconds. The solutions of the stochastic differential equations were obtained using the SDE toolbox for MATLAB[57]. After solving the concentrations of all species, we compute the protrusive force using the concentration for R and use this force to update the potential function in the LSM. The potential function is solved on a Cartesian grid with spatial discretization of 11 points per μm. The level set simulations were carried out using the Level Set Toolbox for MATLAB using the first forward Euler method[58].

Statistics and reproducibility

Statistical significance and P values were determined using the two-tailed t-test in Fig. 5c, 7g,h, and Supplementary Fig. 2c, 3c,d, 6f, 7. One-way ANOVA with post-hoc Tukey HSD test was used in Fig. 3d and 7f. Wilcoxon-Mann-Whitney rank sum test was used in Fig. 4d. Mean values ± either s.d. or s.e.m. were reported as indicated in the figure legends, except in the scatter plots of Fig. 4d where median was shown.

Representative images were presented in Fig. 1a–c, 2b,d,f, 3a–c, 4a–c,e, 5a,b, 7a,c,d,e, and Supplementary Fig. 2a,b,d, 3a,b, 6a which have been repeated in more than 3 independent experiments.

Code availability

Computational codes are available from the corresponding author on reasonable request.

Data availability

Source data for Fig 3d, 7f–h, and Supplementary Fig 2c, 3c,d, 6f have been provided as Supplementary Table 2. All other data supporting the findings of this study are available from the corresponding author on reasonable request.