Publisher Correction: Trapping Phenomenon Attenuates the Consequences of Tipping Points for Limit Cycles

A correction has been published and is appended to both the HTML and PDF versions of this paper. The error has not been fixed in the paper.

with a large transport of heat to the Northern latitudes and the other one corresponding to a shut-down of the circulation consequently ceasing the heat transport 25 . As a part of the thermohaline ocean circulation, local deep ocean convection is also vulnerable to a shut down 26,27 . Both processes, the shut-down of THC on a global scale and deep-ocean convection on a local scale would have a large impact on the climate in the Northern hemisphere leading finally to a cooling in Northern-and Western Europe. The Indian Monsoon is expected to become more wet or more dry depending on which of the processes responsible for such changes like increasing albedo due to aerosol concentration or stronger El Nino's, respectively, are dominant in the future 24 . Another highly debated tipping point relates to the tendency of the arctic ice sheets to become thinner and finally to lead to an ice-free state in summer 28,29 . Finally, we mention the different approaches to study the recurrent switchings of ice ages and warmer climates before the Holocene, which are attributed to several stable states and transitions between them 30 .
Different scenarios are discussed in literature that lead to such critical transitions. On the one hand, changes in the environmental forcing, e.g. atmospheric temperatures or altered precipitation patterns can induce such transitions by reaching a critical threshold at which one of the states loses its stability and the system switches to another state. In mathematical terms, this scenario is related to a bifurcation; the combination of two such bifurcations often comprises a hysteresis 31 allowing for a switching between two alternative states, when a control parameter is varied 15 . When two stable states coexist, then a switching between them is mediated by fluctuations leading to noise-induced transitions 32,33 .
Due to the increasing concerns about such critical transitions on our planet earth, there is an urgent need to identify the approach of a regime shift or a tipping point before its occurrence. Such early warning signals, if applicable, can be used to anticipate the transition and to take measures to slow down the approach in the worst case or to avoid it if the expected alternative state is for some reason undesired. Such avoidance might not always be possible, particularly not in the climate system, but early warnings could be used to take political actions. During the last decade several methods have been developed to gain more insights into how to predict abrupt changes in the system dynamics, induced by its nonlinearity. One of the earliest measures identified is related to the time which is needed to respond to perturbations. While far away from the critical threshold, such perturbations die out quite quickly, this damping becomes significantly slower in the neighborhood of a threshold 32,33 . The perturbation applied can be either a single perturbation or some noise which is inevitable in experiments and in natural processes. In case of a single prescribed perturbation, this measure is easy to implement in experiments and therefore widely used in quantifying the distance to the threshold value.
In case of a noisy system the approach of the critical threshold can be quantified by an increasing variance and autocorrelation 33 . As a second effect noise leads to an irregular switching process between the two (or more) alternative stable states. This process is called flickering 33,34 , attractor hopping 35 or chaotic itinerancy 36,37 depending on the context in which it is studied. Hence, a second indicator has been introduced measuring this flickering or hopping process, which occurs in a bistable (or multistable) region in parameter space in which two or more stable states coexist. It is important to note that the hopping dynamics depends on the noise strength. While a large body of work is devoted to the impact of additive noise, many processes in nature, particularly in ecosystems, are affected by multiplicative noise, which has only rarely been considered. It has to be emphasized that environmental noise in ecosystem dynamics is always multiplicative and plays by far the more important role 38,39 . However, most of the previous work on indicators for critical thresholds is restricted to additive noise describing the impact of fluctuations on physical processes in the climate system, but being of limited value for ecological problems.
Many bistable systems considered in nature possess two stable equilibria, i.e. the system is stationary. For the above mentioned example of shallow lakes, the water in the lake loses transparency by shifting abruptly from a clear to a turbid state when a threshold in the level of nutrients is reached. As a result of this eutrophication process, submerged plants dramatically disappear beyond a tipping point 15 . For the example of desertification, rainfall patterns are the essential environmental conditions determining the shift from a perennial vegetation via localized vegetation patterns to the state of bare soil. Moreover, taking the grazing pressure by livestock in the Sahel zone into account, leads again to a shift from a perennial to an annual vegetation 21 . However, in analyzing the regime shift in the respective ecosystems, the periodic input of nutrients and precipitation due to the seasonal cycle has been neglected, but could have an important impact too. The same argument applies to the analysis of physical processes in the climate system driven by periodic or quasi-periodic changes in the orbital parameters of the sun leading to a variation in solar insolation with periods of about 23,000, 44,000 or 100,000 years, the well-known Milankovitch cycles 40 . Particularly the latter are assumed to be the major drivers for the development of the ice-ages before the Holocene 30 . Recently, these periodic forcings affecting the Albedo of the earth are studied to evaluate the impact of this variation on the Arctic and Antarctic ice cover 29,28 .
In this paper we focus on tipping points and regime shifts of periodically forced systems. In this class of systems, abrupt changes of the dynamics at critical transitions do not occur between steady states (equilibria), but between oscillating states (limit cycles). Though at first sight, the hysteresis curves which are usually drawn to discuss critical transitions look similar, however, the dynamics is quite different. For this scenario we show that a transient structure, a so-called channel, occurring in the system's state space beyond the tipping point, creates a short-term dynamical regime with specific properties which attenuates the criticality of the transition. The smoothing of the transition is demonstrated by computing a finite-time measure of the twisting behavior (rotation property) of the state space surrounding trajectories while inside the channel domain. This measure indicates that, the trajectories passing the channel beyond the tipping point have residual system properties of the limit cycle destroyed at that tipping point. Hence, the channel acts as a "ghost" of the destroyed limit cycle, retaining system trajectories in a very similar fashion. This fact is shown by statistical analysis of the intervals of time spent by noisy trajectories in the neighborhood of the limit cycle (pre-tipping) and the channel (post-tipping). Furthermore, we attribute to the ghost state the inconclusive diagnostic provided by variance and autocorrelation in anticipating tipping points.
Let us now indicate the differences in the dynamical approach to deal with limit cycles instead of asymptotic equilibria. Figure 1 shows the typical bifurcation diagram used to explain the appearance and disappearance of the coexistence of two alternative states. In contrast to the usual diagrams, the lines denote here one point of a limit cycle instead of stationary points. Therefore, the y-axis does not show one coordinate of the stationary state of the system, but one coordinate of the Poincaré section, a special construction which is very useful in analyzing periodic solutions of nonlinear dynamical systems (see Methods section). The Poincaré section in a periodically forced system defines a stroboscopic map in which the system is always analyzed at the same phase of the forcing, i.e. at times t, t + T, t + 2T, … . Hence, a limit cycle corresponds to a fixed point in this stroboscopic map which makes the similarity between Fig. 1 and the well-known sketches of bistability in the stationary case obvious.
There are two saddle-node or fold bifurcations of limit cycles denoted by F 1 and F 2 , at which two limit cycles, a stable and an unstable one, emerge or disappear. Crossing those tipping points the dynamics will change dramatically. Increasing the parameter value p, the continuation of the limit cycle on the upper branch will stop abruptly at F 2 and switch to a periodic behavior corresponding to the lower branch, while decreasing the parameter p and continuing the lower branch will result in a transition to the upper one at the critical threshold F 1 . Our main focus lies in the analysis of the regions around those critical transitions. Firstly, we address the question to what extent the usual criteria of critical slowing down and flickering will signal the approach to the transition (blue region). Secondly, we show that the critical transitions are hidden due to particular structures in state space, so-called channels, which appear in the neighborhood of the fold bifurcations of limit cycles, preventing a clear identification of the critical transitions.

Results
To illustrate our results, we employ a paradigmatic model system, the well-known Duffing oscillator 41 and apply a periodic forcing with amplitude A and frequency ω. In mathematical terms, this simple dynamical system reads: The parameter d is the amplitude of the system damping. The parameter σ controls the noise intensity given by the stochastic forcing ξ(t). The function ξ(t) represents the usual zero mean and unit variance noise with We use a fourth-order Runge-Kutta method to integrate Eq. (1), in this process, the time is measured as a function of the period of the external forcing, i.e., T = 2π/ω.
In a certain parameter range, the system described by Eq. (1) exhibits a generic scenario of bistability between two different limit cycles, i.e. two stable periodic solutions exist separated by an unstable one of saddle character. The corresponding bifurcation diagram is shown in the upper panel of Fig. 2. Though this diagram looks very similar to the general diagram depicted in Fig. 1, it shows only a Poincaré section of the stable limit cycles occurring for the system described by Eq. (1). To characterize those limit cycles in more detail, we compute the generalized winding numbers (GWN) for each of them along the bistable parameter range, the results are depicted in the bottom panel of Fig. 2. In this panel, the GWN is represented by w ∞ , this measure quantifies the asymptotic twisting of the local in neighborhood of limit cycles, a better description of this measure is given in the Methods.
The bifurcation diagram of Fig. 2 (Upper) shows the dependence of the noise-free Duffing oscillator on the forcing amplitude A. Two stable limit cycles, S 1 (blue) and S 2 (yellow), coexist for a range of parameters bounded by two fold bifurcations of limit cycles at the points F 1 and F 2 (tipping points). So, the system is subject to catastrophic shifts, tipping points, as the parameter A reaches the points A 1c or A 2c . Let us now check whether the autocorrelation coefficient at lag-1 and the variance of the system indicate the approach of the critical transition and can serve as early warnings. To this end, we apply now noise to the system and show the resulting behavior  Fig. 3. We focus our analysis on the parameter region close to A ~ A 1c . Hence, in Fig. 3, we reverse the x-axis to better investigate the critical transition at the parameter A 1c . In Fig. 3 (Upper panel), the black line indicates the time evolution of a noisy trajectory with the parameter A varying in the same interval of the bifurcation diagram also indicated in this panel. In Fig. 3 (Middle panel), we show the autocorrelation at lag-1, a measure that usually increases with the approach of critical transitions of systems in equilibrium. In this panel, for limit cycles, we verify a sudden increase in the autocorrelation coefficient as soon as the noisy trajectory starts flickering between the stable limit cycles. But subsequently, it decreases as the system approaches the critical threshold and does not exhibit any significant change as the critical threshold A 1c is passed. Similar behavior is observed for the standard deviation of the noisy trajectory, shown in Fig. 3 (Bottom panel). Therefore, we find that the usual indicators of critical transitions between equilibria may not work for limit cycles. Instead we observe a continuous decreasing of the autocorrelation coefficient and the variance, not suitable to serve as an early warning signal. To explain this behavior, we investigate in more detail the state space structure occurring for parameters succeeding the fold bifurcation at A 1c .
Fold bifurcations of limit cycles are accompanied by the formation of channels in state space through which the trajectory has to go after entering it. To illustrate this behavior which has been first described by Manneville 42,43  in the context of intermittency in turbulence, we show in Fig. 4 the generic principle behind the formation of that channel. As outlined above, a limit cycle corresponds to a fixed point in the Poincaré section. In our case, one point x n in the Poincaré section is mapped onto the next point x n+1 by mapping the limit cycle stroboscopically every period of the forcing, so the fixed points are mapped onto the diagonal x n+1 = x n of the diagram x n+1 × x n shown in Fig. 4. In the bistable region we have three fixed points, two stable and one unstable separating the former two (Fig. 4, left panel). When the fold bifurcation is reached the stable and the unstable limit cycle merge and form an elliptic point (middle panel), while beyond the fold bifurcation a channel appears in state space through which the trajectory moves when it comes close to the region in state space where previously the two limit cycles have been located. Without noise, the trajectory would finally converge to the only stable limit cycle left in the system, denoted by S 1 . Due to the noise, the trajectory is kicked back to the channel from time to time and moving through it again and again. As a consequence of this behavior we observe even beyond the fold bifurcation, that the dynamics returns to the "ghost" of the limit cycle manifested as the channel. The resulting dynamics contains phases where the trajectory is close to the "ghost" and phases where is comes close to the only stable limit cycle but kicked away again by the noise. This way, the flickering dynamics goes on even though the system is well beyond the critical transition. For the very same reason the widely used indicators for critical transitions such as lag-1 autocorrelation function, variance as well as flickering can not signal the approach to the critical transition and the shift or tipping points happens with no warning. In our case, the characteristics of the critical slowing down indicators resemble the case of a smooth transition as analyzed in ref. 44. Additionally, we note that the dynamics before and beyond the critical transition is essentially the same, characterized by the hopping between the two stable states before and between the single stable state and the "ghost" beyond the tipping point. This behavior is generic and will occur for all fold bifurcations of limit cycles forming a channel after the bifurcation.
Let us now discuss the post-tipping behavior in more detail. To demonstrate that indeed the channel is the most essential structure in the state space deforming the noisy system beyond crossing the critical threshold, we analyze the scaling behavior of the length of the transient time with the distance from the threshold, Therefore, we define ε as a parameter measuring the distance from the critical threshold A 1c , i.e., ε = A 1c − A. Then, as a function of this distance ε, we measure the transient time, τ(ε), for trajectories to reach the remaining stable limit cycle (yellow). For these trajectories, we choose a set of initial conditions in the state space region previously occupied by the basin of attraction of the limit cycle S 2 (blue) destroyed in F 1 . In Fig. 5, we show the results for an ensemble of 300 random initial conditions for each distance ε. We find the time that trajectories spend to leave the channel scales as a power-law with the distance ε. The characteristic exponent is equal 0.5, and hence, it corresponds to the value known from the type I intermittency 42,43 . Hence, the characteristic time τ(ε) enables us to quantify the influence of the channel in the time evolution of trajectories as a function of the parameter distance ε.
This equality verifies that trajectories starting in the former basin of attraction of the limit cycle destroyed in the fold bifurcation are in fact being trapped in the channel associated to this bifurcation for a characteristic time, τ(ε). In order to obtain the twisting behavior of trajectories just during the time trapped in the channel, we introduce a finite-time version of the winding number (FTWN) represented by 〈 w(t, ε)〉 . A complete description of this definition is given in section Methods. In the diagram shown in Fig. 6, the color code indicates the FTWN given by 〈 w(t, ε)〉 , in the y-axis we represent the time evolution, t, in units of the period of the forcing, while in the x-axis we show the distance ε from the bifurcation point. The red line represents the function τ(ε) obtained from the adjustment in Fig. 5 for the characteristic time for trajectories to cross the channel. We observe that regardless of the parameter distance ε, the finite-time winding number has a defined value equal to 7.0 (blue in Fig. 6) for times lower than the corresponding τ(ε). Hence, from Fig. 6, we conclude that the post-tipping trajectories, while crossing the channel, conserve the twisting behavior (rotation properties) of the stable limit cycle destroyed in the tipping point.
In the following, we confirm the existence of the residual twisting behavior of the destroyed limit cycle by obtaining the FTWN of sets of initial conditions crossing the channel. Firstly, we choose the parameter A such Figure 4. (a) Two fixed points of the node type (S 1 and S 2 ) are coexisting with a saddle. The black arrows indicate how initial conditions dynamically behave in the neighborhood of each fixed points. (b) As a system parameter is varied the node S 2 and the saddle collide forming an elliptic point E (saddle-node or fold bifurcation). (c) As the parameter crosses the critical bifurcation parameter, the initial conditions (arrows) that used to belong to the attraction domain of the extinct node S 2 are now converging to the node S 1 through a channel formed in the mapping.
that the dynamics will take place in the channel, i.e., A 1c minus a small distance ε = 0.0095, then we compute the FTWN during the time τ(0.0095) = 33.34 (T) for a grid of initial conditions. Attributing different colors to the FTWN obtained for the trajectories corresponding to each initial condition, we clearly distinguish, in the grid of Fig. 7(a), two types of dynamic behavior. (i) the FTWN corresponding to trajectories that cross the channel (initial conditions in blue, fast twisting) and (ii) the the trajectories converging directly to the remaining stable state (initial conditions in yellow, slow twisting). In order to compare the twisting properties of the channel, measured by FTWN, to the twisting behavior around the stable states in the bistable region, we characterize the twisting around the two stable states, S 1 and S 2 , by the asymptotic generalized winding number (GWN). Figure 7(b) shows those twisting properties of trajectories starting on the same grid as in Fig. 7(a) but computed by the asymptotic (GWN) for a forcing amplitude A before the tipping. We notice the similarity between twisting of trajectories crossing the channel beyond tipping (blue in Fig. 7(a)) and around the stable state before tipping (blue in Fig. 7(b)).
To illustrate further that observations of the system's trajectories are insufficient to determine whether the system is bistable (pre-tipping) or has a dynamical channel (post-tipping), we show in Fig. 8(a) the temporal evolution of a trajectory of the noisy the Duffing oscillator as the parameter A increases with time in the same interval as in Fig. 8(a). We notice that, even after the limit cycle S 2 marked in blue disappears in F 1 , the noisy trajectory (black line) is still flickering into the state space region previously occupied by the extinct limit cycle around = −  x 4. This becomes even more obvious when comparing two noisy trajectories with fixed forcing at an amplitude A in the bistable region (pre-tipping) to a trajectory with a forcing amplitude beyond the tipping point (red line in Fig. 8(a). In a statistical sense those two trajectories are indistinguishable, indicating that the pre-tipping and the post-tipping behavior are very similar, with flickering between two distinct state space regions of S 1 and S 2 or the "ghost" of S 2 respectively.
As a consequence, time series as the main window to observations in nature, would show the flickering phenomenon before and beyond the tipping points making the transition in the observed data to appear smooth  instead of abrupt. In order to verify this statement, we investigate the intervals of time, θ, that a noisy trajectory elapses in the neighborhood of the stable state (before the tipping point), and in the channel (beyond the tipping point). The idea behind this study is to extend the notion of escape times to characterize the dynamics beyond the tipping point. In bistable systems one usually computes the mean escape time or mean first passage time to identify the stability of each stable state in a stochastic sense. While for systems possessing a double well potential, it is possible to compute those escape times analytically 45 , one has to rely on numerical estimations for arbitrary multistable systems 35 . Though the vast majority of nonlinear dynamical systems do not possess a potential, the scaling of the escape rates remains valid. Specifically, in Fig. 9, we obtain the distribution of time intervals spent by trajectories in the neighborhood of the stable states and in the channel. The time interval, θ is also expressed in units of the period of the forcing T.
Specifically, we show in Fig. 9 the distributions of the time intervals spent by trajectories in the neighborhood of the stable limit cycle S 2 which will go extinct at the tipping point ( Fig. 9(a)) with the distribution of those time intervals spent in the neighborhood of the "ghost" of S 2 beyond the tipping point ( Fig. 9(b)). Both distributions are exponential distributions, so that the probability density function can be approximated by  where μ = 〈 θ〉 > 0 is the mean value of the distribution. While in the bistable parameter region the mean time spent close to the limit cycle S 2 is 〈 θ〉 ≈ 10 periods of the forcing, it is only slightly shorter (〈 θ〉 ≈ 6 periods of the forcing) beyond the tipping point. However, the density function for the dynamics close to the channel is narrower and higher than in the bistable region, indicating that the shorter intervals of time are more frequent. Hence, even for the parameter A lower than A 1c (beyond tipping point), the frequency with which trajectories visit the neighborhood of the extinct state is not zero, i.e., the flickering phenomenon still occurs after the tipping point. It means that even after the extinction of the limit cycle in A 1c , the state space channel keeps retaining trajectories, avoiding their abrupt definitive transition to the unique survival stable state.
To emphasize that, the characteristics of the dynamics changes smoothly and not abrupt when crossing the tipping point, we show in Fig. 9(c) the changes in the distribution function when decreasing the forcing A from the bistable to the monostable region. In this figure, the color code indicates the probability densities for each parameter A shown in the x-axis. We observe that the distribution of time intervals smoothly changes as the tipping point is approached and passed, indicating that there is no abrupt transition crossing the threshold. However, as the parameter A is passed through the tipping point, a considerable increasing in the density of time interval values around the expected value is observed making the distribution narrower for parameters well beyond the tipping point.

Discussion
In summary, we address tipping points of systems subjected to periodic external forcing. The asymptotic solutions of this class of systems inherently settle into oscillating stable states (limit cycles), a more complex dynamics than the stable steady states (equilibria) for which the tipping points have been widely studied. In nature, the most noticeable occurrences of such oscillating attractors are found in ecology and climate sciences where periodic and quasi-periodic variabilities arise from external factors such as seasonality and astronomical forcing. Here, for a generic periodically forced system that generates such oscillations, we consider the typical hysteretic scenario to investigate tipping points, i.e., a bistable parameter region where the tipping is represented by fold bifurcations of limit cycles rather than steady states. As the parameters are varied and the system reaches a fold bifurcation, in which a stable limit cycle is destroyed leaving a transient structure, a so-called channel, in the state space of the system. Hence, for parameters beyond this tipping point, the channel gives rise to a short-term dynamics which possesses similar properties than the destroyed limit cycle and can therefore be attributed to a "ghost" of the latter. We find that a residual dynamical property of the limit cycle destroyed in the tipping point, namely its twisting behavior, occurs in the short-term dynamics for parameters in the post-tipping region. This finding indicates that the short-term behavior carries dynamical information of the destroyed oscillating stable state.
For system parameters fixed in the post-tipping region, we obtain the time evolution of the system subject to a stochastic noise. With this, we show that the "ghost" attractor retains systems trajectories in a very similar fashion of the stable limit cycle destroyed in the tipping point. Additionally, by obtaining the statistics of the time intervals that noise trajectories spend in the neighborhood of the stable limit cycle and in the neighborhood of the "ghost", we find that the PDFs of waiting times in both regions have the same exponential profile and do not differ much in their expected values. Therefore, the "ghost" dynamics plays an essential role in attenuating the critical transition in a way that it may be seen as a smooth transition when trying to diagnose it from real-world data. Hence, none of the well-known methods like autocorrelation function, variance or flickering are suitable to identify this particular transition properly. The emergence of the "ghost" masks the transition until the system is well beyond the tipping point and makes it to appear smooth instead of catastrophic. Finally, we would like to note that, despite of our choice of dealing with periodically forced systems due to their importance to ecological and climate systems, the "ghost" state masking the critical transition should also be observed in parameter regions succeeding fold bifurcations of limit cycles in autonomous systems.

Methods
Poincaré Section. As we consider systems whose asymptotic behavior are limit cycles, the final dynamics are oscillations rather than equilibria. Bifurcation analysis are performed by defining a Poincaré section, which usually is a hyper-surface arranged transversally to the limit cycle, where the whole system dynamics is described by a discrete system. Letting f be the function that describes the intersections of the limit cycle with the section, for trajectories in three dimensional space, results that = n n n n 1( 1) 2( 1) 1 2 , where x 1n , x 2n are the coordinates of the n th crossing. Consequently, on the surface of section, limit cycles are represented by fixed points of f. Then, states, as shown in bifurcation diagrams such as of Fig. 1, are defined in the surface of section, and in case the section is chosen to be a plane, they are denoted by the (x 1 , x 2 ) plane coordinates.
For the Duffing oscillator described by Eq. (1), a suitable Poincaré section is the so-called T-shift. The dynamics over the section is represented by discrete variables ( x nT ( ), x(nT)) defined as the solution pair ( x, x) collected every period, T = 2π/ω. Fold Bifurcation of Limit Cycles and State Space Channel. In case of limit cycles, a risky bistable configuration occurs when two stable limit cycles are coexisting with one unstable cycle of saddle type. The emergence of a dynamical channel at this scenario can be described on a suitable Poincaré section. We show in Fig. 10(a) that the stable limit cycles yield two fixed points of the node type, S 1 and S 2 in the surface of section, while the unstable limit cycle generates a fixed point of saddle type. The stable manifold of the saddle separates the initial conditions attracted by each node (blue and red in Fig. 10(a)). As the control parameters are varied approaching the fold bifurcation that delimitate the bistability region, such as F 1 and F 2 in Fig. 1, one of the node fixed points approaches the saddle. When the system is set to the parameters at the fold bifurcation point, for instance, the node S 2 and the saddle collide, and they both disappear forming an elliptic fixed point denoted by E in Fig. 10(b). For post-bifurcation parameters, the schematic of the Poincaré section is shown in Fig. 10(c), the node and the saddle no longer exist in the Poincaré section. However, trajectories starting in the space state region, which used to be the basin of attraction of the destroyed node S 2 , converge to the remaining node but not before being attracted by the stable manifold of the unstable elliptic point. Effectively, the system's trajectory behaves as if there existed a channel constraining the trajectory and leading it to the remaining stable state.
The occurrence of these dynamical channels related to fold bifurcations of limit cycles has been first discussed by Pomeau and Manneville, and has been argued to be the mechanism responsible for the laminar phase in the type-I intermittency scenario 42,43,46 . However, in type-I intermittency, chaotic bursts re-inject the trajectory in the dynamical channel. The trajectory spends long time intervals to cross the channel (the laminar phase) eventually escaping to the chaotic phase space region, producing another chaotic burst. In this work, there is no chaotic process to re-inject the trajectory into the channel, so we introduce a Gaussian noise which resets the trajectory to a random configuration belonging to the basin of attraction of the stable state extinct in the fold bifurcation. This procedure successively ejects the trajectory off the neighborhood of the survival stable state, forcing it to successively cross the channel along its time evolution. Regardless of the mechanism used to re-inject the trajectory Figure 10. Two dimensional representation of the bifurcation scenario. (a) Two fixed points of the node type (S 1 and S 2 ) are coexisting with a saddle. The black lines indicate how initial conditions dynamically behave in the neighborhood of the fixed points. The colors indicate the domain of attraction of each node. (b) As a system parameter is varied the node S 2 and the saddle collide forming an elliptic point E (saddle-node or fold bifurcation). (c) As the parameter crosses the critical bifurcation parameter, the initial conditions that used to belong to the attraction domain of the extinct node S 2 are now converging to the node S 1 through a channel formed in the plane.