Cells use molecular working memory to navigate in changing chemoattractant fields

  1. Akhilesh Nandan
  2. Abhishek Das
  3. Robert Lott
  4. Aneta Koseska  Is a corresponding author
  1. Department of Systemic Cell Biology, Max Planck Institute of Molecular Physiology, Germany

Abstract

In order to migrate over large distances, cells within tissues and organisms rely on sensing local gradient cues which are irregular, conflicting, and changing over time and space. The mechanism how they generate persistent directional migration when signals are disrupted, while still remaining adaptive to signal’s localization changes remain unknown. Here, we find that single cells utilize a molecular mechanism akin to a working memory to satisfy these two opposing demands. We derive theoretically that this is characteristic for receptor networks maintained away from steady states. Time-resolved live-cell imaging of Epidermal growth factor receptor (EGFR) phosphorylation dynamics shows that cells transiently memorize position of encountered signals via slow-escaping remnant of the polarized signaling state, a dynamical ‘ghost’, driving memory-guided persistent directional migration. The metastability of this state further enables migrational adaptation when encountering new signals. We thus identify basic mechanism of real-time computations underlying cellular navigation in changing chemoattractant fields.

Editor's evaluation

This paper addresses how cells can robustly maintain direction during movement by ignoring noise in concentration gradients while also being able to adapt to new signals in those gradients. The authors study this tension in EGFR signaling by postulating a form of cellular memory in a theoretical framework based on dynamical systems and bifurcation theory. The authors also carry out experiments that raise further interesting questions. This paper will be of interest to scientists of all stripes working on cell motility and for theorists who take a dynamical systems view of biological phenomena.

https://doi.org/10.7554/eLife.76825.sa0

eLife digest

If we are injured, or fighting an infection, cells in our body migrate over large distances to the site of the wound or infection to act against any invading microbes or repair the damage. Cells navigate to the damaged site by sensing local chemical cues, which are irregular, conflicting and change over time and space. This implies that cells can choose which direction to travel, and stick to it even if the signals around them are disrupted, while still retaining the ability to alter their direction if the location of the signal changes. However, how cells are able to effectively navigate their way through this field of complex chemical cues is poorly understood.

To help resolve this mystery, Nandan, Das et al. studied the epidermal growth factor receptor (EGFR) signaling network which controls how some cells in the body change shape and migrate. The network is activated by specific chemical cues, or ligands, binding to EGFR proteins on the cell surface. The receptors then join together to form pairs, and several tags known as phosphate groups are added to each molecule. This process (known as phosphorylation) switches the receptor pair to an active state, allowing EGFR to relay signals to other proteins in the cell and promote the activity of receptors not bound to a ligand. The phosphorylation state of EGFRs is then modulated over time and across the cell by a network of enzymes called phosphatases which can remove the phosphate groups and switch off the receptor.

To study EGFR phosphorylation dynamics in human cells, Nandan, Das et al. imaged individual cells over time using a microscope. This data was then combined with a mathematical model describing the EGFR signaling network and how cells change their shape over time.

The experiment revealed that the phosphate groups attached to EGFR are not removed immediately when the chemical cue is gone. Instead, the active state is transiently maintained before complete inactivation. This had the effect of encoding a short-lived memory in the signaling network that allowed the cells to continue to migrate in a certain direction even when chemical cues were disrupted. This memory state is dynamic, enabling cells to adapt direction when the cue changes location.

The findings of Nandan, Das et al. reveal the underlying mechanism for how cells decipher complex chemical cues to migrate to where they are needed most. The next steps to follow on from this work will be to understand if other receptors involved in migration work in a similar way.

Introduction

Directed chemotactic behavior relies on generating polarized signaling activity at the plasma membrane of the cell that is translated to an elongated cell shape, and subsequent persistent migration in the direction of the signal. Experimental observations have shown that cells as diverse as social amoeba, neutrophils, leukocytes, fibroblasts, and nerve cells maintain the acquired orientation even when signals are disrupted or noisy (Parent and Devreotes, 1999; Foxman et al., 1999; Ridley et al., 2003). However, not only do they respond robustly to dynamic gradients, they can also adapt the migrational direction by integrating and resolving competing spatial signals, or prioritizing newly encountering attractants (Jilkine and Edelstein-Keshet, 2011; Skoge et al., 2014; Albrecht and Petty, 1998). This suggests that cells likely memorize their recent environment. Numerous models based on positive feedbacks, incoherent feed-forward, excitable or Turing-like networks have been proposed to describe how polarized signaling activity of cell-surface receptors and/or downstream signaling component such as members of the Rho GTPase family can arise (Levchenko and Iglesias, 2002; Levine et al., 2002; Mori et al., 2008; Goryachev and Pokhilko, 2008; Beta et al., 2008; Xiong et al., 2010; Trong et al., 2014; Halatek and Frey, 2018). This polarized activity in turn controls actin and myosin dynamics, and thereby cell migration. Conceptually, the underlying dynamical principles of the proposed models are similar, and can be understood as switching from the stable state of basal- to the stable polarized-signaling steady state in presence of guiding external cues. However, they can account either for sensing and adaptation to non-stationary stimuli or for long-term maintenance of polarized signaling activity, but not both. Thus, how cells process the information from a changing chemoattractant field in real time for long-range navigation remains unknown.

We propose a shift in the conceptual framework, describing theoretically that efficient navigation can be achieved when the polarized signaling state of the receptor network is transiently stable. This is fulfilled in the presence of dynamical ‘ghosts’ at a unique dynamical transition, which we demonstrate in the EGFR signaling network dynamics using a mathematical model, as well as quantitative live-cell imaging of polarized EGFR signaling. We show with a physical model of the cell and migration experiments using microfluidics, that cells generate memory of encountered signals through the ‘ghost’ state, translating it to memory in polarized shape changes and directional migration. Due to the metastability of the ‘ghost’ state, cells can also easily adapt their migration direction depending on the changes in signal localization. We therefore describe a basic mechanism of real-time cellular navigation in complex chemoattractant fields.

Results

Dynamical mechanism of navigation in non-stationary environments

We conjectured that only dynamically metastable receptor signaling states can enable both transient stability of polarized signaling as necessary for robust, memory-guided migration in noisy fields, as well as rapid adaptation of its direction when signals vary in space and time. Our hypothesis is that this can be achieved if biochemical systems are maintained outside, but in the vicinity of the polarization steady state. We therefore approached the problem using the abstract language of dynamical systems theory, where the characteristics of any process directly follow from the type of dynamical transitions, called bifurcations, through which they emerge (Strogatz, 2018).

Directed migration relies on a polarized representation of the directional signal, requiring a reliable mechanism for signal-induced transition from a non-polarized symmetric, to a polarized receptor signaling state, and subsequently polarized cell shape. This transition is thus a symmetry-breaking transition, and we propose that a pitchfork bifurcation (PB, Koseska et al., 2013; Strogatz, 2018) satisfies the necessary dynamical conditions (Figure 1A, Figure 1—figure supplement 1A). Transient memory on the other hand is a unique characteristic of another bifurcation, a saddle-node (SN) bifurcation, that characterizes a transition between stable and unstable steady states. When the SN and thereby a stable steady-state is lost, for example upon signal removal, a remnant or a dynamical ‘ghost’ of the stable state emerges (Strogatz, 2018). These ‘ghost’ states are dynamically metastable and transiently maintain the system in the vicinity of the steady state (Figure 1A, Figure 1—figure supplement 1A). Necessary for manifestation of the ‘ghost’ state is organization at criticality, before the SN. We have previously examined both theoretically and experimentally, the response of receptor networks under uniform growth factor stimulation and determined that the concentration of receptors on the cell membrane regulate the organization of the system at criticality (Stanoev et al., 2018; Stanoev et al., 2020). The features of both bifurcations, cell polarization under spatial cues and a transient memory of this polarization in absence of the cue, will be unified for a sub-critical PB, as it is stabilized via SNPB s. We thus propose that organization at criticality - in the vicinity of a SNPB (gray shaded area in Figure 1—figure supplement 1A; details discussed in Materials and methods), renders a minimal mechanism for cellular responsiveness in changing environments.

Figure 1 with 2 supplements see all
In silico manifestation of metastable polarized membrane signaling, as a mechanism for sensing changing spatial-temporal signals.

(A) Dynamical mechanism: sub-critical pitchfork bifurcation (PB) determines stimulus-induced transition (arrow) between basal unpolarized and polarized receptor signaling state, whereas the associated saddle-node through which the PB is stabilized (SNPB) gives rise to a ‘ghost’ memory state upon signal removal for organization at criticality (before the SNPB). See Figure 1—figure supplement 1A and Methods for detailed description of these transitions. (B) Scheme of the EGFR-PTP interaction network. Ligandless EGFR (Ep) interacts with PTPRG (PRG) and PTPN2 (PN2). Liganded EGFR (E-Ep) promotes autocatalysis of Ep. Causal links - solid black lines; curved arrow lines - diffusion, PM - plasma membrane, ER- endoplasmic reticulum. See also Figure 1—figure supplement 1B. (C) Signal-induced shape-changes during cell polarization. Arrows: local edge velocity direction. Zoom: Viscoelastic model of the cell - parallel connection of an elastic and a viscous element. Ptotal: total pressure; v: local membrane velocity; l: viscoelastic state. Bold letters: vectors. Cell membrane contour: [0,2π]. (D) Top: In silico evolution of spatial EGF distribution. Bottom: Kymograph of Ep for organization at criticality from reaction-diffusion simulations of the network in (B). Triangle - gradient duration. (E) Corresponding exemplary cell shapes with color coded Ep, obtained with the model in (C). (F) Top: Temporal profiles Ep (black) and E-Ep (gray). Green shaded area: EGF gradient presence. Bottom: State-space trajectory of the system with denoted trapping state-space areas (colored) and respective time-scales. See also Figure 1—video 1. Thick/thin line: signal presence/absence. (G) Quantification of in silico cell morphological changes from the example in (E). Triangle - gradient duration. (H) Left: same as in (G), only when stimulated with two consecutive dynamic gradients (triangles) from same direction. Second gradient within the memory phase of the first. See also Figure 1—figure supplement 1D. Right: the second gradient (orange triangle) has opposite direction. See also Figure 1—figure supplement 1E. Dashed line: curve from (G). Mean ± s.d. from n=3 is shown. Parameters: Materials and methods. In (D-H), green(orange)/red lines: stimulus presence/absence.

We described this conjecture mathematically for a general reaction-diffusion model representing the signaling activity on the plasma membrane of a cell, U(x,t)t=F(U)+D2U(x,t), with U being the vector of local densities of active signaling components, D - diffusion constants and F accounting for all chemical reactions. Our theoretical analysis shows that a PB exists if, for a spatial perturbation of the symmetric steady state (Us) of the form U(x,t)=Us+δU(x)eλt, the conditions δU(-x)=-δU(x) and the limit limλ0Fλ=det(J)=0 are simultaneously fulfilled (Materials and methods). This implies that the linearized system has zero-crossing eigenvalues (λ) associated with the odd mode of the perturbation (Paquin-Lefebvre et al., 2020). To probe the sub-critical transition and therefore the necessary organization at criticality, a reduced description in terms of an asymptotic expansion of the amplitude of the polarized state (ϕ) must yield the Landau equation dϕdt=c1ϕ+c2ϕ3-c3ϕ5, guaranteeing the existence of SNPB (see Materials and methods for derivation).

These abstract dynamical transitions can be realized in receptor signaling networks with different topologies and are best analyzed using computational models, whose predictions are then tested in quantitative experiments on living cells. To exemplify the above-mentioned principle, we use the well-characterized Epidermal growth factor receptor (EGFR) sensing network (Reynolds et al., 2003; Baumdick et al., 2015; Stanoev et al., 2018). It constitutes of double negative and negative feedback interactions of the receptor, EGFR (Ep) with two enzymes, the phosphatases PTPRG (PRG) and PTPN2 (PN2, Figure 1B, Figure 1—figure supplement 1B), respectively. Ep and PRG laterally diffuse on the membrane and inhibit each-other’s activities (see Materials and methods for the molecular details of the network). The bidirectional molecular interactions between EGFR and the phosphatases can be mathematically represented using mass action kinetics, giving a system of partial differential equations (PDE) that describes how the dynamics of the constituents evolves in time and space (Equation 14 in Materials and methods). Applying a weakly nonlinear stability analysis (Becherer et al., 2009) to this system of equations shows that the EGFR phosphorylation dynamics undergoes a symmetry-breaking transition (PB) as outlined above (proof in Materials and methods, Figure 1—figure supplement 1C). The PB generates a polarized state that is represented as a inhomogeneous steady state (IHSS) - a combination of a high receptor phosphorylation at the cell front and low in the back of the cell (schematically shown in Figure 1A, Figure 1—figure supplement 1A). This is contrary to a bistable system, where the polarized signaling state would be manifested by two steady states, high and low protein phosphorylation in the front and back of the cell, respectively (Beta et al., 2008). This profiles PB as a robust mechanism of cell polarization. Polarized EGFR signaling on the other hand, will lead to reorganization of the cortical actomyosin cytoskeleton by regulating members of the Rho GTPase family, thereby inducing signal-dependent cell shape changes and subsequent migration (Chiasson-MacKenzie and McClatchey, 2018; Ridley and Hall, 1992). In order to link signaling activity with morphodynamics, we modeled the cell as a viscoelastic cortex surrounding a viscous core (Yang et al., 2008) (Materials and methods), where EGFR signaling dynamics affects cell shape changes through the protrusion/retraction stress and the viscoelastic nature of the cell membrane (Figure 1C).

We first fixed the total EGFR concentration on the cell membrane to a value that corresponds to organization at criticality, and investigated the response of the in silico cell to gradient stimulus. In the absence of stimulus, basal EGFR phosphorylation is uniformly distributed along the cell membrane rendering a symmetrical cell shape (Figure 1D and E). Introducing dynamic gradient stimulus in the simulation (slope changes from steep to shallow over time, Figure 1D, top) led to rapid polarization of EGFR phosphorylation in the direction of the maximal chemoattractant concentration, generating a cell shape with a clear front and back. The polarized signaling state was maintained for a transient period of time after removal of the gradient, corresponding to manifestation of memory of the localization of the previously encountered signal (Figure 1D and E; temporal profile Figure 1F, top). The prolonged polarized state does not result from remnant ligand-bound receptors (E-Ep) on the plasma membrane, as they exponentially decline after signal removal (Figure 1F, top). The memory in polarized signaling was also reflected on the level of the cell morphology, as shown by the difference of normalized cell protrusion area in the front and the back of the cell over time (Figure 1G). Plotting the trajectory that describes the change of the state of the system over time (state-space trajectory, Figure 1F bottom) shows that the temporal memory in EGFR phosphorylation polarization is established due to transient trapping of the signaling state trajectory in state-space, a property of the metastable ‘ghost’ state (Stanoev et al., 2020; Strogatz, 2018) through which the system is maintained away from the steady state. The simulations show that there are two characteristic time-scales present in the system: slow evolution of the system’s dynamics in the ‘ghost’ state due to the trapping, and fast transitions between the steady states (Figure 1—video 1). This emergence of the slow time-scale is another hallmark of systems organized at criticality. What is crucial here however, is that the trapping in the dynamically metastable memory state does not hinder sensing of, and adapting to subsequent signals. The cell polarity is sustained even when the EGF signal is briefly disrupted (Figure 1H left, Figure 1—figure supplement 1D), but also, the cell is able to rapidly reverse direction of polarization when the signal direction is inverted (Figure 1H right, Figure 1—figure supplement 1E).

We next chose in the simulations a higher EGFR concentration on the membrane, such that the system moves from criticality to organization in the stable polarization state (magenta lines, Figure 1—figure supplement 1C). In this scenario, even a transient signal induces switching to the polarized state that is permanently maintained, generating a long-term memory of the direction on the initial signal. Thus, the cell is insensitive to subsequent stimuli from the same direction, whereas consecutive gradients from opposite directions generate conflicting information that cannot be resolved (Figure 1—figure supplement 1F). Organization in the homogeneous, symmetric steady states on the other hand renders cells insensitive to the extracellular signals (Figure 1—figure supplement 1G,H). These response features for organization in the stable steady state regimes resemble the finding of the previously published models: such models cannot simultaneously capture memory in polarization along with continuous adaptation to novel signals, or require fine-tuning of kinetic parameters to explain the experimentally observed cell behavior (Levchenko and Iglesias, 2002; Levine et al., 2002; Mori et al., 2008; Goryachev and Pokhilko, 2008; Beta et al., 2008; Xiong et al., 2010; Trong et al., 2014). This demonstrates that organization at criticality, in a vicinity of a SNPB, is a unique mechanism for processing changing signals.

Cells display temporal memory in polarized receptor phosphorylation resulting from a dynamical ‘ghost’

To test experimentally whether cells maintain memory of the direction of previously encountered signals through prolonged EGFR phosphorylation polarization, and what is the duration of this effect, epithelial breast cancer-derived MCF7 cells were subjected for 1 hr to a stable gradient of fluorescently tagged EGF-Alexa647 (EGF647) with a maximal amplitude of 10 ng/ml applied from the top of the chamber in a computer-programmable microfluidic device (Figure 2A and B). EGFR phosphorylation at the plasma membrane was quantified during and for 3 hr after gradient wash-out (gradient wash-out established in 4–5 min) by determining the rapid translocation of mCherry-tagged phosphotyrosine-binding domain (PTBmCherry) to phosphorylated tyrosines 1086/1148 of ectopically expressed EGFR-mCitrine (EGFRmCitrine) using ratiometric imaging (Offterdinger et al., 2004; Materials and methods). Due to the low endogenous EGFR levels in MCF7 cells, the expression range of EGFRmCitrine was set to mimic the endogenous receptor range in the related MCF10A cell line, such that both cell lines have equivalent signaling properties of downstream effector molecules (Stanoev et al., 2018), and were therefore used in a complementary way in this study.

Figure 2 with 5 supplements see all
Molecular memory in polarized EGFRmCitrine phosphorylation resulting from dynamical state-space trapping is translated to memory in polarized cell shape.

(A) Scheme of microfluidic EGF647-gradient experiment; Zoom: single-cell measurables. Cell membrane contour [0,2π] (20 segments). PTB - phosphotyrosine binding domain, FP/star symbol - fluorescent protein, EGFRp- phosphorylated EGFRmCitrine. Remaining symbols as in Figure 1B. (B) Quantification of EGF647 gradient profile (at 60min, green) and after gradient wash-out (at 65min, red). Mean ± s.d., N=4. (C) Exemplary quantification of, Top: Spatial projection of EGF647 around the cell perimeter. Gaussian fit of the spatial projection is shown. Middle: single-cell EGFRp kymograph. Data was acquired at 1 min intervals in live MCF7EGFRmCitrine cells subjected for 60min to an EGF647 gradient. Other examples in Figure 2—figure supplement 1D. Bottom: respective spatial projection of EGFRp. Gaussian fit of the spatial projection is shown. Mean ± s.d. from n=20 cells, N=7 experiments in Figure 2—figure supplement 1C. (D) Average fraction of polarized plasma membrane area (mean ± s.d.). Single cell profiles in Figure 2—figure supplement 1G. (E) Quantification of memory duration in single cells (median±C.I.). In (D) and (E), n=20, N=7. (F) Top: Exemplary temporal profiles of phosphorylated EGFRmCitrine (black) and EGF647-EGFRmCitrine (gray) corresponding to (C). Bottom: Corresponding reconstructed state-space trajectory (Figure 2—video 1) with denoted trapping state-space areas (colored). Thick/thin line: signal presence/absence. d - embedding time delay. (G) Equivalent as in (F), only in live MCF7-EGFR SNPB cell subjected to 1 hr EGF647 gradient (green shading), and 3 hr after wash-out with 1 µM Lapatinib. Corresponding kymograph shown in Figure 2—figure supplement 2A. Mean ± s.d. temporal profile from n=9, N=2 in Figure 2—figure supplement 2B. Bottom: Corresponding reconstructed state-space trajectory with state-space trapping (colored) (Methods, Figure 2—video 2). (H) Averaged single-cell morphological changes (solidity, mean ± s.d. from n=20, N=7). Average identified memory duration (blue arrow): 40min. Top insets: representative cell masks at distinct time points. (I) Average solidity in MCF7EGFRmCitrine cells subjected to experimental conditions as in (G). Mean ± s.d. from n=9, N=2. Top insets: representative cell masks at distinct time points. In (F-I) green shaded area: EGF647 gradient duration; green/red lines: stimulus presence/absence. Brown line: Lapatinib stimulation. See also Figure 2—figure supplements 1 and 2.

Kymograph analysis of EGFRmCitrine phosphorylation at the plasma membrane of single cells showed polarization in a shallow gradient of EGF647 (as shallow as 10% between front and back of the cell; Figure 2C, Figure 2—figure supplement 1A-D). The direction of EGFRmCitrine phosphorylation polarization coincided with the direction of maximal EGF647 concentration around each cell (π/4 on average, Figure 2—figure supplement 1F). Only few cells manifested basal or symmetric EGFRmCitrine phosphorylation distribution upon gradient stimulation (Figure 2—figure supplement 1A, B, E). Plotting the fraction of plasma membrane area with polarized EGFRmCitrine phosphorylation showed cell-to-cell variability in the polarization kinetics, as well as the maximal amplitude of polarized EGFRmCitrine phosphorylation (Figure 2—figure supplement 1G), in contrast to the rapid EGFR polarization in the numerical simulations (Figure 1D). These differences likely results from the variable positioning of the cells along the gradient in the microfluidic chamber, as well as the variability of total EGFR concentrations in single cells. However, quantification of the polarization duration revealed that, similarly to the numerical predictions, the polarization persisted ~ 40 min on average after gradient removal ([4–159 min], Figure 2D and E).

The memory in EGFRmCitirne phosphorylation was also reflected in the respective single-cell temporal profiles (exemplary profile shown in Figure 2F, top). Reconstructing the state-space trajectory from this temporal profile using Takens’s delay embedding theorem (Takens, 1980) (Materials and methods) showed that before the fast transition to the basal state, the trajectory of the system was trapped in the vicinity of the polarized state (Figure 2F bottom, Figure 2—video 1). Despite the biological and technical noise that affect the measurement of the temporal EGFRmCitrine phosphorylation profile, and thereby the reconstruction of the state-space trajectory, both qualitatively resemble the equivalent numerical profiles (compare Figures 1F and 2F). In contrast, when cells were subjected to an ATP analog EGFR inhibitor (Björkelund et al., 2013) during gradient wash-out, the EGFRmCitrine phosphorylation response exponentially decayed, resulting in a clear absence of transient memory and respective state-space trapping (Figure 2G, Figure 2—figure supplement 2A, B, Figure 2—video 2). Since Lapatinib inhibits the kinase activity of the receptor, the dynamics of the system in this case is mainly guided by the dephosphorylating activity of the phosphatases. Implementing an equivalent of the Lapatinib inhibition in the numerical simulations by decreasing the autocatalytic EGFR activation rate constant after gradient removal verifies that the presence of memory in EGFR phosphorylation cannot be explained only by a dephosphorylation process (Figure 2—figure supplement 2C). This is also evident from the respective state-space trajectory, where the system directly transits from the polarized to the basal state, without intermediate state-space trapping (Figure 2—figure supplement 2D, Figure 2—video 3).

Fitting the experimentally measured single-cell temporal EGFRmCitrine phosphorylation profiles after gradient wash-out using an inverse sigmoid function (Methods) further corroborated that under Lapatinib treatment, phosphorylated EGFRmCitrine exponentially relaxed from the polarized to the basal state (Hill coefficient ≈ 1.28), with a half-life of approx. 10 min (Figure 2—figure supplement 2E, G). Under normal conditions however, the half-life was 30 min on average, reflecting that the phosphorylated EGFRmCitrine is transiently maintained in the metastable signaling state after gradient removal, before rapidly switching to the basal state (Hill coefficient ≈2.88, Figure 2—figure supplement 2F, G). Taken together, this analysis suggests that the memory in polarized EGFRmCitrine phosphorylation results from a dynamically metastable ‘ghost’ state, and not a slow dephosphorylation process.

In order to identify whether the memory in polarized EGFRmCitrine phosphorylation also enables maintaining memory of polarized cell morphology after gradient removal, we quantified the cellular morphological changes using solidity, which is the ratio between the cell’s area and the area of the convex hull. The average single-cell solidity profile over time showed that epithelial cells maintained the polarized cell shape for ~ 40 min after signal removal (Figure 2H, Materials and methods), which directly corresponds to the average memory duration in polarized EGFRmCitrine phosphorylation (Figure 2E). The exemplary quantification of the temporal evolution of the cell protrusion area in direction of the gradient showed equivalent results (Figure 2—figure supplement 2H corresponding to the profile in Figure 2C; memory duration ~ 43 min). In contrast, the absence of memory in EGFRmCitrine phosphorylation under Lapatinib treatment also resulted in absence of transient memory in polarized morphology after stimulus removal (Figure 2I). This establishes a direct link between memory in polarized receptor activity and memory in polarized cell shape.

Transient memory in cell polarization is translated to transient memory in directional migration

To test the phenotypic implications of the transient memory in cell polarization, we analyzed the motility features of the engineered MCF7EGFRmCitrine, as well as of MCF10A cells at physiological EGF concentrations. Cells were subjected to a 5 hr dynamic EGF647 gradient that was linearly distributed within the chamber, with EGF647 ranging between 25-0 ng/ml, allowing for optimal cell migration (Figure 3—figure supplement 1A, B). The gradient steepness was progressively decreased in a controlled manner, rendering an evolution towards a ∼50% shallower gradient over time (Figure 3—figure supplement 1B). Automated tracking of single-cell’s motility trajectories was performed for 14 hr in total. MCF7EGFRmCitrine, as well as MCF10A cells migrated in a directional manner toward the EGF647 source (Figure 3A- and Figure 3—figure supplement 1C, D - left, green trajectory parts). This directed migration persisted for transient period of time after the gradient wash-out (Figure 3A- and Figure 3—figure supplement 1C, D - left, red trajectory parts, Figure 3—video 1), indicating that cells maintain memory of the location of previously encountered source. After the memory phase, the cells transitioned to a migration pattern equivalent to that in the absence of a stimulus (Figure 3A right, Figure 3—figure supplement 1C, D middle). Uniform stimulation with 20 ng/ml EGF647 did not induce directed migration in either of the cell lines, although the overall migration distance was increased in accordance with previous findings (Brüggemann et al., 2021; Figure 3—figure supplement 1C, D, right). Quantification of the directionality of single cells’ motion, that is defined as the displacement over travelled distance, showed that for MCF10A cells, it was significantly higher during the gradient stimulation (5 hr) as compared to no- or uniform-stimulation case (Figure 3B). Moreover, the directionality estimated in the 9 hr time-frame after the gradient removal was greater than the one in continuous stimulus absence, corroborating that cells transiently maintain memory of the previous direction of migration.

Figure 3 with 5 supplements see all
Cells display memory in directional migration toward recently encountered signals.

(A) Left: representative MCF10A single-cell trajectories. Green - 5 hr during and red line - 9 hr after dynamic EGF647 gradient (shaded). Exemplary cell in Figure 3—video 1. Right: Same as in (A), only 14 hr in continuous EGF647 absence. Black dots: end of tracks. (B) Directionality (displacement/distance) in MCF10A single-cell migration during 14 hr absence (0 ng/ml; n=245, N=3) or uniform 20 ng/ml EGF647 stimulation (n=297, N=3); 5 hr dynamic EGF647 gradient (green) and 9 hr during wash-out (red; n=23, N=5). p-values: ***p≤ 0.001, two-sided Welch’s t-test. Error bars: median±95%C.I. (C) Top: Projection of the cells’ relative displacement angles (mean ± s.d.; n=23, N=5) during (green shaded) and after 5 hr dynamic EGF647 gradient. Green/red lines: stimulus presence/absence. Bottom: Kolmogorov-Smirnov (KS) test p-values depicting end of memory in directional migration (blue arrow, t=350min). KS-test estimated using 5 time points window. For (A-C), data sets in Figure 3—figure supplements 1D and 2A. (D) Representative in silico single-cell trajectories. Left: PB(t)RW: Persistent biased random walk, bias is a function of time (green/blue trajectory part - bias on). Right: RW: random walk. (E) Corresponding directionality estimates from n=50 realizations, data in Figure 3—figure supplement 2D. PRW: persistent random walk. p-values: ***p≤ 0.001, two-sided Welch’s t-test. Error bars: median±95%C.I. (F) Same as in (C), top, only from the synthetic PB(t)RW trajectories. (G) MCF10A single-cell trajectories quantified 5 hr during (green) and 9 hr after (brown) dynamic EGF647 gradient (shading) wash-out with 3 µM Lapatinib. n=12, N=5. See also Figure 3—video 2. (H) Directionality in single-cell MCF10A migration after gradient wash-out with (brown, n=12, N=5) and without Lapatinib (red, n=23, N=5). p-values: **p≤ 0.01, KS-test. Error bars: median±95%C.I. (I) Same as in (C), only for the cells in (G). See also Figure 3—figure supplement 2H.

This was also reflected in the projection of the cell’s relative displacement angles (cosθ) estimated along the gradient direction (π) at each time point (Figure 3—figure supplement 2A), representing the angular alignment of the cells to the source direction. The cellular migration trajectories aligned with the source direction (cosθ approached 1) during, and maintained this temporally after gradient removal, before returning to a migration pattern characteristic for stimulus absence or during uniform stimulation (cosθ ~ 0, Figure 3C top, Figure 3—figure supplement 2B). Calculating the similarity between the kernel density distribution estimate (KDE) of the angular alignment distributions at each point in the gradient series with that in continuous stimulus absence, showed that the distributions approach each other only ~ 50 min after the gradient removal (Figure 3C, bottom; Figure 3—figure supplement 2C). Additionally, the calculated similarity between the KDE distributions during the gradient (5 hr) and the 50min memory period further corroborated this finding (Figure 3—figure supplement 2C). The average memory phase in directional motility thus corresponds to the time-frame in which the memory in polarized EGFRmCitrine phosphorylation and cell shape is maintained (Figures 2E and 3C), indicating that the metastable signaling state is translated to a stable prolonged directed migration response after gradient removal.

To investigate whether the motility patterns during the gradient and the memory phase have equivalent characteristics, we fitted the motility data using a modified Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein, 1930; Svensson et al., 2017) and used the extracted migration parameters to generate synthetic single-cell trajectories (Materials and methods). In absence of stimulus, the cellular motion resembled a random walk process (RW: Figure 3D right, Figure 3—figure supplement 2D, E middle), persistent random walk (PRW) was characteristic for the uniform stimulation case (Figure 3—figure supplement 2D, E right), whereas biased PRW described the migration in gradient presence (PBRW, Figure 3D- and Figure 3—figure supplement 2D, left, green trajectory part). Extending the bias duration during the interval of the experimentally observed memory phase (PB(t)RW) was necessary to reproduce the transient persistent motion after gradient removal (Figure 3D- and Figure 3—figure supplement 2D, left, blue trajectory part; Figure 3E and F; Figure 3—figure supplement 2F).

To corroborate the link between memory in polarized receptor activity, memory in polarized cell shape and memory in directional migration, we also quantified the directional migration of MCF10A cells when subjected to Lapatinib during gradient wash-out (Figure 3G). The directionality after gradient removal was significantly lower than in the case without Lapatinib (Figure 3H), suggesting that cells rapidly switch to a RW migration pattern upon gradient wash-out due to the absence of memory in polarized EGFRmCitrine phosphorylation (Figure 2G,I). Thus, single-cell motility trajectories that closely resembled the experimentally observed ones could be mimicked with the PB(t)RW simulation, where the bias duration corresponded to the duration of the gradient (Figure 3—figure supplement 2E left, 2G). Quantification of the average cells’ relative displacement angles showed as well that cosθ approaches 0 exponentially after gradient removal (Figure 3I, Figure 3—figure supplement 2G), suggesting that majority of cells display absence of memory in directional migration under Lapatinib treatment.

In order to dissect better the cell-to-cell variability in this case, we also calculated memory duration from single cell cosθ profiles. For this, single-cell trajectories were first smoothed using Kalman filter (Materials and methods). The quantification showed that majority of the cells displayed absence of or shorter memory in directional migration, with a mean value of ∼25 min (Figure 3—figure supplement 3A, B, D). Since under Lapatinib treatment, EGFR phosphorylation rapidly decays (Figure 2G), this residual memory in some cells likely results from memory in cytoskeletal asymmetries, as previously suggested (Prentice-Mott et al., 2016). Without Lapatinib treatment however, the duration of memory estimated from single-cell cosθ profiles was of the order of 90 min (Figure 3—figure supplement 3A, C, E). If we therefore account in this case also the contribution of cytoskeletal memory, then the memory in directional migration which results from memory in polarized EGFR phosphorylation is on average ∼50 min, similar to the deduced values from the single-cell kymograph quantification (Figure 2E).

Molecular working memory enables cells to navigate in dynamic chemoattractant fields

To test whether the identified memory enables cellular navigation in environments where signals are disrupted but also change over time and space, we subjected cells in the simulations and experiments to a changing growth factor field. The field was generated by a sequence of signals, starting with a dynamic gradient whose steepness changed over time, and was temporary disrupted for a time interval shorter than the interval of memory in cell polarization. This was followed by a second static gradient in the same direction, that after an equivalent disruption period was followed by a third dynamic gradient in the opposite direction (Figure 4A). The in silico migration simulations showed that the cell can sense the initial dynamic gradient and polarizes in the direction of maximal attractant concentration, resulting in directed migration (Figure 4B, Figure 4—figure supplement 1A, Figure 4—video 1). The simulations also predicted that the memory of the previously encountered signal localization enables maintaining robust directional migration even when the signal was disrupted, while still remaining sensitive to the newly emerging signal from the opposite direction. The in silico cell rapidly adapted the orientation when encountering the third signal, demonstrating that the proposed mechanism can also account for prioritizing newly encountered signals. Such a dynamic memory which enables information of previous signals to be temporally maintained while retaining responsiveness to upcoming signals, and thereby manipulate the stored information, in neuronal networks is described as a working memory (Atkinson and Shiffrin, 1968).

Figure 4 with 8 supplements see all
Working memory enables history-dependent single-cell migration in changing chemoattractant field.

(A) Scheme of dynamic spatial-temporal growth factor field implemented in the simulations and experiments. Green(orange)/red: gradient presence/absence. (B) In silico cellular response to the sequence of gradients as depicted in (A), showing changes in EGFR activity, cellular morphology and respective motility trajectory over time. Trajectory color coding corresponding to that in (A), cell contour color coding with respective Ep values as in Figure 1E. Cell size is magnified for better visibility. See also Figure 4—figure supplement 1A, Figure 4—video 1. (C) Representative MCF10A single-cell trajectory and cellular morphologies at distinct time-points, when subjected to dynamic EGF647 gradient field as in (A) (gradient quantification in Figure 4—figure supplement 1E). Trajectory color coding corresponding to that in (A). See also Figure 4—video 4. Full data set in Figure 4—figure supplement 1F. (D) Projection of cells’ relative displacement angles (cosθ) depicting their orientation towards the respective localized signals. Mean ± s.d. from n=12, N=5 is shown. (E) Corresponding kernel density estimates (intervals and color coding in legend). p-values: ***, p≤0.001, ns: not significant, KS-test.

If the signal disruption is however longer than the duration of the working memory, the simulations demonstrated that cells cannot integrate the signals. In turn, cells respond to each signal individually, as the directional migration after the memory is lost, resulting in a shorter range migration trajectory (Figure 4—figure supplement 1B, Figure 4—video 2). On the other hand, if the system has a long-term memory, as resulting from organization in the stable polarized regime, the simulations showed that cellular adaptation to a changing gradient field is hindered (Figure 4—figure supplement 1C,D Figure 4—video 3). The initial dynamic gradient shifted the system to the stable polarization steady state where it was maintained on a long-term, such that sensitivity to upcoming signals from the same direction was hindered. Even more, the cell could not resolve the conflicting information from a subsequent gradient from the opposite direction, as the signals induced high receptor activity on the opposed cell sides, resulting in halted migration. These results therefore highlight the importance of working memory for generating memory-guided migration over long trajectories.

We next tested these predictions experimentally by establishing an equivalent dynamic EGF647 spatial-temporal field in a controlled manner in the microfluidic chamber, and quantified the migratory profile of MCF10A cells (Figure 4—figure supplement 1E). The MCF10A cells sensed the initial dynamic gradient field and migrated in the direction of increasing chemoattractant concentration, maintaining the directionality even when the signal was temporary disrupted. Despite the memory in cell polarization, cells remained responsive and adapted the duration of directional migration when presented with a second static gradient from the same direction, and subsequently prioritized the third, newly encountered signal with opposed orientation (exemplary trajectory in Figure 4C, Figure 4—video 4, Figure 4—figure supplement 1F,G). Thus, the predictions derived by the numerical simulations quantitatively captured that the proposed mechanism of navigation enables integration of, and adaptation to changes in signal localization. The distinction between the simulations and the experiments (Figure 4B and C) is only in the details of the migration pattern, since the PBRW migration mode was not included in the physical model of the cell for simplicity. The temporal memory in directional migration as well as the continuous adaptation of MCF10A cells to novel cues was also reflected in the projection of the cell’s relative displacement angles (Figure 4D). The thereby derived KDE distributions during the first and second gradient (5–245 min; 275–335 min, respectively), as well as the corresponding intervals in which the gradient has been disrupted (245–275 min; 335–365 min, respectively) were statistically similar to each other, demonstrating that cells maintain the direction of migration in the intermittent intervals when the gradient was interrupted (Figure 4E). Moreover, these distributions statistically differed from the one characterizing cellular migration in continuous EGF647 absence (w/o EGF647, distribution symmetrically distributed around cosθ=0). The presence of the third gradient from the opposite direction (365–605 min) on the other hand, induced a shift in the respective KDE distribution to negative cosθ values, reflecting that cells revert the direction of migration (established in ∼10 min). Furthermore, the reverse migration was maintained for approximately 20 min after wash-out of the third gradient (KDE 605–625 min). The statistical similarity between these two distributions demonstrates that cells also establish transient memory of the last detected signal, before reverting to a random walk migration mode (KDE 625–900 min similar to KDE w/o EGF647). These results therefore demonstrate that cells utilize molecular working memory to navigate in changing gradient fields.

Navigation in non-stationary fields, however, also necessitates integration of information, requiring active comparison during migration task execution. We therefore tested next numerically whether the identified organization at criticality enables resolving simultaneous gradients with different amplitudes from opposite sides, that temporally vary in time. In the simulations, the cell sensed the presence of both signals, as reflected in the respective increase in EGFR phosphorylation. However, the net polarization towards the higher amplitude gradient was dominant, resulting in a clear directional migration toward this signal (Figure 4—figure supplement 2A, B). After the gradient removal, the EGFR phosphorylation and the cell shape remained transiently polarized, manifesting memory of the recently encountered stronger signal that was translated to memory in directional migration, before the cell reverted to a random walk migration (Figure 4—video 5). In contrast, if the system has a long-term memory as resulting from organization in the stable polarized state, the simulations showed that EGFR phosphorylation increased almost equivalently with respect to both signals, despite the difference in signal amplitudes. This hindered the responsiveness of the cell such that migration could not be effectively exhibited (Figure 4—figure supplement 1C, D; Figure 4—video 6). These simulations therefore suggest that critical organization of receptor networks is in general crucial for performing complex cellular behavior that goes beyond simple stimulus-response associations.

Discussion

Our data establishes that mammalian cells use a mechanism of working memory to navigate in complex environments where the chemical signals are disrupted or vary over time and space. Previous observations of memory in directed migration have been explained through the presence of bistable dynamics, where the transition from the basal to the polarized steady state and vice versa (after a memory phase) is regulated by two finely tuned thresholds. The authors however did not identify potential molecular elements that store this information, or regulate the thresholds (Skoge et al., 2014). Similarly, the remaining proposed models of polarization also rely on steady-state description of the basal and polarized states (Levine et al., 2002; Mori et al., 2008; Goryachev and Pokhilko, 2008; Beta et al., 2008; Trong et al., 2014), and thereby cannot account for the rapid adaptation to changes in signal localization.

The mechanism of transient memory we report here is realized on a molecular level by a prolonged polarized phosphorylation state of a receptor tyrosine kinase. Dynamically, this state emerges for organization at criticality, where a slow-escaping remnant from the polarized state or a dynamically metastable ‘ghost’ state is generated, and endows cells with robust transient maintenance of directional migration after signal removal. Although the observed memory in directional migration is in part supported by the memory in cytoskletal asymmetries as previously suggested (Prentice-Mott et al., 2016), the memory in receptor signaling we identify here provides a crucial bridge between the rapid receptor phosphorylation/dephosphorylation events and the long-range cellular migration. In particular, the organization at criticality endows the system with a slow time-scale through which the prolonged receptor phosphorylation state can be maintained on average for ∼40–50 min after signal removal, which in turn maintains the polarized cell shape, and thereby directional migration in absence of a signal. Moreover, we have demonstrated that this memory arising from a metastable state uniquely ensures the ability of cells to quickly adapt to changes in the external environment.

Thus, our results suggest that in order to balance between a robust response and adaptation to novel signals, cell utilize an optimal receptor amount at the plasma membrane that corresponds to organization at criticality. The theoretical analysis suggest that the closeness of the receptor amount to the one corresponding to the critical transition is reflected in the memory duration. It can be therefore suggested that the observed variability in the experimentally identified memory length likely results from cell-to-cell variability in receptor concentration at the plasma membrane. Moreover, these results also suggest that a higher number of sensory units at the plasma membrane does not necessarily imply improved sensitivity of cells, but rather counterintuitively , leads to permanent memory of the initially encountered signal. This in turn will limit the cellular responsiveness to upcoming signal changes. It would be therefore of interest to study whether receptor networks are self-organized at criticality through an active sensing mechanism, or this feature has been fine-tuned through evolution, as a means for optimizing sensing and computational capabilities of cells.

Our work furthermore suggest that this general mechanism of a system poised at criticality can explain a wide range of biologically relevant scenarios, from the integration of temporally and spatially varying signals, to how extracellular information is transformed into guidance cues for memory-directed migration. Such memory-guided navigation is advantageous when migration must be realized over long and complex trajectories through dense tissues where the chemical cues are disrupted or only locally organized (Lämmermann et al., 2013). We have demonstrated here that the molecular working memory in cell polarization and therefore the capabilities of cells to navigate in a complex environment are an emergent feature of receptor networks.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Cell line (Homo sapiens)MCF-7ECACCCat.No.86012803
Cell line (Homo sapiens)MCF10AATCCCRL-10317
Recombinant DNA reagentEGFR-mCitrineBaumdick et al., 2015
Recombinant DNA reagentPTB-mCherryFueller et al., 2015
Recombinant DNA reagentcCbl-BFPFueller et al., 2015
Peptide, recombinant proteinFibronectinSigma-AldrichF0895-1MG
Peptide, recombinant proteinCollagenSigma-AldrichC9791-50MG
Chemical compound, drugLapatinibCayman chemicalsCay11493-10
Chemical compound, drugHoechst 33,342Thermo Fisher Sc.62,249
Chemical compound, drugDulbecco’s modified Eagle’s medium (DMEM)PAN BiotechCat. # P04-01500
Chemical compound, drugMEM Amino Acids Solution (50 x)PAN BiotechCat. # P08 32,100
Chemical compound, drugPenicillin- StreptomycinPAN BiotechCat. # P06 07100
Chemical compound, drugFetal Bovine SerumSigma-AldrichCat. # F7524
Chemical compound, drugEGFSigma-AldrichCat. # E9644
Chemical compound, drugHydrocortisoneSigma-AldrichCat. #H-0888
Chemical compound, drugCholera toxinSigma-AldrichCat. #C-8052
Chemical compound, drugInsulinSigma-AldrichCat. #I-1882
Chemical compound, drugHorse SerumInvitrogen26050088
Chemical compound, drugFuGENE6PromegaE2691
Software, algorithmPythonPython software foundationRRID:SCR_008394
Software, algorithmMatlabMathWorksRRID:SCR_001622
Software, algorithmXPPAUThttp://www.math.pitt.edu/~bard/xpp/xpp.html
Software, algorithmTrackmatehttps://doi.org/10.1016/j.ymeth.2016.09.016
Software, algorithmFiji, ImageJhttps://doi.org/10.1038/nmeth.2019
OtherEGF-Alexa647Sonntag et al., 2014Prof. Luc Brunsveld, University of Technology, EindhovenMethods
OtherCellasic ONIX platesMerck ChemicalsM04G-02-5PKMethods

Cell culture

Request a detailed protocol

MCF7 cells (sex: female, ECACC, Cat. No. 86012803) were grown at 37 °C and 5% CO2 in Dulbecco’s Eagle’s medium (DMEM) (PAN-Biotech, Germany), supplemented with 10% inactivated Fetal Calf Serum (FCS) (Sigma-Aldrich), 100 ng ml–1 L-Glutamine, 0.5 mg ml–1 non-essential amino acids, 100 µg ml-1 penicillin and 100 µg ml-1 streptomycin (PAN-Biotech, Germany). Serum starvation was performed by culturing the cells in DMEM supplemented with 0.5% FCS, 100 µg ml-1 penicillin and 100 µg ml-1 streptomycin (PAN-Biotech, Germany). MCF10A cells (sex: female, ATCC-CRL 10317) were grown at 37 °C and 5% CO2 in Mammary Epithelial Cell Growth Basal medium (MEBM from Lonza Pharma & Biotech), supplemented with 5% Horse Serum (HS) (Invitrogen), 20 ng ml–1 EGF (Sigma-Aldrich), 0.5 mg ml–1 hydrocortisone (Sigma-Aldrich), 100 ng ml–1 cholera toxin (Sigma-Aldrich), 10 µg ml-1insulin (Sigma-Aldrich), 100 µg ml-1 penicillin and 100 µg ml-1 streptomycin. Serum starvation was performed by culturing the cells in the DMEM supplemented with 0.5% HS, 0.5 mg ml–1 hydrocortisone (Sigma-Aldrich), 100 ng ml–1 cholera toxin (Sigma-Aldrich), 100 µg ml-1 penicillin and 100 µg ml-1 streptomycin. MCF7 and MCF10A cells were authenticated by Short Tandem Repeat (STR) analysis and did not contain DNA sequences from mouse, rat and hamster (Leibniz-Institut DSMZ). Cells were regularly tested for mycoplasma contamination using MycoAlert Mycoplasma detection kit (Lonza).

Transfection and cell seeding

Request a detailed protocol

For EGFRmCitrine polarization experiments, 2.5 × 105 MCF7 cells were seeded per well in a six-well Lab-Tek chamber (Nunc) until 80% confluence was reached. After 9–10 hr of seeding, transient transfection was performed with a total of 1 µg of plasmids (EGFRmCitrine, PTBmCherry and cCblBFP at ratio 4:3:4 by mass) using FUGENE6 (Roche Diagnostics) transfection reagent and Opti-MEM (Gibco - Thermo Fisher Scientific) according to manufacturer’s procedure. All plasmids were generously provided by Prof. P. Bastiaens, MPI of Molecular Physiology, Dortmund. Cells were incubated for 7–8 hr to allow the expression of the transfected proteins prior to experiments. To detach the cells, the growth media was discarded and cells were washed once with DPBS (PAN Biotech) before adding 100 µl Accutase (Sigma-Aldrich). After 10min incubation period at 37 °C and 5% CO2, fresh growth media was added, and the cell density and viability was measured using cell counter (Vi-CELL XR Cell Viability Analyzer System). After spinning down, the cells were diluted to 10 × 106 cells/ml. The M04-G02 microfluidic gradient plates (Merck Chemicals) were primed for usage by flowing cell culture growth media through the cell chamber for 5min and cells were subsequently seeded according to manufacturer’s instructions.

For migration experiments with uniform EGF647 stimulation, 6-well Lab-Tek plates were coated with Collagen (Sigma-Aldrich) in 0.1 M Acetic acid (Sigma-Aldrich) for MCF7 (100 µg cm-2), and Fibronectin (Sigma-Aldrich) in Phosphate-Buffered Saline (DPBS) (PAN-Biotech) for MCF10A cells (2 µg ml-1), and stored in incubator at 37 °C overnight for evaporation. Excessive media was removed and the wells were washed with DPBS before seeding cells. MCF7 cells were seeded and transfected as described above. In the case of MCF10A cells, 1 × 105 cells per well were used for seeding. For migration experiments with gradient EGF647 stimulation, MCF7 cells were transferred to the coated M04-G02 microfluidic gradient plates as described above. Before seeding, MCF10A cells were detached from 6 well Lab-Teks by discarding the growth media and washing once with DPBS (PAN Biotech) before adding 100 µl Accutase (Sigma-Aldrich). After 20–30 min incubation period at 37 °C and 5% CO2, fresh cell growth media was added, and the cell density and viability were measured using a cell counter (Vi-CELL XR Cell Viability Analyzer System). After spinning down, the cells were diluted to 2 × 106 cells/ml, and subsequently seeded in the microfluidic plates according to manufacturer’s instructions.

Reagents

For gradient quantification, Fluorescein (Sigma Aldrich) was dissolved in Dulbecco’s modified Eagle’s medium (with 25 mM HEPES, without Phenol Red) (PAN Biotech). Imaging media: DMEM without Phenol Red was mixed with 25mM HEPES. For nuclear staining, 20 mM Hoechst 33342 (Thermo Fisher Scientific) was mixed with DPBS and diluted to 2 µM working concentration. EGFR inhibitor Lapatinib (Cayman Chemical, Ann Arbor, MI) was solubilized in DMSO (Thermo Fisher Scientific) to a stock concentration of 5 mM and stored at –20 °C.

Confocal and wide-field microscopy

Request a detailed protocol

Confocal images were recorded using a Leica TCS SP8i confocal microscope (Leica Microsystems) with an environment-controlled chamber (Life Imaging Services) maintained at 37 °C and HC PL APO 63 x/1.2 N.A / motCORR CS2 water objective (Leica Microsystems) or a HC PL FLUOTAR 10 x/0.3 N.A. dry objective (Leica Microsystems). mCitrine, mCherry and Alexa647 were excited with a 470nm–670nm pulsed white light laser (Kit WLL2, NKT Photonics) at 514 min, 561 nm, and 633 nm, respectively. BFP and Hoechst 33342 (Thermo Fisher Scientific) were excited with a 405nm diode laser. The detection of fluorescence emission was restricted with an Acousto-Optical Beam Splitter (AOBS): BFP (425nm–448nm), Hoechst 33342 (425nm–500nm), mCitrine (525nm–551nm), mCherry (580nm–620nm), and Alexa647 (655nm–720nm). Transmission images were recorded at a 150–200% gain. To suppress laser reflection, Notch filter 488/561/633 was used whenever applicable. When using the dry objective for migration experiments, the pinhole was set to 3.14 airy units and 12-bit images of 512 × 512 pixels were acquired in frame sequential mode with 1 x frame averaging. When using the water objective for polarization experiments, the pinhole was fixed (1.7 airy units) for all channels. The Leica Application Suite X (LAS X) software was used.

Wide field images were acquired using an Olympus IX81 inverted microscope (Olympus Life Science) equipped with a MT20 illumination system and a temperature controlled C02 incubation chamber at 37 °C and 5% CO2. Fluorescence and transmission images were collected via a 10 x/0.16 NA air objective and an Orca CCD camera (Hamamatsu Photonics). Hoechst 33342 fluorescence emission was detected between 420nm–460nm via DAPI filter, mCitrine fluorescence emission between 495nm–540nm via YFP filter and Alexa647 fluorescence emission between 705nm–745nm via Cy5 filter. The xCellence (Olympus) software was used.

Gradient establishment for polarization and migration experiments

Request a detailed protocol

The CellAsic Onix Microfluidic Platform (EMD Millipore) was used for gradient cell migration and EGFRmCitrine phosphorylation polarization experiments. For EGFRmCitrine phosphorylation polarization experiments, 1 hr gradient stimulation was established using CellASIC ONIX2 software as follows. (i) Pre-stimulus: Imaging media was flowed from well groups 3 and 4 (CellAsic Onix Manual - https://www.merckmillipore.com/) at low pressure (2.5 kPa) for 5 min. (ii) Gradient establishment: After closing well group 3, pre-loaded EGF647 (10 ng mL–1) was flowed through well group 2 and imaging media from well group 4 at high pressure (15 kPa) for 15 min (iii) Gradient maintenance: The pressure was reduced to 10 kPa for 45 min. (iv) Washout: After closing well groups 2 and 4, imaging media was flowed from well groups 3 and 5 at high pressure (15 kPa) for 15 min and maintained at low pressure (7 kPa) for 165 min. For single gradient migration experiments, this protocol was modified as follows: in step (iii), gradient maintenance was done for 285 min. In step (iv), maintenance was at low pressure for 585 min. 30 ng mL–1 EGF647 was used. For polarization experiments with inhibitor, the same protocol as for polarization experiments was used, except well group 3 and 5 were filled with 1 µM Lapatinib solution and in step (i) well group 3 was kept closed. For single cell gradient migration experiment with inhibitor, 3 µM Lapatinib was used.

For migration experiments under subsequent gradient stimuli / gradient quantification, the following changes in the steps were used: (ii) well group 2 with 30 ng mL–1 EGF647/ 2.5 µM Fluorescein was used. (iii) The gradient maintenance was done for 225 min. (iv) Washout: imaging media was flowed from well groups 3 and 4 at high pressure (15 kPa) for 15 min and maintained at low pressure (7 kPa) for 15 min. (v) Second gradient establishment: After closing well group 3, EGF647 (30 ng ml-1)/ 2.5 µM Fluorescein was flowed from well group 2 and imaging media from well group 4 at high pressure (15 kPa) for 15 min. (vi) The second gradient formed was maintained by reducing the pressure to 10 kPa for 45 min. (vii) Washout: imaging media was flowed from well groups 3 and 4 at high pressure (15 kPa) for 15 min and maintained at low pressure (7 kPa) for 15 min. (viii) Third gradient establishment: After closing well group 4, EGF647 (30 ng mL–1) / 2.5 µM Fluorescein was flowed from well group 5 and imaging media from well group 3 at high pressure (15 kPa) for 15 min. (ix) The third reversed gradient was maintained by reducing the pressure to 10 kPa for 225 min. (x) Washout: imaging media was flowed from well groups 3 and 4 at high pressure (15 kPa) for 15 min and maintained at low pressure (7 kPa) for 285min.

Imaging EGFRmCitrine phosphorylation polarization and single cell migration

Request a detailed protocol

Transfected MCF7EGFRmCitrine cells transferred to M04G-02 gradient plates as described above were incubated for at least 3 hr, followed by serum starvation for at least 6 hr before imaging. Existing cell media was substituted right before imaging with imaging media. Confocal imaging for multiple positions at 1min time interval using adaptive auto-focus system and the water objective was performed concurrently during the duration of the experiment using the Leica TCS SP8i.

For migration experiments under uniform EGF647 stimulation, confocal laser scanning microscopy / transmission imaging of live MCF7EGFRmCitrine / MCF10A cells was done on a Leica TCS SP8i or Olympus IX81 for multiple positions at 3 min and 2 min time interval respectively, using the 10 × dry objective for 14 hr.

EGF647 / Fluorescein gradient quantification

Request a detailed protocol

hEGF647 was generated in the lab of Prof. P. Bastiaens, MPI of molecular Physiology, Dortmund, using the His-CBD-Intein-(Cys)-hEGF-(Cys) plasmid (Sonntag et al., 2014), kindly provided by Prof. Luc Brunsveld, University of Technology, Eindhoven. Human EGF was purified from E. coli BL21 (DE3), N-terminally labeled with Alexa647-maleimide as described previously (Sonntag et al., 2014) and stored in PBS at –20 °C. To quantify the spatial extent of the EGF647 /Fluorescein gradient, gradients were generated following the protocol described in sub-section 5.6 in plates without cells or matrix coating. Confocal images of Alexa647 /GFP channel were acquired at 1min interval. A rectangular region of interest (including the perfusion channels and the culture chamber) was used to obtain an averaged pixel intensity profile using FIJI at each time point. This spatial profile was averaged across multiple experiments and then scaled with the mean intensity value in the perfusion channel, which corresponds to the applied EGF647 /Fluorescein concentration.

Quantifying EGFRmCitrine phosphorylation in single cells

Request a detailed protocol

To quantify plasma membrane EGFRmCitrine phosphorylation in live MCF7EGFRmCitrine cells, single-cell masks were obtained from the EGFRmCitrine channel at each time-point using FIJI (https://imagej.net/Fiji). All pixels within the obtained boundary were radially divided into two segments of equal areas (Stanoev et al., 2018), and the outer segment was taken to represent the plasma membrane. For the kymograph analysis, at each time point, the plasma membrane segment was divided into 4 quadrants in anti-clockwise direction, and each was divided into 5 spatial bins (Figure 2A). The fraction of phosphorylated EGFRmCitrine in each bin, i was estimated as:

(1) EGFRpi(t)=PTBPMi(t)/(PTBT(t)-PTBendo(t))EGFRPMi(t)/EGFRT(t)

where PTBPMi(t) and EGFRPMi(t) are respectively the PTBmCherry and EGFRmCitrine fluorescence at ith plasma membrane bin, PTBT(t) and EGFRT(t) - respective total fluorescence in the whole cell, PTBendo(t) – the PTBmCherry fluorescence on vesicular structures in the cytoplasm. Endosomal structures were identified from the cytosol by intensity thresholding (1.5 s.d. percentile) and PTBmCherry fluorescence from these structures was subtracted from the PTBT(t), to correct for the PTBmCherry fraction bound to the phosphorylated EGFRmCitrine on endosomes.

Temporal profile of the fraction of phosphorylated EGFRmCitrine on the plasma membrane was obtained using:

(2) EGFRp(t)=i=120PTBPMi(t)(PTBT(t)-PTBendo(t))i=120EGFRPMi(t)(EGFRT(t))

and then normalized as:

(3) EGFRp(t)=EGFRp(t)<EGFRp>t[0,5min]maxt(EGFRp(t))<EGFRp>t[0,5min]

with <> being the temporal average in the pre-stimulation interval t[0,5min]. The fraction of liganded receptor was calculated using:

(4) EGF-EGFR(t)=EGFPMEGFRPM(t)

To classify single cells into non-activated, activated (polarized EGFRmCitrine phosphorylation) and pre-activated (uniformly distributed EGFRmCitrine phosphorylation) upon gradient EGF647 stimulation (Figure 2—figure supplement 2A, B), the following method was applied. To identify pre-activated cells, a Gaussian Mixture Model (GMM) was fitted to the histogram of (EGFRpi)t[0,5min] values from all the analysed cells, and the intersection point between the two normal distributions was identified. If more than 30% of the (EGFRpi)t[0,5min] pixel intensity values for any cell lie above the intersection point, the cell is classified as pre-activated. To distinguish between the non-activated and activated cells in the remaining population, average EGFRmCitrine phosphorylation value (EGFRp) per cell was estimated during the pre-stimulation (t[0,5min]) and the stimulation period (t[5min,65min]) (<EGFRp>t[0,65]) from the temporal EGFRmCitrine phosphorylation profiles. Histogram of the respective EGFRp values was again fitted with a GMM model. All cells with an average <EGFRp>t[0,65] value lying below the intersection point were considered to be non-activated, whereas those above - activated.

The average of the spatial projection of the fraction of phosphorylated EGFRmCitrine from single-cell kymographs (Figure 2—figure supplement 1C) was generated from the 20 cells that were polarized in the direction of the EGF647 gradient. For each cell, a temporal average of EGFRp per bin was calculated for the duration of the gradient (t[5min,65min]) and the bin with the maximal EGFRp value was translated to π. The profiles were then smoothened using a rolling average with a window of 7 bins. The resulting profiles were then averaged over all cells and mean ± s.d. is shown.

The local spatial EGF647 distribution around single cells (Figure 2—figure supplement 1F) was estimated as follows: the cell mask obtained using the EGFRmCitrine images were dilated outwards by 8 pixels to account for possible ruffles, and then by additional 15 pixels. The secondary rim of 15 pixels around the cell mask was used to calculate the spatial distribution of EGF647 outside single cells. This outer contour was divided in 20 bins as for the kymographs, and EGF647 intensity was quantified in each bin. The angle between the direction of EGF647 and the direction of EGFR phosphorylation was calculated as the amount of radial bins between the maxima in the spatial projections. This bin-distance was then translated into an angle under the assumption of a circular perimeter.

In order to identify the characteristic features of the EGFRmCitrine phosphorylation profile during the transition from polarized to unpolarized state, the single-cell EGFRp(t) profiles with and without Lapatinib treatment after gradient wash-out were fitted to an inverse sigmoid function given by,

(5) f(t)=a0an+tn

were a0, a are constants and n is the Hill-coefficient (examples in Figure 2—figure supplement 2E, F). Non-linear least square method (python package curve fit) was used to perform the fitting. Under normal conditions (w/o Lapatinib), a10, a0103 and n2.88 fitted well the data (R20.79). The same function however could not describe the EGFRp profiles in the Lapatinib treatment experiment (median R20.33). The Lapatinib treatment profiles were therefore fitted by fixing a=10, and leaving a0 and n as free parameters, as they determine the upper plateau and the steepness of the drop to the basal level. In this case, a019 and n1.28 were identified from the fitting (median R20.84, Figure 2—figure supplement 2E, G). From the fitted profiles in both cases, half-life was estimated to be the time frame in which 50% of EGFRmCitrine phosphorylation is lost after EGF647 removal.

Estimating memory duration in EGFRmCitrine phosphorylation polarization

Request a detailed protocol

The duration of memory in EGFRmCitrine phosphorylation polarization in single cells was estimated from the temporal profile of the fraction of plasma membrane area with high EGFRmCitrine phosphorylation during and after gradient removal (Figure 2D and E). For this, the single-cell kymographs were normalized to a maximal value of 1 using

(6) EGFRpi(t)=EGFRpi(t)<EGFRp>t[0,5min]maxt(EGFRp(t))<EGFRp>t[0,5min]

yielding the value of phosphorylated EGFRmCitrine per bin i per time point t. Using the mean of EGFRp+s.d. over the whole experiment duration as a threshold, all EGFRpi(t) lying above the threshold were taken to constitute the area of polarized EGFRmCitrine phosphorylation. To account for different bin sizes, at each timepoint, the area of all bins with EGFRp above the threshold was summed and divided by the respective total cell area, yielding the temporal evolution of the fraction of polarized cell area (FPA) (Figure 2D). The end of the memory duration per cell was identified as the time point at which FPApercell<(FPAaverages.d.) in 3 consecutive time points (Figure 2E).

Quantifying morphological changes in response to EGF647 in experiments and simulations

Request a detailed protocol

Morphological changes of polarized cells were quantified using the solidity (Figure 2H1) of each cell at each time point and the directed protrusive area towards and away from the gradient (Figure 1G and H; Figure 2—figure supplement 2H). The solidity σ is the ratio between the cell’s area Acell and the area of the convex hull Aconvex (σ=AcellAconvex). The memory duration in cell morphology was calculated from the single-cell solidity profiles, and corresponds to the time-point at which the solidity is below mean-s.d. estimated during gradient presence. The directed cell protrusion area was estimated by comparing single cell masks at two consecutive time points. To reduce noise effects, the masks were first subjected to a 2D Gaussian filtering using the filters.gaussian function from the scipy python package. Protrusions were considered if the area change was greater than 10 pixels or 1.2 µm2 per time point. The front and the back of the cells were determined by identifying an axis that runs perpendicular to the gradient and through the cell nucleus of the initial time point. The directed cell protrusion area was then obtained using Aprot,frontAfront-Aprot,backAback. The final profiles of directed protrusive area were smoothed using 1D Gaussian filtering with the filters.gaussian_filter1d function from the scipy python package. For the equivalent quantification from the simulations, the same procedures were applied without an area threshold. The memory duration was estimated as the time point at which the directed protrusive area crosses zero after the gradient removal.

Quantification of single-cell migration and duration of memory in directed cell migration

Request a detailed protocol

Single cell migration trajectories were extracted using Trackmate (Tinevez et al., 2017) in Fiji (Schindelin et al., 2012) using Hoechst 33342/transmission channel. From the positional information (x and y coordinates) of individual cell tracks, quantities such as Motility, Directionality and cosθ were extracted using custom made Python code (Python Software Foundation, versions 3.7.3, https://www.python.org/). Directionality was calculated as displacement over total distance and statistical significance was tested using two-sided Welch’s t-test. To quantify the memory duration in directed single-cell migration, the Kernel Density Estimate (KDE) from cosθ quantification in the continuous absence of EGF647 (uniform case, between 250–300 min) was compared with a moving window KDE (size of 5 time points) from the gradient migration profile, using two sided Kolmogorov-Smirnov test. To verify the absence of memory when cells were treated with Lapatinib during gradient wash-out, a moving window KDE (5 time points) from cosθ obtained in this case was compared to the KDE in continuous absence of EGF647 (uniform case Figure 3—figure supplement 2B, between 250–300 min) using two sided Kolmogorov-Smirnov test (Figure 3I). Furthermore, the KDE between 300–350 min and 350–840 min (after gradient removal) was statistically equivalent to the KDE in continuous absence of EGF647, confirming the rapid switch from directed to random-walk migration in the Lapatinib case (Figure 3—figure supplement 2H). To estimate the time required for complete reversal of cell migration direction when the cells were subjected to a gradient from opposite direction, KDE distributions were compared between the following time windows: 275–335 min (second gradient), 335–365 min, 365–385 min, 375–385 min, and 365–605 min (third gradient).

To quantify the motility patterns of MCF10A cells in absence, uniform or gradient EGF647 stimulation, we fitted the experimentally obtained single cell migration trajectories using modified Ornstein-Uhlenbeck process (mOU) (Uhlenbeck and Ornstein, 1930) that is defined by the Langevin equation for the velocity vector ν:

(7) dν(t)dt=-1τν(t)+2Dτ(ξ(t)+b(t))

where ξ(t) represents a white noise component, D is a diffusion coefficient characteristic of a Brownian motion, τ is the persistence time and b(t) models the contribution of the time-dependent bias. The experimental data was fitted to obtain values of D and τ. In order to estimate D, Mean Square Displacement (MSD) was calculated from the single cell tracks using MSD(t)=<|xi(t)xi(0)|2>, where xi(t) is the tracked position of i-th cell in the 2D plane, <> is the average across all single cell tracks, and |.| is the Euclidean distance (Selmeczi et al., 2005). To estimate D, the obtained MSD profile was fitted with a linear function (=4Dt). Goodness of Fit for the different experimental conditions: 0 ng/ml EGF647, R2=0.975; for uniform 20 ng/ml EGF647 stimulation, R2=0.995. In order to estimate τ, Velocity Auto-Correlation Function VACF(t)=<νi(t)νi(0)>, where νi(t) is the measured velocity of i-th cell at time t, was fitted with a mono exponential function (=ϕ0e-tτ). Goodness of Fit: for 0 ng/ml EGF647 case - Standard Error of Estimate SEOE=0.0261; for uniform 20 ng/ml EGF647 stimulation case, SEOE=0.0570. Fitted values: for 0 ng/ml EGF647 case, τ=11.105, D=0.425; for uniform 20 ng/ml EGF647 stimulation case, τ=38.143, D=2.207; bias b(t)=0.134.

To compute the duration of memory in directional migration after gradient removal for individual cells (Figure 3—figure supplement 3), single cell migration tracks were first smoothened using a Kalman-filter (python package filterpy.kalman) by predicting the cell position and velocity. The cell’s displacement angles relative to the gradient direction (cosθ) were calculated for each cell at each timepoint, rendering single-cell cosθ plots (Figure 3—figure supplement 3B,C). The memory duration was then calculated as the point where three consecutive timepoints in the cosθ profiles fall below a threshold cosθ value of 0.75.

Reconstructing state-space trajectories from temporal EGFRmCitrine phosphorylation profiles

Request a detailed protocol

The state-space reconstruction in Figure 2F and G was performed using the method of time-delay. For a time series of a scalar variable, a vector x(ti), i=1,...N in state-space in time ti can be constructed as following

(8) X(ti)=[x(ti),x(ti+d),..,x(ti+(m-1)d)]

where i=1 to N(m1)d, d is the embedding delay, m - is a dimension of reconstructed space (embedding dimension). Following the embedding theorems by Takens (Takens, 1980; Sauer et al., 1991), if the sequence X(ti) consists of scalar measurements of the state of a dynamical system, then under certain genericity assumptions, the time delay embedding provides a one-to-one image of the original set, provided m is large enough. The embedding delay was identified using the timeLag function (based on autocorrelation), the embedding dimension using the estimateEmbeddingDims function (based on the nearest-neighbours method), and the state-space reconstruction using the buildTakens function, all from the nonlinearTseries package in R (https://cran.r-project.org/web/packages/nonlinearTseries/index.html). Before state-space reconstructions, time series were smoothened using the Savitzky-Golay filter function in Python. For Figure 2F, d=26, de=3; for Figure 2G, d=50, de=3.

Theoretical consideration of the navigation mechanism in a generalized reaction-diffusion signaling model

Request a detailed protocol

We consider a generalized form of a (mass-conserved) reaction-diffusion (RD) model of an M (URM) component system in N (xRN) dimensional space

(9) U(x,t)t=F(U(x,t))+D2U(x,t)

where FRM is the reaction term, D is a M×M diagonal matrix of diffusion constants Dj,j=1,,M, and 2 is the Laplacian operator. Standard analysis of such models relies on linear stability analysis to find the conditions for a Turing-type instability (Turing, 1952), such that the symmetric steady state becomes unstable and an asymmetric polarized state is stabilized. By its nature, the linear stability analysis makes no prediction about the transition process itself, and thereby the type of bifurcation that underlies it. To provide quantitative description of the symmetry breaking transition in reaction-diffusion models, local perturbation analysis can be applied (Holmes et al., 2015). However, this analysis is mainly restricted to models characterized with large diffusion discrepancy between the signaling components. The conditions for a pitchfork bifurcation (PB)-induced transition in a generic RD model therefore have to be formally defined. Let Us=(uis) for i=1,...,M, be the stable homogeneous symmetric steady state of the RD system. Consider a linear perturbation of the form

(10) U(x,t)=Us+δU(x)e(λt),     δU(x)RM

where δU(x) is the spatial and e(λt) is the temporal part of the perturbation. Substituting Equation 10 in Equation 9 yields a linearized eigenvalue equation whose solution can be determined by solving the characteristic equation, Fλ=det(λIM×MJM×M)=0. J is the Jacobian matrix of the system defined by Jij=Fi(U(x,t))Uj,i=1,.,M,j=1,..,M.

The system exhibits a PB if, an odd eigenfunction δU(x) such that δU(-x)=-δU(x), taken in the limit λ0, fulfills the following condition (Paquin-Lefebvre et al., 2020):

(11) limλ0Fλ=det(J)=0.

When this conditions is satisfied, the symmetric, homogeneous steady state of the system undergoes a pitchfork bifurcation and an inhomogeneous steady state (IHSS) with two branches of asymmetric steady states emerges. In terms of polarization, these branches correspond to front-back-polarized states, where the orientation depends on the direction of the external signal (Figure 1A, Figure 1—figure supplement 1A).

To identify whether the PB is of sub-critical type, and thereby identify the presence of an SNPB, a weakly nonlinear analysis of Equation 9 must be performed to obtain description of the amplitude dynamics of the inhomogeneous state. This can be achieved using an approximate analytical description of the perturbation dynamics based on the Galerkin method (Becherer et al., 2009; Rubinstein et al., 2012; Bozzini et al., 2015). For simplicity, we outline the steps for a one-dimensional system (N=1). As we are interested in the description of a structure of finite spatial size (i.e. finite wavelength k), the final solution of the PDE is expanded around the fastest growing mode, km into a superposition of spatially periodic waves. That means that u(x,t)U can be written as:

(12) u(x,t)n=-+(un(t)enikmx+un*(t)e-nikmx)

where un(t) is the complex amplitude of the nth harmonics. Let the amplitude corresponding to the leading harmonics (n=1) is ϕ(t). After assuming that the amplitude of every other harmonics can be written as a power series of ϕ(t), substituting Equation 12 into Equation 9 allows to write an equation that describes the evolution of ϕ(t). In the case when the resulting equation is of Stuart-Landau type:

(13) dϕdt=c1ϕ+c2ϕ3-c3ϕ5

with c1,c2,c3>0, this corresponds to the normal form of a sub-critical pitchfork bifurcation (Strogatz, 2018). Together with the condition given by Equation 11, the existence of a sub-critical PB for the full system (Equation 9) is guaranteed. A numerical or analytical analysis of Equation 13 enables the identification of the position of the SNPB.

Modeling EGFR phosphorylation polarization dynamics

Request a detailed protocol

The dynamics of the experimentally identified spatially distributed EGFR sensing network (Figure 1B, Figure 1—figure supplement 1B) is described using the following one-dimensional system of partial differential equations (PDEs):

(14) [Ep]t= f1([Ep],[EEp],[RGa],[N2a],[EGFt])+DEp2[Ep]x2[EEp]t= f2([Ep],[EEp],[EGFt])+DEEp2[EEp]x2[RGa]t= f3([Ep],[EEp],[RGa])+DRGa2[RGa]x2[N2a]t= f4([Ep],[EEp],[N2a])

with

 f1=([Et][Ep][EEp])(α1([Et][Ep][EEp])+α2[Ep]+α3[EEp])γ1[RGa][Ep]γ2[N2a][Ep]kon([EGFt][EEp])[Ep]2+1/2koff[EEp];f2=kon([EGFt][EEp])([Ep]2+([Et][Ep][EEp])2)koff[EEp];f3=k1([RGt][RGa])k2[RGa]β1[RGa]([Ep]+[EEp]);

and

f4=ϵ(k1([N2t]-[N2a])-k2[N2a]+β2([Ep]+[E-Ep])([N2t]-[N2a])).

The reaction terms are described in details in Stanoev et al., 2018. In brief, [E-Ep] is the phosphorylated ligand-bound dimeric EGFR, [Ep] - ligandless phosphorylated EGFR, [Et] - total amount of EGFR, [RGa],[RGt] and [N2a],[N2t] - the active and total amount of the membrane localized PTPRG and the ER-bound PTPN2, respectively. Both, the receptor and the deactivating enzymes have active and inactive states, and the model equations describe their state transition rates. Therefore, mass is conserved in the system and the total protein concentrations of the three species ([Et], [RGt] and [N2t]) are constant parameters. Autonomous, autocatalytic and ligand-bound-induced activation of ligandless EGFR ensue from bimolecular interactions with distinct rate constants α1-3, respectively. Other parameters are as follows: k1/k2 — activation/inactivation rate constants of the phosphatases, β1/β2 - receptor-induced regulation rate constants of PTPRG/PTPN2, γ1/γ2 - specific reactivity of the enzymes (PTPRG/PTPN2) towards the receptor. The EGFR-PTPN2 negative feedback is on a time scale (ϵ) approximately two orders of magnitude slower than the phosphorylation-dephosphorylation reaction, as estimated from the ~ 4 min recycling time of EGFRp (Stanoev et al., 2018). This enables, when necessary, to consider a quasi-steady state approximation for the dynamics of PTPN2 for simplicity:

(15) [N2a]qss=[N2t](k1+β2([Ep]+[E-Ep]))k1+k2+β2([Ep]+[E-Ep])

[EGFt] denotes the total ligand concentration. Assuming that at low, physiologically relevant EGF doses, the ligand will be depleted from the solution due to binding to EGFR (Lauffenburger and Linderman, 1996), ligand-binding unbinding was explicitly modeled (kon, koff) in Equation 14, with values corresponding to the experimentally identified ones.

The diffusion terms model the lateral diffusion of the EGFR and PTPRG molecules on the plasma membrane, whereas PTPN2 is ER-bound and does not diffuse. Single particle tracking studies have demonstrated that EGFR molecules on the plasma membrane occupy three distinct mobility states, free, confined and immobile, with the occupations of the free and immobile states decreasing and increasing significantly after EGF stimulation, respectively (2 min after EGF stimulation, corresponding with the time-scale of EGF binding) (Ibach et al., 2015). In the reaction-diffusion (RD) simulations therefore for simplicity, it is assumed that DE-Ep0, whereas diffusion constants of same order are assumed for the ligandless EGFR and PTPRG (DEpDRGa).

Analytical consideration for an SNPB existence in the EGFR network

Request a detailed protocol

To identify analytically the existence of a SNPB in the EGFR receptor network, we performed a weakly nonlinear analysis as described in the general consideration (Section. Theoretical consideration of the navigation mechanism in a generalized reaction-diffusion signaling model). For this, we considered the system Equation 14, where the dynamics of PTPN2 is at quasi-steady state (Equation 15), [E-Ep]=0, and rest of the dependent and independent variables were scaled to have a dimensionless form. Let [Ep~]=[Ep]/E0, [RGa~]=[RGa]/RG0, x~=x/x0, τ=t/t0, such that t0=1/(k1+k2), E0=k1/β2, RG0=(k1+k2)/γ1 and t0/x02=1/DEp. Substituting these into Equation 14 yields the system of dimensionless equations:

(16) [Ep~]τ=q1+q2[Ep~]+q3[Ep~]2-[RGa~][Ep~]-q4(1+[Ep~])[Ep~](1+k+[Ep~])+2[Ep~]x~2[RGa~]τ=r1-[RGa~]-r2[RGa~][Ep~]+D2[RGa~]x~2

with q1=α1[Et]2(k1+k2)β2, q2=(α22α1)[Et]k1+k2, q3=(α1α2)kt(k1+k2)β2, q4=γ2[N2t]k1+k2, k=k2/k1, r1=k1[RGt]γ1(k1+k2)2, r2=β1k1(k1+k2)β2 and D=DRGαDEp.

We further simplify the system Equation 16 by taking the Talyor series expansion of the quasi-steady state approximation of [N2a] around Es, the steady state of [Ep]~

(17) q4(1+[Ep~])[Ep~]1+k+[Ep~]= q7+q8[Ep~]+q9[Ep~]2+o([Ep~]2)

with q7=Esq41+k+EsEsq4(1+k)(1+k+Es)2, q8=Esq41+k+Es+q4(1+k)(1+k+Es)2(1Es), and q9=q4(1+k)(1+k+Es)2, thus yielding:

(18) [E~p]τ=q9+q10[E~p]+q11[E~p]2[RG~a][E~p]+2[E~p]x~2[RG~a]τ=r1[RG~a]r2[RG~a][E~p]+D2[RG~a]x~2

with q9=q1q7,q10=q2q8 and q11=q3q9.

To avoid long expression in the further analysis, we re-name the dependent variables as u1=[Ep~] and u2=[RGa~], and the independent variables as x~=x, τ=t. The system Equation 16 therefore obtains the generic form:

(19) u1t=F1(u1,u2)+2u1x2u2t=F2(u1,u2)+D2u2x2.

In order to perform linear stability analysis, a one-dimensional projection of Equation 19 is considered,

(20) du1fdt=F1(u1f,u2f)-(u1f-u1b)=G1(u1f,u2f,u1b)du2fdt=F2(u1f,u2f)-D(u2f-u2b)=G2(u1f,u2f,u2b)du1bdt=F1(u1b,u2b)-(u1b-u1f)=G3(u1b,u2b,u1f)du2bdt=F2(u1b,u2b)-D(u2b-u2f)=G4(u1b,u2b,u2f)

The simplified one-dimensional geometry assumes a model composed of two compartments (front and back), resembling a projection of the membrane along the main diagonal of the cell. The standard approach of modeling the diffusion along the membrane in this case is simple exchange of the diffusing components. The one-dimensional projection, as demonstrated below, preserves all of the main features of the PDE model.

Let, Us=(u1fsu2fsu1bsu2bs) be the stable symmetric steady state of the system (u1fs=u1bs, u2fs=u2bs). A small amplitude perturbation on this symmetric steady state of the form,

(21) (u1f(t)u2f(t)u1b(t)u2b(t))=(u1fsu2fsu1bsu2bs)+(δu1fδu2fδu1bδu2b)eλt

yields a linearized equation,

(22) λ(dδu1fdtdδu2fdtdδu1bdtdδu2bdt)=J(δu1fδu2fδu1bδu2b)

where

J=(G1u1fG1u2fG1u1b0G2u1fG2u2f0G2u2bG3u1f0G3u1bG3u2b0G4u2fG4u1bG4u2b)

is the Jacobian of the system evaluated at the symmetric steady state. In order to identify existence of PB in the system, the condition given in Equation 11 should be satisfied for an odd mode of the perturbation. For the one-dimensional projection (Equation 20), the odd mode of the perturbation (δU(x)=δU(x)) must yield: δu1f=-δu1b and δu2f=-δu2b. Substituting this into Equation 22 to obtain F-(λ), in the limit λ0 renders:

(23) limλ0F(λ)=det((G1u1f+G3u1b)(G1u1b+G3u1f)(G1u2f+G2u2b)(G2u1f+G4u1b)(G2u2f+G4u2b)(G2u2b+G4u2f))=0

Thus, there exists parameter set for which existence of PB in the system Equation 20 is guaranteed.

To identify whether the PB is sub-critical and thereby identify existence of a SNPB, the solution of the system Equation 19 is approximated as in Equation 12:

(24) u(x,t)=ϕ(t)eikmx+ϕ*(t)e-ikmx+u0(t)+n=23(un(t)enikmx+un*(t)e-nikmx)v(x,t)=ϕ(t)eikmx+ϕ*(t)e-ikmx+v0(t)+n=23(vn(t)enikmx+vn*(t)e-nikmx)

The expansion is taken to n=3rd order, rendering an amplitude equation of 5th order. As described in Becherer et al., 2009, the complex coefficients of the n=0th,n=2nd and n=3rd harmonics can be approximated as power series of ϕ(t). Substituting into Equation 19 allows to derive these coefficients. This yields a system of coupled ODEs representing the time evolution of the complex amplitudes, in this case, for ϕ(t), u0(t), v0(t), u1(t), v1(t), u2(t), v2(t), u3(t) and v3(t). Assuming that the dynamics of the higher order harmonics reaches their steady state much faster than the leading perturbation does, the derivatives of their amplitudes can be set to zero. This allows to obtain expressions of the amplitudes purely as functions ϕ and the parameters of the system as:

(25) u0(ϕ)=(1q10(2(1-q11)-q9|ϕ|2))|ϕ|2v0(ϕ)=(r1|ϕ|2-2r2)|ϕ|2u2(ϕ)=u2(2)ϕ2v2(ϕ)=v2(2)ϕ2u3(ϕ)=u3(3)ϕ3v3(ϕ)=v3(3)ϕ3

where u2(2)=1-q11q10-4km2, v2(2)=-r21+4Dkm2, u3(3)=u2(2)+v2(2)-2q11u2(2)q10-9km2 and v3(3)=-r2(u2(2)+v2(2))1+9Dkm2. The dynamics of the leading harmonics (n=1) can be written as:

(26) dϕdt=c1ϕ+c2ϕ3-c3ϕ5

where c1=q10-km2-r1+q9(1-2q11)q10, c2=(1-q11)(2q11-1)(2q10-1q10-4km2)+r2(2+11+4Dkm2) and c3=2q11u2(2)u3(3)-u2(2)v3(3). Equation 26 is of Stuart-Landau type and represents a normal form of a sub-critical pitchfork bifurcation. This shows the existence of SNPB in the EGFR network.

To corroborate this, we also performed numerical bifurcation analysis on one-dimensional projection (Equation 20) where the reaction terms have the form as defined in Equation 14, including the full form for [N2a], when [E-Ep]=0. The bifurcation analysis (Figure 1—figure supplement 1C) was obtained using the Xppaut software package (Ermentrout, 2016). The parameters in the model Equation 14 have been described in Stanoev et al., 2018, where they were calibrated with experimental data: α1=0.001, α2=0.3, α3=0.7, β1=11, β2=1.1, k1=0.5, k2=0.5, g1=1.9, g2=0.1, kon=0.05, koff=0.28, ϵ=0.01, RGt=1, N2t=1; and the diffusion-like terms have been scaled from the values derived in Orr et al., 2005: DEp~=0.02, DRGa~=0.02 (see also Supplementary file 1).

The bifurcation analysis is performed with respect to total EGFR concentration at the plasma membrane in order to reveal all possible dynamical regimes of the system. This analysis demonstrates that for the spatially distributed EGFR network, the homogeneous steady state (HSS, gray solid line, Figure 1—figure supplement 1C) representing basal non-polarized state losses stability via a symmetry-breaking pitchfork bifurcation (PB), which gives rise to a polarized state represented via an inhomogeneous steady states (IHSS). The polarized state is stabilised via saddle-node bifurcations (SNPB) (Figure 1—figure supplement 1C, magenta branched lines). There is a coexistence between the HSS and the IHSS before the PB, rendering it sub-critical. The IHSS (Koseska et al., 2013) that gives rise to the stable polarized state is a single attractor that describes a heterogeneous state with two branches corresponding to orientation of the front-back-polarized state. The IHSS solution is therefore fundamentally distinct from a bistable system where the high and the low phosphorylation states correspond to two different homogeneous steady states. As the IHSS is a single attractor, the high and low phopshorylation state are interdependent, rendering the PB a unique mechanism for generating robust front-back polarization.

We next describe the dynamical basis of the polarization and memory of polarization in details. We assume that the steady state EGFR concentration at the plasma membrane corresponds to organization at criticality, before the SNPB. For this receptor concentration, only the basal unpolarized state (HSS) is stable (Figure 1—figure supplement 1A, top left, schematic representation). In the presence of a spatially inhomogeneous EGF signal however, the system undergoes a series of complex transitions through which the topology of the phase space changes. In particular, the inhomogeneity introduced by the localized signal leads to unfolding of the pitchfork bifurcation, such that for the same organization (the given EGFR concentration), only the polarized state (the IHSS) is stable (Figure 1—figure supplement 1A, top right). This unfolding of the PB therefore enables robust transition from basal to polarized state. When the EGF signal is removed, the system undergoes again topological phase space changes. However, in this transition, the system does not revert back to the unpolarized state immediately, but rather it is transiently maintained in the ”ghost” of the SNPB that is lost in this transition (Figure 1—figure supplement 1A, low). This is manifested as a transient memory of the polarized state, after which the system rapidly reverts to the basal state.

The reaction diffusion simulations were performed by assuming PTPN2 at quasi-steady state. The cell boundary was represented with a 1D circular domain of length L=2πR (where R=2μm) which was then divided into 20 equal bins. The diffusion terms were approximated by central difference method, enabling for conversion of the PDE system to a system of ordinary differential equations (ODEs). Stochastic simulations with additive white noise were implemented by adding σdWt (σ=0.02, dWt is sampled from a normal distribution with mean 0 and variance 0.01) in the equation for [Ep]. The stochastic sdeint Python package was used. Parameters: DEp=DRGa=0.008μm2/min. DEp was taken from Orr et al., 2005 and scaled to correspond to a cell with perimeter L in the simulations. For organization in the homogenous symmetric steady states (the basal and pre-activated states), organization at criticality or in the stable polarized state (IHSS), Et{1.1,1.85,1.26,1.35} respectively, time step was set to 0.01 min, other parameters as above. Periodic boundary conditions were used. To mimic the dynamic nature of EGF647 gradient, a Gaussian function on a periodic window with varying amplitude and standard deviation was used (shape shown in Figure 1D, top). To represent the state-space trajectory (Figure 1F, bottom), stochastic realization of the one-dimensional projection of the full system (as for the bifurcation analysis) was used.

Physical model of single-cell chemotaxis

Request a detailed protocol

To describe signal-induced cell shape changes and subsequent cell migration, we combined the dynamical description of the gradient sensing capability of the EGFR network (Equation 14, Figure 1B) together with a physical model for cellular migration, thereby implicitly modeling the signal-induced cell shape changes (Figure 1C). In order to couple a mechanical model of the cell with the biochemical EGFR signaling model as a means to simulate large cellular deformations, we utilized the Level Set Method (LSM) (Osher and Sethian, 1988) as described in Yang et al., 2008. Briefly, the cell boundary at time t is described on a two-dimensional Cartesian grid by the closed-contour Γ(t)={x|Ψ(x,t)=0}, that represent the zero-level set of the potential function Ψ(x,t), taken to have an initial form:

(27) Ψ(x,0)={d(x,Γ),if xSd(x,Γ),if xS0,if xΓ

where S identifies the area occupied by the cell and d(x,Γ) is the distance of position x to the curve Γ. Thus, the cell membrane is represented implicitly through the potential function which is defined on the fixed Cartesian grid, eliminating the need to parameterize the boundary, and thereby enabling to handle complex cell boundary geometries.

The shape of the cell (Γ(x,t)) evolves according to the Hamilton-Jacobi equation:

(28) Ψ(x,t)t+v(x,t)Ψ(x,t)=0

The vector v(x,t) is the velocity of the level set moving in the outward direction, thereby intrinsically describing the cell’s membrane protrusion and retraction velocities that are driven by internally generated mechanical forces (e.g. actin polymerization or myosin-II retraction, Bray, 2000). To determine how these forces translate to membrane velocity, a mechanical model that describes the viscoelastic behavior of the cell represented as a viscoelastic cortex surrounding a viscous core, is implemented. Following Yang et al., 2008, the cortex connecting the cell membrane and the cytoplasm is represented by a Voigt model parallel connection of an elastic element kc and a viscous element τc, whereas the cytoplasm is modeled as a purely viscous element, τa, which is placed in series with the Voigt model.

Let l(x,t), xΓ(t) be the viscoelastic state of the cell at time t and at a position x on the membrane, such that |l| represents the length of the numerous parallel unconnected spring-damper systems. The viscoelastic state of the cell then evolves according to:

(29) -kcτcl(t)+1τcPtotal(t)=lv(t)+l(t)t

where ∇ is the gradient operator, the pressure Ptotal(t)=Pprot(t)+Pretr(t)+Parea(t)-Pten(t) is sum of the protrusion, retraction, area conservation, and cortical tension pressures, respectively. The EGFR signaling state ([Ep]) directly determines the protrusion/retraction pressure, since high/low signaling activity triggers actin polymerization / myosin-II retraction following: Pprot(t)= Kprot(([Ep](t)<[Ep](t)>)/([Ep]max(t)<[EP](t)>))n and Pret(t)= Kretr((<[Ep]>[Ep])/(<[EP]>[Ep]max))n, where <.> denotes mean at the membrane, Kprot, Kretr - proportionality constants. The cell is assumed to be flat with uniform thickness, such that the 2D area (A(t)) of the cell is conserved (Parea(t)=Karea(A(0)-A(t))n), Karea - proportionality constant. The pressure generated by the cortical tension therefore depends only on the 2D local surface curvature and the 2D equilibrium pressure, rendering the rounding pressure due to cortical tension to be Pten(t)=Kten(κ(Γ)-1/R)n, with κ(x) being the local membrane curvature, R - initial cell radius, was set to 2 μm, and Kten - proportionality constant. The local membrane velocity v(x), xΓ(t) depends both on the viscoelastic nature of the cell and on the effective pressure profile (Ptotal(t)) and is given by,

(30) v=kcτcl+(1τc+1τa)Ptotal

For the simulations in Figures 1 and 4 and Figure 4—figure supplements 1 and 2 first the stochastic PDEs (Equation 14) are solved and the kymographs of the signalling ([Ep]) activity are generated. The viscoelastic state is initialized with zero value on the membrane, l(x,0)=0. At each time point, Ptotal is estimated, as well as the local membrane velocity using Equation 30. This velocity is then used to evolve both the viscoelastic state (Equation 29) and the potential function (Equation 27).

The spatial discretization of these advection equations (Equations 28 and 29) was performed using the upwindENO2 scheme, as described in the Level Set Toolbox (Mitchell, 2007) and was integrated with first-order forward Euler method. The time step was set to 0.01min and the potential function was solved on a 2D Cartesian grid with spatial discretization of 5 points per µm. All the codes were custom implemented in Python. Parameters: kc=0.1nN/μm3, τc=0.08nNmin/μm3, τa=0.1nNmin/μm3, Kprot=0.08nN/μm2, Kretr=0.05nN/μm2, Karea=0.02nN/μm4, Kten=0.1nN/μm. Kten was taken from the literature, corresponding to an experimentally measured range of cell cortical tension values (Cartagena-Rivera et al., 2016). The rest of the parameters were selected to match the cell migration speed during gradient and memory phase, estimated from the experiments (Figure 3A, v=0.49±0.173μm/min).

Data availability

Source data is provided with the submission. The numerical data used to generate the corresponding figures can be obtained from the codes deposited in https://github.com/akhileshpnn/Cell-memory, (copy archived at swh:1:rev:288921244e5042922e1bbddcf5037a5e87e78723).

The following data sets were generated
    1. Akhilesh N
    (2022) GitHub
    ID Cell-memory. Cells use molecular working memory to navigate in changing chemoattractant fields.

References

  1. Book
    1. Ermentrout B
    (2016)
    XPPAUT
    pitt.
  2. Book
    1. Lauffenburger D
    2. Linderman J
    (1996)
    Receptors: Model for Binding, Trafficking and Signaling
    Oxford University Press.
    1. Sauer T
    2. Yorke JA
    3. Casdagli M
    (1991) Embedology
    Journal of Statistical Physics 65:579–616.
    https://doi.org/10.1007/BF01053745
  3. Book
    1. Takens F
    (1980) Detecting strange attractors in turbulence
    In: Rand D, Young LS, editors. Dynamical Systems and Turbulence, Warwick 1980. Berlin, Heidelberg: Springer Nature. pp. 366–381.
    https://doi.org/10.1007/BFb0091924
    1. Turing A
    (1952) The chemical basis of morphogenesis
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237:37–72.
    https://doi.org/10.1098/rstb.1952.0012

Decision letter

  1. Arvind Murugan
    Reviewing Editor; University of Chicago, United States
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel
  3. Elizabeth R Jerison
    Reviewer; Stanford University, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Cells use molecular working memory to navigate in changing chemoattractant fields" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Elizabeth R Jerison (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The reviewers broadly appreciate the theoretical framework here and the new experimental data gathered. However, a few critical changes are required before publication.

1) The most critical change suggested by the reviewers is tempering and clarifying the relationship between experimental results and theory here. The theoretical work here is already sufficiently interesting, and the experimental results add to it, but the current statements imply a consistency between theory and experiments that is not fully supported.

The authors should temper their claims and clearly delineate the ways in which the experimental results agree and disagree with the theoretical model here. (`All models are wrong, but some are useful.') Such a clear statement would allow future work to build on the results here.

2) The authors should discuss other predictions of their model for future experiments (e.g., see reviewer #2's and #3's comments).

3) Simplify/reduce technical jargon – e.g., see reviewer #1's comments.

Reviewer #1 (Recommendations for the authors):

Appraisal of results supporting Figure #1 and theoretical background:

The theoretical background is sound and detailed. However, much of it rests on the choice of parameters as identified by the authors. It would be beneficial to see the dependence of the existence of saddle-node pitchfork bifurcations (and similarly the stable polarisation state, and homogeneous symmetric steady states) on key parameters (especially the total EGFR on the plasma membrane, and the total amount of the membrane-localized PTPRG and the ER-bound PTPN2).

Also, the authors claim that in their model, "The cell polarity is sustained even when the EGF signal is briefly disrupted, but also, the cell is able to rapidly reverse the direction of polarization when the signal direction is inverted". It would be helpful if the authors indicated (either by a numerical plot, or analytical calculations if possible) the time scales of such responses in their model. Chiefly, it would help to know how long the signal would need to be disrupted before the polarity is no longer sustained, how long the signal direction needs to be inverted before the cell reverses polarity, and how long it takes for the cell to either reverse or lose polarity.

Appraisal of results supporting conclusion #1:

The authors claim that MCF7 cells maintain a memory of the direction of previously encountered signals through prolonged EGFR phosphorylation polarisation and temporal evolution of the cell protrusion that emerges from a saddle-node pitchfork bifurcation which maintains the system away from the steady states. The key pieces of evidence present claim that after subjecting the cells to a stable gradient of EGF for one hour, the cells remain polarised for ~40 minutes. While the experimental evidence shows that the cells are partially polarised for ~40 minutes after gradient washout, the polarisation is less than it is in the presence of the signal and steadily approaches the unpolarised state (this is most clearly seen in the single-cell trajectories supplied in Supp. Figure 2-1F). This is unlike the distinct persistence of polarisation shown in the numerical results (Figure 1F and G). To me, the experimental data suggest that the timescale of relaxation from the polarized state to the unpolarized state due to slow dephosphorylation kinetics is ~40 minutes. The authors should describe what the characteristic chemical timescales of the system are, and if persistence of ~40 minutes is uncharacteristically large. How do the experimental results reflect the existence of a metastable "ghost" state as in saddle-node bifurcations, rather than two steady states that take ~40 minutes to move between them as in a subcritical pitchfork bifurcation without a saddle-node? The authors attempt to answer this by using Lapatinib later on (which should be brought into this section), but the difference due to Lapatinib isn't established as data is shown only for one cell, and further, the unintended effects of Lapatinib on the cells hasn't been shown (a possible way to demonstrate this would be a control experiment with cells treated with Lapatinib for a short amount of time and then exposed to a gradient of EGF to show that there is a memory in such a case).

Appraisal of results supporting conclusion #2:

The authors claim that the transient memory encoded in the EGFR phosphorylation polarisation and cell protrusion translates to memory in directional migration up a linear EGF gradient that persisted after the gradient washout. The authors demonstrate that there is a transient memory in directional migration after subjecting the cells to 5 hours of a linear gradient (however, they do not initially identify the period of this memory in Figure 3A and 3B. They later suggest that it is 50 minutes after gradient washout); however, it is difficult to see if the memory encoded in the EGFR phosphorylation polarisation and cell protrusion is accurately reflected in the memory in directional migration due to the vast differences in the timescales of the experimental evidence for the two (the former is in a constant exposure for 1 hour, while the latter is for a constantly decreasing exposure for 5 hours). Further, the evidence using Lapatinib is not conclusive, as even if we believe that there is a difference between the directionality with and without Lapatinib (Figure 3J), this does not establish the link between transient molecular states and transient directionality. The evidence in Figure 3G is inconclusive as it is not sufficiently different from the case without Lapatinib and Takens's delay-embedding trajectory is based on just one cell.

A key piece whose repercussions were not discussed much by the authors is that the gradient steepness was progressively decreased. The authors treated the entire period before the gradient washout as being the same, and I am still unclear as to what the effect of decreasing the gradient was.

Appraisal of results supporting conclusion #3:

The authors claim that MCF10A cells can sense dynamic gradients while maintaining robust directional migration even when the signal was disrupted. This is seen very clearly in silico in Figure 4B, but the experimental results in Figure 4C clearly differ in that the cells do not seem to have any significant memory following the gradient washout. Further, the ability of the cells to respond to changes in the environment should be contrasted to the case of slow dephosphorylation kinetics to highlight the role of the transient memory state. In general, it would be helpful to understand what the effect of a shorter time scale of the transition from a polarised to unpolarised state would be on the robustness of the system, and what the effect of a longer time scale would be on the adaptation of the system.

Note regarding the Discussion:

The authors do not distinguish their contributions from those made by other research groups that are addressing similar questions. For example, in line 310, the authors cited a set of papers as relevant experimentally but dismissed them as not providing a mechanism even though some of them did. For instance, Skoge et al. (PNAS, 2014), discussed a possible mechanism by which cells exhibit memory and maintain directed motion in the classical back-of-the-wave problem for Dictyostelium, and Prentice-Mott et al. did so for chemotactic neutrophil-like cells. Similarly, in lines 28-30 and lines 148-149 they cite a host of different papers that also model directional sensing in eukaryotic chemotaxis that they dismiss but should address. They state that these models cannot capture memory in polarization along with continuous adaptation or require fine-tuning. But it is not apparent that this is the case. The work by the authors is sufficiently different from many of these other papers and their contribution is significant that it merits dissemination (even if as a different framework for the same problem). But they should acknowledge clearly what their specific contribution is, or how it performs better than other models. If possible, they should cite experimental evidence that distinguishes their predictions from those of other models.

Detailed notes regarding figures and text:

The figures and the captions are very unclear, and the text can be bigger. The captions do not properly explain what is happening (please use complete sentences for captions). A few (but not exhaustive) notes:

Figure 1: Please label the color notation for the bars indicating the direction of the gradient (red, green, and orange) or indicate it somewhere in the caption. A plot similar to Figure 2B would be very helpful.

Figure 1A: This figure is highly confusing without the text. As such it served no purpose for me, and even after reading the text and understanding what it was trying to say, I didn’t understand the figure itself. The caption labels are inconsistent with the coloring in the plot, and there’s a notion of time that was not conveyed in the figure. It might serve the reader to split this into two figures and make the figure captions more explanatory.

1B: Not all the different elements are not explained in the caption (basically everything inside the cell is unlabelled).

1C: The circuit is too small to see.

Supp Figure 1-1: Please label the figures in the correct order (the panels labeled E should alphabetically appear after F and G). Also, it might be beneficial to keep the same scale for Ep (or switch to Ep/Et while identifying Et in the figure itself).

It was difficult to find the parameters used for the IHSS. I would recommend a table of parameters with columns for the different sets of simulations performed, with descriptions of each of the parameters/variables.

1-1B The current description gives the impression that Et is titratable (but as I understand it, the total protein concentration was conserved in the mass-action kinetics) especially as it is placed along with Ep which does change in the simulations. Please elaborate on the description in the figure caption and distinguish Ep and Et as a variable and parameter respectively.

Figure 2A: Please indicate a color bar for EGF just as was provided for EGFRp.

2C: During what time was the spatial projection generated?

D, E, H: Please spell out that n and N refer to the sample size and the number of replicates for the general readership.

2F: The figure is barely legible. Please magnify, and consider normalizing E-Ep.

2G: What are the units for the area? Is it a fraction of the total area?

2H: What does the horizontal line at the solidity of 0.65 represent? Also, how is the "end of memory" determined?

Figure 3A: Please identify the memory phase and label it in blue.

Figure 4D Please add a horizontal line for cos\theta = 0.

4E This figure was very difficult to read. Please consider breaking it into multiple 2D plots as one of the dimensions in the plot currently serves no purpose

Note for Section 1:

This section is fairly heavy on jargon and technical language. While I appreciate the transparency, it might benefit the general reader to reduce the number of terms used. For example, stable inhomogeneous state regime, stable polarization, and inhomogenous steady-state are used interchangeably in the text and the figures. Sticking to a small set of terminology would be beneficial to the reader.

Notes for Section 2:

I would define solidity early on (as of now, the first definition appears in line 557, but solidity appears as a quantity as early as Figure 2). In line 168, I would clarify if it is the polarisation that is shallow or the gradient of EGF. Also, a key piece of evidence that is currently buried is the timescale of gradient washout. It would serve the reader to highlight it in the text.

Please also clarify the relevance of the Takens delay-embedding trajectory. Currently, I found it misleading since it merely states that is a period of transition between two distinguishable states: the polarised state and the unpolarised state, which is already captured in the time plot of the EGFR dynamics.

Notes for Section 3:

It is currently stated that "directed migration persisted for a transient period of time after the gradient wash-out". This doesn't seem to be quantified, and what constitutes a "transient" period is unclear. Thus, the statement ,"The directionality estimated in the 9h time-frame after the gradient removal was greater than the one in continuous stimulus absence" seems to contradict the statement that "After the memory phase, the cells transitioned to a migration pattern equivalent to that in the absence of a stimulus" as the duration of the transient period of transition is unspecified. Also, this statement is not adequately reflected in the statistics of the directionality (Figure 3B).

Notes for Section 4:

This section was well-written, with the exception of describing what the KDE distributions reveal. This is a key piece of statistical evidence that must be more clearly shown and its relevance discussed.

Reviewer #2 (Recommendations for the authors):

– In general, the experiments conducted support the key features of the dynamical model. It would strengthen the authors' conclusions if the effects of perturbations (e.g., by Lapatinib) could be clarified in the main text within the context of the model. For instance, is only the memory lost, or is the cell's ability to polarize in the presence of a gradient also disrupted?

– It would be helpful if one did not have to refer to the caption to read several of the figure panels. (An example: the color-coding in 3A, D).

Reviewer #3 (Recommendations for the authors):

Timescale of memory

In my understanding, the lifetime of a 'ghost' state near a saddle-node bifurcation goes as 1/r^(1/2), where r is the distance to the critical point. This suggests that the timescale associated with the memory state is sensitive to how close the system is to the critical point, and, at least formally, diverges as the system actually approaches this point. This would seem to present a fine-tuning problem: not only does the system need to have parameters tuned to be near criticality, but in fact they have to be exactly the right distance away to achieve a physiologically reasonable intermediate memory timescale. It would be useful for the authors to discuss this: how is the memory timescale controlled? How sensitive is it to parameter changes? Is this somehow a feature, in that the cell can potentially physiologically change its memory timescale?

Suggestions regarding the theoretical exposition:

1. The authors should clarify the dimensional reduction that led to the bifurcation diagrams in Figure 1A and Figure 1—figure supplement 1B. Absent additional justification of this reduction, it would also be clearer to describe this as an approximate treatment of the system (whose behavior is born out by the reaction-diffusion simulations), rather than a proof (line 100).

2. It would be useful to readers to include a more concrete discussion of the EGFR model in the main text, including which features of it drive the behavior of the system and which biochemical parameters control location in the phase space. Additionally, how does the magnitude of the external gradient affect the cell? As a suggestion, the authors could consider moving Methods 5.15, Equation 17, to the main text, with a description of what the important variables are, and what features are important to the presence of the pitchfork bifurcation.

Suggestions regarding the presentation of the measurements:

1. Related to the public review, the authors could strengthen the paper by explicitly discussing the discrepancies between the measurements and the expectations from the model, and potential explanations. Does experimental noise interfere with the EGFR phosphorylation profile measurements? Do the authors believe that some of the cells are not at criticality?

2. The fact that the model depends on cells being biochemically poised near a critical point suggests a variety of stronger experimental tests of the framework. As one example, the authors' analysis suggests that overexpressing EGFR should push the cells away from the critical point, into the inhomogeneous steady-state regime, where they would break symmetry according to the first-encountered gradient and no longer be capable of adapting to a new gradient. This would be quite surprising, as it would correspond to breaking a sensing system by increasing the number of sensors. While performing this experiment is likely beyond the scope of this paper, the authors could strengthen the presentation by discussing this and/or other more counterintuitive predictions of their model in light of existing empirical data and/or future experiments.

3. Figure 2—figure supplement 1C: given that the direction of the polarization relative to the gradient is important, it would be interesting to see all the polarization profiles (and the variability from cell to cell with respect to the direction relative to the gradient).

https://doi.org/10.7554/eLife.76825.sa1

Author response

Essential revisions:

The reviewers broadly appreciate the theoretical framework here and the new experimental data gathered. However, a few critical changes are required before publication.

1) The most critical change suggested by the reviewers is tempering and clarifying the relationship between experimental results and theory here. The theoretical work here is already sufficiently interesting, and the experimental results add to it, but the current statements imply a consistency between theory and experiments that is not fully supported.

The authors should temper their claims and clearly delineate the ways in which the experimental results agree and disagree with the theoretical model here. (`All models are wrong, but some are useful.') Such a clear statement would allow future work to build on the results here.

2) The authors should discuss other predictions of their model for future experiments (e.g., see reviewer #2's and #3's comments).

3) Simplify/reduce technical jargon – e.g., see reviewer #1's comments.

We thank the referees and the editor for their insightful and detailed suggestions that we directly followed in order to generate the amended version of the manuscript. We now particularly discuss the similarities and possible differences between the theory and the experimental observations, more strongly highlight several proposed predictions, but also provide additional predictions of our model with respect to cellular sensing and responsiveness to simultaneous stimuli. We also aimed to reduce the technical jargon and in general, restructure the results as suggested by the referees, as well as provide additional conceptual schemes in order to present the proposed dynamical mechanism and the results in a clear manner to improve the readability of our manuscript. We next provide detailed response to the comments and concerns raised by the referees.

Reviewer #1 (Recommendations for the authors):

Appraisal of results supporting Figure #1 and theoretical background:

The theoretical background is sound and detailed. However, much of it rests on the choice of parameters as identified by the authors. It would be beneficial to see the dependence of the existence of saddle-node pitchfork bifurcations (and similarly the stable polarisation state, and homogeneous symmetric steady states) on key parameters (especially the total EGFR on the plasma membrane, and the total amount of the membrane-localized PTPRG and the ER-bound PTPN2).

Indeed, as the referee suggested, the two key parameters that determine the dynamical behavior of the system are the total EGFR and PTPRG concentrations on the plasma membrane of the cell, and we would like to include in this response also the response regarding the Figure 1—figure supplement 1C (bifurcation analysis of the EGFR network) that is raised later by the referee.

The dynamical behavior of the system is dependent on the ratio between EGFR and PTPRG concentrations at the plasma membrane, as PTPRG exhibits the main dephosphorylating activity of the receptor at the plasma membrane. In Figure 1—figure supplement 1B (Figure 1—figure supplement 1C in the amended version of the manuscript), we have chosen to show a bifurcation diagram with respect to the total EGFR concentration at the plasma membrane (EGFRt), since the structure of the bifurcation diagram remains the same whether one or the other parameter is varied. In order to substantiate this fact, we provide here a two-parameter bifurcation diagram of EGFRt vs. PTPRGt (Author response image 1A) where we demonstrate that there is a large parameter range where pitchfork bifurcation (PB) can be observed in the system. Additionally, we also provide a two-parameter bifurcation diagram of EGFRt vs. PTPN2t (Author response image 1B). This demonstrates that PB always exist above a threshold PTPN2t concentration.

Author response image 1
In-depth bifurcation analysis of the EGFR network.

(A) Two-parameter bifurcation analysis (EGFRt vs. PTPRGt) characterizing the range where pitchfork bifurcations (black line) emerge. (B) Same as in a, only when PTPN2t instead of PTPRGt is considered.

As the referee notes, the total EGFR concentration is conserved in the system, and our hypothesis suggests that the EGFRt concentrations in cells correspond to organization at criticality (before the SNPB, gray shaded region in Figure 1—figure supplement 1C in the amended version). In order to understand how the signal induced transition from stable unpolarized to stable polarized state emerges, we have used EGFRt as a bifurcation parameter to reveal all possible dynamical regimes of the system.

For an EGFRt concentration that correspond to the organization at criticality, the system displays a basal, unpolarized state, as it is organized in the homogenous steady state (HSS) regime (Figure1—figure supplement 1C). Upon presentation of a localized signal however, a characteristic unfolding of the Pitchfork bifurcation [1] occurs (Author response image 2). The HSS are pushed to the sides and the IHSS regime dominates (a general characteristic of IHSS solutions,[2]). Thus, for the same EGFRt concentrations, the basal unpolarized state is now unstable and the system populates the IHSS branch – polarizes. Upon the signal removal, the system will return to the initial HSS regime, however in this case, transiting via the metastable “ghost” state, as schematically shown in Figure 1—figure supplement 1A in the amended version of the manuscript.

Author response image 2
Generic principles of the symmetry breaking mechanism.

Unfolding of the pitchfork bifurcation – the bifurcation diagram of the full system as in a., only when a local EGF gradient signal is present.

Due to the extent of this theoretical analysis and additional concepts necessary to explain the dynamical features of the system, we have chosen only to briefly discuss this in the methods (Section 5.16) and provide a schematic representation of the unfolding in the current version of Figure 1—figure supplement 1A, as well as a simplified version of the transitions in the current version of Figure1A. We strongly believe that including the discussion on these transitions in greater details will distract the reader from the main message of the paper. However, we are currently preparing an additional theoretical manuscript where we explain these details, as well as how it differs from the current models in the field (i.e. Turing, LEGI and Wave-pining) that will be available on bioRxiv in the upcoming month.

The model parameters used result from a detailed experimental characterization of the EGFR network from our previous, as well as experimental work of other groups [3,4,5], including the ligand binding/unbinding rates, recycling time of unphosphorylated receptors, the rates of autonomous/autocatalytic activation of the receptors, the specific phosphatases’ reactivity towards the receptor etc.

Also, the authors claim that in their model, "The cell polarity is sustained even when the EGF signal is briefly disrupted, but also, the cell is able to rapidly reverse the direction of polarization when the signal direction is inverted". It would be helpful if the authors indicated (either by a numerical plot, or analytical calculations if possible) the time scales of such responses in their model. Chiefly, it would help to know how long the signal would need to be disrupted before the polarity is no longer sustained, how long the signal direction needs to be inverted before the cell reverses polarity, and how long it takes for the cell to either reverse or lose polarity.

a) When the disruption of the signal is longer than the memory duration in single cells, the cells will not be able to integrate the information and establish long-term migration. In this case, a so called “stop-andgo” migration pattern will be established. In the amended version of the manuscript, we provide numerical simulations (Figure 4—figure supplement 1B and Figure 4 – video 2), as well as include corresponding discussion in the text.

b) In the case of reverse polarization, cell lose the polarity after the memory in polarization is lost (~20min on average after the removal of the third signal). The duration of the memory in the reversed polarization is quantified in Figure 4E.

c) We thank the referee for this suggestion. We have quantified the time necessary for cells to reverse polarity after presentation of the reversed signal. As shown in Author response image 3, from the KDE of the projection of the cell’s relative displacement angles during the presentation of the third, reversed signal (time point 365min), it can be seen that on average, majority of the cells reverse their direction in migration within 10 min after the reversal signal is presented (yellow, cosθ takes negative values). This information is now included in the manuscript.

Author response image 3
Quantification of the time intevbal in which cells revert their direction of migration, as extracted from KDE distributions of reversal of gradient migration data in Figure 4 —figure supplement 1F.

Appraisal of results supporting conclusion #1:

The authors claim that MCF7 cells maintain a memory of the direction of previously encountered signals through prolonged EGFR phosphorylation polarisation and temporal evolution of the cell protrusion that emerges from a saddle-node pitchfork bifurcation which maintains the system away from the steady states. The key pieces of evidence present claim that after subjecting the cells to a stable gradient of EGF for one hour, the cells remain polarised for ~40 minutes. While the experimental evidence shows that the cells are partially polarised for ~40 minutes after gradient washout, the polarisation is less than it is in the presence of the signal and steadily approaches the unpolarised state (this is most clearly seen in the single-cell trajectories supplied in Supp. Figure 2-1F). This is unlike the distinct persistence of polarisation shown in the numerical results (Figure 1F and G).

In order to address the referee’s comments on the evidence of existence of a SNPB “ghost”, we combine here the answer to this, with the answer to the next comment. The “ghost” state is manifested as transient trapping of the trajectory of the system (depicted through the EGFRp measurements) in the vicinity of the previously stable polarized steady-state (Figure 1F for numerical, and 2F for experimental results), before rapidly resetting to the basal state. In the amended version of the manuscript we provide several additional EGFRp time series profiles from single cell measurements (Figure 2—figure supplement 2 E, F). Note that the shape of the EGFRp profiles w/o Lapatinib treatment is similar to the one shown from the numerical simulation. The absolute levels of EGFRp in the presence of the signal (IHSS state) and after signal removal (memory state, SNPB “ghost”) slightly differ, as also reflected in change of the polarized area after signal removal (Figure 2 —figure supplement 1G).

However, the EGFRp levels at the “ghost” state are transiently steadily maintained, before transiting to the basal unpolarized state. This is significantly different in the case when the kinase activity of the receptor is inhibited using Lapatinib and thereby the “ghost” state is lost (Figure 2G, average from single cells in Figure 2 —figure supplement 2B).

To provide a systematic analysis of the EGFR phosphorylation profiles and their features with and without Lapatinib washout (Figures 2 F, G), we fitted first single-cell EGFRp measurements in both cases with an inverse sigmoidal function (x)=aoan+xn, and quantified the Hill coefficient n, as well as the time frame in which 50% of EGFR phosphorylation is lost (half-life) after EGF signal removal (Figure 2 —figure supplement 2E-G). Under normal conditions, ~30min were necessary to decrease EGFRp values by 50%, corroborating the presence of memory in EGFRp after signal removal, whereas the Hill coefficient was ~2.88, demonstrating that after the memory is lost, there is a rapid transition to the basal state (median R2 of the fits ~0.79). The same function (with a, ao, n chosen from the fit corresponding to the previous case) however could not describe the EGFRp profiles in the Lapatinib treatment experiment, as manifested by the poor goodness-of-fit values (median R2 ~ 0.33). This indicates that the shape of the EGFRp profiles in both cases is different. We therefore fitted the data in the Lapatinib case fixing only a to the values obtained by the control case, and leaving ao and n as free parameters, as they determine the upper plateau and the steepness of the drop to the basal level. In this case, a Hill coefficient of n ~ 1.28 was identified from the fitting (median R2 ~ 0.84), whereas EGFRp levels dropped by 50% within 10min after EGF removal. These results clearly demonstrate that EGFRp is transiently maintained at steady levels after gradient removal, before rapidly transiting to the basal state (property of a “ghost” state), whereas in the Lapatinib case, the system gradually transits to the basal state, reflecting a continuous dephosphorylation process.

Lapatinib is an ATP-analog EGFR inhibitor that locks the tyrosine kinase domain into an inactive conformation, however it does not prevent EGF-induced dimerization of the receptor, as it has been demonstrated using a conformational EGFR sensor [6]. Thus, Lapatinib only affects the kinase activity and therefore, the EGFR phosphorylation profile after Lapatinib stimulation will be determined mainly by the catalytic activity of the phosphatases.

The catalytic activity of fully active PTPs is however two to three orders of magnitude higher than that of tyrosine kinases [18]. Thus, when the system is governed mainly by the phosphatase activity under Lapatinib treatment, EGFR phosphorylation will be reduced to 50% from the maximum level in ~10min. In contrast, the memory in EGFR phosphorylation polarization after removal of the EGF signal is maintained for ~40min on average (up to ~150min in single cells Figure 2—figure supplement 1G), indicating that the transient maintenance of EGFR phosphorylation cannot be explained with slow dephosphorylation.

To me, the experimental data suggest that the timescale of relaxation from the polarized state to the unpolarized state due to slow dephosphorylation kinetics is ~40 minutes. The authors should describe what the characteristic chemical timescales of the system are, and if persistence of ~40 minutes is uncharacteristically large. How do the experimental results reflect the existence of a metastable "ghost" state as in saddle-node bifurcations, rather than two steady states that take ~40 minutes to move between them as in a subcritical pitchfork bifurcation without a saddle-node? The authors attempt to answer this by using Lapatinib later on (which should be brought into this section), but the difference due to Lapatinib isn't established as data is shown only for one cell, and further, the unintended effects of Lapatinib on the cells hasn't been shown (a possible way to demonstrate this would be a control experiment with cells treated with Lapatinib for a short amount of time and then exposed to a gradient of EGF to show that there is a memory in such a case).

A) With respect to the first part of the comment, we would like to refer back to the answer on the previous comment. Moreover, the presence of the “ghost” state in the control and the respective absence in the Lapatinib treated cells is also reflected through the exemplary state-space trajectory reconstructions, where a three- vs. two-state transitions respectively are identified. We also thank the referee for the suggestion to include the Lapatinib characterization in Figure 2, which we followed in the amended version of the manuscript.

B) We thank the referee for the suggestion regarding the possible control experiment. However, as Lapatinib completely blocks the kinase activity of the receptor, pre-stimulation with Lapatinib renders the cells unresponsive to growth factor stimulation for ~96h [7] and therefore this experiment would not yield a possible control in our case. To address this point raised by the referee, we therefore resorted to quantification of the EGFRp temporal profile characteristics in the control and Lapatinib cases, as well as the time interval in which EGFRp is reduced by 50% in single cells after EGF removal, as discussed above.

Appraisal of results supporting conclusion #2:

The authors claim that the transient memory encoded in the EGFR phosphorylation polarisation and cell protrusion translates to memory in directional migration up a linear EGF gradient that persisted after the gradient washout. The authors demonstrate that there is a transient memory in directional migration after subjecting the cells to 5 hours of a linear gradient (however, they do not initially identify the period of this memory in Figure 3A and 3B. They later suggest that it is 50 minutes after gradient washout); however, it is difficult to see if the memory encoded in the EGFR phosphorylation polarisation and cell protrusion is accurately reflected in the memory in directional migration due to the vast differences in the timescales of the experimental evidence for the two (the former is in a constant exposure for 1 hour, while the latter is for a constantly decreasing exposure for 5 hours).

We would like to note here that the memory in directional migration has been estimated from the cell migration traces in Figure 3A (Figure 3 —figure supplement 1D), by quantifying the projection of the cell’s relative displacement angle along the gradient direction (quantification in Figure 3C). Although it has not been technically possible to quantify EGFR phosphorylation polarization and cell migration simultaneously due to the necessity to use different magnification during imaging in both cases, we have addressed the issue whether the memory in EGFR phosphorylation polarization results in memory in directional migration by performing migration experiments under Lapatinib addition, as Lapatinib completely hinders EGFR kinase activity (Figure 3G-I). In addition to the provided directionality quantifications (Figure 3H), in the amended version of the manuscript we provided as well quantification of the projection of the relative displacement angles in the Lapatinib condition (Figure 3I, Figure 3 —figure supplement 2H), as well as additional quantification of memory in directional migration from single-cell trajectories in both cases (Figure 3 —figure supplement 3).

Moreover, we would like to note that regardless of the difference in the time for gradient stimulation, once the system reaches stable polarized state, the transitions that follow after the gradient removal to the basal state through the metastable state remain equivalent, as these are generic properties of the system.

Further, the evidence using Lapatinib is not conclusive, as even if we believe that there is a difference between the directionality with and without Lapatinib (Figure 3J), this does not establish the link between transient molecular states and transient directionality.

Lapatinib blocks the EGFR kinase activity without any effect on its structural properties [6], thus the difference in memory in directional migration after gradient wash-out in the case with and without Lapatinib establishes the link between memory in polarized EGFR signaling and memory in directional migration.

The evidence in Figure 3G is inconclusive as it is not sufficiently different from the case without Lapatinib and Takens's delay-embedding trajectory is based on just one cell.

We note that the averaged single-cell EGFRp profiles upon Lapatinib treatment are shown in the current Figure 2—figure supplement 2B. Together with the additional quantifications on the shape of EGFRp profile with and without Lapatinib, the results show that in the Lapatinib case, EGFR phosphorylation decreases after EGF removal in an exponential fashion. This corresponds to direct transition between the polarized and the basal state without intermediate phase-space trapping, as observed in the exemplary reconstructed state-space trajectory (now in Figure 2G).

In order to further verify that such transient trapping in phase space does not occur in the case when EGFR phosphorylation dynamics is guided mainly by the dephosphorylation of phosphates, such as under the Lapatinib treatment, we have performed additionally numerical simulations where we reduced the EGFR kinase activity through its autocatalysis rate after gradient removal, to mimic the effect of Lapatinib (α2 = 0.25). The results included in Figure 2 —figure supplement 2C, D as well as Figure 2 – video 3, show that equivalently as in the experiments, EGFRp gradually transits from the polarized to the basal state, without intermediate state-space trapping.

A key piece whose repercussions were not discussed much by the authors is that the gradient steepness was progressively decreased. The authors treated the entire period before the gradient washout as being the same, and I am still unclear as to what the effect of decreasing the gradient was.

One of the basic features for efficient polarization mechanism is the capability of sensing steep and shallow gradients with an offset. In order to demonstrate that the proposed mechanism fulfills these conditions, both in the numerical simulations (Figure 1D) as well as in the migration experiments (Figure 3,4), we established a dynamic gradient whose steepness changes over time. The results show that the SNPB mechanism enables sensing gradients that dynamically change in time. This is a more complex form of a spatial signal as compared to the static gradient, for which cells also establish directional polarization and migration. We have used a static gradient as a second gradient in the experiments (and simulations) shown in Figure 4.

Appraisal of results supporting conclusion #3:

The authors claim that MCF10A cells can sense dynamic gradients while maintaining robust directional migration even when the signal was disrupted. This is seen very clearly in silico in Figure 4B, but the experimental results in Figure 4C clearly differ in that the cells do not seem to have any significant memory following the gradient washout.

Although we fully agree with the referee that there is discrepancy between the simulations and experiments, we would like to underline here that the discrepancy refers to the migration distance between the signals, and not a lack of memory in the experimental case. In order to verify this, we have calculated the displacement of the cells in the intervals when the signal is disrupted from the experimental tracks in Figure 4 —figure supplement 1F and this ranges on average between 9-30μm, in some cells even 45μm. That cells migrate and maintain the directionality in these intervals can be also seen from the respective KDE distributions (Figure 4E). We have also included several magnified single cell traces in the amended version of Figure 4 —figure supplement 1G to support these statements.

Further, the ability of the cells to respond to changes in the environment should be contrasted to the case of slow dephosphorylation kinetics to highlight the role of the transient memory state. In general, it would be helpful to understand what the effect of a shorter time scale of the transition from a polarised to unpolarised state would be on the robustness of the system, and what the effect of a longer time scale would be on the adaptation of the system.

As already described, due to the inhibition effect of Lapatinib on the kinase activity of the receptor, the case when cells were subjected to Lapatinib during gradient wash-out describes the behavior of the system when guided mainly by dephosphorylation (absence of a “ghost” state). That in this case, the temporal EGFRp profile has a distinct shape, as well as there is no intermediate state-space trapping was verified both experimentally, as well as with numerical simulations (new Figure 2 —figure supplement 2C, D as well as Figure 2 – video 3). However, in order to test what will be the effect on cell migration if hypothetically dephosphorylation would be slower, we consider here additionally an exponentially decaying model where the dynamics of receptor phosphorylation Rp is given by Equation 1:

Rpt=l(x,t)α(RtRp)βRp+DRpRpx2(1)

I(x,t) is the spatial-temporal signal input, α is an activation term, β – the dephosphorylation constant, DRp – the diffusion constant, and Rt – the total receptor amount. The value of β was chosen to capture a slow dephosphorylation profile after the signal removal (Author response image 4A, left), in contrast to transition through a metastable state in the EGFR-PTP model (Figure 4A in the response, right). In order to demonstrate the difference this brings to spatial-temporal signal processing, two distinct sequences of gradient fields were used (Author response image 4B), where the first two gradients (green) the same timing and duration, whereas the timing of the third gradient (red and blue) differs.

Author response image 4
Optimal response to changing environments is rendered by metastable dynamics.

(A) Response to single gradient stimulation using an exponential decay model for receptor activity dynamics (Equation 1) with slow dephosphorylation kinetics (left); and using the EGFR-PTP model with metastable “ghost” state. (B) Two distinct sequences of spatial-temporal gradient fields. (C) Directed protrusive area estimated from in-silico cell shape changes during stimulus sequence 1 (red) and sequence 2 (blue) for the exponential decay model (left) and EGFR-PTP model (right). Parameters for the exponential decay model: α = 1.5, Rt = 1, β = 0.5 × 10−2.

We performed viscoelastic simulations as described in Figure 1C in the manuscript, using the EGFR-PTP as well as the exponential decay model, to model the receptor activity dynamics, and quantified the morphological changes of the cell in both cases for both signals, using the directed protrusive area as in Figure 1G in the manuscript. When the receptor dynamics is guided by a slow dephosphorylation as in the exponential decay model, the cell morphology remains almost identical for the two distinct spatial-temporal signals (red and blue lines in Author response image 4C left correspond to morphological changes to signals 1 and 2 respectively), suggesting that the cell cannot distinguish between the two distinct signals. In contrast, when the receptor dynamics is guided by a “ghost” memory state in the EGFR-PTP model, the cell morphology directly follows the dynamics of the input signals (red and blue lines in Author response image 4C right correspond to morphological changes to signals 1 and 2 respectively), suggesting that this is a unique mechanism for cellular adaptation to varying spatial-temporal signals.

Note regarding the Discussion:

The authors do not distinguish their contributions from those made by other research groups that are addressing similar questions. For example, in line 310, the authors cited a set of papers as relevant experimentally but dismissed them as not providing a mechanism even though some of them did. For instance, Skoge et al. (PNAS, 2014), discussed a possible mechanism by which cells exhibit memory and maintain directed motion in the classical back-of-the-wave problem for Dictyostelium, and Prentice-Mott et al. did so for chemotactic neutrophil-like cells.

We thank the referee for this suggestion. In the amended version of the manuscript, we have included a more detailed discussion how the mechanism of working memory proposed here differs with the work on the “back-of the-wave” problem in Dictyostelium and chemotactic behavior of neutrophil-like cells. In brief, in Skoge et al., 2014 [8], the authors propose mechanism that relies on LEGI model coupled to a memory-module (positive feedback leading to bistability). The duration of the memory is controlled (memory is turned on and off) through two independent thresholds. The authors however have not proposed a biochemical mechanism how the thresholds or their regulation can be realized, as they are crucial to enable a temporal memory in the polarized signaling activity. We have previously demonstrated that a bistability mechanism, when irreversible, will result in a permanent temporal memory. If the system displays reversible bistability, only memory of the input signal concentration can be exhibited, but not a temporal memory in the activity [9]. Moreover, in Skoge et al., 2014, potential molecular elements that store the information have also not been identified.

In the work of Prentice-Mott et al. [10], a novel mechanism of directional memory in migrating cells resulting from distinct molecular time scales has been proposed. The phenomenological model of transient memory in the cytoskeleton has however not been directly linked to a detailed biochemical mechanism. In the manuscript we present, we provide direct measurements of memory in receptor phosphorylation polarization which directs the memory in directional migration. Moreover, in the amended version of the manuscript, we included single-cell quantification of the duration of memory in directional migration in the case where the kinase activity has been blocked using Lapatinib, which effectively diminishes the memory in receptor phosphorylation. The obtained results are now discussed in relation to the results described in Prentice-Mott et al.

Similarly, in lines 28-30 and lines 148-149 they cite a host of different papers that also model directional sensing in eukaryotic chemotaxis that they dismiss but should address. They state that these models cannot capture memory in polarization along with continuous adaptation or require fine-tuning. But it is not apparent that this is the case. The work by the authors is sufficiently different from many of these other papers and their contribution is significant that it merits dissemination (even if as a different framework for the same problem). But they should acknowledge clearly what their specific contribution is, or how it performs better than other models. If possible, they should cite experimental evidence that distinguishes their predictions from those of other models.

We have decided to omit a detailed discussion on the differences between the SNPB mechanism and previously proposed polarization mechanisms, as in our opinion, this requires an extensive analysis and a detailed discussion which would have significantly extended the scope of the manuscript. We have however performed such analysis by comparing several polarization features (i.e. capability for sensing steep and shallow gradients, threshold activation / spurious polarization to noise, signal amplification, time of polarization etc.) as well as information processing features (sensing subsequent signals, rapid repolarization, signal integration etc.) between the proposed SNPB mechanism and several of the most common polarization models. i.e. Turing model [11, 12], LEGI [13] and Wave-pinning model [14]. The features used for the comparison were drawn from commonly observed features of migrating cells, i.e. leukocytes, neutrophils, Dictiostelium etc. as summarized in an extensive review by the Keshet lab [19]. The main difference between previously published, and the mechanism we propose here, is the organization of the system in parameter space: whereas all previous models are organized in stable attractor regimes, the SNPB mechanism relies on organization at criticality and thereby utilizes metastable states for sensing and adaptation to changing signals. We present an example of these features in this response for the referee’s perusal, however, due to the extent of the analysis, we have chosen to summarize this in a related theoretical manuscript which will be posted on bioRxiv in the upcoming months.

As example we choose to present here the difference between the proposed SNPB mechanism, LEGI, Turing and Wave-pinning models for adaptation to changing signals. First, we quantified the efficiency of the 4 models to respond to signal reversal (Author response image 5A). Estimation of polarization amplification, which is the ratio of response amplitude at the leading edge of the cell before and after the reversal of gradient (RrevRfront) and the time required for reversing the polarization for each of the models demonstrated that rapid reversal is a feature that is unique for the SNPB mechanism (Author response image 5B). Neither the Wave-pinning model nor the Turing model were able to adapt the direction of polarization to the reversed gradient localization in physiologically relevant time, due to the permanent memory resulting from stable attractor organization. We also tested capability of the different models to integrate information from time-varying signals. Due to the same reason, Wave-pinning and Turing models show hindered responsiveness (response is equivalent whether one or two signals are presented, Author response image 5C, left). The LEGI-model, due to the lack of memory, fails to integrate the information and responds to each input individually, and only the SNPB model, due to the presence of metastable dynamics, shows capability to generate a response by integrating both signals (Author response image 5C, right).

Author response image 5
Comparison of information processing features between existing polarisation models (LEGI, Turing, Wavepinning) and the proposed SNPB mechanism.

(A) Schematic of gradient reversal. (B) Quantification of the time for polarization reversal and corresponding polarization amplification (RrevRfront) at stimulus ratio, (SrevSfront)=2. (C) Time series of signaling response at the leading edge (Rfront1) of the cell upon two consecutive gradients from the same direction for Wave-pinning and Turing models (left), and LEGI and SNPB model (right).

Detailed notes regarding figures and text:

The figures and the captions are very unclear, and the text can be bigger. The captions do not properly explain what is happening (please use complete sentences for captions).

We have now revised the figure legends, including all of the comments raised by the referee in order to improve the readability.

A few (but not exhaustive) notes:

Figure 1: Please label the color notation for the bars indicating the direction of the gradient (red, green, and orange) or indicate it somewhere in the caption. A plot similar to Figure 2B would be very helpful.

We have included schematic representation for the gradient direction in the relevant figures.

Figure 1A: This figure is highly confusing without the text. As such it served no purpose for me, and even after reading the text and understanding what it was trying to say, I didn’t understand the figure itself. The caption labels are inconsistent with the coloring in the plot, and there’s a notion of time that was not conveyed in the figure. It might serve the reader to split this into two figures and make the figure captions more explanatory.

Following the referee’s suggestion, we have provided a more schematic representation of the figure in the amended version of the manuscript in order to explain the basis of the dynamical mechanism discussed, included a more detailed schematic representation in the current version of (Figure 1—figure supplement 1A), and improved the explanations in the figure captions.

1B: Not all the different elements are not explained in the caption (basically everything inside the cell is unlabelled).

1C: The circuit is too small to see.

We have improved the size, added an additional legend and described in more details in the figure legend.

Supp Figure 1-1: Please label the figures in the correct order (the panels labeled E should alphabetically appear after F and G). Also, it might be beneficial to keep the same scale for Ep (or switch to Ep/Et while identifying Et in the figure itself).

It was difficult to find the parameters used for the IHSS. I would recommend a table of parameters with columns for the different sets of simulations performed, with descriptions of each of the parameters/variables.

We have revised the labeling in the figure and also added a table of the parameters/ variables for the signaling, as well as the physical model of the cell.

1-1B The current description gives the impression that Et is titratable (but as I understand it, the total protein concentration was conserved in the mass-action kinetics) especially as it is placed along with Ep which does change in the simulations. Please elaborate on the description in the figure caption and distinguish Ep and Et as a variable and parameter respectively.

The answer to this question we have addressed in the answer to comment 1 by the referee and provided additional bifurcation analysis (Figure 1 in this response) where we clarify the raised issues.

Figure 2A: Please indicate a color bar for EGF just as was provided for EGFRp.

Amended.

2C: During what time was the spatial projection generated?

The spatial projection is calculated in the gradient duration period (5min – 65 min), as a temporal average per spatial bin. The basis of this calculation is also given in the methods section.

D, E, H: Please spell out that n and N refer to the sample size and the number of replicates for the general readership.

This is amended.

2F: The figure is barely legible. Please magnify, and consider normalizing E-Ep.

We have now enlarged the figure, however we option to keep the E-Ep profile normalized to the maximal fraction of ligand-bound receptors that can be achieved for the given EGF concentration – the position of the cell within the gradient, as estimated from the experimental data. As shortly also noted in the manuscript, the ligand bound fraction exponentially decays after EGF removal [3] in contrast to the EGFRp profile, demonstrating that the memory in EGFR phosphorylation is not due to residual EGF-bound receptor.

2G: What are the units for the area? Is it a fraction of the total area?

The directed cell protrusion area was estimated by comparing single cell masks at two consecutive time points. Protrusions were considered if the area change was greater than 10 pixels or 1.2µm2 per time point. The directed cell protrusion area was then obtained using Aprot,frontAfrontAprot,backAback and thereby it is a fraction. The extended information is provided in the Methods section.

2H: What does the horizontal line at the solidity of 0.65 represent? Also, how is the "end of memory" determined?

The horizontal line is just a visual guide of the value of the solidity of cells before establishing the gradient. The memory duration in cell morphology was calculated from the single-cell solidity profiles, and corresponds to the time-point at which the solidity is below mean-s.d. estimated during gradient presence, and the relevant description is in Methods, section 5.11.

Figure 3A: Please identify the memory phase and label it in blue.

In the previous version of the manuscript, the duration of memory in directional migration was calculated from the KDE distribution of cosθ, that reflects the population distribution. In order to calculate the duration of memory in directed migration for single-cell tracks and color the tracks accordingly, in the amended version of the manuscript we have first smoothened the trajectories using the Kalman filter and identified the memory duration from the respective cosθ distribution for single cells (Methods, Section 5.12). The quantifications are given in Figure 3 —figure supplement 3, along with few color-coded trajectories. We have however chosen to keep the population averages in the main figure as before.

Figure 4D Please add a horizontal line for cos\theta = 0.

Amended.

4E This figure was very difficult to read. Please consider breaking it into multiple 2D plots as one of the dimensions in the plot currently serves no purpose

The unlabeled dimension in the figure denotes time and we have used the successive positioning of KDE distributions derived from the gradient migration experiment in different time intervals to show how they shift relative to the distribution from normal migration experiment which represent time = 0 min or pre-stimulation scenario. We apologize for failing to label this axis in the previous version of the manuscript, this is now amended in all KDE distribution plots, and the direction of time remains as previously, denoted above the distributions with an arrow.

Note for Section 1:

This section is fairly heavy on jargon and technical language. While I appreciate the transparency, it might benefit the general reader to reduce the number of terms used. For example, stable inhomogeneous state regime, stable polarization, and inhomogenous steady-state are used interchangeably in the text and the figures. Sticking to a small set of terminology would be beneficial to the reader.

We thank the referee for his suggestions. In the revised version, we have tried to minimize the number of jargon specific words and also describe the results in a clearer manner to improve the readability.

Notes for Section 2:

I would define solidity early on (as of now, the first definition appears in line 557, but solidity appears as a quantity as early as Figure 2). In line 168, I would clarify if it is the polarisation that is shallow or the gradient of EGF. Also, a key piece of evidence that is currently buried is the timescale of gradient washout. It would serve the reader to highlight it in the text.

We thank the referee for this suggestion. Solidity is now defined in the paper in relation to Figure 2, we have also specified in the text that the gradient of EGF is shallow (line 178), and added the information that the gradient is washed-out in 4-5min (line 169).

Please also clarify the relevance of the Takens delay-embedding trajectory. Currently, I found it misleading since it merely states that is a period of transition between two distinguishable states: the polarised state and the unpolarised state, which is already captured in the time plot of the EGFR dynamics.

A unique characteristic of the presence of a “ghost” state is a transient trapping of a phase space trajectory near the polarized steady state (Figure 1F, bottom). Reconstructing the phase-space trajectory from EGFRp traces in the control case indeed demonstrates a presence of a “ghost” state in phase space (Figure 2F, bottom), which is lost upon treatment with Lapatinib (Figure 2G, bottom). In order to further verify that such transient trapping in phase space does not occur in the case of a slow dephosphorylation, such as under Lapatinib treatment, as noted to one of the comments above, we also performed additionally numerical simulations where we reduced the EGFR kinase activity through its autocatalysis rate after gradient removal to mimic the effect of Lapatinib (α2 = 0.25). The results in Figure 2 —figure supplement 2C, D show that similarly to the Lapatinib experiments, the system directly transits from the polarized to the basal state without intermittent phase-space trapping.

Notes for Section 3:

It is currently stated that "directed migration persisted for a transient period of time after the gradient wash-out". This doesn't seem to be quantified, and what constitutes a "transient" period is unclear.

In the amended version of the manuscript we tried to improve the readability of the text in general, and highlighted the fact that the transient persistence of the directional migration (memory) after gradient wash-out is quantified in Figure 3C and is ~50min. We have additionally quantified the duration of memory in directional migration from single cell cos θ profiles (Figure 3 —figure supplement 3).

Thus, the statement ,"The directionality estimated in the 9h time-frame after the gradient removal was greater than the one in continuous stimulus absence" seems to contradict the statement that "After the memory phase, the cells transitioned to a migration pattern equivalent to that in the absence of a stimulus" as the duration of the transient period of transition is unspecified. Also, this statement is not adequately reflected in the statistics of the directionality (Figure 3B).

We apologize if the wording has caused confusion, which we try to amend in the revised version. To clarify, that directionality calculated after gradient removal is significantly different (p<0.001) from the one calculated under continuous EGF absence suggests that cell after gradient removal do not immediately revert to the random walk migration as characteristic in the continuous gradient absence, but transiently persist with directional migration (memory), leading to the increased directionality values. In Figure 3C we have quantified that the memory in directional migration after gradient wash out is ~50min. Additionally, from the KDE estimations in Figure 3 —figure supplement 2C, it can also be seen that the KDE estimated after gradient removal approaches the one estimated in continuous absence only after the ~50min of memory phase (compare red and grey distributions, and blue and green).

Notes for Section 4:

This section was well-written, with the exception of describing what the KDE distributions reveal. This is a key piece of statistical evidence that must be more clearly shown and its relevance discussed.

We have now included a more detailed description of Figure 4E in this section.

Reviewer #2 (Recommendations for the authors):

– In general, the experiments conducted support the key features of the dynamical model. It would strengthen the authors' conclusions if the effects of perturbations (e.g., by Lapatinib) could be clarified in the main text within the context of the model. For instance, is only the memory lost, or is the cell's ability to polarize in the presence of a gradient also disrupted?

Corresponding to the referee’s suggestion, we have re-structured and better clarified, especially the results regarding the Lapatinib experiments. From the solidity quantification under Lapatinib condition (now in Figure 2I) and the included exemplary cell masks at distinct time-points, it is evident that the polarized cell shape is quickly lost after EGF removal and Lapatinib administration, in contrast to the normal case (Figure 2H). This suggests that loss of memory in EGFR phosphorylation polarization is reflected in loss of memory in polarized cell shape, and subsequently loss of memory in directional cell migration (Figure 3G), for which we also included additional quantification (Figure 3I, Figure 3 —figure supplement 2H). Moreover, as discussed above, we performed additional numerical simulations to mimic the effect of Lapatinib and show that in absence of memory, the system directly transits from the polarized to the basal state (Figure 2—figure supplement 2D), as experimentally identified in the reconstructed state-space trajectory.

– It would be helpful if one did not have to refer to the caption to read several of the figure panels. (An example: the color-coding in 3A, D).

We hope that the amended version has improved readability.

Reviewer #3 (Recommendations for the authors):

Timescale of memory

In my understanding, the lifetime of a 'ghost' state near a saddle-node bifurcation goes as 1/r^(1/2), where r is the distance to the critical point. This suggests that the timescale associated with the memory state is sensitive to how close the system is to the critical point, and, at least formally, diverges as the system actually approaches this point. This would seem to present a fine-tuning problem: not only does the system need to have parameters tuned to be near criticality, but in fact they have to be exactly the right distance away to achieve a physiologically reasonable intermediate memory timescale. It would be useful for the authors to discuss this: how is the memory timescale controlled? How sensitive is it to parameter changes? Is this somehow a feature, in that the cell can potentially physiologically change its memory timescale?

Indeed, plotting the EGFRmCitine intensity at the plasma membrane vs. the total memory duration in single-cell EGFRmCitrine phosphorylation polarization shows that higher the receptor duration, higher the memory length (Figure 6 in the response). Experimentally, we find a broad distribution of memory duration both, in polarized EGFRmCitrine phosphorylation in MCF7-EGFRmCitrine cells (Figure 2E), as well as memory in directional migration in MCF10A cells (Figure 3 —figure supplement 3A), suggesting that cell-to-cell variability in EGFR concentration generates a span of memory durations.

Due to the characteristics of the EGFR system – ligand-bound receptors are unidirectionally removed from the PM and undergo degradation, whereas ligandless receptor recycle, the total concentration of receptors on the PM will progressively decrease and thereby the system will move away from criticality. However, how long cells will maintain memory in directional migration and thereby organization at criticality depends on the duration and the concentration of the growth factor signal. We have extended the discussion on these points in the manuscript, and also suggest that it would be interesting to study whether receptor networks are self-organized at criticality, or the system evolved to this critical point as a means to have optimal information processing capabilities.

Suggestions regarding the theoretical exposition:

1. The authors should clarify the dimensional reduction that led to the bifurcation diagrams in Figure 1A and Figure 1—figure supplement 1B. Absent additional justification of this reduction, it would also be clearer to describe this as an approximate treatment of the system (whose behavior is born out by the reaction-diffusion simulations), rather than a proof (line 100).

In the revised version of the manuscript, we have added the justification of the one-dimensional projection, as outlined in the response to the first comment by the referee above.

2. It would be useful to readers to include a more concrete discussion of the EGFR model in the main text, including which features of it drive the behavior of the system and which biochemical parameters control location in the phase space. Additionally, how does the magnitude of the external gradient affect the cell? As a suggestion, the authors could consider moving Methods 5.15, Equation 17, to the main text, with a description of what the important variables are, and what features are important to the presence of the pitchfork bifurcation.

We thank the referee for this valuable suggestion. However, due to the extent of the EGFR model/parameter descriptions, we have decided to maintain the equations in the Methods section, but we have extended the description of the model and the dynamical transitions crucial to explain the mechanism of cellular navigation we propose here.

Suggestions regarding the presentation of the measurements:

1. Related to the public review, the authors could strengthen the paper by explicitly discussing the discrepancies between the measurements and the expectations from the model, and potential explanations. Does experimental noise interfere with the EGFR phosphorylation profile measurements? Do the authors believe that some of the cells are not at criticality?

In the amended version of the manuscript, we discuss in greater details the discrepancy between the polarization kinetics and amplitude of polarization in the experiments and simulations, the influence of biochemical and/or technical noise for phase space reconstruction, as well as the main differences between the migration trajectories in Figure 4. We underline however that the numerical simulations predict qualitatively the dynamics and behavior of the system.

2. The fact that the model depends on cells being biochemically poised near a critical point suggests a variety of stronger experimental tests of the framework. As one example, the authors' analysis suggests that overexpressing EGFR should push the cells away from the critical point, into the inhomogeneous steady-state regime, where they would break symmetry according to the first-encountered gradient and no longer be capable of adapting to a new gradient. This would be quite surprising, as it would correspond to breaking a sensing system by increasing the number of sensors. While performing this experiment is likely beyond the scope of this paper, the authors could strengthen the presentation by discussing this and/or other more counterintuitive predictions of their model in light of existing empirical data and/or future experiments.

As the referee suggested, for higher EGFR concentration on the membrane, the system will generate a permanent memory of the direction of the first encountered signal, which will render the cell unresponsive to further changes in signal localization (Figure 4 —figure supplement 1C,D and Figure 4 – video 3). In the amended version of the manuscript, we extended the discussion on this case, and we also provide an additional prediction – how cells resolve simultaneous signals. The numerical simulations predict that the dynamical memory that emerges from organization at criticality enables cells to compare and thereby resolve simultaneous signals with different amplitudes from different spatial localizations, whereas this feature is lost for organization in the IHSS regime (Figure 4 —figure supplement 2).

3. Figure 2—figure supplement 1C: given that the direction of the polarization relative to the gradient is important, it would be interesting to see all the polarization profiles (and the variability from cell to cell with respect to the direction relative to the gradient).

In the amended Figure 2 —figure supplement 1F, we provide a summary of the direction of polarization with respect to the gradient direction.

References:

1. Golubitsky, M., & Schaeffer, D. G. (1985). Singularities and Groups in Bifurcation Theory. 51. https://doi.org/10.1007/978-1-4612-5034-0

2. Koseska, A., Volkov, E., & Kurths, J. (2010). Parameter mismatches and oscillation death in coupled oscillators. Chaos: An Interdisciplinary Journal of Nonlinear Science, 20(2), 023132. https://doi.org/10.1063/1.3456937

3. Stanoev, A., Mhamane, A., Schuermann, K. C., Grecco, H. E., Stallaert, W., Baumdick, M., Brüggemann, Y., Joshi, M. S., Roda-Navarro, P., Fengler, S., Stockert, R., Roßmannek, L., Luig, J., Koseska, A., & Bastiaens, P. I. H. (2018). Interdependence between EGFR and Phosphatases Spatially Established by Vesicular Dynamics Generates a Growth Factor Sensing and Responding Network. Cell Systems, 7(3), 295-309.e11. https://doi.org/10.1016/J.CELS.2018.06.006

4. Reynolds, A. R., Tischer, C., Verveer, P. J., Rocks, O., & Bastiaens, P. I. H. (2003). EGFR activation coupled to inhibition of tyrosine phosphatases causes lateral signal propagation. Nature Cell Biology, 5(5), 447–453. https://doi.org/10.1038/NCB981

5. Baumdick, M., Brüggemann, Y., Schmick, M., Xouri, G., Sabet, O., Davis, L., Chin, J. W., & Bastiaens, P. I. H. (2015). EGF-dependent re-routing of vesicular recycling switches spontaneous phosphorylation suppression to EGFR signaling. eLife, 4. https://doi.org/10.7554/ELIFE.12223

6. Baumdick, M., Gelléri, M., Uttamapinant, C., Beránek, V., Chin, J. W., & Bastiaens, P. I. H. (2018). A conformational sensor based on genetic code expansion reveals an autocatalytic component in EGFR activation. Nature Communications 2018 9:1, 9(1), 1–13. https://doi.org/10.1038/s41467-018-06299-7

7. Wood, E. R., Truesdale, A. T., McDonald, O. B., Yuan, D., Hassell, A., Dickerson, S. H., Ellis, B., Pennisi, C., Horne, E., Lackey, K., Alligood, K. J., Rusnak, D. W., Gilmer, T. M., & Shewchuk, L. (2004). A unique structure for epidermal growth factor receptor bound to GW572016 (Lapatinib): relationships among protein conformation, inhibitor off-rate, and receptor activity in tumor cells. Cancer Research, 64(18), 6652–6659. https://doi.org/10.1158/0008-5472.CAN-04-1168

8. Skoge, M., Yue, H., Erickstad, M., Bae, A., Levine, H., Groisman, A., Loomis, W. F., & Rappel, W. J. (2014). Cellular memory in eukaryotic chemotaxis. Proceedings of the National Academy of Sciences of the United States of America, 111(40), 14448–14453. https://doi.org/10.1073/PNAS.1412197111

9. Stanoev, A., Nandan, A. P., & Koseska, A. (2020). Organization at criticality enables processing of time-varying signals by receptor networks. Molecular Systems Biology, 16(2), e8870. https://doi.org/10.15252/MSB.20198870

10. Prentice-Mott, H. v., Meroz, Y., Carlson, A., Levine, M. A., Davidson, M. W., Irimia, D., Charras, G. T., Mahadevan, L., & Shah, J. v. (2016). Directional memory arises from long-lived cytoskeletal asymmetries in polarized chemotactic cells. Proceedings of the National Academy of Sciences of the United States of America, 113(5), 1267–1272. https://doi.org/10.1073/PNAS.1513289113/SUPPL_FILE/PNAS.1513289113.SM15.MOV

11. Goryachev, A. B., & Pokhilko, A. v. (2008). Dynamics of Cdc42 network embodies a Turing-type mechanism of yeast cell polarity. FEBS Letters, 582(10), 1437–1443. https://doi.org/10.1016/J.FEBSLET.2008.03.029

12. Otsuji, M., Ishihara, S., Co, C., Kaibuchi, K., Mochizuki, A., & Kuroda, S. (2007). A Mass Conserved Reaction–Diffusion System Captures Properties of Cell Polarity. PLOS Computational Biology, 3(6), e108. https://doi.org/10.1371/JOURNAL.PCBI.0030108

13. Parent, C. A., & Devreotes, P. N. (1999). A cell’s sense of direction. Science, 284(5415), 765–769. https://doi.org/10.1126/SCIENCE.284.5415.765

14. Mori, Y., Jilkine, A., & Edelstein-Keshet, L. (2008). Wave-Pinning and Cell Polarity from a Bistable Reaction-Diffusion System. Biophysical Journal, 94(9), 3684. https://doi.org/10.1529/BIOPHYSJ.107.120824

15. Koseska, A., & Bastiaens, P. I. H. (2020). Processing Temporal Growth Factor Patterns by an Epidermal Growth Factor Receptor Network Dynamically Established in Space. Https://Doi.Org/10.1146/Annurev-Cellbio-013020-103810, 36, 359–383. https://doi.org/10.1146/ANNUREV-CELLBIO-013020-103810

16. Ibach J., et al.,(2015). Single particle tracking reveals that EGFR signaling activity is amplified in clathrin coated pits. PloS One 10, e0143162.

17. Stoscheck, C. M., & Carpenter, G. (1984). Characterization of the metabolic turnover of epidermal growth factor receptor protein in A-431 cells. Journal of Cellular Physiology, 120(3), 296–302. https://doi.org/10.1002/JCP.1041200306

18. Fischer, E. H., Charbonneau, H., & Tonks, N. K. (1991). Protein tyrosine phosphatases: a diverse family of intracellular and transmembrane enzymes. Science, 253(5018), 401–406. https://doi.org/10.1126/SCIENCE.1650499

19. Jilkine, A., Edelstein-Keshet, L., (2011). A comparison of mathematical models for polarization of single eukaryotic cells in response to guided cues. PloS Computational Biology. https://doi.org/10.1371/journal.pcbi.1001121

https://doi.org/10.7554/eLife.76825.sa2

Article and author information

Author details

  1. Akhilesh Nandan

    Department of Systemic Cell Biology, Max Planck Institute of Molecular Physiology, Dortmund, Germany
    Present address
    Cellular Computations and Learning, Max Planck Institute for Neurobiology of Behavior - caesar, Bonn, Germany
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review and editing
    Contributed equally with
    Abhishek Das
    Competing interests
    No competing interests declared
  2. Abhishek Das

    Department of Systemic Cell Biology, Max Planck Institute of Molecular Physiology, Dortmund, Germany
    Present address
    Cellular Computations and Learning, Max Planck Institute for Neurobiology of Behavior - caesar, Bonn, Germany
    Contribution
    Data curation, Formal analysis, Investigation, Software, Validation, Visualization, Writing – review and editing
    Contributed equally with
    Akhilesh Nandan
    Competing interests
    No competing interests declared
  3. Robert Lott

    Department of Systemic Cell Biology, Max Planck Institute of Molecular Physiology, Dortmund, Germany
    Present address
    Cellular Computations and Learning, Max Planck Institute for Neurobiology of Behavior - caesar, Bonn, Germany
    Contribution
    Formal analysis, Investigation, Software, Validation, Visualization, Writing – review and editing
    Competing interests
    No competing interests declared
  4. Aneta Koseska

    Department of Systemic Cell Biology, Max Planck Institute of Molecular Physiology, Dortmund, Germany
    Present address
    Cellular Computations and Learning, Max Planck Institute for Neurobiology of Behavior - caesar, Bonn, Germany
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Writing - original draft, Writing – review and editing
    For correspondence
    aneta.koseska@mpinb.mpg.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4263-2340

Funding

Max Planck Society

  • Aneta Koseska

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors thank Angel Stanoev for initial analysis of the EGFR polarity model, critical discussion during the project, as well as valuable comments on the manuscript, Manish Yadav and Monika Scholz for critical reading of the manuscript and valuable suggestions, and Frédéric Paquin-Lefebvre for valuable suggestions on the realization of the reaction-diffusion simulations. All of the experiments were carried in the lab of Philippe Bastiaens, and we are particularly grateful for the opportunity to be part of that engaging and critical community where we could learn and develop this project. We especially thank P Bastiaens for numerous critical discussion and suggestions that were crucial throughout the project, as well as for detailed comments that significantly helped us to improve the manuscript. Data and code availability: All data generated or analysed during this study is included in the manuscript and supporting files. The codes for the numerical simulations are available under https://github.com/akhileshpnn/Cell-memory. Funding: The project was funded by the Max Planck Society, partially through the Lise Meitner Excellence Program.

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Arvind Murugan, University of Chicago, United States

Reviewer

  1. Elizabeth R Jerison, Stanford University, United States

Publication history

  1. Preprint posted: November 12, 2021 (view preprint)
  2. Received: January 6, 2022
  3. Accepted: June 3, 2022
  4. Accepted Manuscript published: June 6, 2022 (version 1)
  5. Version of Record published: July 14, 2022 (version 2)

Copyright

© 2022, Nandan, Das et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 977
    Page views
  • 369
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Akhilesh Nandan
  2. Abhishek Das
  3. Robert Lott
  4. Aneta Koseska
(2022)
Cells use molecular working memory to navigate in changing chemoattractant fields
eLife 11:e76825.
https://doi.org/10.7554/eLife.76825

Further reading

    1. Physics of Living Systems
    Ye Li, Shiqi Liu ... Yilin Wu
    Research Article

    Long-range material transport is essential to maintain the physiological functions of multicellular organisms such as animals and plants. By contrast, material transport in bacteria is often short-ranged and limited by diffusion. Here we report a unique form of actively regulated long-range directed material transport in structured bacterial communities. Using Pseudomonas aeruginosa colonies as a model system, we discover that a large-scale and temporally evolving open channel system spontaneously develops in the colony via shear-induced banding. Fluid flows in the open channels support high-speed (up to 450 µm/s) transport of cells and outer membrane vesicles over centimeters, and help to eradicate colonies of a competing species Staphylococcus aureus. The open channels are reminiscent of human-made canals for cargo transport, and the channel flows are driven by interfacial tension mediated by cell-secreted biosurfactants. The spatial-temporal dynamics of fluid flows in the open channels are qualitatively described by flow profile measurement and mathematical modeling. Our findings demonstrate that mechanochemical coupling between interfacial force and biosurfactant kinetics can coordinate large-scale material transport in primitive life forms, suggesting a new principle to engineer self-organized microbial communities.

    1. Physics of Living Systems
    Nicola Rigolli, Gautam Reddy ... Massimo Vergassola
    Research Article Updated

    Foraging mammals exhibit a familiar yet poorly characterized phenomenon, ‘alternation’, a pause to sniff in the air preceded by the animal rearing on its hind legs or raising its head. Rodents spontaneously alternate in the presence of airflow, suggesting that alternation serves an important role during plume-tracking. To test this hypothesis, we combine fully resolved simulations of turbulent odor transport and Bellman optimization methods for decision-making under partial observability. We show that an agent trained to minimize search time in a realistic odor plume exhibits extensive alternation together with the characteristic cast-and-surge behavior observed in insects. Alternation is linked with casting and occurs more frequently far downwind of the source, where the likelihood of detecting airborne cues is higher relative to ground cues. Casting and alternation emerge as complementary tools for effective exploration with sparse cues. A model based on marginal value theory captures the interplay between casting, surging, and alternation.