From random Poincaré maps to
stochastic mixedmodeoscillation
patterns
Revised, November 14, 2014
Abstract
We quantify the effect of Gaussian white noise on fast–slow dynamical systems with one fast and two slow variables, which display mixedmode oscillations owing to the presence of a foldednode singularity. The stochastic system can be described by a continuousspace, discretetime Markov chain, recording the returns of sample paths to a Poincaré section. We provide estimates on the kernel of this Markov chain, depending on the system parameters and the noise intensity. These results yield predictions on the observed random mixedmode oscillation patterns. Our analysis shows that there is an intricate interplay between the number of smallamplitude oscillations and the global return mechanism. In combination with a local saturation phenomenon near the folded node, this interplay can modify the number of smallamplitude oscillations after a largeamplitude oscillation. Finally, sufficient conditions are derived which determine when the noise increases the number of smallamplitude oscillations and when it decreases this number.
Mathematical Subject Classification. 37H20, 34E17 (primary), 60H10 (secondary)
Keywords and phrases. Singular perturbation, fast–slow system, dynamic bifurcation, folded node, canard, mixedmode oscillation, return map, random dynamical system, firstexit time, concentration of sample paths, Markov chain.
1 Introduction
Oscillation patterns with large variations in amplitude occur frequently in dynamical systems, differential equations and their applications. A class of particular interest are mixedmode oscillations (MMOs) which are patterns consisting of alternating structures of smallamplitude and largeamplitude oscillations. Typical applications arise from chemical systems such as the Belousov–Zhabotinskii reaction [47], the peroxidase–oxidase reaction [25] and autocatalytic reactions [71] as well as from neuroscience, e.g. stellate cells [28], Hodgkin–Huxleytype neurons [73] and pituitary cells [80]. A remarkable number of models for these phenomena lead to differential equations with multiple timescales see e.g. [66, 27, 84, 40]. Frequently, it suffices to consider two timescales and study fast–slow ordinary differential equations (ODEs) which already provide many generic mechanisms leading to MMOs. For a detailed review of the topic we refer to the survey [26], the special issue [23], and references therein. The basic idea is that a local mechanism induces the smallamplitude oscillations (SAOs) while a global return mechanism leads to largeamplitude oscillations (LAOs). In this introduction, we shall just outline the main ideas; the precise development of our setup and results starts in Section 2.
For a deterministic trajectory, we can symbolically write an MMO as a sequence
(1.1) 
where denotes LAOs followed by SAOs, etc. For example, a periodic solution alternating between SAOs and LAO would be or simply with the periodicity understood. A prototypical mechanism to generate SAOs are foldednode singularities [24] which are generic in threedimensional ODEs with one fast and two slow variables [77, 7]. For the global return mechanism, one frequently encounters a relaxationtype structure induced by a cubic (or Sshaped) fastvariable nullcline, also called the critical manifold, which was studied extensively already in the context of van der Poltype oscillators; see e.g. [51, 35, 20, 24] and the references therein. Nondegenerate folds, foldednode singularities and Sshaped critical manifolds form the basic deterministic building blocks for the work in this paper. However, let us mention already here that the stochastic techniques we develop in this paper could potentially be adapted to other cases such as singular Hopf bifurcation and folded saddlenodes [38, 58], bursting oscillations [48, 30], tourbillon structures [81, 26], and other global return mechanisms [40, 61]. Although it is certainly of high interest to study all these cases, it seems to us that the combination of folded singularities and relaxation oscillations is a natural first step as both components are basic elements which occur in a large variety of different models [26].
While in some experiments, remarkably clear MMO patterns have been observed [42, 47], in many other cases the SAOs in the patterns appear noisy [28, 65]. Weak noise acting on a dynamical system is known to induce a variety of phenomena, ranging from small fluctuations around deterministic solutions to large excursions in phase space, as shown, e.g., in stochastic resonance [12, 34, 64] and transitions near tipping points [74, 60, 4]. In the context of oscillatory patterns, the effect of noise on MMO patterns in lowdimensional prototypical models has been studied, for instance, in [52, 76, 69, 87, 43, 62], using numerical simulations, bifurcation theory and asymptotic descriptions of the Fokker–Planck equation.
This work concerns the effect of noise on fast–slow differential equations with one fast and two slow variables, containing a foldednode singularity and an Sshaped critical manifold responsible for the global return mechanism. The resulting stochastic differential equations (SDEs) show a subtle interplay between noise, local and global dynamics, which requires a careful analysis of the behaviour of stochastic sample paths. Our approach builds upon our earlier work [17] which in turn was based upon a pathwise approach to fast–slow SDEs [14, 15].
Our main focus is the derivation of estimates for the Poincaré (or return) map of the stochastic system, for a conveniently chosen twodimensional section . Deterministic return maps in the presence of foldednode singularities have been analyzed, e.g., in [37, 54]. Although the twodimensional Poincaré map is invertible, the strong contraction near attracting critical manifolds implies that it is close to a onedimensional, usually noninvertible map. Figure 1 (a) shows an example of such a onedimensional deterministic return map . The apparent discontinuities are in fact points where the map’s graph displays very narrow dips, due to the presence of socalled canard orbits. Canards are particular solutions of the system staying close to both the attracting and repelling parts of the critical manifold [10, 7, 11], which separate the phase space into sectors of rotation characterized by different numbers of SAOs [24].
The concept of return maps has been extended to stochastic systems, see for instance [85, 18, 44]. This requires some care, because the rapid fluctuations of stochastic sample paths prevent one from using their very first return to to define the map. Instead, one has to consider the first return after a suitably defined, sufficiently large excursion in phase space has taken place. With these precautions, successive intersections of sample paths with define a continuousspace, discretetime Markov chain. The distribution of is obtained from the distribution of via integration with respect to a transition kernel . Under suitable regularity assumptions [6], the theory of harmonic measures ensures that the kernel admits a smooth density , so that the evolution of is specified by an integral equation, namely
(1.2) 
holds for all Borel sets ; see for instance [16, Sections 5.2 and 5.3]. The main aim of the present work is to provide estimates on the kernel . Part of the mathematical challenge is due to the fact that the deterministic flow is not a gradient flow, and thus the stochastic system is irreversible.
Figure 1 (b)–(d) shows simulated stochastic return maps for increasing noise intensity. For each value of , the red points indicate the value of for different realizations of the noise. The deterministic return map is plotted in blue for comparison. Several interesting phenomena can be observed:

The size of fluctuations increases with the noise intensity;

Orbits in sectors with a small number of SAOs (inner sectors) are less affected by noise than those in sectors with a large number of SAOs (outer sectors);

There is a saturation effect, in the sense that for large enough SAO numbers, the typical value of the stochastic return map and its spreading become independent of the sector;

The saturation effect sets in earlier for larger noise intensities.
While the first phenomenon is not surprising, the other observed features are remarkable, and can lead to nonintuitive effects. In the example shown in Figure 1, the deterministic map has a stable fixed point in the 11th sector, so that the deterministic system will display a stable MMO pattern . For sufficiently strong noise, the stochastic system will asymptotically operate in the 12th sector, with occasional transitions to neighbouring sectors such as sectors 11 and 13. Hence, the noise shifts the global return to a higher rotation sector. However, two more noiseinduced effects may also affect the number of observed SAOs. First, the noise may alter the number of SAOs for orbits starting in a given sector, by causing earlier escapes away from the critical manifold. Second, it may produce fluctuations large enough to mask small oscillations. All these effects must be quantified and compared to determine which oscillatory pattern we expect to observe.
The estimates on the kernel we provide in this work yield quantitative information on the above phenomena. In particular, we obtain estimates on the typical size of fluctuations as a function of noise intensity and sector number, and on the onset of the saturation phenomenon. These results complement those already obtained in [17] on the size of fluctuations near a foldednode singularity.
The structure of this article is the following: After introducing the deterministic setup in Section 2, we provide first estimates on noiseinduced fluctuations in Section 3. Sections 4 and 5 extend the analysis to a neighbourhood of the regular fold and of the folded node, respectively. Section 6 combines all the local estimates to provide quantitative results on the kernel. The main results are:
A short discussion of the consequences of these results on the observed MMO patterns is given in Section 6.3. Finally, Section 7 illustrates these results with numerical simulations for the Koper model.
The results obtained here are a first step towards the understanding of stochastic MMOs, that calls for further work. In particular, it would be desirable to obtain a more precise description of the possible MMO patterns. Let us mention two possible ways to achieve this:

Singularly perturbed Markov chains: Consider the ideal case where the sectors of rotation form a Markov partition, meaning that the image of each sector is entirely contained in a sector. Then the dynamics can be described by a topological Markov chain between sectors, see Figure 2. Such a chain will in general not be irreducible, for instance, for the chain shown in Figure 2, state is an absorbing state. In the presence of noise, however, new transitions between sectors appear, typically yielding an irreducible chain. In this sense, the chain for the stochastic system is a singular perturbation of its deterministic limit. For weak, nonvanishing noise, transitions between all states may become possible, but transition times diverge as the noise intensity goes to zero. Methods allowing to determine transition rates in singularly perturbed Markov chains for small positive noise have been developed, for instance, in [75, 41, 5, 86].

Metastable transitions between periodic orbits: Consider a situation where the deterministic system admits several stable periodic orbits, each corresponding to an MMO pattern. Weak noise will induce rare transitions between these orbits. The theory of large deviations [33] provides a way to estimate the exponential asymptotics of transition rates (Arrhenius’ law [3]), via a variational problem. In the reversible case, Kramers’ law [31, 53] provides a more precise expression for transition rates, which are related to exponentially small eigenvalues of the diffusion’s infinitesimal generator, see for instance [21, 22], and [13] for a recent review. For irreversible systems, such precise expressions for transition rates are not available. However, the spectraltheoretic approach may still yield useful information, as in similar irreversible problems involving random Poincaré maps [18, 16].
Notations: We write to denote the absolute value and for the Euclidean norm. For we write for the smallest integer not less than . Furthermore, for we use and . Regarding asymptotics, we use in the usual way, i.e., we write as if and only if . The shorthand is used whenever and hold simultaneously. Furthermore, by we indicate that . Vectors are assumed to be column vectors and denotes the transpose of a vector .
Acknowledgements: C.K. would like to thank the Austrian Academy of Sciences (ÖAW) for support via an APART fellowship as well as the European Commission (EC/REA) for support by a MarieCurie International Reintegration Grant. B.G. and C.K. thank the MAPMO at the Université d’Orléans, N.B. and C.K. thank the CRC 701 at the University of Bielefeld for kind hospitality and financial support. Last but not least, we would like to thank an anonymous referee for constructive comments that led to improvements in the presentation.
2 Mixedmode oscillations – The setup
In this section, we shall outline a typical setup for deterministic mixedmode oscillations based upon threedimensional fast–slow systems of the form
(2.1)  
where and is a small parameter. Throughout, we shall make the following assumption:

The functions are of class .
In particular, (A0) implies that on a fixed compact set there exist uniform bounds on . We remark that the system (2.1) is allowed to depend smoothly upon further system parameters although we do not indicate this dependence in the notation. The critical set of (2.1) is
(2.2) 
Motivated by several applications, such as the Hodgkin–Huxley model [45, 73], the Koper model [51, 59], the forced van der Pol equation [79, 36] and the Rössler model [72], we will assume that the geometric structure of the critical set is an Sshaped smooth manifold; see also Figure 3. More precisely, this assumption can be stated as follows:

Suppose is a smooth manifold composed of five smooth submanifolds,
(2.3) where the twodimensional submanifolds are normally hyperbolic attracting, while the twodimensional submanifold is normally hyperbolic repelling, i.e.,
(2.4) and are onedimensional smooth fold curves consisting of generic fold points
(2.5) Without loss of generality we assume from now on that for all .
Fenichel theory [32] shows that for , the critical submanifolds and perturb to invariant manifolds and , which are close to and in points bounded away from the fold curves .
Setting in (2.1) leads to the slow subsystem
(2.6)  
which is solved by the socalled slow flow. Differentiating implicitly with respect to yields
(2.7) 
for the slow flow, so that the slow subsystem (2.6) can be written as
(2.8) 
where it is understood that all functions are evaluated at points with . One may use that (2.8) can locally be written as a closed system by applying the implicitfunction theorem to express as a graph, e.g. , near the fold as .
Observe that (2.8) is singular on the fold curves as on . The desingularized slow subsystem is obtained by multiplying the righthand side of (2.8) by and applying a rescaling of time. It reads
(2.9) 
We make the following further assumptions:
Let us stress that the above geometric assumptions (A1)–(A3), as well as several further assumptions to follow, provide a convenient framework but that the deterministic and stochastic techniques we present here apply to a much wider range of multiscale systems displaying oscillatory patterns.
On the fast timescale the limit of (2.1) leads to the fast subsystem
(2.11)  
which is solved by the fast flow. It is helpful to decompose the singular limit flows and their perturbations into several parts; see Figure 3 for an illustration. In particular, we consider the sections of the form
(2.12)  
for , , suitably chosen to capture the return map. For an appropriate choice of the constants and (see below or consider the approach in [59]), there are welldefined maps from to .

The geometry of the flowinduced maps and sections is as shown in Figure 3.
In particular, Assumption (A4) implies that there is an transition time on the slow timescale from to as well as from to . (A4) incorporates that there is an spatial separation between each pair of fold/drop curves and it guarantees there is an transition time on the fast timescale from to as well as from to . Furthermore, we exclude the case of a singular Hopf bifurcation [38, 39], where an equilibrium of the full system (2.1) may occur in the neighbourhood of a folded node.
There are four distinct important parts of the flow to analyze:

the flow near the folded node ,

the fast segment ,

the slowflow region near , and

the nondegenerate fold via .
The map can be covered by the same techniques as , and is similar to .
The geometry of flow maps and the possible generation mechanisms for mixedmode oscillations under the assumptions (A0)–(A4) are wellknown; see for example [24, 26]. A main idea is that twisting of the slow manifolds and near a folded node generates SAOs and the global return mechanism via the Sshaped critical manifold induces the LAOs. Fixed points of a full return map, say , correspond to MMOs with a certain pattern
(2.13) 
The main question we address in this paper is how noise influences the patterns (2.13). We are going to split the analysis into two main parts. In Section 3 we provide basic estimates and consider the global part of the return map. Sections 4–5 address local dynamics in the regions near the regular fold and the folded node.
3 The stochastic system
3.1 Estimating stochastic deviations
As a stochastic extension to (2.1) we consider the fast–slow SDE
(3.1)  
where is a dimensional standard Brownian motion on a probability space . The maps
(3.2) 
may depend on , and are assumed to be and to satisfy the usual boundedgrowth condition guaranteing existence of a unique strong solution of (3.1). We shall adopt the shorthand notation to write just instead of .
We will assume that the diffusion coefficients satisfy the following uniform ellipticity assumption:

Let
(3.3) be the diffusion matrix. There exist constants such that
(3.4)
Remark 3.1.
In fact, most of our results remain valid under a weaker hypoellipticity assumption (cf. [6, p. 175] – this weaker condition is needed for the random Poincaré map to have a smooth density). The only result that requires the lower bound in (3.4) is Theorem 6.4, which relies on the earlyescape result [17, Theorem 6.4]. See [46] for recent work under weaker assumptions.
Finally we make the following assumption on the noise intensities:

Assume and .
In fact, in the course of the analysis, we will encounter more restrictive conditions of the form , with . The most stringent of these conditions will be needed for the analysis near the folded node, and requires .
The main goal is to establish bounds on the noiseinduced deviation from a deterministic solution. In [15, Theorem 5.1.18], rather precise bounds for the deviation near normally hyperbolic critical manifolds are derived. We want to adapt these to the other phases of motion. As it turns out, the leadingorder effect of noise occurs near the foldednode singularity. Therefore, it will be sufficient to determine the order of magnitude of noiseinduced deviations during other phases of the dynamics, as a function of and .
We fix a deterministic reference solution and set
(3.5) 
As initial condition we choose as it corresponds to . Substituting in (3.1) and Taylorexpanding, we obtain a system of the form
(3.6) 
where
(3.7) 
The nonlinear terms and satisfy as . The matrix of the system linearized around the deterministic solution has the structure
(3.8) 
where and so on, so that in particular , and . Let
(3.9) 
denote the principal solution of the linear system . Then the solution of (3.6) can be written in the form
(3.10) 
and
(3.11) 
In both equations, we expect the stochastic integrals to give the leading contribution to the fluctuations. They can be estimated by the Bernsteintype inequality Lemma A.2. The magnitude of the other integrals can then be shown to be smaller, using a direct estimate which is valid as long as the system does not exit from the region where the nonlinear terms are negligible; see e.g. [17, p. 4826] or [14, Theorem 2.4].
In order to carry out this program, we need estimates on the elements of the principal solution . Note that the components are in principle larger than the components, but this is compensated by the fact that spends most of the time in the vicinity of stable critical manifolds. The following ODEs will play an important rôle:
(3.12) 
Here and . If is bounded away from , standard singular perturbation theory implies that these ODEs admit solutions and of order (and in fact close to ). If approaches or changes sign, this need no longer be the case, but there may still be solutions such that remains small.
Lemma 3.2.
Assume and that the ODEs (3.12) admit solutions such that is bounded for by a function satisfying . Let . Then for sufficiently small ,
(3.13)  
where is the principal solution of the system
(3.14) 
Proof.
Consider the matrix
(3.15) 
Then the equations (3.12) imply
(3.16) 
where the blocks and are given by
(3.17) 
Consider now the variable . If , then (3.16) implies
(3.18) 
The principal solution of this equation is blockdiagonal, with blocks and , where is the principal solution of . The principal solution of the original equation is then given by
(3.19) 
Furthermore, we have
(3.20) 
Computing the matrix product in (3.19) yields the result. Note that more precise expressions for the matrix elements can be obtained if needed. ∎
To describe the size of fluctuations, for given we introduce stopping times
(3.21) 
Proposition 3.3.
Suppose the assumptions of Lemma 3.2 are satisfied with , bounded uniformly in . Given a finite time horizon of order on the slow timescale, there exist constants such that whenever , and ,
(3.22) 
holds for all .
Proof.
Denote by , , the four terms on the righthand side of (3.10). We will start by estimating and . Since , are assumed to be bounded, we may choose of order in (3.13), and is of order , while the other elements of are of order at most. By Lemma A.2 and the bounds on and , there exists a constant such that
(3.23) 
Indeed, to estimate we first use that on any short time interval with , the stochastic process is close to the martingale , defined by
(3.24) 
since remains of order 1 on these time intervals. First using (3.13) and (A0) and then our choice of and Lemma A.1, we see that the martingale’s variance is bounded by
(3.25) 
for some positive constant . Thus the variance is at most of order for all . Now the first inequality in (3.23) follows immediately from the Bernsteintype estimate Lemma A.2. The prefactor in (3.23) simply counts the number of intervals needed to cover , see e.g. [15, Proposition 3.15] for a detailed proof in a simpler, onedimensional setting. The second inequality in (3.23) is shown similarly.
Furthermore, we have for a constant and . From this, together with Gronwall’s lemma, we deduce that there exists a constant such that
(3.26) 
In a similar way, we find
(3.27) 
Choosing small enough, we can ensure that the terms are negligible, and the result follows by taking the sum of the last two estimates. ∎
The size of typical fluctuations is given by the values of for which the probability (3.22) starts getting small, namely and . We conclude that fluctuations have size in the fast direction, and in the slow direction. Note that for simplicity we ignore the logarithmic contribution arising from the prefactor .
We now want to estimate the noiseinduced spreading for the Poincaré map, starting on the section after the folded node, and arriving on the section before the folded node. As described in Section 2 we decompose the map into several maps, see Figure 3, and estimate the spreading for each map separately. This means that we fix an initial condition on each section, and estimate the deviation of the stochastic sample paths from the deterministic solution when it first hits the next section.
3.2 The fast segments
The fast segments are given by and . By Assumption (A4) there exists a slow time of order in which the deterministic solution starting on reaches a neighbourhood of order of the stable critical manifold. In this neighbourhood, the linearization is negative and of order . To reach an neighbourhood of the critical manifold, an additional slow time of order is required. By the dropcurve transversality assumption (A2) and using (A4), it takes another slow time of at most order to reach the section . For we thus have
(3.28) 
for some positive constants . This implies that whenever ,
(3.29) 
for a constant , and furthermore is negative as soon as is larger than a constant times .
Consider now the equations (3.12) for and . We will show that remains bounded on and that there exists a particular solution which also remains bounded on . For , we proceed in two steps:

For , can grow at most by an amount of order .

For , since is negative, we can use standard singular perturbation theory to show that remains of order , and in fact approaches .
For , we change the direction of time and consider the equation
(3.30) 
We know that is negative, bounded away from , except for a time interval of length near . Thus we conclude that there exists a particular solution which remains bounded, of order , on the whole time interval. Therefore Lemma 3.2 shows that , and remain bounded (in norm), of order , and that remains of order for . As a consequence, we can apply Proposition 3.3 as is, with the result that on the section ,

the spreading in the fast direction is of order ,

the spreading in the slow direction is of order .
3.3 The slow segments
The slow segments are given by and . The analysis of the previous subsection can actually be extended to these segments, because is always negative, bounded away from . The conclusions on typical spreading are the same:

the spreading in the fast direction is of order ,

the spreading in the slow direction is of order .
Note that [15, Theorem 5.1.18] provides a more precise description of the dynamics, by constructing more precise covariance tubes. The qualitative conclusion on typical spreading is the same as above.
4 The regular fold
4.1 Approach
The regular fold corresponds to the transition . We again fix a deterministic solution, now starting on . We choose the origin of the coordinate system on the regular fold and the origin of time in such a way that at time , .
Recall from the deterministic analysis (see e.g. [55, 68] for the twodimensional and [78, 67] for the threedimensional case) that, given of order ,

for , the distance of to the critical manifold grows like ;

there exists a such that for ;

there exists a such that reaches order before time .
In this section, we consider the transition , where is a section on which . In this region, the linearization