Mixed-Mode Oscillations in a Stochastic, Piecewise-Linear System

# Mixed-Mode Oscillations in a Stochastic, Piecewise-Linear System

D.J.W. Simpson and R. Kuske
Department of Mathematics
University of British Columbia
Vancouver, BC, V6T1Z2
The authors acknowledge support from an NSERC Discovery Grant.
###### Abstract

We analyze a piecewise-linear FitzHugh-Nagumo model. The system exhibits a canard near which both small amplitude and large amplitude periodic orbits exist. The addition of small noise induces mixed-mode oscillations (MMOs) in the vicinity of the canard point. We determine the effect of each model parameter on the stochastically driven MMOs. In particular we show that any parameter variation (such as a modification of the piecewise-linear function in the model) that leaves the ratio of noise amplitude to time-scale separation unchanged typically has little effect on the width of the interval of the primary bifurcation parameter over which MMOs occur. In that sense, the MMOs are robust. Furthermore we show that the piecewise-linear model exhibits MMOs more readily than the classical FitzHugh-Nagumo model for which a cubic polynomial is the only nonlinearity. By studying a piecewise-linear model we are able to explain results using analytical expressions and compare these with numerical investigations.

## 1 Introduction

Oscillatory dynamics involving oscillations with greatly differing amplitudes, known as mixed-mode oscillations (MMOs), see Fig. 1, are important in neuron models [1] and in a multitude of chemical reactions [2, 3], refer to [4] for a recent review. Yet there are many open questions regarding the creation, robustness and bifurcations of MMOs. A variety of mechanisms generate MMOs in deterministic systems. Alternatively MMOs may be noise-induced; there are also several scenarios by which this may occur.

We study the following form of the FitzHugh-Nagumo (FHN) model with small, additive, white noise:

 dv=(f(v)−w) dt,dw=ε(αv−σw−λ) dt+D dW, (1)

where represents a potential, is a recovery variable and is a standard Brownian motion. The FHN model is used as a prototypical model of excitable dynamics in a range of scientific fields [5, 6]. Here is a positive constant and , which is regarded as the main bifurcation parameter, controls the growth of oscillations, as seen below. The small parameter , represents the time-scale separation and is the noise amplitude (). Values of and used in, for instance [7, 8], are no larger than the values considered here. By scaling we may assume , except in the special case which corresponds to the van der Pol model (and in this case we may further assume ). We assume that is continuous and roughly of cubic shape. For simplicity, we assume that has a local minimum at and a local maximum at , regardless of the precise function chosen.

If is a cubic, as originally taken by FitzHugh [9] and Nagumo et. al. [10], then, by the above requirements, the cubic must be

 f(v)=3v2−2v3. (2)

Fig. 2-A illustrates the role of the parameter for (1) with (2) in the absence of noise. A small amplitude periodic orbit is created in a Hopf bifurcation at . For the parameters used in Fig. 2, this periodic orbit is stable and its amplitude increases with . Near the amplitude increases exponentially. This rapid growth is known as a canard explosion and is due to time-scale separation and global dynamics [11, 12, 13, 14]. The value of the canard point, , which is well-defined for smooth systems [15, 16], decreases to zero with , as shown in Fig. 3-A. Over an order range of values, (1) with (2) may either settle to equilibrium, exhibit small amplitude oscillations, or exhibit large amplitude oscillations (relaxation oscillations).

As in [17, 18], here we study a piecewise-linear (PWL) FHN model so that, in the presence of noise, the system is amenable to a rigorous analysis without the need for an approximation or limiting scenario. PWL models are commonly used in circuit systems [19, 20, 21]. A PWL version of a driven van der Pol oscillator is studied in [22] to explain the breakdown of canards in experiments. We consider the continuous, PWL function

 f(v)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ηLv,v≤0η1v,01, (3)

where , , and

 η1=w1v1,η2=1−w11−v1. (4)

We state the particular form here in order to briefly illustrate key differences between the smooth and PWL FHN models. Further motivation for (3) is given in §2 and shown in Fig. 5, As shown in Fig. 2-B, (1) with (3) may exhibit a canard explosion. The canard point, , is not well-defined for this system because it lacks global differentiability. Instead we consider the values, and , at which the maximum -value of the periodic orbit of (1) with (3) in the absence of noise is and respectively. The piecewise nature of (3) leads to a natural classification of periodic orbits and oscillations of (1) with (3). (Typically we refer to one complete revolution about the equilbrium as a single oscillation.) With Fig. 2-B in mind, if is the maximum -value of a periodic orbit or single oscillation we declare that the orbit or oscillation is

 small if 01. (5)

Fig. 3-B illustrates typical dependence of and on . In particular we notice that for a fixed choice of the slopes, , in (3), the PWL version of the FHN model does not exhibit small oscillations for arbitrarily small . This is because the two eigenvalues associated with the equilibrium for small are real-valued for sufficiently small negating the possibility of small oscillations, see §2. Unlike for the van der Pol model, values of that are relevant for the FHN model are usually sufficiently large for small oscillations to be important in the PWL model.

The effect of noise in (1) has seen significant recent attention, see for instance [23, 24, 25]. Noise may induce regular oscillations in (1) when in the absence of noise there are no oscillations. There is more than one mechanism that may cause this, most notably stochastic resonance [26] (when a small periodic forcing term is present in addition to noise), coherence resonance [25] (usually when the system is quiescent in the absence of noise), and self-induced stochastic resonance [27] (involving relatively large noise that drives oscillations of periods different from that of the deterministic system).

If is tuned to values near the canard explosion, in the presence of noise the system may exhibit both small amplitude and large amplitude oscillations, i.e. MMOs, as shown in Figs. 4 and 1. Similar MMOs are described in [28] for very small noise by a careful choice of parameter values. For a version of (1) that contains nonlinearity in the equation to better mimic neural behaviour, it has been observed that when is chosen to be just prior to the canard point the frequency of relaxation oscillations increases with noise amplitude [8]. Noise-induced MMOs have been described for three coupled FHN systems near a canard [29]. A signal-to-noise ratio may be defined to quantitatively determine dominant frequencies [30]. Noise-induced MMOs may arise via a different mechanism in the case that the Hopf bifurcation is subcritical [31]. There are a variety of mechanisms for MMOs in three-dimensional systems that we do not consider here, see for instance [32, 33] and references in [4].

In this paper we study noise-driven MMOs in (1) with (3). We use analytical methods to identify parameter values for which MMOs occur and describe the dependence of each model parameter on MMOs. For typical values of the noise amplitude, , MMOs occur over some interval of positive -values. In order to find such intervals we determine exit distributions for forward orbits of (1) with (3) through various cross-sections of phase space. The exit distributions allow us to deduce the amplitude of oscillations and consequently find intervals of MMOs. We show that MMOs are robust in the sense that large variations in other model parameters can have minimal effect on the width of the -intervals.

We note that the model we consider has additive noise in the -equation only, as in, for instance, [8, 25]. This choice allows some simplifications in demonstrating the analytical method, while still capturing qualitatively the behavior that would be observed for more general additive noise. Throughout the paper we indicate where this assumption allows some simplification in the analysis, and we indicate the differences that would need to be addressed for the case of noise also in the -equation.

The remainder of the paper is organized as follows. Section 2 briefly overviews PWL FHN models and provides an analysis of (1) with (3) in the absence of noise. Here we explain the Hopf-like bifurcation at that creates stable oscillations and describe equations for and , Fig. 3. Calculations of exit distributions are detailed in §3. Here we also describe the method by which we use these distributions to find parameter values corresponding to MMOs. Section 4 combines the analysis of the previous sections to determine the effect of each model parameter on MMOs. Finally conclusions are presented in §5.

## 2 Properties of the deterministic system

Analytical results may be derived for (1) when is a PWL function. Arguably the simplest continuous, PWL function that one can use for consists of three line segments (one of them being the straight connection between and ). The FHN model with this function is well-studied [34, 35], refer to [36] for the van der Pol system. However, with this three-piece PWL function, (1) does not exhibit a canard, as shown in [37], and so we do not consider it further. Consequently, as in [22, 38], we use two line segments between and denoting the intermediate point by and the slopes by , specifically (3), as shown in Fig. 5. If instead contains multiple line segments left of such that the slopes of the two lines meeting at are , multiple coexisting attractors commonly exist for small which leads to complications that we do not study here. For simplicity we do not consider comprised of more than four line segments. For canards in PWL FHN models with many segments we refer to reader to the recent work of Rotstein et. al. [38].

In the absence of noise (i.e. when ), (1) with (3) is a continuous, two-dimensional, PWL, ordinary differential equation system:

 ˙v=f(v)−w,˙w=ε(αv−σw−λ). (6)

The phase space, , is divided into four regions

 (7)

by the three switching manifolds, , and , on which the system is non-differentiable.

Each linear component of (6) with (3) has a unique equilibrium, , see Fig. 5 (unless in which case the relevant - and -nullclines are parallel). In the terminology of piecewise-smooth dynamical systems, each is either admissible (lies in the closure of ) or virtual (lies outside the closure of ). The Jacobian, , and the eigenvalues, , associated with each are

 Aj=[ηj−1εα−εσ], (8) ρj=12(ηj−εσ±√(ηj+εσ)2−4εα). (9)

We assume

 ηL<−εσ−2√εα,εσ<η1<ασ, (10)

such that is an attracting node and is either a repelling node or a repelling focus as determined by the sign of . The restriction (10) ensures that stable oscillations are created at , as shown below.

The bifurcation at that results from the interaction of an equilibrium with the switching manifold, , is an example of a discontinuous bifurcation [39, 40, 41]. Effectively, eigenvalues that determine the stability of the admissible equilibrium change discontinuously as the equilibrium crosses the switching manifold at . In general, a bifurcation is expected to occur if one or more eigenvalues “jump” across the imaginary axis at the crossing. Such a bifurcation may be analogous to a smooth bifurcation or it may be unique to piecewise-smooth systems [40]. For two-dimensional systems, codimension-one, discontinuous bifurcations involving a single smooth switching manifold have been completely classified [41, 42].

For the PWL system (6) with (3), an attracting periodic orbit is born at the discontinuous bifurcation, . The relative size of the periodic orbit for small is dependent upon whether the equilibrium, , is a node or a focus. If is an attracting node and is a repelling node, invariant lines corresponding to eigenvectors prevent the creation of a local periodic orbit corresponding to a small oscillation [41]. The periodic orbit generated at has large amplitude (corresponding to a relaxation oscillation). Specifically, as , the maximum value of of the periodic orbit limits on a value greater than .

If instead is a repelling focus, then the bifurcation is a discontinuous analogue of a Hopf bifurcation in that a periodic orbit is created locally. Unlike for a classical Hopf bifurcation, the periodic orbit grows in size linearly with respect to (see Fig. 2-B) which is typical for piecewise-smooth systems.

The value of for which the square-root term in (9) vanishes is the critical value of (see Fig. 3-B) above which the periodic orbit created at is small and below which this orbit is large, and is given by

 εcrit=1σ2(2α−ση1−2√α(α−ση1)). (11)

The curves and , Fig. 3-B, which bound the region of medium oscillations, emanate from . Since the underlying system is PWL, we may obtain analytical expressions relating to these curves by deriving the explicit solution to the flow of each linear component of (1) with (3). We let denote the solution to the linear component of (6) with (3) corresponding to , for an arbitrary initial condition, . For instance:

 [v(1)(t;v0,w0)w(1)(t;v0,w0)] = (12) −1ω1sin(ω1t)cos(ω1t)−η1+εσ2ω1sin(ω1t)⎤⎦[v0−v∗1w0−w∗1]+[v∗1w∗1],

which equals the solution to (6) with (3) for the same initial condition whenever lies in the closure of at all times between and , and where

 ωj=12√∣∣(ηj+εσ)2−4εα∣∣. (13)

Unfortunately we cannot in general explicitly solve (12) for (in particular solve: for ). Consequently we are unable to extract or explicitly in terms of the parameters of the system. For brevity we omit the details and simply note that for the figures in this paper we determine and by numerically solving transcendental expressions. This may be accomplished to any desired accuracy quickly and does not require the use of a differential equation solving method.

The following two approximations are used in the analysis of later sections. For a wide range of parameter values the attracting periodic orbit passes close to the origin. If we approximate by finding where the next intersection of the forward orbit of with the -nullcline is , then we obtain

 λv1≈αv1−σw11+e(η1−εσ)π2ω1, (14)

which is particularly accurate for , as shown in Fig. 3.

Second, the two eigenvalues associated with (9) are and . Within , trajectories rapidly approach the associated slow eigenvector. This eigenvector intersects the switching manifold, , at

 ^wL=λρL,slowα−σηL. (15)

Consequently, trajectories such as large oscillations that spend a relatively long period of continuous time in , exit this region extremely close to the point . This point is important below in the discussion of stochastic dynamics. It is usually sufficient to approximate by considering and the subsequent and finding the value of where intersects . This is because corresponds to the existence of a periodic orbit with a maximum -value of , which must intersect , see (5) and the surrounding discussion.

## 3 Exit distributions

To analyze noise-driven MMOs we consider solutions to (1) with (3) in the presence of noise over long time frames such that transient behaviour has decayed. In this context we determine the fraction of oscillations that are small, the fraction that are medium, and the fraction that are large, referring to (5). One method is to simply solve the system for a long time and count the number of different oscillations. This Monte-Carlo approach is useful for obtaining a basic understanding of the system but poor for an accurate quantitative analysis because the system must be solved accurately for many parameter combinations requiring considerable computation time. Instead, since the system under consideration is PWL, we are able to use exit distributions for the regions (7) to approximate these fractions. This approach does not necessitate arbitrarily small . In contrast, Muratov and Vanden-Eijnden [7] applied stochastic methods to (1) with (2) by considering the system asymptotically (i.e. with arbitrarily small and ) which essentially reduces the problem to one dimension. In [43], the same system is considered but in the limit which also reduces mathematical calculations to one dimension.

Here we describe the exit distributions for forward orbits of (1) with (3) along various cross-sections of phase space. In the following section we use these exit distributions to identify MMOs. The four cross-sections we consider are:

 Σ1={(0,w) | w<0}∪{(v,η1v) | 0≤vw1}∪{(v,η1v) | v∗10}∪{(v,η1v) | 0≤v

as depicted in Fig. 6. We exclude the switching manifold, , from calculations because large oscillations follow a sufficiently predictable path back to when .

One method for computing a first exit distribution is to solve the Fokker-Planck equation for the probability density of the process (1) with (3) and absorbing boundary conditions [44, 45]. Integration of the solution to this boundary value problem at the boundaries in an appropriate manner and over all positive time, may yield the desired exit distribution. However we dismiss this approach as it necessitates extensive numerical computations, in part because drift dominates the diffusion which typically requires extra attention [46, 47, 48]. Instead we utilize the fact that within each region, , the system is Ornstein-Uhlenbeck and, ignoring switching manifolds, has a known explicit solution [45].

The transitional probability density, , i.e. , for the solution to the Ornstein-Uhlenbeck process of , i.e. (1) with , after a time is the Gaussian

 p(1)t(v,w|v0,w0)=12π√det(Θ(t))exp(−12ΔzTΘ(t)−1Δz), (17)

where

 Δz(t;v0,w0)=[v−v(1)(t;v0,w0)w−w(1)(t;v0,w0)]. (18)

The mean, , is the solution to the system in absence of noise (12) and is the covariance matrix given by:

 Θ(t)=D2⎡⎢ ⎢⎣∫t0(eA1s12)2ds∫t0eA1s12eA1s22ds∫t0eA1s12eA1s22ds∫t0(eA1s22)2ds⎤⎥ ⎥⎦=D2[[]ccθ11(t)θ12(t)θ12(t)θ22(t)], (19)

where denotes the -component of the matrix exponential of (8) and we have introduced the for convenience. (Note that (19) would contain more terms if (1) also included noise in the equation.) The probability density (17) obeys the Fokker-Planck equation

 ∂pt(v,w|v0,w0)∂t=−∇⋅Jt(v,w|v0,w0), (20)

where

 Jt(v,w|v0,w0)=[(η1v−w)ptε(αv−σw−λ)pt−D22∂pt∂w], (21)

is the probability current [44, 45]. By integrating (20) and applying the divergence theorem, it follows that the net flow of probability across, say, the -nullcline between and , is given by

 ∫v1v∗1n⋅Jt(v,η1v|v0,w0)dv,

where is the normal vector of the -nullcline pointing outwards [44, 45], i.e. here . If trajectories were unable to cross the -nullcline more than once, then the integral

 ∫∞0n⋅Jt(v,η1v|v0,w0)dt, (22)

would be equal to the density of the first (and last) exit points for escape from below the -nullcline. However, trajectories have multiple intersections with the -nullcline due to the presence of noise. By considering two different time frames, we now show that these multiple intersections have a negligible effect and that (22) represents an exit distribution suitable for our analysis. Specifically we first show that the probability of return to the -nullcline after a short time is small. Then we show that within this short time frame points of multiple intersections are clustered. Finally we show that for longer time intervals after an intersection with the -nullcline, trajectories are far from the nullcline, assuming small noise levels.

We first look at return times for the -nullcline. Closed form expressions for first passage problems of multi-dimensional Ornstein-Uhlenbeck processes are not straight-forward [49, 50]; for this reason we simplify to a one-dimensional problem. Consider the forward orbit of a point on the -nullcline with . Using to represent the distance from the nullcline, (1) with may be written as

 dv=−y dt,dy=((η1−σε)y+(α−ση1)(v−v∗1)ε) dt+D dW, (23)

where we have substituted . Intersections of the orbit with the nullcline are determined by the equation of (23) which we conservatively reduce to

 dy=c dt+D dW, (24)

where, for , the magnitude of the drift has a lower bound:

 c≥(α−ση1)(vmin−v∗1)ε,

assuming for some . For any , we are interested in , i.e. the probability that a solution to (24) with satisfies at some . To calculate this probability we let denote the transitional probability density for (24) with , and condition over the event that , for all :

Notice,

 Pr(y(t)=0 for some t≥δ ∣∣ y(δ)=z)=Pr(y(t)=−z for some t≥0 ∣∣ y(0)=0),

because (24) has no explicit dependence on and . This enables us to write

 (25)

where

 (26)

can be calculated from the density of the first hitting time of to (refer to [51, 52] for more details) producing

 G(z)=c∫∞0p(z,t)dt. (27)

By using (27) and evaluating the integral on the right-hand side of (25), we obtain

 (28)

For instance with , and , whenever the probability of return to the -nullcline after a time of is less than .

Second, for the system (1) with we look at the distribution of future -nullcline intersections up to a time . The solution (17) with evaluated on the -nullcline and normalized is a Gaussian with mean and variance:

 ~v(t) = −(η1θ12−θ22)v(1)+(η1θ11−θ12)w(1)θ22−2η1θ12+η21θ11, (29) ~σ(t)2 = detΘθ22−2η1θ12+η21θ11, (30)

respectively, where the were defined in (19). We observe that (22) undergoes negligible change when convolved by the Gaussian with (29) and (30) evaluated at . For this reason we use (22) to compute exit distributions on the -nullclines. The absence of noise in the equation of (1) ensures multiple rapid crossings through the switching manifolds are not permitted. Consequently we use an integral similar to (22) for exit distributions across the other switching manifolds also. This accounts for all components of each (16).

When the equilibrium, , is admissible, we expect the forward orbit of any point on to escape the lower half of (below the -nullcline) through . We calculate the exit distribution of the orbit through with (22). (Note, for simplicity we omit the term in (21) when using (22) because it is dominated by the other terms in .) Using equally spaced data points and performing this calculation repeatedly, we determine the exit distribution on for any probability density of points on . From the exit distribution on we continue in a similar fashion and compute the exit distributions on , and lastly . Note these calculations use analytical expressions like (22) and not Monte-Carlo simulations. Numerically we observe that the iterative procedure of mapping a distribution on to itself (through , and ) approaches the limiting distribution of the intersection of an arbitrary forward orbit of the system with , Fig. 7. We use the limiting distributions on and to calculate the probability that an arbitrary oscillation is small, medium or large. The results for a range of parameter values are given in the next section.

## 4 Mixed-mode oscillations

In order to understand MMOs quantitatively, we say that (1) with (3) exhibits MMOs whenever both small and large oscillations occur at least of the time. Specifically we find where exit densities corresponding to small and large oscillations both integrate to a value greater than . Fig. 8 illustrates the dependence of MMOs on the primary bifurcation parameter, , and the noise amplitude, . Roughly the range of values which permit MMOs increases with . This matches our intuition, more noise allows for a wider variety of oscillations. We compute Fig. 8 using the iterative scheme described in §3; Monte-Carlo simulations (not shown) give good agreement.

MMOs exist in a region bounded on the left by the curve along which large oscillations occur of the time and on the right by the curve along which small oscillations occur of the time. When the former curve has the value , and the latter curve has the value . This is because the periodic orbit of the system in the absence of noise, (6), changes from small to medium at , and from medium to large at , §2.

From Fig. 8, we see that MMOs do not occur for arbitrarily small even near the canard explosion in contrast to what may be expected. This is because for very small and , medium oscillations dominate. Medium oscillations occur less frequently with increasing . Note also that the MMO regions appear relatively symmetric with respect to .

For the smooth system (1) with (2), we may roughly compute the region of MMOs, by the above definition, from Monte-Carlo simulations, Fig. 8-A. We see that MMOs occur over a smaller parameter range for the smooth version of the FHN model. This distinction is possibly explained by Fig. 4. For the PWL model, all oscillations (including small oscillations) spend sufficient time in to be drawn into the slow eigenvector of this region. Small and large oscillations are intertwined in on their approach to ; the amplitude of one oscillation is practically independent of the previous oscillation. (From a numerical viewpoint, in this situation fewer data points are required than in general.) In contrast, for the smooth system there is a significant distance between small and large oscillations and therefore more noise is required for, say, a large oscillation to follow a small oscillation.

Noting this difference between the smooth and PWL models, we considered another parameter range for the PWL system that has a different exit distribution near the origin. For small values of the slope, , still respecting (10), the MMOs may include small oscillations that do not enter , so that the exit distribution across may be bimodal, as shown in Fig. 7-B. Here large oscillations intersect near (15) whereas the majority of small oscillations intersect on the -nullcline. We considered whether this type of bimodal exit distribution on plays a role analogous to distance between small and large oscillations in the smooth model, but we did not see any evidence of this effect. Specifically the MMO region, Fig. 8-B, has a similar size and shape to the region in Fig. 8-A for which the corresponding value of is an order of magnitude larger.

The boundaries of the MMO regions shown in Fig. 8 are relatively linear, hence we perform an analytical calculation of the slopes at . Let [] denote the slope, , at , of the curve along which of oscillations are small [large].

Let us begin with the curve along which exactly 10% of oscillations are small. This curve intersects at at which the attracting periodic orbit created at intersects at . Here we can focus on small oscillations only, so it suffices to consider the linear systems of and , i.e. (1) with

 f(v)={ηLv,v≤0η1v,v>0. (31)

As is increased, the maximum -value of the deterministic periodic orbit increases at a rate, say, . Due to linearity, this rate is given simply by

 κ1=v1λv1. (32)

When , the periodic orbit intersects at some point with and the line at . If we now consider small but leave all other parameters unchanged, over a long time frame trajectories intersect at points approximately normally distributed about . Since small oscillations neglect switching of (1) at , intersection points on are similarly approximately normally distributed about , as shown in Fig. 9, with a standard deviation of say, , where is a constant that we compute below. The -value of intersection points on then have the distribution , using (32). That is, if denotes the probability density for these -values, then

 qsmall(v)=1√2πγ1De−12γ21D2(v−v1−κ1(λ−λv1))2. (33)

Then of oscillations are small when

 ∫v1−∞qsmall(v)dv=12(1−erf(κ1(λ−λv1)√2γ1D))=0.1. (34)

By rearranging the previous equation we deduce that the slope of the curve at is

 ssmall=dDdλ=κ1√2γ1 erf−1(0.8). (35)

From §2, may be accurately calculated by solving transcendental equations.

We obtain a good approximation to as follows. Due to strong contraction in , the distribution of points on has a standard deviation that is much smaller than the standard deviation of points on . Than it is reasonable to approximate the distribution on by the single value . Then is equivalent to the exit distribution along the -nullcline, thus by (22),

 qsmall(v)=˙w(v,η1v)∫∞0p(1)t(v,η1v|0,wλv1)dt+O(D2), (36)

where refers to (6). By (17) and (19),

 qsmall(v)=˙w(v,η1v)2πD2∫∞01√θ11θ22−θ212e−ϕ(v,t)D2dt+O(D2), (37)

where

 ϕ(v,t)=12(θ11θ22−θ212)[v−v(1),η1v−w(1)][θ22(t)−θ12(t)−θ12(t)θ11(t)][v−v(1)η1v−w(1)]. (38)

In the limit , the asymptotic approximation to integral in (37) is determined from the main contribution of , which is its maximum value; formally this is achieved by Watson’s lemma [53]. We omit the details of this calculation which produces

 γ1=√θ1,1. (39)

We calculate the slope of the curve on which of oscillations are large at in a similar fashion. Again we approximate the density of intersection points on by a point mass but this time we use the value (15) and compute the density of intersection points on the switching manifold, . For small this density is approximately Gaussian, i.e. , where when the deterministic trajectory passes through the points (or rather very near to this point), and . If denotes this probability density, then

 qlarge(w)=1√2πγ2De−12γ22D2(w−wλ1−κ2(λ−λ1))2. (40)

The constant, , may be computed from (12) using the chain rule for differentiation:

 κ2=∂w(1)∂λ=∂w(1)∂t∣∣∣t=tint∂tint∂λ+∂w(1)∂w0∂w0∂λ+∂w(1)∂v∗1∂v∗1∂λ+∂w(1)∂w∗1∂w∗1∂λ. (41)

where is the time taken for the trajectory to go from to and

 ∂tint∂λ=(∂v(1)∂w0∂w0∂λ+∂v(1)∂v∗1∂v∗1∂λ+∂v(1)∂w∗1∂w∗1∂λ)/∂v(1)∂t∣∣∣t=tint. (42)

By a calculation similar to that for described above, we obtain

 γ2=√θ11(˙w˙v)2−2θ12(˙w˙v)+θ22, (43)

using (6).

Unlike small oscillations, large oscillations traverse and so we must also consider the flow in these regions. For and near , we let denote the first intersection of the backwards orbit from with so that we may distinguish medium and large oscillations on when . We write

 ^w2=wλ1+κ3(λ−λ1), (44)

ignoring higher order terms and where may be calculated in a manner similar to . For small , the forward orbit of any point , with , has the probability, , of undergoing a large, rather than medium, oscillation before returning to . We find that for small there is a very sharp transition of at so that it suffices to use the approximation , where for and for . Consequently, of oscillations are large when

 ∫wλ1+κ3(λ−λ1)−∞qlarge(w)dw=12(1−erf((κ2−κ3)(λ−λ1)√2γ2D))=0.1, (45)

where is given by (40). By rearranging this expression we arrive at

 slarge=dDdλ=κ2−κ3√2γ2 erf−1(0.8). (46)

We have verified that computation of and by the expressions (35) and (46) matches with computation by exit distributions (22) as in Fig. 8.

The variation of the slopes (35) and (46) with respect to , and is shown in Fig. 10. Both (35) and (46) approach zero as , but accounting for the fact that the noise amplitude, , is not multiplied by in (1), the scaled values and vary relatively slightly. With increasing , and increase slightly and approach the same value; with increasing , and decrease slightly. Due to linearity, if is held constant and is increased, is unchanged. We do not need to consider variation in and because these values may be written in terms of and (4). Additionally, the slopes are not strongly affected by (as long as ) and . Therefore, for the most part, and lie in, say, the interval . Hence for intermediate values of although the interval of values which permit MMOs varies widely with the system parameters, the width of this interval is robust with respect to parameter change. For small values of , MMOs may not occur at all. For very large values of , multiple crossings on the may generate different results.

## 5 Conclusions

We have studied MMOs in a PWL version of the FHN model, (1) with (3). To obtain quantitative results we have defined oscillations as small, medium, or large by the maximum -value attained (5). Furthermore we define MMOs by where at least of oscillations are small and at least are large (this approach may be applied to any preferred values of the percentages). We incorporate noise additively in one equation for transparency of analysis. Numerically we have observed that additive noise in both equations yields similar mixed-mode dynamics.

The existence of a canard explosion in the system with no noise is a consequence of incorporating a four-piece PWL function into the model. We have introduced an analogy for canard points of smooth systems, specifically we identify two values, and , at which the periodic orbit of the deterministic system changes from small to medium, and from medium to large. Both values increase with increasing , and , and decrease with increasing as shown in Fig. 10.

Near the canard explosion noise drives MMOs. The boundaries of the regions of MMOs shown in Fig. 8 are approximately linear for small and we have calculated their slopes, , at . For small values of we have shown that the probability current provides an approximation for exit distributions; for large values of this approximation may no longer be valid. Unless the noise amplitude is extremely small, MMOs exist over some interval of values. We have illustrated that for constant typically the width of this interval changes minimally with a relatively large variation in the values of the other system parameters.

For the results in this paper we have used the value for the slope of the -nullcline for . With similar or more negative values of the system exhibits the same qualitative behaviour. However, with larger values of , say , (1) with (3) may have multiple attracting solutions in the absence of noise. The coexistence of attracting small and large periodic orbits, naturally produces MMOs in the presence of noise though this is via a different mechanism than the one studied here.

The FHN model of Makarov et. al. [8], which includes additional nonlinearity, readily exhibits MMOs. It is possible that this addition strengthens the attraction of the periodic orbit for , mimicking the slow eigenvector discussed above and causing MMOs to be more robust.

## References

• [1] I. Erchova and D.J. McGonigle. Rhythms of the brain: An examination of mixed mode oscillation approaches to the analysis of neurophysiological data. Chaos, 18:015115, 2008.
• [2] D. Barkley. Slow manifolds and mixed-mode oscillations in the Belousov-Zhabotinskii reaction. J. Chem. Phys., 89(9):5547–5559, 1988.
• [3] V. Petrov, S.K. Scott, and K. Showalter. Mixed-mode oscillations in chemical systems. J. Chem. Phys., 97(9):6191–6198, 1992.
• [4] M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H.M. Osinga, and M. Wechselberger. Mixed-mode oscillations with multiple time scales. Under review, 2010.
• [5] C. Rocsoreanu, A. Georgescu, and N. Giurgiteanu. The FitzHugh-Nagumo Model: Bifurcation and Dynamics. Kluwer, Norwell, MA, 2000.
• [6] J. Keener and J. Sneyd. Mathematical Physiology. Springer-Verlag, New York, 1998.
• [7] C.B. Muratov and E. Vanden-Eijnden. Noise induced mixed mode oscillations in a relaxation oscillator near the onset of a limit cycle. Chaos, 18:015111, 2008.
• [8] V.A. Makarov, V.I. Nekorkin, and M.G. Velarde. Spiking behavior in a noise-driven system combining oscillatory and excitatory properties. Phys. Rev. E, 86(15):3431–3434, 2001.
• [9] R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophys. J., 1(6):445–466, 1961.
• [10] J. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simulating nerve axon. Proc. Inst. Radio Eng., 50(10):2061–2070, 1962.
• [11] E. Benot, J.F. Callot, F. Diener, and M. Diener. Chasse au canard. Collect. Math., 32(1-2):37–119, 1981. (in French).
• [12] W. Eckhaus. Relaxation oscillations including a standard chase on french ducks. In Asymptotic Analysis II., volume 985 of Lecture Notes in Mathematics, pages 449–494. Springer-Verlag, New York, 1983.
• [13] S.M. Baer and T. Erneux. Singular Hopf bifurcation to relaxation oscillations. SIAM J. Appl. Math., 46(5):721–739, 1986.
• [14] S.M. Baer and T. Erneux. Singular Hopf bifurcation to relaxation oscillations. II. SIAM J. Appl. Math., 52(6):1651–1664, 1992.
• [15] M. Krupa and P. Szmolyan. Relaxation oscillation and canard explosion. J. Diff. Eq., 174:312–368, 2001.
• [16] M. Krupa and P. Szmolyan. Extending geometric singular perturbation theory to nonhyperbolic points - fold and canard points in two dimensions. SIAM J. Math. Anal., 33(2):286–314, 2001.
• [17] S. Coombes. Neuronal networks with gap junctions: A study of piecewise linear planar neuron models. SIAM J. Appl. Dyn. Sys., 7(3):1101–1129, 2008.
• [18] A. Tonnelier and W. Gerstner. Piecewise linear differential equations and integrate-and-fire neurons: Insights from two-dimensional membrane models. Phys. Rev. E, 67:021908, 2003.
• [19] S. Banerjee and G.C. Verghese, editors. Nonlinear Phenomena in Power Electronics. IEEE Press, New York, 2001.
• [20] Z.T. Zhusubaliyev and E. Mosekilde. Bifurcations and Chaos in Piecewise-Smooth Dynamical Systems. World Scientific, Singapore, 2003.
• [21] C.K. Tse. Complex Behavior of Switching Power Converters. CRC Press, Boca Raton, FL, 2003.
• [22] M. Sekikawa, N. Inaba, and T. Tsubouchi. Chaos via duck solution breakdown in a piecewise linear van der Pol oscillator driven by an extremely small periodic perturbation. Phys. D, 194:227–249, 2004.
• [23] B. Lindner, J. Garcia-Ojalvo, A. Neiman, and L. Schimansky-Geier. Effects of noise in excitable systems. Phys. Reports, 392:321–424, 2004.
• [24] R.E. Lee DeVille, E. Vanden-Eijnden, and C.B. Muratov. Two distinct mechanisms of coherence in randomly perturbed dynamical systems. Phys. Rev. E, 72:031105, 2005.
• [25] A.S. Pikovsky and J. Kurths. Coherence resonance in a noise-driven excitable system. Phys. Rev. Lett., 78(5):775–778, 1997.
• [26] L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni. Stochastic resonance. Rev. Modern Phys., 70(1):223–287, 1998.
• [27] G. Hu, T. Ditzinger, C.Z. Ning, and H. Haken. Stochastic resonance without external periodic force. Phys. Rev. Lett., 71(6):807–810, 1993.
• [28] J. Durham and J. Moehlis. Feedback control of canards. Chaos, 18:015110, 2008.
• [29] X. Li, J. Wang, and W. Hu. Effects of chemical synapses on the enhancement of signal propagation in coupled neurons near the canard regime. Phys. Rev. E, 76:041902, 2007.
• [30] G. Zhao, Z. Hou, and H. Xin. Canard explosion and coherent biresonance in the rate oscillation of CO oxidation on platinum surface. J. Phys. Chem. A, 109:8515–8519, 2005.
• [31] N. Yu, R. Kuske, and Y.X. Li. Stochastic phase dynamics and noise-induced mixed-mode oscillations in coupled oscillators. Chaos, 18:015112, 2008.
• [32] M. Desroches, B. Krauskopf, and H.M. Osinga. Mixed-mode oscillations and slow manifolds in the self-coupled Fitzhugh-Nagumo system. Chaos, 18:015107, 2008.
• [33] H.G. Rotstein, M. Wechselberger, and N. Kopell. Canard induced mixed-mode oscillations in a medial entorhinal cortex layer II stellate cell model. SIAM J. Appl. Dyn. Sys., 7(4):1582–1611, 2008.
• [34] B. Lindner and L. Schimansky-Geier. Coherence and stochastic resonance in a two-state system. Phys. Rev. E, 61(6):61