Stochastic Resonance in Nonpotential Systems

We propose a method to analytically show the possibility for the appearance of a maximum in the signal-to-noise ratio in nonpotential systems. We apply our results to the FitzHugh-Nagumo model under a periodic external forcing, showing that the model exhibits stochastic resonance. The procedure that we follow is based on the reduction to a one-dimensional dynamics in the adiabatic limit, and in the topology of the phase space of the systems under study. Its application to other nonpotential systems is also discussed.


I. INTRODUCTION
Stochastic resonance ( SR ) [1]- [8] is a phenomenon in which an enhancement of the response of a non-linear system is observed when this system is yielded to an external forcing at some optimized nonzero noise level. Since the original proposition of SR as a possible explanation for periodic recurrences in the global climate dynamics, SR has become the object of copious theoretical and experimental research in a wide variety of systems in physics, biology and chemistry. In all these works the possibility of noise having beneficial effects in the dynamics of nonlinear systems has been pointed out. The original formulation of the problem, in terms of a bistable system and a periodic forcing has been extended to systems under the action of aperiodic forcing [9] and nondynamical systems [7], [10].
In the present work we focus our attention on nonpotential systems. Non-potential systems correspond to systems far from equilibrium for which the principle of detailed balance does not hold. There are abundant examples of such systems arising from biological, chemical and physical problems. Our contribution in this paper is to develop a formalism which allows us to analytically treat a wide class of nonpotential systems among which one can include excitable and threshold systems as well as, for example, symmetric double-well models [11]. In particular, we apply our approach to continue the work undertaken by some authors in studying the stochastic properties of the FHN model. This is a well known model with wide application in the field of neuronal research [13], [5]. Apart from several numerical simulations done on this subject, Collins et al. [9] have carried out some analytical work on this matter in the area of aperiodic stochastic resonance . Some experimental research has been performed to show the existence of SR in this model [24]. The results obtained were compared to the predictions of the theory of SR in nondynamical systems [10]. Our scheme allows to analytically approach this problem in a simple way by using a generalization of the kinetic equations approach used in the case of potential systems ( see [2], [14], [15]).
All of the aforementioned models have a common feature; their dynamics exhibit three fixed points: an unstable point between two stable points. This feature established some resemblance between the processes described by these models and the hopping through a potential barrier. There are a great variety of systems that contains these features. The FitzHugh-Nagumo model, in its bistable regime, belongs to them. It is worth pointing out that this is not the regime in which this model is used to model neural activity. In this context the FHN model is taken to be in the excitable regime where only one globally attractor exists. As it is pointed out by Wiesenfeld et al. [6], a simple model of excitable system consists, among other things, of a threshold or potential barrier. Our theory provides a way to compute the escape rates from the attractors of a type of two-dimensional non-potential systems, and therefore it furnishes us with the main ingredient to apply the theory developed by Wiesenfeld et al.. In this fact you can find the relevance of our work to the field of excitable systems. Another example which fit the characteristics we are asking for is the class of symmetric double well systems [11]. The Sel'kov model [16] for autocatalytic systems gathers these features, too. This paper is organized as follows. In Section II we precisely define the range of applicability of our theory. We study the dynamics of the fluctuations and compute the kinetic equations. In Section III we introduce the FitzHugh-Nagumo model. We analyze its stochastic properties and show the existence of stochastic resonance. Finally in Section IV we discuss our main results.

II. DYNAMICS OF THE FLUCTUATIONS
In this work we study the stochastic properties of a class of two-dimensional nonpotential noisy dynamical systems. These system may be characterized by the peculiar topology of the phase space of their corresponding deterministic underlying versions. For the case we are concerned the dynamics is characterized by the presence of three aligned fixed points; an unstable point between two stable points. An example of this kind was given by Maier and Stein [11] in the context of the escape problem. They studied a system with a symmetric phase space consisting of an hyperbolic point between two sinks whose attraction basins were separate by the separatrix of the hyperbolic point.
To begin with consider a general two-dimensional noisy dynamical system [12] where x is the vector whose components are state variables, the field u( x) is the drift, g is the noise matrix and ξ(t) is a gaussian white noise of zero mean and correlation function given by with D being the noise intensity. For the case discussed in [11] the components of the drift are where P 3 (x) is a third order polynomial, and m is an integer such that 0 ≤ m ≤ 2. From the second of these equations, it is easy to check that all the fixed points are aligned. By equating eq. (3) to zero, one obtains a third order equation with three real roots for the proper values of the parameters. The corresponding Fokker-Planck equation for the probability density, ρ(x, t), is where D = D g · g T is the diffusion tensor.
Let us now assume that the system is potential. In such a case it is possible to write the Fokker-Planck equation as a continuity equation: where J is the diffusion current given by To obtain this expression, we have defined U and µ ( a generalized chemical potential) as follows For a potential system, the function U as defined in (8) is simply the potential energy. In the nonpotential case, however, the value of U will depend on the path of integration we choose and, in general, we can not achieve eqs. (6) and (7). Now consider the weak noise limit and think about some general characteristics of the probability distribution. If we let the system evolve during a sufficiently long time the probability distribution will have two maxima at the two stable fixed points ( SFP ) of the deterministic dynamics and a minimum in the unstable fixed point ( USFP ) [17]. On the other hand, in the weak noise limit the probability distribution will be very narrow around the line on which the fixed point lies. Thus we can assume that the fluctuations run over this line and, therefore, their dynamics are practically one-dimensional. This approximation can be justified accounting for the assumption of low noise level and the adiabatic limit. As well as in the one-dimensional double-well problem the adiabatic hypothesis implies that the representative point of the system is, in this long time limit, in one of the two wells [2], [14], in our two dimensional problem implies that the dynamics will be restricted to be on the nullcline. On the other hand, as it was found by Maier and Stein in [11], the distance to the nullcline, in the limit of weak noise intensity, is normally distributed with variance equal to the square root of noise level. Therefore, the fluctuations will run in a very narrow region around the nullcline. Consequently, although the whole system is nonpotential, we can reduce the dynamics being one-dimensional whose potential is given by the function U taking as the integration path the null cline. For the case discussed previously from eqs. (3), (4) and (8), under these circumstances ,one then obtains the drift corresponding to the dynamics of the fluctuations in the weak noise limit and Our next step will be to discretize the dynamics on the null cline . In particular we will obtain the kinetic equations. To this end we define the populations n + ( n − ) as the population on the right ( left ) of the USFP [18].
is the portion of the phase space on the right( left ) of a line orthogonal to the line which contains the fixed points and passing through the USFP .
In the adiabatic limit we can assume that the population is strongly concentrated in a small region around the SFP, as suggested by the picture of the probability density that we have profiled when the maxima in this long time limit is very high. This corresponds, in the assumption of one-dimensional fluctuations dynamics, having a bistable potential with two deep minima at the SFP and maximum at the USFP, or equivalently a high potential barrier.
In the present context we understand by adiabatic aproximation a long time limit in which all the system has arrived to a quasistationary state such that the probability of the system to be in a state different from the stationary stable states is practically zero.
So in this limit we can suppose that the system reaches a quasi-stationary state in which a quasi-stationary diffusion current is established. In addition, is assumed to be uniform between the two maxima of the probability density and, in the weak noise limit this current is concentrated in the line joining the three fixed points without loss of generality, the system can be taken to lie in the axis y = 0.
where x + ( x − ) is the coordinate of the fixed point on the right ( left ) of the USFP. The kinetic equation for n + is given by By using the divergence theorem and the assumptions about the form of the diffusion current, we have and, proceeding in the same way for n − , we obtain In addition, due to the height of the maxima in the probability density and the weakness of the noise, we can also consider that equilibrium in each side of the USFP is reached independently, thus, the generalized chemical potential is given by where x 0 is the coordinate of the USFP. The tensorial character of the generalized potential has been removed due to the dynamics reduction to one dimension. In equation (18) µ ≡ µ 11 has been defined. By using equation (18) in (7) we obtain where U corresponds to the integral over the adequate path of (10), U + and U − are its values at the SFP ( its minima ), Ψ ≡ g 11 ρ and Ψ ± = g 11 ρ(x ± , t).
In order to obtain the expression for the quasi-stationary current J 1 (t) we follow the same procedure as in [14]. From the definition of the probability current and the adiabatic hypothesis we have where a diagonal diffusion tensor is assumed. Integrating now over x and taking into account that, due to the height of the barrier, the main contribution to these integrals is around the maximum of the potential x 0 , one obtains where and µ + , µ − are the values of the generalized chemical potential at SFP. By using eq. (19) we can rewrite (21) in the following way where U + and U − are the values of the potential at the SFP. On the other hand, by using eq. (19) in eq. (12) one obtains the following relations These expressions allow us to write the kinetic equations for the populations n + and n − in the form where the kinetic coefficients are given by With this result we have obtained of the kinetic equations for the nonpotential system.

III. THE FITZHUGH-NAGUMO MODEL
In this section, we will apply the results of the previous section to the study of the FitzHugh-Nagumo ( FHN ) neural model [13], [19]- [21]. This model is a variant of the Hodgkin-Huxley model [21]- [23] which accounts for the essentials of the regenerative firing mechanisms in nerve cells. The FHN equations correspond to an excitable threshold model but, as will be seen briefly, due to their cubic non-linearity, they exhibit the characteristic behavior of a bistable system. Our main objective is to show analytically the appearance of SR in this model under a periodic external forcing.
The non-dimensional equations of the FHN model are [21] dv where 0 < a < 1 is essentially the threshold value, b and γ are positive constants and I a is the applied current. For the sake of simplicity, and without loss of generality, we will take I a = 0. The drift field for this model is given by As can be seen from eq. (28) the null cline of the deterministic dynamics of this equations is the line v = γ b w. By substitution on the right hand side of equation (27) This is a third order equation, which for certain values of the parameters has three roots ( see Figure 1 ). Among these three fixed points two are stable: F − and F + , and one unstable: F 0 , situated between the other two [21]. Thus, in this case, this system fulfills the conditions in order to our theory to be applied. The function U defined in Section II has to be integrated in this case over line v = γ b w. Performing this integration one obtains On the other hand, when this system is in a noisy environment, in the limit of weak noise, we can approximate the dynamics of the fluctuations by the one-dimensional Langevin equation that is, the fluctuations run along the line v = γ b w. As can be easily checked, U is the potential for the deterministic part of this equation. The two SFP of this one-dimensional dynamics are the two minima of (32) and the USFP is its maximum. Collins et al . [9] have arrived previously to similar conclusions by another approach in the context of the study of aperiodic stochastic resonance. Fig. 2 shows the asymmetric shape of U . Before going any further, it would be interesting to summarize what this picture has to say to us about the physics of the problem [13], [9]. In the FHN model there is a fast variable, v(t), and a recovery-like variable, w(t). After the barrier threshold has been crossed, i.e. the system has gone to an "excited" state, the system returns ( even in the deterministic case ) to the "rest" state. As can be seen in Fig. 2, there is one of the stable states for which the potential is larger than for the other stable state. Therefore, there is a more stable state, which corresponds to the rest state to which the system after some time , under the action of the noise, returns. In our scheme the presence of noise is necessary in order to return to the rest state, because of the elimenation of w(t).
In order to show how this scheme can account for the existence of stochastic resonance in the FHN model, we will suppose that the system is under the action of periodic forcing. For simplicity's sake, we will assume that the parameter a is a periodic function: a = a 0 (1 + ǫ 0 sinω s t) where ǫ 0 is a small parameter and a 0 (1 + α) < 1.
To take a as an oscillatory factor implies that the positions of F 0 ( the USFP ) and F + ( one of the SFP ) as well as the values of the potential at these points become oscillatory functions, too. The position of F − (the other SFP ) remains constant. Let v 0 , v − and v + be the v-coordinate of the maximum and the minima, respectively, of the potential; one has To compute the moments and the power spectrum, we assume that, in the limit of weak noise, the probability density ( in one dimension ) can be written as [2] p(v, t) = n where n + and n − come from the kinetic equations (25). The formal solution to these equations is found to be with In order to calculate n ± we perform a Taylor expansion of the transition rates up to the first order in respect to the parameter ǫ(t).
where φ 0 has been defined as φ 0 ≡ ǫ 0 /D and α ± 0 and α ± 1 are given by with d ± 0 being the zero order contribution of |U ′′ 0 |U ′′ ± , ξ 0 and η 0 the zero and first order contribution of U at the USFP and ξ ± and η ± the zero and first order contribution of U at the SFPs. Its first order contribution can be neglected in the limit of weak noise. By introducing these Taylor expansions in equations (36) and (37), we get the populations up to the first order in the parameter φ 0 .
where Φ ≡ arctg(α 0 /ω 0 ) and α 0 = α + 0 + α − 0 . The quantity n ± (t|v t0 , t 0 ) is the conditional probability that v(t) is in the + state at time t given that the state at time t 0 was v t0 . From this equation it is possible to compute the statistical properties of the process v(t). Of particular interest to our purposes is to find its autocorrelation function which is given by ( The conditional correlation function is given by Let us make some considerations which allows us to simplify the computation of the autocorrelation function. It is clear from equation (43) that the Fourier transform of the autocorrelation function will depend on t as well as in the frequency. This dependency is avoided by taking its average over the period of the external forcing [2]. The autocorrelation function is to be computed up to the second order in parameter φ 0 ∼ D −1 .Thus, in the limit of the weak noise D −1 , can be neglected in comparision to D −2 . Therefore, the only contribution of the first order term of v ± to the autocorrelation function comes from its product with the zero order term of the product of n's. But, on doing the average this term vanishes. So, finally, and taking into account that v − = 0, we arrive to v(t)v(t + τ )|v t0 , t 0 = λ 2 + n + (t + τ |v + , t)n + (t|v t0 , t 0 ), where the overlined indicates an average over t and λ + is the position of F + up to order zero. From equations (43) and (44) taking the average and the limit t 0 → −∞, we finally obtain the following expression for the autocorrelation function.
With this results we can now compute the averaged power spectrum given by where the last equality follows from the commutative character of averaging and Fourier transform. We obtain, after Fourier transform (50) In this expression the fraction of the total power in the broadband noisy part of the spectrum, which usually is a small fraction of the total power, has been neglected [2]. Note that in the power spectrum there is a term proportional to δ(Ω). This is due to the asymmetry of the potential that originates a mean probability current between one stable state and the more stable one is more stable than the other one. This term has not been obtained in previous approaches to this problem [9], [10] From equation (47) the signal to noise ratio, R, can be obtained as a function of the noise level D, by making Ω = ω 0 .
We have plotted R in Figure 3 as a function of the noise level. It exhibits a maximum for a non-zero noise level and therefore this model shows stochastic resonance.

IV. CONCLUSIONS
In this paper we have proposed a method to analyze the possibility for the appearance of SR in nonpotential systems. In particular we have treated the FHN model. A lot of numerical work on this subject has been performed but only a few papers have treated this problem from an analytical perspective due to the difficulties inherent to the nature of nonpotential system. The work by Collins et al [9] where aperiodic SR is discussed, and the paper by Pei et al. [24], in which the theoretical framework developed by Gingl et al. [10] was used to interpret some experimental results, are among the ones discussing analytical aspects.
Far from being specific of this problem, the method we have proposed can be extended to a wide variety of different nonpotential systems ranging from threshold systems to some autocatalytic models and symmetric double well systems from fields so diverse as chemical kinetics and population dynamics. Although, for the sake of brevity, we do not quote our results here we have applied our approach to that of Collins et al. [9] with a periodic input and we obtain similar results in both cases, not the same because in the case of these authors the input is additive and the threshold is maintained constant. We have applied our theory to the standard double-well model [11] and we have obtain that this model exhibits SR. The result in this case is equivalent to the one obtained in [2] for a symmetric quartic potential.
The scheme we present in this paper allows us to treat potential and the class of aforementioned nonpotential systems in an unified way. Moreover, it reproduces the essentials of the physics of the problem as can be seen from the obtaining of the refractory current in eq. (52) [13].
Finally, it is worth pointing out that our method, which constitutes essentially an extension of the Kramers rate theory to this kind of nonpotential systems, enables one to compute the kinetic coefficients in a simple and direct way.
• Figure 2.-Asymmetric one-dimensional potential as a function of v. Non-dimensional variables. •