Optimal biochemical information processing at criticality

How cells utilize surface receptors for chemoreception is a recurrent question spanning between physics and biology over the past few decades. However, the dynamical mechanism for processing time-varying signals is still unclear. Using dynamical systems formalism to describe criticality in non-equilibrium systems, we propose generic principle for temporal information processing through phase-space trajectories using dynamic transient memory. In contrast to short-term memory, dynamic memory generated via ghost attractor enables signal integration depending on stimulus history, and thus balance between stability and plasticity in receptor responses. We propose that self-organization at criticality can arise through fluctuation-sensing mechanism, illustrated for the experimentally established epidermal growth factor sensing system. This framework applies irrespective of the intrinsic node dynamics or network size, as we show using also a basic neuronal model. Processing of non-stationary signals, a feature previously attributed only to neuronal networks, thus uniquely emerges for biochemical networks organized at criticality.


I. INTRODUCTION
In a wide variety of biological processes including embryogenesis, immune cells motility, wound healing or cancer metastasis [1][2][3], cells sense and respond to timevarying chemical signals that reflect the non-stationary environment to which they readily adapt. Sensing of chemical signals occurs through receptor-coupled enzymatic activity that transmits the information about the environment through chemical modifications. The receptor activity dynamics on the other hand emerges from the biochemical network in which the receptor is embedded. To perceive and process continuous streams of multimodal inputs, the cellular sensing system must satisfy two general, but seemingly opposed demands. Sensing requires plasticity in the receptor activity responses to enable rapid responsiveness to continuous changes in the environment, whereas robust signaling responses require prolonged receptor activity through a transient memory to be maintained on the cell surface after signal removal [4]. Such a memory is necessary to process the temporal dependencies that are inherent in time-varying signals [5]. However, a theoretical framework that explains how these opposed features emerge on the level of single cells has been missing, since common models for computation such as Turing machines or attractor networks cannot capture this dynamics. While Turing machines describe off-line computation on (static) discrete inputs [6], attractor networks that use multiistability for memory realization settle into a fixed-point and can thereby temporally limit the perception of upcoming stimuli [4].
Processing of time-varying chemical signals by cell surface receptors on the other hand dynamically resembles the real-time computations of sensory stimuli carried out by neural microcircuits in the cerebral cortex [7,8]. To * Correspondence to: aneta.koseska@mpi-dortmund.mpg.de overcome the limitations of Turing computation in the latter case, universal frameworks using transient dynamics and state-dependent trajectories have been proposed as principal forms of computations that rely on the highdimensionality of the networks with complex intrinsic neuronal dynamics on several temporal scales [7,9,10]. These formalisms however cannot be directly translated to the equivalent computational problem in biochemical networks. First, receptors do not exhibit intrinsic complex activity dynamics such as neurons, and second, cortical microcircuits are usually consisted of several hundreds of neurons, whereas the core receptor networks can have as few as four nodes [4].
We propose here a theoretical framework for cellular processing of time-varying signals in a vicinity of saddlenode (SN) bifurcation, where the information is encoded through metastable state generated via ghost attractor, and interpreted using phase-space trajectories. In contrast to networks with short-term memory, we demonstrate that the dynamical transient memory emerging due to this critical organization enables signal integration and optimal information processing even with twocomponent sensing systems. Using a single spiking leaky integrate-and-fire (LIF) neuron model [11] we demonstrate that these features are exhibited for a different type of intrinsic dynamics, and are thereby generic. Drawing a parallel to epidemic spreading on networks, we explain how such a critical behavior can take place on a molecular level of cell surface receptors. Using the example of the epidermal growth factor receptor (EGFR) system, we propose that a simple fluctuation-sensing mechanism could enable information processing capabilities to emerge in a self-organized manner, as a generic means to balance between plasticity and robustness in cellular sensing.

II. DYNAMICAL INFORMATION PROCESSING AT CRITICALITY USING STATE-DEPENDENT TRAJECTORIES
Computational outcomes of systems that sense timevarying signals depend not only on the current stimulus, but also on memory of the previous sequence(s) in order to integrate the signal information [5]. The necessity for such a memory in receptor activity on the cell surface is even more pronounced to enable robust signaling response, since the majority of the receptors are rapidly internalized upon ligand binding and unidirectionally trafficked towards the cell interior where they are degraded [12,13]. In broad range of biological systems, attractor networks are generically thought to underlie these memory functions [14][15][16]. A minimal cellular sensing network that accounts for memory in receptor activity (R a ) is a two-component toggle-switch (Fig. 1a), where the double-negative feedback interaction [4,17,18] between the active receptor and an inactivating enzyme, a protein P DNF,a , follows the law of mass action: The system integrates the R a activation-and the mutual inhibition mechanisms (Fig. 1a) that govern the protein state transitions between their active (R a , P DN F,a ) and inactive (R i , P DN F,i ) states. They are described in further detail in Appendix A with the corresponding parameters.
Bistability is exhibited between two saddle-node bifurcations for a broad range of the bifurcation parameterthe P DN F,T /R T concentration ratio -in absence of chemical stimulus input (Fig. 1b), and it is also maintained for a certain input range (Fig. 1c). The effective input for the cells in this case is the fraction of ligand-bound receptors that reflects the extracellular ligand concentration (Appendix A and [4]). Processing time-varying signals however, necessitates that the robustness that arises in attractor networks via bistability must be balanced with plasticity in the receptor responses in order to maintain sensitivity to upcoming stimuli.
To understand how this can occur in bistable systems, we studied qualitatively the dynamical R a -P DNF,a behavior by analyzing how the phase-space trajectories evolve in relation to the changes in the geometry of the underlying phase space as a function of a pulsed stimulus. Generally, the relative positioning of the nullclines, which are determined by the system parameters, shapes the phase space geometry. In non-autonomous or inputdriven systems, the underlying phase space can be altered either through its geometry (change in the positioning, shape and size of the attractors) or topology (change in the number or stability of the attractors) [19,20]. We therefore also estimated the associated quasi-potential landscapes [21] (Figs. 1d to 1f, Appendix B and [19]). When starting from the valley of basal receptor activity in the double-well quasi-potential landscape characteristic for the bistable organization (Fig. 1b, green) in absence of stimulus (Fig. 1e middle), a topological phase space change that reflects the vanishing of this state occurs at a threshold signal concentration. This is manifested through a transition to the monostable state with high receptor activity (Fig. 1e, i →ii green transition). Upon signal removal, the reverse topological change leads to re-establishing of bistability (ii →iii green transition). However the trajectory remains in the occupied high activity steady state (green circle in Fig. 1e middle, and top). Thus, the first pulse will activate the receptors and this will hinder further responsiveness to upcoming stimuli due to the long-term memory that results from this stable attractor organization ( Fig. 1e top, inset).
In contrast, organization in the monostable regime ( Fig. 1b, blue) does not result in memory in receptor activity ( Fig. 1d top and inset). The continuous and reversible repositioning of the single steady-state attractor that captivates the state trajectory through the induced changes in the phase space geometry (adding/removing stimulus: i →ii, ii →iii, blue transitions, Fig. 1d) in this case leads to receptor response that closely follows that of the input (Fig. 1d top, inset). This indicates that neither bi-nor monostability can simultaneously account for plasticity and robustness in cellular responses to timevarying chemical cues.
For positioning in the vicinity of the saddle-node bifurcation point however (Fig. 1b, magenta), there is only one stable attractor -the basal activity state (Fig. 1f middle). However, due to the closeness to the bifurcation point, the dynamics of the system here is qualitatively different than the remaining the monostable region. For a supra-threshold input pulse, the transition from the basal monostable -to the high activity monostable state (i →ii magenta transition) via the bistable region, determines robust activation of the receptor, whereas upon input removal, these consecutive topological transitions are reversed. During the reversal, there is a delay between establishing the single stable attractor (magenta state iii ), and the system trajectory converging to it (magenta state iv ), resulting in prolonged receptor activity before relaxation to basal level (Fig. 1f top). This delay causes a transient memory in receptor activity that does not hinder further responsiveness of the system (Fig. 1f top, inset).
The transient memory is a consequence of the critical dynamical behavior near the SN bifurcation. In this organization, the nullclines intersect only once, indicating a single low stable steady state. However, they are positioned very close to one another in the phase space area where the high steady state would be stable (compare Magenta arrows: signal administration (i →ii ) and removal (ii →iii →iv ). The iii →iv transition and the associated phase space plot demonstrate the existence of a ghost attractor. (g) Bifurcation diagram of the LIF neuron model (Appendix C), depicting membrane voltage as a function of the input current (I app ). HB: Hopf bifurcation, SNLC: saddle-node on limit cycle. Solid/dashed lines: stable/unstable steady state (black) and limit cycle (red). Shading same as in b. (h) Firing of a single LIF neuron during an input current increase (yellow) followed by transient memory after input resetting. Magenta: envelope of the neuronal response. Inset: Response upon consecutive inputs.
Parameters: Appendix C.
Figs. 1f and 1e, middle), resulting in a quasi-potential landscape with a very shallow slope ( Fig. 1f middle, inset). Thus, when the system transits back from the bistable to the monostable region, the remnant of the saddle-node that disappeared in this transition generates a metastable state that continues to capture the incoming trajectories (See Supplemental Video [22]). Such delayed dynamics referred to as a feature of a ghost attractor [23], has been previously shown in some driven dynamical systems such as ferroelectrics or semiconductor lasers [24]. In this organization, the system does not process information using the stable attractors associated with specific states (0-1 case), but rather the information is maintained via the transient memory and interpreted through the phase-space trajectories. The trajectories are in turn navigated through the dynamic quasi-potential landscape by the input-driven topological transitions. This is a generic feature of systems organized at a critical proximity to a saddle-node bifurcation point. We demonstrate this on a system with different intrinsic dynamics -a single leaky integrate-and-fire (LIF) neuron model (Eq. (C1) in Appendix C). Bifurcation analysis showed that in this case, bistability between a resting and a spiking state is marked by a saddle-node on a limit cycle (SNLC) bifurcation (Fig. 1g). When organized in the vicinity of the SNLC, a continuous supra-threshold current input induces repetitive neuronal firing, which is transiently maintained after input resetting, due to the presence of the ghost of the SNLC (Fig. 1h). The envelope of the neuronal pulsing activity closely resembles the receptor activity profile at criticality (compare Fig. 1h and 1f top). Similar neuronal behavior where switching from a low-activity (baseline) state into a stimulusselective activity state that is maintained throughout a delayed period after signal removal has been experimentally related to working memory processes in the large-scale neuronal networks in primate prefrontal cortex [25,26]. It has also been demonstrated that stable attractor networks, as well as feedforward chain and random chaotic networks, cannot account for both stable coding as well as strong temporal dynamics of the neuronal population during working memory [26]. Thus, the critical behavior emerging for organization in the vicinity of the saddle-node bifurcation, although here simplified for a single neuron, could possibly provide the underling dynamical mechanism that supports working memory activity in the prefrontal cortex.

III. TRANSIENT MOLECULAR MEMORY INTERPRETED THOUGH CRITICAL EPIDEMIC SPREADING
To understand how transient memory can be realized on a molecular level, we studied how transient receptor activity can be generated and maintained in absence of stimulus using single molecule reaction-diffusion simulation framework on a two-dimensional surface (Ap-pendix D). Microscopic, or single-molecule receptor dynamics can be regarded as analogous to the susceptibleinfected-susceptible (SIS) epidemic spreading system [27,28]. In this analogy, single receptor molecules are susceptible to activation, equivalently to the S-state of an agent in the SIS model, and once active they can then propagate their state via direct contact (infected state of the SIS model). Here, receptor molecules can become susceptible again by interacting with active P molecules. The basic reproduction number R 0 that determines the transmission potential of an infection [29] corresponds to the average number of newly 'infected' molecules by a single active receptor molecule in the course of its active lifetime, i.e. before its deactivation. If R 0 < 1, the overall activity in the system decays (Fig. 2a, top), whereas if R 0 > 1, the system exhibits supercritical behavior and the activity propagates in an epidemic-like branching fashion (Fig. 2a, bottom).
To relate the specific R 0 realization when crossing the saddle-node bifurcations, we derived analytically the dependence of R 0 to the main parameter that determines the positioning of the system in the vicinity of the SN bifurcation. For simplicity, we use γ DN F as a bifurcation parameter, which is proportional to the specific reactivity of P DNF,a to R a (γ DN F ∝ P DN F,T /R T , Fig. 1b). The activation transmission potential in every point of the phase space can be described as are the average molecular activation rate and lifetime, respectively (kinetic parameters given in Appendix D). When initiated at the basal state, hence with a fully susceptible population (R a = 0, P DN F,a = k 1 /(k 1 + k 2 ), Appendix D), due to R 0 = α 2 R T (k 1 + k 2 )/(γ DN F P DN F,T k 1 ) ∝ 1/γ DN F the critical threshold R 0 = 1 is crossed at the γ DN F value corresponding to SN1 (Fig. 2b bottom and middle; see Appendix D). For γ DN F smaller than SN1, activity propagation is ensued as R 0 > 1. However, once the system reaches the high activity steady state, R 0 is maintained at 1 (Fig. 2b top and middle; Appendix D), because mass conservation limits the number of susceptible molecules. The system loses the transmission potential for γ DN F values higher than SN2.
Equivalently to the macroscopic case (Figs. 1d to 1f), the phase space trajectories showed a clear attraction towards a single- (Fig. 2c top) or two coexisting attractors for different initial conditions (Fig. 2c middle) and different γ DN F values. In the latter case, dynamical switching between the two steady states was observed. In the vicinity of SN2 however, the system exhibits critical behavior R 0 →1 (Fig. 2b, magenta region). The proximity of R 0 to 1 together with the diffusion-induced spatial inhomogeneities in P DN F,T /R T concentration can effectively increase the local transmission potential above the critical value 1, thereby generating local pockets of active receptor that transiently sustain and further propagate this state across the surface. This results in interchanging periods of inactivation and re-activation bursts around the ghost attractor state (e.g. green profile in Fig. 2c, lower-right panel), manifested as prolonged receptor activity before the system settles to basal state ( Fig. 2c bottom equivalent to Fig. 1f).

IV. TRANSIENT VS. SHORT-TERM MEMORY: SIGNAL INTEGRATION AND OPTIMAL INFORMATION PROCESSING
Coming back to the macroscopic description of the system, the transient memory does not only enable plasticity and robustness in receptor responses, but it is also a unique manifestation of a dynamic memory. In other words, the total duration of the transient memory in receptor activity depends on the previous stimulus history, thus enabling signal integration (Fig. 3a, magenta). This feature arises from the metastability of the ghost state where the phase space trajectory gets transiently captured after stimulus removal, thus permitting responsiveness to subsequent signals. Equivalent behavior was observed also for the neuronal model (Fig. 3a, inset).
The stable attractor realization of a simple shortterm memory (reversible bistability) on the other hand emerges through the hysteresis of signal activa-tion/deactivation levels (Fig. 3b). This memory is not manifested through prolonged temporal activation as for the transient memory (Fig. 3c, magenta), shown through step-wise input modulation (Fig. 3c, orange). This in turn excludes processing temporal dependencies that are inherent in time-varying signals, such that signal integration does not take place (Fig. 3a, orange). The other stable attractor realization, the long-term memory (irreversible bistability, Fig. 3c, green), does not allow responsiveness to upcoming signals to take place (Fig. 1e).
This crucial feature of dynamic memory is further complemented with additional information processing capabilities to enable optimal computation of time-varying signals. In the case of cell surface receptor networks for example, organization at criticality also endows cells with robust sensing via receptor activation in a switchlike manner (Fig. 3d), as a consequence of the topological phase space transitions (Fig. 1f, i →ii ). Even more, the dynamic range of the response amplitude rapidly increases when transiting from the monostable towards the bistable regime, with a clear peak at the SN2 bifurcation (right to left, Fig. 3e). Using the single-molecule reaction-diffusion simulation framework, it can be additionally demonstrated that robustness to noise is also optimal at criticality: the probability for spurious activa- tion in absence of stimulus highly increases in the bistable region, whereas in the vicinity of the SN2 bifurcation, this probability is close to zero (Fig. 3f). These results therefore show that critical organization is crucial for optimal information processing of time-varying signals.

V. SELF-ORGANIZED POSITIONING AT CRITICALITY BY FLUCTUATION SENSING
We next investigated how such positioning can be realized for receptor networks. For this, we use the example of the proto-oncogenic epidermal growth factor receptor (EGFR), for which a critical organization between a monostable and bistable mode of operation was recently experimentally demonstrated [4]. The dynamics of EGF sensing in this case is regulated by a fourcomponent spatially-distributed network, where the interactions of EGFR (R) with three specific membraneassociated enzymes -protein tyrosine phosphatases via a double negative (P DNF , PTPRG), a negative feedback (P NF , PTPN2) and a negative regulation (P NR , PTPRJ) are coupled via the vesicular trafficking of the receptor (Fig. 4a). Ligand-bound EGF receptors (LR) promote autocatalytic activation of ligandless receptors (red arrows) [12,17,18,30], and thereby transfer information about the extracellular environment before they are internalized and degraded [12,13]. The vesicular recycling on the other hand, brings the internalized and deactivated ligandless EGFR back to the plasma membrane, thereby establishing the EGF signal processing network. Numerical simulations using a two-compartment model that takes the trafficking-induced redistribution of receptors explicitly into account (Eq. (E1) in Appendix E) showed that input-induced decrease in plasma membrane receptor concentration rapidly shifts the operation of the system into the monostable regime. This results in decreasing dynamic range of the activation response amplitude and loss of transient memory of the recent stimuli (Fig. 4b, compare blue to magenta lines). Thus, the EGFR concentration at the membrane, previously a determining bifurcation parameter, must be dynamically maintained at the critical transition. Two-parameter bifurcation analysis showed how the positioning of the SN bifurcation point depends on the total receptor concentration R T and its recycling rate constant k rec (Fig. 4c).   Therefore, a self-organizing mechanism by which k rec is up-regulated as a result of the decrease in total receptor concentration would effectively retain the system in the vicinity of the SN for several subsequent growth factor pulses (Fig. 4d). It should be noted however, that the SN positioning asymptotically approaches a minimal receptor concentration below which the receptor recycling rate can no longer sufficiently compensate for the loss of receptors from the membrane (Fig. 4c, dashed line). Additionally, the receptor recycling machinery may also impose an upper bound on k rec , further limiting the resetting capacity. Such a dynamically-maintained organization would require a mechanism for sensing receptor abundance to estimate the divergence from the saddle-node bifurcation point, and an actuating mechanism to translate this positioning into a corresponding recycling rate. Information about the amount of receptors, especially in noisy settings, could be encoded in the fluctuations of the ac-tivity state. The temporal signature of these fluctuations depends on the alignment of the nullclines -hence on the positioning of the system in parameter space [8], which is in this case determined by the concentration of receptors. The dominant frequency in basal EGFR activity fluctuations that we estimated as an average of multiple stochastic activity profiles (Appendix G) was lowest around the SN bifurcation [24] (Fig. 4e). It is therefore possible to postulate an actuating molecular mechanism where cells use an effector protein downstream of EGFR that is embedded in a simple low-pass filter to up-regulate EGFR recycling when low frequency fluctuations of EGFR activity are present and thereby maintain positioning in the SN bifurcation vicinity (Fig. 4f). The previously identified characteristics of Akt [31][32][33], a serine-threonine kinase downstream of EGFR, qualify it as a candidate that could induce this prolonged maintenance of the system at the critical transition. However, once the system escapes the low-frequency valley around the SN and enters the monostable region, the cell will rely on EGFR synthesis to re-establish the organization at criticality in the long term. This example thus demonstrates that even a simple low-pass filter can help cells to maintain positioning at criticality from which optimal information processing capabilities in growth factor sensing emerge.

VI. DISCUSSION
The physics of how cells perceive time-varying chemical signals is a fundamentally important question that dates back to the work of Berg and Purcell [34]. In the past decades, the limits of biochemical sensing have been explored using equilibrium and non-equilibrium descriptions of sensing through ligand binding/unbinding events to stationary receptors [34][35][36][37][38][39]. The ability of most receptors to detect low ligand concentrations however necessitates high affinity binding [40][41][42], limiting the detection of rapidly varying signals. Furthermore, the localization of many cell surface receptors is dynamically maintained through vesicular recycling, whereas upon ligand binding, unidirectional internalization and receptor degradation occurs [12,13,43,44]. Although a rapid removal of ligand-bound receptors from the cell surface would in principle enable sensing of non-stationary cues via the remaining ligandless receptors, it would not allow for signal integration and robust signaling responses. From the dynamical system formalism it follows that the topology of the receptor network determines the dynamical possibilities of the system. However, as we showed here, the optimal positioning in parameter space is what enables biochemical sensing of non-stationary environments via robust but plastic non-linear activity receptor responses. Such a critical self-organized behavior [45] resembles the one of the mammalian hearing system [46], flocks of birds [47] or biofilms [48].
We demonstrated that processing non-stationary inputs in real-time as ubiquitous during chemical sensing, necessitates information to be encoded through metastable states and interpreted using phase-space trajectories. The minimal requirements for such processing mechanism are met by the existence of a ghost of a saddle-node. This mechanism can be realized for various intrinsic dynamics and network sizes, as we have shown for two-component cell surface receptor network and single LIF neuron. In both cases, the dynamic transient memory that emerges from the critical organization allows for history-dependent signal integration and generally, optimized information processing -computational features that were previously only attributed to largescale neuronal networks [49][50][51][52][53][54].
Conceptually, the framework of information processing with state-dependent trajectories thus suggests that the notion of computation with stable attractors [55,56], resembling the legacy of von Neumann and Turing [6], likely should be adapted for cellular processing of nonstationary signals. In this case, the dynamic transient memory and thereby related information processing features that emerge at criticality ensure optimal balance between stability and overall responsiveness to upcoming cues, as pervasive for systems that operate in a continuously changing environment.

ACKNOWLEDGMENTS
The authors thank Philippe Bastiaens for numerous discussions and suggestions that were crucial for the development of this work, as well as for critically reading the manuscript.
A.K. conceptualized the study. A.S. performed most of the analytical derivations, numerical simulations and bifurcation analysis with help of A.P.N. and A.K. A.S. and A.K. interpreted the results and wrote the manuscript with help of A.P.N.
The authors declare no conflict of interests. (1) we model a minimal network motif that exhibits bistability, a double-negative feedback, using law of mass action. Both, the receptor, as well as the deactivating enzyme have active (R a , P DNF,a ) and inactive (R i , P DNF,i ) states, and their state transition rates are described by the model equations. Therefore, mass is conserved in the system and the total protein concentrations of both species (R T , P DN F,T ) are constant parameters. This allows R i = 1 − R a and P DN F,i = 1 − P DN F,a , expressed as fractions from the total protein concentrations. Since the fraction of ligand-bound receptors (LR a ) is mapped to the cell from the ligand concentration in the environment [4], it is considered as an input parameter in the system.
Autonomous, autocatalytic and ligand-bound-induced activation of ligandless R i ensue from bimolecular interactions with distinct rate constants α 1−3 , respectively. Other parameters: k 2/1 -P DNF activation/inactivation rate constant ratio, k 1 -kinetic constant that does not influence the steady state values of the system,β DN F = β DN F R T /k 1 -receptor-induced regulation rate constant of P DNF ,γ DN F = γ DN F P DN F,T /R T -specific reactivity of the enzyme towards the receptor, thus proportional to the local effective P DN F,T /R T ratio. In the analysis we refer to changes of P DN F,T /R T when numericallyγ DN F is varied.
For simulations with growth factor pulses in Figs. 1d to 1f and Fig. 3a, binding/unbinding of ligand to modulate LR a was introduced. Thus, −k on R a L T + 1 2 k of f LR a was added as additional term to the differential equation of R a , and the dynamics of LR a was modeled with For the simulations in the Supplementary Video [22] a stochastic dif-ferential equation model was constructed from Eq. (1) by adding a multiplicative noise term σX i (1 − X i )dW t , where σ = 0.05, dW t is the Brownian motion term and X i (1 − X i ) is the state-dependent function for each variable i that accounts for mass conservation and normalization of the variables. The model was solved with ∆t = 0.01 using the Euler solver from the Financial Toolbox in MATLAB. The parameters corresponding to Figs. 1a to 1f and Figs. 3a to 3e are: α 1 = 0.0017, α 2 = 0.3, α 3 = 1.0,β DN F = 36.0558,γ DN F ∈ {2.5, 2.957, 3.5, 4.3} (bistability, memory, reversible bistability, monostability), k 1 = 0.01, k 2/1 = 0.5, R T = 0.8, k on = 0.003, k of f = 0.01668. The L T amplitude during the pulse was set to produce 0.15 of LR a in steady-state.

Appendix B: Quasi-potential landscape computation
The numerical computation of the quasi-potential landscapes, corresponding to the phase space diagrams in Figs. 1d to 1f was conducted using an approach adapted from the one in [57]. Multiple trajectories were calculated starting from initial conditions distributed on a grid in the phase space. Rate of change in the quasipotential (initiated arbitrarily to 0) was calculated by in every iteration and the quasi-potential was integrated by the ODE solver together with the system trajectory integration. Quasipotentials of trajectories converging to the same attractor were aligned to match at the steady-state level. Quasi-potentials of different attractors were subsequently aligned at the initial points of the neighboring trajectories that converge to the different respective attractors, i.e. at the separatrix points. Additionally, neighboring separatrix pairs were weighted by the angle between their derivatives (θ), according to the formula: 1 2 (1 − cos θ). This gives greater weight to diverging pairs, effectively aligning the separatrix quasi-potential values at the saddle steady state point. The quasi-potential landscape at every point in phase space was finally estimated by interpolation from the aligned quasi-potential values of all of the trajectory points.

Appendix C: Model of a leaky integrate-and-fire (LIF) neuron
We used a well-established model of a LIF neuron [11], described with the following set of equations: where V is the membrane potential and n -activation variable for K + .
Appendix D: Estimation of the basic reproduction number using single-molecule reaction diffusion framework We considered a two-dimensional domain representing the plasma membrane containing reacting and diffusing single molecules. The spatial coordinates of the molecules were updated using Brownian dynamics and time was discretized to intervals of length ∆t. First-order unimolecular reactions occur spontaneously with proba-bilityk∆t, wherek is the intrinsic reaction rate constant. Second-order bimolecular reactions on the other hand are modeled using the Doi method [58], following the Smoluchowski single-particle framework for describing diffusion influenced reactions [59]. An interaction takes place between two molecules that have diffused within a proximity distance σ of each other, and a reaction ensues with a probabilityg∆t, whereg is the microscopic bimolecular reaction rate constant. σ is of order of the molecule radius. In the rare event that a substrate molecule is in proximity of n > 1 other enzyme molecules, reaction takes place with probability 1 − (1 −g∆t) n , assuming any of the enzyme molecules can affect the state of the given substrate molecule. Formation of the product proceeds immediately upon successful bimolecular enzyme-substrate interaction, i.e. the state of the substrate molecule is directly changed. To model the interactions between R and P DNF (Fig. 1a), we assumed that both particles diffuse across the 2D domain with equal diffusion rates D R = D P DN F = D = 0.1 µm 2 /s. The interaction radius σ was set to 2ρ, where ρ = 10 nm is the molecule radius. 500 receptor molecules and 700 P DN F molecules were randomly deployed on a 5 µm × 5 µm surface and allowed to diffuse using Brownian dynamics with periodic boundary conditions. ∆t was set to 1 × 10 −3 s to ensure that √ 4D∆t ≤ σ, i.e. any interaction between two molecules that come in proximity is detected, and also to ensure that no reaction probability is greater than one. The state transitions of R and P DNF occur in accordance with the macroscopic description -Eq. (1), in absence of external input (LR a = 0). The microscopic rate constants are therefore proportional to the ones in our main ODE model:α 1 = 0.0017/(σ 2 π),α 2 = 0.3/(σ 2 π),β DN F = 0.360558/(σ 2 π), k 1 = 0.005/∆t,k 2 = 0.0025/∆t. They were set to produce faster kinetics, due to numerical and data storage constraints. By analogy to the macroscopic bifurcation analysis (γ DN F ∝ P DN F /R T ),γ DN F was varied between 0 and 0.92/(σ 2 π) (P DN F,T and R T were kept constant) to modulate the specific reactivity of P DNF,a towards R a . To calculate the basic reproduction number R 0 , the number of substrate receptor molecules R 0,j (t) that each R a,j molecule successfully activated within its activation lifetime was recorded, after that molecule has been previously activated itself at time t. R 0 was calculated as an average of these figures within a certain time interval. Theoretically R 0 depends on the probability of activating a susceptible receptor molecule and the duration of the activity lifetime, in analogy to the SIS epidemic spreading model. Neglecting the effect from autonomous activation of R (with low rate α 1 ), we arrive at the following definition R 0 ≡ α2R T (1−Ra) γ DN F P DN F,T P DN F,a . From Eq. (1) it is straightforward to determine the basal activity stable steady state as R a = 0, P DN F,a = k1 k1+k2 , and thus R 0 = α2R T (k1+k2) γ DN F P DN F,T k1 . By employing linear stability analysis we could determine that the condition for stability of this steady state is indeed R 0 < 1. On the other hand, it follows readily from the first equation that R 0 = 1 around the high activity stable steady state, where R a = 0. There we could also find that γ DN F = α2R T 2 k1P DN F,T (1−R a ) k1+k2 R T β DN F +R a = C 1 (1−R a )(C 2 +R a ), thus there is approximately a quadratic dependence betweenγ DN F and R a . This form was used in Fig. 2b, middle, to estimate the bifurcation profile: C 1 = 1.388, C 2 = 0.385, which translates toγ DN F,SN 1 = 0.5344, γ DN F,SN 2 = 0.66. To estimate the extensively occupied high and basal receptor activity states from the trajectories in R a -P DNF,a phase space (Fig. 2c), Gaussian mixture distribution was fitted to the data with two components and data points were pruned iteratively with 90th percentile cut-off until convergence. The trajectories, as well as temporal profiles on Fig. 2c were visualized with downsampling to every fifth frame due to size constraints, without loss of visual quality.
Appendix E: Compartmental model of spatial-temporal EGFR activation regulation The experimentally derived EGFR-PTP network [4] was implemented using a two-compartmental model that includes explicitly the vesicular trafficking between the plasma membrane and the endosomal compartments (Fig. 4a), as described using the following system of ODEs: P DNF (PTPRG), P NF (PTPN2), P NR (PTPRJ) represent the major protein tyrosine phosphatases that regulate EGFR (R, LR) activity,γ X = γ X P X,T /R T -specific reactivity that each P X ∈ {P DNF (PTPRG), P NF (PTPN2), P NR (PTPRJ)} has towards EGFR and is therefore proportional to local effective P X,T /R T ratio. k in , k rec , k deg denote receptor internalization, recycling and degradation rate constants, respectively, i,a -inactive and active species, E, P M -endosomal and plasma membrane species. For the self-organizing criticality model (Fig. 4d), the dynamically maintained k dyn rec = R T −R T ,asymp R T (1−LR E i )−R T ,asymp k rec , where R T,asymp = 1.086 is the lower bound asymptotic value of R T in dependence to k rec (dashed line, Fig. 4c). Saturation level of 2.5 for the multiplier term is also assumed, beyond which the recycling rate can no longer increase. Parameters: β DN F =β DN F , γ DN F = 3.0, γ N F = 3.0, γ N R = 0.001, k in = 0.02, k rec = 0.042, k i in = 0.2, k deg = 0.2, R T = P DN F,T = P N F,T = P N R,T = 1.0. Other parameters are the same as in Fig. 1.
the single cell dose-response data from where the topology of the sensing network (Fig. 4a) was derived. We convert from dimensionless time to minutes by equating the EGFR phosphorylation kinetics and duration in the simulations using the kinetic parameters to the experimental values in [4]. The parameters for the microscopic dynamics in the single-molecule reaction diffusion simulations were set to scale the macroscopic ODE parameters and set to produce faster kinetics due to numerical reasons, as described in the corresponding section.