NOISE-SENSITIVE MEASURE FOR STOCHASTIC RESONANCE IN BIOLOGICAL OSCILLATORS

There has been ample experimental evidence that a variety of biological systems use the mechanism of stochastic resonance for tasks such as prey capture and sensory information processing. Traditional quantities for the characterization of stochastic resonance, such as the signal-to-noise ratio, possess a low noise sensitivity in the sense that they vary slowly about the optimal noise level. To tune to this level for improved system performance in a noisy environment, a high sensitivity to noise is required. Here we show that, when the resonance is understood as a manifestation of phase synchronization, the average synchronization time between the input and the output signal has an extremely high sensitivity in that it exhibits a cusp-like behavior about the optimal noise level. We use a class of biological oscillators to demonstrate this phenomenon, and provide a theoretical analysis to establish its generality. Whether a biological system actually takes advantage of phase synchronization and the cusp-like behavior to tune to optimal noise level presents an interesting issue of further theoretical and experimental research.

1. Introduction.One of the remarkable nonlinear phenomena in biological systems is stochastic resonance (SR), which can roughly be characterized as the optimization of certain system performance by noise.Early theoretical evidence suggests that noise is essential for the transmission of sensory information, possibly through the mechanism of SR [1,2].There is also speculation that internal noise of a biological oscillator may play a constructive role in information transfer through SR [3].So far, there have been experimental demonstrations of SR in a variety of biological systems [4,5,6,7,8,9,10], including the interesting discovery that SR enhances the electrosensory information available to paddlefish for prey capture [10].There has even been a psychophysical experiment demonstrating that SR can be used as a measuring tool to characterize the ability of the human brain to interpret visual patterns immersed in noise [11].
SR was introduced in 1981 [12,13] as a plausible mechanism to account for the periodic occurrence of global glaciation (ice ages).It has stimulated a large amount of research, has been identified in a variety of natural and engineering systems, and 584 YING-CHENG LAI, AND KWANGHO PARK has continued to be an interesting topic in nonlinear science [14,15].Given a nonlinear system, its response to a weak signal is generally influenced by noise but, when SR occurs, noise can enhance the response.For periodic signals, a signalto-noise ratio (SNR) can be defined in terms of the dominant spectral peak in the frequency domain and can thus characterize the resonance in a natural way [16].For an aperiodic signal, there may be no well-defined peaks in the Fourier spectrum.In this case, the correlation between the input and the output signal [17,18,19], entropies, and other quantities derived from the information theory [20,21,22,23,24,25] can be used for characterization, where SR means the optimization of such a measure by noise.There are also nonlinear systems, in particular excitable systems, for which the performance optimization can occur in a range of the noise level.This is referred to as aperiodic stochastic resonance (ASR) [17,18,19].
A general feature associated with existing measures for characterization of SR is that either they vary slowly with noise about the optimal value, exhibiting a "bell-shape" behavior as typically seen in the SNR, or they are insensitive to noise variation (e.g., the correlation measure in the case of ASR [17,18,19]).Because of the ubiquity of SR in biological systems [4,5,6,7,8,9,10], a natural question is how a biological oscillator tunes to the optimal noise level to realize SR.For this purpose a measure that is highly sensitive to noise variation is desired.There are also potential technological applications where one might be interested in such a noise-sensitive measure.An example is to develop a device to assess the working environment based on the principle of SR [26,27,28,29].Specifically, for various types of measuring devices in a noisy environment, it is desirable to have the signal spectral peak as pronounced as possible with respect to the broad, noisy background.The principle of SR can naturally be used to detect the optimal noise level (or the optimal working condition).Andó and Graziani recognized that the SNR is in general not suitable for this purpose, as it does not allow for online tuning of the noise variance because of its insensitivity to noise variation about the optimal level.In a series of papers [26,27,28,29], they developed mathematical models utilizing feed-forward estimation theory and tested experimental devices based on the Schmitt trigger to overcome the difficulty.
In this paper, we present a measure of SR that has an extremely high sensitivity in the sense that, as a function of the noise amplitude, it exhibits a cusp-like behavior about the optimal noise level.In particular, recently it has been shown that SR can be understood as a manifestation of phase synchronization between the input and the output signal.The relationship between SR and phase synchronization has been demonstrated in noisy bistable systems [30,31,32,33,34,35,36,37] and in excitable systems with periodic [38,39] and with aperiodic signals [40].To see how SR can be studied through phase synchronization, imagine an input signal x in (t) that oscillates in time.A phase variable φ in (t) can be defined where one cycle of oscillation in x in (t) corresponds to an increase of 2π in φ in (t).A similar phase variable φ out (t) can be defined for the output signal x out (t).There is phase synchronization [41] if the phase difference satisfies ∆φ(t) ≡ |φ out (t) − φ in (t)| ≤ 2π for all t [42].In the presence of noise or due to the lack of coupling, phase synchronization can occur only in finite time intervals.That is, ∆φ(t) can remain bounded within 2π for a finite amount of time before a phase slip, typically 2π, occurs.Given a noise amplitude D, one can then measure the average time τ (D) for phase synchronization.The principal result of this paper is that, as the optimal noise level associated with SR is approached, this time increases so rapidly that mathematically, its behavior can be described as cusp-like.There is thus an extremely high sensitivity of τ (D) to noise variation, making it appealing for biological systems to achieve optimal-noise tuning or in device applications for detecting and realizing optimal working conditions1 .
We emphasize that, while many quantities have been proposed to characterize SR [14,15], all of them exhibit a bell-shape type of variation about the resonance and, hence, they are unsuitable for the type of applications mentioned above.In this sense, the average phase-synchronization time stands out as a measure that is characteristically different from the existing ones.A brief account of the cusp-like behavior of this time and a physical theory was published recently [44].The purpose of this paper is to extend this phenomenon to biological oscillators and provide more extensive theoretical and numerical support.In particular, we shall use one of the paradigmatic models for biological oscillators, the FitzHugh-Nagumo (FHN) system, and provide numerical evidence for the cusp-like behavior in the average phase-synchronization time (Sec.2).To explain the numerical finding (Sec.3), we shall use the theoretical approach in references [17,18,19] in the study of ASR to reduce the system of FHN oscillators to a minimal model (Sec.3.1): the double-well potential system that has been a paradigm to address many fundamental issues in SR.Using this idealized model with a simple periodic input signal, we can analyze the dynamics of phase synchronization by examining the stochastic transitions between the potential wells, leading to analytic formulas for τ (D) (Sec.3.2).Near the optimal noise level, τ (D) exhibits a cusp-like maximum with distinct values of derivative depending on whether the optimal level is approached from below or above.Although the specifics of τ (D) depend on the details of the system and the input signal, our analysis and further numerical evidence using the double-well potential model (Sec.4) suggest that the cusp-like behavior be general.
2. Array of FitzHugh-Nagumo oscillators.We consider the following array of FitzHugh-Nagumo (FHN) oscillators [45,46], a paradigmatic model for studying SR in biological oscillators: where 1, 0 < b < 1/2 are parameters, and S(t) is the input signal.To be as general as possible, we assume that the input signal is noisy: the term D J η(t) in equation (1) thus models this external noise term.The term Dξ i (t) simulates an adjustable noise source, suggesting that this is a type of internal noise source used by the biological oscillator for tuning to optimal noise level.In equation (1), η(t) and ξ i (t) (i = 1, . . ., N ) are independent Gaussian random processes of zero mean and unit variance: η(t)η(t ) = δ(t − t ) and ξ i (t)ξ i (t ) = δ(t − t ).The output from an FHN oscillator is typically a spike train.It is convenient to use the instantaneous firing rate, which is the number of spikes per unit time, to represent a smooth output signal for comparing with the input signal.For an array of FHN oscillators, the ensemble average firing rate can be used.A reasonably good match between S(t) and the firing pattern is achieved for moderate noise level (c).(e-h) The corresponding set of plots for rectangular input signal.
We present a series of numerical plots leading to the cusp-like behavior in the average synchronization time.We fix the following set of parameters in the FHN equation: = 0.005 and b = 0.15 (so that each FHN oscillator is subthreshold).To illustrate, we use two types of periodic input signals: (1) sinusoidal signal S(t) = Asin(ω 0 t) and ( 2) rectangular signal S(t) described by S(t) = A and −A for 0 ≤ t < T 0 /2 and for T 0 /2 ≤ t < T 0 , respectively, with the period given by T 0 = 2π/ω 0 .In our simulations we choose (arbitrarily) A = 0.045 and ω 0 = 0.15.The stochastic differential-equation system is numerically solved by using a standard second-order routine [47].Figures 1 (a-d) show, for the sinusoidal signal, the noisy input signal S(t) + 0.0015η(t) and the output spike trains x(t) for a single FHN oscillator for D = 0, 0.025, and 0.075, respectively.We see that without the internal noise (Fig. 1(b)), the output instantaneous firing rate r E (t) is low and it cannot represent the input signal S(t).For strong noise (Fig. 1(d)), the firing rate is high but its temporal variation does not match that of S(t), either.A reasonably good match  between r E (t) and S(t) is apparently achieved for moderate noise (Fig. 1(c)), where the firing behavior is significantly stronger in time intervals where the input signal is near its maxima.Similar behaviors occur for the rectangular input signal, as shown in a corresponding set of plots in figures 1 (e-h).For an array of FHN oscillators, for some appropriate noise level, there can be a nearly perfect match in the wave forms between the input signal and the average firing rate of the oscillators, as shown in figures 2(a-d), where the number of oscillators is N = 100.The phase variables associated with the noisy input signal S(t) + D J η(t) and the average output firing rate r E (t) can be calculated by using the standard Hilbert transform method.Figure 3 , respectively, where we observe an apparent cusp-like feature near the optimal noise level.

3.1.
Reduction of FHN model to double-well potential system.Heuristically, the dynamics of a single FHN oscillator can be reduced to the motion of a classical mechanical particle in a double-well potential system [17,18,19].For pedagogical purposes, we outline the major steps in the reduction process.Using the change of variables, x → x + 1/2 and y → y − b + 1/2, we can convert the dynamical equations for a single FHN oscillator to where A ≡ b − 1/2.For 1, the time rate of change of x(t) is much greater than that of y(t) and, hence, x(t) and y(t) can be regarded as a fast and a slow variable, respectively.Using the approximation ẏ ≈ 0 or y(t) ≈ x f (t) + Dξ(t), where x f (t) is a solution of the FHN system in the presence of signal S(t) but in the absence of noise, we can simplify the x-equation as where D ζ(t) ≡ D J η(t)−Dξ(t) represents the combined noise and D = D 2 J + D 2 is its amplitude.We have where the time-dependent potential function is given by which is a tilted double-well potential.The solution of the single FHN equations can then be interpreted as describing the motion of a heavily damped particle in the potential, with time-dependent slope of tilting.A firing event in a single FHN oscillator is equivalent to a crossing of the particle through the barrier.The ensemble-averaged firing rate r(t) is determined by the Kramers formula [48,49].Using perturbative analyses for x f (t) and to find the locations of the local minima of the potential well as well as the maximum of the barrier, one can obtain the following ensemble-averaged firing rate [17,18,19]: where B is a constant which is the "distance" of the system's excitation level to the threshold [17,18,19].
For the array of FHN oscillators in equation ( 1), the mean firing rate is the same as the ensemble-averaged firing rate of a single FHN oscillator.Fluctuations of the firing rate are determined by noise of the following form: D J η(t)+(D/N ) where σ(D ) is positive and proportional to D , and κ(t) is a Gaussian random signal.

Theoretical formulas of τ (D)
for the double-well potential system.We now consider particle motion in an idealized double-well potential in the presence of external driving and noise, subject to strong damping.The Langevin equation can be written as where ) is the potential function, D is the noise amplitude, and ξ(t) is the white noise term that satisfies ξ(t) = 0 and ξ(t)ξ(t ) = δ(t − t ).The potential has two wells, one at x l = −1 and another at x r = 1, and a barrier at x = 0. To make analysis feasible, we consider the case where the external driving F (t) is a periodically rectangular signal of period T 0 = 1, The tilted potential V (x, t) can assume one of two forms in the first and second half of a period, as shown schematically in figure 5. Due to the external forcing, the two wells become asymmetric with respect to each other.At a given time, one of the wells is deep and another is shallow, and they switch periodically in time with period T 0 .It is thus natural to assign a phase variable φ(t) for a particle in a well: φ(t) = 0 if the particle is in the well at x r and φ(t) = π if it is in the well at x l .Due to noise, a particle initially in one well can overcome the potential barrier to go to another well, vice versa.The rate of this transition (the probability of transition per unit time) is determined by the Kramers formula [48,49]: where E b is the barrier height.Let E d and E s be the barrier heights when the particle is in the deep and shallow well, respectively, where E d > E s .The rates of transition from the deep to the shallow well and the opposite are R d ∼ exp (−E d /D) and R s ∼ exp (−E s /D), respectively.For different noise strength, the response of the particle to the external forcing, as measured by the transitions between the two wells, determines the extent of phase synchronization between the input and the output.
Imagine that for some noise strength, in one period of time the probability of transition from the shallow to the deep well is appreciable but that of the opposite transition is negligible.Assume initially a particle is in the right well.During the first half-period the particle moves to the deep well, generating a π phase change.In the second half-period, the deep well becomes shallow and vice versa, and hence the particle moves to the right well again, as shown in figure 5.There is then another π phase change in the second half-period.The total phase change in one period is then 2π, which matches exactly the phase change associated with the external forcing (input signal).That is, the phase of the particle can be locked with respect to that of the input signal, giving rise to perfect phase synchronization.Because of noise, such a perfect synchronization cannot be achieved indefinitely.Let D opt be the noise amplitude for which the average synchronization time reaches a maximum value τ max 1.Let Ψ 2π (D) be the probability for a 2π change in one driving period.We have Consider first the case D < D opt .In the extreme case where D ≈ 0, the Kramers rates are essentially zero so that a particle initially in one potential well will remain there for a long time.Due to the 2π phase change in the input signal in one period, there will be a corresponding 2π change in the phase difference ∆φ between the input and the output signal.We have Ψ 2π (0) ≈ 1.As D is increased from zero, it becomes possible for a particle in the shallow well to move to the deep well so that R s will increase, but if D is small, we expect R d to remain negligible because of the higher potential barrier, as shown in figure 6.This will reduce Ψ 2π (D) from unity.The amount of reduction is given by the Kramers rate R s .The probability for 2π phase change is thus (1 − C 0 R s ), where C 0 is a constant that can be determined by the condition Ψ 2π (D opt ) ≈ δ.We obtain The average phase synchronization time for D < D opt is given by τ < (D) ≈ 1/Ψ < 2π (D).We have and . (Color online) For D > D opt , transition from the deep to the shallow well is possible but if this happens, there is a high probability that the particle will move back to the deep well in relatively short time.
For small value of δ, we thus expect to observe that τ < (D) increases rapidly as D → D opt from below.
For D > D opt , it is possible for a particle in the deep well to overcome the high potential barrier to move to the shallow well.Thus both transition rates become important.Because of the large noise, in one driving period a particle in the deep well can jump to the shallow well and then moves rapidly to the deep well again because of the relatively low potential barrier for this transition (R s ≈ 1).This process induces a 2π phase difference between the particle and the input signal.For larger noise, multiple transitions in one driving period are possible so that the phase difference would increase continuously with time and, as a practical matter, no phase synchronization can be observed.Since our interest is in the average synchronization time, only the noise range in which a single transition can possibly occur in one driving period is relevant.We thus have where the constant C 1 is given by The average phase synchronization time for D > D opt is given by τ > (D) ≈ 1/Ψ > 2π (D).Taking the derivative we obtain Again we observe that Equations ( 16) and ( 18) thus indicate a cusp-like behavior in τ (D) about D opt .Moreover, we have For small value of δ, for D → D − opt , the rise of τ < (D) can be pronounced than that of τ > (D) for D → D + opt .In general, we expect to see an asymmetric behavior in τ (D) near D opt .Note also that equation (17) implies that for D > D opt , the average time τ > (D) obeys the following scaling law with the noise amplitude: We emphasize that equations ( 16) and ( 18) indicate a fast rising and a fast falling behavior in the average synchronization time only for noise amplitude below and above the optimal value, respectively.Our argument is not applicable when the noise amplitude is in the infinitesimal vicinity of the optimal value.Thus, our heuristic theory cannot predict whether there is a cusp behavior in the mathematical sense of discontinuity in the derivative.A recent analytic expression [50,51] for the instantaneous phase-diffusion coefficient in a periodically driven system suggests, however, a smooth behavior in the phase-synchronization time about the optimal noise level.In particular, the diffusion coefficient shows a sharp but smooth peak at the optimal noise level.Since the average phase-synchronization time can be related to the inverse of the diffusion coefficient [52], it is reasonable that the behavior of this time also be smooth.

4.
Numerical results with the double-well system.Here we present numerical support for the cusp-like behavior in the double-well system with both periodic rectangular and sinusoidal driving.4.1.Periodic rectangular driving.By examining the average synchronization time, the optimal noise amplitude is determined to be D opt ≈ 0.033.Figures 8  (a-c) show, for driving amplitude F 0 = 0.18, the output signal x(t) with respect to the input for three values of the noise amplitude.For D = 0.02 < D opt (a), there are infrequent mismatches between the phases.For D = 0.033 ≈ D opt , we observe almost a perfect phase match between the input and the output signal in the time interval displayed.For D = 0.1 > D opt , phase mismatches occur quite often.These suggest that noise of amplitude D near D opt results in a maximal degree of phase synchronization between the input and the output signal.
Figure 9 shows, for the same driving amplitude, evolutions of the phase difference between the input and the output signal.Here the phase variable associated with the output is defined to be where x(t) is the Hilbert transform of x(t): and D = 0.025, respectively.Within the time considered (600 driving periods), we observe 2π phase slips for all noise levels except for D = 0.033, indicating that relatively long phase synchronization has been achieved and, hence, this is approximately the optimal noise amplitude, which has been observed to yield a maximum in the SNR.In this sense the measure of average phase-synchronization time is consistent with the traditional measures for characterizing SR.As D deviates away from D opt , 2π phase slips occur more often.
To verify the cusp-like behavior in the average phase-synchronization time τ , we choose a number of values of the noise amplitude about D opt .For each value, we use 20 realizations of the stochastic system to calculate the average value of the time between successive 2π phase slips.The result is shown in figure 10 (a), where we see that τ exhibits apparently a cusp-like behavior about the optimal noise amplitude D opt .As the noise amplitude is increased from zero and approaches the optimal value, there are four orders of magnitude of increase in the average phase-synchronization time, indicating an extremely high sensitivity to noise as compared with the traditional measures.The asymmetric behavior in τ about D opt , as predicted by our theoretical analysis, is shown in figure 10      nearly horizontal trace without large steps correspond to the case of D ≈ D opt .The cusp-like behavior in the average phase-synchronization time is shown in figure 14.
The cusp-like behavior in the average phase-synchronization time can also occur with respect to variation in the frequency of the external signal, the so-called bonafide stochastic-resonance phenomenon [53,54,55].To demonstrate this, we fix D = 0.031 ≈ D opt and calculate the average synchronization time τ as a function of the frequency of the input signal.The result is shown in figure 15, where τ exhibits a similar cusp-like behavior as in figures 10 and 14.Thus, the synchronization time is sensitive not only to noise variation, but also to other parameters such as the frequency of the driving signal.This may be interesting from the standpoint of frequency tuning in biological systems or in device applications.5. Discussion.Characterization of SR by using the phase-synchronization time may be of fundamental interest because this represents an alternative way to study SR.This approach can also be practically useful because the synchronization time depends on the noise level much more sensitively than the traditional measures such as the SNR.This may provide insights into the mechanism for biological systems to tune noise to achieve optimal performance through SR.In terms of technological applications, suppose an instrument is to be built based on the phenomenon of SR.Using this time measure can be more advantageous because of the higher precision it can potentially offer.In this paper, we have presented numerical evidence for the high noise sensitivity in a paradigmatic model of biological oscillators.To be able to obtain analytic understanding, we have used the standard double-well potential model with a periodic input signal.The dynamics of phase synchronization is then analyzed based on the transitions between the potential wells, with the help of the Kramers formula.Our principal finding is that, near the optimal noise level, the function τ (D) exhibits a cusp-like maximum with distinct values of derivative depending on whether the optimal level is approached from below or above.Although the specifics of τ (D) depend on the details of the system and the input signal, our analysis and numerical computations indicate that the cusp-like behavior is general.While our analysis is heuristic, a more rigorous treatment may be possible using a recently proposed two-state, discrete phase model for SR [36,37].Although there is a huge body of literature on SR, to our knowledge, the cusp-like behavior in the synchronization time has not been noticed previously.
Our approach to understanding stochastic resonance may also be useful for the phenomenon of resonant activation [56,57,58,59] where, for a particle in a potential well with a time-varying barrier, in the presence of noise the average crossing time, which is the average of times required to diffuse over each of barriers, can exhibit a minimum as a parameter controlling the barrier height varies.In previous works, the reported resonance peak is typically broad [56,57,58,59].Our results here imply that if resonant activation is treated using phase synchronization, it is possible that the average synchronization time can exhibit a cusp-like, sharp maximum.A possible setting to establish this is to assume that the barrier height is controlled by a time-varying signal (e.g., chaotic) for which a phase variable can be defined.The relative phase difference between the particle and this signal, and consequently phase synchronization, can then be investigated as in this paper.

Figure 1 .
Figure 1.(a) Noisy sinusoidal input signal S(t) + 0.0015η(t).(bd) Output spike trains x(t) from a single FHN oscillator for internal noise amplitude D = 0 (b), 0.025 (c), and 0.075 (d), respectively.A reasonably good match between S(t) and the firing pattern is achieved for moderate noise level (c).(e-h) The corresponding set of plots for rectangular input signal.

Figure 2 .
Figure 2.For an array of 100 FHN oscillators, match between the noisy sinusoidal input signal (a) and the average firing rate (b) for D = 0.025.Similar match occurs for the rectangular signal, as in (c) and (d).

Figure 3 .
Figure 3.Time evolutions of the phase difference between the input signal and the output average firing rate from a system of N FHN oscillators: (a) N = 3 for sinusoidal input signal and D = 0.06, 0.055, 0.025, 0.015, and 0.01 (from top down), (b) N = 4 for periodic rectangular input signal and D = 0.085, 0.07, 0.04, 0.025, and 0.02 (from top down).

Figure 4 .
Figure 4. Cusp-like behavior in the average phasesynchronization time about the optimal noise level: (a) for N = 3 and sinusoidal input signal, and (b) for N = 4 and periodic rectangular input signal.The synchronization times are expressed in units of the periods of the respective input signals.

N
i=1 ξ i (t), which can be written as D ζ (t), where D = D 2 J + D 2 /N and ζ (t) is also a Gaussian random signal of zero mean and unit variance.Taking into account random fluctuations in the ensemble-averaged firing rate, we have

Figure 5 .
Figure 5. (Color online) Tilted double-well potential as a result of external forcing.

Figure 6 .
Figure 6.(color online) For D < D opt , transition from the shallow to the deep well can occur but the opposite is unlikely.

Figure 8 .Figure 9 .Figure 10 .
Figure 8. (Color online) For the double-well potential system under periodic rectangular driving of amplitude F 0 = 0.18, the output signal x(t) in relation to the input driving for three values of the noise amplitude: (a) D = 0.02 < D opt , (b) D = 0.033 ≈ D opt , and (c) D = 0.1 > D opt .

4 . 2 .
(b), which is a blowup of part of figure 10 (a) near D opt .Support for the predicted scaling law (20) for D > D opt is shown in figure 11.Sinusoidal driving.We have obtained similar results when the input driving is a smooth sinusoidal signal A sin t.In particular, figures 12 (a-c) show, for A = 0.225, the output signal x(t) in relation to the input for three values of the noise amplitude: (a) D = 0.0175 < D opt , (b) D = 0.031 ≈ D opt , and (c) D = 0.075 > D opt .We see nearly perfect phase match between the input and the output signal for D ≈ D opt (panel (b)).
Figure 13  shows the evolutions of the phase difference between the input and the output signal for five different values of noise level.The

Figure 11 .
Figure 11.(Color online) For the double-well potential system under periodic rectangular driving, evidence for the scaling law(20).

Figure 12 .Figure 13 .
Figure 12. (Color online) For the double-well potential system under sinusoidal driving of amplitude A = 0.225, the output signal x(t) in relation to the input driving for three values of the noise amplitude: (a) D = 0.0175 < D opt , (b) D = 0.031 ≈ D opt , and (c) D = 0.075 > D opt .

Figure 14 .Figure 15 .
Figure 14.(Color online) For the double-well potential system under sinusoidal driving, cusp-like behavior in the dependence of the average phase-synchronization time on noise amplitude.