Early-warning signals for bifurcations in random dynamical systems with bounded noise
aa r X i v : . [ m a t h . D S ] A p r Early-warning signals for bifurcations inrandom dynamical systems with bounded noise
Christian Kuehn ∗ , Giuseppe Malavolta † , and Martin Rasmussen ‡ April 9, 2018
Abstract
We consider discrete-time one-dimensional random dynamical systemswith bounded noise, which generate an associated set-valued dynamical sys-tem. We provide necessary and sufficient conditions for a discontinuous bifur-cation of a minimal invariant set of the set-valued dynamical system in termsof the derivatives of the so-called extremal maps. We propose an algorithmfor reconstructing the derivatives of the extremal maps from a time seriesthat is generated by iterations of the original random dynamical system. Wedemonstrate that the derivative reconstructed for different parameters canbe used as an early-warning signal to detect an upcoming bifurcation, andapply the algorithm to the bifurcation analysis of the stochastic return mapof the Koper model, which is a three-dimensional multiple time scale ordi-nary differential equation used as prototypical model for the formation ofmixed-mode oscillation patterns. We apply our algorithm to data generatedby this map to detect an upcoming transition.
Keywords:
Bifurcation, bounded noise, early-warning signal, fast-slow sys-tem, Koper model, mixed-mode oscillations, random dynamical system, set-valueddynamical system.
There has been a steadily increasing interest into dealing with sudden changesin the behaviour of dynamical systems, sometimes referred to as critical transi-tions and tipping points in the applied sciences. Many critical transitions canbe modelled mathematically using differential equations involving bifurcations .For instance, medical conditions such as asthma attacks and epileptic seizures ∗ Technical University of Munich, Faculty of Mathematics, Boltzmannstr. 3, 85748 Garchingbei M¨unchen, Germany † Department of Mathematics, Imperial College London, 180 Queen’s Gate, London SW7 2AZ,United Kingdom ‡ Department of Mathematics, Imperial College London, 180 Queen’s Gate, London SW7 2AZ,United Kingdom +
09, Kue13].In low-dimensional autonomous (time-independent) dynamical systems, crit-ical transitions have been extensively studied in the context of bifurcation the-ory, by explaining and classifying ways in which attractors can lose stability andgive rise to new types of behaviours, from simple transitions between stationarityand oscillations, to transitions between predictable and intrinsically unpredictable(chaotic) dynamics [CH96, Kuz95]. There is a central assumption of very slow(adiabatic) variation of parameters in the classical bifurcation theory, and theproblem of critical transitions in complex systems and their early-warning signalsmotivates the need to deal with transitions in nonautonomous (time-dependent)and random dynamical systems [KR11, Arn98]. Critical transitions in form of rate-induced tipping in nonautonomous dynamical systems have recently been studiedin [APW17, RS16]In this paper, we consider discrete-time random dynamical systems withbounded noise, which give rise to a set-valued dynamical system. We study bi-furcations of the set-valued dynamical system in form of discontinuous changesin the minimal invariant sets (which correspond to critical transitions), and weaddress the question of an early warning to predict such transitions. We pro-vide an algorithm that approximates the derivatives of the extremal mappings ofone-dimensional set-valued dynamical systems from a time series generated by thecorresponding random dynamical system with bounded noise. The use of deriva-tives of extremal maps as early-warning signals has not been established before.We show that the derivative is equal to 1 at a bifurcation point necessarily, andprovides a universal threshold for this reason, which does not depend on the typeof system or the type of bounded noise.We apply our theoretical results to the return map in the Koper model [Kop95].This model is a prototypical example exhibiting mixed-mode oscillations (MMOs)[DGK + Acknowledgements.
Christian Kuehn was partially supported by a LichtenbergProfessorship of VolkswagenStiftung. Giuseppe Malavolta and Martin Rasmussenwere supported by an EPSRC Career Acceleration Fellowship EP/I004165/1. Mar-tin Rasmussen was supported by funding from the European Union’s Horizon 2020research and innovation programme for the ITN CRITICS under Grant AgreementNumber 643073.
Consider the dynamics of continuous mappings f : R → R perturbed by additivebounded noise, given by the random difference equation x n +1 = f ( x n ) + ξ n , where { ξ n } n ∈ N is a sequence of i.i.d. random variables with values in a compact set K ⊂ R . The collective behavior of all future trajectories can then be studied viathe set-valued mapping F : R → K ( R ), defined by F ( x ) := { f ( x )+ y ∈ R : y ∈ K } ,where K ( R ) is the set of all nonempty compact subsets of R , equipped with theHausdorff distance h : K ( R ) × K ( R ) → R +0 .Note that for two nonempty sets A, B ⊂ R and x ∈ R , the distance of x to A is given by dist( x, A ) := inf {| x − y | : y ∈ A } , and the Hausdorff semi-distance of A and B is defined by d ( A | B ) := sup { dist( x, B ) : x ∈ A } . Then the Hausdorffdistance of A and B is given by h ( A, B ) := max { d ( A | B ) , d ( B | A ) } .In this paper, we study set-valued mappings F : R → K ( R ) without makingreference to a deterministic mapping f as above. We assume throughout thatset-valued mappings F : R → K ( R ) satisfy the following conditions. Hypothesis 1 (on the set-valued mapping F ) . We assume the following hypothe-ses:(H1) F ( x ) is an interval for all x ∈ R .(H2) There exists an ε > such that F ( x ) contains an interval of length ε for all x ∈ R .(H3) The so-called extremal maps f − : R → R and f + : R → R , given by f − ( x ) := min F ( x ) and f + ( x ) = max F ( x ) , (1) are continuously differentiable. Fundamental objects of interest for set-valued mappings are minimal invariantsets. 3 efinition 2 (Minimal invariant set) . A set E ∈ K ( R ) is called a minimal invari-ant for a set-valued mapping F : R → K ( R ) if F ( E ) = E and E does not containany proper nonempty subset with this property.The importance of minimal invariant sets is due to the fact that they are thesupports of stationary densities of the corresponding Markov semi-group [ZH07].To study bifurcations of set-valued mappings, we consider families of set-valueddynamical mappings { F α } α ∈ A , where A ⊂ R is open. Hypothesis 3 (on the family of set-valued mappings { F α } α ∈ A ) . We assume thefollowing hypotheses:(H4) The mapping ( α, x ) F α ( x ) is continuous for all α ∈ A and x ∈ R .(H5) The mapping α F α ( x ) is continuous for every α ∈ A , uniformly in x ∈ R . Let M α be the union of all minimal invariant sets of the set-valued mapping F α , α ∈ A . We are interested in discontinuous changes of M α under variation of α , which we will call discontinuous topological bifurcations . Definition 4 (Discontinuous topological bifurcation of minimal invariant sets) . We say that a family of set-valued mappings { F α } α ∈ A admits a discontinuoustopological bifurcation of minimal invariant sets at a point α ∈ A if the mapping α
7→ M α is discontinuous at α in Hausdorff distance.It will be convenient to define when a single minimal invariant set topologicallybifurcates. We say that a family of set-valued mappings { F α } α ∈ A topologicallybifurcates at a minimal invariant set E α of F α if for any neighborhood U of E α ,the mapping α U ∩ M α is discontinuous at α in Hausdorff distance.We aim at providing necessary conditions for a discontinuous bifurcation ofminimal invariant sets in this paper. As a first step towards this goal, we showthat a minimal invariant set satisfying certain properties on the derivatives of theextremal maps persists under continuous perturbations. Proposition 5.
Let { F α } α ∈ A be a family of set-valued mappings satisfying theassumptions (H1)–(H5). Suppose that there exists a minimal invariant set E α =[ e , e ] of F α such that (cid:12)(cid:12)(cid:12)(cid:12) d f − α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) < (cid:12)(cid:12)(cid:12)(cid:12) d f + α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) < , (cid:12)(cid:12)(cid:12)(cid:12) d f − α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) < , (cid:12)(cid:12)(cid:12)(cid:12) d f + α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) < , where f − α and f + α are the extremal mappings of F α as defined in (H3). Thenthere exists a non-bifurcating family of minimal invariant sets in a neighborhoodof α , i.e. there exists δ > such that for every α ∈ B δ ( α ) , there exists aminimal invariant set E α of F α , and { F α } α ∈ A does not topologically bifurcate atthe minimal invariant set E α of F α . Note that in [ZH07], a similar result is proved for random diffeomorphisms un-der the assumption that the minimal invariant set is isolated . An isolated minimalinvariant set E of a set-valued mapping F satisfies the property that for any openset W ⊃ E , there exists an open set U with E ⊂ U ⊂ W such that F ( U ) ⊂ U . Itis easy to show that under the assumptions of Proposition 5, the minimal invariantset E α of F α is isolated. 4 roof. We show that there exist two neighborhoods W and W of the points e and e , respectively, such that for every x ∈ W ∪ W one haslim n →∞ d ( F nα ( x ) | E α ) = 0 . (2)Since (cid:12)(cid:12)(cid:12) d f − α d x ( e ) (cid:12)(cid:12)(cid:12) < (cid:12)(cid:12)(cid:12) d f + α d x ( e ) (cid:12)(cid:12)(cid:12) <
1, there exist and L < W of e such that for every x, y ∈ W , we have | f − α ( x ) − f − α ( y ) | ≤ L | x − y | and | f + α ( x ) − f + α ( y ) | ≤ L | x − y | . Similarly, since (cid:12)(cid:12)(cid:12) d f − α d x ( e ) (cid:12)(cid:12)(cid:12) < (cid:12)(cid:12)(cid:12) d f + α d x ( e ) (cid:12)(cid:12)(cid:12) <
1, there exist and L < W of e such that for every x, y ∈ W , we have | f − α ( x ) − f − α ( y ) | ≤ L | x − y | and | f + α ( x ) − f + α ( y ) | ≤ L | x − y | . With L := max { L , L } <
1, it follows that for every x ∈ W i , i ∈ { , } , we have | f − α ( x ) − f − α ( e i ) | ≤ L | x − e i | = L dist( x, E α ) , (3)and | f + α ( x ) − f + α ( e i ) | ≤ L | x − e i | = L dist( x, E α ) . (4)Define the set W := W ∪ W ∪ E α and note that F ( W ) ⊂ W . It follows that d ( F nα ( x ) | E α ) ≤ L n d ( W | E α ) for all x ∈ W and n ∈ N , where in addition to (3) and (4), we also use the invariance of E α . This establishes(2). The conclusion follows using [LRR15, Proposition 4.1, Theorem 6.1], whichmakes full use of (H1)–(H5). Any open neighborhood B δ ( E α ) with δ > B δ ( E α ) ⊂ W separates the minimal invariant set from its dual according to[LRR15, Proposition 4.1]. By [LRR15, Theorem 6.1], the necessary condition fora discontinuous topological bifurcation is not satisfied, which means that of F α does not topologically bifurcate at the minimal invariant set E α . We note thatthe results in [LRR15] are formulated for compact spaces, but the results can beapplied if our family of set-valued mappings is truncated outside a compact setthat contains the minimal invariant set E α in its interior.It is possible to formulate a necessary condition for having a bifurcation that isa direct consequence of Proposition 5 and is expressed in terms of the derivativesof the extremal maps. Corollary 6 (Necessary condition for a discontinuous topological bifurcation) . Let { F α } α ∈ A be a family of set-valued mappings satisfying (H1)–(H5). For afixed α ∈ A , let E α = [ e , e ] be a minimal invariant set of F α . If { F α } α ∈ A topologically bifurcates at E α , then at least one of the four conditions (cid:12)(cid:12)(cid:12)(cid:12) d f − α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) ≥ , (cid:12)(cid:12)(cid:12)(cid:12) d f + α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) ≥ , (cid:12)(cid:12)(cid:12)(cid:12) d f − α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) ≥ , (cid:12)(cid:12)(cid:12)(cid:12) d f + α d x ( e ) (cid:12)(cid:12)(cid:12)(cid:12) ≥ Example 7.
Consider the family of set-valued mappings F α : R → K ( R ), α > F α ( x ) := B σ (cid:0) α arctan( x ) + x (cid:1) for all x ∈ R , where σ >
0. For a given and fixed σ >
0, there exists α > • For α ≥ α , there exist exactly two minimal invariant sets E α ⊂ ( −∞ , E α ⊂ (0 , ∞ ). They are symmetric with respect to the origin, since F α is symmetric as well. For α > α , each extremal map f + α and f − α of the set-valued mapping F α has a pair of attracting fixed points, which are boundarypoints of the minimal invariant sets E α and E α . • For α < α , there exists only one minimal invariant set E α , that is nota continuation of the two minimal invariant sets from above, since it alsocontains the interval between these two sets. Each extremal mapping f + α and f − α now has only one fixed point that bounds the new minimal invariantset (see Figure 1). αx Figure 1: Dependence of the minimal invariant sets on α for the family of set-valued mappings F α , α >
0, from Example 7: there are two minimal invariantsets for α > α , and after the discontinuous bifurcation, there is only one biggerminimal invariant set for α ≤ α .Importantly, the value α depends on σ and can be observed as the value suchthat one of the extremal mappings f + α or f − α is tangent to the identity line (see6igure 2), meaning the derivative in this intersection point is equal to 1. We willuse this later to develop an early-warning signal for this bifurcation. x − − x − − x − − f + α ( x ) and f − α ( x ) from Example 7. Fromleft to right: α > α (two minimal invariants sets), α = α (two minimal invariantsets, α < α (one minimal invariant set).Note that this bifurcation can not be detected in the limit α ր α by consid-ering derivatives of the extremal maps at the boundary of E α , α < α . This doesnot contradict the statement of Proposition 5 and Corollary 6, since at α = α ,there are two minimal invariant sets, and the necessary condition for a bifurcationestablished in Corollary 6 holds for both of these invariant sets.In the following proposition we provide a sufficient condition for a discontin-uous topological bifurcation. We require that the extremal mappings are strictlymonotone increasing. Proposition 8 (Sufficient condition for a discontinuous topological bifurcation) . Consider open intervals
I, U ⊂ R and family of set-valued mappings F α : I →K ( I ) , α ∈ U . We assume that the extreme mapping f + α as defined in (1) istwo times continuously differentiable. In addition, we require that α f + α ( x ) isdifferentiable for all x ∈ R . Assume that for a fixed α ∈ U , there exists a minimalinvariant set E α = [ e ( α ) , e ( α )] . Furthermore, we suppose that the followingfour conditions hold.(i) f + α : I → I is strictly monotone increasing for every α ∈ U .(ii) dd x f + α ( x ) (cid:12)(cid:12) x = e ( α ) = 1 .(iii) d d x f + α ( x ) (cid:12)(cid:12) x = e ( α ) > .(iv) dd α f + α ( e ( α )) (cid:12)(cid:12) α = α > .Then the family of set-valued mappings ( F α ) α ∈ U admits a discontinuous topologicalbifurcations at the minimal invariant set E α . Remark 9.
Analogous conditions can be formulated for e ( α ) ∈ E α andthe extremal mappings f − α , α ∈ U . In this case we would then require that d d x f − α ( x ) (cid:12)(cid:12) x = e ( α ) < dd α f − α ( e ( α )) (cid:12)(cid:12) α = α < roof of Proposition 8. Assume that the family E α does not admit a discontinuoustopological bifurcation at α . Then there exists a neighborhood W ⊂ U of α anda continuous family of minimal invariant sets E α = [ e ( α ) , e ( α )] for F α , α ∈ W .Because of (i) and the minimal invariance of E α = [ e ( α ) , e ( α )] for F α , it followsthat f + α ( e ( α )) = e ( α ) for all α ∈ W . (5)Conditions (ii), (iii), (iv) imply that the mappings f + α , α ∈ W , admit a saddlenode bifurcation at ( α , e ( α )), which means in particular, that there exists an¯ α > α and a neighborhood V ⊂ I of x such that the mapping f + α has no fixedpoints in V for any α ∈ ( α , ¯ α ). This contradicts (5) and finishes the proof of thisproposition.We now consider a topological bifurcation for a non-invertible mapping. Example 10.
Consider the quadratic map f α : R → R , f α ( x ) = − x + α . Weconsider the corresponding set-valued mapping F α : R → K ( R ), given by F α ( x ) := B σ ( − x + α ) , for a fixed noise level σ >
0. In the following, we determine the values of σ > α > α > E α = [ a, b ]of F α ( x ) such that f ± α are strictly monotonically decreasing on E α . Note that f + α ( a ) = b and f − α ( b ) = a , since f + α is order-reversing. We get a = f + α ( f − α ( a )) and b = f − α ( f + α ( b )) . Consider now an α > E α = E α ∪ E α = [ a, c ] ∪ [ d, b ] such that f ± α are strictly monotonically decreasing on E α .Since F α ( E α ) = E α and F α ( E α ) = E α , we have f + α ( a ) = b , f − α ( b ) = a , f − α ( c ) = d and f + α ( d ) = c . It follows that a = f − α ( f + α ( a )) , b = f + α ( f − α ( b )) , c = f + α ( f − α ( c )) , d = f − α ( f + α ( d )) . It follows in both cases that the boundary points of E are fixed points of f + α ◦ f − α or f − α ◦ f + α , so these mappings play a crucial role in the bifurcation analysis. Wedescribe the bifurcation occurring for F α as illustrated in Figure 3 for a minimalinvariant set located in [0 , ∞ ). For a choice of a fixed and small enough σ > α > B δ ( α ), δ >
0, of α , both depending on σ , such that the following statements hold.8 .775 0.800 0.825 0.850 0.875 0.900 0.9250.00.20.40.60.81.0 αx Figure 3: Dependence of the minimal invariant sets on α for the family of set-valuedmappings F α from Example 10 with σ = 0 . α < α , and after the discontinuous bifurcation, there exists adisconnected minimal invariant given as union of two intervals for α ≥ α . x .
40 0 . . . x .
40 0 . . . x .
40 0 . . . f + α and f + α (in black), f − α ◦ f + α and f + α ◦ f − α (ingray) from Example 10 with σ = 0 . Left: α < α ; one connected minimalinvariant set. Middle: α = α ; one minimal invariant set with two connectedcomponents; the second iteration mappings are tangent to the identity. Right: α > α ; one minimal invariant set with two connected components.(i) For α ∈ ( α − δ, α ), there exists a single connected minimal invariant set E α . The map f − α ◦ f + α (bottom map in gray in Figure 4) has one (attractive)fixed points q − ( α ) and the map f + α ◦ f − α (top map in gray in Figure 4) hasalso one (attractive) fixed point q + ( α ), with q − ( α ) < q + ( α ). Due to theabove analysis, we obtain E α = [ q − ( α ) , q + ( α )] for all α ∈ ( α − δ, α ) . Note that the fixed points q − ( α ) and q + ( α ) can be continued on the interval( α − δ, α + δ ).(ii) At α = α , both maps f − α ◦ f + α and f + α ◦ f − α undergo a saddle node bifurcation,and for each mapping, a new pair of fixed points is created: r − ( α ) and r − ( α )9or the mapping f − α ◦ f + α (note that r − ( α ) = r − ( α )), and r +1 ( α ) and r +2 ( α )for the mapping f + α ◦ f − α (note that r +1 ( α ) = r +2 ( α )). Both pairs are definedfor α ∈ [ α , α + δ ).For α ∈ [ α , α + δ ), the four fixed points of the second iterates q − ( α )
0, andthe map f − α ◦ f + α has two fixed points 0 < q − ( α ) < r − ( α ) = r − ( α ). Thendefine the set E α := [ r − ( α ) , q + ( α )]. Since the extremal mappings f + α and f − α are strictly monotonically decreasing on [ q − ( α ) , q + ( α )], we have E α := F α ( E α ) = [ f − α ( q + ( α )) , f + α ( r − ( α ))] . (7)It is easy to see that E α ∩ E α = ∅ . We iterate again and use the fact that q + ( α )and r − ( α ) are fixed points of f + α ◦ f − α and f − α ◦ f + α , respectively. We obtain F α ( E α ) = [ f − α ( f + α ( r − ( α ))) , f + α ( f − α ( q + ( α )))] = [ r − ( α ) , q + ( α )] = E α . (8)The sets E α and E α are cyclically mapped into each other by F α and by com-bining (7) and (8) we have[ f − α ( q + ( α )) , f + α ( r − ( α ))] = E α == F α ( E α ) = [ f − α ( f + α ( f − α ( q + ( α )))) , f + α ( f − α ( f + α ( r − ( α ))))] . Then one has f + α ( f − α ( f + α ( r − ( α )))) = f + α ( r − ( α )) , i.e. f + α ( r − ( α )) is a fixed point of f + α ◦ f − α . It is straightforward to calculatethat f + α ( r − ( α )) > q − ( α ) and f + α ( r − ( α )) < r − ( α ). This contradicts theassumption that f + α ◦ f − α has only one fixed point for x >
0. It follows that10 − α ◦ f + α and f + α ◦ f − α both have two exactly different fixed points at α = α and x >
0, one of which is tangent to the identity.We note that the extremal maps at the boundary of the sets still satisfy thenecessary condition of Corollary 6 (for any x > , x ∈ E α , one has | ( f ± α ) ′ ( x ) | > f − α ◦ f + α and f + α ◦ f − α must be considered instead here, although the bifurcationdiagram is very similar to the the pitchfork bifurcation described in Example 7. Set-valued mappings describe the topological part of a discrete-time random dy-namical system with bounded noise, in the sense that the support of the ergodicstationary measures correspond to the minimal invariant sets of the set-valuedmapping [ZH07]. This means that the results of the previous section extend tothe context of random dynamical systems. We have seen in the previous sectionthat that the derivative of the extremal mappings at the boundary of the minimalinvariant sets can be used to indicate a discontinuous topological bifurcation. Inthis section, we present an algorithm for estimating the derivative of the extremalmaps from a time series generated by iterations of a random dynamical systemwith bounded noise.Consider a time series S = { x , x , x , . . . , x n } as a realization of iterations of x i +1 = h ( x i , ξ i ) for all i ∈ { , . . . , n − } , (9)where h : R × [0 , → R , and { ξ i } i ∈ N is a sequence of i.i.d. random variables withvalues in [0 , h has an associated set-valued mapping F defined by F ( x ) = h ( x, [0 , x ∈ R . We suppose in the following that the extremal maps f + and f − of F , as definedin (1), are continuously differentiable.Let E = [ e , e ] be a minimal invariant set for F such that f − ( e ) = e and f + ( e ) = e . We aim at providing an approximation of the derivative of f − at e ,depending on the choice of two intervals of radius δ and a choice of two intervalsin the image of the map h of size ε . Note that there is a analogous algorithm forestimating the derivative of f + at e . Algorithm 11.
We approximate d f − d x ( e ) as follows.(A1) Consider two subintervals of E of length δ >
0, given by I := ( m , m + δ ) and I := ( m , m + δ ) , δ := m − m − δ > I and I are separated), and define I := [ m , m + δ ] . We require that I is close e .(A2) Assume that there exist w , . . . , w k ∈ S ∩ I and z , . . . , z k ∈ S ∩ I suchthat w i = x γ i and z i = x η i for some sequences { γ i } i =1 ,...,k and { η i } i =1 ,...,k ,and each consecutive iteration w + i = x γ i +1 is contained in˜ I := h min x ∈ I f − ( x ) , min x ∈ I f − ( x ) + ε i , and z + i = x η i +1 is contained in˜ I := h min x ∈ I f − ( x ) , min x ∈ I f − ( x ) + ε i for a given and fixed ε > f − in e is then givenby D := k P k i =1 z + i − k P k i =1 w + i ∆ , where ∆ = m + δ − m is the length of the interval I .We note that the approximation D depends on the time series S , on ε and onthe choices of the intervals I and I .It is well known that the Markov semi-group generated by (9) admits an ergodicstationary measure that is supported on the minimal invariant set E (see [ZH07]).We assume that the stationary measure is absolutely continuous with respect tothe Lebesgue measure on E , and the Lebesgue density does not vanish in theinterior of E . In addition, (9) induces a Markov chain { X i } i ∈ N that is strictlystationary and ergodic.We rewrite the map h as h ( x, ξ ) = f − ( x ) + g ( x, ξ ) , where g : R × [0 , → R is a map of the same regularity as h . We introduce therandom variable X + := f − ( X ) + g ( X, ξ ) , (10)where X is distributed according to the stationary distribution of the Markovsemi-group generated by (9). Note that X + has the same distribution as X , butwe will mainly be interested in the distribution of ( X, X + ), and it can be shownthat the delayed Markov chain { ( X i , X i +1 ) } i ∈ N is stationary and ergodic withthis distribution.We rewrite the difference quotient introduced in (A3) in terms of the aboverandom variables, and we show that the quantity calculated in (A3) is an estimatefrom below of the derivative of the extremal map f − in the interval I . Note that12e require I and I , and hence I , to be close to e , as stated in (A1), if we wantto approximate d f − d x ( e ).In (A3), two time averages are performed in order to approximate the deriva-tive. We now use the ergodic theorem to justify this approach. Proposition 12.
Consider the ergodic stationary process { X i } i ∈ N , and let I , I , ˜ I and ˜ I be the intervals introduced in (A1) and (A2). Define for each i ∈ N the random variables χ ( i ) := I ( X i ) ˜ I ( X i +1 ) and χ ( i ) := I ( X i ) ˜ I ( X i +1 ) , where A denotes the indicator function, i.e. A ( x ) = 1 if x ∈ A , and A ( x ) = 0 if x / ∈ A . Define for each n ∈ N the random variables k ( n ) := n X i =1 χ ( i ) and k ( n ) := n X i =1 χ ( i ) , (11) and consider the random variables X and X + as introduced in (10) . Then thefollowing limits exist almost surely, and the two inequalities lim n →∞ k ( n ) n X i =1 χ ( i ) h ( X i , ξ i ) ≥ E (cid:2) I ( X ) ˜ I ( X + ) f − ( X ) (cid:3) P (cid:0) X ∈ I , X + ∈ ˜ I (cid:1) (12) and lim n →∞ k ( n ) n X i =1 χ ( i ) h ( X i , ξ i ) ≤ E (cid:2) I ( X ) ˜ I ( X + ) f − ( X ) (cid:3) P (cid:0) X ∈ I , X + ∈ ˜ I (cid:1) + ε (13) hold. Remark 13.
Denote the joint Lebesgue density of X and X + by φ : R × R → R .We aim at rewriting the quantity E (cid:2) I ( X ) ˜ I ( X + ) f − ( X ) (cid:3) P (cid:0) X ∈ I , X + ∈ ˜ I (cid:1) . Firstly, the numerator reads as E (cid:2) I ( X ) ˜ I ( X + ) f − ( X ) (cid:3) = Z R Z R I ( x ) ˜ I ( y ) f − ( x ) φ ( x, y ) d x d y . Note that the function x R ˜ I φ ( x, y )d y is strictly positive whenever x ∈ I , since P (cid:0) X i +1 ∈ ˜ I | X i ∈ I (cid:1) >
0. The extra term I ( X ) ˜ I ( X + ) inside the integralrestricts the domain of integration, and hence, φ is not a probability density on I × ˜ I . However, P ( X ∈ I , X ∈ ˜ I ) is the desired normalization factor, as wehave P ( X ∈ I , X ∈ ˜ I ) = Z I d x Z ˜ I φ ( x, y ) d y > . Hence, the function x R I d x R ˜ I φ ( x, y )d y Z ˜ I φ ( x, y ) d y
13s a probability density describing the statistics of those points in I for which theconsecutive iteration is in ˜ I . Proof.
We prove (13). Note thatlim n →∞ k ( n ) n X i =1 χ ( i ) h ( X i , ξ i ) = lim n →∞ k ( n ) nn n X i =1 χ ( i ) (cid:0) f − ( X i ) + g ( X i , ξ i ) (cid:1) . We consider the limitslim n →∞ n k ( n ) and lim n →∞ n n X i =1 χ ( i )( f − ( X i ) + g ( X i , ξ i )) (14)separately. For both limits the ergodic theorem is applied to the delayed Markovchain { ( X i , X i +1 ) } n ∈ N , and we consider the observables O ( X i , X i +1 ) = I ( X i ) ˜ I ( X i +1 ) ,O ( X i , X i +1 ) = I ( X i ) ˜ I ( X i +1 )( f − ( X i ) + g ( X i , ξ i )) = I ( X i ) ˜ I ( X i +1 ) X i +1 . The first observable leads tolim n →∞ n k ( n ) = lim n →∞ n n X i =1 I ( X i ) ˜ I ( X i +1 ) = P ( X ∈ I , X + ∈ ˜ I ) . The second observable yieldslim n →∞ n n X i =1 χ ( i )( f − ( X i ) + g ( X i , ξ i )) ≤ lim n →∞ n n X i =1 χ ( i )( f − ( X i ) + ε )= E [ I ( X ) ˜ I ( X + ) f − ( X )] + ε P ( X ∈ I , X + ∈ ˜ I ) . combining these two limit relations giveslim n →∞ k ( n ) nn n X i =1 χ ( i ) h ( X i , ξ i ) ≤ E [ I ( X ) ˜ I ( X + ) f − ( X )] P ( X ∈ I , X + ∈ ˜ I ) + ε . The proof for (12) is analogous, with the difference that the perturbation g ( x, ξ )is estimated from below by zero.The following corollary is a direct consequence of Proposition 12 and it showsthat D in (A3) is an estimate from below of the derivative of f − in the interval I . Corollary 14.
Under the conditions of Proposition 12, almost surely, we havelim n →∞ k ( n ) n X i =1 χ ( i ) h ( X i , ξ i ) − k ( n ) n X i =1 χ ( i ) h ( X i , ξ i ) ! ≤ max x ∈ I ( f − ) ′ ( x ) + ε ∆ . roof. Proposition 12 yieldslim n →∞ k ( n ) n X i =1 χ ( i ) h ( X i , ξ i ) − k ( n ) n X i =1 χ ( i ) h ( X i , ξ i ) ! ≤ E [ I ( X ) ˜ I ( X + ) f − ( X )] P ( X ∈ I , X + ∈ ˜ I ) − E [ I ( X ) ˜ I ( X + ) f − ( X )] P ( X ∈ I , X + ∈ ˜ I ) ! + ε ∆ (15) ≤ max x ∈ I ,y ∈ I (cid:0) f ( y ) − f ( x ) (cid:1) + ε ∆ ≤ max x ∈ I ,y ∈ I f ( y ) − f ( x ) y − x + ε ∆ ≤ max x ∈ I ( f − ) ′ ( x ) + ε ∆ . We approximated the expectation in (15) by considering the maximum. Next weused ∆ > y − x and the intermediate value theorem. Remark 15.
The error ε ∆ depends on the size of the intervals I , I , ˜ I and ˜ I .In order to reduce this error, one could shrink the size of the intervals, but as aconsequence, we will need to wait a longer time to obtain a sample of points insatisfying the conditions in (A2). Remark 16.
In Corollary 14, the estimate of the derivative takes two typesof errors into account, one coming from the noise, and the other given by theunderestimation of the difference quotient D (note that ∆ serves as upper boundin the denominator). We want to compare the derivative of the extremal map f − ( x ) = h ( x,
0) with the difference quotient introduced in (A3) for k = k = 1.For any x ∈ I , y ∈ I and ξ ∈ [0 , (cid:12)(cid:12)(cid:12)(cid:12) max c ∈ [ m ,m + δ ] ∂h∂x ( c, − h ( y, ξ ) − h ( x, ξ )∆ (cid:12)(cid:12)(cid:12)(cid:12) == (cid:12)(cid:12)(cid:12)(cid:12) max c ∈ [ m ,m + δ ] ∂h∂x ( c, − y − xy − x h ( y, ξ ) − h ( x, ξ )∆ (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) max c ∈ [ m ,m + δ ] ∂h∂x ( c, − y − x ∆ ∂h∂x (¯ c, ξ ) (cid:12)(cid:12)(cid:12)(cid:12) == (cid:12)(cid:12)(cid:12)(cid:12) max c ∈ [ m ,m + δ ] ∂h∂x ( c, (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − y − x ∆ ∂h∂x (¯ c, ξ )max c ∈ [ m ,m + δ ] ∂h∂x ( c, ξ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . Here, we used the intermediate value theorem to obtain ¯ c ∈ [ m , m + δ ]. Notethat if the map h is at least continuously differentiable, the term on the left handside is bounded, and then the error can be controlled by the quantity (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − y − x ∆ ∂h∂x (¯ c, ξ )max c ∈ [ m ,m + δ ] ∂h∂x ( c, ξ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . The error depends on the size of the intervals, since y − x ∆ appears and on thevariation of the function h as the ratio of the derivatives appears. We want toemphasize that the shorter the intervals I , I , ˜ I , ˜ I are, the smaller is the error is.On the other hand, from the application point of view, the probability of having15onsecutive points in I i and ˜ I is also lower, since P ( X ∈ I i , X + ∈ ˜ I i ) for i ∈ { , } ,tends to zero as the size of the intervals tend to zero. When applying Algorithm 11to a time series, it is then crucial to consider the size of the intervals in order tohave a suitable sample, and the variation of the extremal function in the intervals,in order to have a meaningful estimate. Remark 17.
In the case of additive noise, Algorithm 11 can be simplified signif-icantly. In (A2), the intervals ˜ I and ˜ I , and the condition that the next iterationmust be contained in ˜ I and ˜ I can be dropped. In order to show this we state andsketch the proof of Proposition 12 and Corollary 14 in the case of additive noise.Consider the stationary process { X i } i ∈ N as previously defined, and let ξ i be asequence of i.i.d. bounded random variables distributed according to a randomvariable ξ such that X i is independent of ξ i for every i ∈ N . In the additive noisecase, the map h reads as h ( X i , ξ i ) = f ( X i ) + ξ i . Define for i ∈ N the random variablesΓ ( i ) := I ( X i ) and Γ ( i ) := I ( X i ) . Then the statements (12) and (13) of Proposition 12 read aslim n →∞ k ( n ) n X i =1 Γ ( i )( f ( X i ) + ξ i ) = E [ I ( X ) f ( X )] P ( X ∈ I ) + E [ ξ ] (16)and lim n →∞ k ( n ) n X i =1 Γ ( i )( f ( X i ) + ξ i ) = E [ I ( X ) f ( X )] P ( X ∈ I ) + E [ ξ ] , (17)where k , k are defined in (11). The calculation of the limits (16) and (17) isanalogous to the calculation performed in the proof of Proposition 12. Whencompared to (12) and (13), note that in the right hand side of the equalities (16)and (17), a different contribution of the noise appears. For (16), we havelim n →∞ k ( n ) n X i =1 Γ ( i ) ξ i = 1 P ( X ∈ I ) E [ ξ ] P ( X ∈ I ) = E [ ξ ] , and similarly for (17), we havelim n →∞ k ( n ) n X i =1 Γ ( i ) ξ i = E [ ξ ] . We do not need an estimate of the limits, since the effect of the noise will cancelout. On the other hand, the additive version of Corollary 14 reads aslim n →∞ k ( n ) n X i =1 Γ ( i )( f ( X i ) + ξ i ) − k ( n ) n X i =1 Γ ( i )( f ( X i ) + ξ i ) ! ≤ max x ∈ I ( f − ) ′ ( x ) .
16n fact, one haslim n →∞ k ( n ) n X i =1 Γ ( i )( f ( X i ) + ξ i ) − k ( n ) n X i =1 Γ ( i )( f ( X i ) + ξ i ) ! ≤ (cid:20) E [ I ( X ) h ( X )] P ( X ∈ I ) + E [ ξ ] − E [ I ( X ) h ( X )] P ( X ∈ I ) − E [ ξ ] (cid:21) (18) ≤ max x ∈ I ,y ∈ I (cid:0) f ( y ) − f ( x ) (cid:1) ≤ max x ∈ I ,y ∈ I f ( y ) − f ( x ) y − x ≤ max x ∈ I ( f − ) ′ ( x ) . The proof is the same as in Corollary 14, except for step (18), where the contri-bution of the noise cancels out. Hence we have that the intervals ˜ I and ˜ I andthe condition on the next iterated in (A2) can be dropped.The absence of the intervals ˜ I and ˜ I implies that the algorithm is moreeffective in the additive case for two reasons. Firstly, there are less conditions onthe points of the time series, implying that one has generally a greater number ofpoints eligible for the calculating the estimate when compared to the general case.Secondly, we obtain a better estimate from below of the derivative, since the term ε ∆ does not appear in the additive case. In this section, we numerically study the stochastic return map of the Kopermodel with additive bounded noise. The Koper model is a prototypical exampleof a system exhibiting mixed-mode oscillations (MMOs), i.e. periodic orbits withwidely varying amplitudes in one period. Depending upon the value of a parameter λ ∈ R the system exhibits different stable modes of oscillations. Bifurcations ofthe minimal invariant set of the set-valued mapping associated to the stochasticreturn map correspond to the system shifting from an oscillation mode to another.With the algorithm of Section 3 we reconstruct the derivative of the extremal map f − λ for different values of λ to find an early-warning sign for upcoming bifurcations.We consider the Koper model given by the ordinary differential equation ε d x d τ = y − x + 3 x =: f ( x, y, z ) , d y d τ = kx − y + λ ) + z, (19)d z d τ = λ + y − z, where k, λ ∈ R are real parameters and ε is assumed to be small and positive.We fix k = −
10 and study the system for a special parameter range λ ∈ I =( − , − / + ε with one fast variable x , and two slow variables y, z . Therefore,17he analysis of MMOs can be carried out using the framework of multiple timescale dynamical systems [DGK +
12, Kue15]. We provide a brief description of themultiscale geometry of the system. The critical manifold is given by C := { ( x, y, z ) ∈ R | y = x − x } . Note that C can be viewed as the algebraic constraint for the slow subsystemobtained from (19) in the limit ε →
0. Furthermore, we can also re-scale time t := τ /ε and on the time scale t , the critical manifold consists of equilibrium pointfor the fast subsystem d x d t = y − x + 3 x, d y d t = 0 , (20)d z d t = 0 , obtained again by a formal singular limit ε → C is normally hyperbolic (whichjust means here that ( ∂ x f ) | C is nonzero) away from the two fold lines F + := { x = 1 } ∩ C and F − := { x = − } ∩ C . The folds F ± separate the normally hyperbolic repelling manifold C r = C ∩ { x < } ∩ { x > − } from the two normally hyperbolic attracting manifolds C a +0 = C ∩ { x > } C a − = C ∩ { x < − } . In particular, C r consists of repellingequilibria of (20) while C a ± consist of attracting equilibria. There are two foldedsingularities for the system located at q ± = ( ± , ± , λ ± , which can be viewed as equilibria of a suitable desingularized slow subsystem; seee.g. [DGK +
12, Kue15, Wec05] for more details. The main point for us here is thatit is well-understood how the interplay between folded singularities and a globalfast-slow return mechanism can generate MMOs as shown in Figure 5. Due to thesystem symmetry ( x, y, z, λ, k ) → ( − x, − y, − z, − λ, k ) , we can consider q + only. It can be shown that for λ ∈ I = ( − , − /
6) the foldedsingularity q + is a so-called folded node [Kue11]. It is known [BKW06, SW01]that in the vicinity of a folded node various small-amplitude oscillations (SAOs)can be generated. The idea is that a set of special canard orbits partitions theneighborhood of q + in the ( x, z ) plane into certain sectors of rotation. An orbitof system (19) entering a sector will perform s small oscillations before making atransition to a fast segment. To each sector corresponds a different number s ofoscillations.Due to the geometry of the S -shaped critical manifold and its attractivity prop-erties it follows that we also expect large-amplitude oscillations (LAOs) in the sys-tem due to a relaxation-oscillation type switching between O ( ε )-neighbourhoods18f the two attracting branches C a ± (see also Figure 5). The combination of theSAO and the LAO mechanisms can then yield MMOs.It is possible to define a return map P λ on a suitable codimension-one sectionΣ to study the MMOs via the return map P λ : Σ → Σ . Here it will be convenient to fix Σ asΣ = (cid:26) ( x, y, z ) ∈ R | x = − , y ∈ ( − , − , z ∈ ( z , z ) (cid:27) . The points z , z will be chosen according to λ and the noise amplitude. Due tothe fast-slow structure, the map P λ is almost one-dimensional [Guc08] and it iscommon to project the map P λ further to reduce it to a one-dimensional map p λ : Σ ∩ { y = y } → R , where we define the map as z p λ ( z ) := ( π z ◦ P λ )( y , z )with π z denoting the projection onto the z axis. The projection is made in orderto simplify the analysis even if the resulting p λ is non-invertible. Sometimes it iseven possible to calculate an explicit (asymptotic) expression for the nearly one-dimensional map in certain parameter regimes as explained for a slightly different,but related, three-time scale model in [KPK08]. −2 −1 0 1 2 3−8.2−8−7.8−7.6−7.4−7.2−7−2.5−2−1.5−1−0.500.511.522.5 xzy ty Figure 5: A trajectory of a solution of system (19) and its projection on the y component, corresponding to the initial condition ( x, y, z ) = ( − . , − , −
8) and λ = − . z p λ ( z ). The jumpscorrespond to the location of canard orbits; see Figure 6. An attracting fixed19 zp λ Figure 6: Deterministic return map for the Koper system λ = − .
9. The jumpscorrespond to the location of the canard orbits. In fact, the map is continuousand the gaps represent steep jumps [Guc08, DGK + p λ corresponds to a stable mode of oscillations. It is keythat by varying λ the fixed point of the return map can change and the systemcan transition to a different stable mode of oscillation. The goal is to introduce(bounded) noise in the system of ordinary differential equations to predict the shiftto a new mode of oscillations when λ is slowly increased towards a bifurcation.Consider again the system (19).We introduce bounded additive noise in the model and consider the family ofrandom maps x ( t ) y ( t ) z ( t ) = x ( s ) y ( s ) z ( s ) + y ( s ) − x ( s ) + 3 x ( s ) ε (cid:0) kx ( s ) − y ( s ) + λ ) + z ( s ) (cid:1) ε (cid:0) λ + y ( s ) − z ( s ) (cid:1) ( t − s ) + A ( ξ ( t ) − ξ ( s )) , (21)where ξ = ( ξ ( t ) , ξ ( t ) , ξ ( t )) ⊤ ∈ [ − , is a three-dimensional bounded stochasticprocess with independent increments, t > s , and we fix the time step ( t − s ). Eachincrement ξ ( t ) − ξ ( s ) is uniformly distributed in [ − , A = σ . . . . . . . . . . Note that in the limit ( t − s ) → stochastic return map , p stλ defined on the samesection Σ of the map p λ . The definition of a stochastic return map is more prob-lematic, since the noise can destroy the return mechanism or we could have anearly return to the section, i.e. before the trajectory goes through the global re-turn mechanism. A formal construction is presented in [WK90]. In the numericalreconstruction performed in this paper Σ and the initial conditions are chosen insuch a way that a realization of the orbit of the random differential equation doesmake a global return to Σ. In practice we consider initial conditions ( x, y, z ) with x < − , y < z ∈ ( z , z ) and we stop the iteration once the orbit ( x ( t ) , y ( t ) , z ( t ))is such that y ( t ) < x ( t ) crosses x = − . -6.89 -6.885 -6.88 -6.875 -6.87 -6.865 -6.86 -6.8550.20.40.60.811.21.41.61.8 λDp λ -6.89 -6.885 -6.88 -6.875 -6.87 -6.865 -6.860.20.30.40.50.60.70.80.911.1 λDf − Figure 7: (a) Derivative of the deterministic return map p λ evaluated at the fixedpoint x λ for λ ≤ λ = − .
86. For λ > − .
86 the derivative is evaluated in x λ asthe fixed point disappears and the map has a periodic orbit. (b) Derivative of f − atthe lower boundary of the minimal invariant set. For any given λ the lower boundof the minimal invariant set is retrieved from Figure 8. We considered n = 3000realizations of the stochastic return map, say z i = p stλ ( ξ i , m λ ), w i = p stλ ( ξ i m λ + ε ),for i = 1 , . . . , n and approximated the derivative by d λ = min i { w i }− min i { z i } ε . Weconsidered ε = 0 . λ and a fixed noise amplitude σ = 0 .
01. Notice that at theboundary between two sectors a trajectory of the random differential equationis more likely to have two possible oscillations modes according to the noise re-alization. This is a consequence of the fact that also the stochastic map is aprojection onto the z axis and the noise affecting the other components may alsocause a jump to a different sector. The graph of the stochastic return map ap-proximates the graph of the set-valued map. From the reconstruction we find thatthere exists an isolated minimal invariant set E λ for λ ∗ < λ < λ for some λ ∗ .In λ ≈ − .
86, we observe that the set E λ satisfies the necessary conditions forhaving a bifurcation. Since it does not disappear for λ > λ the set E λ blows uplower semi-continuously in the sense presented in [LRR15]. The bifurcation canbe observed in the reconstruction of the minimal invariant set in Figure 8. Weevaluated the derivative of the extremal map f − associated to the stochastic re-turn map p stλ with the algorithm presented in Section 3. The methodology and theresults are presented in Figure 7. The derivative of f − λ evaluated at the boundary21 λz Figure 8: Minimal invariant set of the stochastic return map reconstructed for λ ∈ [ − . , − . λ ≈ − . E λ crosses 1 from below as λ crosses λ . This can be used as early-warningsignal due to Proposition 8. −8.2 −8.1 −8 −7.9 −7.8 −7.7−8.25−8.2−8.15−8.1−8.05−8−7.95−7.9−7.85−7.8−7.75 (a) zp λ −8.2 −8.1 −8 −7.9 −7.8 −7.7−8.2−8.15−8.1−8.05−8−7.95−7.9−7.85−7.8−7.75−7.7 (b) zp λ Figure 9: (a) λ = − .
9. Approximation of the set-valued dynamical system as-sociated to the stochastic return map before bifurcation. The minimal invariantset is easily identified by the intersection of the identity line and the map. (b) λ = − .
85 Approximation of the set-valued dynamical system associated to thestochastic return map after bifurcation. Notice that the system now can randomlyjump between two different modes of oscillations.22 eferences [APW17] P. Ashwin, C. Perryman, and S. Wieczorek,
Parameter shifts fornonautonomous systems in low dimension: bifurcation-and rate-induced tipping , Nonlinearity (2017), no. 6, 2185–2210.[Arn98] Ludwig Arnold, Random Dynamical Systems , Springer Monographs inMathematics, Springer-Verlag, Berlin, 1998.[BKW06] M. Brøns, M. Krupa, and M. Wechselberger,
Mixed mode oscillationsdue to the generalized canard phenomenon , Fields Institute Communi-cations (2006), 39–63.[CH96] S.-N. Chow and J.K. Hale, Methods of Bifurcation Theory ,Grundlehren der mathematischen Wissenschaften, vol. 251, Springer,1996.[DGK +
12] M. Desroches, J. Guckenheimer, C. Kuehn, B. Krauskopf, H. Osinga,and M. Wechselberger,
Mixed-mode oscillations with multiple timescales , SIAM Review (2012), no. 2, 211–288.[Guc08] J. Guckenheimer, Return maps of folded nodes and folded saddle-nodes ,Chaos (2008), 015108.[Kop95] M.T.M. Koper, Bifurcations of mixed-mode oscillations in a three-variable autonomous Van der Pol-Duffing model with a cross-shapedphase diagram , Physica D (1995), 72–94.[KPK08] M. Krupa, N. Popovic, and N. Kopell, Mixed-mode oscillations in threetime-scale systems: A prototypical example , SIAM Journal of AppliedDynamical Systems (2008), no. 2, 361–420.[KR11] P.E. Kloeden and M. Rasmussen, Nonautonomous Dynamical Systems ,Mathematical Surveys and Monographs, vol. 176, American Mathe-matical Society, Providence, Rhodes Island, 2011.[Kue11] C. Kuehn,
On decomposing mixed-mode oscillations and their returnmaps , Chaos (2011), no. 3, 033107.[Kue13] , A mathematical framework for critical transitions: normalforms, variance and applications , Journal of Nonlinear Science (2013), no. 3, 457–510.[Kue15] , Multiple Time Scale Dynamics , Applied Mathematical Sci-ences, vol. 191, Springer, Cham, 2015.[Kuz95] Y.A. Kuznetsow,
Elements of Applied Bifurcation Theory , AppliedMathematical Sciences, vol. 112, Springer, 1995.[LRR15] J.S.W. Lamb, M. Rasmussen, and C.S. Rodrigues,
Topological bifur-cations of minimal invariant sets for set-valued dynamical systems ,Proceedings of the American Mathematical Society (2015), no. 9,3927–3937. 23RS16] P. Ritchie and J. Sieber,
Early-warning indicators for rate-induced tip-ping , Chaos (2016), no. 9, 093116, 13.[SBB +
09] M. Scheffer, J. Bascompte, W.A. Brock, V. Brovkhin, S.R. Carpenter,V. Dakos, H. Held, E.H. van Nes, M. Rietkerk, and G. Sugihara,
Early-warning signals for critical transitions , Nature (2009), 53–59.[SW01] P. Szmolyan and M. Wechselberger,
Canards in R , Journal of Differ-ential Equations (2001), 419–453.[Wec05] M. Wechselberger, Existence and bifurcation of canards in R in thecase of a folded node , SIAM Journal of Applied Dynamical Systems (2005), no. 1, 101–139.[WK90] J.B. Weiss and E. Knobloch, A stochastic return map for stochasticdifferential equations , Journal of Statistical Physics (1990), no. 5,863–883.[ZH07] H. Zmarrou and A.J. Homburg, Bifurcations of stationary densities ofrandom diffeomorphisms , Ergodic Theory and Dynamical Systems27