## Abstract

We consider the noise-induced escapes in an excitable system possessing a quasi-threshold manifold, along which there exists a certain point of minimal quasi-potential. In the weak noise limit, the optimal escaping path turns out to approach this particular point asymptotically, making it analogous to an ordinary saddle. Numerical simulations are performed and an elaboration on the effect of small but finite noise is given, which shows that the ridges where the prehistory probability distribution peaks are located mainly within the region where the quasi-potential increases gently. The cases allowing anisotropic noise are discussed and we found that varying the noise term in the slow variable would dramatically raise the whole level of quasi-potentials, leading to significant changes in both patterns of optimal paths and exit locations.

## 1. Introduction

The evolution of many dynamical systems involves a strong separation of time scales, leading to the interplay of fast and slow variables. Examples of interest arise in a wide variety of areas including planetary motions, geophysical flows, climate–weather interaction models and neural systems. A distinctive feature shared by many fast–slow systems, especially in neuroscience, is excitability. An excitable system may exhibit qualitatively disparate responses to slightly different perturbations in magnitudes, which suggests the existence of a threshold determining the generation of large-amplitude excursions. In other words, the essential features of excitability and occurrence of action potentials in many neurons lie in traversing through a threshold in a broad sense.

From the viewpoint of dynamical systems, a threshold, if one exists, is never a fixed number but a manifold manifesting itself as an invariant set, e.g. a single curve in two-dimensional systems. It was first noticed by FitzHugh [1] when he formulated a mathematical definition of firing thresholds using geometric analysis of neural models. However, in many cases, there is no well-defined threshold manifold due to the absence of saddle states. For many neural systems near a bifurcation from the resting state to sustained spiking activity, there may exist a region consisting of many trajectories corresponding to intermediate-amplitude responses ranging between subthreshold responses and stereotypical spikes. This was referred to as a fuzzy threshold set by Izhikevich [2], which, however, could be quite narrow in some models, including the Hodgkin–Huxley model [3]. In particular, the differences between stimuli engendering the respective responses of small potentials and action potentials could be so small that for any practical purposes such models exhibit all-or-none behaviours. In this sense, the threshold set could be regarded as a quasi-threshold manifold.

On the other hand, small random perturbations may have a profound impact on the dynamics and lead to the emergence of novel behaviours which could be rather deterministic in suitable limits [4]. Noise-driven excitable systems can produce responses with a high degree of coherence, with the leading phenomenon referred to as coherence resonance, which was first proposed by Pikovsky & Kurths [5]. They studied the extent of coherence resonance by a relation between the noise-induced activation time and the excursion time. Another well-known example describing the dramatic influence of small noise is stochastic resonance [6], which is when a system driven by a weak periodic excitation will oscillate exactly with the frequency of a driving force in the presence of vanishingly small noise without which no such oscillation would be expected. The immanent mechanism of stochastic resonance is also related to a match between the driving frequency and the escaping rate, which depends on the barrier height of a potential well [7]. Moreover, spatial coherence resonance extracted by additive Gaussian noise in excitable media was found to be best pronounced at an optimal noise level, which leads to an appropriate relation between the noise robust excursion time and the activation time, which is heavily subject to noise intensity [8]. Besides, a similar mechanism underlies other noise-induced phenomena such as spatial decoherence [9] and delay-aided stochastic multi-resonances [10].

The Freidlin–Wentzell large deviation theory [7] provides a proper framework to depict long-term effects of small random perturbations. Briefly speaking, the theory builds on the fact that almost unlikely events, when they occur, do so with an overwhelming probability in the way that is least improbable. This makes the prehistory of noise-induced transitions predictable by employing a certain action functional, whose minimum serves as an exponential rate of the stationary probability density in the Wentzel–Kramers–Brillouin (WKB) approximation [11] and gives the exponential order of the transition rate between metastable states. Moreover, the minimizing path of the action functional gives the pathway in whose vicinity the transitions occur with maximal probability. In past decades, a great deal of work has been devoted to problems of noise-induced transitions or escapes from equilibria [12], limit cycles [13,14] and chaotic attractors [15,16]. However, most works studied cases where a well-defined threshold manifold can be found. Recently, noise-induced escapes in an excitable neuronal model were studied by Khovanov *et al*. [17], who found that the quasi-threshold structure was revealed via the patterns of optimal fluctuational paths. They also compared the responses of a monostable resonator and integrator with stochastic input signals and a mixture of periodic and random stimuli. Based on that, Franovic *et al*. [18] investigated the activation process of both a single unit and multiple interacting excitable units by analysing the effects of two noise sources on statistical features of the action process and the modifications due to different forms of coupling. Meanwhile, cases of large assemblies of excitable units [19] were also considered by the mean-field approach and macroscopic excitable behaviours exhibited by the assembly were demonstrated explicitly.

In this paper, we investigate the excitable dynamics under weak noise in the FitzHugh–Nagumo (FHN) model in which there exists no well-defined threshold manifold. Compared with previous research, we attempt to provide some insights into the mechanism of excitability by using the concept of quasi-potential. Both the quasi-threshold manifold and the right branch of the N-shaped nullcline are regarded as boundary conditions for the problem of noise-induced escapes. We found that optimal paths exit through a certain point on the quasi-threshold manifold. Statistical analyses are then performed and we analyse the effect of finite noise by taking into account the distribution of the quasi-potential. We also change the noise terms added to the fast and slow variable and study the influence of various noise combinations. This paper is organized as follows. In §2, we formulate the problem of the FHN model and discuss some properties of its quasi-threshold manifold. Theoretical results of the large deviation theory providing foundations for our study are given in §3. Escaping scenarios under isotropic noise are investigated in §4, where the effect of small but finite noise is also discussed. The influence of anisotropic noise on the exit location along the quasi-threshold manifold is discussed in §5. Conclusions are then given in §6.

## 2. Formulation

We investigate the FHN model driven under a Gaussian white noise (*ξ*_{1}(*t*),*ξ*_{2}(*t*)) in the following form:
*x* and *y* represent the fast and slow variables, which correspond to the membrane voltage and the slow activation of an outward current, respectively. The parameter *ε* having small values delineates the separation in time scales between the fast and slow dynamics, and *a* controls the bifurcation of system (2.1), whereby the critical value is |*a*|=1. For |*a*|>1, system (2.1) has a unique stable equilibrium (*x*_{eq}, *y*_{eq})=(−*a*,−*a*+*a*^{3}/3), whereas for |*a*|<1 a supercritical Hopf bifurcation occurs, giving birth to a small-amplitude limit cycle attractor. In this paper, we consider the dynamics of system (2.1) with the parameters fixed: *a*=1.05, *ε*=0.05. The noise *ξ*_{1}(*t*) may account for synaptic noises due to random external inputs and the noise *ξ*_{2}(*t*) can be interpreted as internal thermal fluctuations in the opening of the iron-gating channels. They satisfy 〈*ξ*_{i}(*t*)〉=0 and 〈*ξ*_{i}(*t*)*ξ*_{j}(*s*)〉=*Dδ*_{ij}*δ*(*t*−*s*), where 〈 ⋅ 〉 denotes the mean value. The overall intensity of the noise is set by the small parameter *D*≪1 and the coefficients *r*_{x} and *r*_{y} determine the relative noise intensities of individual sources.

The phase space of system (2.1) in the absence of noise is shown in figure 1. The red and blue dashed lines denote the *x*- and *y*-nullcline, across which the derivatives *x*_{eq},*y*_{eq}) determined by the parameter *a*. The *x*-nullcline has a cubic N-shape and can be divided into three parts: the left, middle and right branch. As one can see from the small arrows representing the vector field **F** in figure 1, the left and right branches are attracting and the middle branch is repelling. In the singular limit *ε*→0, (2.1) turns into a one-dimensional system with *y* fixed, which means the middle branch of the *x*-nullcline forms a well-defined threshold manifold consisting of repellers.

However, for small but finite *ε* it has been shown [20] that a single line forming the threshold manifold does not exist, but a foliation of lines can be specified. That is to say, a layer, made up of infinite canard trajectories, plays the role of the threshold set in this case. Moreover, the layer is found to be quite thin and tends to one line outside the close vicinity of the top point P in figure 1. According to the discussion in §1, it is this quasi-threshold manifold [1] for finite *ε* that corresponds to the ‘ghost separatrix’ [17,21] in the singular limit. The quasi-threshold manifold is related to a special canard trajectory, depicted in figure 1 by a green solid line. This trajectory follows near the middle branch of the cubic nullcline and leads its way to the top point P. Note the highly unstable flow on both sides of this trajectory: a small perturbation pushing the initial state above or below results in a subthreshold or superthreshold response; see the thin cyan and black lines in figure 1. They are obtained by direct integration of system (2.1) with different initial values *y*_{0} and the differences between those labelled by 1 and 2 are smaller than 10^{−4}. Throughout this paper, we refer to this special canard trajectory as a ‘separatrix’. An easy way to compute the separatrix in two-dimensional relaxation oscillators is to integrate the system (2.1) backwards (

## 3. Large deviation theory

The Fokker–Plank equation describing the probability density for system (2.1) can be solved by the WKB approximation [7],
*P*(**x**) is the stationary probability density, *C*(**x**) is a prefactor beyond the scope of this paper and *S*(**x**) is the action calculated along a certain trajectory terminating at point **x**. To the leading order of intensity *D*, we obtain the following auxiliary Hamiltonian system:
*x*,*y*) and momentum (*p*_{x},*p*_{y}), the Hamiltonian flow given by system (3.2) determines the most probable (optimal) paths of system (2.1) connecting the initial **x**(0) and the final **x**(*T*) points. In fact, the large deviation theory gives a rough estimate for the probability that a trajectory *ψ*∈*C*[0,*T*], i.e.
*δ* and *D* sufficiently small. Here *P*_{x} denotes the probability under the condition **x**(0)=*x* and we assume that *ψ*(0)=*x*. If *B* is a Borel subset of **R**^{2}, we then have
*ε*→0. That is to say, the most probable paths are essentially minimizers to the action functional [7]:
*H*(**x**,**p**) via a Legendre transformation [7], **F** denotes the deterministic vector field and **A**=*Λ**Λ*^{T} (see ** Λ** in equation (2.1)). Accordingly, this functional extremal problem

*δS*/

*δx*=0 with given boundary conditions can be solved by the steepest descent methods, and through applying proper relations given in [22] we can also obtain the Hamiltonian system (3.2).

Two remarks should be made concerning equation (3.6). First, the physical meaning of the conjugate momentum **p** introduced in the Hamiltonian system (3.2) is specified. In fact, its absolute value measures the divergence of the most probable paths from the deterministic dynamics [23]. An approximation to zero implies that the fluctuational paths nearly follow the motions in the absence of noise. Thus, the momentum **p** can also be treated as an indicator signifying the effect of noise during the escaping process.

Second, *T* in equation (3.6) is finite. However, we remark that the infimum in equation (3.5) is usually not attained since the initial or terminal states may be equilibria or limit cycles, so that a path connecting to an equilibrium at one end takes infinite time to approach or deviate from it in theory [7]. In other words, the global minimum of the action functional is only achieved at *T*) denotes the space of all absolutely continuous functions *ψ*:[0,*T*]→**R**^{2} so that *ψ*(0)=**x**_{1} and *ψ*(*T*)=**x**_{2}. Similar to gradient systems, it is found that the expected exit time from the domain of attraction of a metastable state exponentially depends on the depth of the potential well [7], i.e.
*U*_{x1} denotes the basin boundary of the state **x**_{1}. In other words, the quasi-potential can measure the relative stability of each metastable state in non-gradient systems [24]. For more details one can refer to Freidlin & Wentzell [7] and, in this paper, it turns out to be a cardinal tool for our study.

## 4. Phenomena caused by isotropic noise

### (a) Escaping from the separatrix

In this section, we are going to consider the noise-induced phenomena of the FHN model (2.1) by isotropic sources, i.e. *r*_{x}=*r*_{y}=1. We consider two types of boundaries here, that is, the separatrix and the right branch of the *x*-nullcline (SR). The separatrix represents the quasi-threshold manifold of the system, which results in the extended problem of noise-induced escapes. The consideration of SR is related to the noise-induced spiking process, which is of great interest in neuroscience.

We calculate the quasi-potential of system (2.1) by the adjusted ordered upwind method (OUM) [25], which considers and accepts the values of adjacent points in ascending orders. A regular mesh of 30 000×30 000 points is used to cover the region [−4,4]×[−3,3], which is sufficiently accurate for our purposes. From the contour plot shown in figure 2*a*, one can see that there exists a large region between the left and right branches of the *x*-nullcline, where the values of the quasi-potential change extremely slowly and remain nearly constant. An elaboration on this will be given later, but here we only note that the contour lines start to accumulate sharply outside the *x*-nullcline. In figure 2*b*, we give a three-dimensional plot of the quasi-potential near the equilibrium to visualize the escape scenario from the potential well.

We first consider the extended escape problem by treating the separatrix as the boundary. It is parameterized by discretized points, at which the values of the quasi-potential are computed through cubic interpolations, as depicted in figure 3*a*. As one can see from the inset plot in figure 3*a*, the values of the quasi-potential remain around zero in comparison with those of the left end of the separatrix, i.e. those near point Q in figure 1, where they begin to increase dramatically. However, taking a closer look, we found that the quasi-potential has a non-monotone change hidden in the flat zero part, giving rise to a minimum value. The point corresponding to this minimum is then indicated by a black square on the separatrix in figure 3*b*, denoted as *S*_{0}. Equation (3.4) implies that this point determines the expected exit time and, furthermore, the escaping path fluctuating to it should be the optimal one in theory [7].

In order to verify the validity of our calculations, the geometric minimum action method [26] (gMAM) is used. Recalling the double minimization in equation (3.7), this method is based on a reformulation of the action functional *S*_{T}(*ψ*) on the space of curves so that the minimizing curve **x** may still exist. As proved in [26], the value of the action functional integrated along a certain curve is independent of the choice of its reparametrization *α*. We choose several points on the separatrix and find the minimizing curve connecting each one of them to the equilibrium **x**_{eq}=(−*a*,−*a*+*a*^{3}/3) by gMAM. Then the computed action corresponding to each minimizer is plotted by a green triangle in figure 3*a*, showing a rather good agreement with values given by the OUM.

Once we have found the quasi-potential, the velocity vector of the optimal path and the gradient of the quasi-potential are related via
**F** denotes the deterministic vector field of system (2.1). Subsequently, it leads us to shoot a flowline of −**G** from the first arrival at any point **x** back to **x**_{eq}. This flowline is necessarily the optimal path providing a global minimum to the action functional with given boundary conditions: **x**_{eq} and **x**. A rigorous proof can be found in [25].

The flowline from *S*_{0} back to the equilibrium is shown by a solid blue line in figure 3*b*. Meanwhile, the curve connecting both ends is also obtained through gMAM and plotted by a red solid line. Evidently, they can be hardly distinguishable from each other. However, one may query that these two paths are given by fixing both end points. Next, we cover this shortage by using the method of the action plot, which finds solutions to different initial value problems of the Hamiltonian system (3.2). Some 10^{5} initial points are uniformly distributed on a small circle of radius 10^{−5} centred at **x**_{eq}. The initial conditions of the four-dimensional system (3.2) at these points are constructed by the proper relations given in [27]. The Hamiltonian system is then integrated until its flow in the (*x*,*y*)-plane crosses the separatrix, and by that moment the cumulative action is conserved. We stress that we do not specify explicitly which point on the separatrix the projections of the Hamiltonian flow would traverse. In figure 3*b*, the Hamiltonian path corresponding to the global minimum action is plotted by a black dashed line, and it is in accord with the results of OUM and gMAM, escaping from the separatrix after traversing *S*_{0}.

### (b) Arriving at the nullcline

Previously, we have discussed the case of escaping from the separatrix along which there exists a point of minimum quasi-potential. Regarding it as an extended escape problem, we found that a certain point *S*_{0} lying on the quasi-threshold manifold assumes the role of saddles in the standard escaping scenarios. In these terms, the quasi-threshold manifold (separatrix) turns out to be an analogue of the well-defined threshold (stable manifold), with account taken of the rapid divergences of trajectories close by. Now let us go further and allow the escaping trajectories to reach the *x*-nullcline SR to see if our assertion holds.

To begin with, we must notice an essential difference between the separatrix and SR, namely the separatrix is an invariant set of the original deterministic system (2.1) while the SR is not. It follows that the optimal path would approach the SR neither tangentially nor asymptotically (i.e. *p*_{x},*p*_{y}→0). It is apparent if one recalls the fact that the momentum **p** actually measures the divergence from the deterministic dynamics. However, a distinguished feature of the action plot is that each point of discontinuity corresponds to a heteroclinic Hamiltonian trajectory connecting two stationary states of the original system. Usually one is an attractor and the other is saddle-type and lying on the basin boundary. We remark that the lowest action heteroclinic trajectory is the optimal path [13,14].

Keeping this in mind, we obtain the minimum action path by taking the SR as the boundary of the action plot; see the green solid line in figures 3*b* and 4*a*. It approaches the separatrix, passes through *S*_{0} and subsequently moves to the SR. We also found that these two optimal paths with different boundary conditions completely coincide with each other; see figure 3*b* for an illustration. More importantly, this coincidence does not depend on the radius of the small circle or number of initial points used in the action plot, therefore it is technically robust. Regarding the patterns, obviously the green line is only tangent to the separatrix rather than the SR. Its momenta *p*_{x} and *p*_{y} are shown in figure 5*a* by two thick solid lines in red and yellow, respectively, and both display the same tendency to 0. However, tending to 0 actually occurs while the optimal path is approaching *S*_{0}, since the plus sign ‘+’ in figure 5*a* corresponds to the point labelled by the same sign in figure 4*b*. After traversing *S*_{0}, the momenta *p*_{x} and *p*_{y} remain as zero and it then takes little time to fulfil the excursion and reach the SR.

To further demonstrate our discussions, we depict the vector field **G** given by equation (4.1) with green arrows and a flowline from *S*_{0} towards SR is also shown. Indicated by a thin blue line in figure 4*a*, it shows a good agreement with the optimal path given by the action plot. In conclusion, whether the boundary is the separatrix or SR, as *D*→0 the optimal path approaches the separatrix tangentially and asymptotically until exiting from a certain point *S*_{0} with the lowest quasi-potential along the separatrix.

Before ending this subsection, let us briefly discuss the time series of the momentum **p**. It can be found that *p*_{x} and *p*_{y} alternate to achieve their maxima, as shown by the squares and crosses in figure 5*a*. In fact, it stems from the spiral dynamics near **x**_{eq} if we recall the roles that the momenta play. The deterministic vector field **F** spirals into the equilibrium while **G** spirals out. Both of them are anticlockwise. When *p*_{x} attains its maxima, *F*_{y} and *G*_{y} are in the same direction (both upwards or downwards) but *F*_{x} and *G*_{x} are opposite (one left and the other right); see the red squares in figure 4*b*. Thus, the maxima of *p*_{x} correspond approximately to the zeros of *p*_{y}. The same is true for the maxima of *p*_{y}, denoted by the red crosses in figure 4*b*. Figure 5*b* describes the cumulative process of the action of the optimal path, the crosses in which correspond to those in figures 5*a* and 4*b*.

### (c) Effect of finite noise intensity

The Freidlin–Wentzell large deviation theory provides results in the weak noise limit, i.e. *D*→0. However, the noise intensity *D* in numerical simulations cannot be arbitrarily small. Now we are in a position to consider the effect of finite intensity on the escaping scenarios for the FHN model (2.1).

Statistical results are obtained by means of tracing out ridges of the prehistory probability distribution [28] *P*_{h}, which turns out to be an appropriate method to describe the distribution of fluctuational paths, both theoretically and experimentally [29]. This distribution contains all the information on the temporal evolution of the system before arriving at a given point, and it should peak sharply around the optimal path as *D*→0. Setting *D*=10^{−5}, we integrated system (2.1) by the second-order stochastic Runge–Kutta method [30], with initial conditions at **x**_{eq} and stopping criteria *x*>1.75. The time step for numerical simulations is Δ*t*=10^{−4}. The *P*_{h} is displayed in figure 6 with its ridges indicated by a yellow dashed line. The momenta *p*_{x} and *p*_{y} during the escape are also simulated by averaging the noise realizations conserved simultaneously with the paths. After performing a low-pass filtering process, we obtain the two zigzag lines as shown in figure 5*a*. The low-pass filtering process helps to reveal the underlying patterns of noise by eliminating their high-frequency components. As can be seen, they are in good agreement with *p*_{x} and *p*_{y} given by the action plot.

For comparison, we also plot the statistical result as a yellow line in figure 4*a*,*b*, which clearly shows that the statistical path crosses the separatrix and spikes before reaching *S*_{0}. This discrepancy from our theoretical analysis mainly derives from two aspects. First, as described by equation (3.4), the action functional essentially estimates the probability or difficulty that **x**(*t*) passes through the neighbourhoods of various paths *ψ*. Since the quasi-potential given in figure 2*a* has a large flat region between the left and right branches of the *x*-nullcline, the activation energies of different paths within this region are almost equal. Furthermore, due to the highly divergent dynamics around the separatrix, most realizations move immediately to the SR once they have crossed the separatrix. These contribute to the sharp ridges of *P*_{h}. Nevertheless, we find that the distribution appears not to change substantially for a weaker intensity. We should remark that further decreasing the noise intensity *D* would tremendously increase the mean exit time and the inducing escapes can hardly be observed within a reasonable length of time, so that constructing a reliable *P*_{h} is hardly practical.

A natural question arises: Is there anything special about the path along which the distribution peaks? A careful observation of the vector field **G** in figure 4*a* suggests that, after the Hamiltonian trajectories arrive at SR, they have two options: either moving right to infinity or moving left to the left branch of the *x*-nullcline. As a consequence, a critical trajectory that moves upwards vertically must exist and is found approximatively by both the method of the action plot and the OUM; see the thick black line and the thin red line in figure 4. The black line is given by a direct integration of the Hamiltonian system (3.2) from a certain initial condition, and the red line is given by shooting a flowline of −**G** from the small red triangle above the SR where the differences *p*_{x} and *p*_{y} of this particular path are plotted as two thin black lines in figure 5*a*. They evolve near the optimal path but do not tend to 0 at the end. Then the critical path has a slightly larger action than the optimal one, as shown in figure 5*b*. Its intersection points with the separatrix and SR are both indicated by black triangles in figure 4 and a Hamiltonian trajectory that crosses the SR below *N**_{0} will diverge right to infinity. It is found from figure 4 that under small but finite noise most realizations of system (2.1) concentrate around this critical path.

Next, we try to give some qualitative discussions regarding the behaviours of the statistical path, which crosses the separatrix at *S**_{0} below *S*_{0} shown in figure 4*b*. The quasi-potential at *S**_{0} equals 1.0600×10^{−4}, which is 6.7×10^{−7} larger than that of *S*_{0}. It is such a delicate difference that only sufficiently weak noise could trigger distinct dynamics. For the current intensity *D*=10^{−5}, they are hardly distinguishable. As we have mentioned above, for a smaller noise *D*≪10^{−5} physical observations would become impractical due to the astronomical order of the exit times. However, if we slightly increase the intensity, such as *D*=5×10^{−5}, the ridges of *P*_{h} shift downwards, shown as the magenta line in figure 6, in which we omit its statistical distribution for clarity. Further increasing *D* would not obviously change where the distribution peaks again, but only disperse the distribution, making its ridges unsharp. Both the black and magenta curves in figure 6 are plotted in figure 2*a* by two dash-dotted lines, which can be regarded as a rough partition of the contour plot. The quasi-potential above the black dash-dotted line, which corresponds to the critical Hamiltonian trajectory, remains nearly constant. Between the black and magenta lines the quasi-potential rises gently, and below the magenta line it increases rapidly. Therefore, under finite but small noise, in which case the large deviation theory applies, the exit locations on the separatrix depend on the noise intensity *D* and the Hamiltonian trajectory passing through that exit point determines the distributions of most realizations then.

Finally, due to the existence of a large flat region in the quasi-potential and the inherent features of the FHN system, the optimal path, as a global minimizer of the action functional, is no longer the most probable one in the sense of statistics under finite noise. However, it is this quasi-potential given by the large deviation theory that in turn provides a reasonable interpretation of the statistical results under small noise.

## 5. Influence of anisotropic noise

In this section, we make allowance for the cases that *r*_{x} and *r*_{y} are not equal and discuss the influence of this variation in brief. In order to calculate the quasi-potential in such a case, a transformation of coordinates is necessary since the OUM in [25] is proposed only for identical noise terms. Through a linear transformation *X* and *Y* denote stochastic processes and *W*_{1,2} represent two independent Wiener processes. Two remarks are given here. First, in the new coordinates the drift term of the original system (2.1) is changed. After computing the quasi-potential by the OUM, we can switch back to the original coordinates for subsequent studies. Second, due to the transformation we applied, the singular cases when one of *r*_{x} and *r*_{y} is zero is beyond the ability of this method.

We have calculated quasi-potentials for several combinations of intensity terms: fixing *r*_{x}=1, *r*_{y}=0.5,0.1,0.01 and fixing *r*_{y}=1, *r*_{x}=0.5,0.1,0.01, and plotted the points where the minima along the separatrix are located in figure 7. It can be seen that, as *r*_{y}→0, the positions of the minimum quasi-potentials are moving up markedly, but are moving down only slightly as *r*_{x}→0. Furthermore, as *r*_{y} decreases, the levels of the quasi-potential of the whole plane rise significantly, which suggests that the required actions for escapes from the potential well, in other words the difficulty, are highly raised. This can also be verified by the order of the magnitude of the momentum *p*_{x} for *r*_{x}=1,*r*_{y}=0, which is 10^{−2} as shown in the inset plot. By contrast, decreasing *r*_{x} while holding *r*_{y}=1 would not obviously change the values of the quasi-potential and *p*_{y} for *r*_{x}=0,*r*_{y}=1 is found to have the same order as the case *r*_{x}=1,*r*_{y}=1. Optimal paths for each case can also be obtained by similar procedures to those described in the previous section. In figure 7, we only show the results of gMAM by thin solid lines in corresponding colours. As *r*_{y} decreases with *r*_{x} fixed, the patterns of the optimal paths are found to change obviously while nothing occurs in the other case.

Therefore, we may draw the conclusion that the noise source added to the slow variable *y* plays an ascendant role in the escaping scenarios, thereby also in excitable ones, by a remarkable change of the magnitudes of the quasi-potential. We are convinced that this may provide certain insight into the occurrence of coherence resonance [5,31] when noises are typically added to the slow variable and the bifurcation parameter *a* tends to 1 from above. Nevertheless, substantial works are in order to fully verify this assertion.

## 6. Conclusion

In this paper, we studied the noise-induced escape in the FHN model, which has a quasi-threshold manifold consisting of a layer of infinite canard trajectories. We found that this quasi-threshold manifold, referred to as a separatrix, serves as a well-defined threshold often occupied by the stable manifold. Through computing the quasi-potential defined in the large deviation theory, a certain point is found to attain the minimum potential value along the separatrix. The optimal path in the limit *D*→0 approaches the separatrix tangentially and asymptotically until it reaches this minimum potential point. Then it moves to the right branch of the *x*-nullcline, spiking an action voltage. Furthermore, we demonstrated that this escaping scenario is independent of the choice of boundary condition. Various methods were employed, such as the gMAM, the OUM and the classical method of the action plot. All of them provide consistent results, which guarantee the validity of our discussions.

Numerical simulations under small but finite noise were then performed, giving rise to the prehistory probability distribution whose ridges describe the behaviours of most noise-driven realizations. Owing to the nearly constant values of the quasi-potentials and the highly diverging dynamics around the quasi-threshold manifold, the statistical path does not agree with the theoretical one in the limit of weak noise. In fact, the difference in the costed actions between them is much lower than the magnitude of the noise intensity used in the simulation, suggesting that such an order of intensity does not have a sufficiently high ‘resolution’ to distinguish these two paths, so that most realizations spike before reaching the minimum potential point in theory. Moreover, we found a critical path by both the method of the action plot and the OUM, which fits in well with the statistical path and manifests itself as a partition boundary in the contour of the quasi-potential. We found that the prehistory probability distributions of several noise intensities mainly concentrate within the region where the quasi-potential accumulates gently. As the noise intensity is raised the sharp ridges shift from the critical path down to the lower boundary, after reaching which the distributions flatten and the peaks unsharpen.

In addition, we briefly considered the cases where noises of different intensities were added to fast and slow variables. We found that varying the noise term in the slow variable would lead to a significant change in both the patterns of the optimal escaping paths and the locations of the minimum quasi-potentials along the quasi-threshold manifold. More importantly, decreasing the noise terms in the slow variable would dramatically raise the whole level of quasi-potentials, making noise-induced escape more difficult. As a consequence, it is the noise added to the slow variable that plays a dominant role in the course of excitability.

As a heuristic discussion to finish our paper, some potential to extend our approaches to coupled excitable units could be expected if certain global characteristics representing the ensemble states or patterns of large units are defined. Fortunately, it is rather promising as demonstrated in [19] that an assembly of excitable units may itself exhibit macroscopic excitable behaviours as a single unit. A low-dimensional approximate model of the original interactive system could be obtained by the mean-field approach and this then makes subsequent analysis of large deviation theory possible. However, to what extent this derivative model can approximate the original one, or whether it eliminates certain inner delicacy, requires further study.

## Data accessibility

The algorithm to calculate the quasi-potential can be found in [25] and necessary parameters for all the numerical results are shown in the captions of figures. Readers can perform the method to reproduce the present numerical results by any program language such as C, Fortran or Matlab.

## Competing interests

We have no competing interests.

## Authors' contributions

Z.C., J.Z. and X.L. designed the study and organized the entire research programme. Z.C. wrote the manuscript, Matlab code and performed the simulations. J.Z. offered considerable advice during the research. X.L. corrected and revised the manuscript, supervised the modelling and simulation. All authors gave final approval for publication.

## Funding

This work was supported by the National Natural Science Foundation of China (nos. 11472126, 11232007), the Research Fund of the State Key Laboratory of Mechanics and Control of Mechanical Structures (Nanjing University of Aeronautics and Astronautics) (grant no. MCMS-0116G01).

## Acknowledgements

We gratefully thank Lanmei Li and Ping Chen for their generous support during the writing of this paper.

- Received January 27, 2017.
- Accepted April 13, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.