From random Poincaré maps tostochastic mixed-mode-oscillation patterns

From random Poincaré maps to
stochastic mixed-mode-oscillation patterns

Nils Berglund MAPMO, CNRS – UMR 7349, Université d’Orléans, Fédération Denis Poisson – FR 2964, B.P. 6759, 45067 Orléans Cedex 2, France.    Barbara Gentz Faculty of Mathematics, University of Bielefeld, Universitätsstr. 25, 33615 Bielefeld, Germany.Research supported by CRC 701 “Spectral Structures and Topological Methods in Mathematics”.    Christian Kuehn Vienna University of Technology, Institute for Analysis and Scientific Computing, 1040 Vienna, Austria.
December 22, 2013
Revised, November 14, 2014

We quantify the effect of Gaussian white noise on fast–slow dynamical systems with one fast and two slow variables, which display mixed-mode oscillations owing to the presence of a folded-node singularity. The stochastic system can be described by a continuous-space, discrete-time 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 mixed-mode oscillation patterns. Our analysis shows that there is an intricate interplay between the number of small-amplitude oscillations and the global return mechanism. In combination with a local saturation phenomenon near the folded node, this interplay can modify the number of small-amplitude oscillations after a large-amplitude oscillation. Finally, sufficient conditions are derived which determine when the noise increases the number of small-amplitude 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, mixed-mode oscillation, return map, random dynamical system, first-exit 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 mixed-mode oscillations (MMOs) which are patterns consisting of alternating structures of small-amplitude and large-amplitude 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–Huxley-type 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 small-amplitude oscillations (SAOs) while a global return mechanism leads to large-amplitude oscillations (LAOs). In this introduction, we shall just outline the main ideas; the precise development of our set-up and results starts in Section 2.

For a deterministic trajectory, we can symbolically write an MMO as a sequence


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 folded-node singularities [24] which are generic in three-dimensional ODEs with one fast and two slow variables [77, 7]. For the global return mechanism, one frequently encounters a relaxation-type structure induced by a cubic (or S-shaped) fast-variable nullcline, also called the critical manifold, which was studied extensively already in the context of van der Pol-type oscillators; see e.g. [51, 35, 20, 24] and the references therein. Non-degenerate folds, folded-node singularities and S-shaped 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 saddle-nodes [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 low-dimensional 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 folded-node singularity and an S-shaped 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 two-dimensional section . Deterministic return maps in the presence of folded-node singularities have been analyzed, e.g., in [37, 54]. Although the two-dimensional Poincaré map is invertible, the strong contraction near attracting critical manifolds implies that it is close to a one-dimensional, usually non-invertible map. Figure 1 (a) shows an example of such a one-dimensional deterministic return map . The apparent discontinuities are in fact points where the map’s graph displays very narrow dips, due to the presence of so-called 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].

Figure 1: -coordinate of the first-return map on the section  for the Koper model (7.1). Parameter values are , , , , with noise intensities (a) , (b) , (c) , and (d) . The horizontal and vertical lines mark the location of canards.

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 continuous-space, discrete-time 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


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 non-intuitive 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 noise-induced 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 folded-node singularity.

The structure of this article is the following: After introducing the deterministic set-up in Section 2, we provide first estimates on noise-induced 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:

  • Theorem 6.1 (global return map) quantifies the effect of noise during the global return phase;

  • Theorem 6.2 (local map for inner sectors) provides estimates on noise-induced fluctuations for orbits starting near a folded node in sectors with small SAO number; together with Theorem 6.1 it yields bounds on the size of fluctuations of the Poincaré map in all inner sectors;

  • Theorem 6.4 (local map for outer sectors) gives similar estimates for orbits starting in sectors with a large SAO number; in particular, it proves the saturation effect.

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.











Figure 2: A singularly perturbed Markov chain.

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, non-vanishing 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 spectral-theoretic 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 Marie-Curie International Re-integration 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 Mixed-mode oscillations – The setup

In this section, we shall outline a typical setup for deterministic mixed-mode oscillations based upon three-dimensional fast–slow systems of the form


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


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 S-shaped 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,


    where the two-dimensional submanifolds are normally hyperbolic attracting, while the two-dimensional submanifold is normally hyperbolic repelling, i.e.,


    and are one-dimensional smooth fold curves consisting of generic fold points


    Without loss of generality we assume from now on that for all .

Figure 3: Sketch illustrating the definition of the different sections. The horizontal coordinate is , the vertical one is , and points out of the plane.

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


which is solved by the so-called slow flow. Differentiating implicitly with respect to yields


for the slow flow, so that the slow subsystem (2.6) can be written as


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 implicit-function 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 right-hand side of (2.8) by and applying a rescaling of time. It reads


We make the following further assumptions:

  • Suppose all fold points on satisfy the normal switching condition [67, 78]


    Furthermore, assume that the projections of along the -coordinate onto , which are also called the drop curves, are transverse to the slow flow.

  • Assume that the normal switching condition fails only at a unique singularity and is a node equilibrium point of (2.9); in this case, is called a folded node (or folded-node singularity) [9, 77].

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


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


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 well-defined maps from to .

  • The geometry of the flow-induced 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 slow-flow region near , and

  • the non-degenerate 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 mixed-mode oscillations under the assumptions (A0)–(A4) are well-known; 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 S-shaped critical manifold induces the LAOs. Fixed points of a full return map, say , correspond to MMOs with a certain pattern


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 45 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


where is a -dimensional standard Brownian motion on a probability space . The maps


may depend on , and are assumed to be and to satisfy the usual bounded-growth 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


    be the diffusion matrix. There exist constants such that

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 early-escape 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 noise-induced 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 leading-order effect of noise occurs near the folded-node singularity. Therefore, it will be sufficient to determine the order of magnitude of noise-induced deviations during other phases of the dynamics, as a function of and .

We fix a deterministic reference solution and set


As initial condition we choose as it corresponds to . Substituting in (3.1) and Taylor-expanding, we obtain a system of the form




The nonlinear terms and satisfy as . The matrix of the system linearized around the deterministic solution has the structure


where and so on, so that in particular , and . Let


denote the principal solution of the linear system . Then the solution of (3.6) can be written in the form




In both equations, we expect the stochastic integrals to give the leading contribution to the fluctuations. They can be estimated by the Bernstein-type 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:


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 ,


where is the principal solution of the system


Consider the matrix


Then the equations (3.12) imply


where the blocks and are given by


Consider now the variable . If , then (3.16) implies


The principal solution of this equation is block-diagonal, with blocks and , where is the principal solution of . The principal solution of the original equation is then given by


Furthermore, we have


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

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 ,


holds for all .


Denote by , , the four terms on the right-hand 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


Indeed, to estimate we first use that on any short time interval with , the stochastic process is close to the martingale , defined by


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


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 Bernstein-type 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, one-dimensional 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


In a similar way, we find


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 noise-induced 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 drop-curve transversality assumption (A2) and using (A4), it takes another slow time  of at most order  to reach the section . For we thus have


for some positive constants . This implies that whenever ,


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


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 two-dimensional and [78, 67] for the three-dimensional 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