Cells use molecular working memory to navigate in changing chemoattractant fields
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. Timeresolved livecell imaging of Epidermal growth factor receptor (EGFR) phosphorylation dynamics shows that cells transiently memorize position of encountered signals via slowescaping remnant of the polarized signaling state, a dynamical ‘ghost’, driving memoryguided persistent directional migration. The metastability of this state further enables migrational adaptation when encountering new signals. We thus identify basic mechanism of realtime 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.sa0eLife 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 shortlived 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 EdelsteinKeshet, 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 feedforward, excitable or Turinglike networks have been proposed to describe how polarized signaling activity of cellsurface 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 polarizedsignaling steady state in presence of guiding external cues. However, they can account either for sensing and adaptation to nonstationary stimuli or for longterm maintenance of polarized signaling activity, but not both. Thus, how cells process the information from a changing chemoattractant field in real time for longrange 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 livecell 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 realtime cellular navigation in complex chemoattractant fields.
Results
Dynamical mechanism of navigation in nonstationary environments
We conjectured that only dynamically metastable receptor signaling states can enable both transient stability of polarized signaling as necessary for robust, memoryguided 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 signalinduced transition from a nonpolarized symmetric, to a polarized receptor signaling state, and subsequently polarized cell shape. This transition is thus a symmetrybreaking 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 saddlenode ($SN$) bifurcation, that characterizes a transition between stable and unstable steady states. When the $SN$ and thereby a stable steadystate 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 subcritical $PB$, as it is stabilized via $S{N}_{PB}$ s. We thus propose that organization at criticality  in the vicinity of a $S{N}_{PB}$ (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.
We described this conjecture mathematically for a general reactiondiffusion model representing the signaling activity on the plasma membrane of a cell, $\frac{\partial \mathbf{U}(\mathbf{x},\mathbf{t})}{\partial t}=\mathbf{F}(\mathbf{U})+\mathbf{D}{\nabla}^{2}\mathbf{U}(\mathbf{x},\mathbf{t})$, with $\mathbf{U}$ being the vector of local densities of active signaling components, $\mathbf{D}$  diffusion constants and $\mathbf{F}$ accounting for all chemical reactions. Our theoretical analysis shows that a $PB$ exists if, for a spatial perturbation of the symmetric steady state (${\mathbf{U}}_{\mathbf{s}}$) of the form $\mathbf{U}(\mathbf{x},t)={\mathbf{U}}_{\mathbf{s}}+\delta \mathbf{U}(\mathbf{x}){e}^{\lambda t}$, the conditions $\delta \mathbf{U}(\mathbf{x})=\delta \mathbf{U}(\mathbf{x})$ and the limit ${lim}_{\lambda \to 0}{F}_{\lambda}=det(J)=0$ are simultaneously fulfilled (Materials and methods). This implies that the linearized system has zerocrossing eigenvalues ($\lambda $) associated with the odd mode of the perturbation (PaquinLefebvre et al., 2020). To probe the subcritical transition and therefore the necessary organization at criticality, a reduced description in terms of an asymptotic expansion of the amplitude of the polarized state ($\varphi $) must yield the Landau equation $\frac{d\varphi}{dt}={c}_{1}\varphi +{c}_{2}{\varphi}^{3}{c}_{3}{\varphi}^{5}$, guaranteeing the existence of $S{N}_{PB}$ (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 abovementioned principle, we use the wellcharacterized 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 (${E}_{p}$) with two enzymes, the phosphatases PTPRG (${P}_{RG}$) and PTPN2 (${P}_{N2}$, Figure 1B, Figure 1—figure supplement 1B), respectively. ${E}_{p}$ and ${P}_{RG}$ laterally diffuse on the membrane and inhibit eachother’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 symmetrybreaking 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 signaldependent cell shape changes and subsequent migration (ChiassonMacKenzie 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 ligandbound receptors ($E{E}_{p}$) 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 (statespace 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 statespace, 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 timescales 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 timescale 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 longterm 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 finetuning 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 $S{N}_{PB}$, 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 cancerderived MCF7 cells were subjected for 1 hr to a stable gradient of fluorescently tagged EGFAlexa647 (EGF^{647}) with a maximal amplitude of 10 ng/ml applied from the top of the chamber in a computerprogrammable microfluidic device (Figure 2A and B). EGFR phosphorylation at the plasma membrane was quantified during and for 3 hr after gradient washout (gradient washout established in 4–5 min) by determining the rapid translocation of mCherrytagged phosphotyrosinebinding domain (PTB^{mCherry}) to phosphorylated tyrosines 1086/1148 of ectopically expressed EGFRmCitrine (EGFR^{mCitrine}) using ratiometric imaging (Offterdinger et al., 2004; Materials and methods). Due to the low endogenous EGFR levels in MCF7 cells, the expression range of EGFR^{mCitrine} 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.
Kymograph analysis of $EGFR}^{mCitrine$ phosphorylation at the plasma membrane of single cells showed polarization in a shallow gradient of EGF^{647} (as shallow as 10% between front and back of the cell; Figure 2C, Figure 2—figure supplement 1AD). The direction of $EGFR}^{mCitrine$ phosphorylation polarization coincided with the direction of maximal EGF^{647} concentration around each cell ($\pi /4$ on average, Figure 2—figure supplement 1F). Only few cells manifested basal or symmetric $EGFR}^{mCitrine$ phosphorylation distribution upon gradient stimulation (Figure 2—figure supplement 1A, B, E). Plotting the fraction of plasma membrane area with polarized $EGFR}^{mCitrine$ phosphorylation showed celltocell variability in the polarization kinetics, as well as the maximal amplitude of polarized $EGFR}^{mCitrine$ 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 $EGFR}^{mCitirne$ phosphorylation was also reflected in the respective singlecell temporal profiles (exemplary profile shown in Figure 2F, top). Reconstructing the statespace 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 $EGFR}^{mCitrine$ phosphorylation profile, and thereby the reconstruction of the statespace 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 washout, the $EGFR}^{mCitrine$ phosphorylation response exponentially decayed, resulting in a clear absence of transient memory and respective statespace 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 statespace trajectory, where the system directly transits from the polarized to the basal state, without intermediate statespace trapping (Figure 2—figure supplement 2D, Figure 2—video 3).
Fitting the experimentally measured singlecell temporal $EGFR}^{mCitrine$ phosphorylation profiles after gradient washout using an inverse sigmoid function (Methods) further corroborated that under Lapatinib treatment, phosphorylated $EGFR}^{mCitrine$ exponentially relaxed from the polarized to the basal state (Hill coefficient ≈ 1.28), with a halflife of approx. 10 min (Figure 2—figure supplement 2E, G). Under normal conditions however, the halflife was 30 min on average, reflecting that the phosphorylated $EGFR}^{mCitrine$ 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 $EGFR}^{mCitrine$ phosphorylation results from a dynamically metastable ‘ghost’ state, and not a slow dephosphorylation process.
In order to identify whether the memory in polarized $EGFR}^{mCitrine$ 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 singlecell 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 $EGFR}^{mCitrine$ 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 $EGFR}^{mCitrine$ 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 $MCF7EGFR}^{mCitrine$, as well as of MCF10A cells at physiological EGF concentrations. Cells were subjected to a 5 hr dynamic EGF^{647} gradient that was linearly distributed within the chamber, with EGF^{647} ranging between 250 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 singlecell’s motility trajectories was performed for 14 hr in total. $MCF7EGFR}^{mCitrine$, as well as MCF10A cells migrated in a directional manner toward the EGF^{647} 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 washout (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 EGF^{647} 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 uniformstimulation case (Figure 3B). Moreover, the directionality estimated in the 9 hr timeframe 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.
This was also reflected in the projection of the cell’s relative displacement angles ($\mathrm{cos}\theta $) estimated along the gradient direction ($\pi $) 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 ($\mathrm{cos}\theta $ approached 1) during, and maintained this temporally after gradient removal, before returning to a migration pattern characteristic for stimulus absence or during uniform stimulation ($\mathrm{cos}\theta \approx $ ~ 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 timeframe in which the memory in polarized $EGFR}^{mCitrine$ 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 OrnsteinUhlenbeck process (Uhlenbeck and Ornstein, 1930; Svensson et al., 2017) and used the extracted migration parameters to generate synthetic singlecell 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 washout (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 washout due to the absence of memory in polarized $EGFR}^{mCitrine$ phosphorylation (Figure 2G,I). Thus, singlecell 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 $\mathrm{cos}\theta $ 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 celltocell variability in this case, we also calculated memory duration from single cell $\mathrm{cos}\theta $ profiles. For this, singlecell 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 (PrenticeMott et al., 2016). Without Lapatinib treatment however, the duration of memory estimated from singlecell $\mathrm{cos}\theta $ 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 singlecell 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).
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 longterm 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 longterm, 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 memoryguided migration over long trajectories.
We next tested these predictions experimentally by establishing an equivalent dynamic EGF^{647} spatialtemporal 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 EGF^{647} absence (w/o EGF^{647}, distribution symmetrically distributed around $\mathrm{cos}\theta =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 $\mathrm{cos}\theta $ values, reflecting that cells revert the direction of migration (established in ∼10 min). Furthermore, the reverse migration was maintained for approximately 20 min after washout 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 EGF^{647}). These results therefore demonstrate that cells utilize molecular working memory to navigate in changing gradient fields.
Navigation in nonstationary 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 longterm 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 stimulusresponse 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 steadystate 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 slowescaping 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 (PrenticeMott et al., 2016), the memory in receptor signaling we identify here provides a crucial bridge between the rapid receptor phosphorylation/dephosphorylation events and the longrange cellular migration. In particular, the organization at criticality endows the system with a slow timescale 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 celltocell 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 selforganized at criticality through an active sensing mechanism, or this feature has been finetuned 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 memorydirected migration. Such memoryguided 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
Cell culture
Request a detailed protocolMCF7 cells (sex: female, ECACC, Cat. No. 86012803) were grown at 37 °C and 5% CO_{2} in Dulbecco’s Eagle’s medium (DMEM) (PANBiotech, Germany), supplemented with 10% inactivated Fetal Calf Serum (FCS) (SigmaAldrich), 100 ng ml^{–1} LGlutamine, 0.5 mg ml^{–1} nonessential amino acids, 100 µg ml^{1} penicillin and 100 µg ml^{1} streptomycin (PANBiotech, 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 (PANBiotech, Germany). MCF10A cells (sex: female, ATCCCRL 10317) were grown at 37 °C and 5% CO_{2} in Mammary Epithelial Cell Growth Basal medium (MEBM from Lonza Pharma & Biotech), supplemented with 5% Horse Serum (HS) (Invitrogen), 20 ng ml^{–1} EGF (SigmaAldrich), 0.5 mg ml^{–1} hydrocortisone (SigmaAldrich), 100 ng ml^{–1} cholera toxin (SigmaAldrich), 10 µg ml^{1}insulin (SigmaAldrich), 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 (SigmaAldrich), 100 ng ml^{–1} cholera toxin (SigmaAldrich), 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 (LeibnizInstitut DSMZ). Cells were regularly tested for mycoplasma contamination using MycoAlert Mycoplasma detection kit (Lonza).
Transfection and cell seeding
Request a detailed protocolFor $EGFR}^{mCitrine$ polarization experiments, 2.5 × 10^{5} MCF7 cells were seeded per well in a sixwell LabTek 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 ($EGF{R}^{mCitrine}$, $PT{B}^{mCherry}$ and $cCb{l}^{BFP}$ at ratio 4:3:4 by mass) using FUGENE6 (Roche Diagnostics) transfection reagent and OptiMEM (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 (SigmaAldrich). After $10\text{}\mathrm{min}$ incubation period at 37 °C and 5% CO_{2}, fresh growth media was added, and the cell density and viability was measured using cell counter (ViCELL XR Cell Viability Analyzer System). After spinning down, the cells were diluted to 10 × 10^{6} cells/ml. The M04G02 microfluidic gradient plates (Merck Chemicals) were primed for usage by flowing cell culture growth media through the cell chamber for $5\text{}\mathrm{min}$ and cells were subsequently seeded according to manufacturer’s instructions.
For migration experiments with uniform $EG{F}^{647}$ stimulation, 6well LabTek plates were coated with Collagen (SigmaAldrich) in 0.1 M Acetic acid (SigmaAldrich) for MCF7 (100 µg cm^{2}), and Fibronectin (SigmaAldrich) in PhosphateBuffered Saline (DPBS) (PANBiotech) 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 × 10^{5} cells per well were used for seeding. For migration experiments with gradient EGF^{647} stimulation, MCF7 cells were transferred to the coated M04G02 microfluidic gradient plates as described above. Before seeding, MCF10A cells were detached from 6 well LabTeks by discarding the growth media and washing once with DPBS (PAN Biotech) before adding 100 µl Accutase (SigmaAldrich). After 20–30 min incubation period at 37 °C and 5% CO_{2}, fresh cell growth media was added, and the cell density and viability were measured using a cell counter (ViCELL XR Cell Viability Analyzer System). After spinning down, the cells were diluted to 2 × 10^{6} 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 $25\text{}\mathrm{m}\mathrm{\text{M}}$ 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 widefield microscopy
Request a detailed protocolConfocal images were recorded using a Leica TCS SP8i confocal microscope (Leica Microsystems) with an environmentcontrolled 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 $405\text{}\mathrm{nm}$ diode laser. The detection of fluorescence emission was restricted with an AcoustoOptical 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 12bit 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 C0_{2} incubation chamber at 37 °C and 5% CO_{2}. 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 protocolThe CellAsic Onix Microfluidic Platform (EMD Millipore) was used for gradient cell migration and $EGFR}^{mCitrine$ phosphorylation polarization experiments. For $EGFR}^{mCitrine$ phosphorylation polarization experiments, 1 hr gradient stimulation was established using CellASIC ONIX2 software as follows. (i) Prestimulus: 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, preloaded EGF^{647} (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} EGF^{647} 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} EGF^{647}/ 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, EGF^{647} (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, EGF^{647} (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 $285\text{}\mathrm{min}$.
Imaging $EGF{R}^{mCitrine}$ phosphorylation polarization and single cell migration
Request a detailed protocolTransfected $MCF7EGFR}^{mCitrine$ cells transferred to M04G02 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 $1\text{}\mathrm{min}$ time interval using adaptive autofocus system and the water objective was performed concurrently during the duration of the experiment using the Leica TCS SP8i.
For migration experiments under uniform EGF^{647} stimulation, confocal laser scanning microscopy / transmission imaging of live $MCF7EGFR}^{mCitrine$ / 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.
EGF^{647} / Fluorescein gradient quantification
Request a detailed protocolhEGF^{647} was generated in the lab of Prof. P. Bastiaens, MPI of molecular Physiology, Dortmund, using the HisCBDIntein(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), Nterminally labeled with Alexa647maleimide as described previously (Sonntag et al., 2014) and stored in PBS at –20 °C. To quantify the spatial extent of the EGF^{647} /Fluorescein gradient, gradients were generated following the protocol described in subsection 5.6 in plates without cells or matrix coating. Confocal images of Alexa647 /GFP channel were acquired at $1\text{}\mathrm{min}$ 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 EGF^{647} /Fluorescein concentration.
Quantifying $EGFR}^{mCitrine$ phosphorylation in single cells
Request a detailed protocolTo quantify plasma membrane $EGFR}^{mCitrine$ phosphorylation in live $MCF7EGFR}^{mCitrine$ cells, singlecell masks were obtained from the $EGFR}^{mCitrine$ channel at each timepoint 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 anticlockwise direction, and each was divided into 5 spatial bins (Figure 2A). The fraction of phosphorylated $EGFR}^{mCitrine$ in each bin, $i$ was estimated as:
where $PT{B}_{PM}^{i}(t)$ and $EGF{R}_{PM}^{i}(t)$ are respectively the $PTB}^{mCherry$ and $EGF{R}^{mCitrine}$ fluorescence at ${i}^{th}$ plasma membrane bin, $PT{B}_{T}(t)$ and $EGF{R}_{T}(t)$  respective total fluorescence in the whole cell, $PT{B}_{endo}(t)$ – the $PTB}^{mCherry$ fluorescence on vesicular structures in the cytoplasm. Endosomal structures were identified from the cytosol by intensity thresholding (1.5 s.d. percentile) and $PTB}^{mCherry$ fluorescence from these structures was subtracted from the $PT{B}_{T}(t)$, to correct for the $PTB}^{mCherry$ fraction bound to the phosphorylated $EGFR}^{mCitrine$ on endosomes.
Temporal profile of the fraction of phosphorylated $EGFR}^{mCitrine$ on the plasma membrane was obtained using:
and then normalized as:
with <> being the temporal average in the prestimulation interval $t\in [0,5min]$. The fraction of liganded receptor was calculated using:
To classify single cells into nonactivated, activated (polarized $EGFR}^{mCitrine$ phosphorylation) and preactivated (uniformly distributed $EGFR}^{mCitrine$ phosphorylation) upon gradient EGF^{647} stimulation (Figure 2—figure supplement 2A, B), the following method was applied. To identify preactivated cells, a Gaussian Mixture Model (GMM) was fitted to the histogram of ${(EGF{R}_{p}^{i})}_{t\in [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 ${(EGF{R}_{p}^{i})}_{t\in [0,5min]}$ pixel intensity values for any cell lie above the intersection point, the cell is classified as preactivated. To distinguish between the nonactivated and activated cells in the remaining population, average $EGFR}^{mCitrine$ phosphorylation value ($EGF{R}_{p}$) per cell was estimated during the prestimulation ($t\in [0,5min]$) and the stimulation period ($t\in [5min,65min]$) ($<EGF{R}_{p}{>}_{t\in [0,65]}$) from the temporal $EGFR}^{mCitrine$ phosphorylation profiles. Histogram of the respective $EGF{R}_{p}$ values was again fitted with a GMM model. All cells with an average $<EGF{R}_{p}{>}_{t\in [0,65]}$ value lying below the intersection point were considered to be nonactivated, whereas those above  activated.
The average of the spatial projection of the fraction of phosphorylated $EGFR}^{mCitrine$ from singlecell kymographs (Figure 2—figure supplement 1C) was generated from the 20 cells that were polarized in the direction of the EGF^{647} gradient. For each cell, a temporal average of $EGF{R}_{p}$ per bin was calculated for the duration of the gradient ($t\in [5min,65min]$) and the bin with the maximal $EGF{R}_{p}$ value was translated to $\pi $. 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 EGF^{647} distribution around single cells (Figure 2—figure supplement 1F) was estimated as follows: the cell mask obtained using the $EGFR}^{mCitrine$ 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 EGF^{647} outside single cells. This outer contour was divided in 20 bins as for the kymographs, and EGF^{647} intensity was quantified in each bin. The angle between the direction of EGF^{647} and the direction of EGFR phosphorylation was calculated as the amount of radial bins between the maxima in the spatial projections. This bindistance was then translated into an angle under the assumption of a circular perimeter.
In order to identify the characteristic features of the $EGFR}^{mCitrine$ phosphorylation profile during the transition from polarized to unpolarized state, the singlecell $EGF{R}_{p}(t)$ profiles with and without Lapatinib treatment after gradient washout were fitted to an inverse sigmoid function given by,
were a_{0}, $a$ are constants and $n$ is the Hillcoefficient (examples in Figure 2—figure supplement 2E, F). Nonlinear least square method (python package curve fit) was used to perform the fitting. Under normal conditions (w/o Lapatinib), $a\sim 10$, ${a}_{0}\sim {10}^{3}$ and $n\sim 2.88$ fitted well the data (${R}^{2}\sim 0.79$). The same function however could not describe the EGFRp profiles in the Lapatinib treatment experiment (median ${R}^{2}\sim 0.33$). The Lapatinib treatment profiles were therefore fitted by fixing $a=10$, and leaving a_{0} 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}_{0}\sim 19$ and $n\sim 1.28$ were identified from the fitting (median ${R}^{2}\sim 0.84$, Figure 2—figure supplement 2E, G). From the fitted profiles in both cases, halflife was estimated to be the time frame in which 50% of $EGFR}^{mCitrine$ phosphorylation is lost after EGF^{647} removal.
Estimating memory duration in $\mathrm{E}\mathrm{G}\mathrm{F}\mathrm{R}}^{mCitrine$ phosphorylation polarization
Request a detailed protocolThe duration of memory in $EGFR}^{mCitrine$ phosphorylation polarization in single cells was estimated from the temporal profile of the fraction of plasma membrane area with high $EGFR}^{mCitrine$ phosphorylation during and after gradient removal (Figure 2D and E). For this, the singlecell kymographs were normalized to a maximal value of 1 using
yielding the value of phosphorylated $EGFR}^{mCitrine$ per bin $i$ per time point $t$. Using the mean of $EGF{R}_{p}+s.d.$ over the whole experiment duration as a threshold, all $EGF{R}_{p}^{i}(t)$ lying above the threshold were taken to constitute the area of polarized $EGFR}^{mCitrine$ phosphorylation. To account for different bin sizes, at each timepoint, the area of all bins with $EGF{R}_{p}$ 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 $FP{A}_{percell}<(FP{A}_{average}s.d.)$ in 3 consecutive time points (Figure 2E).
Quantifying morphological changes in response to EGF^{647} in experiments and simulations
Request a detailed protocolMorphological 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 $\sigma $ is the ratio between the cell’s area ${A}_{cell}$ and the area of the convex hull ${A}_{convex}$ ($\sigma =\frac{{A}_{cell}}{{A}_{convex}}$). The memory duration in cell morphology was calculated from the singlecell solidity profiles, and corresponds to the timepoint at which the solidity is below means.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 µm^{2} 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 $\frac{{A}_{prot,front}}{{A}_{front}}\frac{{A}_{prot,back}}{{A}_{back}}$. The final profiles of directed protrusive area were smoothed using 1D Gaussian filtering with the $filters.gaussian\text{\_}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 singlecell migration and duration of memory in directed cell migration
Request a detailed protocolSingle 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 $\mathrm{cos}\theta $ 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 twosided Welch’s ttest. To quantify the memory duration in directed singlecell migration, the Kernel Density Estimate (KDE) from $\mathrm{cos}\theta $ quantification in the continuous absence of EGF^{647} (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 KolmogorovSmirnov test. To verify the absence of memory when cells were treated with Lapatinib during gradient washout, a moving window KDE (5 time points) from $\mathrm{cos}\theta $ obtained in this case was compared to the KDE in continuous absence of EGF^{647} (uniform case Figure 3—figure supplement 2B, between 250–300 min) using two sided KolmogorovSmirnov 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 EGF^{647}, confirming the rapid switch from directed to randomwalk 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 $EG{F}^{647}$ stimulation, we fitted the experimentally obtained single cell migration trajectories using modified OrnsteinUhlenbeck process (mOU) (Uhlenbeck and Ornstein, 1930) that is defined by the Langevin equation for the velocity vector $\nu $:
where $\xi (t)$ represents a white noise component, D is a diffusion coefficient characteristic of a Brownian motion, $\tau $ is the persistence time and $b(t)$ models the contribution of the timedependent bias. The experimental data was fitted to obtain values of D and $\tau $. In order to estimate D, Mean Square Displacement (MSD) was calculated from the single cell tracks using $MSD(t)=<{\mathbf{x}}_{i}(t){\mathbf{x}}_{i}(0){}^{2}>$, where ${\mathbf{x}}_{i}(t)$ is the tracked position of ith 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 EGF^{647}, ${R}^{2}=0.975$; for uniform 20 ng/ml EGF^{647} stimulation, ${R}^{2}=0.995$. In order to estimate $\tau $, Velocity AutoCorrelation Function $VACF(t)=<{\nu}_{i}(t)\cdot {\nu}_{i}(0)>$, where ${\nu}_{i}(t)$ is the measured velocity of $i$th cell at time t, was fitted with a mono exponential function ($={\varphi}_{0}\cdot {e}^{\frac{t}{\tau}}$). Goodness of Fit: for 0 ng/ml EGF^{647} case  Standard Error of Estimate $SEOE=0.0261$; for uniform 20 ng/ml EGF^{647} stimulation case, $SEOE=0.0570$. Fitted values: for 0 ng/ml EGF^{647} case, $\tau =11.105$, $D=0.425$; for uniform 20 ng/ml EGF^{647} stimulation case, $\tau =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 Kalmanfilter (python package filterpy.kalman) by predicting the cell position and velocity. The cell’s displacement angles relative to the gradient direction ($\mathrm{cos}\theta $) were calculated for each cell at each timepoint, rendering singlecell $\mathrm{cos}\theta $ plots (Figure 3—figure supplement 3B,C). The memory duration was then calculated as the point where three consecutive timepoints in the $\mathrm{cos}\theta $ profiles fall below a threshold $\mathrm{cos}\theta $ value of 0.75.
Reconstructing statespace trajectories from temporal $EGFR}^{mCitrine$ phosphorylation profiles
Request a detailed protocolThe statespace reconstruction in Figure 2F and G was performed using the method of timedelay. For a time series of a scalar variable, a vector $x({t}_{i})$, $i=1,...N$ in statespace in time t_{i} can be constructed as following
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({t}_{i})$ consists of scalar measurements of the state of a dynamical system, then under certain genericity assumptions, the time delay embedding provides a onetoone 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 nearestneighbours method), and the statespace reconstruction using the buildTakens function, all from the nonlinearTseries package in R (https://cran.rproject.org/web/packages/nonlinearTseries/index.html). Before statespace reconstructions, time series were smoothened using the SavitzkyGolay filter function in Python. For Figure 2F, $d=26$, ${d}_{e}=3$; for Figure 2G, $d=50$, ${d}_{e}=3$.
Theoretical consideration of the navigation mechanism in a generalized reactiondiffusion signaling model
Request a detailed protocolWe consider a generalized form of a (massconserved) reactiondiffusion (RD) model of an $M$ ($\mathbf{U}\in {\mathbf{R}}^{M}$) component system in $N$ ($\mathbf{x}\in {\mathbf{R}}^{N}$) dimensional space
where $\mathbf{F}\in {\mathbf{R}}^{M}$ is the reaction term, $\mathbf{D}$ is a $M\times M$ diagonal matrix of diffusion constants ${D}_{j},j=1,\mathrm{\dots},M$, and ${\nabla}^{2}$ is the Laplacian operator. Standard analysis of such models relies on linear stability analysis to find the conditions for a Turingtype 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 reactiondiffusion 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 ${\mathbf{U}}_{\mathbf{s}}=({u}_{is})$ for $i=1,...,M$, be the stable homogeneous symmetric steady state of the RD system. Consider a linear perturbation of the form
where $\delta \mathbf{U}(\mathbf{x})$ is the spatial and ${e}^{(\lambda 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}_{\lambda}=det(\lambda {I}_{M\times M}{J}_{M\times M})=0$. J is the Jacobian matrix of the system defined by ${J}_{ij}=\frac{\partial {\mathbf{F}}_{i}(\mathbf{U}(x,t))}{\partial {U}_{j}},i=1,\mathrm{\dots}.,M,j=1,\mathrm{\dots}..,M$.
The system exhibits a $PB$ if, an odd eigenfunction $\delta \mathbf{U}(\mathbf{x})$ such that $\delta \mathbf{U}(\mathbf{x})=\delta \mathbf{U}(\mathbf{x})$, taken in the limit $\lambda \to 0$, fulfills the following condition (PaquinLefebvre et al., 2020):
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 frontbackpolarized 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 subcritical type, and thereby identify the presence of an $S{N}_{PB}$, 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 onedimensional 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, k_{m} into a superposition of spatially periodic waves. That means that $u(x,t)\in \mathbf{U}$ can be written as:
where ${u}_{n}(t)$ is the complex amplitude of the ${n}^{th}$ harmonics. Let the amplitude corresponding to the leading harmonics $(n=1)$ is $\varphi (t)$. After assuming that the amplitude of every other harmonics can be written as a power series of $\varphi (t)$, substituting Equation 12 into Equation 9 allows to write an equation that describes the evolution of $\varphi (t)$. In the case when the resulting equation is of StuartLandau type:
with ${c}_{1},{c}_{2},{c}_{3}>0$, this corresponds to the normal form of a subcritical pitchfork bifurcation (Strogatz, 2018). Together with the condition given by Equation 11, the existence of a subcritical 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 $S{N}_{PB}$.
Modeling EGFR phosphorylation polarization dynamics
Request a detailed protocolThe dynamics of the experimentally identified spatially distributed EGFR sensing network (Figure 1B, Figure 1—figure supplement 1B) is described using the following onedimensional system of partial differential equations (PDEs):
with
and
The reaction terms are described in details in Stanoev et al., 2018. In brief, $[E{E}_{p}]$ is the phosphorylated ligandbound dimeric EGFR, $[{E}_{p}]$  ligandless phosphorylated EGFR, $[{E}_{t}]$  total amount of EGFR, $[R{G}_{a}],[R{G}_{t}]$ and $[N{2}_{a}],[N{2}_{t}]$  the active and total amount of the membrane localized PTPRG and the ERbound 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 ($[{E}_{t}]$, $[R{G}_{t}]$ and $[N{2}_{t}]$) are constant parameters. Autonomous, autocatalytic and ligandboundinduced activation of ligandless EGFR ensue from bimolecular interactions with distinct rate constants ${\alpha}_{13}$, respectively. Other parameters are as follows: ${k}_{1}/{k}_{2}$ — activation/inactivation rate constants of the phosphatases, ${\beta}_{1}/{\beta}_{2}$  receptorinduced regulation rate constants of $PTPRG/PTPN2$, ${\gamma}_{1}/{\gamma}_{2}$  specific reactivity of the enzymes ($PTPRG/PTPN2$) towards the receptor. The EGFRPTPN2 negative feedback is on a time scale ($\u03f5$) approximately two orders of magnitude slower than the phosphorylationdephosphorylation reaction, as estimated from the ~ 4 min recycling time of $EGF{R}_{p}$ (Stanoev et al., 2018). This enables, when necessary, to consider a quasisteady state approximation for the dynamics of PTPN2 for simplicity:
$[EG{F}_{t}]$ 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), ligandbinding unbinding was explicitly modeled (${k}_{on}$, ${k}_{off}$) 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 ERbound 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 timescale of EGF binding) (Ibach et al., 2015). In the reactiondiffusion (RD) simulations therefore for simplicity, it is assumed that ${D}_{E{E}_{p}}\approx 0$, whereas diffusion constants of same order are assumed for the ligandless EGFR and PTPRG (${D}_{{E}_{p}}\sim {D}_{R{G}_{a}}$).
Analytical consideration for an $S{N}_{PB}$ existence in the EGFR network
Request a detailed protocolTo identify analytically the existence of a $S{N}_{PB}$ 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 reactiondiffusion signaling model). For this, we considered the system Equation 14, where the dynamics of PTPN2 is at quasisteady state (Equation 15), $[E{E}_{p}]=0$, and rest of the dependent and independent variables were scaled to have a dimensionless form. Let $[\stackrel{~}{{E}_{p}}]=[{E}_{p}]/{E}_{0}$, $[\stackrel{~}{R{G}_{a}}]=[R{G}_{a}]/R{G}_{0}$, $\stackrel{~}{x}=x/{x}_{0}$, $\tau =t/{t}_{0}$, such that ${t}_{0}=1/({k}_{1}+{k}_{2})$, ${E}_{0}={k}_{1}/{\beta}_{2}$, $R{G}_{0}=({k}_{1}+{k}_{2})/{\gamma}_{1}$ and ${t}_{0}/{x}_{0}^{2}=1/{D}_{{E}_{p}}$. Substituting these into Equation 14 yields the system of dimensionless equations:
with $q}_{1}=\frac{{\alpha}_{1}\cdot [{E}_{t}{]}^{2}}{({k}_{1}+{k}_{2})\cdot {\beta}_{2}$, $q}_{2}=\frac{({\alpha}_{2}2\cdot {\alpha}_{1})\cdot [{E}_{t}]}{{k}_{1}+{k}_{2}$, $q}_{3}=\frac{({\alpha}_{1}{\alpha}_{2})\cdot {k}_{t}}{({k}_{1}+{k}_{2})\cdot {\beta}_{2}$, $q}_{4}=\frac{{\gamma}_{2}\cdot [N{2}_{t}]}{{k}_{1}+{k}_{2}$, $k={k}_{2}/{k}_{1}$, $r}_{1}=\frac{{k}_{1}\cdot [R{G}_{t}]\cdot {\gamma}_{1}}{({k}_{1}+{k}_{2}{)}^{2}$, $r}_{2}=\frac{{\beta}_{1}\cdot {k}_{1}}{({k}_{1}+{k}_{2})\cdot {\beta}_{2}$ and $D=\frac{{D}_{R{G}_{\alpha}}}{{D}_{{E}_{p}}}$.
We further simplify the system Equation 16 by taking the Talyor series expansion of the quasisteady state approximation of $[N{2}_{a}]$ around $E}_{s$, the steady state of $[\stackrel{~}{{E}_{p}]}$
with $q}_{7}=\frac{{E}_{s}{q}_{4}}{1+k+{E}_{s}}\frac{{E}_{s}{q}_{4}(1+k)}{(1+k+{E}_{s}{)}^{2}$, ${q}_{8}=\frac{{E}_{s}{q}_{4}}{1+k+{E}_{s}}+\frac{{q}_{4}(1+k)}{(1+k+{E}_{s}{)}^{2}}(1{E}_{s})$, and $q}_{9}=\frac{{q}_{4}(1+k)}{(1+k+{E}_{s}{)}^{2}$, thus yielding:
with $q}_{9}={q}_{1}{q}_{7},\phantom{\rule{thinmathspace}{0ex}}{q}_{10}={q}_{2}{q}_{8$ and $q}_{11}={q}_{3}{q}_{9$.
To avoid long expression in the further analysis, we rename the dependent variables as ${u}_{1}=[\stackrel{~}{{E}_{p}}]$ and ${u}_{2}=[\stackrel{~}{R{G}_{a}}]$, and the independent variables as $\stackrel{~}{x}=x$, $\tau =t$. The system Equation 16 therefore obtains the generic form:
In order to perform linear stability analysis, a onedimensional projection of Equation 19 is considered,
The simplified onedimensional 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 onedimensional projection, as demonstrated below, preserves all of the main features of the PDE model.
Let, ${\mathbf{U}}_{\mathbf{s}}=\left(\begin{array}{c}\hfill {u}_{1fs}\hfill \\ \hfill {u}_{2fs}\hfill \\ \hfill {u}_{1bs}\hfill \\ \hfill {u}_{2bs}\hfill \end{array}\right)$ be the stable symmetric steady state of the system (${u}_{1fs}={u}_{1bs}$, ${u}_{2fs}={u}_{2bs}$). A small amplitude perturbation on this symmetric steady state of the form,
yields a linearized equation,
where
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 onedimensional projection (Equation 20), the odd mode of the perturbation $(\delta \mathbf{U}(\mathbf{}\mathbf{x})=\delta \mathbf{U}(\mathbf{x}))$ must yield: $\delta {u}_{1f}=\delta {u}_{1b}$ and $\delta {u}_{2f}=\delta {u}_{2b}$. Substituting this into Equation 22 to obtain ${F}_{}(\lambda )$, in the limit $\lambda \to 0$ renders:
Thus, there exists parameter set for which existence of PB in the system Equation 20 is guaranteed.
To identify whether the PB is subcritical and thereby identify existence of a $S{N}_{PB}$, the solution of the system Equation 19 is approximated as in Equation 12:
The expansion is taken to $n={3}^{rd}$ order, rendering an amplitude equation of ${5}^{th}$ order. As described in Becherer et al., 2009, the complex coefficients of the $n={0}^{th},n={2}^{nd}$ and $n={3}^{rd}$ harmonics can be approximated as power series of $\varphi (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 $\varphi (t)$, ${u}_{0}(t)$, ${v}_{0}(t)$, ${u}_{1}(t)$, ${v}_{1}(t)$, ${u}_{2}(t)$, ${v}_{2}(t)$, ${u}_{3}(t)$ and ${v}_{3}(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 $\varphi $ and the parameters of the system as:
where ${u}_{2}^{(2)}=\frac{1{q}_{11}}{{q}_{10}4{k}_{m}^{2}}$, ${v}_{2}^{(2)}=\frac{{r}_{2}}{1+4D{k}_{m}^{2}}$, ${u}_{3}^{(3)}=\frac{{u}_{2}^{(2)}+{v}_{2}^{(2)}2{q}_{11}{u}_{2}^{(2)}}{{q}_{10}9{k}_{m}^{2}}$ and ${v}_{3}^{(3)}=\frac{{r}_{2}({u}_{2}^{(2)}+{v}_{2}^{(2)})}{1+9D{k}_{m}^{2}}$. The dynamics of the leading harmonics $(n=1)$ can be written as:
where ${c}_{1}={q}_{10}{k}_{m}^{2}{r}_{1}+\frac{{q}_{9}(12{q}_{11})}{{q}_{10}}$, ${c}_{2}=(1{q}_{11})(2{q}_{11}1)(\frac{2}{{q}_{10}}\frac{1}{{q}_{10}4{k}_{m}^{2}})+{r}_{2}(2+\frac{1}{1+4D{k}_{m}^{2}})$ and ${c}_{3}=2{q}_{11}{u}_{2}^{(2)}{u}_{3}^{(3)}{u}_{2}^{(2)}{v}_{3}^{(3)}$. Equation 26 is of StuartLandau type and represents a normal form of a subcritical pitchfork bifurcation. This shows the existence of $S{N}_{PB}$ in the EGFR network.
To corroborate this, we also performed numerical bifurcation analysis on onedimensional projection (Equation 20) where the reaction terms have the form as defined in Equation 14, including the full form for $[N{2}_{a}]$, when $[E{E}_{p}]=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: ${\alpha}_{1}=0.001$, ${\alpha}_{2}=0.3$, ${\alpha}_{3}=0.7$, ${\beta}_{1}=11$, ${\beta}_{2}=1.1$, ${k}_{1}=0.5$, ${k}_{2}=0.5$, ${g}_{1}=1.9$, ${g}_{2}=0.1$, ${k}_{on}=0.05$, ${k}_{off}=0.28$, $\u03f5=0.01$, $R{G}_{t}=1$, $N{2}_{t}=1$; and the diffusionlike terms have been scaled from the values derived in Orr et al., 2005: $\stackrel{~}{{D}_{{E}_{p}}}=0.02$, $\stackrel{~}{{D}_{R{G}_{a}}}=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 nonpolarized state losses stability via a symmetrybreaking pitchfork bifurcation ($PB$), which gives rise to a polarized state represented via an inhomogeneous steady states (IHSS). The polarized state is stabilised via saddlenode bifurcations ($S{N}_{PB}$) (Figure 1—figure supplement 1C, magenta branched lines). There is a coexistence between the HSS and the IHSS before the $PB$, rendering it subcritical. 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 frontbackpolarized 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 frontback 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 $S{N}_{PB}$. 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 $S{N}_{PB}$ 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 quasisteady state. The cell boundary was represented with a 1D circular domain of length $L=2\pi R$ (where $R=2\mu 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 $\sigma \cdot d{W}_{t}$ ($\sigma =0.02$, $d{W}_{t}$ is sampled from a normal distribution with mean 0 and variance 0.01) in the equation for $[{E}_{p}]$. The stochastic sdeint Python package was used. Parameters: ${D}_{{E}_{p}}={D}_{R{G}_{a}}=0.008\phantom{\rule{thinmathspace}{0ex}}\mu {m}^{2}/min$. D_{Ep} 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 preactivated states), organization at criticality or in the stable polarized state (IHSS), ${E}_{t}\in \{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 $EG{F}^{647}$ 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 statespace trajectory (Figure 1F, bottom), stochastic realization of the onedimensional projection of the full system (as for the bifurcation analysis) was used.
Physical model of singlecell chemotaxis
Request a detailed protocolTo describe signalinduced 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 signalinduced 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 twodimensional Cartesian grid by the closedcontour $\mathrm{\Gamma}(t)=\{\mathbf{x}\mathrm{\Psi}(\mathbf{x},t)=0\}$, that represent the zerolevel set of the potential function $\mathrm{\Psi}(\mathbf{x},t)$, taken to have an initial form:
where $S$ identifies the area occupied by the cell and $d(\mathbf{x}\mathbf{,}\mathbf{\Gamma}\mathbf{)}$ is the distance of position $x$ to the curve $\mathrm{\Gamma}$. 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 ($\mathrm{\Gamma}(\mathbf{x},t)$) evolves according to the HamiltonJacobi equation:
The vector $\mathbf{v}(\mathbf{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 myosinII 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 k_{c} and a viscous element ${\tau}_{c}$, whereas the cytoplasm is modeled as a purely viscous element, ${\tau}_{a}$, which is placed in series with the Voigt model.
Let $\mathbf{l}(\mathbf{x},t)$, $\mathbf{x}\in \mathrm{\Gamma}(t)$ be the viscoelastic state of the cell at time $t$ and at a position $\mathbf{x}$ on the membrane, such that $\mathbf{l}$ represents the length of the numerous parallel unconnected springdamper systems. The viscoelastic state of the cell then evolves according to:
where ∇ is the gradient operator, the pressure ${\mathbf{P}}_{\mathbf{total}}(t)={\mathbf{P}}_{\mathbf{prot}}(t)+{\mathbf{P}}_{\mathbf{retr}}(t)+{\mathbf{P}}_{\mathbf{area}}(t){\mathbf{P}}_{\mathbf{ten}}(t)$ is sum of the protrusion, retraction, area conservation, and cortical tension pressures, respectively. The EGFR signaling state ($[{E}_{p}]$) directly determines the protrusion/retraction pressure, since high/low signaling activity triggers actin polymerization / myosinII retraction following: ${\mathbf{P}}_{\mathbf{p}\mathbf{r}\mathbf{o}\mathbf{t}}(t)=$ ${K}_{prot}(([{E}_{p}](t)<[{E}_{p}](t)>)$$/([{E}_{p}{]}_{max}(t)<[{E}_{P}](t)>))\mathbf{n}$ and ${\mathbf{P}}_{\mathbf{r}\mathbf{e}\mathbf{t}}(t)=$ ${K}_{retr}((<[{E}_{p}]>[{E}_{p}])$$/(<[{E}_{P}]>[{E}_{p}{]}_{max}))\mathbf{n}$, where <.> denotes mean at the membrane, ${K}_{prot}$, ${K}_{retr}$  proportionality constants. The cell is assumed to be flat with uniform thickness, such that the 2D area ($A(t)$) of the cell is conserved (${\mathbf{P}}_{\mathbf{area}}(\mathbf{t})={K}_{area}(A(0)A(t))\mathbf{n}$), ${K}_{area}$  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 ${\mathbf{P}}_{\mathbf{ten}}(t)={K}_{ten}(\kappa (\mathrm{\Gamma})1/R)\mathbf{n}$, with $\kappa (x)$ being the local membrane curvature, R  initial cell radius, was set to 2 $\mu m$, and ${K}_{ten}$  proportionality constant. The local membrane velocity $\mathbf{v}(\mathbf{x})$, $\mathbf{x}\in \mathrm{\Gamma}(t)$ depends both on the viscoelastic nature of the cell and on the effective pressure profile (${\mathbf{P}}_{\mathbf{total}}(t)$) and is given by,
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 ($[{E}_{p}]$) activity are generated. The viscoelastic state is initialized with zero value on the membrane, $\mathbf{l}(\mathbf{x},0)=0$. At each time point, ${\mathbf{P}}_{\mathbf{total}}$ 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 firstorder 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: $k}_{c}=0.1\phantom{\rule{thinmathspace}{0ex}}nN/\mu {m}^{3$, $\tau}_{c}=0.08\phantom{\rule{thinmathspace}{0ex}}nN\phantom{\rule{thinmathspace}{0ex}}min/\mu {m}^{3$, $\tau}_{a}=0.1\phantom{\rule{thinmathspace}{0ex}}nN\phantom{\rule{thinmathspace}{0ex}}min/\mu {m}^{3$, $K}_{prot}=0.08\phantom{\rule{thinmathspace}{0ex}}nN/\mu {m}^{2$, $K}_{retr}=0.05\phantom{\rule{thinmathspace}{0ex}}nN/\mu {m}^{2$, $K}_{area}=0.02\phantom{\rule{thinmathspace}{0ex}}nN/\mu {m}^{4$, ${K}_{ten}=0.1\phantom{\rule{thinmathspace}{0ex}}nN/\mu m$. $K}_{ten$ was taken from the literature, corresponding to an experimentally measured range of cell cortical tension values (CartagenaRivera 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\pm 0.173\mu 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/Cellmemory, (copy archived at swh:1:rev:288921244e5042922e1bbddcf5037a5e87e78723).

GitHubID Cellmemory. Cells use molecular working memory to navigate in changing chemoattractant fields.
References

Human memory: A proposed system and its control processesPsychology of Learning and Motivation 2:89–195.https://doi.org/10.1016/S00797421(08)604223

A bistable mechanism for directional sensingNew Journal of Physics 10:083015.https://doi.org/10.1088/13672630/10/8/083015

Weakly nonlinear analysis of Turing patterns in a morphochemical model for metal growthComputers & Mathematics with Applications 70:1948–1969.https://doi.org/10.1016/j.camwa.2015.08.019

BookCell Movements: From Molecules to MotilityGarland Science, Taylor and Francis Group.https://doi.org/10.4324/9780203833582

Integrating conflicting chemotactic signals. The role of memory in leukocyte navigationThe Journal of Cell Biology 147:577–588.https://doi.org/10.1083/jcb.147.3.577

Rethinking pattern formation in reaction–diffusion systemsNature Physics 14:507–514.https://doi.org/10.1038/s4156701700405

Transition from amplitude to oscillation death via Turing bifurcationPhysical Review Letters 111:024103.https://doi.org/10.1103/PhysRevLett.111.024103

BookReceptors: Model for Binding, Trafficking and SignalingOxford University Press.

Wavepinning and cell polarity from a bistable reactiondiffusion systemBiophysical Journal 94:3684–3697.https://doi.org/10.1529/biophysj.107.120824

Imaging phosphorylation dynamics of the epidermal growth factor receptorThe Journal of Biological Chemistry 279:36972–36981.https://doi.org/10.1074/jbc.M405830200

Fronts propagating with curvaturedependent speed: Algorithms based on HamiltonJacobi formulationsJournal of Computational Physics 79:12–49.https://doi.org/10.1016/00219991(88)900022

Pattern formation in a coupled membranebulk reactiondiffusion model for intracellular polarization and oscillationsJournal of Theoretical Biology 497:110242.https://doi.org/10.1016/j.jtbi.2020.110242

Fiji: an opensource platform for biologicalimage analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019

Cell motility as persistent random motion: theories from experimentsBiophysical Journal 89:912–931.https://doi.org/10.1529/biophysj.105.061150

Cellular memory in eukaryotic chemotaxisPNAS 111:14448–14453.https://doi.org/10.1073/pnas.1412197111

Organization at criticality enables processing of timevarying signals by receptor networksMolecular Systems Biology 16:e8870.https://doi.org/10.15252/msb.20198870

BookDetecting strange attractors in turbulenceIn: Rand D, Young LS, editors. Dynamical Systems and Turbulence, Warwick 1980. Berlin, Heidelberg: Springer Nature. pp. 366–381.https://doi.org/10.1007/BFb0091924

Parameterspace topology of models for cell polarityNew Journal of Physics 16:065009.https://doi.org/10.1088/13672630/16/6/065009

The chemical basis of morphogenesisPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237:37–72.https://doi.org/10.1098/rstb.1952.0012

On the Theory of the Brownian MotionPhysical Review 36:823–841.https://doi.org/10.1103/PhysRev.36.823

Modeling cellular deformations using the level set formalismBMC Systems Biology 2:1–16.https://doi.org/10.1186/17520509268
Article and author information
Author details
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 PaquinLefebvre for valuable suggestions on the realization of the reactiondiffusion 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/Cellmemory. Funding: The project was funded by the Max Planck Society, partially through the Lise Meitner Excellence Program.
Version history
 Preprint posted: November 12, 2021 (view preprint)
 Received: January 6, 2022
 Accepted: June 3, 2022
 Accepted Manuscript published: June 6, 2022 (version 1)
 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

 2,052
 views

 487
 downloads

 10
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading

 Cell Biology
 Physics of Living Systems
The organelles of eukaryotic cells maintain distinct protein and lipid compositions required for their specific functions. The mechanisms by which many of these components are sorted to their specific locations remain unknown. While some motifs mediating subcellular protein localization have been identified, many membrane proteins and most membrane lipids lack known sorting determinants. A putative mechanism for sorting of membrane components is based on membrane domains known as lipid rafts, which are laterally segregated nanoscopic assemblies of specific lipids and proteins. To assess the role of such domains in the secretory pathway, we applied a robust tool for synchronized secretory protein traffic (RUSH, Retention Using Selective Hooks) to protein constructs with defined affinity for raft phases. These constructs consist solely of singlepass transmembrane domains (TMDs) and, lacking other sorting determinants, constitute probes for membrane domainmediated trafficking. We find that while raft affinity can be sufficient for steadystate PM localization, it is not sufficient for rapid exit from the endoplasmic reticulum (ER), which is instead mediated by a short cytosolic peptide motif. In contrast, we find that Golgi exit kinetics are highly dependent on raft affinity, with raft preferring probes exiting the Golgi ~2.5fold faster than probes with minimal raft affinity. We rationalize these observations with a kinetic model of secretory trafficking, wherein Golgi export can be facilitated by protein association with raft domains. These observations support a role for raftlike membrane domains in the secretory pathway and establish an experimental paradigm for dissecting its underlying machinery.

 Microbiology and Infectious Disease
 Physics of Living Systems
Bacteria in biofilms secrete potassium ions to attract free swimming cells. However, the basis of chemotaxis to potassium remains poorly understood. Here, using a microfluidic device, we found that Escherichia coli can rapidly accumulate in regions of high potassium concentration on the order of millimoles. Using a bead assay, we measured the dynamic response of individual flagellar motors to stepwise changes in potassium concentration, finding that the response resulted from the chemotaxis signaling pathway. To characterize the chemotactic response to potassium, we measured the dose–response curve and adaptation kinetics via an Förster resonance energy transfer (FRET) assay, finding that the chemotaxis pathway exhibited a sensitive response and fast adaptation to potassium. We further found that the two major chemoreceptors Tar and Tsr respond differently to potassium. Tar receptors exhibit a biphasic response, whereas Tsr receptors respond to potassium as an attractant. These different responses were consistent with the responses of the two receptors to intracellular pH changes. The sensitive response and fast adaptation allow bacteria to sense and localize small changes in potassium concentration. The differential responses of Tar and Tsr receptors to potassium suggest that cells at different growth stages respond differently to potassium and may have different requirements for potassium.