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

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

Figure 1 with 2 supplements see all

Download asset Open asset

We described this conjecture mathematically for a general reaction-diffusion 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 zero-crossing eigenvalues ($\lambda$) associated with the odd mode of the perturbation (Paquin-Lefebvre et al., 2020). To probe the sub-critical transition and therefore the necessary organization at criticality, a reduced description in terms of an asymptotic expansion of the amplitude of the polarized state ($\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 above-mentioned principle, we use the well-characterized Epidermal growth factor receptor (EGFR) sensing network (Reynolds et al., 2003; Baumdick et al., 2015; Stanoev et al., 2018). It constitutes of double negative and negative feedback interactions of the receptor, EGFR ($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 each-other’s activities (see Materials and methods for the molecular details of the network). The bidirectional molecular interactions between EGFR and the phosphatases can be mathematically represented using mass action kinetics, giving a system of partial differential equations (PDE) that describes how the dynamics of the constituents evolves in time and space (Equation 14 in Materials and methods). Applying a weakly nonlinear stability analysis (Becherer et al., 2009) to this system of equations shows that the EGFR phosphorylation dynamics undergoes a symmetry-breaking transition ($PB$) as outlined above (proof in Materials and methods, Figure 1—figure supplement 1C). The $PB$ generates a polarized state that is represented as a inhomogeneous steady state (IHSS) - a combination of a high receptor phosphorylation at the cell front and low in the back of the cell (schematically shown in Figure 1A, Figure 1—figure supplement 1A). This is contrary to a bistable system, where the polarized signaling state would be manifested by two steady states, high and low protein phosphorylation in the front and back of the cell, respectively (Beta et al., 2008). This profiles $PB$ as a robust mechanism of cell polarization. Polarized EGFR signaling on the other hand, will lead to reorganization of the cortical actomyosin cytoskeleton by regulating members of the Rho GTPase family, thereby inducing signal-dependent cell shape changes and subsequent migration (Chiasson-MacKenzie and McClatchey, 2018; Ridley and Hall, 1992). In order to link signaling activity with morphodynamics, we modeled the cell as a viscoelastic cortex surrounding a viscous core (Yang et al., 2008) (Materials and methods), where EGFR signaling dynamics affects cell shape changes through the protrusion/retraction stress and the viscoelastic nature of the cell membrane (Figure 1C).

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

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

To test experimentally whether cells maintain memory of the direction of previously encountered signals through prolonged EGFR phosphorylation polarization, and what is the duration of this effect, epithelial breast cancer-derived MCF7 cells were subjected for 1 hr to a stable gradient of fluorescently tagged EGF-Alexa647 (EGF⁶⁴⁷) with a maximal amplitude of 10 ng/ml applied from the top of the chamber in a computer-programmable microfluidic device (Figure 2A and B). EGFR phosphorylation at the plasma membrane was quantified during and for 3 hr after gradient wash-out (gradient wash-out established in 4–5 min) by determining the rapid translocation of mCherry-tagged phosphotyrosine-binding domain (PTB^mCherry) to phosphorylated tyrosines 1086/1148 of ectopically expressed EGFR-mCitrine (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.

Figure 2 with 5 supplements see all

Download asset Open asset

Figure 2—source data 1

Source data for Figure 2.

https://cdn.elifesciences.org/articles/76825/elife-76825-fig2-data1-v2.xlsx

Download elife-76825-fig2-data1-v2.xlsx

Kymograph analysis of ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation at the plasma membrane of single cells showed polarization in a shallow gradient of EGF⁶⁴⁷ (as shallow as 10% between front and back of the cell; Figure 2C, Figure 2—figure supplement 1A-D). The direction of ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation polarization coincided with the direction of maximal EGF⁶⁴⁷ concentration around each cell ($\pi /4$ on average, Figure 2—figure supplement 1F). Only few cells manifested basal or symmetric ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation distribution upon gradient stimulation (Figure 2—figure supplement 1A, B, E). Plotting the fraction of plasma membrane area with polarized ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation showed cell-to-cell variability in the polarization kinetics, as well as the maximal amplitude of polarized ${\displaystyle {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 ${\displaystyle {EGFR}^{mCitirne}}$ phosphorylation was also reflected in the respective single-cell temporal profiles (exemplary profile shown in Figure 2F, top). Reconstructing the state-space trajectory from this temporal profile using Takens’s delay embedding theorem (Takens, 1980) (Materials and methods) showed that before the fast transition to the basal state, the trajectory of the system was trapped in the vicinity of the polarized state (Figure 2F bottom, Figure 2—video 1). Despite the biological and technical noise that affect the measurement of the temporal ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation profile, and thereby the reconstruction of the state-space trajectory, both qualitatively resemble the equivalent numerical profiles (compare Figures 1F and 2F). In contrast, when cells were subjected to an ATP analog EGFR inhibitor (Björkelund et al., 2013) during gradient wash-out, the ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation response exponentially decayed, resulting in a clear absence of transient memory and respective state-space trapping (Figure 2G, Figure 2—figure supplement 2A, B, Figure 2—video 2). Since Lapatinib inhibits the kinase activity of the receptor, the dynamics of the system in this case is mainly guided by the dephosphorylating activity of the phosphatases. Implementing an equivalent of the Lapatinib inhibition in the numerical simulations by decreasing the autocatalytic EGFR activation rate constant after gradient removal verifies that the presence of memory in EGFR phosphorylation cannot be explained only by a dephosphorylation process (Figure 2—figure supplement 2C). This is also evident from the respective state-space trajectory, where the system directly transits from the polarized to the basal state, without intermediate state-space trapping (Figure 2—figure supplement 2D, Figure 2—video 3).

Fitting the experimentally measured single-cell temporal ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation profiles after gradient wash-out using an inverse sigmoid function (Methods) further corroborated that under Lapatinib treatment, phosphorylated ${\displaystyle {EGFR}^{mCitrine}}$ exponentially relaxed from the polarized to the basal state (Hill coefficient ≈ 1.28), with a half-life of approx. 10 min (Figure 2—figure supplement 2E, G). Under normal conditions however, the half-life was 30 min on average, reflecting that the phosphorylated ${\displaystyle {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 ${\displaystyle {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 ${\displaystyle {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 single-cell solidity profile over time showed that epithelial cells maintained the polarized cell shape for ~ 40 min after signal removal (Figure 2H, Materials and methods), which directly corresponds to the average memory duration in polarized ${\displaystyle {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 ${\displaystyle {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.

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

Figure 3 with 5 supplements see all

Download asset Open asset

Figure 3—source data 1

Source data for Figure 3.

https://cdn.elifesciences.org/articles/76825/elife-76825-fig3-data1-v2.xlsx

Download elife-76825-fig3-data1-v2.xlsx

This was also reflected in the projection of the cell’s relative displacement angles ($\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 ($\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 ($\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 time-frame in which the memory in polarized ${\displaystyle {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 Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein, 1930; Svensson et al., 2017) and used the extracted migration parameters to generate synthetic single-cell trajectories (Materials and methods). In absence of stimulus, the cellular motion resembled a random walk process (RW: Figure 3D right, Figure 3—figure supplement 2D, E middle), persistent random walk (PRW) was characteristic for the uniform stimulation case (Figure 3—figure supplement 2D, E right), whereas biased PRW described the migration in gradient presence (PBRW, Figure 3D- and Figure 3—figure supplement 2D, left, green trajectory part). Extending the bias duration during the interval of the experimentally observed memory phase (PB(t)RW) was necessary to reproduce the transient persistent motion after gradient removal (Figure 3D- and Figure 3—figure supplement 2D, left, blue trajectory part; Figure 3E and F; Figure 3—figure supplement 2F).

To corroborate the link between memory in polarized receptor activity, memory in polarized cell shape and memory in directional migration, we also quantified the directional migration of MCF10A cells when subjected to Lapatinib during gradient wash-out (Figure 3G). The directionality after gradient removal was significantly lower than in the case without Lapatinib (Figure 3H), suggesting that cells rapidly switch to a RW migration pattern upon gradient wash-out due to the absence of memory in polarized ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation (Figure 2G,I). Thus, single-cell motility trajectories that closely resembled the experimentally observed ones could be mimicked with the PB(t)RW simulation, where the bias duration corresponded to the duration of the gradient (Figure 3—figure supplement 2E left, 2G). Quantification of the average cells’ relative displacement angles showed as well that $\cos \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 cell-to-cell variability in this case, we also calculated memory duration from single cell $\cos \theta$ profiles. For this, single-cell trajectories were first smoothed using Kalman filter (Materials and methods). The quantification showed that majority of the cells displayed absence of or shorter memory in directional migration, with a mean value of ∼25 min (Figure 3—figure supplement 3A, B, D). Since under Lapatinib treatment, EGFR phosphorylation rapidly decays (Figure 2G), this residual memory in some cells likely results from memory in cytoskeletal asymmetries, as previously suggested (Prentice-Mott et al., 2016). Without Lapatinib treatment however, the duration of memory estimated from single-cell $\cos \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 single-cell kymograph quantification (Figure 2E).

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

Figure 4 with 8 supplements see all

Download asset Open asset

Figure 4—source data 1

Source data for Figure 4.

https://cdn.elifesciences.org/articles/76825/elife-76825-fig4-data1-v2.xlsx

Download elife-76825-fig4-data1-v2.xlsx

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

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

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

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

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

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

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

Wide field images were acquired using an Olympus IX81 inverted microscope (Olympus Life Science) equipped with a MT20 illumination system and a temperature controlled C0₂ incubation chamber at 37 °C and 5% CO₂. 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.

The CellAsic Onix Microfluidic Platform (EMD Millipore) was used for gradient cell migration and ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation polarization experiments. For ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation polarization experiments, 1 hr gradient stimulation was established using CellASIC ONIX2 software as follows. (i) Pre-stimulus: Imaging media was flowed from well groups 3 and 4 (CellAsic Onix Manual - https://www.merckmillipore.com/) at low pressure (2.5 kPa) for 5 min. (ii) Gradient establishment: After closing well group 3, pre-loaded EGF⁶⁴⁷ (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⁶⁴⁷ 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⁶⁴⁷/ 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⁶⁴⁷ (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⁶⁴⁷ (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\min$.

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

For migration experiments under uniform EGF⁶⁴⁷ stimulation, confocal laser scanning microscopy / transmission imaging of live ${\displaystyle {MCF7-EGFR}^{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.

hEGF⁶⁴⁷ was generated in the lab of Prof. P. Bastiaens, MPI of molecular Physiology, Dortmund, using the His-CBD-Intein-(Cys)-hEGF-(Cys) plasmid (Sonntag et al., 2014), kindly provided by Prof. Luc Brunsveld, University of Technology, Eindhoven. Human EGF was purified from E. coli BL21 (DE3), N-terminally labeled with Alexa647-maleimide as described previously (Sonntag et al., 2014) and stored in PBS at –20 °C. To quantify the spatial extent of the EGF⁶⁴⁷ /Fluorescein gradient, gradients were generated following the protocol described in sub-section 5.6 in plates without cells or matrix coating. Confocal images of Alexa647 /GFP channel were acquired at $1\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⁶⁴⁷ /Fluorescein concentration.

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

where $PT{B}_{PM}^i(t)$ and $EGF{R}_{PM}^i(t)$ are respectively the ${\displaystyle {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 ${\displaystyle {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 ${\displaystyle {PTB}^{mCherry}}$ fluorescence from these structures was subtracted from the $PT{B}_T(t)$, to correct for the ${\displaystyle {PTB}^{mCherry}}$ fraction bound to the phosphorylated ${\displaystyle {EGFR}^{mCitrine}}$ on endosomes.

Temporal profile of the fraction of phosphorylated ${\displaystyle {EGFR}^{mCitrine}}$ on the plasma membrane was obtained using:

and then normalized as:

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

To classify single cells into non-activated, activated (polarized ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation) and pre-activated (uniformly distributed ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation) upon gradient EGF⁶⁴⁷ stimulation (Figure 2—figure supplement 2A, B), the following method was applied. To identify pre-activated cells, a Gaussian Mixture Model (GMM) was fitted to the histogram of ${(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 pre-activated. To distinguish between the non-activated and activated cells in the remaining population, average ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation value ($EGF{R}_p$) per cell was estimated during the pre-stimulation ($t\in [0,5min]$) and the stimulation period ($t\in [5min,65min]$) (${\displaystyle <EGF{R}_p{>}_{t\in [0,65]}}$) from the temporal ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation profiles. Histogram of the respective $EGF{R}_p$ values was again fitted with a GMM model. All cells with an average ${\displaystyle <EGF{R}_p{>}_{t\in [0,65]}}$ value lying below the intersection point were considered to be non-activated, whereas those above - activated.

The average of the spatial projection of the fraction of phosphorylated ${\displaystyle {EGFR}^{mCitrine}}$ from single-cell kymographs (Figure 2—figure supplement 1C) was generated from the 20 cells that were polarized in the direction of the EGF⁶⁴⁷ 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⁶⁴⁷ distribution around single cells (Figure 2—figure supplement 1F) was estimated as follows: the cell mask obtained using the ${\displaystyle {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⁶⁴⁷ outside single cells. This outer contour was divided in 20 bins as for the kymographs, and EGF⁶⁴⁷ intensity was quantified in each bin. The angle between the direction of EGF⁶⁴⁷ and the direction of EGFR phosphorylation was calculated as the amount of radial bins between the maxima in the spatial projections. This bin-distance was then translated into an angle under the assumption of a circular perimeter.

In order to identify the characteristic features of the ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation profile during the transition from polarized to unpolarized state, the single-cell $EGF{R}_p(t)$ profiles with and without Lapatinib treatment after gradient wash-out were fitted to an inverse sigmoid function given by,

were a₀, $a$ are constants and $n$ is the Hill-coefficient (examples in Figure 2—figure supplement 2E, F). Non-linear least square method (python package curve fit) was used to perform the fitting. Under normal conditions (w/o Lapatinib), $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₀ 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, half-life was estimated to be the time frame in which 50% of ${\displaystyle {EGFR}^{mCitrine}}$ phosphorylation is lost after EGF⁶⁴⁷ removal.

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

yielding the value of phosphorylated ${\displaystyle {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 ${\displaystyle {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 ${\displaystyle FP{A}_{per-cell}<(FP{A}_{average}-s.d.)}$ in 3 consecutive time points (Figure 2E).

Morphological changes of polarized cells were quantified using the solidity (Figure 2H1) of each cell at each time point and the directed protrusive area towards and away from the gradient (Figure 1G and H; Figure 2—figure supplement 2H). The solidity $\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 single-cell solidity profiles, and corresponds to the time-point at which the solidity is below mean-s.d. estimated during gradient presence. The directed cell protrusion area was estimated by comparing single cell masks at two consecutive time points. To reduce noise effects, the masks were first subjected to a 2D Gaussian filtering using the $filters.gaussian$ function from the $scipy$ python package. Protrusions were considered if the area change was greater than 10 pixels or 1.2 µm² 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.

Single cell migration trajectories were extracted using Trackmate (Tinevez et al., 2017) in Fiji (Schindelin et al., 2012) using Hoechst 33342/transmission channel. From the positional information (x and y coordinates) of individual cell tracks, quantities such as Motility, Directionality and $\cos \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 two-sided Welch’s t-test. To quantify the memory duration in directed single-cell migration, the Kernel Density Estimate (KDE) from $\cos \theta$ quantification in the continuous absence of EGF⁶⁴⁷ (uniform case, between 250–300 min) was compared with a moving window KDE (size of 5 time points) from the gradient migration profile, using two sided Kolmogorov-Smirnov test. To verify the absence of memory when cells were treated with Lapatinib during gradient wash-out, a moving window KDE (5 time points) from $\cos \theta$ obtained in this case was compared to the KDE in continuous absence of EGF⁶⁴⁷ (uniform case Figure 3—figure supplement 2B, between 250–300 min) using two sided Kolmogorov-Smirnov test (Figure 3I). Furthermore, the KDE between 300–350 min and 350–840 min (after gradient removal) was statistically equivalent to the KDE in continuous absence of EGF⁶⁴⁷, confirming the rapid switch from directed to random-walk migration in the Lapatinib case (Figure 3—figure supplement 2H). To estimate the time required for complete reversal of cell migration direction when the cells were subjected to a gradient from opposite direction, KDE distributions were compared between the following time windows: 275–335 min (second gradient), 335–365 min, 365–385 min, 375–385 min, and 365–605 min (third gradient).

To quantify the motility patterns of MCF10A cells in absence, uniform or gradient $EG{F}^{647}$ stimulation, we fitted the experimentally obtained single cell migration trajectories using modified Ornstein-Uhlenbeck process (mOU) (Uhlenbeck and Ornstein, 1930) that is defined by the Langevin equation for the velocity vector $\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 time-dependent 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 ${\displaystyle MSD(t)=<|{\mathbf{x}}_i(t)-{\mathbf{x}}_i(0){|}^2>}$, where ${\mathbf{x}}_i(t)$ is the tracked position of i-th cell in the 2D plane, <> is the average across all single cell tracks, and $|.|$ is the Euclidean distance (Selmeczi et al., 2005). To estimate D, the obtained MSD profile was fitted with a linear function ($=4Dt$). Goodness of Fit for the different experimental conditions: 0 ng/ml EGF⁶⁴⁷, $R^2=0.975$; for uniform 20 ng/ml EGF⁶⁴⁷ stimulation, $R^2=0.995$. In order to estimate $\tau$, Velocity Auto-Correlation Function ${\displaystyle 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⁶⁴⁷ case - Standard Error of Estimate $SEOE=0.0261$; for uniform 20 ng/ml EGF⁶⁴⁷ stimulation case, $SEOE=0.0570$. Fitted values: for 0 ng/ml EGF⁶⁴⁷ case, $\tau =11.105$, $D=0.425$; for uniform 20 ng/ml EGF⁶⁴⁷ 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 Kalman-filter (python package filterpy.kalman) by predicting the cell position and velocity. The cell’s displacement angles relative to the gradient direction ($\cos \theta$) were calculated for each cell at each timepoint, rendering single-cell $\cos \theta$ plots (Figure 3—figure supplement 3B,C). The memory duration was then calculated as the point where three consecutive timepoints in the $\cos \theta$ profiles fall below a threshold $\cos \theta$ value of 0.75.

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

where $i=1$ to ${\displaystyle N-(m-1)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 one-to-one image of the original set, provided $m$ is large enough. The embedding delay was identified using the timeLag function (based on autocorrelation), the embedding dimension using the estimateEmbeddingDims function (based on the nearest-neighbours method), and the state-space reconstruction using the buildTakens function, all from the nonlinearTseries package in R (https://cran.r-project.org/web/packages/nonlinearTseries/index.html). Before state-space reconstructions, time series were smoothened using the Savitzky-Golay filter function in Python. For Figure 2F, $d=26$, $d_e=3$; for Figure 2G, $d=50$, $d_e=3$.

We consider a generalized form of a (mass-conserved) reaction-diffusion (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 Turing-type instability (Turing, 1952), such that the symmetric steady state becomes unstable and an asymmetric polarized state is stabilized. By its nature, the linear stability analysis makes no prediction about the transition process itself, and thereby the type of bifurcation that underlies it. To provide quantitative description of the symmetry breaking transition in reaction-diffusion models, local perturbation analysis can be applied (Holmes et al., 2015). However, this analysis is mainly restricted to models characterized with large diffusion discrepancy between the signaling components. The conditions for a pitchfork bifurcation ($PB$)-induced transition in a generic RD model therefore have to be formally defined. Let ${\mathbf{U}}_{\mathbf{s}}=(u_{is})$ for ${\displaystyle 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, ${\displaystyle 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 (Paquin-Lefebvre 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 front-back-polarized states, where the orientation depends on the direction of the external signal (Figure 1A, Figure 1—figure supplement 1A).

To identify whether the PB is of sub-critical type, and thereby identify the presence of an $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 one-dimensional system ($N=1$). As we are interested in the description of a structure of finite spatial size (i.e. finite wavelength $k$), the final solution of the PDE is expanded around the fastest growing mode, 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 Stuart-Landau type:

with ${\displaystyle c_1,c_2,c_3>0}$, this corresponds to the normal form of a sub-critical pitchfork bifurcation (Strogatz, 2018). Together with the condition given by Equation 11, the existence of a sub-critical PB for the full system (Equation 9) is guaranteed. A numerical or analytical analysis of Equation 13 enables the identification of the position of the $S{N}_{PB}$.

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

with

and

The reaction terms are described in details in Stanoev et al., 2018. In brief, $[E-E_p]$ is the phosphorylated ligand-bound 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 ER-bound PTPN2, respectively. Both, the receptor and the deactivating enzymes have active and inactive states, and the model equations describe their state transition rates. Therefore, mass is conserved in the system and the total protein concentrations of the three species ($[E_t]$, $[R{G}_t]$ and $[N{2}_t]$) are constant parameters. Autonomous, autocatalytic and ligand-bound-induced activation of ligandless EGFR ensue from bimolecular interactions with distinct rate constants ${\alpha }_{1-3}$, respectively. Other parameters are as follows: $k_1/k_2$ — activation/inactivation rate constants of the phosphatases, ${\beta }_1/{\beta }_2$ - receptor-induced regulation rate constants of $PTPRG/PTPN2$, ${\gamma }_1/{\gamma }_2$ - specific reactivity of the enzymes ($PTPRG/PTPN2$) towards the receptor. The EGFR-PTPN2 negative feedback is on a time scale ($ϵ$) approximately two orders of magnitude slower than the phosphorylation-dephosphorylation reaction, as estimated from the ~ 4 min recycling time of $EGF{R}_p$ (Stanoev et al., 2018). This enables, when necessary, to consider a quasi-steady 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), ligand-binding 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 ER-bound and does not diffuse. Single particle tracking studies have demonstrated that EGFR molecules on the plasma membrane occupy three distinct mobility states, free, confined and immobile, with the occupations of the free and immobile states decreasing and increasing significantly after EGF stimulation, respectively (2 min after EGF stimulation, corresponding with the time-scale of EGF binding) (Ibach et al., 2015). In the reaction-diffusion (RD) simulations therefore for simplicity, it is assumed that $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}$).

To 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 reaction-diffusion signaling model). For this, we considered the system Equation 14, where the dynamics of PTPN2 is at quasi-steady state (Equation 15), $[E-E_p]=0$, and rest of the dependent and independent variables were scaled to have a dimensionless form. Let $[\widetilde{E_p}]=[E_p]/E_0$, $[\widetilde{R{G}_a}]=[R{G}_a]/R{G}_0$, $\tilde{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 ${\displaystyle q_1=\frac{{\alpha }_1\cdot [E_t{]}^2}{(k_1+k_2)\cdot {\beta }_2}}$, ${\displaystyle q_2=\frac{({\alpha }_2-2\cdot {\alpha }_1)\cdot [E_t]}{k_1+k_2}}$, ${\displaystyle q_3=\frac{({\alpha }_1-{\alpha }_2)\cdot k_t}{(k_1+k_2)\cdot {\beta }_2}}$, ${\displaystyle q_4=\frac{{\gamma }_2\cdot [N{2}_t]}{k_1+k_2}}$, ${\displaystyle k=k_2/k_1}$, ${\displaystyle r_1=\frac{k_1\cdot [R{G}_t]\cdot {\gamma }_1}{(k_1+k_2{)}^2}}$, ${\displaystyle r_2=\frac{{\beta }_1\cdot k_1}{(k_1+k_2)\cdot {\beta }_2}}$ and ${\displaystyle D=\frac{D_{R{G}_{\alpha }}}{D_{E_p}}}$.

We further simplify the system Equation 16 by taking the Talyor series expansion of the quasi-steady state approximation of $[N{2}_a]$ around ${\displaystyle E_s}$, the steady state of ${\displaystyle [\widetilde{E_p]}}$

with ${\displaystyle q_7=\frac{E_s{q}_4}{1+k+E_s}-\frac{E_s{q}_4(1+k)}{(1+k+E_s{)}^2}}$, ${\displaystyle 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 ${\displaystyle q_9=\frac{q_4(1+k)}{(1+k+E_s{)}^2}}$, thus yielding:

with ${\displaystyle q_9=q_1-q_7,q_{10}=q_2-q_8}$ and ${\displaystyle q_{11}=q_3-q_9}$.

To avoid long expression in the further analysis, we re-name the dependent variables as $u_1=[\widetilde{E_p}]$ and $u_2=[\widetilde{R{G}_a}]$, and the independent variables as $\tilde{x}=x$, $\tau =t$. The system Equation 16 therefore obtains the generic form:

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

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

Let, ${\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